From 34c9d8444829ad1f27946d5c014a0c2d91a7e66e Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Thu, 11 Nov 2021 09:36:46 -0700 Subject: [PATCH 01/12] Add SPP option to several physics parameterizations --- compns_stochy.F90 | 62 +++++++++++++++++++++++++++++++++ get_stochy_pattern.F90 | 30 +++++++++++++--- stochastic_physics.F90 | 76 +++++++++++++++++++++++++++++++++++------ stochy_data_mod.F90 | 71 ++++++++++++++++++++++++++++++++++++-- stochy_namelist_def.F90 | 11 +++++- 5 files changed, 232 insertions(+), 18 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 61c7f7cf..977cdaeb 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -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_spperts/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 @@ -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 @@ -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: @@ -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) @@ -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_spperts) +#else + rewind (nlunit) + open (unit=nlunit, file=fn_nml, action='READ', status='OLD', iostat=ios) + read(nlunit,nam_spperts) +#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 @@ -223,6 +247,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 @@ -310,6 +335,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. ! @@ -320,6 +380,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 ! diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index 3d235bee..81d1b791 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -5,12 +5,13 @@ module get_stochy_pattern_mod lats_node_a, lon_dim_a, len_trie_ls, & len_trio_ls, ls_dim, nodes, stochy_la2ga, & coslat_a, latg, latg2, levs, lonf, skeblevs - 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 @@ -364,18 +365,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 :: 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_master()) 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 @@ -389,6 +392,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 @@ -417,6 +426,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/), varid4a) + ierr=NF90_PUT_ATT(ncid,varid4a,"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/), varid4b) + ierr=NF90_PUT_ATT(ncid,varid4b,"long_name","spectral cofficients SPP") + endif ierr=NF90_ENDDEF(ncid) if (ierr .NE. 0) then write(0,*) 'error creating stochastic restart file' @@ -448,6 +463,13 @@ subroutine write_stoch_restart_atm(sfile) enddo enddo endif + if (nspp > 0) then + do n=1,nspp + do k=1,n_var_spp + call write_pattern(rpattern_spp(n),ncid,k,n,varid4a,varid4b,.true.,ierr) + enddo + enddo + endif if (is_master() ) then ierr=NF90_CLOSE(ncid) if (ierr .NE. 0) then diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 67e4033f..7c36ed4f 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -22,11 +22,15 @@ module stochastic_physics subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in, fn_nml, nlunit, & xlon,xlat, & do_sppt_in, do_shum_in, do_skeb_in, lndp_type_in, n_var_lndp_in, use_zmtnblck_out, skeb_npass_out, & - lndp_var_list_out, lndp_prt_list_out, ak, bk, nthreads, mpiroot, mpicomm, iret) + lndp_var_list_out, lndp_prt_list_out, & + n_var_spp_in, spp_var_list_out, spp_prt_list_out, do_spp_in, & + ak, bk, nthreads, mpiroot, mpicomm, iret) !\callgraph !use stochy_internal_state_moa use stochy_data_mod, only : init_stochdata,gg_lats,gg_lons,nsppt, & - rad2deg,INTTYP,wlon,rnlat,gis_stochy,vfact_skeb,vfact_sppt,vfact_shum,skeb_vpts,skeb_vwts,sl + rad2deg,INTTYP,wlon,rnlat,gis_stochy, & + vfact_skeb,vfact_sppt,vfact_shum,skeb_vpts,skeb_vwts,sl, & + nspp, vfact_spp use stochy_namelist_def use spectral_layout_mod,only:me,master,nodes,colrad_a,latg,lonf,skeblevs use mpi_wrapper, only : mpi_wrapper_initialize,mype,npes,is_master @@ -44,13 +48,16 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in character(len=*), intent(in) :: fn_nml real(kind=kind_dbl_prec), intent(in) :: xlon(:,:) real(kind=kind_dbl_prec), intent(in) :: xlat(:,:) -logical, intent(in) :: do_sppt_in, do_shum_in, do_skeb_in +logical, intent(in) :: do_sppt_in, do_shum_in, do_skeb_in ,do_spp_in integer, intent(in) :: lndp_type_in, n_var_lndp_in +integer, intent(in) :: n_var_spp_in real(kind=kind_dbl_prec), intent(in) :: ak(:), bk(:) logical, intent(out) :: use_zmtnblck_out integer, intent(out) :: skeb_npass_out character(len=3), dimension(max_n_var_lndp), intent(out) :: lndp_var_list_out real(kind=kind_dbl_prec), dimension(max_n_var_lndp), intent(out) :: lndp_prt_list_out +character(len=3), dimension(max_n_var_spp), intent(out) :: spp_var_list_out +real(kind=kind_dbl_prec), dimension(max_n_var_spp), intent(out) :: spp_prt_list_out ! Local variables @@ -58,7 +65,7 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in integer :: nblks,len real*8 :: PRSI(levs),PRSL(levs),dx real, allocatable :: skeb_vloc(:) -integer :: k,kflip,latghf,blk,k2 +integer :: k,kflip,latghf,blk,k2,v character*2::proc ! Initialize MPI and OpenMP @@ -118,13 +125,26 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in & ' namelist settings n_var_lndp in physics nml, and lndp_* in nam_sfcperts' iret = 20 return +else if (n_var_spp_in .ne. n_var_spp) then + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings n_var_spp in physics nml, and spp_* in nam_spperts' + write(0,*) 'n_var_spp, n_var_spp_in', n_var_spp, n_var_spp_in + iret = 20 + return +else if (do_spp_in.neqv.do_spp) then + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_spp and spp' + iret = 20 + return end if ! update remaining model configuration parameters from namelist use_zmtnblck_out=use_zmtnblck skeb_npass_out=skeb_npass lndp_var_list_out=lndp_var_list lndp_prt_list_out=lndp_prt_list -if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0) ) return +spp_var_list_out=spp_var_list +spp_prt_list_out=spp_prt_list +if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0) .AND. (.NOT. do_spp)) return allocate(sl(levs)) do k=1,levs sl(k)= 0.5*(ak(k)/101300.+bk(k)+ak(k+1)/101300.0+bk(k+1)) ! si are now sigmas @@ -205,6 +225,19 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in if (is_master()) print *,'shum vert profile',k,sl(k),vfact_shum(k) enddo endif +if (do_spp) then + allocate(vfact_spp(levs)) + do k=1,levs + if (sl(k) .lt. spp_sigtop1(1) .and. sl(k) .gt. spp_sigtop2(1)) then + vfact_spp(k) = (sl(k)-spp_sigtop2(1))/(spp_sigtop1(1)-spp_sigtop2(1)) + else if (sl(k) .lt. spp_sigtop2(1)) then + vfact_spp(k) = 0.0 + else + vfact_spp(k) = 1.0 + endif + if (is_master()) print *,'spp vert profile',k,sl(k),vfact_spp(k) + enddo +endif ! get interpolation weights ! define gaussian grid lats and lons latghf=latg/2 @@ -302,16 +335,17 @@ end subroutine init_stochastic_physics_ocn !allocates and polulates the necessary arrays subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, skebu_wts, & - skebv_wts, sfc_wts,nthreads) + skebv_wts, sfc_wts, spp_wts, nthreads) !\callgraph !use stochy_internal_state_mod use stochy_data_mod, only : nshum,rpattern_shum,rpattern_sppt,nsppt,rpattern_skeb,nskeb,& - gis_stochy,vfact_sppt,vfact_shum,vfact_skeb, rpattern_sfc, nlndp + gis_stochy,vfact_sppt,vfact_shum,vfact_skeb, rpattern_sfc, nlndp, & + rpattern_spp, nspp, vfact_spp use get_stochy_pattern_mod,only : get_random_pattern_scalar,get_random_pattern_vector, & get_random_pattern_sfc use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,sppt_logit, & - lndp_type, n_var_lndp + lndp_type, n_var_lndp, n_var_spp, do_spp, spp_stddev_cutoff, spp_prt_list use mpi_wrapper, only: is_master use spectral_layout_mod,only:me implicit none @@ -325,18 +359,19 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s real(kind=kind_dbl_prec), intent(inout) :: skebu_wts(:,:,:) real(kind=kind_dbl_prec), intent(inout) :: skebv_wts(:,:,:) real(kind=kind_dbl_prec), intent(inout) :: sfc_wts(:,:,:) +real(kind=kind_dbl_prec), intent(inout) :: spp_wts(:,:,:,:) integer, intent(in) :: nthreads real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:), tmpl_wts(:,:,:) !D-grid -integer :: k +integer :: k,v integer j,ierr,i integer :: nblks, blk, len, maxlen character*120 :: sfile character*6 :: STRFH logical :: do_advance_pattern -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. (n_var_spp .le. 0)) return ! Update number of threads in shared variables in spectral_layout_mod and set block-related variables nblks = size(blksz) @@ -410,6 +445,22 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s ENDDO if (is_master()) print*,'sfc_wts=',sfc_wts(1,1,:) deallocate(tmpl_wts) +endif +if (n_var_spp .GE. 1) then + DO v=1,n_var_spp + tmp_wts=0.0 + call get_random_pattern_scalar(rpattern_spp(v),nspp,gis_stochy,tmp_wts) + DO blk=1,nblks + len=blksz(blk) + DO k=1,levs + if (spp_stddev_cutoff(v).gt.0.0) then + spp_wts(blk,1:len,k,v)=MAX(MIN(tmp_wts(1:len,blk)*vfact_spp(k),spp_stddev_cutoff(v)),-1.0*spp_stddev_cutoff(v)) * spp_prt_list(v) + else + spp_wts(blk,1:len,k,v)=tmp_wts(1:len,blk)*vfact_spp(k) * spp_prt_list(v) + endif + ENDDO + ENDDO + ENDDO endif deallocate(tmp_wts) deallocate(tmpu_wts) @@ -457,6 +508,7 @@ end subroutine run_stochastic_physics_ocn subroutine finalize_stochastic_physics() use stochy_data_mod, only : nshum,rpattern_shum,rpattern_sppt,nsppt,rpattern_skeb,nskeb,& vfact_sppt,vfact_shum,vfact_skeb, skeb_vwts,skeb_vpts, & + rpattern_spp, vfact_spp, nspp, & rpattern_sfc, nlndp,gg_lats,gg_lons,sl,skebu_save,skebv_save,gis_stochy use spectral_layout_mod, only : lat1s_h,lat1s_a ,lon_dims_a,wgt_a,sinlat_a,coslat_a,colrad_a,wgtcs_a,rcs2_a,lats_nodes_h,global_lats_h implicit none @@ -483,6 +535,10 @@ subroutine finalize_stochastic_physics() if (nlndp > 0) then if (allocated(rpattern_sfc)) deallocate(rpattern_sfc) endif + if (nspp > 0) then + if (allocated(rpattern_spp)) deallocate(rpattern_spp) + if (allocated(vfact_spp)) deallocate(vfact_spp) + endif deallocate(lat1s_a) deallocate(lon_dims_a) diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index 9247c4f9..3fd4ce01 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -23,16 +23,17 @@ module stochy_data_mod public :: init_stochdata,init_stochdata_ocn type(random_pattern), public, save, allocatable, dimension(:) :: & - rpattern_sppt,rpattern_shum,rpattern_skeb, rpattern_sfc,rpattern_epbl1,rpattern_epbl2,rpattern_ocnsppt + rpattern_sppt,rpattern_shum,rpattern_skeb, rpattern_sfc,rpattern_epbl1,rpattern_epbl2,rpattern_ocnsppt,rpattern_spp integer, public :: nepbl=0 integer, public :: nocnsppt=0 integer, public :: nsppt=0 integer, public :: nshum=0 integer, public :: nskeb=0 integer, public :: nlndp=0 ! this is the number of different patterns (determined by the tau/lscale input) + integer, public :: nspp =0 ! this is the number of different patterns (determined by the tau/lscale input) real*8, public,allocatable :: sl(:) - real(kind=kind_dbl_prec),public, allocatable :: vfact_sppt(:),vfact_shum(:),vfact_skeb(:) + real(kind=kind_dbl_prec),public, allocatable :: vfact_sppt(:),vfact_shum(:),vfact_skeb(:),vfact_spp(:) real(kind=kind_dbl_prec),public, allocatable :: skeb_vwts(:,:),skeb_vpts(:,:) real(kind=kind_dbl_prec),public, allocatable :: gg_lats(:),gg_lons(:) real(kind=kind_dbl_prec),public :: wlon,rnlat,rad2deg @@ -75,7 +76,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) call compns_stochy (me,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,delt,iret) if (iret/=0) return ! need to make sure that non-zero irets are being trapped. if(is_master()) print*,'in init stochdata',nodes,lat_s - 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 ! initialize the specratl pattern generatore (including gaussian grid decomposition) ! if (nodes.GE.lat_s/2) then ! lat_s=(int(nodes/12)+1)*24 @@ -123,12 +124,15 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) !enddo if (n_var_lndp>0) nlndp=1 if (is_master()) print *,' nlndp = ', nlndp + if (n_var_spp>0) nspp=n_var_spp + if (is_master()) print *,' nspp = ', nspp if (nsppt > 0) allocate(rpattern_sppt(nsppt)) if (nshum > 0) allocate(rpattern_shum(nshum)) if (nskeb > 0) allocate(rpattern_skeb(nskeb)) ! mg, sfc perts if (nlndp > 0) allocate(rpattern_sfc(nlndp)) + if (nspp > 0) allocate(rpattern_spp(nspp)) ! if stochini is true, then read in pattern from a file if (is_master()) then @@ -433,6 +437,67 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) enddo ! k, n_var_lndp enddo ! n, nlndp endif ! nlndp > 0 + if (nspp > 0) then + if (is_master()) then + print *, 'Initialize random pattern for SPP-PERTS' + if (stochini) then + ierr=NF90_INQ_VARID(stochlun,"spppert_seed", varid1) + if (ierr .NE. 0) then + write(0,*) 'error inquring SPP-PERTS seed' + iret = ierr + return + end if + ierr=NF90_INQ_VARID(stochlun,"ppcpert_spec", varid2) + if (ierr .NE. 0) then + write(0,*) 'error inquring SPP-PERTS spec' + iret = ierr + return + end if + endif + endif + ones = 1. + call patterngenerator_init(spp_lscale(1:nspp),delt,spp_tau(1:nspp),ones(1:nspp),iseed_spp,rpattern_spp, & + lonf,latg,jcap,gis_stochy%ls_node,nspp,n_var_spp,0,new_lscale) + do n=1,nspp + if (is_master()) print *, 'Initialize random pattern for SPP PERTS' + do k=1,n_var_spp + nspinup = spinup_efolds*spp_tau(n)/delt + if (stochini) then + call read_pattern(rpattern_spp(n),jcapin,stochlun,k,n,varid1,varid2,.true.,ierr) + if (ierr .NE. 0) then + write(0,*) 'error reading SPP pattern' + iret = ierr + return + endif + if (is_master()) print *, 'spp pattern read',n,k,minval(rpattern_spp(n)%spec_o(:,:,k)), maxval(rpattern_spp(n)%spec_o(:,:,k)) + else + call getnoise(rpattern_spp(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_spp(n)%spec_e(nn,1,k)=noise_e(nn,1) + rpattern_spp(n)%spec_e(nn,2,k)=noise_e(nn,2) + nm = rpattern_spp(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_spp(n)%spec_e(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_e(nn,1,k)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_e(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_e(nn,2,k)*rpattern_spp(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_spp(n)%spec_o(nn,1,k)=noise_o(nn,1) + rpattern_spp(n)%spec_o(nn,2,k)=noise_o(nn,2) + nm = rpattern_sfc(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_sfc(n)%spec_o(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,1,k)*rpattern_spp(n)%varspectrum(nm) + rpattern_sfc(n)%spec_o(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,2,k)*rpattern_spp(n)%varspectrum(nm) + enddo + do nn=1,nspinup + call patterngenerator_advance(rpattern_spp(n),k,.false.) + enddo + if (is_master()) print *, 'spp pattern initialized, ',n, k, minval(rpattern_spp(n)%spec_o(:,:,k)), maxval(rpattern_spp(n)%spec_o(:,:,k)) + endif ! stochini + enddo ! k, n_var_spp + enddo ! n, nspp + endif ! nspp > 0 + + if (is_master() .and. stochini) CLOSE(stochlun) deallocate(noise_e,noise_o) end subroutine init_stochdata diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index adc3f022..f5b31ece 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -11,6 +11,7 @@ module stochy_namelist_def public integer, parameter :: max_n_var_lndp = 6 ! must match value used in GFS_typedefs + integer, parameter :: max_n_var_spp = 6 ! must match value used in GFS_typedefs integer nssppt,nsshum,nsepbl,nsocnsppt,nsskeb,lon_s,lat_s,ntrunc ! pjp stochastic phyics @@ -30,7 +31,7 @@ module stochy_namelist_def integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_epbl,iseed_ocnsppt,iseed_epbl2 logical stochini,sppt_logit,new_lscale logical use_zmtnblck - logical do_shum,do_sppt,do_skeb,pert_epbl,do_ocnsppt + logical do_shum,do_sppt,do_skeb,pert_epbl,do_ocnsppt,do_spp real(kind=kind_dbl_prec), dimension(5) :: lndp_lscale,lndp_tau integer n_var_lndp @@ -39,4 +40,12 @@ module stochy_namelist_def character(len=3), dimension(max_n_var_lndp) :: lndp_var_list real(kind=kind_dbl_prec), dimension(max_n_var_lndp) :: lndp_prt_list + real(kind=kind_dbl_prec), dimension(max_n_var_spp) :: spp_lscale & + & , spp_tau,spp_stddev_cutoff & + & , spp_sigtop1, spp_sigtop2 + integer n_var_spp + integer(8),dimension(max_n_var_spp) ::iseed_spp + character(len=3), dimension(max_n_var_spp) :: spp_var_list + real(kind=kind_dbl_prec), dimension(max_n_var_spp) :: spp_prt_list + end module stochy_namelist_def From 7f1159e3c0e90f3c51c9e08d11225cbd57eaaa6d Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Mon, 22 Nov 2021 09:37:31 -0700 Subject: [PATCH 02/12] bug fix --- stochy_data_mod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index 3fd4ce01..e22eeb38 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -483,10 +483,10 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) do nn=1,len_trio_ls rpattern_spp(n)%spec_o(nn,1,k)=noise_o(nn,1) rpattern_spp(n)%spec_o(nn,2,k)=noise_o(nn,2) - nm = rpattern_sfc(n)%idx_o(nn) + nm = rpattern_spp(n)%idx_o(nn) if (nm .eq. 0) cycle - rpattern_sfc(n)%spec_o(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,1,k)*rpattern_spp(n)%varspectrum(nm) - rpattern_sfc(n)%spec_o(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,2,k)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_o(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,1,k)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_o(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,2,k)*rpattern_spp(n)%varspectrum(nm) enddo do nn=1,nspinup call patterngenerator_advance(rpattern_spp(n),k,.false.) From 0b868b1592475ffa804cb3006a4070350b214ba2 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Tue, 23 Nov 2021 12:59:27 -0700 Subject: [PATCH 03/12] merge fix --- stochy_data_mod.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index 6782fad3..f80f1cd0 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -436,7 +436,6 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) do n=1,nspp if (is_rootpe()) print *, 'Initialize random pattern for SPP PERTS' do k=1,n_var_spp - nspinup = spinup_efolds*spp_tau(n)/delt if (stochini) then call read_pattern(rpattern_spp(n),jcapin,stochlun,k,n,varid1,varid2,.true.,ierr) if (ierr .NE. 0) then @@ -463,9 +462,6 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) rpattern_spp(n)%spec_o(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,1,k)*rpattern_spp(n)%varspectrum(nm) rpattern_spp(n)%spec_o(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,2,k)*rpattern_spp(n)%varspectrum(nm) enddo - do nn=1,nspinup - call patterngenerator_advance(rpattern_spp(n),k,.false.) - enddo if (is_rootpe()) print *, 'spp pattern initialized, ',n, k, minval(rpattern_spp(n)%spec_o(:,:,k)), maxval(rpattern_spp(n)%spec_o(:,:,k)) endif ! stochini enddo ! k, n_var_spp From 40ec485d17f1cd3045edaa87897e2c65a2260c68 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Tue, 23 Nov 2021 13:34:53 -0700 Subject: [PATCH 04/12] merge fixes --- stochastic_physics.F90 | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 69a12972..57d0327a 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -93,7 +93,6 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in ! replace INTTYP=0 ! bilinear interpolation call init_stochdata(levs,dtp,input_nml_file_in,fn_nml,nlunit,iret) -print*,'back from init stochdata',iret if (iret .ne. 0) return ! check namelist entries for consistency if (do_sppt_in.neqv.do_sppt) then @@ -235,7 +234,7 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in else vfact_spp(k) = 1.0 endif - if (is_master()) print *,'spp vert profile',k,sl(k),vfact_spp(k) + if (is_rootpe()) print *,'spp vert profile',k,sl(k),vfact_spp(k) enddo endif ! get interpolation weights @@ -342,14 +341,8 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s use get_stochy_pattern_mod,only : get_random_pattern_scalar,get_random_pattern_vector, & get_random_pattern_sfc use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,sppt_logit, & -<<<<<<< HEAD lndp_type, n_var_lndp, n_var_spp, do_spp, spp_stddev_cutoff, spp_prt_list -use mpi_wrapper, only: is_master -use spectral_layout_mod,only:me -======= - lndp_type, n_var_lndp use mpi_wrapper, only: is_rootpe ->>>>>>> 4e8a7cff92ee158104963dddea5ae335071cbbb9 implicit none ! Interface variables @@ -372,12 +365,8 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s character*120 :: sfile character*6 :: STRFH logical :: do_advance_pattern -<<<<<<< HEAD if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0 ) .AND. (n_var_spp .le. 0)) return -======= -if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0 ) ) return ->>>>>>> 4e8a7cff92ee158104963dddea5ae335071cbbb9 ! Update number of threads in shared variables in spectral_layout_mod and set block-related variables nblks = size(blksz) From 86755acf295cda84687ea87f0145af4d0065e7f3 Mon Sep 17 00:00:00 2001 From: llpcarson Date: Tue, 30 Nov 2021 21:14:41 +0000 Subject: [PATCH 05/12] add iostat to new nml read --- compns_stochy.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index cc358782..177937c2 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -163,7 +163,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) #endif #ifdef INTERNAL_FILE_NML - read(input_nml_file, nml=nam_spperts) + read(input_nml_file, nml=nam_spperts, iostat=ios) #else rewind (nlunit) open (unit=nlunit, file=fn_nml, action='READ', status='OLD', iostat=ios) From d4446d9a8030ae69b041d7d9bc3da9b5759a6338 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Fri, 11 Feb 2022 02:52:31 +0000 Subject: [PATCH 06/12] Fix application of pattern generator when using SPP. --- get_stochy_pattern.F90 | 70 +++++++++++++++++++++++++++++++++++++----- stochastic_physics.F90 | 15 ++++----- stochy_data_mod.F90 | 25 +++++++-------- 3 files changed, 81 insertions(+), 29 deletions(-) diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index b493f8c6..82fdcd69 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -21,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. @@ -275,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 @@ -327,7 +381,7 @@ subroutine write_stoch_restart_atm(sfile) 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 @@ -385,10 +439,10 @@ subroutine write_stoch_restart_atm(sfile) 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/), varid4a) - ierr=NF90_PUT_ATT(ncid,varid4a,"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/), varid4b) - ierr=NF90_PUT_ATT(ncid,varid4b,"long_name","spectral cofficients SPP") + 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 @@ -424,7 +478,7 @@ subroutine write_stoch_restart_atm(sfile) if (nspp > 0) then do n=1,nspp do k=1,n_var_spp - call write_pattern(rpattern_spp(n),ncid,k,n,varid4a,varid4b,.true.,ierr) + call write_pattern(rpattern_spp(n),ncid,k,n,varid5a,varid5b,.true.,ierr) enddo enddo endif diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 57d0327a..b1306cba 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -339,7 +339,7 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s gis_stochy,vfact_sppt,vfact_shum,vfact_skeb, rpattern_sfc, nlndp, & rpattern_spp, nspp, vfact_spp use get_stochy_pattern_mod,only : get_random_pattern_scalar,get_random_pattern_vector, & - get_random_pattern_sfc + get_random_pattern_sfc,get_random_pattern_spp use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,sppt_logit, & lndp_type, n_var_lndp, n_var_spp, do_spp, spp_stddev_cutoff, spp_prt_list use mpi_wrapper, only: is_rootpe @@ -357,7 +357,7 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s real(kind=kind_dbl_prec), intent(inout) :: spp_wts(:,:,:,:) integer, intent(in) :: nthreads -real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:), tmpl_wts(:,:,:) +real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:),tmpl_wts(:,:,:),tmp_spp_wts(:,:,:) !D-grid integer :: k,v integer j,ierr,i @@ -441,20 +441,21 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s deallocate(tmpl_wts) endif if (n_var_spp .GE. 1) then - DO v=1,n_var_spp - tmp_wts=0.0 - call get_random_pattern_scalar(rpattern_spp(v),nspp,gis_stochy,tmp_wts) + allocate(tmp_spp_wts(gis_stochy%nx,gis_stochy%ny,n_var_spp)) + call get_random_pattern_spp(rpattern_spp,nspp,gis_stochy,tmp_spp_wts) + DO v=1,n_var_spp DO blk=1,nblks len=blksz(blk) DO k=1,levs if (spp_stddev_cutoff(v).gt.0.0) then - spp_wts(blk,1:len,k,v)=MAX(MIN(tmp_wts(1:len,blk)*vfact_spp(k),spp_stddev_cutoff(v)),-1.0*spp_stddev_cutoff(v)) * spp_prt_list(v) + spp_wts(blk,1:len,k,v)=MAX(MIN(tmp_spp_wts(1:len,blk,v)*vfact_spp(k),spp_stddev_cutoff(v)),-1.0*spp_stddev_cutoff(v)) * spp_prt_list(v) else - spp_wts(blk,1:len,k,v)=tmp_wts(1:len,blk)*vfact_spp(k) * spp_prt_list(v) + spp_wts(blk,1:len,k,v)=tmp_spp_wts(1:len,blk,v)*vfact_spp(k) * spp_prt_list(v) endif ENDDO ENDDO ENDDO + deallocate(tmp_spp_wts) endif deallocate(tmp_wts) deallocate(tmpu_wts) diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index f80f1cd0..837587e3 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -435,40 +435,37 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) lonf,latg,jcap,gis_stochy%ls_node,nspp,n_var_spp,0,new_lscale) do n=1,nspp if (is_rootpe()) print *, 'Initialize random pattern for SPP PERTS' - do k=1,n_var_spp if (stochini) then - call read_pattern(rpattern_spp(n),jcapin,stochlun,k,n,varid1,varid2,.true.,ierr) + call read_pattern(rpattern_spp(n),jcapin,stochlun,1,n,varid1,varid2,.true.,ierr) if (ierr .NE. 0) then write(0,*) 'error reading SPP pattern' iret = ierr return endif - if (is_rootpe()) print *, 'spp pattern read',n,k,minval(rpattern_spp(n)%spec_o(:,:,k)), maxval(rpattern_spp(n)%spec_o(:,:,k)) + if (is_rootpe()) print *, 'spp pattern read',n,1,minval(rpattern_spp(n)%spec_o(:,:,1)), maxval(rpattern_spp(n)%spec_o(:,:,1)) else call getnoise(rpattern_spp(n),noise_e,noise_o) do nn=1,len_trie_ls - rpattern_spp(n)%spec_e(nn,1,k)=noise_e(nn,1) - rpattern_spp(n)%spec_e(nn,2,k)=noise_e(nn,2) + rpattern_spp(n)%spec_e(nn,1,1)=noise_e(nn,1) + rpattern_spp(n)%spec_e(nn,2,1)=noise_e(nn,2) nm = rpattern_spp(n)%idx_e(nn) if (nm .eq. 0) cycle - rpattern_spp(n)%spec_e(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_e(nn,1,k)*rpattern_spp(n)%varspectrum(nm) - rpattern_spp(n)%spec_e(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_e(nn,2,k)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_e(nn,1,1) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_e(nn,1,1)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_e(nn,2,1) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_e(nn,2,1)*rpattern_spp(n)%varspectrum(nm) enddo do nn=1,len_trio_ls - rpattern_spp(n)%spec_o(nn,1,k)=noise_o(nn,1) - rpattern_spp(n)%spec_o(nn,2,k)=noise_o(nn,2) + rpattern_spp(n)%spec_o(nn,1,1)=noise_o(nn,1) + rpattern_spp(n)%spec_o(nn,2,1)=noise_o(nn,2) nm = rpattern_spp(n)%idx_o(nn) if (nm .eq. 0) cycle - rpattern_spp(n)%spec_o(nn,1,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,1,k)*rpattern_spp(n)%varspectrum(nm) - rpattern_spp(n)%spec_o(nn,2,k) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,2,k)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_o(nn,1,1) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,1,1)*rpattern_spp(n)%varspectrum(nm) + rpattern_spp(n)%spec_o(nn,2,1) = rpattern_spp(n)%stdev*rpattern_spp(n)%spec_o(nn,2,1)*rpattern_spp(n)%varspectrum(nm) enddo - if (is_rootpe()) print *, 'spp pattern initialized, ',n, k, minval(rpattern_spp(n)%spec_o(:,:,k)), maxval(rpattern_spp(n)%spec_o(:,:,k)) + if (is_rootpe()) print *, 'spp pattern initialized, ',n, 1, minval(rpattern_spp(n)%spec_o(:,:,1)), maxval(rpattern_spp(n)%spec_o(:,:,1)) endif ! stochini - enddo ! k, n_var_spp enddo ! n, nspp endif ! nspp > 0 - if (is_rootpe() .and. stochini) CLOSE(stochlun) deallocate(noise_e,noise_o) end subroutine init_stochdata From 7a74cddedb10334a1fc43a4e25c2e361017fab1d Mon Sep 17 00:00:00 2001 From: jeff beck Date: Fri, 11 Feb 2022 04:18:59 +0000 Subject: [PATCH 07/12] Fix write call. --- get_stochy_pattern.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index 82fdcd69..1fd00bec 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -477,9 +477,7 @@ subroutine write_stoch_restart_atm(sfile) endif if (nspp > 0) then do n=1,nspp - do k=1,n_var_spp - call write_pattern(rpattern_spp(n),ncid,k,n,varid5a,varid5b,.true.,ierr) - enddo + call write_pattern(rpattern_spp(n),ncid,1,n,varid5a,varid5b,.true.,ierr) enddo endif if (is_rootpe() ) then From 095cd2885008db05a8adef1643bcd9282158ff0c Mon Sep 17 00:00:00 2001 From: jeff beck Date: Sat, 12 Feb 2022 18:17:15 +0000 Subject: [PATCH 08/12] Fix to stochastic pattern scaling - match WRF-based SPP code --- stochastic_physics.F90 | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index b1306cba..7dde1398 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -65,7 +65,7 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in integer :: nblks,len real*8 :: PRSI(levs),PRSL(levs),dx real, allocatable :: skeb_vloc(:) -integer :: k,kflip,latghf,blk,k2,v +integer :: k,kflip,latghf,blk,k2,v,i character*2::proc ! Initialize MPI and OpenMP @@ -447,11 +447,16 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s DO blk=1,nblks len=blksz(blk) DO k=1,levs - if (spp_stddev_cutoff(v).gt.0.0) then - spp_wts(blk,1:len,k,v)=MAX(MIN(tmp_spp_wts(1:len,blk,v)*vfact_spp(k),spp_stddev_cutoff(v)),-1.0*spp_stddev_cutoff(v)) * spp_prt_list(v) - else - spp_wts(blk,1:len,k,v)=tmp_spp_wts(1:len,blk,v)*vfact_spp(k) * spp_prt_list(v) - endif + spp_wts(blk,1:len,k,v)=tmp_spp_wts(1:len,blk,v)*vfact_spp(k) + DO i=1,len + if (spp_wts(blk,i,k,v) .GT. spp_stddev_cutoff(v)*spp_prt_list(v)) then + spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*spp_stddev_cutoff(v)*spp_prt_list(v) + endif + if (spp_wts(blk,i,k,v) .LT. -1.0*spp_stddev_cutoff(v)*spp_prt_list(v)) then + spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*spp_stddev_cutoff(v)*spp_prt_list(v) + endif + !spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*vfact_spp(k) + ENDDO ENDDO ENDDO ENDDO From 5d0d2391d5690d038edb76a8fec88c7821e1e049 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Sat, 12 Feb 2022 21:53:35 +0000 Subject: [PATCH 09/12] Fix cutoff thresholding --- stochastic_physics.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 7dde1398..4027495f 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -450,10 +450,10 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s spp_wts(blk,1:len,k,v)=tmp_spp_wts(1:len,blk,v)*vfact_spp(k) DO i=1,len if (spp_wts(blk,i,k,v) .GT. spp_stddev_cutoff(v)*spp_prt_list(v)) then - spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*spp_stddev_cutoff(v)*spp_prt_list(v) + spp_wts(blk,i,k,v)=spp_stddev_cutoff(v)*spp_prt_list(v) endif if (spp_wts(blk,i,k,v) .LT. -1.0*spp_stddev_cutoff(v)*spp_prt_list(v)) then - spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*spp_stddev_cutoff(v)*spp_prt_list(v) + spp_wts(blk,i,k,v)=-1.0*spp_stddev_cutoff(v)*spp_prt_list(v) endif !spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*vfact_spp(k) ENDDO From 659c621ffbbc3802e4603bd9e3afb7aece3655c7 Mon Sep 17 00:00:00 2001 From: jeff beck Date: Sun, 13 Feb 2022 01:12:15 +0000 Subject: [PATCH 10/12] Final spp_wts configuration. --- stochastic_physics.F90 | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 4027495f..83ab324c 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -447,16 +447,11 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s DO blk=1,nblks len=blksz(blk) DO k=1,levs - spp_wts(blk,1:len,k,v)=tmp_spp_wts(1:len,blk,v)*vfact_spp(k) - DO i=1,len - if (spp_wts(blk,i,k,v) .GT. spp_stddev_cutoff(v)*spp_prt_list(v)) then - spp_wts(blk,i,k,v)=spp_stddev_cutoff(v)*spp_prt_list(v) - endif - if (spp_wts(blk,i,k,v) .LT. -1.0*spp_stddev_cutoff(v)*spp_prt_list(v)) then - spp_wts(blk,i,k,v)=-1.0*spp_stddev_cutoff(v)*spp_prt_list(v) - endif - !spp_wts(blk,i,k,v)=spp_wts(blk,i,k,v)*vfact_spp(k) - ENDDO + if (spp_stddev_cutoff(v).gt.0.0) then + spp_wts(blk,1:len,k,v)=MAX(MIN(tmp_spp_wts(1:len,blk,v)*vfact_spp(k),spp_stddev_cutoff(v)),-1.0*spp_stddev_cutoff(v))*spp_prt_list(v) + else + spp_wts(blk,1:len,k,v)=tmp_spp_wts(1:len,blk,v)*vfact_spp(k)*spp_prt_list(v) + endif ENDDO ENDDO ENDDO From 079820e43fe5a64a42b736b03b4e34559906d8ee Mon Sep 17 00:00:00 2001 From: jeff beck Date: Tue, 15 Feb 2022 23:34:26 +0000 Subject: [PATCH 11/12] Update SPP namelist entry --- compns_stochy.F90 | 2 +- stochastic_physics.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 177937c2..74e7ee32 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -66,7 +66,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & lndp_tau,lndp_lscale ! For SPP physics parameterization perterbations - namelist /nam_spperts/spp_var_list, spp_prt_list, iseed_spp, & + 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) diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 83ab324c..70ba6037 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -122,7 +122,7 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in return else if (n_var_spp_in .ne. n_var_spp) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & - & ' namelist settings n_var_spp in physics nml, and spp_* in nam_spperts' + & ' namelist settings n_var_spp in physics nml, and spp_* in nam_sppperts' write(0,*) 'n_var_spp, n_var_spp_in', n_var_spp, n_var_spp_in iret = 20 return From 450508d1c9ca5110da9c4079b54e8515eca2fa3c Mon Sep 17 00:00:00 2001 From: jeff beck Date: Wed, 16 Feb 2022 00:44:17 +0000 Subject: [PATCH 12/12] Update SPP namelist entry --- compns_stochy.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 74e7ee32..15b68bd3 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -163,11 +163,11 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) #endif #ifdef INTERNAL_FILE_NML - read(input_nml_file, nml=nam_spperts, iostat=ios) + 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_spperts) + read(nlunit,nam_sppperts) #endif if (me == 0) then