Skip to content

Commit

Permalink
Land stochastic perturbations (wrapper PR for NCAR#65) (NCAR#68)
Browse files Browse the repository at this point in the history 
* Move initialization of stochastic physics after the physics
initialization in CCPP.
* Add albedo variables to land perturbations with lndp_type=2 option. Change to accommodate soil perturbations with RUC LSM.

* Max/min soil moisture variables are introduced via GFS_Control_type
variables instead of through the use of namelist_soilveg*. This is a
more flexible way for different LSMs.

* Added pores and resid variables for max/min soil moisture to GFS_typedefs.f90.

* Remove tracer_sanitizer from all suites and from CCPP prebuild config

* Add namelist option to apply land surface perturbations at every time step, clean up stochastic_physics/stochastic_physics_wrapper.F90

Co-authored-by: tanyasmirnova <[email protected]>
  • Loading branch information
@climbfuji @tanyasmirnova
climbfuji and tanyasmirnova authored Jan 9, 2021
1 parent b1e98cf commit 6ecee94
Show file tree
Hide file tree
Showing 8 changed files with 139 additions and 30 deletions.
9 changes: 6 additions & 3 deletions atmos_model.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -627,9 +627,9 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
IPD_Interstitial, commglobal, mpp_npes(), Init_parm)

!--- 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, ierr)
! if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
!
#else
call IPD_initialize (IPD_Control, IPD_Data, IPD_Diag, IPD_Restart, Init_parm)
#endif
Expand Down Expand Up @@ -684,6 +684,9 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
! Initialize the CCPP physics
call CCPP_step (step="physics_init", nblks=Atm_block%nblks, ierr=ierr)
if (ierr/=0) call mpp_error(FATAL, 'Call to CCPP physics_init step failed')
!--- 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')
#endif

!--- set the initial diagnostic timestamp
Expand Down
1 change: 0 additions & 1 deletion ccpp/config/ccpp_prebuild_config.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,6 @@
'ccpp/physics/physics/ozphys_2015.f',
'ccpp/physics/physics/precpd.f',
'ccpp/physics/physics/phys_tend.F90',
'ccpp/physics/physics/tracer_sanitizer.F90',
'ccpp/physics/physics/radlw_main.F90',
'ccpp/physics/physics/radsw_main.F90',
'ccpp/physics/physics/rascnv.F90',
Expand Down
2 changes: 1 addition & 1 deletion ccpp/physics
Submodule physics updated 8 files
+10 −2 physics/sfc_drv.f
+16 −0 physics/sfc_drv.meta
+7 −2 physics/sfc_drv_ruc.F90
+16 −0 physics/sfc_drv_ruc.meta
+10 −2 physics/sfc_noahmp_drv.f
+16 −0 physics/sfc_noahmp_drv.meta
+0 −113 physics/tracer_sanitizer.F90
+0 −124 physics/tracer_sanitizer.meta
1 change: 0 additions & 1 deletion ccpp/suites/suite_FV3_GSD_noah.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,6 @@
<scheme>mp_thompson_post</scheme>
<scheme>GFS_MP_generic_post</scheme>
<scheme>cu_gf_driver_post</scheme>
<!-- <scheme>tracer_sanitizer</scheme> -->
<scheme>maximum_hourly_diagnostics</scheme>
</subcycle>
</group>
Expand Down
1 change: 0 additions & 1 deletion ccpp/suites/suite_FV3_GSD_v0.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,6 @@
<scheme>mp_thompson_post</scheme>
<scheme>GFS_MP_generic_post</scheme>
<scheme>cu_gf_driver_post</scheme>
<!-- <scheme>tracer_sanitizer</scheme> -->
<scheme>maximum_hourly_diagnostics</scheme>
</subcycle>
</group>
Expand Down
26 changes: 21 additions & 5 deletions gfsphysics/GFS_layer/GFS_typedefs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -811,7 +811,9 @@ module GFS_typedefs
integer :: lsoil_lsm !< number of soil layers internal to land surface model
integer :: lsnow_lsm !< maximum number of snow layers internal to land surface model
integer :: lsnow_lsm_lbound!< lower bound for snow arrays, depending on lsnow_lsm
real(kind=kind_phys), pointer :: zs(:) => null() !< depth of soil levels for land surface model
real(kind=kind_phys), pointer :: zs(:) => null() !< depth of soil levels for land surface model
real(kind=kind_phys), pointer :: pores(:) => null() !< max soil moisture for a given soil type for land surface model
real(kind=kind_phys), pointer :: resid(:) => null() !< min soil moisture for a given soil type for land surface model
logical :: rdlai !< read LAI from input file (for RUC LSM or NOAH LSM WRFv4)
logical :: ua_phys !< flag for using University of Arizona? extension to NOAH LSM WRFv4
logical :: usemonalb !< flag to read surface diffused shortwave albedo from input file for NOAH LSM WRFv4
Expand Down Expand Up @@ -1110,6 +1112,8 @@ module GFS_typedefs
integer :: skeb_npass
integer :: lndp_type
integer :: n_var_lndp
logical :: lndp_each_step ! flag to indicate that land perturbations are applied at every time step,
! otherwise they are applied only after gcycle is run
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
Expand Down Expand Up @@ -3469,8 +3473,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
logical :: do_shum = .false.
logical :: do_skeb = .false.
integer :: skeb_npass = 11
integer :: lndp_type = 0
integer :: n_var_lndp = 0
integer :: lndp_type = 0
integer :: n_var_lndp = 0
logical :: lndp_each_step = .false.

