Skip to content
Merged
2 changes: 1 addition & 1 deletion compns_stochy.F90
Original file line number Diff line number Diff line change
Expand Up @@ -294,7 +294,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
'land perturbations will be applied to selected paramaters, using newer scheme designed for DA ens spread'
do k =1,n_var_lndp
select case (lndp_var_list(k))
case('vgf','smc','stc')
case('vgf','smc','stc','alb', 'sal')
if (me==0) print*, 'land perturbation will be applied to ', lndp_var_list(k)
case default
print*, 'ERROR: land perturbation requested for new parameter - will need to be coded in lndp_apply_pert', lndp_var_list(k)
Expand Down
124 changes: 105 additions & 19 deletions lndp_apply_perts.F90
Original file line number Diff line number Diff line change
Expand Up @@ -15,17 +15,22 @@ module lndp_apply_perts_mod
!====================================================================
! Driver for applying perturbations to sprecified land states or parameters
! Draper, July 2020.
! Note on location: requires access to namelist_soilveg
! Note on location: requires access to namelist_soilveg or namelist_soilveg_ruc
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated

subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
lndp_prt_list, sfc_wts, xlon, xlat, stype, maxsmc,param_update_flag, &
smc, slc, stc, vfrac, ierr)
subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, &
dtf, n_var_lndp, lndp_var_list, &
lndp_prt_list, sfc_wts, xlon, xlat, stype, maxsmcnoah, maxsmc, &
param_update_flag, &
smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, &
snoalb, ierr)
!snoalb, semis, ierr)

implicit none

! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, lsm
integer, intent(in) :: lsoil_lsm, lsm_ruc
character(len=3), intent(in) :: lndp_var_list(:)
real(kind=kind_dbl_prec), intent(in) :: lndp_prt_list(:)
real(kind=kind_dbl_prec), intent(in) :: dtf
Expand All @@ -35,41 +40,82 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
logical, intent(in) :: param_update_flag
! true = parameters have been updated, apply perts
real(kind=kind_dbl_prec), intent(in) :: stype(:,:)
real(kind=kind_dbl_prec), intent(in) :: maxsmcnoah(:)
real(kind=kind_dbl_prec), intent(in) :: maxsmc(:)
real(kind=kind_dbl_prec), intent(in) :: zs_lsm(:)

! intent(inout)
real(kind=kind_dbl_prec), intent(inout) :: smc(:,:,:)
real(kind=kind_dbl_prec), intent(inout) :: slc(:,:,:)
real(kind=kind_dbl_prec), intent(inout) :: stc(:,:,:)
real(kind=kind_dbl_prec), intent(inout) :: vfrac(:,:)
real(kind=kind_dbl_prec), intent(inout) :: snoalb(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alvsf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alnsf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alvwf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alnwf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: facsf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: facwf(:,:)
!real(kind=kind_dbl_prec), intent(inout) :: semis(:,:)

! 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, soiltyp, k, nsoil
logical :: print_flag

real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert
real(kind=kind_dbl_prec), dimension(max(lsoil,lsoil_lsm)) :: zslayer, smc_vertscale, stc_vertscale

! decrease in applied pert with depth
real(kind=kind_dbl_prec), dimension(4), parameter :: smc_vertscale = (/1.0,0.5,0.25,0.125/)
real(kind=kind_dbl_prec), dimension(4), parameter :: stc_vertscale = (/1.0,0.5,0.25,0.125/)

! model-dependent values, hard-wired in noah code.
!-- Noah lsm
real(kind=kind_dbl_prec), dimension(4), parameter :: smc_vertscale_noah = (/1.0,0.5,0.25,0.125/)
real(kind=kind_dbl_prec), dimension(4), parameter :: stc_vertscale_noah = (/1.0,0.5,0.25,0.125/)
real(kind=kind_dbl_prec), dimension(4), parameter :: zs_noah = (/0.1, 0.3, 0.6, 1.0/)
real(kind=kind_dbl_prec), parameter :: minsmc = 0.02
!-- RUC lsm
real(kind=kind_dbl_prec), 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_dbl_prec), 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_dbl_prec), parameter :: minsmc = 0.02
real(kind=kind_dbl_prec), parameter :: zsbot = 6.

