module lndp_apply_perts_mod
use kinddef, only : kind_dbl_prec, kind_phys
use stochy_namelist_def
implicit none
private
public :: lndp_apply_perts
contains
!====================================================================
! lndp_apply_perts
!====================================================================
! Driver for applying perturbations to specified land states or parameters,
! following the LNDP_TYPE=2 scheme.
! Draper, July 2020.
! Can select perturbations by specifying lndp_var_list and lndp_pert_list
! Notes on how the perturbations are added.
! 1. Model prognostic variables.
! If running a long forecast or cycling DA system (as in global UFS),
! perturbing the prognostic variables only at the initial conditions will
! have very limited impact, and they should instead be perturbed at every time step.
! In this case, the pertrubations should be specified as a rate (pert/hr)
! to avoid the ensemble spread being dependent on the model time step.
!
! For a short forecast (~days, as in regional HRRR), can see impact from
! perturbing only the initial conditions. In this case, the perturbation
! is specified as an absolute value (not a rate).
!
! 2. Model parameters:
! The timing of how to perturb the parameters depends on how / whether
! the parameters are updated over time. For the UFS global system, global_cycle
! is periodically called to update the parameters (controlled by FHCYC).
! Each time it's called global_cycle overwrites most of the
! prior parameters (overwriting any perturbations applied to those
! parameters). Hence, the perturbations are applied only immediately after global_cycle
! has been called, and the parameters are not applied as a rate (since they
! don't accumulate).
! For the regional models, FHCYC is 0, and the global_cycle is not called, so
! can perturb parameters every time step. Hence, need to specify the perturbations
! as a rate.
!
! The above cases are controlled by the lndp_model_type variable.
!
! If adding perturbations to new parameters, need to check how/whether
! the parameters are updated by the model.
subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg, &
lsoil, dtf, kdt, n_var_lndp, lndp_var_list, lndp_prt_list, &
sfc_wts, xlon, xlat, stype, smcmax, smcmin, param_update_flag, &
smc, slc, stc, vfrac, alnsf, alnwf, snoalb, semis, zorll, ierr)
implicit none
! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, kdt, iopt_dveg
integer, intent(in) :: lsm, lsm_noah, lsm_ruc, lsm_noahmp
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 just been updated by global_cycle
integer, intent(in) :: stype(:,:)
real(kind=kind_phys), intent(in) :: smcmax(:)
real(kind=kind_phys), intent(in) :: smcmin(:)
! 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(:,:)
real(kind=kind_phys), intent(inout) :: snoalb(:,:)
real(kind=kind_phys), intent(inout) :: alnsf(:,:)
real(kind=kind_phys), intent(inout) :: alnwf(:,:)
real(kind=kind_phys), intent(inout) :: semis(:,:)
real(kind=kind_phys), intent(inout) :: zorll(:,:)
! intent(out)
integer, intent(out) :: ierr
! local
integer :: nblks, print_i, print_nb, i, nb
!integer :: this_im, v, soiltyp, k
integer :: this_im, v, k
logical :: print_flag, do_pert_state, do_pert_param
real(kind=kind_phys) :: p, min_bound, max_bound, pert
real(kind=kind_phys) :: tmp_smc
real(kind=kind_phys) :: conv_hr2tstep, tfactor_state, tfactor_param
real(kind=kind_phys), dimension(lsoil) :: zslayer, smc_vertscale, stc_vertscale
! decrease in applied pert with depth
!-- Noah lsm
real(kind=kind_phys), dimension(4), parameter :: smc_vertscale_noah = (/1.0,0.5,0.25,0.125/)
real(kind=kind_phys), dimension(4), parameter :: stc_vertscale_noah = (/1.0,0.5,0.25,0.125/)
real(kind=kind_phys), dimension(4), parameter :: zs_noah = (/0.1, 0.3, 0.6, 1.0/)
!-- RUC lsm
real(kind=kind_phys), dimension(9), parameter :: smc_vertscale_ruc = (/1.0,0.9,0.8,0.6,0.4,0.2,0.1,0.05,0./)
real(kind=kind_phys), dimension(9), parameter :: stc_vertscale_ruc = (/1.0,0.9,0.8,0.6,0.4,0.2,0.1,0.05,0./)
real(kind=kind_phys), dimension(9), parameter :: zs_ruc = (/0.05, 0.15, 0.20, 0.20, 0.40, 0.60, 0.60, 0.80, 1.00/)
ierr = 0
if (lsm/=lsm_noah .and. lsm/=lsm_ruc .and. lsm/=lsm_noahmp) then
write(6,*) 'ERROR: lndp_apply_pert assumes LSM is Noah, Noah-MP, or RUC,', &
' may need to adapt variable names for a different LSM'
ierr=10
return
endif
if (lsm==lsm_noahmp) then
do v = 1,n_var_lndp
select case (trim(lndp_var_list(v)))
case ('alb','sal','emi','zol')
print*, &
'ERROR: lndp_prt_list option in lndp_apply_pert', trim(lndp_var_list(v)) , &
' has not been checked for Noah-MP. Please check how the parameter is set/updated ', &
' before applying. Note: in Noah-MP many variables that have traditionally been', &
' externally specified parameters are now prognostic. Also, many parameters are', &
' set at each Noah-MP model call from data tables'
ierr = 10
return
case ('veg')
if (iopt_dveg .NE. 4 ) then
print*, &
'ERROR: veg perturbations have not yet been coded for dveg options other than 4'
ierr = 10
return
endif
end select
enddo
endif
! for perturbations applied as a rate, lndp_prt_list input is per hour. Converts to per timestep
conv_hr2tstep = dtf/3600. ! conversion factor from per hour to per tstep.
! determine whether updating state variables and/or parameters
do_pert_state=.false.
do_pert_param=.false.
select case (lndp_model_type)
case(1) ! global, perturb states every time step (pert applied as a rate)
! perturb parameters only when they've been update (pert is not a rate)
do_pert_state=.true.
tfactor_state=conv_hr2tstep
if (param_update_flag) then
do_pert_param=.true.
tfactor_param=1.
endif
case(2) ! regional, perturb states only at first time step (pert is not a rate)
! perurb parameters at every time step (pert is a rate)
if ( kdt == 2 ) then
do_pert_state=.true.
tfactor_state=1.
endif
do_pert_param = .true.
tfactor_param = conv_hr2tstep
case(3) ! special case to apply perturbations at initial time step only (pert is not a rate)
if ( kdt == 2 ) then
do_pert_state=.true.
tfactor_state=1.
do_pert_param=.true.
tfactor_param=1.
endif
case default
print*, &
'ERROR: unrecognised lndp_model_type option in lndp_apply_pert, exiting', trim(lndp_var_list(v))
ierr = 10
return
end select
if (lsm == lsm_noah .or. lsm == lsm_noahmp) then
do k = 1, lsoil
zslayer(k) = zs_noah(k)
smc_vertscale(k) = smc_vertscale_noah(k)
stc_vertscale(k) = stc_vertscale_noah(k)
enddo
elseif (lsm == lsm_ruc) then
do k = 1, lsoil
zslayer(k) = zs_ruc(k)
smc_vertscale(k) = smc_vertscale_ruc(k)
stc_vertscale(k) = stc_vertscale_ruc(k)
enddo
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 ( 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')
if (do_pert_state) then
p=5.
min_bound = smcmin(stype(nb,i))
max_bound = smcmax(stype(nb,i))
! with RUC LSM perturb smc only at time step = 2, as in HRRR
do k=1,lsoil
! apply perts to a copy of smc, retain original smc
! for later update to liquid soil moisture.
! note: previously we were saving the ice water content
! (smc-slc) and subtracting this from the perturbed smc to
! get the perturbed slc. This was introducing small errors in the slc
! when passing back to the calling program, I think due to precision issues,
! as the ice content is typically zero. Clara Draper, March, 2022.
tmp_smc = smc(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)/(zslayer(k)*1000.)
pert = pert*tfactor_state
call apply_pert('smc',pert,print_flag, tmp_smc,ierr,p,min_bound, max_bound)
! assign all of applied pert to the liquid soil moisture
slc(nb,i,k) = slc(nb,i,k) + tmp_smc - smc(nb,i,k)
smc(nb,i,k) = tmp_smc
enddo
endif
case('stc')
if (do_pert_state) then
do k=1,lsoil
pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v)
pert = tfactor_state
call apply_pert('stc',pert,print_flag, stc(nb,i,k),ierr)
enddo
endif
!=================================================================
! Parameter updates
!=================================================================
! are all of the params below included in noah?
case('vgf') ! vegetation fraction
if (do_pert_param) then
p =5.
min_bound=0.
max_bound=1.
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*tfactor_param
call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound)
endif
case('alb') ! albedo
if (do_pert_param) then
p =5.
min_bound=0.0
max_bound=0.4
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*tfactor_param
call apply_pert ('alnsf',pert,print_flag, alnsf(nb,i), ierr,p,min_bound, max_bound)
call apply_pert ('alnwf',pert,print_flag, alnwf(nb,i), ierr,p,min_bound, max_bound)
endif
case('sal') ! snow albedo
if (do_pert_param) then
p =5.
min_bound=0.3
max_bound=0.85
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*tfactor_param
call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound)
endif
case('emi') ! emissivity
if (do_pert_param) then
p =5.
min_bound=0.8
max_bound=1.
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*tfactor_param
call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound)
endif
case('zol') ! land roughness length
if (do_pert_param) then
p =5.
min_bound=0.
max_bound=300.
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*tfactor_param
call apply_pert ('zol',pert,print_flag, zorll(nb,i), ierr,p,min_bound, max_bound)
endif
case default
print*, &
'ERROR: unrecognised lndp_prt_list option in lndp_apply_pert, exiting', trim(lndp_var_list(v))
ierr = 10
return
end select
enddo
enddo
enddo
end subroutine lndp_apply_perts
!====================================================================
! 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 )
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 lndp_apply_perts_mod