Skip to content
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

Add support for Stochastically Perturbed Parameterizations (SPP) in FV3 #51

Merged
merged 15 commits into from
Feb 23, 2022
Merged

Add support for Stochastically Perturbed Parameterizations (SPP) in FV3 #51

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
34c9d84
Add SPP option to several physics parameterizations
Nov 11, 2021
7f1159e
bug fix
Nov 22, 2021
d2e1b95
Merge branch 'master' of https://github.com/noaa-psd/stochastic_physi…
Nov 23, 2021
0b868b1
merge fix
Nov 23, 2021
40ec485
merge fixes
Nov 23, 2021
86755ac
add iostat to new nml read
llpcarson Nov 30, 2021
912c52b
Merge branch 'master' of https://github.com/noaa-psd/stochastic_physi…
JeffBeck-NOAA Jan 20, 2022
d4446d9
Fix application of pattern generator when using SPP.
JeffBeck-NOAA Feb 11, 2022
049984c
Merge branch 'master' of https://github.com/noaa-psd/stochastic_physi…
JeffBeck-NOAA Feb 11, 2022
7a74cdd
Fix write call.
JeffBeck-NOAA Feb 11, 2022
095cd28
Fix to stochastic pattern scaling - match WRF-based SPP code
JeffBeck-NOAA Feb 12, 2022
5d0d239
Fix cutoff thresholding
JeffBeck-NOAA Feb 12, 2022
659c621
Final spp_wts configuration.
JeffBeck-NOAA Feb 13, 2022
079820e
Update SPP namelist entry
JeffBeck-NOAA Feb 15, 2022
450508d
Update SPP namelist entry
JeffBeck-NOAA Feb 16, 2022
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
62 changes: 62 additions & 0 deletions compns_stochy.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,9 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt
namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, &
lndp_tau,lndp_lscale
! For SPP physics parameterization perterbations
namelist /nam_sppperts/spp_var_list, spp_prt_list, iseed_spp, &
spp_tau,spp_lscale,spp_sigtop1, spp_sigtop2,spp_stddev_cutoff

