diff --git a/.gitmodules b/.gitmodules index 789f67889..d253f6966 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,5 +8,5 @@ branch = master [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/NCAR/ccpp-physics.git - branch = master + url = https://github.com/NCAR/ccpp-physics + branch = master diff --git a/atmos_model.F90 b/atmos_model.F90 index 900c9143e..35b2c4ceb 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -297,7 +297,9 @@ subroutine update_atmos_radiation_physics (Atmos) #endif !--- call stochastic physics pattern generation / cellular automata - call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block) + call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr) + if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') + !--- if coupled, assign coupled fields @@ -628,7 +630,8 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step) #endif !--- Initialize stochastic physics pattern generation / cellular automata for first time step - call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block) + call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr) + if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') Atmos%Diag => IPD_Diag diff --git a/ccpp/framework b/ccpp/framework index 209f1c92d..f5d4cd2bf 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit 209f1c92d99b7d4cc63e0d41c652fcfd730bd9fa +Subproject commit f5d4cd2bf7752ebf1e4ed16dfdfae71dbfabfb76 diff --git a/ccpp/physics b/ccpp/physics index 8617587ed..4c17ff716 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 8617587edb95aa097b7bbc2735990393bc6d9b90 +Subproject commit 4c17ff716a92d8ac0261d6cc2365bbd7a752b74a diff --git a/gfsphysics/CMakeLists.txt b/gfsphysics/CMakeLists.txt index 34faf6f4c..2c5040588 100644 --- a/gfsphysics/CMakeLists.txt +++ b/gfsphysics/CMakeLists.txt @@ -17,7 +17,6 @@ endif() set(CCPP_SOURCES physics/mersenne_twister.f physics/namelist_soilveg.f - physics/physparam.f physics/set_soilveg.f physics/noahmp_tables.f90 diff --git a/gfsphysics/GFS_layer/GFS_driver.F90 b/gfsphysics/GFS_layer/GFS_driver.F90 index 4abde43d0..e32ec4c26 100644 --- a/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_driver.F90 @@ -633,6 +633,7 @@ subroutine GFS_time_vary_step (Model, Statein, Stateout, Sfcprop, Coupling, & if (mod(Model%kdt,Model%nscyc) == 1) THEN call gcycle (nblks, Model, Grid(:), Sfcprop(:), Cldprop(:)) endif + ! if not updating surface params through fcast, perturb params once at start of fcast endif !--- determine if diagnostics buckets need to be cleared diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 02eb00e00..899955f03 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -660,6 +660,7 @@ subroutine GFS_physics_driver & real :: pshltr,QCQ,rh02 real(kind=kind_phys), allocatable, dimension(:,:) :: den + real(kind=kind_phys) :: lndp_vgf !! Initialize local variables (for debugging purposes only, !! because the corresponding variables Interstitial(nt)%... !! are reset to zero every time). @@ -928,34 +929,28 @@ subroutine GFS_physics_driver & ! alb1d(i) = zero vegf1d(i) = zero enddo - if (Model%do_sfcperts) then - if (Model%pertz0(1) > zero) then - z01d(:) = Model%pertz0(1) * Coupling%sfc_wts(:,1) -! if (me == 0) print*,'Coupling%sfc_wts(:,1) min and max',minval(Coupling%sfc_wts(:,1)),maxval(Coupling%sfc_wts(:,1)) -! if (me == 0) print*,'z01d min and max ',minval(z01d),maxval(z01d) - endif - if (Model%pertzt(1) > zero) then - zt1d(:) = Model%pertzt(1) * Coupling%sfc_wts(:,2) - endif - if (Model%pertshc(1) > zero) then - bexp1d(:) = Model%pertshc(1) * Coupling%sfc_wts(:,3) - endif - if (Model%pertlai(1) > zero) then - xlai1d(:) = Model%pertlai(1) * Coupling%sfc_wts(:,4) - endif -! --- do the albedo percentile calculation in GFS_radiation_driver instead --- ! -! if (Model%pertalb(1) > zero) then -! do i=1,im -! call cdfnor(Coupling%sfc_wts(i,5),cdfz) -! alb1d(i) = cdfz -! enddo -! endif - if (Model%pertvegf(1) > zero) then - do i=1,im - call cdfnor(Coupling%sfc_wts(i,6),cdfz) - vegf1d(i) = cdfz - enddo - endif + lndp_vgf=-999. + + if (Model%lndp_type==1) then + do k =1,Model%n_var_lndp + select case(Model%lndp_var_list(k)) + case ('rz0') + z01d(:) = Model%lndp_prt_list(k)* Coupling%sfc_wts(:,k) + case ('rzt') + zt1d(:) = Model%lndp_prt_list(k)* Coupling%sfc_wts(:,k) + case ('shc') + bexp1d(:) = Model%lndp_prt_list(k) * Coupling%sfc_wts(:,k) + case ('lai') + xlai1d(:) = Model%lndp_prt_list(k)* Coupling%sfc_wts(:,k) + case ('vgf') +! note that the pertrubed vegfrac is being used in sfc_drv, but not sfc_diff + do i=1,im + call cdfnor(Coupling%sfc_wts(i,k),cdfz) + vegf1d(i) = cdfz + enddo + lndp_vgf = Model%lndp_prt_list(k) + end select + enddo endif !*## CCPP ## ! @@ -1856,6 +1851,7 @@ subroutine GFS_physics_driver & ! &,' pgr=',pgr(ipr),' sfcemis=',sfcemis(ipr) !## CCPP ##* sfc_drv.f/lsm_noah_run + call sfc_drv & ! --- inputs: (im, lsoil, Statein%pgr, & @@ -1867,7 +1863,8 @@ subroutine GFS_physics_driver & Sfcprop%shdmin, Sfcprop%shdmax, Sfcprop%snoalb, & Radtend%sfalb, flag_iter, flag_guess, Model%lheatstrg, & Model%isot, Model%ivegsrc, & - bexp1d, xlai1d, vegf1d, Model%pertvegf, & + bexp1d, xlai1d, vegf1d,lndp_vgf, & + ! --- input/output: weasd3(:,1), snowd3(:,1), tsfc3(:,1), tprcp3(:,1), & Sfcprop%srflag, smsoil, stsoil, slsoil, Sfcprop%canopy, & diff --git a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 index ebec30c4d..da5078f7b 100644 --- a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 @@ -1221,6 +1221,7 @@ subroutine GFS_radiation_driver & ! mg, sfc perts real(kind=kind_phys), dimension(size(Grid%xlon,1)) :: alb1d + real(kind=kind_phys) :: lndp_alb real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+ltp) :: cldtausw real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+ltp) :: cldtaulw @@ -1846,14 +1847,18 @@ subroutine GFS_radiation_driver & ! --- scale random patterns for surface perturbations with ! perturbation size ! --- turn vegetation fraction pattern into percentile pattern - alb1d(:) = zero - if (Model%do_sfcperts) then - if (Model%pertalb(1) > zero) then - do i=1,im - call cdfnor(Coupling%sfc_wts(i,5),alb1d(i)) - enddo + alb1d(:) = 0. + lndp_alb = -999. + if (Model%lndp_type ==1) then + do k =1,Model%n_var_lndp + if (Model%lndp_var_list(k) == 'alb') then + do i=1,im + call cdfnor(Coupling%sfc_wts(i,k),alb1d(i)) + lndp_alb = Model%lndp_prt_list(k) + enddo + endif + enddo endif - endif ! mg, sfc-perts !*## CCPP ## @@ -1870,7 +1875,7 @@ subroutine GFS_radiation_driver & Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf, & Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice, & Sfcprop%tisfc, IM, & - alb1d, Model%pertalb, & ! mg, sfc-perts + alb1d, lndp_alb, & ! mg, sfc-perts sfcalb) ! --- outputs !> -# Approximate mean surface albedo from vis- and nir- diffuse values. diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 3c1252b1e..4d3140d24 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -517,7 +517,6 @@ module GFS_typedefs real (kind=kind_phys), pointer :: skebu_wts (:,:) => null() ! real (kind=kind_phys), pointer :: skebv_wts (:,:) => null() ! real (kind=kind_phys), pointer :: sfc_wts (:,:) => null() ! mg, sfc-perts - integer :: nsfcpert=6 !< number of sfc perturbations !--- aerosol surface emissions for Thompson microphysics real (kind=kind_phys), pointer :: nwfa2d (:) => null() !< instantaneous water-friendly sfc aerosol source @@ -1045,14 +1044,14 @@ module GFS_typedefs logical :: do_shum logical :: do_skeb integer :: skeb_npass - logical :: do_sfcperts - integer :: nsfcpert=6 - real(kind=kind_phys) :: pertz0(5) ! mg, sfc-perts - real(kind=kind_phys) :: pertzt(5) ! mg, sfc-perts - real(kind=kind_phys) :: pertshc(5) ! mg, sfc-perts - real(kind=kind_phys) :: pertlai(5) ! mg, sfc-perts - real(kind=kind_phys) :: pertalb(5) ! mg, sfc-perts - real(kind=kind_phys) :: pertvegf(5) ! mg, sfc-perts + integer :: lndp_type + integer :: n_var_lndp + character(len=3) :: lndp_var_list(6) ! dimension here must match n_var_max_lndp in stochy_nml_def + real(kind=kind_phys) :: lndp_prt_list(6) ! dimension here must match n_var_max_lndp in stochy_nml_def + ! also previous code had dimension 5 for each pert, to allow + ! multiple patterns. It wasn't fully coded (and wouldn't have worked + ! with nlndp>1, so I just dropped it). If we want to code it properly, + ! we'd need to make this dim(6,5). !--- tracer handling character(len=32), pointer :: tracer_names(:) !< array of initialized tracers from dynamic core integer :: ntrac !< number of tracers @@ -1948,6 +1947,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: uustar_ocean(:) => null() !< real (kind=kind_phys), pointer :: vdftra(:,:,:) => null() !< real (kind=kind_phys), pointer :: vegf1d(:) => null() !< + real (kind=kind_phys) :: lndp_vgf !< + integer, pointer :: vegtype(:) => null() !< real (kind=kind_phys), pointer :: w_upi(:,:) => null() !< real (kind=kind_phys), pointer :: wcbmax(:) => null() !< @@ -2800,9 +2801,9 @@ subroutine coupling_create (Coupling, IM, Model) Coupling%skebv_wts = clear_val endif - !--- stochastic physics option - if (Model%do_sfcperts) then - allocate (Coupling%sfc_wts (IM,Model%nsfcpert)) + !--- stochastic land perturbation option + if (Model%lndp_type .NE. 0) then + allocate (Coupling%sfc_wts (IM,Model%n_var_lndp)) Coupling%sfc_wts = clear_val endif @@ -3314,15 +3315,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: use_zmtnblck = .false. logical :: do_shum = .false. logical :: do_skeb = .false. - integer :: skeb_npass = 11 - logical :: do_sfcperts = .false. ! mg, sfc-perts - integer :: nsfcpert = 6 ! mg, sfc-perts - real(kind=kind_phys) :: pertz0 = -999. - real(kind=kind_phys) :: pertzt = -999. - real(kind=kind_phys) :: pertshc = -999. - real(kind=kind_phys) :: pertlai = -999. - real(kind=kind_phys) :: pertalb = -999. - real(kind=kind_phys) :: pertvegf = -999. + integer :: skeb_npass = 11 + integer :: lndp_type = 0 + integer :: n_var_lndp = 0 !--- aerosol scavenging factors character(len=20) :: fscav_aero(20) = 'default' @@ -3399,7 +3394,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & do_deep, jcap, & cs_parm, flgmin, cgwf, ccwf, cdmbgwd, sup, ctei_rm, crtrh, & dlqf, rbcr, shoc_parm, psauras, prauras, wminras, & - do_sppt, do_shum, do_skeb, do_sfcperts, & + do_sppt, do_shum, do_skeb, lndp_type, n_var_lndp, & !--- Rayleigh friction prslrd0, ral_ts, ldiag_ugwp, do_ugwp, do_tofd, & ! --- Ferrier-Aligo @@ -4007,21 +4002,15 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%e0fac = e0fac !--- stochastic physics options - ! do_sppt, do_shum, do_skeb and do_sfcperts are namelist variables in group + ! do_sppt, do_shum, do_skeb and lndp_type are namelist variables in group ! physics that are parsed here and then compared in init_stochastic_physics ! to the stochastic physics namelist parametersto ensure consistency. Model%do_sppt = do_sppt Model%use_zmtnblck = use_zmtnblck Model%do_shum = do_shum Model%do_skeb = do_skeb - Model%do_sfcperts = do_sfcperts ! mg, sfc-perts - Model%nsfcpert = nsfcpert ! mg, sfc-perts - Model%pertz0 = pertz0 - Model%pertzt = pertzt - Model%pertshc = pertshc - Model%pertlai = pertlai - Model%pertalb = pertalb - Model%pertvegf = pertvegf + Model%lndp_type = lndp_type + Model%n_var_lndp = n_var_lndp !--- cellular automata options Model%nca = nca @@ -5075,7 +5064,8 @@ subroutine control_print(Model) print *, ' do_sppt : ', Model%do_sppt print *, ' do_shum : ', Model%do_shum print *, ' do_skeb : ', Model%do_skeb - print *, ' do_sfcperts : ', Model%do_sfcperts + print *, ' lndp_type : ', Model%lndp_type + print *, ' n_var_lndp : ', Model%n_var_lndp print *, ' ' print *, 'cellular automata' print *, ' nca : ', Model%nca @@ -7015,6 +7005,7 @@ subroutine interstitial_phys_reset (Interstitial, Model) Interstitial%uustar_ocean = huge Interstitial%vdftra = clear_val Interstitial%vegf1d = clear_val + Interstitial%lndp_vgf = clear_val Interstitial%vegtype = 0 Interstitial%wcbmax = clear_val Interstitial%weasd_ice = huge diff --git a/gfsphysics/GFS_layer/GFS_typedefs.meta b/gfsphysics/GFS_layer/GFS_typedefs.meta index 4c5d75267..dab6a5c17 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.meta +++ b/gfsphysics/GFS_layer/GFS_typedefs.meta @@ -1900,7 +1900,7 @@ standard_name = weights_for_stochastic_surface_physics_perturbation long_name = weights for stochastic surface physics perturbation units = none - dimensions = (horizontal_dimension,number_of_surface_perturbations) + dimensions = (horizontal_dimension,number_of_land_surface_variables_perturbed) type = real kind = kind_phys [dqdti] @@ -3691,8 +3691,8 @@ type = real kind = kind_phys [do_sppt] - standard_name = flag_for_stochastic_surface_physics_perturbations - long_name = flag for stochastic surface physics perturbations + standard_name = flag_for_stochastic_physics_perturbations + long_name = flag for stochastic physics perturbations units = flag dimensions = () type = logical @@ -3714,60 +3714,32 @@ units = flag dimensions = () type = logical -[do_sfcperts] - standard_name = flag_for_stochastic_surface_perturbations - long_name = flag for stochastic surface perturbations option - units = flag +[lndp_type] + standard_name = index_for_stochastic_land_surface_perturbation_type + long_name = index for stochastic land surface perturbations type + units = index dimensions = () - type = logical -[nsfcpert] - standard_name = number_of_surface_perturbations - long_name = number of surface perturbations + type = integer +[n_var_lndp] + standard_name = number_of_land_surface_variables_perturbed + long_name = number of land surface variables perturbed units = count dimensions = () type = integer -[pertz0] - standard_name = magnitude_of_perturbation_of_momentum_roughness_length - long_name = magnitude of perturbation of momentum roughness length - units = frac - dimensions = (5) - type = real - kind = kind_phys -[pertzt] - standard_name = magnitude_of_perturbation_of_heat_to_momentum_roughness_length_ratio - long_name = magnitude of perturbation of heat to momentum roughness length ratio - units = frac - dimensions = (5) - type = real - kind = kind_phys -[pertshc] - standard_name = magnitude_of_perturbation_of_soil_type_b_parameter - long_name = magnitude of perturbation of soil type b parameter - units = frac - dimensions = (5) - type = real - kind = kind_phys -[pertlai] - standard_name = magnitude_of_perturbation_of_leaf_area_index - long_name = magnitude of perturbation of leaf area index - units = frac - dimensions = (5) - type = real - kind = kind_phys -[pertalb] - standard_name = magnitude_of_surface_albedo_perturbation - long_name = magnitude of surface albedo perturbation - units = frac - dimensions = (5) - type = real - kind = kind_phys -[pertvegf] - standard_name = magnitude_of_perturbation_of_vegetation_fraction - long_name = magnitude of perturbation of vegetation fraction - units = frac - dimensions = (5) +[lndp_prt_list] + standard_name =magnitude_of_perturbations_for_landperts + long_name = magnitude of perturbations for landperts + units = variable + dimensions = (number_of_land_surface_variables_perturbed) type = real kind = kind_phys +[lndp_var_list] + standard_name = variables_to_be_perturbed_for_landperts + long_name = variables to be perturbed for landperts + units = none + dimensions = (number_of_land_surface_variables_perturbed) + type = character + kind = len=3 [ntrac] standard_name = number_of_tracers long_name = number of tracers @@ -8591,6 +8563,13 @@ dimensions = (horizontal_dimension,vertical_dimension,number_of_vertical_diffusion_tracers) type = real kind = kind_phys +[lndp_vgf] + standard_name = magnitude_of_perturbation_of_vegetation_fraction + long_name = magnitude of perturbation of vegetation fraction + units = frac + dimensions = () + type = real + kind = kind_phys [vegf1d] standard_name = perturbation_of_vegetation_fraction long_name = perturbation of vegetation fraction diff --git a/gfsphysics/physics/GFS_debug.F90 b/gfsphysics/physics/GFS_debug.F90 index 2887d6e64..e1953a9ff 100644 --- a/gfsphysics/physics/GFS_debug.F90 +++ b/gfsphysics/physics/GFS_debug.F90 @@ -468,7 +468,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts ) call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts ) end if - if (Model%do_sfcperts) then + if (Model%lndp_type .NE. 0) then call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts' , Coupling%sfc_wts ) end if if (Model%do_ca) then diff --git a/gfsphysics/physics/radiation_surface.f b/gfsphysics/physics/radiation_surface.f index 99f0ebc2f..9ae258a0c 100644 --- a/gfsphysics/physics/radiation_surface.f +++ b/gfsphysics/physics/radiation_surface.f @@ -382,7 +382,7 @@ subroutine setalb & & slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif, & & alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, & & sncovr, snoalb, albPpert ! sfc-perts, mgehne - real (kind=kind_phys), dimension(5), intent(in) :: pertalb ! sfc-perts, mgehne + real (kind=kind_phys), intent(in) :: pertalb ! sfc-perts, mgehne ! --- outputs real (kind=kind_phys), dimension(IMAX,NF_ALBD), intent(out) :: & @@ -620,12 +620,12 @@ subroutine setalb & ! sfc-perts, mgehne *** ! perturb all 4 kinds of surface albedo, sfcalb(:,1:4) - if (pertalb(1) > 0.0) then + if (pertalb>0.0) then do i = 1, imax do kk=1, 4 ! compute beta distribution parameters for all 4 albedos m = sfcalb(i,kk) - s = pertalb(1)*m*(1.-m) + s = pertalb*m*(1.-m) alpha = m*m*(1.-m)/(s*s)-m beta = alpha*(1.-m)/m ! compute beta distribution value corresponding diff --git a/gfsphysics/physics/sfc_drv.f b/gfsphysics/physics/sfc_drv.f index 84b4b84d5..80e081909 100644 --- a/gfsphysics/physics/sfc_drv.f +++ b/gfsphysics/physics/sfc_drv.f @@ -182,7 +182,7 @@ subroutine sfc_drv & ! --- input: integer, intent(in) :: im, km, isot, ivegsrc - real (kind=kind_phys), dimension(5), intent(in) :: pertvegf + real (kind=kind_phys), intent(in) :: pertvegf integer, dimension(im), intent(in) :: soiltyp, vegtype, slopetyp @@ -383,10 +383,10 @@ subroutine sfc_drv & ! perturb vegetation fraction that goes into sflx, use the same ! perturbation strategy as for albedo (percentile matching) vegfp = vegfpert(i) ! sfc-perts, mgehne - if (pertvegf(1) > zero) then + if (pertvegf > zero) then ! compute beta distribution parameters for vegetation fraction mv = shdfac - sv = pertvegf(1)*mv*(1.-mv) + sv = pertvegf*mv*(one-mv) alphav = mv*mv*(one-mv)/(sv*sv)-mv betav = alphav*(one-mv)/mv ! compute beta distribution value corresponding diff --git a/stochastic_physics/makefile b/stochastic_physics/makefile index eb721c8c4..c841571a4 100644 --- a/stochastic_physics/makefile +++ b/stochastic_physics/makefile @@ -18,7 +18,7 @@ $(info $$ESMF_INC is [${ESMF_INC}]) LIBRARY = libstochastic_physics_wrapper.a -FFLAGS += -I$(FMS_DIR) -I ../../stochastic_physics -I../ccpp/physics -I../atmos_cubed_sphere +FFLAGS += -I$(FMS_DIR) -I ../../stochastic_physics -I../ccpp/physics -I../ccpp/build/physics -I../atmos_cubed_sphere SRCS_f = diff --git a/stochastic_physics/stochastic_physics_wrapper.F90 b/stochastic_physics/stochastic_physics_wrapper.F90 index 2fb8cafd1..d4803c639 100644 --- a/stochastic_physics/stochastic_physics_wrapper.F90 +++ b/stochastic_physics/stochastic_physics_wrapper.F90 @@ -13,6 +13,12 @@ module stochastic_physics_wrapper_mod real(kind=kind_phys), dimension(:,:,:), allocatable, save :: skebv_wts real(kind=kind_phys), dimension(:,:,:), allocatable, save :: sfc_wts + real(kind=kind_phys), dimension(:,:,:), allocatable, save :: smc + real(kind=kind_phys), dimension(:,:,:), allocatable, save :: stc + real(kind=kind_phys), dimension(:,:,:), allocatable, save :: slc + real(kind=kind_phys), dimension(:,:), allocatable, save :: vfrac + real(kind=kind_phys), dimension(:,:), allocatable, save :: stype + ! For cellular automata real(kind=kind_phys), dimension(:,:,:), allocatable, save :: ugrs real(kind=kind_phys), dimension(:,:,:), allocatable, save :: qgrs @@ -37,7 +43,7 @@ module stochastic_physics_wrapper_mod !------------------------------- ! CCPP step !------------------------------- - subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) + subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) #ifdef OPENMP use omp_lib @@ -49,47 +55,56 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) use atmosphere_mod, only: Atm, mygrid use stochastic_physics, only: init_stochastic_physics, run_stochastic_physics - use stochastic_physics_sfc, only: run_stochastic_physics_sfc use cellular_automata_global_mod, only: cellular_automata_global use cellular_automata_sgs_mod, only: cellular_automata_sgs + use lndp_apply_perts_mod, only: lndp_apply_perts + use namelist_soilveg, only: maxsmc implicit none type(GFS_control_type), intent(inout) :: GFS_Control type(GFS_data_type), intent(inout) :: GFS_Data(:) type(block_control_type), intent(inout) :: Atm_block + integer, intent(out) :: ierr integer :: nthreads, nb + logical :: param_update_flag #ifdef OPENMP nthreads = omp_get_max_threads() #else nthreads = 1 #endif + ierr = 0 ! Initialize initalize_stochastic_physics: if (GFS_Control%kdt==0) then - if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. GFS_Control%do_sfcperts) then + if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .GT. 0) ) then ! Initialize stochastic physics call init_stochastic_physics(GFS_Control%levs, GFS_Control%blksz, GFS_Control%dtp, & GFS_Control%input_nml_file, GFS_Control%fn_nml, GFS_Control%nlunit, GFS_Control%do_sppt, GFS_Control%do_shum, & - GFS_Control%do_skeb, GFS_Control%do_sfcperts, GFS_Control%use_zmtnblck, GFS_Control%skeb_npass, GFS_Control%nsfcpert, & - GFS_Control%pertz0, GFS_Control%pertzt, GFS_Control%pertshc, GFS_Control%pertlai, GFS_Control%pertalb, GFS_Control%pertvegf, & - GFS_Control%ak, GFS_Control%bk, nthreads, GFS_Control%master, GFS_Control%communicator) + GFS_Control%do_skeb, GFS_Control%lndp_type, GFS_Control%n_var_lndp, GFS_Control%use_zmtnblck, GFS_Control%skeb_npass, & + GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list, & + GFS_Control%ak, GFS_Control%bk, nthreads, GFS_Control%master, GFS_Control%communicator, ierr) + if (ierr/=0) then + write(6,*) 'call to init_stochastic_physics failed' + return + endif end if - ! Get land surface perturbations here (move to "else" block below if wanting to update each time-step) - if (GFS_Control%do_sfcperts) then + if ( GFS_Control%lndp_type .EQ. 1 ) then ! this scheme sets perts once ! Copy blocked data into contiguous arrays; no need to copy sfc_wts in (intent out) allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz))) allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz))) - allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs)) + allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%n_var_lndp)) do nb=1,Atm_block%nblks xlat(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlat(:) xlon(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlon(:) end do - call run_stochastic_physics_sfc(GFS_Control%blksz, xlat=xlat, xlon=xlon, sfc_wts=sfc_wts) + call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%phour, GFS_Control%blksz, xlat=xlat, xlon=xlon, & + sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, sfc_wts=sfc_wts, & + nthreads=nthreads) ! Copy contiguous data back; no need to copy xlat/xlon, these are intent(in) - just deallocate do nb=1,Atm_block%nblks GFS_Data(nb)%Coupling%sfc_wts(:,:) = sfc_wts(nb,1:GFS_Control%blksz(nb),:) @@ -97,8 +112,7 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) deallocate(xlat) deallocate(xlon) deallocate(sfc_wts) - end if - + end if ! Consistency check for cellular automata if(GFS_Control%do_ca)then ! DH* The current implementation of cellular_automata assumes that all blocksizes are the @@ -111,7 +125,7 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) else initalize_stochastic_physics - if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb) then + if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .EQ. 2) ) then ! Copy blocked data into contiguous arrays; no need to copy weights in (intent(out)) allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz))) allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz))) @@ -129,11 +143,14 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) allocate(skebu_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs)) allocate(skebv_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs)) end if + if ( GFS_Control%lndp_type .EQ. 2 ) then ! this scheme updates through forecast + allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%n_var_lndp)) + end if + call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%phour, GFS_Control%blksz, xlat=xlat, xlon=xlon, & - sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, nthreads=nthreads) + sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, sfc_wts=sfc_wts, & + nthreads=nthreads) ! Copy contiguous data back; no need to copy xlat/xlon, these are intent(in) - just deallocate - deallocate(xlat) - deallocate(xlon) if (GFS_Control%do_sppt) then do nb=1,Atm_block%nblks GFS_Data(nb)%Coupling%sppt_wts(:,:) = sppt_wts(nb,1:GFS_Control%blksz(nb),:) @@ -154,6 +171,52 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) deallocate(skebu_wts) deallocate(skebv_wts) end if + if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme + do nb=1,Atm_block%nblks + GFS_Data(nb)%Coupling%sfc_wts(:,:) = sfc_wts(nb,1:GFS_Control%blksz(nb),:) + end do + + allocate(smc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil)) + allocate(slc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil)) + allocate(stc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil)) + allocate(stype(1:Atm_block%nblks,maxval(GFS_Control%blksz))) + allocate(vfrac(1:Atm_block%nblks,maxval(GFS_Control%blksz))) + do nb=1,Atm_block%nblks + stype(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%stype(:) + smc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%smc(:,:) + slc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%slc(:,:) + stc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%stc(:,:) + vfrac(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%vfrac(:) + end do + + ! determine whether land paramaters have been over-written + if (mod(GFS_Control%kdt,GFS_Control%nscyc) == 1) then ! logic copied from GFS_driver + param_update_flag = .true. + else + param_update_flag = .false. + endif + call lndp_apply_perts( GFS_Control%blksz, GFS_Control%lsm, GFS_Control%lsoil, GFS_Control%dtf, & + GFS_Control%n_var_lndp, GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list, & + sfc_wts, xlon, xlat, stype, maxsmc,param_update_flag, smc, slc,stc, vfrac, ierr) + if (ierr/=0) then + write(6,*) 'call to GFS_apply_lndp failed' + return + endif + deallocate(stype) + deallocate(sfc_wts) + do nb=1,Atm_block%nblks + GFS_Data(nb)%Sfcprop%smc(:,:) = smc(nb,1:GFS_Control%blksz(nb),:) + GFS_Data(nb)%Sfcprop%slc(:,:) = slc(nb,1:GFS_Control%blksz(nb),:) + GFS_Data(nb)%Sfcprop%stc(:,:) = stc(nb,1:GFS_Control%blksz(nb),:) + GFS_Data(nb)%Sfcprop%vfrac(:) = vfrac(nb,1:GFS_Control%blksz(nb)) + enddo + deallocate(smc) + deallocate(slc) + deallocate(stc) + deallocate(vfrac) + endif ! lndp block + deallocate(xlat) + deallocate(xlon) end if endif initalize_stochastic_physics