diff --git a/.gitmodules b/.gitmodules index d253f6966..789f67889 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,5 +8,5 @@ branch = master [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/NCAR/ccpp-physics - branch = master + url = https://github.com/NCAR/ccpp-physics.git + branch = master diff --git a/atmos_model.F90 b/atmos_model.F90 index 35b2c4ceb..900c9143e 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -297,9 +297,7 @@ subroutine update_atmos_radiation_physics (Atmos) #endif !--- call stochastic physics pattern generation / cellular automata - call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block, ierr) - if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') - + call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block) !--- if coupled, assign coupled fields @@ -630,8 +628,7 @@ 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, ierr) - if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') + call stochastic_physics_wrapper(IPD_Control, IPD_Data, Atm_block) Atmos%Diag => IPD_Diag diff --git a/ccpp/framework b/ccpp/framework index f5d4cd2bf..209f1c92d 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit f5d4cd2bf7752ebf1e4ed16dfdfae71dbfabfb76 +Subproject commit 209f1c92d99b7d4cc63e0d41c652fcfd730bd9fa diff --git a/ccpp/physics b/ccpp/physics index 4c17ff716..8617587ed 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 4c17ff716a92d8ac0261d6cc2365bbd7a752b74a +Subproject commit 8617587edb95aa097b7bbc2735990393bc6d9b90 diff --git a/gfsphysics/CMakeLists.txt b/gfsphysics/CMakeLists.txt index 2c5040588..34faf6f4c 100644 --- a/gfsphysics/CMakeLists.txt +++ b/gfsphysics/CMakeLists.txt @@ -17,6 +17,7 @@ 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 e32ec4c26..4abde43d0 100644 --- a/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_driver.F90 @@ -633,7 +633,6 @@ 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 899955f03..02eb00e00 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -660,7 +660,6 @@ 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). @@ -929,28 +928,34 @@ subroutine GFS_physics_driver & ! alb1d(i) = zero vegf1d(i) = zero enddo - 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 + 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 endif !*## CCPP ## ! @@ -1851,7 +1856,6 @@ 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, & @@ -1863,8 +1867,7 @@ 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,lndp_vgf, & - + bexp1d, xlai1d, vegf1d, Model%pertvegf, & ! --- 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 da5078f7b..ebec30c4d 100644 --- a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 @@ -1221,7 +1221,6 @@ 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 @@ -1847,18 +1846,14 @@ subroutine GFS_radiation_driver & ! --- scale random patterns for surface perturbations with ! perturbation size ! --- turn vegetation fraction pattern into percentile pattern - 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 + 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 endif + endif ! mg, sfc-perts !*## CCPP ## @@ -1875,7 +1870,7 @@ subroutine GFS_radiation_driver & Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf, & Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice, & Sfcprop%tisfc, IM, & - alb1d, lndp_alb, & ! mg, sfc-perts + alb1d, Model%pertalb, & ! 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 4d3140d24..3c1252b1e 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -517,6 +517,7 @@ 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 @@ -1044,14 +1045,14 @@ module GFS_typedefs logical :: do_shum logical :: do_skeb integer :: skeb_npass - 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). + 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 !--- tracer handling character(len=32), pointer :: tracer_names(:) !< array of initialized tracers from dynamic core integer :: ntrac !< number of tracers @@ -1947,8 +1948,6 @@ 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() !< @@ -2801,9 +2800,9 @@ subroutine coupling_create (Coupling, IM, Model) Coupling%skebv_wts = clear_val endif - !--- stochastic land perturbation option - if (Model%lndp_type .NE. 0) then - allocate (Coupling%sfc_wts (IM,Model%n_var_lndp)) + !--- stochastic physics option + if (Model%do_sfcperts) then + allocate (Coupling%sfc_wts (IM,Model%nsfcpert)) Coupling%sfc_wts = clear_val endif @@ -3315,9 +3314,15 @@ 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 - integer :: lndp_type = 0 - integer :: n_var_lndp = 0 + 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. !--- aerosol scavenging factors character(len=20) :: fscav_aero(20) = 'default' @@ -3394,7 +3399,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, lndp_type, n_var_lndp, & + do_sppt, do_shum, do_skeb, do_sfcperts, & !--- Rayleigh friction prslrd0, ral_ts, ldiag_ugwp, do_ugwp, do_tofd, & ! --- Ferrier-Aligo @@ -4002,15 +4007,21 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%e0fac = e0fac !--- stochastic physics options - ! do_sppt, do_shum, do_skeb and lndp_type are namelist variables in group + ! do_sppt, do_shum, do_skeb and do_sfcperts 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%lndp_type = lndp_type - Model%n_var_lndp = n_var_lndp + 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 !--- cellular automata options Model%nca = nca @@ -5064,8 +5075,7 @@ subroutine control_print(Model) print *, ' do_sppt : ', Model%do_sppt print *, ' do_shum : ', Model%do_shum print *, ' do_skeb : ', Model%do_skeb - print *, ' lndp_type : ', Model%lndp_type - print *, ' n_var_lndp : ', Model%n_var_lndp + print *, ' do_sfcperts : ', Model%do_sfcperts print *, ' ' print *, 'cellular automata' print *, ' nca : ', Model%nca @@ -7005,7 +7015,6 @@ 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 519f532f5..28c1963a5 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.meta +++ b/gfsphysics/GFS_layer/GFS_typedefs.meta @@ -2038,7 +2038,7 @@ standard_name = weights_for_stochastic_surface_physics_perturbation long_name = weights for stochastic surface physics perturbation units = none - dimensions = (horizontal_dimension,number_of_land_surface_variables_perturbed) + dimensions = (horizontal_dimension,number_of_surface_perturbations) type = real kind = kind_phys active = (index_for_stochastic_land_surface_perturbation_type .ne. 0) @@ -3834,8 +3834,8 @@ type = real kind = kind_phys [do_sppt] - standard_name = flag_for_stochastic_physics_perturbations - long_name = flag for stochastic physics perturbations + standard_name = flag_for_stochastic_surface_physics_perturbations + long_name = flag for stochastic surface physics perturbations units = flag dimensions = () type = logical @@ -3857,32 +3857,60 @@ units = flag dimensions = () type = logical -[lndp_type] - standard_name = index_for_stochastic_land_surface_perturbation_type - long_name = index for stochastic land surface perturbations type - units = index +[do_sfcperts] + standard_name = flag_for_stochastic_surface_perturbations + long_name = flag for stochastic surface perturbations option + units = flag dimensions = () - type = integer -[n_var_lndp] - standard_name = number_of_land_surface_variables_perturbed - long_name = number of land surface variables perturbed + type = logical +[nsfcpert] + standard_name = number_of_surface_perturbations + long_name = number of surface perturbations units = count dimensions = () type = integer -[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) +[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) 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 @@ -8825,13 +8853,6 @@ 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 e1953a9ff..2887d6e64 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%lndp_type .NE. 0) then + if (Model%do_sfcperts) 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 9ae258a0c..99f0ebc2f 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), intent(in) :: pertalb ! sfc-perts, mgehne + real (kind=kind_phys), dimension(5), 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>0.0) then + if (pertalb(1) > 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*m*(1.-m) + s = pertalb(1)*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 80e081909..84b4b84d5 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), intent(in) :: pertvegf + real (kind=kind_phys), dimension(5), 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 > zero) then + if (pertvegf(1) > zero) then ! compute beta distribution parameters for vegetation fraction mv = shdfac - sv = pertvegf*mv*(one-mv) + sv = pertvegf(1)*mv*(1.-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 c841571a4..eb721c8c4 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../ccpp/build/physics -I../atmos_cubed_sphere +FFLAGS += -I$(FMS_DIR) -I ../../stochastic_physics -I../ccpp/physics -I../atmos_cubed_sphere SRCS_f = diff --git a/stochastic_physics/stochastic_physics_wrapper.F90 b/stochastic_physics/stochastic_physics_wrapper.F90 index d4803c639..2fb8cafd1 100644 --- a/stochastic_physics/stochastic_physics_wrapper.F90 +++ b/stochastic_physics/stochastic_physics_wrapper.F90 @@ -13,12 +13,6 @@ 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 @@ -43,7 +37,7 @@ module stochastic_physics_wrapper_mod !------------------------------- ! CCPP step !------------------------------- - subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) + subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block) #ifdef OPENMP use omp_lib @@ -55,56 +49,47 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) 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%lndp_type .GT. 0) ) then + if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. GFS_Control%do_sfcperts) 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%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 + 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) end if - if ( GFS_Control%lndp_type .EQ. 1 ) then ! this scheme sets perts once + ! Get land surface perturbations here (move to "else" block below if wanting to update each time-step) + if (GFS_Control%do_sfcperts) then ! 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),GFS_Control%n_var_lndp)) + allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs)) 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(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) + call run_stochastic_physics_sfc(GFS_Control%blksz, xlat=xlat, xlon=xlon, sfc_wts=sfc_wts) ! 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),:) @@ -112,7 +97,8 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) 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 @@ -125,7 +111,7 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) else initalize_stochastic_physics - if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .EQ. 2) ) then + if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb) 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))) @@ -143,14 +129,11 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) 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, sfc_wts=sfc_wts, & - nthreads=nthreads) + sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_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),:) @@ -171,52 +154,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) 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