Skip to content

More cleanup for Tanya #2

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.

Sign up for GitHub

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

Merged
tanyasmirnova merged 1 commit into from
Dec 30, 2020
Merged

More cleanup for Tanya #2

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
150 changes: 68 additions & 82 deletions lndp_apply_perts.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,35 +14,33 @@ module lndp_apply_perts_mod
! lndp_apply_perts
!====================================================================
! Driver for applying perturbations to sprecified land states or parameters
! Draper, July 2020.
! Draper, July 2020.

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, smcmax, smcmin, &
param_update_flag, &
smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, &
snoalb, semis, ierr)
subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, &
dtf, 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, ierr)

implicit none

! intent(in)
! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, lsm
integer, intent(in) :: lsoil_lsm, lsm_ruc
integer, intent(in) :: n_var_lndp, lsoil
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(:,:,:)
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)
logical, intent(in) :: param_update_flag
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) :: smcmax(:)
real(kind=kind_dbl_prec), intent(in) :: smcmin(:)
real(kind=kind_dbl_prec), intent(in) :: zs_lsm(:)

! 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(:,:,:)
Expand All @@ -56,16 +54,16 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm,
real(kind=kind_dbl_prec), intent(inout) :: facwf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: semis(:,:)

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

! local
integer :: nblks, print_i, print_nb, i, nb
integer :: this_im, v, soiltyp, k, nsoil
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), dimension(max(lsoil,lsoil_lsm)) :: zslayer, smc_vertscale, stc_vertscale
real(kind=kind_dbl_prec), dimension(lsoil) :: zslayer, smc_vertscale, stc_vertscale

! decrease in applied pert with depth
!-- Noah lsm
Expand All @@ -75,17 +73,15 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm,
!-- 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.
real(kind=kind_dbl_prec), dimension(9), parameter :: zs_ruc = (/0.00, 0.05, 0.20, 0.40, 0.60, 1.00, 1.60, 2.20, 3.00/)

ierr = 0
ierr = 0

if (lsm .NE. 1 .and. lsm .ne. 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
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

!write (0,*) 'Input to lndp_apply_pert'
Expand All @@ -94,21 +90,15 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm,
!write (0,*) 'n_var_lndp, lndp_var_list =', n_var_lndp, lndp_var_list
!write (0,*) 'smcmin =',smcmin

zslayer(:) = 0.
smc_vertscale(:) = 0.
stc_vertscale(:) = 0.
if (lsm == 1) then
nsoil = lsoil
do k = 1, nsoil
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
nsoil = lsoil_lsm
zslayer(nsoil) = zsbot - zs_lsm(nsoil)
do k = 1, nsoil-1
zslayer(k) = zs_lsm(k+1) - zs_lsm(k)
do k = 1, lsoil
zslayer(k) = zs_ruc(k)
smc_vertscale(k) = smc_vertscale_ruc(k)
stc_vertscale(k) = stc_vertscale_ruc(k)
enddo
Expand All @@ -128,42 +118,42 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm,

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 = smcmin(soiltyp)
max_bound = smcmax(soiltyp)
do k=1,nsoil

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
! 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
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

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

do k=1,nsoil
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)
Expand All @@ -173,58 +163,54 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm,
! Parameter updates - only if param_update_flag = TRUE
!=================================================================
case('vgf') ! vegetation fraction
!if (param_update_flag) then
p =5.
if (param_update_flag) then
p =5.
min_bound=0.
max_bound=1.

! lndp_prt_list(v) = 0.1 in input.nml
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
case('alb') ! albedo
!if (param_update_flag) then
if (param_update_flag) then
p =5.
min_bound=0.0
max_bound=0.4

! lndp_prt_list(v) = 0.4 (wrf)
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
endif
case('sal') ! snow albedo
!if (param_update_flag) then
if (param_update_flag) then
p =5.
min_bound=0.3
max_bound=0.85

! lndp_prt_list(v) = 0.4 (wrf)
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
endif
case('emi') ! emissivity
!if (param_update_flag) then
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
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 @@ -241,14 +227,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 @@ -259,9 +245,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 @@ -280,21 +266,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 @@ -317,7 +303,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