diff --git a/cime_config/machines/cmake_macros/gnu_eddi.cmake b/cime_config/machines/cmake_macros/gnu_eddi.cmake index 009de57acb5d..be450901b2cb 100644 --- a/cime_config/machines/cmake_macros/gnu_eddi.cmake +++ b/cime_config/machines/cmake_macros/gnu_eddi.cmake @@ -1,12 +1,12 @@ +string(APPEND CMAKE_CXX_FLAGS " -fpermissive") if (COMP_NAME STREQUAL gptl) - string(APPEND CPPDEFS " -DHAVE_VPRINTF -DHAVE_GETTIMEOFDAY -DHAVE_BACKTRACE") + string(APPEND CPPDEFS " -DHAVE_VPRINTF -DHAVE_TIMES -DHAVE_GETTIMEOFDAY -DHAVE_BACKTRACE") endif() string(APPEND CMAKE_Fortran_FLAGS_RELEASE " -fno-unsafe-math-optimizations ") -string(APPEND CMAKE_Fortran_FLAGS_DEBUG " -g -fbacktrace -fbounds-check -ffpe-trap=invalid,zero,overflow -Wall") +string(APPEND CMAKE_Fortran_FLAGS_DEBUG " -g -fbacktrace -fbounds-check -ffpe-trap=invalid,zero,overflow -Wsurprising -finit-real=snan") #-fsanitize=thread") string(APPEND CMAKE_EXE_LINKER_FLAGS " -lcurl") if (MPILIB STREQUAL mpi-serial) set(SCC "gcc") -endif() -if (MPILIB STREQUAL mpi-serial) set(SFC "gfortran") endif() + diff --git a/cime_config/machines/config_machines.xml b/cime_config/machines/config_machines.xml index 79b1dd55ee13..1c20ba7ea233 100644 --- a/cime_config/machines/config_machines.xml +++ b/cime_config/machines/config_machines.xml @@ -4939,15 +4939,16 @@ small developer workhorse at lbl climate sciences + eddi LINUX gnu openmpi ngeet - /raid1/lbleco/e3sm/ - /home/rgknox/Models/InputDatasets/cesm_input_data/ - /home/rgknox/Models/InputDatasets/cesm_input_data/atm/datm7/ - /home/rgknox/Models//cesm_archive/$CASE - /home/rgknox/Models//cesm_baselines/$COMPILER + /home/rgknox/Models/cime_output + /home/rgknox/Models/InputDatasets/e3sm_input_datasets/ + /home/rgknox/Models/InputDatasets/e3sm_input_datasets/atm/datm7/ + /home/rgknox/LandRuns/E3SM/e3sm_archive/$CASE + /home/rgknox/LandRuns/E3SM/e3sm_baselines/$COMPILER /raid1/lbleco/cesm/cesm_tools/cprnc/cprnc 1 none @@ -4961,14 +4962,14 @@ mpirun - -np {{ total_tasks }} - -npernode $MAX_MPITASKS_PER_NODE + -np {{ total_tasks }} + --map-by ppr:{{ tasks_per_node }}:socket:PE=$ENV{OMP_NUM_THREADS} --bind-to core - /usr/local/share/cmake-3.21/ - $ENV{NETCDF_HOME} + $SHELL{nc-config --prefix} + $SHELL{nf-config --prefix} diff --git a/cime_config/machines/config_pio.xml b/cime_config/machines/config_pio.xml index 6687b995dc21..7d241fb48341 100644 --- a/cime_config/machines/config_pio.xml +++ b/cime_config/machines/config_pio.xml @@ -61,7 +61,7 @@ netcdf netcdf netcdf - netcdf + netcdf netcdf netcdf netcdf diff --git a/components/data_comps/datm/src/datm_comp_mod.F90 b/components/data_comps/datm/src/datm_comp_mod.F90 index 9b33c02190e4..0639e4e86f13 100644 --- a/components/data_comps/datm/src/datm_comp_mod.F90 +++ b/components/data_comps/datm/src/datm_comp_mod.F90 @@ -1081,8 +1081,11 @@ subroutine datm_comp_run(EClock, x2a, a2x, & rtmp = maxval(avstrm%rAttr(stdew,:)) call shr_mpi_max(rtmp,tdewmax,mpicom,'datm_tdew',all=.true.) endif - if (my_task == master_task) & - write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax + if (my_task == master_task) then + write(logunit,*) trim(subname),' tbotmax: ',tbotmax + !write(logunit,*) trim(subname), 'tdewmax: ',tdewmax,stdew + write(logunit,*) trim(subname), 'anidrmax: ',anidrmax + end if endif lsize = mct_avect_lsize(a2x) do n = 1,lsize diff --git a/components/elm/bld/ELMBuildNamelist.pm b/components/elm/bld/ELMBuildNamelist.pm index ebde3a5d0093..f080424d463e 100755 --- a/components/elm/bld/ELMBuildNamelist.pm +++ b/components/elm/bld/ELMBuildNamelist.pm @@ -3764,20 +3764,26 @@ sub add_default { # query the definition to find out if the variable is an input pathname my $is_input_pathname = $definition->is_input_pathname($var); + - # The default values for input pathnames are relative. If the namelist - # variable is defined to be an absolute pathname, then prepend - # the E3SM inputdata root directory. - if (not defined $settings{'no_abspath'}) { - if (defined $settings{'set_abspath'}) { - $val = set_abs_filepath($val, $settings{'set_abspath'}); - } else { - if ($is_input_pathname eq 'abs') { - $val = set_abs_filepath($val, $inputdata_rootdir); - } - } + if ($is_input_pathname eq 'landroot') { + my $landroot = abs_path("$ProgDir/.."); + $val = set_abs_filepath($val,$landroot); + } else { + # The default values for input pathnames are relative. If the namelist + # variable is defined to be an absolute pathname, then prepend + # the E3SM inputdata root directory. + if (not defined $settings{'no_abspath'}) { + if (defined $settings{'set_abspath'}) { + $val = set_abs_filepath($val, $settings{'set_abspath'}); + } else { + if ($is_input_pathname eq 'abs') { + $val = set_abs_filepath($val, $inputdata_rootdir); + } + } + } } - + # query the definition to find out if the variable takes a string value. # The returned string length will be >0 if $var is a string, and 0 if not. my $str_len = $definition->get_str_len($var); diff --git a/components/elm/bld/namelist_files/namelist_defaults.xml b/components/elm/bld/namelist_files/namelist_defaults.xml index 71a404db857f..d056cb8780ab 100644 --- a/components/elm/bld/namelist_files/namelist_defaults.xml +++ b/components/elm/bld/namelist_files/namelist_defaults.xml @@ -131,10 +131,10 @@ attributes from the config_cache.xml file (with keys converted to upper-case). lnd/clm2/paramdata/clm_params_c180524_phs.nc - - - -lnd/clm2/paramdata/fates_params_api.41.0.0_14pft_c250813.nc + + + +src/external_models/fates/parameter_files/fates_params_default.json lnd/clm2/paramdata/CNP_parameters_c131108.nc diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml index 9886773ee45a..6755818936b1 100644 --- a/components/elm/bld/namelist_files/namelist_definition.xml +++ b/components/elm/bld/namelist_files/namelist_definition.xml @@ -531,7 +531,7 @@ constants for biogeochem modules + input_pathname="landroot" group="elm_inparm" valid_values="" > Full pathname datafile with fates parameters diff --git a/components/elm/src/biogeophys/CanopyFluxesMod.F90 b/components/elm/src/biogeophys/CanopyFluxesMod.F90 index 1c278c75ce97..c82ddc7eb436 100644 --- a/components/elm/src/biogeophys/CanopyFluxesMod.F90 +++ b/components/elm/src/biogeophys/CanopyFluxesMod.F90 @@ -42,7 +42,9 @@ module CanopyFluxesMod use ColumnDataType , only : col_es, col_ef, col_ws use VegetationType , only : veg_pp use VegetationDataType , only : veg_es, veg_ef, veg_ws, veg_wf - + use decompMod , only : fates_pproc + use FatesInterfaceTypesMod, only : bc_out_type + use omp_lib !!! using elm_instMod messes with the compilation order use elm_instMod , only : alm_fates, soil_water_retention_curve use TemperatureType , only : temperature_vars @@ -52,10 +54,11 @@ module CanopyFluxesMod ! ! !PUBLIC TYPES: implicit none - save + ! ! !PUBLIC MEMBER FUNCTIONS: public :: CanopyFluxes + public :: PatchLoadBalance contains @@ -130,8 +133,6 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & integer :: time ! ! !LOCAL VARIABLES: - real(r8), pointer :: bsun(:) ! sunlit canopy transpiration wetness factor (0 to 1) - real(r8), pointer :: bsha(:) ! shaded canopy transpiration wetness factor (0 to 1) real(r8), parameter :: btran0 = 0.0_r8 ! initial value real(r8), parameter :: zii = 1000.0_r8 ! convective boundary layer height [m] real(r8), parameter :: beta = 1.0_r8 ! coefficient of convective velocity [-] @@ -210,7 +211,7 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & real(r8) :: wtalq(bounds%begp:bounds%endp) ! normalized latent heat cond. for air and leaf [-] real(r8) :: wtgaq ! normalized latent heat cond. for air and ground [-] real(r8) :: el(bounds%begp:bounds%endp) ! vapor pressure on leaf surface [pa] - real(r8) :: deldT ! derivative of "el" on "t_veg" [pa/K] + real(r8) :: deldT ! derivative of "el" on "t_veg" [pa/K] (dummy) real(r8) :: qsatl(bounds%begp:bounds%endp) ! leaf specific humidity [kg/kg] real(r8) :: qsatldT(bounds%begp:bounds%endp) ! derivative of "qsatl" on "t_veg" real(r8) :: e_ref2m ! 2 m height surface saturated vapor pressure [Pa] @@ -243,11 +244,14 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & real(r8) :: o2(bounds%begp:bounds%endp) ! atmospheric o2 partial pressure (pa) real(r8) :: svpts(bounds%begp:bounds%endp) ! saturation vapor pressure at t_veg (pa) real(r8) :: eah(bounds%begp:bounds%endp) ! canopy air vapor pressure (pa) + real(r8) :: rssun_old(bounds%begp:bounds%endp) ! used for determining convergence via change in resitance + real(r8) :: rssha_old(bounds%begp:bounds%endp) ! from one iteration to the next real(r8) :: s_node ! vol_liq/eff_porosity real(r8) :: smp_node ! matrix potential real(r8) :: smp_node_lf ! F. Li and S. Levis real(r8) :: vol_liq ! partial volume of liquid water in layer integer :: itlef ! counter for leaf temperature iteration [-] + integer :: itstoma ! counter for stoma iteration [-] integer :: iter_final ! number of iterations used integer :: nmozsgn(bounds%begp:bounds%endp) ! number of times stability changes sign real(r8) :: w ! exp(-LSAI) @@ -267,10 +271,8 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & integer :: filterp_noveg(bounds%endp-bounds%begp+1) ! bare ground pft filter integer :: fn ! number of values in vegetated pft filter integer :: filterp(bounds%endp-bounds%begp+1) ! vegetated pft filter - integer :: fnorig ! number of values in pft filter copy - integer :: fporig(bounds%endp-bounds%begp+1) ! temporary filter - integer :: fnold ! temporary copy of pft count integer :: f ! filter index + logical :: in_filterp(bounds%begp:bounds%endp) ! tells us if patch p is in filterp logical :: found ! error flag for canopy above forcing hgt integer :: index ! patch index for error real(r8) :: egvf ! effective green vegetation fraction @@ -318,12 +320,41 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & real(r8) :: slope_rad, deg2rad character(len=64) :: event !! timing event + logical :: converge_stoma ! logical switch that flags if the tveg loop converged + logical :: converge_tveg ! logical swithc that flags if the stomatal loop converged + real(r8) :: del_gs ! The maximum difference in stomatal conductance + ! from current iteration to previous, between sunlit and + ! shaded portions of the leaves [m/s] + + real(r8) :: reldel_gs + real(r8) :: hmean_gs + real(r8) :: t_end,t_start,w_end,w_start + real(r8) :: w_diff(0:fates_pproc-1) + real(r8) :: setup_overhead + + ! Indices for patch parallel threading + integer :: it ! thread iterator + integer :: np ! number of patches on this thread + integer :: totp(0:fates_pproc-1) + integer :: npatch(0:fates_pproc-1) + ! Indices for raw and rah integer, parameter :: above_canopy = 1 ! Above canopy integer, parameter :: below_canopy = 2 ! Below canopy ! Lower bound for VPD (based on CLM) real(r8), parameter :: vpd_min = 50._r8 + + ! We set the minum allowable difference in the conductance iteration + ! to be equal to the maximum allowable stomatal resistance (this number is from fates) + real(r8),parameter :: max_del_gs = 1._r8/2.e8_r8 ! [m/s] + real(r8),parameter :: max_reldel_gs = 0.02 + + integer, parameter :: itmax_stomata = 10 + + logical, parameter :: do_b4b = .false. ! Set this true to reproduce results before + ! refactoring the patch-loops + !------------------------------------------------------------------------------ associate( & @@ -431,6 +462,8 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & rssun => photosyns_vars%rssun_patch , & ! Output: [real(r8) (:) ] leaf sunlit stomatal resistance (s/m) (output from Photosynthesis) rssha => photosyns_vars%rssha_patch , & ! Output: [real(r8) (:) ] leaf shaded stomatal resistance (s/m) (output from Photosynthesis) + !rssun_old => photosyns_vars%rssun_old_patch , & ! Output: [real(r8) (:) ] previous leaf sunlit stomatal resistance (s/m) (output from Photosynthesis) + !rssha_old => photosyns_vars%rssha_old_patch , & ! Output: [real(r8) (:) ] previous leaf shaded stomatal resistance (s/m) (output from Photosynthesis) grnd_ch4_cond => ch4_vars%grnd_ch4_cond_patch , & ! Output: [real(r8) (:) ] tracer conductance for boundary layer [m/s] @@ -465,17 +498,14 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & begp => bounds%begp , & endp => bounds%endp & ) - if (use_hydrstress) then - bsun => energyflux_vars%bsun_patch ! Output:[real(r8) (:) ] sunlit canopy transpiration wetness factor (0 to 1) - bsha => energyflux_vars%bsha_patch ! Output:[real(r8) (:) ] sunlit canopy transpiration wetness factor (0 to 1) - end if + ! Determine step size dtime = dtime_mod !yr = year_curr; mon = mon_curr; day = day_curr; time = secs_curr; irrig_nsteps_per_day = ((irrig_length + (dtime - 1))/dtime) ! round up - ! First - set the following values over points where frac vegetation covered by snow is zero + ! First - set the following values over points where frac vegetation NOT covered by snow is zero ! (e.g. btran, t_veg, rootr, rresis) do fp = 1,num_nolakeurbanp p = filter_nolakeurbanp(fp) @@ -504,11 +534,13 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & ! Filter patches where frac_veg_nosno IS NON-ZERO ! ----------------------------------------------------------------- fn = 0 + in_filterp(:) = .false. do fp = 1,num_nolakeurbanp p = filter_nolakeurbanp(fp) if (frac_veg_nosno(p) /= 0) then fn = fn + 1 filterp(fn) = p + in_filterp(p) = .true. end if end do @@ -734,7 +766,6 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & wind_speed0(p) = max(0.01_r8, hypot(forc_u(t), forc_v(t))) wind_speed_adj(p) = wind_speed0(p) ur(p) = max(1.0_r8, sqrt(wind_speed_adj(p)**2 + ugust(t)**2)) - prev_tau(p) = tau_est(t) else ur(p) = max(1.0_r8,sqrt(forc_u(t)*forc_u(t)+forc_v(t)*forc_v(t)+ugust(t)*ugust(t))) @@ -748,7 +779,6 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & dthv(p) = dth(p)*(1._r8+0.61_r8*forc_q(t))+0.61_r8*forc_th(t)*dqh(p) zldis(p) = forc_hgt_u_patch(p) - displa(p) - ! Check to see if the forcing height is below the canopy height if (zldis(p) < 0._r8) then found = .true. @@ -766,457 +796,477 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & end if end if - do f = 1, fn - p = filterp(f) - c = veg_pp%column(p) - - ! Initialize Obukhov length scale and wind speed - - call MoninObukIni(ur(p), thv(c), dthv(p), zldis(p), z0mv(p), um(p), obu(p)) - num_iter(p) = 0._r8 - - end do - - ! Set counter for leaf temperature iteration (itlef) - - itlef = 1 - fnorig = fn - fporig(1:fn) = filterp(1:fn) - - ! Begin stability iteration event = 'can_iter' call t_start_lnd(event) - ITERATION : do while (itlef <= itmax .and. fn > 0) - - ! Determine friction velocity, and potential temperature and humidity - ! profiles of the surface boundary layer - call FrictionVelocity (begp, endp, fn, filterp, & - displa(begp:endp), z0mv(begp:endp), z0hv(begp:endp), z0qv(begp:endp), & - obu(begp:endp), itlef, ur(begp:endp), um(begp:endp), ugust_total(begp:endp), ustar(begp:endp), & - temp1(begp:endp), temp2(begp:endp), temp12m(begp:endp), temp22m(begp:endp), fm(begp:endp), & - frictionvel_vars) - - do f = 1, fn - p = filterp(f) - c = veg_pp%column(p) - t = veg_pp%topounit(p) - g = veg_pp%gridcell(p) - - tlbef(p) = t_veg(p) - del2(p) = del(p) - - ! Determine aerodynamic resistances - ram1(p) = 1._r8/(ustar(p)*ustar(p)/um(p)) - rah(p,above_canopy) = 1._r8/(temp1(p)*ustar(p)) - raw(p,above_canopy) = 1._r8/(temp2(p)*ustar(p)) - - ! Forbid removing more than 99% of wind speed in a time step. - ! This is mainly to avoid convergence issues since this is such a - ! basic form of iteration in this loop... - if (implicit_stress) then - tau(p) = forc_rho(t)*wind_speed_adj(p)/ram1(p) - call shr_flux_update_stress(wind_speed0(p), wsresp(t), tau_est(t), & - tau(p), prev_tau(p), tau_diff(p), prev_tau_diff(p), & - wind_speed_adj(p)) - ur(p) = max(1.0_r8, sqrt(wind_speed_adj(p)**2 + ugust(t)**2)) - end if - - ! Bulk boundary layer resistance of leaves - - uaf(p) = um(p)*sqrt( 1._r8/(ram1(p)*um(p)) ) - - ! Use pft parameter for leaf characteristic width - ! dleaf_patch if this is not an ed patch. - ! Otherwise, the value has already been loaded - ! during the FATES dynamics and/or initialization call - if(.not.veg_pp%is_fates(p)) then - dleaf_patch(p) = dleaf(veg_pp%itype(p)) - end if - - - cf = 0.01_r8/(sqrt(uaf(p))*sqrt( dleaf_patch(p) )) - rb(p) = 1._r8/(cf*uaf(p)) - rb1(p) = rb(p) - - ! Parameterization for variation of csoilc with canopy density from - ! X. Zeng, University of Arizona - - w = exp(-(elai(p)+esai(p))) - ! changed by K.Sakaguchi from here - ! transfer coefficient over bare soil is changed to a local variable - ! just for readability of the code (from line 680) - csoilb = (vkc/(0.13_r8*(z0mg(c)*uaf(p)/1.5e-5_r8)**0.45_r8)) - - !compute the stability parameter for ricsoilc ("S" in Sakaguchi&Zeng,2008) - - ri = ( grav*htop(p) * (taf(p) - t_grnd(c)) ) / (taf(p) * uaf(p) **2.00_r8) - - !! modify csoilc value (0.004) if the under-canopy is in stable condition - - if ( (taf(p) - t_grnd(c) ) > 0._r8) then - ! decrease the value of csoilc by dividing it with (1+gamma*min(S, 10.0)) - ! ria ("gmanna" in Sakaguchi&Zeng, 2008) is a constant (=0.5) - ricsoilc = csoilc / (1.00_r8 + ria*min( ri, 10.0_r8) ) - csoilcn = csoilb*w + ricsoilc*(1._r8-w) - else - csoilcn = csoilb*w + csoilc*(1._r8-w) - end if - - !! Sakaguchi changes for stability formulation ends here - - rah(p,below_canopy) = 1._r8/(csoilcn*uaf(p)) - raw(p,below_canopy) = rah(p,below_canopy) - if (use_lch4) then - grnd_ch4_cond(p) = 1._r8/(raw(p,above_canopy)+raw(p,below_canopy)) - end if - - ! Stomatal resistances for sunlit and shaded fractions of canopy. - ! Done each iteration to account for differences in eah, tv. - - svpts(p) = el(p) ! Pa - eah(p) = forc_pbot(t) * qaf(p) / mm_epsilon ! Pa - rhaf(p) = eah(p)/svpts(p) - - ! variables for history fields - rah_above(p) = rah(p,above_canopy) - raw_above(p) = raw(p,above_canopy) - rah_below(p) = rah(p,below_canopy) - raw_below(p) = raw(p,below_canopy) - vpd(p) = max((svpts(p) - eah(p)), vpd_min) * pa_to_kpa ! kPa - end do - - ! Modification for shrubs proposed by X.D.Z - ! Equivalent modification for soy following AgroIBIS - ! NOTE: the following block of code was moved out of Photosynthesis subroutine and - ! into here by M. Vertenstein on 4/6/2014 as part of making the photosynthesis - ! routine a separate module. This move was also suggested by S. Levis in the previous - ! version of the code. - ! BUG MV 4/7/2014 - is this the correct place to have it in the iteration? - ! THIS SHOULD BE MOVED OUT OF THE ITERATION but will change answers - - - do f = 1, fn - p = filterp(f) - c = veg_pp%column(p) - if(.not.veg_pp%is_fates(p)) then - ! soybean (crop with N fixation) - if (crop(veg_pp%itype(p)) >= 1 .and. nfixer(veg_pp%itype(p)) == 1) then - - btran(p) = min(1._r8, btran(p) * 1.25_r8) - end if - end if - end do - - - if ( use_fates ) then -#ifndef _OPENACC - call alm_fates%wrap_photosynthesis(bounds, fn, filterp(1:fn), & - svpts(begp:endp), eah(begp:endp), o2(begp:endp), & - co2(begp:endp), rb(begp:endp), dayl_factor(begp:endp), & - atm2lnd_vars, canopystate_vars, photosyns_vars) + totp(:) = 0 + + ! Recommended for this section: + ! export OMP_PROC_BIND=spread + ! export OMP_WAIT_POLICY=PASSIVE +#ifdef _OPENMP + t_start = omp_get_wtime() + w_diff(:) = 0. +#endif + !$OMP PARALLEL num_threads(fates_pproc) PRIVATE (it,np,f,p,c,t,g,itstoma,itlef,converge_stoma, & + !$OMP converge_tveg,cf,w,csoilb,ri,csoilcn, & + !$OMP ricsoilc,wta,wtl,wtshi,wtg0,wtga,rppdry, & + !$OMP efpot,rpp,wtaq,wtlq,snow_depth_c,fsno_dl, & + !$OMP elai_dl,rdl,wtsqi,wtgq0,wtgaq,dc1,dc2,efsh, & + !$OMP efeold,erre,lw_grnd,dels,ecidif,tstar, & + !$OMP qstar,thvstar,wc, iter_final,del_gs,reldel_gs,w_start,w_end, & + !$OMP delt,delt_snow,delt_soil,delt_h2osfc,delq_snow,delq_soil, & + !$OMP delq_h2osfc,slope_rad, e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT) if(fates_pproc>1) + + it = 0 + w_start = 0 + if(fates_pproc>1) then +#ifdef _OPENMP + it = omp_get_thread_num() + w_start = omp_get_wtime() #endif - else ! not use_fates - - if ( use_hydrstress ) then - call PhotosynthesisHydraulicStress (bounds, fn, filterp, & - svpts(begp:endp), eah(begp:endp), o2(begp:endp), co2(begp:endp), rb(begp:endp), bsun(begp:endp), & - bsha(begp:endp), btran(begp:endp), dayl_factor(begp:endp), & - qsatl(begp:endp), qaf(begp:endp), & - atm2lnd_vars, soilstate_vars, surfalb_vars, solarabs_vars, & - canopystate_vars, photosyns_vars) - else - call Photosynthesis (bounds, fn, filterp, & - svpts(begp:endp), eah(begp:endp), o2(begp:endp), co2(begp:endp), rb(begp:endp), btran(begp:endp), & - dayl_factor(begp:endp), atm2lnd_vars, surfalb_vars, solarabs_vars, & - canopystate_vars, photosyns_vars, 'sun') - end if - - if ( use_c13 ) then - call Fractionation (bounds, fn, filterp, & - cnstate_vars, solarabs_vars, surfalb_vars, photosyns_vars, & - 1) - endif - - do f = 1, fn - p = filterp(f) - c = veg_pp%column(p) - ! soybean (crop with N fixation) - if (crop(veg_pp%itype(p)) >= 1 .and. nfixer(veg_pp%itype(p)) == 1) then - btran(p) = min(1._r8, btran(p) * 1.25_r8) - end if - end do - - if ( .not. use_hydrstress ) then - call Photosynthesis (bounds, fn, filterp, & - svpts(begp:endp), eah(begp:endp), o2(begp:endp), co2(begp:endp), rb(begp:endp), btran(begp:endp), & - dayl_factor(begp:endp), atm2lnd_vars,surfalb_vars, solarabs_vars, & - canopystate_vars, photosyns_vars, 'sha') - end if - - - if ( use_c13 ) then - call Fractionation (bounds, fn, filterp, & - cnstate_vars, solarabs_vars, surfalb_vars, photosyns_vars, & - 0) - end if - - end if ! end of if use_fates - - do f = 1, fn - p = filterp(f) - c = veg_pp%column(p) - t = veg_pp%topounit(p) - g = veg_pp%gridcell(p) - - ! Sensible heat conductance for air, leaf and ground - ! Moved the original subroutine in-line... - - wta = 1._r8/rah(p,above_canopy) ! air - wtl = (elai(p)+esai(p))/rb(p) ! leaf - wtg(p) = 1._r8/rah(p,below_canopy) ! ground - wtshi = 1._r8/(wta+wtl+wtg(p)) - wtl0(p) = wtl*wtshi ! leaf - wtg0 = wtg(p)*wtshi ! ground - wta0(p) = wta*wtshi ! air - - wtga = wta0(p)+wtg0 ! ground + air - wtal(p) = wta0(p)+wtl0(p) ! air + leaf - - ! Fraction of potential evaporation from leaf - - if (fdry(p) > 0._r8) then - rppdry = fdry(p)*rb(p)*(laisun(p)/(rb(p)+rssun(p)) + & - laisha(p)/(rb(p)+rssha(p)))/elai(p) - else - rppdry = 0._r8 - end if - - ! Calculate canopy conductance for methane / oxygen (e.g. stomatal conductance & leaf bdy cond) - if (use_lch4) then - canopy_cond(p) = (laisun(p)/(rb(p)+rssun(p)) + laisha(p)/(rb(p)+rssha(p)))/max(elai(p), 0.01_r8) - end if - - efpot = forc_rho(t)*wtl*(qsatl(p)-qaf(p)) - ! When the hydraulic stress parameterization is active calculate rpp - ! but not transpiration - if ( use_hydrstress ) then - if (efpot > 0._r8) then - if (btran(p) > btran0) then - rpp = rppdry + fwet(p) - else - rpp = fwet(p) - end if - !Check total evapotranspiration from leaves - rpp = min(rpp, (qflx_tran_veg(p)+h2ocan(p)/dtime)/efpot) - else - rpp = 1._r8 - end if - else - - if (efpot > 0._r8) then - if (btran(p) > btran0) then - qflx_tran_veg(p) = efpot*rppdry - rpp = rppdry + fwet(p) - else - !No transpiration if btran below 1.e-10 - rpp = fwet(p) - qflx_tran_veg(p) = 0._r8 - end if - !Check total evapotranspiration from leaves - rpp = min(rpp, (qflx_tran_veg(p)+h2ocan(p)/dtime)/efpot) - else - !No transpiration if potential evaporation less than zero - rpp = 1._r8 - qflx_tran_veg(p) = 0._r8 - end if - end if - ! Update conductances for changes in rpp - ! Latent heat conductances for ground and leaf. - ! Air has same conductance for both sensible and latent heat. - ! Moved the original subroutine in-line... - - wtaq = frac_veg_nosno(p)/raw(p,above_canopy) ! air - wtlq = frac_veg_nosno(p)*(elai(p)+esai(p))/rb(p) * rpp ! leaf - - !Litter layer resistance. Added by K.Sakaguchi - snow_depth_c = z_dl ! critical depth for 100% litter burial by snow (=litter thickness) - fsno_dl = snow_depth(c)/snow_depth_c ! effective snow cover for (dry)plant litter - elai_dl = lai_dl*(1._r8 - min(fsno_dl,1._r8)) ! exposed (dry)litter area index - rdl = ( 1._r8 - exp(-elai_dl) ) / ( 0.004_r8*uaf(p)) ! dry litter layer resistance - - ! add litter resistance and Lee and Pielke 1992 beta - if (delq(p) < 0._r8) then !dew. Do not apply beta for negative flux (follow old rsoil) - wtgq(p) = frac_veg_nosno(p)/(raw(p,below_canopy)+rdl) - else - if (do_soilevap_beta()) then - wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,below_canopy)+rdl) - endif - end if - - wtsqi = 1._r8/(wtaq+wtlq+wtgq(p)) - - wtgq0 = wtgq(p)*wtsqi ! ground - wtlq0(p) = wtlq*wtsqi ! leaf - wtaq0(p) = wtaq*wtsqi ! air - - wtgaq = wtaq0(p)+wtgq0 ! air + ground - wtalq(p) = wtaq0(p)+wtlq0(p) ! air + leaf - - dc1 = forc_rho(t)*cpair*wtl - dc2 = hvap*forc_rho(t)*wtlq - - efsh = dc1*(wtga*t_veg(p)-wtg0*t_grnd(c)-wta0(p)*thm(p)) - efe(p) = dc2*(wtgaq*qsatl(p)-wtgq0*qg(c)-wtaq0(p)*forc_q(t)) - - ! Evaporation flux from foliage - - erre = 0._r8 - if (efe(p)*efeb(p) < 0._r8) then - efeold = efe(p) - efe(p) = 0.1_r8*efeold - erre = efe(p) - efeold - end if - ! fractionate ground emitted longwave - lw_grnd=(frac_sno(c)*t_soisno(c,snl(c)+1)**4 & - +(1._r8-frac_sno(c)-frac_h2osfc(c))*t_soisno(c,1)**4 & - +frac_h2osfc(c)*t_h2osfc(c)**4) - - dt_veg(p) = (sabv(p) + air(p) + bir(p)*t_veg(p)**4 + & - cir(p)*lw_grnd - efsh - efe(p)) / & - (- 4._r8*bir(p)*t_veg(p)**3 +dc1*wtga +dc2*wtgaq*qsatldT(p)) - t_veg(p) = tlbef(p) + dt_veg(p) - dels = dt_veg(p) - del(p) = abs(dels) - err(p) = 0._r8 - if (del(p) > delmax) then - dt_veg(p) = delmax*dels/del(p) - t_veg(p) = tlbef(p) + dt_veg(p) - err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + & - 4._r8*dt_veg(p)) + cir(p)*lw_grnd - & - (efsh + dc1*wtga*dt_veg(p)) - (efe(p) + & - dc2*wtgaq*qsatldT(p)*dt_veg(p)) - end if - - ! Fluxes from leaves to canopy space - ! "efe" was limited as its sign changes frequently. This limit may - ! result in an imbalance in "hvap*qflx_evap_veg" and - ! "efe + dc2*wtgaq*qsatdt_veg" - - efpot = forc_rho(t)*wtl*(wtgaq*(qsatl(p)+qsatldT(p)*dt_veg(p)) & - -wtgq0*qg(c)-wtaq0(p)*forc_q(t)) - qflx_evap_veg(p) = rpp*efpot - - ! Calculation of evaporative potentials (efpot) and - ! interception losses; flux in kg m**-2 s-1. ecidif - ! holds the excess energy if all intercepted water is evaporated - ! during the timestep. This energy is later added to the - ! sensible heat flux. - if ( use_hydrstress ) then - ecidif = max(0._r8,qflx_evap_veg(p)-qflx_tran_veg(p)-h2ocan(p)/dtime) - qflx_evap_veg(p) = min(qflx_evap_veg(p),qflx_tran_veg(p)+h2ocan(p)/dtime) - else - - ecidif = 0._r8 - if (efpot > 0._r8 .and. btran(p) > btran0) then - qflx_tran_veg(p) = efpot*rppdry - else - qflx_tran_veg(p) = 0._r8 - end if - ecidif = max(0._r8, qflx_evap_veg(p)-qflx_tran_veg(p)-h2ocan(p)/dtime) - qflx_evap_veg(p) = min(qflx_evap_veg(p),qflx_tran_veg(p)+h2ocan(p)/dtime) - end if - - ! The energy loss due to above two limits is added to - ! the sensible heat flux. - eflx_sh_veg(p) = efsh + dc1*wtga*dt_veg(p) + err(p) + erre + hvap*ecidif - - ! Re-calculate saturated vapor pressure, specific humidity, and their - ! derivatives at the leaf surface - - call QSat(t_veg(p), forc_pbot(t), el(p), deldT, qsatl(p), qsatldT(p)) - - ! Update vegetation/ground surface temperature, canopy air - ! temperature, canopy vapor pressure, aerodynamic temperature, and - ! Monin-Obukhov stability parameter for next iteration. - - taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p) - qaf(p) = wtlq0(p)*qsatl(p) + wtgq0*qg(c) + forc_q(t)*wtaq0(p) - - ! Update Obukhov length scale and wind speed including the - ! stability effect - - dth(p) = thm(p)-taf(p) - dqh(p) = forc_q(t)-qaf(p) - delq(p) = wtalq(p)*qg(c)-wtlq0(p)*qsatl(p)-wtaq0(p)*forc_q(t) + end if - tstar = temp1(p)*dth(p) - qstar = temp2(p)*dqh(p) + np = canopystate_vars%patch_par(it)%npatch + + patch_loop: do f = 1, np - thvstar = tstar*(1._r8+0.61_r8*forc_q(t)) + 0.61_r8*forc_th(t)*qstar - zeta(p) = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c)) + p = canopystate_vars%patch_par(it)%patch_list(f) + c = veg_pp%column(p) + t = veg_pp%topounit(p) + g = veg_pp%gridcell(p) + + if_filterp: if (in_filterp(p)) then + + totp(it) = totp(it) + 1 + + ! Initialize Obukhov length scale and wind speed + call MoninObukIni(ur(p), thv(c), dthv(p), zldis(p), z0mv(p), um(p), obu(p)) + + num_iter(p) = 0 + !rssun(p) is carried over from previous time-step + !rssha(p) is carried over from previous time-step + rssun_old(p) = -100._r8 + rssha_old(p) = -100._r8 + + ! Begin stability iterations + ! We evaluate temperature convergence inside + ! of stomatal conductance convergence. We nest these + ! loops to minimize stomatal conductance calculations + + itstoma = 0 + converge_stoma = .false. + iterate_stoma: do while(.not.converge_stoma) + + ! Set counter for leaf temperature iteration (itlef) + itlef = 1 + converge_tveg = .false. + iterate_tveg: do while(.not.converge_tveg) + + ! Determine friction velocity, and potential temperature and humidity + ! profiles of the surface boundary layer + call FrictionVelocity (begp, endp, 1, [p], & + displa(begp:endp), z0mv(begp:endp), z0hv(begp:endp), z0qv(begp:endp), & + obu(begp:endp), itlef, ur(begp:endp), um(begp:endp), ugust_total(begp:endp), ustar(begp:endp), & + temp1(begp:endp), temp2(begp:endp), temp12m(begp:endp), temp22m(begp:endp), fm(begp:endp), & + frictionvel_vars) + + tlbef(p) = t_veg(p) + del2(p) = del(p) + + ! Determine aerodynamic resistances + ram1(p) = 1._r8/(ustar(p)*ustar(p)/um(p)) + rah(p,above_canopy) = 1._r8/(temp1(p)*ustar(p)) + raw(p,above_canopy) = 1._r8/(temp2(p)*ustar(p)) + + ! Forbid removing more than 99% of wind speed in a time step. + ! This is mainly to avoid convergence issues since this is such a + ! basic form of iteration in this loop... + if (implicit_stress) then + tau(p) = forc_rho(t)*wind_speed_adj(p)/ram1(p) + call shr_flux_update_stress(wind_speed0(p), wsresp(t), tau_est(t), & + tau(p), prev_tau(p), tau_diff(p), prev_tau_diff(p), & + wind_speed_adj(p)) + ur(p) = max(1.0_r8, sqrt(wind_speed_adj(p)**2 + ugust(t)**2)) + end if + + ! Bulk boundary layer resistance of leaves + + uaf(p) = um(p)*sqrt( 1._r8/(ram1(p)*um(p)) ) + + ! Use pft parameter for leaf characteristic width + ! dleaf_patch if this is not an ed patch. + ! Otherwise, the value has already been loaded + ! during the FATES dynamics and/or initialization call + if(.not.veg_pp%is_fates(p)) then + dleaf_patch(p) = dleaf(veg_pp%itype(p)) + end if + + cf = 0.01_r8/(sqrt(uaf(p))*sqrt( dleaf_patch(p) )) + + rb(p) = 1._r8/(cf*uaf(p)) + rb1(p) = rb(p) + + ! Parameterization for variation of csoilc with canopy density from + ! X. Zeng, University of Arizona + + w = exp(-(elai(p)+esai(p))) + + ! changed by K.Sakaguchi from here + ! transfer coefficient over bare soil is changed to a local variable + ! just for readability of the code (from line 680) + csoilb = (vkc/(0.13_r8*(z0mg(c)*uaf(p)/1.5e-5_r8)**0.45_r8)) + + !compute the stability parameter for ricsoilc ("S" in Sakaguchi&Zeng,2008) + + ri = ( grav*htop(p) * (taf(p) - t_grnd(c)) ) / (taf(p) * uaf(p) **2.00_r8) + + !! modify csoilc value (0.004) if the under-canopy is in stable condition + + if ( (taf(p) - t_grnd(c) ) > 0._r8) then + ! decrease the value of csoilc by dividing it with (1+gamma*min(S, 10.0)) + ! ria ("gmanna" in Sakaguchi&Zeng, 2008) is a constant (=0.5) + ricsoilc = csoilc / (1.00_r8 + ria*min( ri, 10.0_r8) ) + csoilcn = csoilb*w + ricsoilc*(1._r8-w) + else + csoilcn = csoilb*w + csoilc*(1._r8-w) + end if + + !! Sakaguchi changes for stability formulation ends here + + rah(p,below_canopy) = 1._r8/(csoilcn*uaf(p)) + raw(p,below_canopy) = rah(p,below_canopy) + if (use_lch4) then + grnd_ch4_cond(p) = 1._r8/(raw(p,above_canopy)+raw(p,below_canopy)) + end if + + ! Stomatal resistances for sunlit and shaded fractions of canopy. + ! Done each iteration to account for differences in eah, tv. + + svpts(p) = el(p) ! Pa + eah(p) = forc_pbot(t) * qaf(p) / mm_epsilon ! Pa + rhaf(p) = eah(p)/svpts(p) + + ! variables for history fields + rah_above(p) = rah(p,above_canopy) + raw_above(p) = raw(p,above_canopy) + rah_below(p) = rah(p,below_canopy) + raw_below(p) = raw(p,below_canopy) + vpd(p) = max((svpts(p) - eah(p)), vpd_min) * pa_to_kpa ! kPa + + ! Modification for shrubs proposed by X.D.Z + ! Equivalent modification for soy following AgroIBIS + ! NOTE: the following block of code was moved out of Photosynthesis subroutine and + ! into here by M. Vertenstein on 4/6/2014 as part of making the photosynthesis + ! routine a separate module. This move was also suggested by S. Levis in the previous + ! version of the code. + ! BUG MV 4/7/2014 - is this the correct place to have it in the iteration? + ! THIS SHOULD BE MOVED OUT OF THE ITERATION but will change answers - + + if(.not.veg_pp%is_fates(p)) then + ! soybean (crop with N fixation) + if (crop(veg_pp%itype(p)) >= 1 .and. nfixer(veg_pp%itype(p)) == 1) then + btran(p) = min(1._r8, btran(p) * 1.25_r8) + end if + end if + + if(do_b4b)then + call WrapPhotosynthesis(bounds,p,svpts,eah,o2,co2,rb,dayl_factor, & + btran,qsatl,qaf,atm2lnd_vars,canopystate_vars,photosyns_vars, & + soilstate_vars, surfalb_vars,solarabs_vars,cnstate_vars,energyflux_vars) + end if + + ! Sensible heat conductance for air, leaf and ground + ! Moved the original subroutine in-line... + + wta = 1._r8/rah(p,above_canopy) ! air + wtl = (elai(p)+esai(p))/rb(p) ! leaf + wtg(p) = 1._r8/rah(p,below_canopy) ! ground + wtshi = 1._r8/(wta+wtl+wtg(p)) + wtl0(p) = wtl*wtshi ! leaf + wtg0 = wtg(p)*wtshi ! ground + wta0(p) = wta*wtshi ! air + + wtga = wta0(p)+wtg0 ! ground + air + wtal(p) = wta0(p)+wtl0(p) ! air + leaf + + ! Fraction of potential evaporation from leaf + + if (fdry(p) > 0._r8) then + rppdry = fdry(p)*rb(p)*(laisun(p)/(rb(p)+rssun(p)) + & + laisha(p)/(rb(p)+rssha(p)))/elai(p) + else + rppdry = 0._r8 + end if + + ! Calculate canopy conductance for methane / oxygen + ! (e.g. stomatal conductance & leaf bdy cond) + if (use_lch4) then + canopy_cond(p) = (laisun(p)/(rb(p)+rssun(p)) + laisha(p) / & + (rb(p)+rssha(p)))/max(elai(p), 0.01_r8) + end if + + efpot = forc_rho(t)*wtl*(qsatl(p)-qaf(p)) + ! When the hydraulic stress parameterization is active calculate rpp + ! but not transpiration + if ( use_hydrstress ) then + if (efpot > 0._r8) then + if (btran(p) > btran0) then + rpp = rppdry + fwet(p) + else + rpp = fwet(p) + end if + !Check total evapotranspiration from leaves + rpp = min(rpp, (qflx_tran_veg(p)+h2ocan(p)/dtime)/efpot) + else + rpp = 1._r8 + end if + else + if (efpot > 0._r8) then + if (btran(p) > btran0) then + qflx_tran_veg(p) = efpot*rppdry + rpp = rppdry + fwet(p) + else + !No transpiration if btran below 1.e-10 + rpp = fwet(p) + qflx_tran_veg(p) = 0._r8 + end if + !Check total evapotranspiration from leaves + rpp = min(rpp, (qflx_tran_veg(p)+h2ocan(p)/dtime)/efpot) + else + !No transpiration if potential evaporation less than zero + rpp = 1._r8 + qflx_tran_veg(p) = 0._r8 + end if + end if + + ! Update conductances for changes in rpp + ! Latent heat conductances for ground and leaf. + ! Air has same conductance for both sensible and latent heat. + ! Moved the original subroutine in-line... + + wtaq = frac_veg_nosno(p)/raw(p,above_canopy) ! air + wtlq = frac_veg_nosno(p)*(elai(p)+esai(p))/rb(p) * rpp ! leaf + + !Litter layer resistance. Added by K.Sakaguchi + snow_depth_c = z_dl ! critical depth for 100% litter burial by snow (=litter thickness) + fsno_dl = snow_depth(c)/snow_depth_c ! effective snow cover for (dry)plant litter + elai_dl = lai_dl*(1._r8 - min(fsno_dl,1._r8)) ! exposed (dry)litter area index + rdl = ( 1._r8 - exp(-elai_dl) ) / ( 0.004_r8*uaf(p)) ! dry litter layer resistance + + ! add litter resistance and Lee and Pielke 1992 beta + if (delq(p) < 0._r8) then !dew. Do not apply beta for negative flux (follow old rsoil) + wtgq(p) = frac_veg_nosno(p)/(raw(p,below_canopy)+rdl) + else + if (do_soilevap_beta()) then + wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,below_canopy)+rdl) + endif + end if + + wtsqi = 1._r8/(wtaq+wtlq+wtgq(p)) + + wtgq0 = wtgq(p)*wtsqi ! ground + wtlq0(p) = wtlq*wtsqi ! leaf + wtaq0(p) = wtaq*wtsqi ! air + + wtgaq = wtaq0(p)+wtgq0 ! air + ground + wtalq(p) = wtaq0(p)+wtlq0(p) ! air + leaf + + dc1 = forc_rho(t)*cpair*wtl + dc2 = hvap*forc_rho(t)*wtlq + + efsh = dc1*(wtga*t_veg(p)-wtg0*t_grnd(c)-wta0(p)*thm(p)) + efe(p) = dc2*(wtgaq*qsatl(p)-wtgq0*qg(c)-wtaq0(p)*forc_q(t)) + + ! Evaporation flux from foliage + erre = 0._r8 + if (efe(p)*efeb(p) < 0._r8) then + efeold = efe(p) + efe(p) = 0.1_r8*efeold + erre = efe(p) - efeold + end if + + ! fractionate ground emitted longwave + lw_grnd=(frac_sno(c)*t_soisno(c,snl(c)+1)**4 & + +(1._r8-frac_sno(c)-frac_h2osfc(c))*t_soisno(c,1)**4 & + +frac_h2osfc(c)*t_h2osfc(c)**4) + + dt_veg(p) = (sabv(p) + air(p) + bir(p)*t_veg(p)**4 + & + cir(p)*lw_grnd - efsh - efe(p)) / & + (- 4._r8*bir(p)*t_veg(p)**3 +dc1*wtga +dc2*wtgaq*qsatldT(p)) + t_veg(p) = tlbef(p) + dt_veg(p) + dels = dt_veg(p) + del(p) = abs(dels) + err(p) = 0._r8 + if (del(p) > delmax) then + dt_veg(p) = delmax*dels/del(p) + t_veg(p) = tlbef(p) + dt_veg(p) + err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + & + 4._r8*dt_veg(p)) + cir(p)*lw_grnd - & + (efsh + dc1*wtga*dt_veg(p)) - (efe(p) + & + dc2*wtgaq*qsatldT(p)*dt_veg(p)) + end if + + ! Fluxes from leaves to canopy space + ! "efe" was limited as its sign changes frequently. This limit may + ! result in an imbalance in "hvap*qflx_evap_veg" and + ! "efe + dc2*wtgaq*qsatdt_veg" + + efpot = forc_rho(t)*wtl*(wtgaq*(qsatl(p)+qsatldT(p)*dt_veg(p)) & + -wtgq0*qg(c)-wtaq0(p)*forc_q(t)) + qflx_evap_veg(p) = rpp*efpot + + ! Calculation of evaporative potentials (efpot) and + ! interception losses; flux in kg m**-2 s-1. ecidif + ! holds the excess energy if all intercepted water is evaporated + ! during the timestep. This energy is later added to the + ! sensible heat flux. + if ( use_hydrstress ) then + ecidif = max(0._r8,qflx_evap_veg(p)-qflx_tran_veg(p)-h2ocan(p)/dtime) + qflx_evap_veg(p) = min(qflx_evap_veg(p),qflx_tran_veg(p)+h2ocan(p)/dtime) + else + + ecidif = 0._r8 + if (efpot > 0._r8 .and. btran(p) > btran0) then + qflx_tran_veg(p) = efpot*rppdry + else + qflx_tran_veg(p) = 0._r8 + end if + ecidif = max(0._r8, qflx_evap_veg(p)-qflx_tran_veg(p)-h2ocan(p)/dtime) + qflx_evap_veg(p) = min(qflx_evap_veg(p),qflx_tran_veg(p)+h2ocan(p)/dtime) + end if + + ! The energy loss due to above two limits is added to + ! the sensible heat flux. + eflx_sh_veg(p) = efsh + dc1*wtga*dt_veg(p) + err(p) + erre + hvap*ecidif + + ! Re-calculate saturated vapor pressure, specific humidity, and their + ! derivatives at the leaf surface + + call QSat(t_veg(p), forc_pbot(t), el(p), deldT, qsatl(p), qsatldT(p)) + + ! Update vegetation/ground surface temperature, canopy air + ! temperature, canopy vapor pressure, aerodynamic temperature, and + ! Monin-Obukhov stability parameter for next iteration. + + taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p) + qaf(p) = wtlq0(p)*qsatl(p) + wtgq0*qg(c) + forc_q(t)*wtaq0(p) + + ! Update Obukhov length scale and wind speed including the + ! stability effect + + dth(p) = thm(p)-taf(p) + dqh(p) = forc_q(t)-qaf(p) + delq(p) = wtalq(p)*qg(c)-wtlq0(p)*qsatl(p)-wtaq0(p)*forc_q(t) + + tstar = temp1(p)*dth(p) + qstar = temp2(p)*dqh(p) + + thvstar = tstar*(1._r8+0.61_r8*forc_q(t)) + 0.61_r8*forc_th(t)*qstar + zeta(p) = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c)) + + if (zeta(p) >= 0._r8) then !stable + zeta(p) = min(2._r8,max(zeta(p),0.01_r8)) + um(p) = max(ur(p),0.1_r8) + else !unstable + zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8)) + if ((.not. atm_gustiness) .or. force_land_gustiness) then + wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8 + ugust_total(p) = sqrt(ugust(t)**2 + wc**2) + um(p) = sqrt(ur(p)*ur(p)+wc*wc) + else + um(p) = max(ur(p),0.1_r8) + end if + end if + obu(p) = zldis(p)/zeta(p) + + if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1 + if (nmozsgn(p) >= 4) obu(p) = zldis(p)/(-0.01_r8) + obuold(p) = obu(p) + + ! laminar boundary resistance for h2o over leaf, + ! should I make this consistent for latent heat calculation? + lbl_rsc_h2o(p) = getlblcef(forc_rho(t),t_veg(p))*uaf(p)/(uaf(p)**2._r8+1.e-10_r8) + + ! Test for convergence + iter_final = itlef + itlef = itlef+1 + + if(itlef>itmin)then + dele(p) = abs(efe(p)-efeb(p)) + efeb(p) = efe(p) + det(p) = max(del(p),del2(p)) + num_iter(p) = real(itlef,r8) + + if ( (.not. (det(p) < dtmin .and. dele(p) < dlemin) .or. & + (implicit_stress .and. abs(tau_diff(p)) >= dtaumin)) .and. & + (itlef <= itmax)) then + converge_tveg = .false. + else + converge_tveg = .true. + end if + end if + + end do iterate_tveg + + ! Evaluate quality of conductance solution + ! + ! Criteria for finding a solution to the outer loop + ! + ! 1) Always make sure that at least one photosynthesis call + ! is made. (ie itstoma>0) + ! 2) Calculate the change in resistance that was made on the + ! last solution. If the difference is negligable, and + ! condition 1 is satisfied, then you have a solution + ! 3) Exit if too many attempts and accept what you have + ! (ie. itstoma>itmax_stomata + + !reldel_rs = 2._r8*max( abs(rssun(p)-rssun_old(p))/(rssun(p)+rssun_old(p)), & + ! abs(rssha(p)-rssha_old(p))/(rssha(p)+rssha_old(p)) ) + + del_gs = max( abs(1._r8/rssun(p)-1._r8/rssun_old(p)), & + abs(1._r8/rssha(p)-1._r8/rssha_old(p)) ) + + ! Let's use the harmonic mean of the conductances + ! which is the inverse of the sum of resistances + hmean_gs = 2._r8/(rssun(p)+rssha(p)) + reldel_gs = del_gs / hmean_gs + + ! (x/1)/(1/y) + ! y + reldel_gs = max( rssun(p)*abs(1._r8/rssun(p)-1._r8/rssun_old(p)), & + rssha(p)*abs(1._r8/rssha(p)-1._r8/rssha_old(p))) + + istoma_converge_if: if( do_b4b .or. & + !(del_gs < max_del_gs ) .or. & + (reldel_gs < max_reldel_gs) .or. & + (itstoma>=itmax_stomata) ) then + converge_stoma = .true. - if (zeta(p) >= 0._r8) then !stable - zeta(p) = min(2._r8,max(zeta(p),0.01_r8)) - um(p) = max(ur(p),0.1_r8) - else !unstable - zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8)) - if ((.not. atm_gustiness) .or. force_land_gustiness) then - wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8 - ugust_total(p) = sqrt(ugust(t)**2 + wc**2) - um(p) = sqrt(ur(p)*ur(p)+wc*wc) else - um(p) = max(ur(p),0.1_r8) - end if - end if - obu(p) = zldis(p)/zeta(p) - if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1 - if (nmozsgn(p) >= 4) obu(p) = zldis(p)/(-0.01_r8) - obuold(p) = obu(p) + ! Update the outer (stomata c) counter + itstoma = itstoma + 1 + num_iter(p) = num_iter(p) + 1 - end do ! end of filtered pft loop + ! Update the previous resistances + rssun_old(p) = rssun(p) + rssha_old(p) = rssha(p) - do f = 1, fn - p = filterp(f) - t = veg_pp%topounit(p) - lbl_rsc_h2o(p) = getlblcef(forc_rho(t),t_veg(p))*uaf(p)/(uaf(p)**2._r8+1.e-10_r8) !laminar boundary resistance for h2o over leaf, should I make this consistent for latent heat calculation? - enddo - - ! Test for convergence - iter_final = itlef - itlef = itlef+1 - if (itlef > itmin) then - do f = 1, fn - p = filterp(f) - dele(p) = abs(efe(p)-efeb(p)) - efeb(p) = efe(p) - det(p) = max(del(p),del2(p)) - num_iter(p) = real(itlef,r8) - end do - fnold = fn - fn = 0 - do f = 1, fnold - p = filterp(f) - if (.not. (det(p) < dtmin .and. dele(p) < dlemin) .or. & - (implicit_stress .and. abs(tau_diff(p)) >= dtaumin)) then - fn = fn + 1 - filterp(fn) = p - end if - end do - end if + ! Call photosynthesis and retrieve + ! updated stomatal conductances - end do ITERATION ! End stability iteration - call t_stop_lnd(event) + ! Instead of updating stomatal conductances + ! we hold the value calculated in the outer loop + ! as constant during this inner loop (tveg) iteration - fn = fnorig - filterp(1:fn) = fporig(1:fn) - - do f = 1, fn - p = filterp(f) - c = veg_pp%column(p) - t = veg_pp%topounit(p) - g = veg_pp%gridcell(p) + call WrapPhotosynthesis(bounds,p,svpts,eah,o2,co2,rb,dayl_factor, & + btran,qsatl,qaf,atm2lnd_vars,canopystate_vars,photosyns_vars, & + soilstate_vars, surfalb_vars,solarabs_vars,cnstate_vars,energyflux_vars) + end if istoma_converge_if + end do iterate_stoma + ! Energy balance check in canopy lw_grnd=(frac_sno(c)*t_soisno(c,snl(c)+1)**4 & @@ -1314,44 +1364,317 @@ subroutine CanopyFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, & ' wind_speed_adj= ',wind_speed_adj(p),' iter_final= ',iter_final end if end if - - end do - - if ( use_fates ) then - + if ( use_fates ) then + #ifndef _OPENACC - call alm_fates%wrap_accumulatefluxes(bounds,fn,filterp(1:fn)) - call alm_fates%wrap_hydraulics_drive(bounds,fn,filterp(1:fn),soilstate_vars, & - solarabs_vars,energyflux_vars) + ! Calculate non-leaf respiration terms + ! and accumulate fluxes + call alm_fates%WrapAccumulateFluxes(bounds,p) #endif - else - - ! Determine total photosynthesis - call PhotosynthesisTotal(fn, filterp, & - atm2lnd_vars, cnstate_vars, canopystate_vars, photosyns_vars) - - ! Filter out patches which have small energy balance errors; report others - - fnold = fn - fn = 0 - do f = 1, fnold - p = filterp(f) - if (abs(err(p)) > 0.1_r8) then - fn = fn + 1 - filterp(fn) = p - end if - end do + end if - do f = 1, fn - p = filterp(f) - write(iulog,*) 'energy balance in canopy ',p,', err=',err(p) - end do + end if if_filterp + + end do patch_loop +#ifdef _OPENMP + w_end = omp_get_wtime() + w_diff(it) = (w_end-w_start) +#endif + + !$OMP END PARALLEL + +#ifdef _OPENMP + t_end = omp_get_wtime() + setup_overhead = maxval(w_diff)/(t_end - t_start) + do f=0,fates_pproc-1 + npatch(f) = canopystate_vars%patch_par(f)%npatch + end do + !write(iulog,'(A)') '---------------------' + ! write(iulog,'(A,2x,F7.5,1x,A)') "Group time: ", t_end - t_start, "s" + !write(iulog,'(A,2x,F7.5)') "Setup Efficiency: ", setup_overhead + !write(iulog,'(A,*(2x,I3))') "Worker patch: ",npatch(:) + !write(iulog,'(A,*(2x,F7.5))') "Balance Efficiency: ",(minval(w_diff(:))+0.00001)/maxval(w_diff(:)) +#endif + + call t_stop_lnd(event) + + if(sum(totp).ne.fn)then + write(iulog,*) 'totp should equal fn...',sum(totp),fn + write(iulog,*) 'filterp:',filterp(1:fn) + do it=0,fates_pproc-1 + np = canopystate_vars%patch_par(it)%npatch + write(iulog,*) 'patch_list: ',canopystate_vars%patch_par(it)%patch_list(1:np) + write(iulog,*) 'totp(it): ',totp(it) + end do + stop + end if + + if ( use_fates ) then +#ifndef _OPENACC + call alm_fates%wrap_hydraulics_drive(bounds,fn,filterp(1:fn),soilstate_vars, & + solarabs_vars,energyflux_vars) +#endif + else + ! Determine total photosynthesis + call PhotosynthesisTotal(fn, filterp, & + atm2lnd_vars, cnstate_vars, canopystate_vars, photosyns_vars) + end if + + ! Report high energy balance errors + do f = 1, fn + p = filterp(f) + if (abs(err(p)) > 0.1_r8) then + write(iulog,*) 'energy balance in canopy ',p,', err=',err(p) end if + end do + + end associate +end subroutine CanopyFluxes + + ! ========================================================================= + + subroutine WrapPhotosynthesis(bounds,p,svpts,eah,o2,co2,rb,dayl_factor, & + btran,qsatl,qaf,atm2lnd_vars,canopystate_vars,photosyns_vars, & + soilstate_vars, surfalb_vars,solarabs_vars,cnstate_vars,energyflux_vars) + + + type(bounds_type) , intent(in) :: bounds + integer,intent(in) :: p ! patch index from begp:endp + real(r8),intent(in) :: svpts(bounds%begp:bounds%endp) + real(r8),intent(in) :: eah(bounds%begp:bounds%endp) + real(r8),intent(in) :: o2(bounds%begp:bounds%endp) + real(r8),intent(in) :: co2(bounds%begp:bounds%endp) + real(r8),intent(in) :: rb(bounds%begp:bounds%endp) + real(r8),intent(in) :: dayl_factor(bounds%begp:bounds%endp) + real(r8),intent(inout) :: btran(bounds%begp:bounds%endp) + real(r8),intent(in) :: qsatl(bounds%begp:bounds%endp) + real(r8),intent(in) :: qaf(bounds%begp:bounds%endp) + + type(atm2lnd_type) , intent(inout) :: atm2lnd_vars + type(canopystate_type) , intent(inout) :: canopystate_vars + type(photosyns_type) , intent(inout) :: photosyns_vars + type(soilstate_type) , intent(inout) :: soilstate_vars + type(surfalb_type) , intent(inout) :: surfalb_vars + type(solarabs_type) , intent(inout) :: solarabs_vars + type(cnstate_type) , intent(inout) :: cnstate_vars + type(energyflux_type) , intent(inout) :: energyflux_vars - end associate + integer :: begp,endp + begp = bounds%begp + endp = bounds%endp + + ! Instead of updating stomatal conductances + ! we hold the value calculated in the outer loop + ! as constant during this inner loop (tveg) iteration + if_fates: if ( use_fates ) then +#ifndef _OPENACC - end subroutine CanopyFluxes + call alm_fates%WrapPatchPhotosynthesis(bounds, p, & + svpts(begp:endp), eah(begp:endp), o2(begp:endp), & + co2(begp:endp), rb(begp:endp), dayl_factor(begp:endp), & + atm2lnd_vars, canopystate_vars, photosyns_vars) +#endif + else ! not use_fates + if ( use_hydrstress ) then + call PhotosynthesisHydraulicStress (bounds, 1, [p], & + svpts(begp:endp), eah(begp:endp), o2(begp:endp), co2(begp:endp), rb(begp:endp), & + energyflux_vars%bsun_patch(begp:endp), energyflux_vars%bsha_patch(begp:endp), & + btran(begp:endp), dayl_factor(begp:endp), & + qsatl(begp:endp), qaf(begp:endp), & + atm2lnd_vars, soilstate_vars, surfalb_vars, solarabs_vars, & + canopystate_vars, photosyns_vars) + else + call Photosynthesis (bounds, 1, [p], & + svpts(begp:endp), eah(begp:endp), o2(begp:endp), co2(begp:endp), rb(begp:endp), btran(begp:endp), & + dayl_factor(begp:endp), atm2lnd_vars, surfalb_vars, solarabs_vars, & + canopystate_vars, photosyns_vars, 'sun') + end if + + if ( use_c13 ) then + call Fractionation (bounds, 1, [p], & + cnstate_vars, solarabs_vars, surfalb_vars, photosyns_vars, 1) + endif + + ! soybean (crop with N fixation) + if (crop(veg_pp%itype(p)) >= 1 .and. nfixer(veg_pp%itype(p)) == 1) then + btran(p) = min(1._r8, btran(p) * 1.25_r8) + end if + + if ( .not. use_hydrstress ) then + call Photosynthesis (bounds, 1, [p], & + svpts(begp:endp), eah(begp:endp), o2(begp:endp), co2(begp:endp), rb(begp:endp), btran(begp:endp), & + dayl_factor(begp:endp), atm2lnd_vars,surfalb_vars, solarabs_vars, & + canopystate_vars, photosyns_vars, 'sha') + end if + + if ( use_c13 ) then + call Fractionation (bounds, 1, [p], & + cnstate_vars, solarabs_vars, surfalb_vars, photosyns_vars, 0) + end if + end if if_fates ! end of if use_fates + + end subroutine WrapPhotosynthesis + + ! ========================================================================================================= + + pure subroutine sort_tasks(n, load_size, indices) + + integer, intent(in) :: n + real(8), intent(inout) :: load_size(n) + integer, intent(inout) :: indices(n) + + integer :: i, j, temp_idx + real(8) :: temp_load + + ! Insertion sort is extremely efficient for N < 30 + do i = 2, n + temp_load = load_size(i) + temp_idx = indices(i) + j = i - 1 + + ! Move elements that are smaller than temp_load to one position ahead + ! to achieve descending order + do while (j >= 1) + if (load_size(j) < temp_load) then + load_size(j+1) = load_size(j) + indices(j+1) = indices(j) + j = j - 1 + else + exit + end if + end do + + load_size(j+1) = temp_load + indices(j+1) = temp_idx + end do + + end subroutine sort_tasks + + !call PatchLoadBalance(bounds_clump, filter(nc)%num_nolakeurbanp, & + ! filter(nc)%nolakeurbanp,canopystate_vars, & + ! alm_fates%fates(nc)%bc_out,alm_fates%f2hmap(nc)%hsites,.false.) + + subroutine PatchLoadBalance(bounds, num_nolakeurbanp, filter_nolakeurbanp, canopystate_vars, fates_bc_out, fates_hsites, is_coldstart) + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_nolakeurbanp ! number of column non-lake, non-urban points in pft filter + integer , intent(in) :: filter_nolakeurbanp(:) ! patch filter for non-lake, non-urban points + type(canopystate_type) , intent(inout) :: canopystate_vars + type(bc_out_type) , intent(in) :: fates_bc_out(:) + integer , intent(in) :: fates_hsites(bounds%begc:) + logical , intent(in) :: is_coldstart + + integer :: filterp(bounds%endp-bounds%begp+1) + real(r8) :: load_size(bounds%endp-bounds%begp+1) + real(r8) :: task_load(0:fates_pproc-1) ! this works because M <= N + integer :: f,fn,ifp,c,p,s + integer :: npt ! number of patches per thread + integer :: ct ! Current thread + integer :: np ! number of patches on this proc + integer :: ctask ! thread with lowest task load + real(r8), parameter :: bl_default = 1._r8 ! The load size for a big-leaf canopy + ! This is roughly how many leaf layers + + ! If there is no patch-parallism + ! or there is no fates, then we have nothing to update here + + np = bounds%endp-bounds%begp+1 + + filterp(:)=0 + fn = 0 + if(is_coldstart)then + do p = bounds%begp,bounds%endp + c = veg_pp%column(p) + if(col_pp%is_soil(c) .and. col_pp%active(c))then + fn = fn + 1 + filterp(fn) = p + load_size(fn) = 1 + end if + end do + else + do p = bounds%begp,bounds%endp + c = veg_pp%column(p) + + if(col_pp%is_soil(c) .and. col_pp%active(c))then + if(.not.veg_pp%is_fates(p))then + fn = fn + 1 + filterp(fn) = p + load_size(fn) = bl_default + else + ! If this patch is under FATES jurisdiction + ! we only include it in the list if this + ! is actually an active FATES patch + ! is_fates determines jurisdiction, is_veg + ! determines active or not active + if(veg_pp%is_veg(p)) then + s = fates_hsites(c) + fn = fn + 1 + filterp(fn) = p + ifp = p - col_pp%pfti(c) + load_size(fn) = fates_bc_out(s)%load_size(ifp) + end if + end if + end if + end do + end if + + !if(fn==0)then + ! write(iulog,*)'fn:',fn + ! write(iulog,*)'is_veg:',veg_pp%is_veg(bounds%begp:bounds%endp) + ! stop + !end if + + !!fn = 0 + !!do f = 1,np + !! p = filter_nolakeurbanp(f) + !! c = veg_pp%column(p) + !! s = fates_hsites(c) + !! fn = fn + 1 + !! filterp(fn) = p + ! This is a predictive filter, we don't load + ! balance every day, so we are using the + ! sno fraction from once a day to guess + ! the sno fraction over the time-steps. + ! The actual snow fraction will indeed be used + ! to filter out patches though + !! if (canopystate_vars%frac_veg_nosno_patch(p) == 0) then + !! load_size(fn) = 0.25*bl_default + !! else + !! if(veg_pp%is_fates(p))then + !! ifp = p - col_pp%pfti(c) + !! load_size(fn) = fates_bc_out(s)%load_size(ifp) + !! else + !! load_size(fn) = bl_default + !! end if + !! end if + !!end do + + do ct = 0,fates_pproc-1 + canopystate_vars%patch_par(ct)%npatch = 0 + canopystate_vars%patch_par(ct)%patch_list(:) = 0 + end do + + if(fates_pproc==1 .or. .not.use_fates)then + do f = 1,fn + canopystate_vars%patch_par(0)%patch_list(f) = filterp(f) + end do + canopystate_vars%patch_par(0)%npatch = fn + else + ! Sorts the load sizes and their indices in decending order + call sort_tasks(fn, load_size, filterp) + task_load(:) = 0._r8 + do f = 1,fn + ct = minloc(task_load,dim=1)-1 + task_load(ct) = task_load(ct) + load_size(f) + canopystate_vars%patch_par(ct)%npatch = canopystate_vars%patch_par(ct)%npatch + 1 + npt = canopystate_vars%patch_par(ct)%npatch + canopystate_vars%patch_par(ct)%patch_list(npt) = filterp(f) + end do + end if + + end subroutine PatchLoadBalance + + end module CanopyFluxesMod diff --git a/components/elm/src/biogeophys/CanopyStateType.F90 b/components/elm/src/biogeophys/CanopyStateType.F90 index b07db4b7ad04..fdfc4088ed0b 100644 --- a/components/elm/src/biogeophys/CanopyStateType.F90 +++ b/components/elm/src/biogeophys/CanopyStateType.F90 @@ -5,7 +5,7 @@ module CanopyStateType use shr_infnan_mod , only : shr_infnan_isnan,nan => shr_infnan_nan, assignment(=) use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun - use decompMod , only : bounds_type + use decompMod , only : bounds_type,fates_pproc use landunit_varcon , only : istsoil, istcrop use elm_varcon , only : spval,ispval use elm_varpar , only : nlevcan, nvegwcs @@ -31,6 +31,13 @@ module CanopyStateType !------------------------------------------------------------------------------ !$acc declare create(perchroot ) !$acc declare create(perchroot_alt) + + ! + type :: PatchPar_type + integer :: npatch + integer,allocatable :: patch_list(:) + end type PatchPar_type + ! ! !PUBLIC TYPES: type, public :: CanopyState_type @@ -76,6 +83,9 @@ module CanopyStateType ! for non-ED/FATES this is the same as pftcon%dleaf() real(r8), pointer :: lbl_rsc_h2o_patch (:) ! laminar boundary layer resistance for water over dry leaf (s/m) real(r8) , pointer :: vegwp_patch (:,:) ! patch vegetation water matric potential (mm) + + type(patchpar_type), pointer :: patch_par(:) ! Patch index data for patch-level parallism + contains procedure, public :: Init @@ -114,6 +124,7 @@ subroutine InitAllocate(this, bounds ) integer :: begp, endp integer :: begc, endc integer :: begg, endg + integer :: m,mt !------------------------------------------------------------------------ @@ -160,7 +171,13 @@ subroutine InitAllocate(this, bounds ) allocate(this%lbl_rsc_h2o_patch (begp:endp)) ; this%lbl_rsc_h2o_patch (:) = spval allocate(this%vegwp_patch (begp:endp,1:nvegwcs)) ; this%vegwp_patch (:,:) = spval - + allocate(this%patch_par(0:fates_pproc-1)) + do m = 0,fates_pproc-1 + allocate(this%patch_par(m)%patch_list(1:endp-begp+1)) + this%patch_par(m)%patch_list(:) = -1 + this%patch_par(m)%npatch = 0 + end do + end subroutine InitAllocate !----------------------------------------------------------------------- diff --git a/components/elm/src/biogeophys/FrictionVelocityMod.F90 b/components/elm/src/biogeophys/FrictionVelocityMod.F90 index 9d2883ffc81c..856ede6ec647 100644 --- a/components/elm/src/biogeophys/FrictionVelocityMod.F90 +++ b/components/elm/src/biogeophys/FrictionVelocityMod.F90 @@ -16,7 +16,6 @@ module FrictionVelocityMod ! ! !PUBLIC TYPES: implicit none - save logical, public :: implicit_stress = .false. logical, public :: atm_gustiness = .false. @@ -35,7 +34,7 @@ module FrictionVelocityMod contains !------------------------------------------------------------------------------ - subroutine FrictionVelocity(lbn, ubn, fn, filtern, & + pure subroutine FrictionVelocity(lbn, ubn, fn, filtern, & displa, z0m, z0h, z0q, & obu, iter, ur, um, ugust, ustar, & temp1, temp2, temp12m, temp22m, fm,frictionvel_vars,landunit_index) @@ -393,10 +392,10 @@ subroutine FrictionVelocity(lbn, ubn, fn, filtern, & end associate - end subroutine FrictionVelocity + end subroutine FrictionVelocity !------------------------------------------------------------------------------ - real(r8) function StabilityFunc1(zeta) + pure real(r8) function StabilityFunc1(zeta) !$acc routine seq ! !DESCRIPTION: ! Stability function for rib < 0. @@ -420,7 +419,7 @@ real(r8) function StabilityFunc1(zeta) end function StabilityFunc1 !------------------------------------------------------------------------------ - real(r8) function StabilityFunc2(zeta) + pure real(r8) function StabilityFunc2(zeta) !$acc routine seq ! !DESCRIPTION: ! Stability function for rib < 0. @@ -442,7 +441,7 @@ real(r8) function StabilityFunc2(zeta) end function StabilityFunc2 !----------------------------------------------------------------------- - subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu) + pure subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu) !$acc routine seq ! !DESCRIPTION: ! Initialization of the Obukhov length scale. diff --git a/components/elm/src/biogeophys/PhotosynthesisType.F90 b/components/elm/src/biogeophys/PhotosynthesisType.F90 index e00e725ef96b..0be8022a766e 100644 --- a/components/elm/src/biogeophys/PhotosynthesisType.F90 +++ b/components/elm/src/biogeophys/PhotosynthesisType.F90 @@ -355,6 +355,13 @@ subroutine InitCold(this, bounds) do p = bounds%begp,bounds%endp l = veg_pp%landunit(p) + ! These nominal resistance values will be overwritten + ! during the canopy flux iteration process upon first + ! use. We just need something that wont break + ! the model on the first guess + this%rssun_patch(p) = 1.e-5 + this%rssha_patch(p) = 1.e-5 + this%lmrcanopy_patch(p) = 0.0_r8 this%alphapsnsun_patch(p) = spval @@ -406,7 +413,16 @@ subroutine Restart(this, bounds, ncid, flag) dim1name='pft', long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%rc13_psnsha_patch) endif - + + call restartvar(ncid=ncid, flag=flag, varname='rssun_patch',xtype=ncd_double, & + dim1name='pft', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=this%rssun_patch) + + call restartvar(ncid=ncid, flag=flag, varname='rssha_patch',xtype=ncd_double, & + dim1name='pft', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=this%rssha_patch) + + end subroutine Restart !------------------------------------------------------------------------------ subroutine TimeStepInit (this, bounds) diff --git a/components/elm/src/biogeophys/QSatMod.F90 b/components/elm/src/biogeophys/QSatMod.F90 index e6636bf6e979..afa989ffaff2 100644 --- a/components/elm/src/biogeophys/QSatMod.F90 +++ b/components/elm/src/biogeophys/QSatMod.F90 @@ -7,7 +7,7 @@ module QSatMod ! !PUBLIC TYPES: use shr_kind_mod , only: r8 => shr_kind_r8 implicit none - save + private ! ! !PUBLIC MEMBER FUNCTIONS: @@ -60,7 +60,7 @@ module QSatMod !----------------------------------------------------------------------- - subroutine QSat (T, p, es, esdT, qs, qsdT) + elemental subroutine QSat (T, p, es, esdT, qs, qsdT) ! ! !DESCRIPTION: ! Computes saturation mixing ratio and the change in saturation diff --git a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 index f36f74411fb5..41898180dfb7 100644 --- a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 +++ b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 @@ -964,7 +964,7 @@ subroutine SurfaceAlbedo(bounds, & if(use_fates)then #ifndef _OPENACC - call alm_fates%wrap_canopy_radiation(bounds,surfalb_vars,nextsw_cday,declinp1) + call alm_fates%wrap_canopy_radiation(bounds,surfalb_vars,canopystate_vars,nextsw_cday,declinp1) #endif else diff --git a/components/elm/src/dyn_subgrid/dynSubgridDriverMod.F90 b/components/elm/src/dyn_subgrid/dynSubgridDriverMod.F90 index 09ce19b7f98d..0a54b71fc310 100644 --- a/components/elm/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/components/elm/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -144,13 +144,19 @@ subroutine dynSubgrid_init(bounds, glc2lnd_vars, crop_vars) call dyncrop_interp(bounds, crop_vars) end if - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + ! THIS SHOULD FASTER THAN PARALLEL DO (RGK) + !! OMP PARALLEL num_threads(nclumps) PRIVATE (nc, bounds_clump) if(nclumps>1) + !! nc = omp_get_thread_num() + + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) call dynSubgrid_wrapup_weight_changes(bounds_clump, glc2lnd_vars) end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + !! OMP END PARALLEL + end subroutine dynSubgrid_init !----------------------------------------------------------------------- @@ -241,7 +247,7 @@ subroutine dynSubgrid_driver(bounds_proc, & ! Do initialization, prior to land cover change ! ========================================================================== - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) @@ -289,7 +295,7 @@ subroutine dynSubgrid_driver(bounds_proc, & ! Do everything else related to land cover change ! ========================================================================== - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90 index 55dd8d76c4f2..25c1ae8d5a30 100644 --- a/components/elm/src/main/controlMod.F90 +++ b/components/elm/src/main/controlMod.F90 @@ -17,7 +17,7 @@ module controlMod use shr_log_mod , only: errMsg => shr_log_errMsg use abortutils , only: endrun use spmdMod , only: masterproc - use decompMod , only: clump_pproc + use decompMod , only: clump_pproc,fates_pproc use elm_varpar , only: maxpatch_pft, maxpatch_glcmec, more_vertlayers,nlevdecomp_full, nsoilorder use histFileMod , only: max_tapes, max_namlen use histFileMod , only: hist_empty_htapes, hist_dov2xy, hist_avgflag_pertape, hist_type1d_pertape @@ -374,12 +374,8 @@ subroutine control_init( ) ! Set clumps per procoessor -#if (defined _OPENMP) - clump_pproc = omp_get_max_threads() -#else - clump_pproc = 1 -#endif + override_nsrest = nsrest if (masterproc) then @@ -586,8 +582,24 @@ subroutine control_init( ) endif endif +#if (defined _OPENMP) + if(use_fates)then + fates_pproc = omp_get_max_threads() + clump_pproc = 1 + else + clump_pproc = omp_get_max_threads() + fates_pproc = 1 + end if +#else + clump_pproc = 1 + fates_pproc = 1 +#endif + endif ! end of if-masterproc if-block + + + ! ---------------------------------------------------------------------- ! Read in other namelists for other modules ! ---------------------------------------------------------------------- @@ -969,7 +981,14 @@ subroutine control_spmd() ! clump decomposition variables call mpi_bcast (clump_pproc, 1, MPI_INTEGER, 0, mpicom, ier) + call mpi_bcast (fates_pproc, 1, MPI_INTEGER, 0, mpicom, ier) +#if (defined _OPENMP) + if(fates_pproc>1)then + call omp_set_max_active_levels(1) + end if +#endif + ! lateral connectivity call mpi_bcast (lateral_connectivity, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (domain_decomp_type, len(domain_decomp_type), MPI_CHARACTER, 0, mpicom, ier) diff --git a/components/elm/src/main/decompInitMod.F90 b/components/elm/src/main/decompInitMod.F90 index c03856404056..b8d52fbeca37 100644 --- a/components/elm/src/main/decompInitMod.F90 +++ b/components/elm/src/main/decompInitMod.F90 @@ -119,6 +119,8 @@ subroutine decompInit_moab(lni,lnj,amask) write(iulog,*) 'decompInit_moab(): allocation error for procinfo%cid' call endrun(msg=errMsg(__FILE__, __LINE__)) endif + + procinfo%nclumps = clump_pproc procinfo%cid(:) = -1 procinfo%ncells = moab_gcell%num_owned ! owned elements in the current task diff --git a/components/elm/src/main/decompMod.F90 b/components/elm/src/main/decompMod.F90 index ff9fc3383b15..1dcf9f105cc5 100644 --- a/components/elm/src/main/decompMod.F90 +++ b/components/elm/src/main/decompMod.F90 @@ -10,13 +10,19 @@ module decompMod ! Must use shr_sys_abort rather than endrun here to avoid circular dependency use shr_sys_mod , only : shr_sys_abort use elm_varctl , only : iulog + use elm_varctl , only : use_fates use elm_varcon , only : grlnd, nameg, namet, namel, namec, namep, nameCohort use mct_mod , only : mct_gsMap ! ! !PUBLIC TYPES: implicit none integer, public :: clump_pproc ! number of clumps per MPI process - + integer, public :: fates_pproc ! number of clumps (for FATES) per MPI proces + ! This is mutually exclusive with clump_pproc + ! If FATES is active, clump_pproc will be forced + ! to 1, and FATES will reserve the omp_num_threads() + ! for its patch-level processes + ! Define possible bounds subgrid levels integer, parameter, public :: BOUNDS_SUBGRID_GRIDCELL = 1 integer, parameter, public :: BOUNDS_SUBGRID_TOPOUNIT = 2 @@ -297,7 +303,7 @@ subroutine get_clump_bounds_new (n, bounds) !write(iulog,*) 'SPM omp debug decompMod 1 ', & !OMP_GET_NUM_THREADS(),OMP_GET_MAX_THREADS(),OMP_GET_THREAD_NUM() - if ( OMP_GET_NUM_THREADS() == 1 .and. OMP_GET_MAX_THREADS() > 1 )then + if (.not.use_fates .and. ( OMP_GET_NUM_THREADS() == 1 .and. OMP_GET_MAX_THREADS() > 1 ))then call shr_sys_abort( trim(subname)//' ERROR: Calling from inside a non-threaded region)') end if #endif diff --git a/components/elm/src/main/elm_driver.F90 b/components/elm/src/main/elm_driver.F90 index e9ee1045739e..81dbea6bc068 100644 --- a/components/elm/src/main/elm_driver.F90 +++ b/components/elm/src/main/elm_driver.F90 @@ -40,6 +40,7 @@ module elm_driver ! use BareGroundFluxesMod , only : BareGroundFluxes use CanopyFluxesMod , only : CanopyFluxes + use CanopyFluxesMod , only : PatchLoadBalance use SedYieldMod , only : SoilErosion use SoilFluxesMod , only : SoilFluxes ! (formerly Biogeophysics2Mod) use UrbanFluxesMod , only : UrbanFluxes @@ -330,7 +331,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) ! the associated filter updates, too (otherwise we get a carbon balance error) ! ================================================================================== - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -363,7 +364,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) ! Zero fluxes for transient land cover ! ============================================================================ - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -481,7 +482,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) else call t_startf('cnbalchk_at_grid') - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -551,7 +552,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) ! changes due to dynamic area adjustments can break column-level conservation ! ============================================================================ - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -677,7 +678,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) ! snow accumulation exceeds 10 mm. ! ============================================================================ - !$OMP PARALLEL DO PRIVATE (nc,l,c, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc,l,c, bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -1315,6 +1316,11 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) call alm_fates%dynamics_driv( bounds_clump, top_as, & top_af, atm2lnd_vars, soilstate_vars, & canopystate_vars, frictionvel_vars, soil_water_retention_curve) + + call PatchLoadBalance(bounds_clump, filter(nc)%num_nolakeurbanp, & + filter(nc)%nolakeurbanp,canopystate_vars, & + alm_fates%fates(nc)%bc_out,alm_fates%f2hmap(nc)%hsites,.false.) + end if end if @@ -1381,7 +1387,8 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) ! for the two-stream radiation scheme. if (use_fates .and. .not.doalb .and. get_nstep() == 1 .and. nsrest == nsrStartup) then if ( trim(finidat) == '' .or. fates_radiation_model=='twostream') then - call alm_fates%wrap_canopy_radiation(bounds_clump,surfalb_vars,nextsw_cday,declinp1) + call alm_fates%wrap_canopy_radiation(bounds_clump,surfalb_vars, & + canopystate_vars,nextsw_cday,declinp1) end if end if @@ -1449,7 +1456,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) if (create_glacier_mec_landunit) then call t_startf('lnd2glc') - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) call lnd2glc_vars%update_lnd2glc(bounds_clump, & @@ -1466,7 +1473,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) if (iac_present) then call t_startf('lnd2iac') - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) call lnd2iac_vars%update_lnd2iac(bounds_clump) diff --git a/components/elm/src/main/elm_initializeMod.F90 b/components/elm/src/main/elm_initializeMod.F90 index b061e481884d..8e23fff07224 100644 --- a/components/elm/src/main/elm_initializeMod.F90 +++ b/components/elm/src/main/elm_initializeMod.F90 @@ -21,7 +21,6 @@ module elm_initializeMod use readParamsMod , only : readSharedParameters, readPrivateParameters use ncdio_pio , only : file_desc_t use ELMFatesInterfaceMod , only : ELMFatesGlobals1,ELMFatesGlobals2 - use ELMFatesParamInterfaceMod, only: FatesReadPFTs use BeTRSimulationELM, only : create_betr_simulation_elm use SoilLittVertTranspMod, only : CreateLitterTransportList use iso_c_binding @@ -37,7 +36,7 @@ module elm_initializeMod use ColumnDataType , only : col_es use VegetationType , only : veg_pp use VegetationDataType , only : veg_es - + use CanopyFluxesMod , only : PatchLoadBalance use elm_instMod use WaterBudgetMod , only : WaterBudget_Reset use CNPBudgetMod , only : CNPBudget_Reset @@ -345,14 +344,6 @@ subroutine initialize1( ) call soilorder_conrd() - ! Read in FATES parameter values early in the call sequence as well - ! The PFT file, specifically, will dictate how many pfts are used - ! in fates, and this will influence the amount of memory we - ! request from the model, which is relevant in set_fates_global_elements() - if (use_fates) then - call FatesReadPFTs() - end if - ! Read surface dataset and set up subgrid weight arrays call surfrd_get_data(begg, endg, ldomain, fsurdat) @@ -453,7 +444,7 @@ subroutine initialize1( ) call t_stopf('init_filters') nclumps = get_proc_clumps() - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) call reweight_wrapup(bounds_clump, & @@ -908,7 +899,7 @@ subroutine initialize2( ) 'finidat and finidat_interp_source cannot both be non-blank') end if - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) call reweight_wrapup(bounds_clump, & @@ -944,7 +935,7 @@ subroutine initialize2( ) end if - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) call reweight_wrapup(bounds_clump, & @@ -1059,7 +1050,7 @@ subroutine initialize2( ) !------------------------------------------------------------ if (create_glacier_mec_landunit) then - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -1072,6 +1063,16 @@ subroutine initialize2( ) !$OMP END PARALLEL DO end if + ! This load balancing is used in canopy fluxes + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) + do nc = 1,nclumps + call get_clump_bounds(nc, bounds_clump) + call PatchLoadBalance(bounds_clump, filter(nc)%num_nolakeurbanp, & + filter(nc)%nolakeurbanp,canopystate_vars, & + alm_fates%fates(nc)%bc_out,alm_fates%f2hmap(nc)%hsites,.true.) + end do + !$OMP END PARALLEL DO + !------------------------------------------------------------ ! Deallocate wt_nat_patch !------------------------------------------------------------ @@ -1089,7 +1090,7 @@ subroutine initialize2( ) ! If fates is using satellite phenology mode, make sure to call the SatellitePhenology ! procedure prior to init_coldstart which will eventually call leaf_area_profile if ( use_fates_sp ) then - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) call SatellitePhenology(bounds_clump, & diff --git a/components/elm/src/main/elmfates_interfaceMod.F90 b/components/elm/src/main/elmfates_interfaceMod.F90 index 7c7643404d56..650cb8292810 100644 --- a/components/elm/src/main/elmfates_interfaceMod.F90 +++ b/components/elm/src/main/elmfates_interfaceMod.F90 @@ -77,6 +77,7 @@ module ELMFatesInterfaceMod use elm_varctl , only : use_lch4 use elm_varctl , only : use_century_decomp use elm_varctl , only : carbon_only,carbonnitrogen_only,carbonphosphorus_only + use elm_varctl , only : fates_paramfile use elm_varcon , only : tfrz use elm_varcon , only : spval use elm_varcon , only : denice @@ -167,12 +168,15 @@ module ELMFatesInterfaceMod use EDPftVarcon , only : EDpftvarcon_inst use FatesRadiationDriveMod, only : FatesSunShadeFracs use FatesRadiationDriveMod, only : FatesNormalizedCanopyRadiation + use FatesRadiationDriveMod, only : FatesNormalizedPatchRadiation use EDBtranMod , only : btran_ed, & get_active_suction_layers use EDCanopyStructureMod , only : canopy_summarization, update_hlm_dynamics use EDCanopyStructureMod , only : UpdateFatesAvgSnowDepth + use FatesPlantRespPhotosynthMod, only : FatesCondPhotoPatch + use FatesPlantRespPhotosynthMod, only : FatesPlantRespPatch use FatesPlantRespPhotosynthMod, only : FatesPlantRespPhotosynthDrive - use EDAccumulateFluxesMod , only : AccumulateFluxes_ED + use EDAccumulateFluxesMod , only : FatesAccumulatePatchFluxes use FatesSoilBGCFluxMod , only : UnPackNutrientAquisitionBCs use FatesSoilBGCFluxMod , only : FluxIntoLitterPools use PRTGenericMod, only : prt_cnp_flex_allom_hyp @@ -206,10 +210,8 @@ module ELMFatesInterfaceMod use dynFATESLandUseChangeMod, only : fates_harvest_luh_mass use FatesInterfaceTypesMod , only : bc_in_type, bc_out_type - - use ELMFatesParamInterfaceMod, only : fates_param_reader_ctsm_impl - use FatesParametersInterface, only : fates_param_reader_type - use FatesParametersInterface, only : fates_parameters_type + use decompMod , only : fates_pproc + use perf_mod , only : t_startf, t_stopf @@ -267,7 +269,8 @@ module ELMFatesInterfaceMod procedure, public :: wrap_sunfrac procedure, public :: wrap_btran procedure, public :: wrap_photosynthesis - procedure, public :: wrap_accumulatefluxes + procedure, public :: WrapPatchPhotosynthesis + procedure, public :: WrapAccumulateFluxes procedure, public :: prep_canopyfluxes procedure, public :: wrap_canopy_radiation procedure, public :: wrap_WoodProducts @@ -331,7 +334,7 @@ subroutine ELMFatesGlobals1() integer :: pass_use_luh2 integer :: pass_masterproc logical :: verbose_output - type(fates_param_reader_ctsm_impl) :: var_reader + if (use_fates) then @@ -378,6 +381,8 @@ subroutine ELMFatesGlobals1() end if call set_fates_ctrlparms('masterproc',ival=pass_masterproc) + call set_fates_ctrlparms('parteh_mode',ival=fates_parteh_mode) + end if ! The following call reads in the parameter file @@ -394,8 +399,8 @@ subroutine ELMFatesGlobals1() ! want fates to handle crops, so again, it should be ignored. ! (RGK 07-2022) - call SetFatesGlobalElements1(use_fates,natpft_size,0,var_reader) - + call SetFatesGlobalElements1(use_fates,natpft_size,0,fates_paramfile) + natpft_size = fates_maxPatchesPerSite max_patch_per_col= max(natpft_size, numcft, maxpatch_urb) @@ -476,7 +481,6 @@ subroutine ELMFatesGlobals2() call set_fates_ctrlparms('hlm_name',cval='ELM') call set_fates_ctrlparms('hio_ignore_val',rval=spval) call set_fates_ctrlparms('soilwater_ipedof',ival=get_ipedof(0)) - call set_fates_ctrlparms('parteh_mode',ival=fates_parteh_mode) call set_fates_ctrlparms('seeddisp_cadence',ival=fates_seeddisp_cadence) call set_fates_ctrlparms('hist_hifrq_dimlevel',ival=fates_history_dimlevel(1)) @@ -842,7 +846,6 @@ subroutine init(this, bounds_proc, flandusepftdat) use spmdMod, only : npes use decompMod, only : procinfo - use FatesInterfaceMod, only : FatesReportParameters use FatesParameterDerivedMod, only : param_derived use FatesInterfaceTypesMod, only : numpft_fates => numpft use elm_varsur, only : wt_nat_patch @@ -907,8 +910,7 @@ subroutine init(this, bounds_proc, flandusepftdat) end if nclumps = get_proc_clumps() - - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,nmaxcol,s,c,l,g,collist,pi,pf,ft) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,nmaxcol,s,c,l,g,collist,pi,pf,ft) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -918,7 +920,7 @@ subroutine init(this, bounds_proc, flandusepftdat) ! Allocate the mapping that points columns to FATES sites, 0 is NA allocate(this%f2hmap(nc)%hsites(bounds_clump%begc:bounds_clump%endc)) - + ! Initialize all columns with a zero index, which indicates no FATES site this%f2hmap(nc)%hsites(:) = 0 @@ -1057,11 +1059,11 @@ subroutine init(this, bounds_proc, flandusepftdat) call InitHydrSites(this%fates(nc)%sites,this%fates(nc)%bc_in) end if - if( this%fates(nc)%nsites == 0 ) then - write(iulog,*) 'Clump ',nc,' had no valid FATES sites' - write(iulog,*) 'This will likely cause problems until code is improved' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + !if( this%fates(nc)%nsites == 0 ) then + ! write(iulog,*) 'Clump ',nc,' had no valid FATES sites' + ! write(iulog,*) 'This will likely cause problems until code is improved' + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + !end if ! Set patch itypes on natural veg columns to nonsense ! This will force a crash if the model outside of FATES tries to think @@ -1079,9 +1081,6 @@ subroutine init(this, bounds_proc, flandusepftdat) call this%init_history_io(bounds_proc) - ! Report Fates Parameters (debug flag in lower level routines) - call FatesReportParameters(masterproc) - ! Fire data to send to FATES call create_fates_fire_data_method( this%fates_fire_data_method ) @@ -1395,11 +1394,13 @@ subroutine dynamics_driv(this, bounds_clump, top_as_inst, & ! Part IV: ! Update history IO fields that depend on ecosystem dynamics ! --------------------------------------------------------------------------------- + call t_startf('hbuf_fatesdyn') call fates_hist%update_history_dyn( nc, & this%fates(nc)%nsites, & this%fates(nc)%sites, & this%fates(nc)%bc_in) - + call t_stopf('hbuf_fatesdyn') + if (masterproc) then write(iulog, *) 'FATES dynamics complete' end if @@ -1711,12 +1712,12 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump, & z0m(p) = this%fates(nc)%bc_out(s)%z0m_pa(ifp) displa(p) = this%fates(nc)%bc_out(s)%displa_pa(ifp) dleaf_patch(p) = this%fates(nc)%bc_out(s)%dleaf_pa(ifp) - - - end do + + + if(abs(areacheck - 1.0_r8).gt.1.e-9_r8)then write(iulog,*) 'area wrong in interface',areacheck - 1.0_r8 call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -1824,7 +1825,7 @@ subroutine restart( this, bounds_proc, ncid, flag, & call this%fates_restart%Init(nclumps, fates_bounds) ! Define the bounds on the first dimension for each thread - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,fates_clump) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,fates_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -1834,7 +1835,7 @@ subroutine restart( this, bounds_proc, ncid, flag, & end do !$OMP END PARALLEL DO - !$OMP PARALLEL DO PRIVATE (nc,s,c,g) + !$OMP PARALLEL DO PRIVATE (nc,s,c,g) if(nclumps>1) do nc = 1,nclumps allocate(this%fates_restart%restart_map(nc)%site_index(this%fates(nc)%nsites)) @@ -1870,7 +1871,7 @@ subroutine restart( this, bounds_proc, ncid, flag, & ! --------------------------------------------------------------------------------- if(flag=='write')then - !$OMP PARALLEL DO PRIVATE (nc) + !$OMP PARALLEL DO PRIVATE (nc) if(nclumps>1) do nc = 1, nclumps if (this%fates(nc)%nsites>0) then call this%fates_restart%set_restart_vectors(nc,this%fates(nc)%nsites, & @@ -1943,7 +1944,7 @@ subroutine restart( this, bounds_proc, ncid, flag, & if(flag=='read')then - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,s) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,s) if(nclumps>1) do nc = 1, nclumps if (this%fates(nc)%nsites>0) then @@ -2115,7 +2116,7 @@ subroutine init_coldstart(this, canopystate_inst, soilstate_inst, frictionvel_in nclumps = get_proc_clumps() - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,s,c,j,vol_ice,eff_porosity) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,s,c,j,vol_ice,eff_porosity) if(nclumps>1) do nc = 1, nclumps if ( this%fates(nc)%nsites>0 ) then @@ -2374,17 +2375,16 @@ subroutine prep_canopyfluxes(this, bounds_clump ) integer :: c,s integer :: nc + if (.not.use_fates_planthydro) return + nc = bounds_clump%clump_index do s = 1, this%fates(nc)%nsites - ! filter flag == 1 means that this patch has not been called for photosynthesis - this%fates(nc)%bc_in(s)%filter_photo_pa(:) = 1 ! set transpiration input boundary condition to zero. The exposed ! vegetation filter may not even call every patch. - if (use_fates_planthydro) then - this%fates(nc)%bc_in(s)%qflx_transp_pa(:) = 0._r8 - end if - + + this%fates(nc)%bc_in(s)%qflx_transp_pa(:) = 0._r8 + end do end subroutine prep_canopyfluxes @@ -2695,46 +2695,113 @@ subroutine wrap_photosynthesis(this, bounds_clump, fn, filterp, & end subroutine wrap_photosynthesis - ! ====================================================================================== + ! ==================================================================================== + + subroutine WrapPatchPhotosynthesis(this, bounds_clump, p, & + esat_tv, eair, oair, cair, rb, dayl_factor, & + atm2lnd_inst, canopystate_inst, photosyns_inst) + + use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun + use decompMod , only : bounds_type + use elm_varcon , only : rgas, tfrz, namep + use elm_varctl , only : iulog + use quadraticMod , only : quadratic + + ! + ! !ARGUMENTS: + class(hlm_fates_interface_type), intent(inout) :: this + type(bounds_type) , intent(in) :: bounds_clump + integer , intent(in) :: p ! patch of interest + real(r8) , intent(in) :: esat_tv(bounds_clump%begp: ) ! saturation vapor pressure at t_veg (Pa) + real(r8) , intent(in) :: eair( bounds_clump%begp: ) ! vapor pressure of canopy air (Pa) + real(r8) , intent(in) :: oair( bounds_clump%begp: ) ! Atmospheric O2 partial pressure (Pa) + real(r8) , intent(in) :: cair( bounds_clump%begp: ) ! Atmospheric CO2 partial pressure (Pa) + real(r8) , intent(in) :: rb( bounds_clump%begp: ) ! boundary layer resistance (s/m) + real(r8) , intent(in) :: dayl_factor( bounds_clump%begp: ) ! scalar (0-1) for daylength + type(atm2lnd_type) , intent(in) :: atm2lnd_inst + type(canopystate_type) , intent(inout) :: canopystate_inst + type(photosyns_type) , intent(inout) :: photosyns_inst + + integer :: nlevsoil + integer :: s,t,c,ifp,j,icp,nc + real(r8) :: dtime + + call t_startf('edpsn') + associate(& + t_soisno => col_es%t_soisno , & + t_veg => veg_es%t_veg , & + tgcm => veg_es%thm , & + forc_pbot => top_as%pbot , & + rssun => photosyns_inst%rssun_patch , & + rssha => photosyns_inst%rssha_patch) + + + nc = bounds_clump%clump_index + c = veg_pp%column(p) + t = col_pp%topounit(c) + s = this%f2hmap(nc)%hsites(c) + ifp = p - col_pp%pfti(c) - subroutine wrap_accumulatefluxes(this, bounds_clump, fn, filterp) + ! This should be done during prep! + nlevsoil = this%fates(nc)%bc_in(s)%nlevsoil + this%fates(nc)%bc_in(s)%t_soisno_sl(1:nlevsoil) = t_soisno(c,1:nlevsoil) ! soil temperature (Kelvin) + this%fates(nc)%bc_in(s)%forc_pbot = forc_pbot(t) ! atmospheric pressure (Pa) + this%fates(nc)%bc_in(s)%dayl_factor_pa(ifp) = dayl_factor(p) ! scalar (0-1) for daylength + this%fates(nc)%bc_in(s)%oair_pa(ifp) = oair(p) ! Atmospheric O2 partial pressure (Pa) + this%fates(nc)%bc_in(s)%cair_pa(ifp) = cair(p) ! Atmospheric CO2 partial pressure (Pa) + + ! These probably all dynamically respond to the LEB calculation + this%fates(nc)%bc_in(s)%esat_tv_pa(ifp) = esat_tv(p) ! saturation vapor pressure at t_veg (Pa) + this%fates(nc)%bc_in(s)%eair_pa(ifp) = eair(p) ! vapor pressure of canopy air (Pa) + this%fates(nc)%bc_in(s)%rb_pa(ifp) = rb(p) ! boundary layer resistance (s/m) + this%fates(nc)%bc_in(s)%t_veg_pa(ifp) = t_veg(p) ! vegetation temperature (Kelvin) + this%fates(nc)%bc_in(s)%tgcm_pa(ifp) = tgcm(p) ! air temperature at agcm + + call FatesCondPhotoPatch(ifp,this%fates(nc)%sites(s),this%fates(nc)%bc_in(s),this%fates(nc)%bc_out(s)) + + rssun(p) = this%fates(nc)%bc_out(s)%rssun_pa(ifp) + rssha(p) = this%fates(nc)%bc_out(s)%rssha_pa(ifp) + + photosyns_inst%psnsun_patch(p) = spval + photosyns_inst%psnsha_patch(p) = spval + + ! This flags that this patch has experienced photosynthesis + this%fates(nc)%bc_in(s)%filter_photo_pa(ifp) = 3 + + end associate + call t_stopf('edpsn') + + end subroutine WrapPatchPhotosynthesis + + ! ====================================================================================== + + subroutine WrapAccumulateFluxes(this, bounds_clump, p) ! !ARGUMENTS: class(hlm_fates_interface_type), intent(inout) :: this type(bounds_type) , intent(in) :: bounds_clump - integer , intent(in) :: fn ! size of pft filter - integer , intent(in) :: filterp(fn) ! pft filter + integer , intent(in) :: p ! Locals - integer :: s,c,p,ifp,icp + integer :: s,c,ifp real(r8) :: dtime integer :: nc + dtime = real(get_step_size(),r8) nc = bounds_clump%clump_index - ! Run a check on the filter - do icp = 1,fn - p = filterp(icp) - c = veg_pp%column(p) - s = this%f2hmap(nc)%hsites(c) - ifp = p-col_pp%pfti(c) - if(this%fates(nc)%bc_in(s)%filter_photo_pa(ifp) /= 3)then - write(iulog,*) 'Not all patches on the natveg column in the canopys' - write(iulog,*) 'filter ran canopy fluxes: s, p, icp, ifp: ',s,p,icp,ifp - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end do - - - dtime = real(get_step_size(),r8) - - call AccumulateFluxes_ED(this%fates(nc)%nsites, & - this%fates(nc)%sites, & - this%fates(nc)%bc_in, & - this%fates(nc)%bc_out, & - dtime) - return - end subroutine wrap_accumulatefluxes + c = veg_pp%column(p) + s = this%f2hmap(nc)%hsites(c) + ifp = p-col_pp%pfti(c) + call FatesPlantRespPatch(ifp,this%fates(nc)%sites(s),this%fates(nc)%bc_in(s),dtime) + + call FatesAccumulatePatchFluxes(ifp,this%fates(nc)%sites(s), dtime) + + + return + end subroutine WrapAccumulateFluxes + ! ====================================================================================== subroutine wrap_WoodProducts(this, bounds_clump, fc, filterc) @@ -2864,18 +2931,21 @@ end subroutine wrap_FatesCarbonStocks ! ====================================================================================== - subroutine wrap_canopy_radiation(this, bounds_clump, surfalb_inst,nextsw_cday,declinp1) + subroutine wrap_canopy_radiation(this, bounds_clump, surfalb_inst,canopystate_vars,nextsw_cday,declinp1) use shr_orb_mod, only: shr_orb_cosz - + use omp_lib + ! Arguments class(hlm_fates_interface_type), intent(inout) :: this type(bounds_type), intent(in) :: bounds_clump type(surfalb_type) , intent(inout) :: surfalb_inst + type(canopystate_type) , intent(in) :: canopystate_vars real(r8),intent(in) :: nextsw_cday,declinp1 ! locals integer :: s,c,p,ifp,icp,nc,g + integer :: it,np,f associate(& albgrd_col => surfalb_inst%albgrd_col , & !in @@ -2899,16 +2969,7 @@ subroutine wrap_canopy_radiation(this, bounds_clump, surfalb_inst,nextsw_cday,de coszen_col(c) = shr_orb_cosz (nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1) this%fates(nc)%bc_in(s)%coszen = coszen_col(c) - - do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno - - p = ifp+col_pp%pfti(c) - if (veg_es%t_veg(p) <= tfrz) then - this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = veg_ws%fwet(p) - else - this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = 0._r8 - end if - end do + this%fates(nc)%sites(s)%coszen = this%fates(nc)%bc_in(s)%coszen if(coszen_col(c) > 0._r8) then @@ -2921,12 +2982,52 @@ subroutine wrap_canopy_radiation(this, bounds_clump, surfalb_inst,nextsw_cday,de this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = spval end if + + do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno + + p = ifp+col_pp%pfti(c) + if (veg_es%t_veg(p) <= tfrz) then + this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = veg_ws%fwet(p) + else + this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = 0._r8 + end if + + !call FatesNormalizedPatchRadiation(ifp, & + ! this%fates(nc)%sites(s), & + ! this%fates(nc)%bc_in(s), & + ! this%fates(nc)%bc_out(s)) + + end do + + end do - - call FatesNormalizedCanopyRadiation( & - this%fates(nc)%sites, & - this%fates(nc)%bc_in, & - this%fates(nc)%bc_out) + + !$OMP PARALLEL num_threads(fates_pproc) PRIVATE (it,np,f,p,c,ifp,s) if(fates_pproc>1) + it = 0 + if(fates_pproc>1) then +#ifdef _OPENMP + it = omp_get_thread_num() +#endif + end if + np = canopystate_vars%patch_par(it)%npatch + patch_loop: do f = 1, np + p = canopystate_vars%patch_par(it)%patch_list(f) + if(veg_pp%is_fates(p))then + c = veg_pp%column(p) + ifp = p - col_pp%pfti(c) + s = this%f2hmap(nc)%hsites(c) + call FatesNormalizedPatchRadiation(ifp, & + this%fates(nc)%sites(s), & + this%fates(nc)%bc_in(s), & + this%fates(nc)%bc_out(s)) + end if + end do patch_loop + !$OMP END PARALLEL + + !call FatesNormalizedCanopyRadiation( & + ! this%fates(nc)%sites, & + ! this%fates(nc)%bc_in, & + ! this%fates(nc)%bc_out) ! Pass FATES BC's back to HLM ! ----------------------------------------------------------------------------------- @@ -3110,13 +3211,14 @@ subroutine wrap_update_hifrq_hist(this, bounds_clump ) end do ! Update history variables that track these variables + call t_startf('hbuf_fateshifrq') call fates_hist%update_history_hifrq(nc, & this%fates(nc)%nsites, & this%fates(nc)%sites, & this%fates(nc)%bc_in, & this%fates(nc)%bc_out, & dtime) - + call t_stopf('hbuf_fateshifrq') end associate end subroutine wrap_update_hifrq_hist @@ -3397,7 +3499,7 @@ subroutine init_history_io(this,bounds_proc) call fates_hist%Init(nclumps, fates_bounds) ! Define the bounds on the first dimension for each thread - !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,fates_clump) + !$OMP PARALLEL DO PRIVATE (nc,bounds_clump,fates_clump) if(nclumps>1) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -3617,7 +3719,8 @@ subroutine wrap_hydraulics_drive(this, bounds_clump, & real(r8) :: dtime integer :: nlevsoil - + !this%fates(nc)%bc_in(s)%qflx_transp_pa(:) = 0._r8 + if ( .not.use_fates_planthydro ) return nc = bounds_clump%clump_index diff --git a/components/elm/src/main/elmfates_paraminterfaceMod.F90 b/components/elm/src/main/elmfates_paraminterfaceMod.F90 deleted file mode 100644 index deb069260835..000000000000 --- a/components/elm/src/main/elmfates_paraminterfaceMod.F90 +++ /dev/null @@ -1,225 +0,0 @@ -module ELMFatesParamInterfaceMod - ! NOTE(bja, 2017-01) this code can not go into the main clm-fates - ! interface module because of circular dependancies with pftvarcon. - - use FatesGlobals, only : fates_log - use FatesParametersInterface, only : fates_parameters_type - use FatesParametersInterface, only : fates_param_reader_type - use shr_kind_mod, only : r8 => shr_kind_r8 - - implicit none - private - - ! Extend the fates parameter reader type with Read parameter function - type, extends(fates_param_reader_type) :: fates_param_reader_ctsm_impl - ! !PRIVATE MEMBER DATA: - contains - ! !PUBLIC MEMBER FUNCTIONS: - procedure, public :: Read ! Read params from disk - end type - - public :: fates_param_reader_ctsm_impl - - ! NOTE(bja, 2017-01) these methods can NOT be part of the hlm-fates - ! interface type because they are called before the instance is - ! initialized. - public :: FatesReadPFTs - - logical :: DEBUG = .false. - - character(len=*), parameter, private :: sourcefile = & - __FILE__ - -contains - - - !----------------------------------------------------------------------- - subroutine FatesReadPFTs() - - use elm_varctl, only : use_fates, paramfile, fates_paramfile - use spmdMod, only : masterproc - - use EDPftvarcon , only : EDPftvarcon_inst - - use fileutils , only : getfil - use ncdio_pio , only : file_desc_t, ncd_pio_closefile, ncd_pio_openfile - - implicit none - - character(len=32) :: subname = 'FatesReadPFTs' - class(fates_parameters_type), allocatable :: fates_params - logical :: is_host_file - - character(len=256) :: locfn ! local file name - type(file_desc_t) :: ncid ! pio netCDF file id - - if (use_fates) then - if (masterproc) then - write(fates_log(), *) 'clmfates_interfaceMod.F90::'//trim(subname)//' :: CLM reading ED/FATES '//' PFTs ' - end if - - allocate(fates_params) - call fates_params%Init() - call EDPftvarcon_inst%Init() - - call EDPftvarcon_inst%Register(fates_params) - - is_host_file = .false. - call ParametersFromNetCDF(fates_paramfile, is_host_file, fates_params) - - is_host_file = .true. - call ParametersFromNetCDF(paramfile, is_host_file, fates_params) - - call EDPftvarcon_inst%Receive(fates_params) - - call fates_params%Destroy() - deallocate(fates_params) - end if - - end subroutine FatesReadPFTs - - !----------------------------------------------------------------------- - subroutine SetParameterDimensions(ncid, is_host_file, fates_params) - ! Get the list of dimensions used by the fates parameters, - ! retreive them from the parameter file, then give the information - ! back to fates. - use FatesParametersInterface, only : fates_parameters_type, param_string_length, max_dimensions, max_used_dimensions - use ncdio_pio , only : file_desc_t - - implicit none - - type(file_desc_t), intent(inout) :: ncid - logical, intent(in) :: is_host_file - class(fates_parameters_type), intent(inout) :: fates_params - - integer :: num_used_dimensions - character(len=param_string_length) :: used_dimension_names(max_used_dimensions) - integer :: used_dimension_sizes(max_used_dimensions) - - call fates_params%GetUsedDimensions(is_host_file, num_used_dimensions, used_dimension_names) - - call GetUsedDimensionSizes(ncid, num_used_dimensions, used_dimension_names, used_dimension_sizes) - - call fates_params%SetDimensionSizes(is_host_file, num_used_dimensions, used_dimension_names, used_dimension_sizes) - - end subroutine SetParameterDimensions - - !----------------------------------------------------------------------- - subroutine GetUsedDimensionSizes(ncid, num_used_dimensions, dimension_names, dimension_sizes) - - use ncdio_pio , only : ncd_inqdid, ncd_inqdlen - use FatesParametersInterface, only : param_string_length - use ncdio_pio, only : file_desc_t - - - implicit none - - type(file_desc_t), intent(inout) :: ncid - integer, intent(in) :: num_used_dimensions - character(len=param_string_length), intent(in) :: dimension_names(:) - integer, intent(out) :: dimension_sizes(:) - - integer :: d, max_dim_size, num_dims - integer :: dim_len, dim_id - - dimension_sizes(:) = 0 - max_dim_size = 0 - - do d = 1, num_used_dimensions - call ncd_inqdid(ncid, dimension_names(d), dim_id) - call ncd_inqdlen(ncid, dim_id, dim_len) - dimension_sizes(d) = dim_len - !write(*, *) '--> ', trim(dimension_names(d)), ' setting size ', dimension_sizes(d) - end do - - end subroutine GetUsedDimensionSizes - - !----------------------------------------------------------------------- - subroutine ParametersFromNetCDF(filename, is_host_file, fates_params) - - use abortutils, only : endrun - use fileutils , only : getfil - use ncdio_pio , only : file_desc_t, ncd_pio_closefile, ncd_pio_openfile - use paramUtilMod, only : readNcdio - - use FatesParametersInterface, only : fates_parameters_type - use FatesParametersInterface, only : param_string_length, max_dimensions, max_used_dimensions - use FatesParametersInterface, only : dimension_shape_scalar, dimension_shape_1d, dimension_shape_2d - - implicit none - - character(len=*), intent(in) :: filename - logical, intent(in) :: is_host_file - class(fates_parameters_type), intent(inout) :: fates_params - - character(len=40) :: subname = 'clmfates_interface::ReadParameters' - character(len=256) :: locfn ! local file name - type(file_desc_t) :: ncid ! pio netCDF file id - integer :: dimid ! netCDF dimension id - integer :: i, num_params, dimension_shape - integer :: max_dim_size - real(r8), allocatable :: data(:, :) - character(len=param_string_length) :: name - integer :: dimension_sizes(max_dimensions) - character(len=param_string_length) :: dimension_names(max_dimensions) - integer :: size_dim_1, size_dim_2 - logical :: is_host_param - - call getfil (filename, locfn, 0) - call ncd_pio_openfile (ncid, trim(locfn), 0) - - call SetParameterDimensions(ncid, is_host_file, fates_params) - max_dim_size = fates_params%GetMaxDimensionSize() - allocate(data(max_dim_size, max_dim_size)) - - num_params = fates_params%num_params() - do i = 1, num_params - call fates_params%GetMetaData(i, name, dimension_shape, dimension_sizes, dimension_names, is_host_param) - if (is_host_file .eqv. is_host_param) then - select case(dimension_shape) - case(dimension_shape_scalar) - size_dim_1 = 1 - size_dim_2 = 1 - case(dimension_shape_1d) - size_dim_1 = dimension_sizes(1) - size_dim_2 = 1 - case(dimension_shape_2d) - size_dim_1 = dimension_sizes(1) - size_dim_2 = dimension_sizes(2) - case default - write(fates_log(),*) 'dimension shape:',dimension_shape - call endrun(msg='unsupported number of dimensions reading parameters.') - end select - if(DEBUG) then - write(fates_log(), *) 'clmfates_interfaceMod.F90:: reading '//trim(name) - end if - call readNcdio(ncid, name, dimension_shape, dimension_names, subname, data(1:size_dim_1, 1:size_dim_2)) - call fates_params%SetData(i, data(1:size_dim_1, 1:size_dim_2)) - end if - end do - deallocate(data) - call ncd_pio_closefile(ncid) - end subroutine ParametersFromNetCDF - - !----------------------------------------------------------------------- - - subroutine Read(this, fates_params ) - ! - ! !DESCRIPTION: - ! Read 'fates_params' parameters from storage. - ! - ! USES - use elm_varctl, only : fname_len, paramfile, fates_paramfile - ! !ARGUMENTS: - class(fates_param_reader_ctsm_impl) :: this - class(fates_parameters_type), intent(inout) :: fates_params - !----------------------------------------------------------------------- - logical :: is_host_file = .false. - - call ParametersFromNetCDF(fates_paramfile, is_host_file, fates_params) - - end subroutine Read - - !----------------------------------------------------------------------- - -end module ELMFatesParamInterfaceMod diff --git a/components/elm/src/main/filterMod.F90 b/components/elm/src/main/filterMod.F90 index 80d376abd889..015e94fa061b 100644 --- a/components/elm/src/main/filterMod.F90 +++ b/components/elm/src/main/filterMod.F90 @@ -178,7 +178,7 @@ subroutine allocFiltersOneGroup(this_filter) ! Loop over clumps on this processor -!$OMP PARALLEL DO PRIVATE (nc,bounds) +!$OMP PARALLEL DO PRIVATE (nc,bounds) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds) diff --git a/components/elm/src/main/initGridCellsMod.F90 b/components/elm/src/main/initGridCellsMod.F90 index a0e99dd2102d..58d52cc42917 100644 --- a/components/elm/src/main/initGridCellsMod.F90 +++ b/components/elm/src/main/initGridCellsMod.F90 @@ -114,7 +114,7 @@ subroutine initGridcells nclumps = get_proc_clumps() ! FIX(SPM,032414) add private vars for cohort and perhaps patch dimension - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump, ti, li, ci, pi, gdc, topounit) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump, ti, li, ci, pi, gdc, topounit) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) diff --git a/components/elm/src/main/restFileMod.F90 b/components/elm/src/main/restFileMod.F90 index 291178b98233..442de0dfbde2 100644 --- a/components/elm/src/main/restFileMod.F90 +++ b/components/elm/src/main/restFileMod.F90 @@ -527,7 +527,7 @@ subroutine restFile_read( bounds, file, & ! The reweight_wrapup call needs to be done inside a clump loop, so we set that up ! here. nclumps = get_proc_clumps() - !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) + !$OMP PARALLEL DO PRIVATE (nc, bounds_clump) if(nclumps>1) do nc = 1, nclumps call get_clump_bounds(nc, bounds_clump) call reweight_wrapup(bounds_clump, glc2lnd_vars%icemask_grc(bounds_clump%begg:bounds_clump%endg))