-
Notifications
You must be signed in to change notification settings - Fork 161
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Updated land perturbation scheme #148
Changes from 28 commits
6587ec6
c136caf
b0315a7
8c7305c
069b58e
5a20c62
1cef958
142a549
472f674
a5c653a
6da74eb
5c2669b
ef26a23
5c57a58
4cc6e04
c2b4d20
03332e0
35df12e
458e5be
f52214d
02160db
4856f8f
2fdca27
bba8815
dee00b4
e268506
25a039e
9aa3e37
60d6b32
e56e33a
3177ce1
6e6658e
71ca8e4
4e6e1b0
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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') | ||
|
||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thanks for adding these error checks. |
||
|
||
!--- 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 | ||
|
||
|
+2 −2 | physics/GFS_MP_generic.meta | |
+1 −1 | physics/GFS_debug.F90 | |
+8 −5 | physics/GFS_rrtmg_pre.F90 | |
+22 −14 | physics/GFS_rrtmgp_sw_pre.F90 | |
+27 −20 | physics/GFS_rrtmgp_sw_pre.meta | |
+2 −2 | physics/GFS_stochastics.meta | |
+32 −39 | physics/GFS_surface_generic.F90 | |
+32 −52 | physics/GFS_surface_generic.meta | |
+3 −3 | physics/radiation_surface.f | |
+14 −2 | physics/rrtmg_sw_pre.F90 | |
+3 −3 | physics/sfc_drv.f | |
+1 −1 | physics/sfc_drv.meta |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,249 @@ | ||
module GFS_land_perts | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Shouldn't this entire code be in the CCPP "stochastic_physics" group, either as an additional scheme or as part of the existing Once we remove IPD, there will be no access to the soilveg data etc. anymore. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Update: see my comment below for file |
||
|
||
use machine, only: kind_phys | ||
|
||
implicit none | ||
|
||
private | ||
|
||
public :: GFS_apply_lndp | ||
|
||
contains | ||
|
||
!==================================================================== | ||
! GFS_apply_lndp | ||
!==================================================================== | ||
! Driver for applying perturbations to sprecified land states or parameters | ||
! Draper, July 2020. | ||
! Note on location: requires access to namelist_soilveg | ||
|
||
subroutine GFS_apply_lndp(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, & | ||
lndp_prt_list, sfc_wts, xlon, xlat, stype, param_update_flag, & | ||
smc, slc, stc, vfrac, ierr) | ||
|
||
use namelist_soilveg ! needed for maxsmc | ||
|
||
implicit none | ||
|
||
! intent(in) | ||
integer, intent(in) :: blksz(:) | ||
integer, intent(in) :: n_var_lndp, lsoil, lsm | ||
character(len=3), intent(in) :: lndp_var_list(:) | ||
real(kind=kind_phys), intent(in) :: lndp_prt_list(:) | ||
real(kind=kind_phys), intent(in) :: dtf | ||
real(kind=kind_phys), intent(in) :: sfc_wts(:,:,:) | ||
real(kind=kind_phys), intent(in) :: xlon(:,:) | ||
real(kind=kind_phys), intent(in) :: xlat(:,:) | ||
logical, intent(in) :: param_update_flag | ||
! true = parameters have been updated, apply perts | ||
real(kind=kind_phys), intent(in) :: stype(:,:) | ||
|
||
! intent(inout) | ||
real(kind=kind_phys), intent(inout) :: smc(:,:,:) | ||
real(kind=kind_phys), intent(inout) :: slc(:,:,:) | ||
real(kind=kind_phys), intent(inout) :: stc(:,:,:) | ||
real(kind=kind_phys), intent(inout) :: vfrac(:,:) | ||
|
||
! intent(out) | ||
integer, intent(out) :: ierr | ||
|
||
! local | ||
integer :: nblks, print_i, print_nb, i, nb | ||
integer :: this_im, v, soiltyp, k | ||
logical :: print_flag | ||
|
||
real(kind=kind_phys) :: p, min_bound, max_bound, tmp_sic, pert | ||
|
||
! decrease in applied pert with depth | ||
real(kind=kind_phys), dimension(4), parameter :: smc_vertscale = (/1.0,0.5,0.25,0.125/) | ||
real(kind=kind_phys), dimension(4), parameter :: stc_vertscale = (/1.0,0.5,0.25,0.125/) | ||
|
||
! model-dependent values, hard-wired in noah code. | ||
real(kind=kind_phys), dimension(4), parameter :: zs_noah = (/0.1, 0.3, 0.6, 1.0/) | ||
real(kind=kind_phys), parameter :: minsmc = 0.02 | ||
|
||
ierr = 0 | ||
|
||
if (lsm .NE. 1 ) then | ||
write(6,*) 'ERROR: GFS_land_pert assumes LSM is noah, ', & | ||
' may need to adapt variable names for a different LSM' | ||
ierr=10 | ||
return | ||
endif | ||
|
||
nblks = size(blksz) | ||
|
||
call set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb) | ||
|
||
do nb =1,nblks | ||
do i = 1, blksz(nb) | ||
|
||
!if ( nint(Sfcprop(nb)%slmsk(i)) .NE. 1) cycle ! skip if not land | ||
|
||
!if ( ((isot == 1) .and. (soiltyp == 16)) & | ||
! .or.( (isot == 0) .and. (soiltyp == 9 )) ) cycle ! skip if land-ice | ||
|
||
if ( smc(nb,i,1) .EQ. 1.) cycle ! skip non-soil (land-ice and non-land) | ||
! set printing | ||
if ( (i==print_i) .and. (nb==print_nb) ) then | ||
print_flag = .true. | ||
else | ||
print_flag=.false. | ||
endif | ||
|
||
do v = 1,n_var_lndp | ||
select case (trim(lndp_var_list(v))) | ||
!================================================================= | ||
! State updates - performed every cycle | ||
!================================================================= | ||
case('smc') | ||
p=5. | ||
soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc | ||
min_bound = minsmc | ||
max_bound = maxsmc(soiltyp) | ||
|
||
do k=1,lsoil | ||
!store frozen soil moisture | ||
tmp_sic= smc(nb,i,k) - slc(nb,i,k) | ||
|
||
! perturb total soil moisture | ||
! factor of sldepth*1000 converts from mm to m3/m3 | ||
pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zs_noah(k)*1000.) | ||
pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep | ||
! (necessary for state vars only) | ||
call apply_pert('smc',pert,print_flag, smc(nb,i,k),ierr,p,min_bound, max_bound) | ||
|
||
! assign all of applied pert to the liquid soil moisture | ||
slc(nb,i,k) = smc(nb,i,k) - tmp_sic | ||
enddo | ||
|
||
case('stc') | ||
|
||
do k=1,lsoil | ||
pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) | ||
pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep | ||
! (necessary for state vars only) | ||
call apply_pert('stc',pert,print_flag, stc(nb,i,k),ierr) | ||
enddo | ||
!================================================================= | ||
! Parameter updates - only if param_update_flag = TRUE | ||
!================================================================= | ||
case('vgf') ! vegetation fraction | ||
if (param_update_flag) then | ||
p =5. | ||
min_bound=0. | ||
max_bound=1. | ||
|
||
pert = sfc_wts(nb,i,v)*lndp_prt_list(v) | ||
call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) | ||
endif | ||
case default | ||
print*, & | ||
'ERROR: unrecognised lndp_prt_list option in GFS_apply_lndp, exiting', trim(lndp_var_list(v)) | ||
ierr = 10 | ||
return | ||
end select | ||
enddo | ||
enddo | ||
enddo | ||
end subroutine GFS_apply_lndp | ||
|
||
!==================================================================== | ||
! apply_perts | ||
!==================================================================== | ||
! Apply perturbations to selected land states or parameters | ||
|
||
subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax) | ||
|
||
! intent in | ||
logical, intent(in) :: print_flag | ||
real(kind=kind_phys), intent(in) :: pert | ||
character(len=*), intent(in) :: vname ! name of variable being perturbed | ||
|
||
real(kind=kind_phys), optional, intent(in) :: p ! flat-top paramater, 0 = no flat-top | ||
! flat-top function is used for bounded variables | ||
! to reduce the magnitude of perturbations near boundaries. | ||
real(kind=kind_phys), optional, intent(in) :: vmin, vmax ! min,max bounds of variable being perturbed | ||
|
||
! intent (inout) | ||
real(kind=kind_phys), intent(inout) :: state | ||
|
||
! intent out | ||
integer :: ierr | ||
|
||
!local | ||
real(kind=kind_phys) :: z | ||
|
||
if ( print_flag ) then | ||
write(*,*) 'LNDP - applying lndp to ',vname, ', initial value', state | ||
endif | ||
|
||
! apply perturbation | ||
if (present(p) ) then | ||
if ( .not. (present(vmin) .and. present(vmax) )) then | ||
ierr=20 | ||
print*, 'error, flat-top function requires min & max to be specified' | ||
endif | ||
|
||
z = -1. + 2*(state - vmin)/(vmax - vmin) ! flat-top function | ||
state = state + pert*(1-abs(z**p)) | ||
else | ||
state = state + pert | ||
endif | ||
|
||
if (present(vmax)) state = min( state , vmax ) | ||
if (present(vmin)) state = max( state , vmin ) | ||
!state = max( min( state , vmax ), vmin ) | ||
|
||
if ( print_flag ) then | ||
write(*,*) 'LNDP - applying lndp to ',vname, ', final value', state | ||
endif | ||
|
||
end subroutine apply_pert | ||
|
||
|
||
!==================================================================== | ||
! set_printing_nb_i | ||
!==================================================================== | ||
! routine to turn on print flag for selected location | ||
! | ||
subroutine set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb) | ||
|
||
implicit none | ||
|
||
! intent (in) | ||
integer, intent(in) :: blksz(:) | ||
real(kind=kind_phys), intent(in) :: xlon(:,:) | ||
real(kind=kind_phys), intent(in) :: xlat(:,:) | ||
|
||
|
||
! intent (out) | ||
integer, intent(out) :: print_i, print_nb | ||
|
||
! local | ||
integer :: nblks,nb,i | ||
real, parameter :: plon_trunc = 114.9 | ||
real, parameter :: plat_trunc = -26.6 | ||
real, parameter :: delta = 1. | ||
|
||
nblks = size(blksz) | ||
|
||
print_i = -9 | ||
print_nb = -9 | ||
do nb = 1,nblks | ||
do i = 1,blksz(nb) | ||
if ( (xlon(nb,i)*57.29578 > plon_trunc) .and. (xlon(nb,i)*57.29578 < plon_trunc+delta ) .and. & | ||
(xlat(nb,i)*57.29578 > plat_trunc ) .and. (xlat(nb,i)*57.29578 < plat_trunc+delta ) ) then | ||
print_i=i | ||
print_nb=nb | ||
write(*,*) 'LNDP -print flag is on', xlon(nb,i)*57.29578, xlat(nb,i)*57.29578, nb, i | ||
return | ||
endif | ||
enddo | ||
enddo | ||
|
||
end subroutine set_printing_nb_i | ||
|
||
end module GFS_land_perts | ||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. At some point we agreed to no longer accept changes to IPD unless they are backports of GFS v16 implementation bug fixes. Just mentioning it here, but I am not going to ask you to change anything. |
||
!! 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, & | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for reverting
.gitmodules
to where it was originally (the new version in this PR is how it should be, no need to add.git
and also good to remove the unnecessary double whitespace and the trailing whitespace below.