!--- aerosol scavenging factors
character(len=20) :: fscav_aero(20) = 'default'
Expand Down Expand Up @@ -3555,7 +3560,8 @@ 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, &
lndp_type, n_var_lndp, lndp_each_step, &
!--- Rayleigh friction
prslrd0, ral_ts, ldiag_ugwp, do_ugwp, do_tofd, &
! --- Ferrier-Aligo
Expand Down Expand Up @@ -3981,6 +3987,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
end if
! Set number of ice model layers
Model%kice = kice

! Allocate variable for min/max soil moisture for a given soil type
allocate (Model%pores(30))
allocate (Model%resid(30))
Model%pores = clear_val
Model%resid = clear_val
!
if (lsnow_lsm /= 3) then
write(0,*) 'Logic error: NoahMP expects the maximum number of snow layers to be exactly 3 (see sfc_noahmp_drv.f)'
Expand Down Expand Up @@ -4216,6 +4228,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%do_skeb = do_skeb
Model%lndp_type = lndp_type
Model%n_var_lndp = n_var_lndp
Model%lndp_each_step = lndp_each_step

!--- cellular automata options
Model%nca = nca
Expand Down Expand Up @@ -5159,6 +5172,8 @@ subroutine control_print(Model)
print *, ' aoasis : ', Model%aoasis
print *, ' fasdas : ', Model%fasdas
print *, ' kice : ', Model%kice
print *, ' shape(pores) : ', shape(Model%pores)
print *, ' shape(resid) : ', shape(Model%resid)
#endif
print *, ' ivegsrc : ', Model%ivegsrc
print *, ' isot : ', Model%isot
Expand Down Expand Up @@ -5316,7 +5331,8 @@ subroutine control_print(Model)
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 *, ' n_var_lndp : ', Model%n_var_lndp
print *, ' lndp_each_step : ', Model%lndp_each_step
print *, ' '
print *, 'cellular automata'
print *, ' nca : ', Model%nca
Expand Down
14 changes: 14 additions & 0 deletions gfsphysics/GFS_layer/GFS_typedefs.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3285,6 +3285,20 @@
type = real
kind = kind_phys
active = (flag_for_land_surface_scheme == flag_for_ruc_land_surface_scheme)
[pores]
standard_name = maximum_soil_moisture_content_for_land_surface_model
long_name = maximum soil moisture for a given soil type for land surface model
units = m
dimensions = (30)
type = real
kind = kind_phys
[resid]
standard_name = minimum_soil_moisture_content_for_land_surface_model
long_name = minimum soil moisture for a given soil type for land surface model
units = m
dimensions = (30)
type = real
kind = kind_phys
[rdlai]
standard_name = flag_for_reading_leaf_area_index_from_input
long_name = flag for reading leaf area index from initial conditions
Expand Down
115 changes: 97 additions & 18 deletions stochastic_physics/stochastic_physics_wrapper.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,23 @@ module stochastic_physics_wrapper_mod
real(kind=kind_phys), dimension(:,:,:), allocatable, save :: skebv_wts
real(kind=kind_phys), dimension(:,:,:), allocatable, save :: sfc_wts

integer, save :: lsoil = -999
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
!albedo
real(kind=kind_phys), dimension(:,:), allocatable, save :: snoalb
real(kind=kind_phys), dimension(:,:), allocatable, save :: alvsf
real(kind=kind_phys), dimension(:,:), allocatable, save :: alnsf
real(kind=kind_phys), dimension(:,:), allocatable, save :: alvwf
real(kind=kind_phys), dimension(:,:), allocatable, save :: alnwf
real(kind=kind_phys), dimension(:,:), allocatable, save :: facsf
real(kind=kind_phys), dimension(:,:), allocatable, save :: facwf
!emissivity
real(kind=kind_phys), dimension(:,:), allocatable, save :: semis

real(kind=kind_phys), dimension(:,:), allocatable, save :: stype

! For cellular automata
Expand Down Expand Up @@ -58,7 +71,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
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