ierr = 0

if (lsm .NE. 1 ) then
write(6,*) 'ERROR: lndp_apply_pert assumes LSM is noah, ', &
if (lsm .NE. 1 .and. lsm .ne. lsm_ruc) then
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated
write(6,*) 'ERROR: lndp_apply_pert assumes LSM is noah or ruc, ', &
' may need to adapt variable names for a different LSM'
ierr=10
return
endif

!write (0,*) 'Input to lndp_apply_pert'
!write (0,*) 'lsm, lsoil, lsm_ruc, lsoil_lsm =', lsm, lsoil, lsm_ruc, lsoil_lsm
!write (0,*) 'zs_lsm =', zs_lsm
!write (0,*) 'n_var_lndp, lndp_var_list =', n_var_lndp, lndp_var_list
!write (0,*) 'maxsmcnoah, maxsmc =', maxsmcnoah, maxsmc

zslayer(:) = 0.
smc_vertscale(:) = 0.
stc_vertscale(:) = 0.
if (lsm == 1) then
nsoil = lsoil
do k = 1, nsoil
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
nsoil = lsoil_lsm
zslayer(nsoil) = zsbot - zs_lsm(nsoil)
do k = 1, nsoil-1
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated
zslayer(k) = zs_lsm(k+1) - zs_lsm(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)
Expand Down Expand Up @@ -99,26 +145,29 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
p=5.
soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc
min_bound = minsmc
max_bound = maxsmc(soiltyp)

do k=1,lsoil
max_bound = maxsmcnoah(soiltyp)
!write (*,*) 'max_bound=', max_bound

do k=1,nsoil
!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.)
! lndp_prt_list(v) = 0.3 in input.nml
pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zslayer(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
!write (0,*) 'nb, i, smc(nb,i,:)', nb, i, smc(nb,i,:)

case('stc')

do k=1,lsoil
do k=1,nsoil
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)
Expand All @@ -128,14 +177,51 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
! Parameter updates - only if param_update_flag = TRUE
Comment thread
DomHeinzeller marked this conversation as resolved.
!=================================================================
case('vgf') ! vegetation fraction
if (param_update_flag) then
!if (param_update_flag) then
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated
p =5.
min_bound=0.
max_bound=1.

! lndp_prt_list(v) = 0.1 in input.nml
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated
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
!endif
!write (0,*) 'nb, i, vfrac(nb,i)', nb, i, vfrac(nb,i)
case('alb') ! albedo
!if (param_update_flag) then
p =5.
min_bound=0.0
max_bound=0.4

! lndp_prt_list(v) = 0.4 (wrf)
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
!call apply_pert ('alvsf',pert,print_flag, alvsf(nb,i), ierr,p,min_bound, max_bound)
call apply_pert ('alnsf',pert,print_flag, alnsf(nb,i), ierr,p,min_bound, max_bound)
!call apply_pert ('alvwf',pert,print_flag, alvwf(nb,i), ierr,p,min_bound, max_bound)
call apply_pert ('alnwf',pert,print_flag, alnwf(nb,i), ierr,p,min_bound, max_bound)
!call apply_pert ('facsf',pert,print_flag, facsf(nb,i), ierr,p,min_bound, max_bound)
!call apply_pert ('facwf',pert,print_flag, facwf(nb,i), ierr,p,min_bound, max_bound)
!endif
case('sal') ! snow albedo
!if (param_update_flag) then
p =5.
min_bound=0.3
max_bound=0.85

! lndp_prt_list(v) = 0.4 (wrf)
Comment thread
DomHeinzeller marked this conversation as resolved.
Outdated
pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound)
!endif
!case('emi') ! emissivity
! if (param_update_flag) then
! p =5.
! min_bound=0.
! max_bound=1.

! lndp_prt_list(v) = 0.1 (wrf)
! pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
! call apply_pert ('semis',pert,print_flag, semis(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))
Expand Down