rerth =6.3712e+6 ! radius of earth (m)
tol=0.01 ! tolerance for calculations
Expand All @@ -79,12 +82,15 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
skeb = -999. ! stochastic KE backscatter amplitude
lndp_var_list = 'XXX'
lndp_prt_list = -999.
spp_var_list = 'XXX'
spp_prt_list = -999.
! logicals
do_sppt = .false.
use_zmtnblck = .false.
new_lscale = .false.
do_shum = .false.
do_skeb = .false.
do_spp = .false.
! C. Draper July 2020.
! input land pert variables:
! LNDP_TYPE = 0
Expand Down Expand Up @@ -125,6 +131,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
skeb_sigtop1 = 0.1
skeb_sigtop2 = 0.025
shum_sigefold = 0.2
spp_sigtop1 = 0.1
spp_sigtop2 = 0.025
! reduce amplitude of sppt near surface (lowest 2 levels)
sppt_sfclimit = .false.
! gaussian or power law variance spectrum for skeb (0: gaussian, 1:
Expand All @@ -133,6 +141,11 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
skeb_varspect_opt = 0
sppt_logit = .false. ! logit transform for sppt to bounded interval [-1,+1]
stochini = .false. ! true= read in pattern, false=initialize from seed
! For SPP perturbations
spp_lscale = -999. ! length scales
spp_tau = -999. ! time scales
spp_stddev_cutoff = 0 ! cutoff/limit for std-dev (zero==no limit applied)
iseed_spp = 0 ! random seeds (if 0 use system clock)

#ifdef INTERNAL_FILE_NML
read(input_nml_file, nml=nam_stochy)
Expand All @@ -149,8 +162,19 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
read(nlunit,nam_sfcperts)
#endif

#ifdef INTERNAL_FILE_NML
read(input_nml_file, nml=nam_sppperts, iostat=ios)
#else
rewind (nlunit)
open (unit=nlunit, file=fn_nml, action='READ', status='OLD', iostat=ios)
read(nlunit,nam_sppperts)
#endif

if (me == 0) then
print *,' in compns_stochy'
print*,'spp_lscale=',spp_lscale
print*,'spp_tau=',spp_tau
print*,'spp_stddev_cutoff=',spp_stddev_cutoff
endif

! PJP stochastic physics additions
Expand Down Expand Up @@ -211,6 +235,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
if (skeb(k).GT.0) l_min=min(skeb_lscale(k),l_min)
enddo
if (lndp_type.GT.0) l_min=min(lndp_lscale(1),l_min)
if (spp_prt_list(1).GT.0) l_min=min(spp_lscale(1),l_min)
!ntrunc=1.5*circ/l_min
ntrunc=circ/l_min
if (me==0) print*,'ntrunc calculated from l_min',l_min,ntrunc
Expand Down Expand Up @@ -298,6 +323,41 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
iret = 10
return
end select
!
! SPP perts - parse nml input
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! count requested pert variables
n_var_spp= 0
do k =1,size(spp_var_list)
if ( (spp_var_list(k) .EQ. 'XXX') .or. (spp_prt_list(k) .LE. 0.) ) then
cycle
else
n_var_spp=n_var_spp+1
spp_var_list( n_var_spp) = spp_var_list(k) !
spp_prt_list( n_var_spp) = spp_prt_list(k)
endif
enddo
IF (n_var_spp > 0 ) THEN
do_spp=.true.
ENDIF
if (n_var_spp > max_n_var_spp) then
print*, 'ERROR: SPP physics perturbation requested for too many parameters', &
'increase max_n_var_spp'
iret = 10
return
endif
if (me==0) print*, &
'SPP physics perturbations will be applied to selected parameters', n_var_spp
do k =1,n_var_spp
select case (spp_var_list(k))
case('pbl','sfc', 'mp','rad','gwd')
if (me==0) print*, 'SPP physics perturbation will be applied to ', spp_var_list(k)
case default
print*, 'ERROR: SPP physics perturbation requested for new parameter - will need to be coded in spp_apply_pert', spp_var_list(k)
iret = 10
return
end select
enddo
!
! All checks are successful.
!
Expand All @@ -308,6 +368,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
print *, ' do_skeb : ', do_skeb
print *, ' lndp_type : ', lndp_type
if (lndp_type .NE. 0) print *, ' n_var_lndp : ', n_var_lndp
print *, ' do_spp : ', do_spp
print *, ' n_var_spp : ', n_var_spp
endif
iret = 0
!
Expand Down
88 changes: 81 additions & 7 deletions get_stochy_pattern.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,13 @@ module get_stochy_pattern_mod
len_trio_ls, ls_dim, stochy_la2ga, &
coslat_a, latg, levs, lonf, skeblevs,&
four_to_grid, spec_to_four, dezouv_stochy,dozeuv_stochy
use stochy_namelist_def, only : n_var_lndp, ntrunc, stochini
use stochy_namelist_def, only : n_var_lndp, ntrunc, stochini,n_var_spp
use stochy_data_mod, only : gg_lats, gg_lons, inttyp, nskeb, nshum, nsppt, &
nocnsppt,nepbl,nlndp, &
rnlat, rpattern_sfc, rpattern_skeb, &
rpattern_shum, rpattern_sppt, rpattern_ocnsppt,&
rpattern_epbl1, rpattern_epbl2, skebu_save, &
nspp,rpattern_spp, &
skebv_save, skeb_vwts, skeb_vpts, wlon
use stochy_patterngenerator_mod, only: random_pattern, ndimspec, &
patterngenerator_advance
Expand All @@ -20,7 +21,7 @@ module get_stochy_pattern_mod
implicit none
private

public get_random_pattern_vector
public get_random_pattern_vector,get_random_pattern_spp
public get_random_pattern_sfc,get_random_pattern_scalar
public write_stoch_restart_atm,write_stoch_restart_ocn
logical :: first_call=.true.
Expand Down Expand Up @@ -274,7 +275,61 @@ subroutine get_random_pattern_scalar(rpattern,npatterns,&
deallocate(workg)

end subroutine get_random_pattern_scalar
!

!>@brief The subroutine 'get_random_pattern_spp' converts spherical harmonics
!to the gaussian grid then interpolates to the target grid
!>@details This subroutine is for a 2-D (lat-lon) scalar field
subroutine get_random_pattern_spp(rpattern,npatterns,&
gis_stochy,pattern_3d)

! generate a random pattern for stochastic physics
implicit none
type(random_pattern), intent(inout) :: rpattern(npatterns)
type(stochy_internal_state) :: gis_stochy
integer,intent(in):: npatterns

integer i,j,lat,n

! logical lprint

real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: workg
real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a)
integer kmsk0(lonf,gis_stochy%lats_node_a)
real(kind=kind_dbl_prec) :: pattern_3d(gis_stochy%nx,gis_stochy%ny,npatterns)
real(kind=kind_dbl_prec) :: pattern_1d(gis_stochy%nx)

allocate(workg(lonf,latg))
do n=1,npatterns
kmsk0 = 0
glolal = 0.
call patterngenerator_advance(rpattern(n),1,.false.)
call scalarspect_to_gaugrid(rpattern(n),gis_stochy, &
glolal,1)

workg = 0.
do j=1,gis_stochy%lats_node_a
lat=gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+j)
do i=1,lonf
workg(i,lat) = glolal(i,j)
enddo
enddo

call mp_reduce_sum(workg,lonf,latg)

