Skip to content

Commit

Permalink
lndp updates for new stochastic_physics_wrapper
Browse files Browse the repository at this point in the history
  • Loading branch information
@ClaraDraper-NOAA
ClaraDraper-NOAA committed Aug 14, 2020
1 parent d7332a5 commit ae4a3bf
Show file tree
Hide file tree
Showing 8 changed files with 347 additions and 77 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,10 @@ add_library(
./cellular_automata_global.F90
./cellular_automata_sgs.F90
./update_ca.F90
./lndp_apply_perts.F90
)

target_link_libraries(stochastic_physics sp::sp_d)
target_link_libraries(stochastic_physics fms)
target_link_libraries(stochastic_physics gfsphysics)

9 changes: 8 additions & 1 deletion compns_stochy.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,6 +267,13 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
lndp_prt_list( n_var_lndp) = lndp_prt_list(k)
endif
enddo
if (n_var_lndp > max_n_var_lndp) then
print*, 'ERROR: land perturbation requested for too many parameters', &
'increase max_n_var_lndp'
iret = 10
return
endif


if (lndp_type==1) then
if (me==0) print*, &
Expand All @@ -290,7 +297,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
case('vgf','smc','stc')
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 GFS_land_pert', lndp_var_list(k)
print*, 'ERROR: land perturbation requested for new parameter - will need to be coded in lndp_apply_pert', lndp_var_list(k)
iret = 10
return
end select
Expand Down
11 changes: 5 additions & 6 deletions get_stochy_pattern.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ end subroutine get_random_pattern_fv3
!>@brief The subroutine 'get_random_pattern_fv3_sfc' converts spherical harmonics to the gaussian grid then interpolates to the cubed-sphere grid once
!>@details This subroutine is for a 2-D (lat-lon) scalar field
subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,&
gis_stochy,xlat,xlon,blksz,nblks,maxlen,do_advance_patterns,pattern_3d)
gis_stochy,xlat,xlon,blksz,nblks,maxlen,pattern_3d)
!\callgraph

! generate a random pattern for stochastic physics
Expand All @@ -112,7 +112,6 @@ subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,&
type(stochy_internal_state), target :: gis_stochy
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:),xlon(:,:)
integer,intent(in) :: npatterns,blksz(:),nblks,maxlen
logical, intent(in) :: do_advance_pattern
real(kind=kind_dbl_prec), intent(out) :: pattern_3d(nblks,maxlen,n_var_lndp)

integer i,j,l,lat,ierr,n,nn,k,nt
Expand All @@ -133,8 +132,8 @@ subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,&
kmsk0 = 0
glolal = 0.
do n=1,npatterns
if (do_advance_pattern) call patterngenerator_advance(rpattern(n),k,.false.)
if (is_master()) print *, 'Random pattern for SFC-PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k))
call patterngenerator_advance(rpattern(n),k,.false.)
if (is_master()) print *, 'Random pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k))
call scalarspect_to_gaugrid( &
rpattern(n)%spec_e(:,:,k),rpattern(n)%spec_o(:,:,k),wrk2d,&
gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,&
Expand All @@ -153,7 +152,7 @@ subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,&
enddo

call mp_reduce_sum(workg,lonf,latg)
if (is_master()) print *, 'workg after mp_reduce_sum for SFC-PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(workg), maxval(workg)
if (is_master()) print *, 'workg after mp_reduce_sum for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(workg), maxval(workg)

! interpolate to cube grid

Expand All @@ -168,7 +167,7 @@ subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,&
pattern_3d(blk,:,k)=pattern_1d(:)
end associate
enddo
if (is_master()) print *, '3D pattern for SFC-PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k))
if (is_master()) print *, '3D pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k))
deallocate(rslmsk)
deallocate(workg)

Expand Down
249 changes: 249 additions & 0 deletions lndp_apply_perts.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,249 @@
module lndp_apply_perts_mod

use kinddef, only : kind_dbl_prec

implicit none

private

public :: lndp_apply_perts

contains

!====================================================================
! lndp_apply_perts
!====================================================================
! Driver for applying perturbations to sprecified land states or parameters
! Draper, July 2020.
! Note on location: requires access to namelist_soilveg

subroutine lndp_apply_perts(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_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
! true = parameters have been updated, apply perts
real(kind=kind_dbl_prec), intent(in) :: stype(:,:)

! 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)
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_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert

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

ierr = 0

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
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 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_dbl_prec), intent(in) :: pert
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
! 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
integer :: ierr

!local
real(kind=kind_dbl_prec) :: 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_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), 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


3 changes: 2 additions & 1 deletion makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,8 @@ SRCS_F90 = \
./initialize_spectral_mod.F90 \
./cellular_automata_global.F90 \
./cellular_automata_sgs.F90 \
./update_ca.F90
./update_ca.F90 \
./lndp_apply_perts.F90

SRCS_c =

Expand Down
Loading

0 comments on commit ae4a3bf

Please sign in to comment.