Expand Down Expand Up @@ -108,11 +120,24 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%n_var_lndp))
end if
if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
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))
if (GFS_Control%lsm == GFS_Control%lsm_noah) then
lsoil = GFS_Control%lsoil
elseif (GFS_Control%lsm == GFS_Control%lsm_ruc) then
lsoil = GFS_Control%lsoil_lsm
endif
allocate(smc(1:Atm_block%nblks,maxval(GFS_Control%blksz),lsoil))
allocate(slc(1:Atm_block%nblks,maxval(GFS_Control%blksz),lsoil))
allocate(stc(1:Atm_block%nblks,maxval(GFS_Control%blksz),lsoil))
allocate(stype(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(vfrac(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(snoalb(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(alvsf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(alnsf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(alvwf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(alnwf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(facsf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(facwf(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(semis(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
endif

do nb=1,Atm_block%nblks
Expand Down Expand Up @@ -171,31 +196,76 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)

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(:)
snoalb(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%snoalb(:)
alvsf(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%alvsf(:)
alnsf(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%alnsf(:)
alvwf(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%alvwf(:)
alnwf(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%alnwf(:)
facsf(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%facsf(:)
facwf(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%facwf(:)
semis(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Radtend%semis(:)
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.
if (GFS_Control%lsm == GFS_Control%lsm_noah) then
do nb=1,Atm_block%nblks
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(:,:)
end do
elseif (GFS_Control%lsm == GFS_Control%lsm_ruc) then
do nb=1,Atm_block%nblks
smc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%smois(:,:)
slc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%sh2o(:,:)
stc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%tslb(:,:)
end do
endif

! determine whether land paramaters have been over-written to
! trigger applying perturbations (logic copied from GFS_driver),
! or if perturbations should be applied at every time step
if (mod(GFS_Control%kdt,GFS_Control%nscyc) == 1 ) then
param_update_flag = .true.
else
param_update_flag = .false.
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)

call lndp_apply_perts(GFS_Control%blksz, GFS_Control%lsm, GFS_Control%lsm_noah, GFS_Control%lsm_ruc, lsoil, &
GFS_Control%dtf, GFS_Control%kdt, GFS_Control%lndp_each_step, &
GFS_Control%n_var_lndp, GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list, &
sfc_wts, xlon, xlat, stype, GFS_Control%pores, GFS_Control%resid,param_update_flag, &
smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, snoalb, semis, ierr)
if (ierr/=0) then
write(6,*) 'call to GFS_apply_lndp failed'
return
endif

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))
GFS_Data(nb)%Sfcprop%vfrac(:) = vfrac(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%snoalb(:) = snoalb(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%alvsf(:) = alvsf(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%alnsf(:) = alnsf(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%alvwf(:) = alvwf(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%alnwf(:) = alnwf(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%facsf(:) = facsf(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Sfcprop%facwf(:) = facwf(nb,1:GFS_Control%blksz(nb))
GFS_Data(nb)%Radtend%semis(:) = semis(nb,1:GFS_Control%blksz(nb))
enddo

if (GFS_Control%lsm == GFS_Control%lsm_noah) then
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),:)
enddo
elseif (GFS_Control%lsm == GFS_Control%lsm_ruc) then
do nb=1,Atm_block%nblks
GFS_Data(nb)%Sfcprop%smois(:,:) = smc(nb,1:GFS_Control%blksz(nb),:)
GFS_Data(nb)%Sfcprop%sh2o(:,:) = slc(nb,1:GFS_Control%blksz(nb),:)
GFS_Data(nb)%Sfcprop%tslb(:,:) = stc(nb,1:GFS_Control%blksz(nb),:)
enddo
endif

endif ! lndp block
end if

Expand Down Expand Up @@ -313,6 +383,7 @@ subroutine stochastic_physics_wrapper_end (GFS_Control)
if (allocated(skebv_wts)) deallocate(skebv_wts)
end if
if ( GFS_Control%lndp_type .EQ. 2 ) then ! this scheme updates through forecast
lsoil = -999
if (allocated(sfc_wts)) deallocate(sfc_wts)
end if
if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
Expand All @@ -321,6 +392,14 @@ subroutine stochastic_physics_wrapper_end (GFS_Control)
if (allocated(stc)) deallocate(stc)
if (allocated(stype)) deallocate(stype)
if (allocated(vfrac)) deallocate(vfrac)
if (allocated(snoalb)) deallocate(snoalb)
if (allocated(alvsf)) deallocate(alvsf)
if (allocated(alnsf)) deallocate(alnsf)
if (allocated(alvwf)) deallocate(alvwf)
if (allocated(alnwf)) deallocate(alnwf)
if (allocated(facsf)) deallocate(facsf)
if (allocated(facwf)) deallocate(facwf)
if (allocated(semis)) deallocate(semis)
endif
call finalize_stochastic_physics()
endif
Expand Down

0 comments on commit 6ecee94

Please sign in to comment.