Skip to content

Commit

Permalink
Merge branch 'master' into ocn_stoch
Browse files Browse the repository at this point in the history
  • Loading branch information
@pjpegion
pjpegion authored Jan 22, 2021
2 parents 8d42b09 + c39bb8a commit d52b5cb
Show file tree
Hide file tree
Showing 2 changed files with 161 additions and 64 deletions.
2 changes: 1 addition & 1 deletion compns_stochy.F90
View file
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','emi','zol')
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
223 changes: 160 additions & 63 deletions lndp_apply_perts.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,58 +14,105 @@ module lndp_apply_perts_mod
! lndp_apply_perts
!====================================================================
! Driver for applying perturbations to sprecified land states or parameters
! Draper, July 2020.
! Note on location: requires access to namelist_soilveg
! Draper, July 2020.

subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, &
dtf, kdt, lndp_each_step, &
n_var_lndp, lndp_var_list, lndp_prt_list, &
sfc_wts, xlon, xlat, stype, smcmax, smcmin, param_update_flag, &
smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, &
snoalb, semis, zorll, ierr)

subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
lndp_prt_list, sfc_wts, stype, maxsmc,param_update_flag, &
smc, slc, stc, vfrac, ierr)
use stochy_data_mod, only : gis_stochy
implicit none

! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, lsm
character(len=3), intent(in) :: lndp_var_list(:)
! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, kdt
logical, intent(in) :: lndp_each_step
integer, intent(in) :: lsm, lsm_noah, 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
real(kind=kind_dbl_prec), intent(in) :: sfc_wts(:,:,:)
logical, intent(in) :: param_update_flag
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)
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) :: maxsmc(:)
real(kind=kind_dbl_prec), intent(in) :: smcmax(:)
real(kind=kind_dbl_prec), intent(in) :: smcmin(:)

! intent(inout)
! 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(:,:)

! intent(out)
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(:,:)
real(kind=kind_dbl_prec), intent(inout) :: zorll(:,:)

! intent(out)
integer, intent(out) :: ierr

! local
integer :: nblks, print_i, print_nb, i, nb
integer :: this_im, v, soiltyp, k
logical :: print_flag
integer :: this_im, v, soiltyp, k
logical :: print_flag

real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert
real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert, factor
real(kind=kind_dbl_prec), dimension(lsoil) :: 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), 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) then
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

ierr = 0
!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,*) 'smcmin =',smcmin

! lndp_prt_list input is per hour, factor converts to per timestep
! Do conversion only when variables are perturbed at every time step
if(lndp_each_step) then
factor = dtf/3600.
else
factor = 1.
endif

if (lsm .NE. 1 ) then
write(6,*) 'ERROR: lndp_apply_pert assumes LSM is noah, ', &
' may need to adapt variable names for a different LSM'
ierr=10
return
if (lsm == lsm_noah) 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)
Expand All @@ -82,39 +129,43 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &

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
if ( (i==print_i) .and. (nb==print_nb) ) then
print_flag = .true.
else
print_flag=.false.
else
print_flag=.false.
endif

do v = 1,n_var_lndp
do v = 1,n_var_lndp
select case (trim(lndp_var_list(v)))
!=================================================================
! State updates - performed every cycle
!=================================================================
case('smc')
p=5.
case('smc')
p=5.
soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc
min_bound = minsmc
max_bound = maxsmc(soiltyp)
min_bound = smcmin(soiltyp)
max_bound = smcmax(soiltyp)

if ((lsm /= lsm_ruc) .or. (lsm == lsm_ruc .and. kdt == 2)) then
! with RUC LSM perturb smc only at time step = 2, as in HRRR
do k=1,lsoil
!store frozen soil moisture
tmp_sic= smc(nb,i,k) - slc(nb,i,k)

! perturb total soil moisture
! 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
! 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
! assign all of applied pert to the liquid soil moisture
slc(nb,i,k) = smc(nb,i,k) - tmp_sic
enddo
endif

case('stc')
case('stc')

do k=1,lsoil
pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v)
Expand All @@ -126,22 +177,68 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
! Parameter updates - only if param_update_flag = TRUE
!=================================================================
case('vgf') ! vegetation fraction
if (param_update_flag) then
p =5.
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.
max_bound=1.

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound)
endif
case default
case('alb') ! albedo
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.0
max_bound=0.4

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
!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 .or. lndp_each_step) then
p =5.
min_bound=0.3
max_bound=0.85

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound)
endif
case('emi') ! emissivity
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.8
max_bound=1.

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound)
endif
case('zol') ! land roughness length
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.
max_bound=300.

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
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
'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

Expand All @@ -158,14 +255,14 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
character(len=*), intent(in) :: vname ! name of variable being perturbed

real(kind=kind_dbl_prec), optional, intent(in) :: p ! flat-top paramater, 0 = no flat-top
! flat-top function is used for bounded variables
! flat-top function is used for bounded variables
! to reduce the magnitude of perturbations near boundaries.
real(kind=kind_dbl_prec), optional, intent(in) :: vmin, vmax ! min,max bounds of variable being perturbed

! intent (inout)
real(kind=kind_dbl_prec), intent(inout) :: state
! intent out

! intent out
integer :: ierr

!local
Expand All @@ -176,9 +273,9 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
endif

! apply perturbation
if (present(p) ) then
if (present(p) ) then
if ( .not. (present(vmin) .and. present(vmax) )) then
ierr=20
ierr=20
print*, 'error, flat-top function requires min & max to be specified'
endif

Expand All @@ -197,21 +294,21 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
endif

end subroutine apply_pert


!====================================================================
! set_printing_nb_i
! 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
implicit none

! intent (in)
integer, intent(in) :: blksz(:)
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)


! intent (out)
Expand All @@ -234,7 +331,7 @@ subroutine set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb)
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
return
endif
enddo
enddo
Expand Down

0 comments on commit d52b5cb

Please sign in to comment.