! interpolate to cube grid
do j=1,gis_stochy%ny
pattern_1d = 0
associate( tlats=>gis_stochy%parent_lats(1:gis_stochy%len(j),j),&
tlons=>gis_stochy%parent_lons(1:gis_stochy%len(j),j))
call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
pattern_3d(:,j,n)=pattern_1d(:)
end associate
enddo
enddo
deallocate(workg)

end subroutine get_random_pattern_spp

!>@brief The subroutine 'scalarspect_to_gaugrid' converts scalar spherical harmonics to a scalar on a gaussian grid
!>@details This subroutine is for a 2-D (lat-lon) scalar field
Expand Down Expand Up @@ -322,18 +377,20 @@ end subroutine scalarspect_to_gaugrid
subroutine write_stoch_restart_atm(sfile)
!\callgraph
use netcdf
use stochy_namelist_def, only : do_sppt,do_shum,do_skeb,lndp_type
use stochy_namelist_def, only : do_sppt,do_shum,do_skeb,lndp_type,do_spp
implicit none
character(len=*) :: sfile
integer :: stochlun,k,n,isize,ierr
integer :: ncid,varid1a,varid1b,varid2a,varid2b,varid3a,varid3b,varid4a,varid4b
integer :: ncid,varid1a,varid1b,varid2a,varid2b,varid3a,varid3b,varid4a,varid4b,varid5a,varid5b
integer :: seed_dim_id,spec_dim_id,zt_dim_id,ztsfc_dim_id,np_dim_id,npsfc_dim_id
integer :: ztspp_dim_id,npspp_dim_id

include 'netcdf.inc'

if ( ( .NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0 ) ) return
if ( ( .NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0 ) .AND. (.NOT. do_spp)) return
stochlun=99
if (is_rootpe()) then
if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0 .OR. nlndp>0 ) then
if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0 .OR. nlndp>0 .OR. nspp>0 ) then
ierr=nf90_create(trim(sfile),cmode=NF90_CLOBBER,ncid=ncid)
ierr=NF90_PUT_ATT(ncid,NF_GLOBAL,"ntrunc",ntrunc)
call random_seed(size=isize) ! get seed size
Expand All @@ -347,6 +404,12 @@ subroutine write_stoch_restart_atm(sfile)
ierr=NF90_DEF_DIM(ncid,"n_var_lndp",n_var_lndp,ztsfc_dim_id)
ierr=NF90_PUT_ATT(ncid,ztsfc_dim_id,"long_name","number of sfc perturbation types")
endif
if (nspp .GT. 0) then
ierr=NF90_DEF_DIM(ncid,"num_patterns_spp",nspp,npspp_dim_id) ! should be 5
ierr=NF90_PUT_ATT(ncid,npspp_dim_id,"long_name","number of random patterns for spp)")
ierr=NF90_DEF_DIM(ncid,"n_var_spp",n_var_spp,ztspp_dim_id)
ierr=NF90_PUT_ATT(ncid,ztspp_dim_id,"long_name","number of spp perturbation types")
endif
ierr=NF90_DEF_DIM(ncid,"ndimspecx2",2*ndimspec,spec_dim_id)
ierr=NF90_PUT_ATT(ncid,spec_dim_id,"long_name","number of spectral cofficients")
if (do_sppt) then
Expand Down Expand Up @@ -375,6 +438,12 @@ subroutine write_stoch_restart_atm(sfile)
ierr=NF90_DEF_VAR(ncid,"sfcpert_spec",NF90_DOUBLE,(/spec_dim_id, ztsfc_dim_id, npsfc_dim_id/), varid4b)
ierr=NF90_PUT_ATT(ncid,varid4b,"long_name","spectral cofficients SHUM")
endif
if (nspp>0) then
ierr=NF90_DEF_VAR(ncid,"spp_seed",NF90_DOUBLE,(/seed_dim_id, ztspp_dim_id, npspp_dim_id/), varid5a)
ierr=NF90_PUT_ATT(ncid,varid5a,"long_name","random number seed for SPP")
ierr=NF90_DEF_VAR(ncid,"spp_spec",NF90_DOUBLE,(/spec_dim_id, ztspp_dim_id, npspp_dim_id/), varid5b)
ierr=NF90_PUT_ATT(ncid,varid5b,"long_name","spectral cofficients SPP")
endif
ierr=NF90_ENDDEF(ncid)
if (ierr .NE. 0) then
write(0,*) 'error creating stochastic restart file'
Expand Down Expand Up @@ -406,6 +475,11 @@ subroutine write_stoch_restart_atm(sfile)
enddo
enddo
endif
if (nspp > 0) then
do n=1,nspp
call write_pattern(rpattern_spp(n),ncid,1,n,varid5a,varid5b,.true.,ierr)
enddo
endif
if (is_rootpe() ) then
ierr=NF90_CLOSE(ncid)
if (ierr .NE. 0) then
Expand Down
Loading