From 5c9c8772049fc7d074bf030f09f3724b9e12632c Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Sat, 25 May 2024 11:03:06 -0600 Subject: [PATCH 1/8] Fix normalization --- compns_stochy.F90 | 199 +++++++++++++++++++++------- stochastic_physics.F90 | 146 +++++++++++++++------ stochy_data_mod.F90 | 151 ++++++++++++++++++++-- stochy_namelist_def.F90 | 21 ++- stochy_patterngenerator.F90 | 249 ++++++++++++++++++++++++++++-------- 5 files changed, 611 insertions(+), 155 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 194730e..06be05c 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -56,15 +56,19 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! namelist /nam_stochy/ntrunc,lon_s,lat_s,sppt,sppt_tau,sppt_lscale,sppt_logit, & - iseed_shum,iseed_sppt,shum,shum_tau,& - shum_lscale,stochini,skeb_varspect_opt,sppt_sfclimit, & - skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & - skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& - shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & - epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt - namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & - lndp_tau,lndp_lscale + iseed_shum,iseed_sppt,shum,shum_tau, & + shum_lscale,stochini,skeb_varspect_opt,sppt_sfclimit, & + skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & + skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2, & + shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & + epbl,epbl_lscale,epbl_tau,iseed_epbl, & + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt, & + ocnskeb,ocnskeb_lscale,ocnskeb_tau,iseed_ocnskeb + namelist /nam_sfcperts/lndp_type,lndp_model_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 @@ -79,12 +83,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 @@ -95,9 +102,18 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! perturbations are assigned once at the start of the forecast ! LNDP_TYPE = 2 ! this is the newer land pert scheme, introduced and tested for impact on UFS/GDAS cycling stsyem -! perturbations are assigned at each time step (for state variables), or each time parameters are updated -! and the perturbations evolve over time. +! see https://journals.ametsoc.org/view/journals/hydr/22/8/JHM-D-21-0016.1.xml lndp_type = 0 ! +! LNDP_MODEL_TYPE +! integer indicating the model type for applying perturbations for lndp_type=2 scheme. +! 1 - global model +! (cycling of prognostic variables between DA cycles, +! parameters may be periodically updated during forecast) +! 2 - regional model +! (short foreast, prognotic variables re-initialized each forecast, +! parameters not updated during forecast) +! 3 - special case to apply perturbations only at start of forecast. + lndp_model_type = 0 ! lndp_lscale = -999. ! length scales lndp_tau = -999. ! time scales iseed_lndp = 0 ! random seeds (if 0 use system clock) @@ -111,7 +127,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) shum_tau = -999. skeb_tau = -999. skeb_vdof = 5 ! proxy for vertical correlation, 5 is close to 40 passes of the 1-2-1 filter in the GFS - skebnorm = 0 ! 0 - random pattern is stream function, 1- pattern is kenorm, 2- pattern is vorticity + skebnorm = 0 ! 0 - random pattern is stream function, 1- pattern is kenorm, 2- pattern is vorticity, 3- normalized following Berner et al. 2009 sppt_lscale = -999. ! length scales shum_lscale = -999. skeb_lscale = -999. @@ -125,6 +141,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 +151,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 +172,20 @@ 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) + close(nlunit) +#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 @@ -179,26 +214,32 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) endif ENDIF ! compute frequencty to estimate dissipation timescale - IF (skebint == 0.) skebint=deltim - nsskeb=nint(skebint/deltim) ! skebint in seconds - IF(nsskeb<=0 .or. abs(nsskeb-skebint/deltim)>tol) THEN - WRITE(0,*) "SKEB interval is invalid",skebint - iret=9 - return - ENDIF - IF (spptint == 0.) spptint=deltim - nssppt=nint(spptint/deltim) ! spptint in seconds - IF(nssppt<=0 .or. abs(nssppt-spptint/deltim)>tol) THEN - WRITE(0,*) "SPPT interval is invalid",spptint - iret=9 - return - ENDIF - IF (shumint == 0.) shumint=deltim - nsshum=nint(shumint/deltim) ! shumint in seconds - IF(nsshum<=0 .or. abs(nsshum-shumint/deltim)>tol) THEN - WRITE(0,*) "SHUM interval is invalid",shumint - iret=9 - return + IF (do_skeb) THEN + IF (skebint == 0.) skebint=deltim + nsskeb=nint(skebint/deltim) ! skebint in seconds + IF(nsskeb<=0 .or. abs(nsskeb-skebint/deltim)>tol) THEN + WRITE(0,*) "SKEB interval is invalid",skebint + iret=9 + return + ENDIF + ENDIF + IF (do_sppt) THEN + IF (spptint == 0.) spptint=deltim + nssppt=nint(spptint/deltim) ! spptint in seconds + IF(nssppt<=0 .or. abs(nssppt-spptint/deltim)>tol) THEN + WRITE(0,*) "SPPT interval is invalid",spptint + iret=9 + return + ENDIF + ENDIF + IF (do_shum) THEN + IF (shumint == 0.) shumint=deltim + nsshum=nint(shumint/deltim) ! shumint in seconds + IF(nsshum<=0 .or. abs(nsshum-shumint/deltim)>tol) THEN + WRITE(0,*) "SHUM interval is invalid",shumint + iret=9 + return + ENDIF ENDIF !calculate ntrunc if not supplied if (ntrunc .LT. 1) then @@ -211,6 +252,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 @@ -245,13 +287,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) cycle else n_var_lndp=n_var_lndp+1 - lndp_var_list( n_var_lndp) = lndp_var_list(k) ! for lndp_type==2: - ! for state variables, unit is pert per hour - ! for parmaters, no time dimension in unit - ! since perturbations do not accumulate - ! (i.e., global_cycle overwrites the paramaters - ! each time it's called, so any previous perturbations - ! are lost). + lndp_var_list( n_var_lndp) = lndp_var_list(k) lndp_prt_list( n_var_lndp) = lndp_prt_list(k) endif enddo @@ -265,7 +301,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) if (lndp_type==1) then if (me==0) print*, & - 'lndp_type=1, land perturbations will be applied to selected paramaters, using older scheme designed for S2S fcst spread' + 'lndp_type=1, land perturbations will be applied to selected paramaters, using older scheme designed for S2S fcst spread with the noah LSM' ! sanity-check requested input do k =1,n_var_lndp select case (lndp_var_list(k)) @@ -280,6 +316,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) elseif(lndp_type==2) then if (me==0) print*, & 'land perturbations will be applied to selected paramaters, using newer scheme designed for DA ens spread' + ! check requested parameters have been coded. + ! note, Noah-MP specific checks will be done later (since need to know lsm type) do k =1,n_var_lndp select case (lndp_var_list(k)) case('vgf','smc','stc','alb', 'sal','emi','zol') @@ -290,6 +328,11 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) return end select enddo + if ( (lndp_model_type < 1) .or. (lndp_model_type > 3) ) then + print*, 'ERROR: for lndp_type=2, must have lndp_model_type = 1,2,3' + iret = 10 + return + endif endif case default @@ -298,6 +341,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. ! @@ -307,7 +385,10 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) print *, ' do_shum : ', do_shum print *, ' do_skeb : ', do_skeb print *, ' lndp_type : ', lndp_type + print *, ' lndp_model_type : ', lndp_model_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 ! @@ -359,15 +440,16 @@ subroutine compns_stochy_ocn (deltim,iret) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! namelist /nam_stochy/ntrunc,lon_s,lat_s,sppt,sppt_tau,sppt_lscale,sppt_logit, & - iseed_shum,iseed_sppt,shum,shum_tau, & - shum_lscale,stochini,skeb_varspect_opt,sppt_sfclimit, & - skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & - skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& - shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & - epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt + iseed_shum,iseed_sppt,shum,shum_tau, & + shum_lscale,stochini,skeb_varspect_opt,sppt_sfclimit, & + skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & + skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2, & + shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & + epbl,epbl_lscale,epbl_tau,iseed_epbl, & + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt, & + ocnskeb,ocnskeb_lscale,ocnskeb_tau,iseed_ocnskeb - namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & + namelist /nam_sfcperts/lndp_type,lndp_model_type,lndp_var_list, lndp_prt_list, iseed_lndp, & lndp_tau,lndp_lscale @@ -382,22 +464,29 @@ subroutine compns_stochy_ocn (deltim,iret) ! (each is an array of length 5) epbl = -999. ! stochastic physics tendency amplitude ocnsppt = -999. ! stochastic physics tendency amplitude + ocnskeb = -999. ! stochastic physics tendency amplitude ! logicals pert_epbl = .false. do_ocnsppt = .false. + do_ocnskeb = .false. new_lscale = .false. epblint = 0 ocnspptint = 0 + ocnskebint = 0 epbl_tau = -999. ! time scales ocnsppt_tau = -999. ! time scales + ocnskeb_tau = -999. ! time scales epbl_lscale = -999. ! length scales ocnsppt_lscale = -999. ! length scales + ocnskeb_lscale = -999. ! length scales iseed_epbl = 0 ! random seeds (if 0 use system clock) iseed_epbl2 = 0 ! random seeds (if 0 use system clock) iseed_ocnsppt = 0 ! random seeds (if 0 use system clock) + iseed_ocnskeb = 0 ! random seeds (if 0 use system clock) rewind (nlunit) open (unit=nlunit, file='input.nml', action='READ', status='OLD', iostat=ios) read(nlunit,nam_stochy) + close(nlunit) if (mpp_pe()==mpp_root_pe()) then print *,' in compns_stochy_ocn' @@ -410,6 +499,9 @@ subroutine compns_stochy_ocn (deltim,iret) IF (ocnsppt(1) > 0 ) THEN do_ocnsppt=.true. ENDIF + IF (ocnskeb(1) > 0 ) THEN + do_ocnskeb=.true. + ENDIF ! compute frequencty to update random pattern IF (epblint == 0.) epblint=deltim nsepbl=nint(epblint/deltim) ! epblint in seconds @@ -421,7 +513,14 @@ subroutine compns_stochy_ocn (deltim,iret) IF (ocnspptint == 0.) ocnspptint=deltim nsocnsppt=nint(ocnspptint/deltim) ! ocnspptint in seconds IF(nsocnsppt<=0 .or. abs(nsocnsppt-ocnspptint/deltim)>tol) THEN - WRITE(0,*) "ePBL interval is invalid",ocnspptint + WRITE(0,*) "ocnsppt interval is invalid",ocnspptint + iret=9 + return + ENDIF + IF (ocnskebint == 0.) ocnskebint=deltim + nsocnskeb=nint(ocnskebint/deltim) ! ocnskebint in seconds + IF(nsocnskeb<=0 .or. abs(nsocnskeb-ocnskebint/deltim)>tol) THEN + WRITE(0,*) "ocnskeb interval is invalid",ocnskebint iret=9 return ENDIF @@ -433,6 +532,7 @@ subroutine compns_stochy_ocn (deltim,iret) do k=1,5 if (epbl(k).GT.0) l_min=min(epbl_lscale(k),l_min) if (ocnsppt(k).GT.0) l_min=min(ocnsppt_lscale(k),l_min) + if (ocnskeb(k).GT.0) l_min=min(ocnskeb_lscale(k),l_min) enddo !ntrunc=1.5*circ/l_min ntrunc=circ/l_min @@ -459,6 +559,7 @@ subroutine compns_stochy_ocn (deltim,iret) print *, 'ocean stochastic physics' print *, ' pert_epbl : ', pert_epbl print *, ' do_ocnsppt : ', do_ocnsppt + print *, ' do_ocnskeb : ', do_ocnskeb endif iret = 0 if (iseed_epbl(1) > 0) iseed_epbl2(1)=iseed_epbl(1)-1234567 diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 28a0841..edabaab 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, spp_stddev_cutoff_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_transforms,only:colrad_a,latg,lonf,skeblevs use mpi_wrapper, only : mpi_wrapper_initialize,mype,npes,is_rootpe @@ -44,13 +48,17 @@ 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(:), intent(out) :: lndp_var_list_out +real(kind=kind_dbl_prec), dimension(:), intent(out) :: lndp_prt_list_out +character(len=3), dimension(:), intent(out) :: spp_var_list_out +real(kind=kind_dbl_prec), dimension(:), intent(out) :: spp_prt_list_out +real(kind=kind_dbl_prec), dimension(:), intent(out) :: spp_stddev_cutoff_out ! Local variables @@ -58,7 +66,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,i character*2::proc ! Initialize MPI and OpenMP @@ -86,7 +94,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 @@ -107,22 +114,38 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in else if (lndp_type_in /= lndp_type) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings lndp_type in physics and nam_sfcperts' - print*,'lndp_type',lndp_type_in,lndp_type iret = 20 return else if (n_var_lndp_in /= n_var_lndp) then - print*,'n_var_lndp',n_var_lndp_in , n_var_lndp write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' 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_sppperts' + 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 +if (n_var_lndp>0) then + lndp_var_list_out=lndp_var_list(1:n_var_lndp) + lndp_prt_list_out=lndp_prt_list(1:n_var_lndp) +endif +if (n_var_spp>0) then + spp_var_list_out=spp_var_list(1:n_var_spp) + spp_prt_list_out=spp_prt_list(1:n_var_spp) + spp_stddev_cutoff_out=spp_stddev_cutoff(1:n_var_spp) +endif +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 @@ -203,6 +226,19 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in if (is_rootpe()) 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_rootpe()) print *,'spp vert profile',k,sl(k),vfact_spp(k) + enddo +endif ! get interpolation weights ! define gaussian grid lats and lons latghf=latg/2 @@ -223,11 +259,10 @@ end subroutine init_stochastic_physics !!!!!!!!!!!!!!!!!!!! subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_in,do_sppt_in, & - mpiroot, mpicomm, iret) + do_skeb_in, mpiroot, mpicomm, iret) use stochy_data_mod, only : init_stochdata_ocn,gg_lats,gg_lons,& rad2deg,INTTYP,wlon,rnlat,gis_stochy_ocn use spectral_transforms , only : latg,lonf,colrad_a -!use MOM_grid, only : ocean_grid_type use stochy_namelist_def use mersenne_twister, only: random_gauss use mpi_wrapper, only : mpi_wrapper_initialize,mype,npes,is_rootpe @@ -236,7 +271,7 @@ subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_i real,intent(in) :: delt integer,intent(in) :: nx,ny,nz real,intent(in) :: geoLonT(nx,ny),geoLatT(nx,ny) -logical,intent(in) :: pert_epbl_in,do_sppt_in +logical,intent(in) :: pert_epbl_in,do_sppt_in,do_skeb_in integer,intent(in) :: mpiroot, mpicomm integer, intent(out) :: iret real(kind=kind_dbl_prec), parameter :: con_pi =4.0d0*atan(1.0d0) @@ -269,6 +304,11 @@ subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_i & ' namelist settings pert_epbl and epbl' iret = 20 return +else if (do_skeb_in.neqv.do_ocnskeb) then + write(0,'(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & + & ' namelist settings do_skeb and skeb' + iret = 20 + return end if ! get interpolation weights @@ -286,6 +326,7 @@ subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_i enddo WLON=gg_lons(1)-(gg_lons(2)-gg_lons(1)) RNLAT=gg_lats(1)*2-gg_lats(2) +print*,'finished ocean init' end subroutine init_stochastic_physics_ocn !!!!!!!!!!!!!!!!!!!! @@ -297,16 +338,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 + 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 + lndp_type, n_var_lndp, n_var_spp, do_spp, spp_stddev_cutoff, spp_prt_list use mpi_wrapper, only: is_rootpe implicit none @@ -319,17 +361,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(:,:,:) +real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:),tmpl_wts(:,:,:),tmp_spp_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) @@ -402,6 +446,23 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s ENDDO ENDDO deallocate(tmpl_wts) +endif +if (n_var_spp .GE. 1) then + 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_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 + deallocate(tmp_spp_wts) endif deallocate(tmp_wts) deallocate(tmpu_wts) @@ -410,45 +471,50 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s end subroutine run_stochastic_physics -subroutine run_stochastic_physics_ocn(sppt_wts,t_rp1,t_rp2) -!use MOM_forcing_type, only : mech_forcing -!use MOM_grid, only : ocean_grid_type +subroutine run_stochastic_physics_ocn(sppt_wts,skeb_wts,t_rp1,t_rp2) use stochy_internal_state_mod -use stochy_data_mod, only : nepbl,nocnsppt,rpattern_epbl1,rpattern_epbl2,rpattern_ocnsppt, gis_stochy_ocn +use stochy_data_mod, only : nepbl,nocnsppt,nocnskeb,rpattern_epbl1,rpattern_epbl2,rpattern_ocnsppt, rpattern_ocnskeb, gis_stochy_ocn use get_stochy_pattern_mod,only : get_random_pattern_scalar use stochy_namelist_def implicit none -!type(ocean_grid_type), intent(in) :: G -real, intent(inout) :: sppt_wts(:,:),t_rp1(:,:),t_rp2(:,:) +real, intent(inout) :: sppt_wts(:,:),t_rp1(:,:),t_rp2(:,:),skeb_wts(:,:) real, allocatable :: tmp_wts(:,:) -if (pert_epbl .OR. do_ocnsppt) then +if (pert_epbl .OR. do_ocnsppt .OR. do_ocnskeb ) then allocate(tmp_wts(gis_stochy_ocn%nx,gis_stochy_ocn%ny)) if (pert_epbl) then call get_random_pattern_scalar(rpattern_epbl1,nepbl,gis_stochy_ocn,tmp_wts) t_rp1(:,:)=2.0/(1+exp(-1*tmp_wts)) call get_random_pattern_scalar(rpattern_epbl2,nepbl,gis_stochy_ocn,tmp_wts) t_rp2(:,:)=2.0/(1+exp(-1*tmp_wts)) - else - t_rp1(:,:)=1.0 - t_rp2(:,:)=1.0 +! else +! t_rp1(:,:)=1.0 +! t_rp2(:,:)=1.0 endif if (do_ocnsppt) then call get_random_pattern_scalar(rpattern_ocnsppt,nocnsppt,gis_stochy_ocn,tmp_wts) sppt_wts=2.0/(1+exp(-1*tmp_wts)) - else - sppt_wts=1.0 +! else +! sppt_wts=1.0 endif + if (do_ocnskeb) then + call get_random_pattern_scalar(rpattern_ocnskeb,nocnskeb,gis_stochy_ocn,skeb_wts,normalize=.true.) +! else +! skeb_wts=1.0 + endif + !print*,'after get_random_pattern_scalar',skeb_wts(1,1) deallocate(tmp_wts) -else - sppt_wts(:,:)=1.0 - t_rp1(:,:)=1.0 - t_rp2(:,:)=1.0 +!else +! sppt_wts(:,:)=1.0 +! skeb_wts(:,:)=1.0 +! t_rp1(:,:)=1.0 +! t_rp2(:,:)=1.0 endif 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_transforms, only : lat1s_a ,lon_dims_a,wgt_a,sinlat_a,coslat_a,colrad_a,rcs2_a implicit none @@ -475,6 +541,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 f844856..3d10b39 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -10,7 +10,8 @@ module stochy_data_mod use constants_mod, only : radius use mpi_wrapper, only: mp_bcst, is_rootpe, mype use stochy_patterngenerator_mod, only: random_pattern, patterngenerator_init,& - getnoise, patterngenerator_advance,ndimspec,chgres_pattern,computevarspec_r + getnoise, getnoise_un, patterngenerator_advance, patterngenerator_advance_jb, ndimspec, & + chgres_pattern, computevarspec_r use stochy_internal_state_mod ! use mersenne_twister_stochy, only : random_seed use mersenne_twister, only : random_seed @@ -21,16 +22,19 @@ 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,rpattern_ocnskeb integer, public :: nepbl=0 integer, public :: nocnsppt=0 + integer, public :: nocnskeb=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(:,:) integer ,public, allocatable :: skeb_vpts(:,:) real(kind=kind_dbl_prec),public, allocatable :: gg_lats(:),gg_lons(:) @@ -56,7 +60,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) integer, intent(out) :: iret real :: ones(5) - real :: rnn1 + real :: rnn1,gamma_sum integer :: nn,k,nm,stochlun,ierr,n integer :: locl,indev,indod,indlsod,indlsev integer :: l,jbasev,jbasod @@ -74,9 +78,9 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) call compns_stochy (mype,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 ( (.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 + call initialize_spectral(gis_stochy) - allocate(noise_e(len_trie_ls,2),noise_o(len_trio_ls,2)) ! determine number of random patterns to be used for each scheme. do n=1,size(sppt) @@ -113,13 +117,16 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) ! endif !enddo if (n_var_lndp>0) nlndp=1 + if (n_var_spp>0) nspp=n_var_spp if (is_rootpe()) print *,' nlndp = ', nlndp + if (is_rootpe()) 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_rootpe()) then @@ -408,6 +415,60 @@ 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_rootpe()) 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_rootpe()) print *, 'Initialize random pattern for SPP PERTS' + if (stochini) then + 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,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,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,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,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,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, 1, minval(rpattern_spp(n)%spec_o(:,:,1)), maxval(rpattern_spp(n)%spec_o(:,:,1)) + endif ! stochini + enddo ! n, nspp + endif ! nspp > 0 + if (is_rootpe() .and. stochini) CLOSE(stochlun) deallocate(noise_e,noise_o) end subroutine init_stochdata @@ -423,9 +484,10 @@ subroutine init_stochdata_ocn(nlevs,delt,iret) real, intent(in) :: delt integer, intent(out) :: iret - integer :: nn,nm,stochlun,n,jcapin + integer :: nn,nm,stochlun,n,jcapin,n2 integer :: l,jbasev,jbasod - integer :: indev,indod,indlsod,indlsev,varid1,varid2,varid3,varid4,ierr + integer :: locl,indev,indod,indlsod,indlsev,varid1,varid2,varid3,varid4,ierr + real :: gamma_sum,pi real(kind_dbl_prec),allocatable :: noise_e(:,:),noise_o(:,:) include 'function_indlsod' @@ -434,10 +496,11 @@ subroutine init_stochdata_ocn(nlevs,delt,iret) stochlun=99 levs=nlevs + pi = 4.0d0*atan(1.0d0) iret=0 call compns_stochy_ocn (delt,iret) if(is_rootpe()) print*,'in init stochdata_ocn' - if ( (.NOT. pert_epbl) .AND. (.NOT. do_ocnsppt) ) return + if ( (.NOT. pert_epbl) .AND. (.NOT. do_ocnsppt) .AND. (.NOT. do_ocnskeb) ) return call initialize_spectral(gis_stochy_ocn) if (iret/=0) return allocate(noise_e(len_trie_ls,2),noise_o(len_trio_ls,2)) @@ -458,12 +521,21 @@ subroutine init_stochdata_ocn(nlevs,delt,iret) endif enddo + do n=1,size(ocnskeb) + if (ocnskeb(n) > 0) then + nocnskeb=nocnskeb+1 + else + exit + endif + enddo + if (nepbl > 0) then allocate(rpattern_epbl1(nepbl)) allocate(rpattern_epbl2(nepbl)) endif if (nocnsppt > 0) allocate(rpattern_ocnsppt(nocnsppt)) + if (nocnskeb > 0) allocate(rpattern_ocnskeb(nocnskeb)) ! if stochini is true, then read in pattern from a file if (is_rootpe()) then @@ -629,6 +701,67 @@ subroutine init_stochdata_ocn(nlevs,delt,iret) endif enddo endif + + if (nocnskeb > 0) then + if (is_rootpe()) then + if (stochini) then + ierr=NF90_INQ_VARID(stochlun,"ocnskeb_seed", varid1) + if (ierr .NE. 0) then + write(0,*) 'error inquring OCNSPPT seed' + iret = ierr + return + end if + ierr=NF90_INQ_VARID(stochlun,"ocnskeb_spec", varid2) + if (ierr .NE. 0) then + write(0,*) 'error inquring OCNSPPT spec' + iret = ierr + return + end if + endif + endif + if (is_rootpe()) print *, 'Initialize random pattern for ocnskeb' + call patterngenerator_init(ocnskeb_lscale(1:nocnskeb),ocnskebint,ocnskeb_tau(1:nocnskeb),ocnskeb(1:nocnskeb),iseed_ocnskeb,rpattern_ocnskeb, & + lonf,latg,jcap,gis_stochy_ocn%ls_node,nocnskeb,1,0,new_lscale,.true.) + + do n=1,nocnskeb + if (stochini) then + call read_pattern(rpattern_ocnskeb(n),jcapin,stochlun,1,n,varid1,varid2,.false.,ierr) + if (ierr .NE. 0) then + write(0,*) 'error reading OCNSKEB pattern' + iret = ierr + return + end if + else + call getnoise_un(rpattern_ocnskeb(n),noise_e,noise_o) + !call getnoise(rpattern_ocnskeb(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_ocnskeb(n)%spec_e(nn,1,1)=noise_e(nn,1) + rpattern_ocnskeb(n)%spec_e(nn,2,1)=noise_e(nn,2) + nm = rpattern_ocnskeb(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_ocnskeb(n)%spec_e(nn,1,1) = rpattern_ocnskeb(n)%stdev*rpattern_ocnskeb(n)%spec_e(nn,1,1)*rpattern_ocnskeb(n)%varspectrum(nm) + rpattern_ocnskeb(n)%spec_e(nn,2,1) = rpattern_ocnskeb(n)%stdev*rpattern_ocnskeb(n)%spec_e(nn,2,1)*rpattern_ocnskeb(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_ocnskeb(n)%spec_o(nn,1,1)=noise_o(nn,1) + rpattern_ocnskeb(n)%spec_o(nn,2,1)=noise_o(nn,2) + nm = rpattern_ocnskeb(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_ocnskeb(n)%spec_o(nn,1,1) = rpattern_ocnskeb(n)%stdev*rpattern_ocnskeb(n)%spec_o(nn,1,1)*rpattern_ocnskeb(n)%varspectrum(nm) + rpattern_ocnskeb(n)%spec_o(nn,2,1) = rpattern_ocnskeb(n)%stdev*rpattern_ocnskeb(n)%spec_o(nn,2,1)*rpattern_ocnskeb(n)%varspectrum(nm) + enddo + print*,'calling patterngenerator_advance norm init' + call patterngenerator_advance_jb(rpattern_ocnskeb(n)) + !call patterngenerator_advance(rpattern_ocnskeb(n)) +! if (is_rootpe()) then +! print*,'after advance' +! do nn=1,81,5 +! print*,nn,rpattern_ocnskeb(n)%spec_o(nn,1,1),noise_o(nn,1) +! enddo +! endif + endif + enddo + endif deallocate(noise_e,noise_o) end subroutine init_stochdata_ocn diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index b9d0ca2..922e2a4 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -10,8 +10,9 @@ module stochy_namelist_def implicit none public - integer, parameter :: max_n_var_lndp = 6 ! must match value used in GFS_typedefs - integer nssppt,nsshum,nsepbl,nsocnsppt,nsskeb,lon_s,lat_s,ntrunc + integer, parameter :: max_n_var_lndp = 20 ! 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,nsocnskeb, nsskeb,lon_s,lat_s,ntrunc ! pjp stochastic phyics integer skeb_varspect_opt,skeb_npass @@ -20,23 +21,33 @@ module stochy_namelist_def real(kind=kind_dbl_prec) :: skeb_sigtop1,skeb_sigtop2, & sppt_sigtop1,sppt_sigtop2,shum_sigefold, & skeb_vdof - real(kind=kind_dbl_prec) skeb_diss_smooth,epblint,ocnspptint,spptint,shumint,skebint,skebnorm + real(kind=kind_dbl_prec) skeb_diss_smooth,epblint,ocnspptint,ocnskebint,spptint,shumint,skebint,skebnorm real(kind=kind_dbl_prec), dimension(5) :: skeb,skeb_lscale,skeb_tau real(kind=kind_dbl_prec), dimension(5) :: sppt,sppt_lscale,sppt_tau real(kind=kind_dbl_prec), dimension(5) :: shum,shum_lscale,shum_tau real(kind=kind_dbl_prec), dimension(5) :: epbl,epbl_lscale,epbl_tau real(kind=kind_dbl_prec), dimension(5) :: ocnsppt,ocnsppt_lscale,ocnsppt_tau + real(kind=kind_dbl_prec), dimension(5) :: ocnskeb,ocnskeb_lscale,ocnskeb_tau integer,dimension(5) ::skeb_vfilt - integer(kind=kind_dbl_prec),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_epbl,iseed_ocnsppt,iseed_epbl2 + integer(kind=kind_dbl_prec),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_epbl,iseed_ocnsppt,iseed_ocnskeb,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,do_ocnskeb real(kind=kind_dbl_prec), dimension(5) :: lndp_lscale,lndp_tau integer n_var_lndp integer(kind=kind_dbl_prec),dimension(5) ::iseed_lndp integer lndp_type + integer lndp_model_type 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 diff --git a/stochy_patterngenerator.F90 b/stochy_patterngenerator.F90 index 7060765..0f0bb3b 100644 --- a/stochy_patterngenerator.F90 +++ b/stochy_patterngenerator.F90 @@ -1,5 +1,7 @@ -!>@brief The module 'stochy_patterngenerator_mod' contains the derived type random_pattern +!>nbrief The module 'stochy_patterngenerator_mod' contains the derived type random_pattern !! which controls the characteristics of the random pattern +!! This is a modified version of the original one stochy_patterngenerator.F90, where the +!! the random patterns are not properly normalized module stochy_patterngenerator_mod !> generate random patterns with specified temporal and spatial auto-correlation @@ -14,9 +16,9 @@ module stochy_patterngenerator_mod implicit none private - public :: computevarspec, setvarspect,& - patterngenerator_init, patterngenerator_destroy, getnoise, & - patterngenerator_advance, ndimspec,& + public :: computevarspec, setvarspect, & + patterngenerator_init, patterngenerator_destroy, getnoise, getnoise_un, & + patterngenerator_advance, patterngenerator_advance_jb, ndimspec, & chgres_pattern,computevarspec_r ! ----------------------------------------------- @@ -24,12 +26,12 @@ module stochy_patterngenerator_mod !>@details A seperate instance of this type is needed for each pattern type,public :: random_pattern ! start type define ! ----------------------------------------------- - real(kind_dbl_prec), public :: lengthscale ! length scale in m - real(kind_dbl_prec), public :: tau - real(kind_dbl_prec), public :: dt - real(kind_dbl_prec), public :: phi - real(kind_dbl_prec), public :: stdev - real(kind_evod), allocatable, dimension(:), public :: varspectrum, varspectrum1d, lap + real(kind_phys), public :: lengthscale ! length scale in m + real(kind_phys), public :: tau + real(kind_phys), public :: dt + real(kind_phys), public :: phi + real(kind_phys), public :: stdev + real(kind_dbl_prec), allocatable, dimension(:), public :: varspectrum, varspectrum1d, lap integer, allocatable, dimension(:), public ::& degree,order,idx_e,idx_o integer, allocatable, dimension(:,:), public :: idx @@ -49,23 +51,32 @@ module stochy_patterngenerator_mod !! temporaral and spatial correlations. subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& nlon, nlat, jcap, ls_nodes, npatterns,& - nlevs, varspect_opt,new_lscale) + nlevs, varspect_opt,new_lscale,bn) !\callgraph - real(kind_dbl_prec), intent(in),dimension(npatterns) :: lscale,tscale,stdev + real(kind_phys), intent(in),dimension(npatterns) :: lscale,tscale,stdev real, intent(in) :: delt integer, intent(in) :: nlon,nlat,jcap,npatterns,varspect_opt integer, intent(in) :: ls_nodes(ls_dim,3),nlevs logical, intent(in) :: new_lscale + logical, optional, intent(in) :: bn ! berner normalization type(random_pattern), intent(out), dimension(npatterns) :: rpattern integer(8), intent(inout) :: iseed(npatterns) integer m,j,l,n,nm,nn,np,indev1,indev2,indod1,indod2 integer(8) count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 integer count4, ierr + logical :: bn_local ! berner normalization ! integer member_id integer indlsod,indlsev,jbasev,jbasod include 'function_indlsod' include 'function_indlsev' + bn_local=.false. + if (present(bn)) bn_local=bn + if (present(bn)) then + print*,'Berner norm=',bn + else + print*,'old norm' + endif nlons = nlon nlats = nlat ntrunc = jcap @@ -139,6 +150,7 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& rpattern(np)%lengthscale = lscale(np) rpattern(np)%dt = delt rpattern(np)%phi = exp(-delt/tscale(np)) + !rpattern(np)%phi = 0 rpattern(np)%stdev = stdev(np) allocate(rpattern(np)%varspectrum(ndimspec)) allocate(rpattern(np)%varspectrum1d(0:ntrunc)) @@ -168,13 +180,14 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& rpattern(np)%seed = count4 ! set seed (to be the same) on all tasks. Save random state. call random_setseed(rpattern(np)%seed,rpattern(np)%rstate) - if (varspect_opt .ne. 0 .and. varspect_opt .ne. 1) then + !print*,'calling setvarspect',bn_local + if (varspect_opt .lt. 0 .or. varspect_opt .gt. 1) then if (is_rootpe()) then - print *,'WARNING: illegal value for varspect_opt (should be 0 or 1), using 0 (gaussian spectrum)...' + print *,'WARNING: illegal value for varspect_opt (should be 0,1,or 2), using 0 (gaussian spectrum)...' endif - call setvarspect(rpattern(np),0,new_lscale) + call setvarspect(rpattern(np),0,new_lscale,bn_local) else - call setvarspect(rpattern(np),varspect_opt,new_lscale) + call setvarspect(rpattern(np),varspect_opt,new_lscale,bn_local) endif enddo ! n=1,npatterns end subroutine patterngenerator_init @@ -200,8 +213,8 @@ end subroutine patterngenerator_destroy subroutine computevarspec(rpattern,dataspec,var) !\callgraph ! compute globally integrated variance from spectral coefficients - complex(kind_evod), intent(in) :: dataspec(ndimspec) - real(kind_evod), intent(out) :: var + complex(kind_dbl_prec), intent(in) :: dataspec(ndimspec) + real(kind_dbl_prec), intent(out) :: var type(random_pattern), intent(in) :: rpattern integer n var = 0. @@ -220,8 +233,8 @@ end subroutine computevarspec subroutine computevarspec_r(rpattern,dataspec,var) !\callgraph ! compute globally integrated variance from spectral coefficients - real(kind_dbl_prec), intent(in) :: dataspec(2*ndimspec) - real(kind_dbl_prec), intent(out) :: var + real(kind_phys), intent(in) :: dataspec(2*ndimspec) + real(kind_phys), intent(out) :: var type(random_pattern), intent(in) :: rpattern integer n var = 0. @@ -239,11 +252,11 @@ end subroutine computevarspec_r !! variance from real spectral c subroutine getnoise(rpattern,noise_e,noise_o) !\callgraph - real(kind_dbl_prec), intent(out) :: noise_e(len_trie_ls,2) - real(kind_dbl_prec), intent(out) :: noise_o(len_trio_ls,2) + real(kind_phys), intent(out) :: noise_e(len_trie_ls,2) + real(kind_phys), intent(out) :: noise_o(len_trio_ls,2) ! generate white noise with unit variance in spectral space type(random_pattern), intent(inout) :: rpattern - real :: noise(2*ndimspec) + real(kind_dbl_prec) :: noise(2*ndimspec) integer nm,nn call random_gauss(noise,rpattern%rstate) noise(1) = 0.; noise(ndimspec+1) = 0. @@ -272,14 +285,49 @@ subroutine getnoise(rpattern,noise_e,noise_o) enddo end subroutine getnoise +!>@brief The subroutine 'getnoise' scales spectral cofficients with +!! white noise to the appropriate amplitude for speherical harmonincs +!! variance from real spectral c + subroutine getnoise_un(rpattern,noise_e,noise_o) +!\callgraph + real(kind_phys), intent(out) :: noise_e(len_trie_ls,2) + real(kind_phys), intent(out) :: noise_o(len_trio_ls,2) + ! generate white noise with unit variance in spectral space + type(random_pattern), intent(inout) :: rpattern + real(kind_dbl_prec) :: noise(2*ndimspec) + integer nm,nn + call random_gauss(noise,rpattern%rstate) + noise(1) = 0.; noise(ndimspec+1) = 0. + noise_e = 0.; noise_o = 0. + ! subset + do nn=1,len_trie_ls + nm = rpattern%idx_e(nn) + if (nm == 0) cycle + noise_e(nn,1) = noise(nm) + noise_e(nn,2) = noise(ndimspec+nm) + if (rpattern%order(nm) .eq. 0) then + noise_e(nn,2) = 0. + endif + enddo + do nn=1,len_trio_ls + nm = rpattern%idx_o(nn) + if (nm == 0) cycle + noise_o(nn,1) = noise(nm) + noise_o(nn,2) = noise(ndimspec+nm) + if (rpattern%order(nm) .eq. 0) then + noise_o(nn,2) = 0. + endif + enddo + end subroutine getnoise_un + !>@brief The subroutine 'patterngenerator_advance' advance 1st-order autoregressive process subroutine patterngenerator_advance(rpattern,k,skeb_first_call) !\callgraph ! advance 1st-order autoregressive process with ! specified autocorrelation (phi) and variance spectrum (spectrum) - real(kind_dbl_prec) :: noise_e(len_trie_ls,2) - real(kind_dbl_prec) :: noise_o(len_trio_ls,2) + real(kind_phys) :: noise_e(len_trie_ls,2) + real(kind_phys) :: noise_o(len_trio_ls,2) type(random_pattern), intent(inout) :: rpattern logical, intent(in) :: skeb_first_call integer j,l,n,nn,nm,k,k2 @@ -307,65 +355,158 @@ subroutine patterngenerator_advance(rpattern,k,skeb_first_call) enddo end subroutine patterngenerator_advance +!>@brief The subroutine 'patterngenerator_advance_jb' advance 1st-order autoregressive process + subroutine patterngenerator_advance_jb(rpattern) +!\callgraph + + ! advance 1st-order autoregressive process with + ! specified autocorrelation (phi) and variance spectrum (spectrum) + real(kind_phys) :: noise_e(len_trie_ls,2) + real(kind_phys) :: noise_o(len_trio_ls,2) + type(random_pattern), intent(inout) :: rpattern + integer j,l,n,nn,nm + call getnoise_un(rpattern,noise_e,noise_o) + !call getnoise(rpattern,noise_e,noise_o) + !print*,'in advance_jb' + !print*, 'random pattern stdev:', rpattern%stdev + do nn=1,len_trie_ls + nm = rpattern%idx_e(nn) + if (nm == 0) cycle + rpattern%spec_e(nn,1,1) = rpattern%phi*rpattern%spec_e(nn,1,1) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*noise_e(nn,1)*rpattern%varspectrum(nm) + rpattern%spec_e(nn,2,1) = rpattern%phi*rpattern%spec_e(nn,2,1) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*noise_e(nn,2)*rpattern%varspectrum(nm) + enddo + do nn=1,len_trio_ls + nm = rpattern%idx_o(nn) + if (nm == 0) cycle + rpattern%spec_o(nn,1,1) = rpattern%phi*rpattern%spec_o(nn,1,1) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*noise_o(nn,1)*rpattern%varspectrum(nm) + rpattern%spec_o(nn,2,1) = rpattern%phi*rpattern%spec_o(nn,2,1) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*noise_o(nn,2)*rpattern%varspectrum(nm) + enddo + end subroutine patterngenerator_advance_jb + !>@brief The subroutine 'setvarspect' calculates the variance spectrum ! from a specified decorrelation lengthscale - subroutine setvarspect(rpattern,varspect_opt,new_lscale) + subroutine setvarspect(rpattern,varspect_opt,new_lscale,berner_normalize) !\callgraph ! define variance spectrum (isotropic covariance) ! normalized to unit global variance type(random_pattern), intent(inout) :: rpattern integer, intent(in) :: varspect_opt logical, intent(in) :: new_lscale - integer :: n - complex(kind_evod) noise(ndimspec) - real(kind_evod) var,rerth,inv_rerth_sq + logical, intent(in) :: berner_normalize + integer :: n, nm + complex(kind_dbl_prec) noise(ndimspec) + real(kind_dbl_prec) var,b_jb,rerth,inv_rerth_sq,pi,gamma_sum,deltaE + pi = 4.0d0*atan(1.0d0) rerth =6.3712e+6 ! radius of earth (m) inv_rerth_sq=1.0/(rerth**2) + !print*,'setvarspect,ntrunc',berner_normalize,ntrunc ! 1d variance spectrum (as a function of total wavenumber) if (varspect_opt == 0) then ! gaussian + print*, 'Gaussian variance spectrum' ! rpattern%lengthscale is interpreted as an efolding length ! scale, in meters. - do n=0,ntrunc - rpattern%varspectrum1d(n) = exp(-rpattern%lengthscale**2*(float(n)*(float(n)+1.))/(4.*rerth**2)) - enddo ! scaling factors for spectral coeffs of white noise pattern with unit variance if (new_lscale) then - !fix for proper lengthscale - rpattern%varspectrum = ntrunc*exp((rpattern%lengthscale*0.25)**2*rpattern%lap*inv_rerth_sq) + !fix for proper lengthscale + print*, 'Proper lengthscale condition' +! print*, 'rpattern lap:',rpattern%lap + rpattern%varspectrum = exp((rpattern%lengthscale*0.25)**2*rpattern%lap*inv_rerth_sq) + do n=0,ntrunc + rpattern%varspectrum1d(n) = exp(-(rpattern%lengthscale*0.25)**2*float(n)*(float(n)+1.)*inv_rerth_sq) + enddo else - rpattern%varspectrum = sqrt(ntrunc*exp(rpattern%lengthscale**2*rpattern%lap/(4.*rerth**2))) + print*, 'Not proper lengthscale condition' + !print*, 'rpattern lap:',rpattern%lap + rpattern%varspectrum = exp(rpattern%lengthscale**2*rpattern%lap/(4.*rerth**2)) + do n=0,ntrunc + rpattern%varspectrum1d(n) = exp(-rpattern%lengthscale**2*(float(n)*(float(n)+1.))/(4.*rerth**2)) + enddo + print*,'Finished' endif else if (varspect_opt == 1) then ! power law ! rpattern%lengthscale is interpreted as a power, not a length. - do n=0,ntrunc + !print*,'Power Law lengthscale',rpattern%lengthscale + rpattern%varspectrum1d(0) = 0.0 + do n=1,ntrunc rpattern%varspectrum1d(n) = float(n)**(rpattern%lengthscale) enddo ! scaling factors for spectral coeffs of white noise pattern with unit variance - rpattern%varspectrum = sqrt(ntrunc*(rpattern%degree**(rpattern%lengthscale))) + rpattern%varspectrum(1) = 0 + rpattern%varspectrum(2:) = rpattern%degree(2:)**(rpattern%lengthscale) + + !do n=0,ntrunc + ! rpattern%varspectrum1d(n) = float(n)**(rpattern%lengthscale) + !enddo + ! scaling factors for spectral coeffs of white noise pattern with unit variance + !rpattern%varspectrum = sqrt(ntrunc*(rpattern%degree**(rpattern%lengthscale))) + !rpattern%varspectrum = rpattern%degree**(rpattern%lengthscale)! modified following equation 4 in Berner et al 2009 endif - noise = 0. - do n=1,ndimspec - if (rpattern%order(n) .ne. 0.) then - noise(n) = cmplx(1.,1.)/sqrt(2.*rpattern%degree(n)+1) - else - noise(n) = sqrt(2.)/sqrt(2.*rpattern%degree(n)+1.) - endif - enddo - noise(1) = 0 ! no global mean. - ! make sure global mean variance is 1. - noise = noise*sqrt(1./ntrunc) - noise = rpattern%varspectrum*noise - call computevarspec(rpattern,noise,var) - rpattern%varspectrum = rpattern%varspectrum/sqrt(var) - rpattern%varspectrum1d = rpattern%varspectrum1d/var + print*, 'Berner Normalization:',berner_normalize + if (.not. berner_normalize) then ! normalize the GEFS way + noise = 0. + do n=1,ndimspec + if (rpattern%order(n) .ne. 0.) then + noise(n) = cmplx(1.,1.)/sqrt(2.*rpattern%degree(n)+1) + else + noise(n) = sqrt(2.)/sqrt(2.*rpattern%degree(n)+1.) + endif + enddo + noise(1) = 0 ! no global mean. + ! make sure global mean variance is 1. + noise = noise*sqrt(1./ntrunc) + noise = rpattern%varspectrum*noise + call computevarspec(rpattern,noise,var) + rpattern%varspectrum = rpattern%varspectrum/sqrt(var) + rpattern%varspectrum1d = rpattern%varspectrum1d/var + + else if (berner_normalize) then ! normalize by Berner et al. 2009 + do n=1,len_trie_ls + nm = rpattern%idx_e(n) + if (nm == 0) cycle + if (rpattern%order(nm) .ne. 0.) then + rpattern%varspectrum(nm) = rpattern%varspectrum(nm)/sqrt(ntrunc*(2.*rpattern%degree(nm)+1)*(1+rpattern%phi)) + else + rpattern%varspectrum(nm) = rpattern%varspectrum(nm)*sqrt(2.)/sqrt(ntrunc*(2.*rpattern%degree(n)+1.)*(1+rpattern%phi)) + endif + enddo + do n=1,len_trio_ls + nm = rpattern%idx_o(n) + if (nm == 0) cycle + if (rpattern%order(nm) .ne. 0.) then + rpattern%varspectrum(nm) = rpattern%varspectrum(nm)/sqrt(ntrunc*(2.*rpattern%degree(nm)+1)*(1+rpattern%phi)) + else + rpattern%varspectrum(nm) = rpattern%varspectrum(nm)*sqrt(2.)/sqrt(ntrunc*(2.*rpattern%degree(n)+1.)*(1+rpattern%phi)) + endif + enddo + do n=0,ntrunc + rpattern%varspectrum1d(n) = rpattern%varspectrum1d(n)/sqrt(ntrunc*(2.*n+1.)) + enddo + gamma_sum=0 + !print*,'summing gamma' + do n=1,ntrunc + !print*,n,gamma_sum,rpattern%varspectrum1d(n)**2 + !gamma_sum = gamma_sum+n*(n+1.)*(2*n+1)*(rpattern%varspectrum1d(n)**2) + gamma_sum = gamma_sum+n*(n+1.)*(2.*n+2.)*(rpattern%varspectrum1d(n)**2) + enddo + !b_jb=sqrt(((1-rpattern%phi)*4*pi*(rerth**2))/(gamma_sum*sqrt(var))) !!!Changed phi to (1-phi) + b_jb = sqrt(4*pi*(rerth**2)*(1-rpattern%phi)/gamma_sum) + !print*,'normalizing by b_jb',gamma_sum,b_jb + rpattern%varspectrum = rpattern%varspectrum*b_jb + rpattern%varspectrum1d = rpattern%varspectrum1d*b_jb + endif end subroutine setvarspect + !>@brief The subroutine 'chgres_pattern' truncates the spherical harmonics if !! restarting from a higher-resolution pattern subroutine chgres_pattern(pattern2din,pattern2dout,ntruncin,ntruncout) !\callgraph - real(kind_dbl_prec), intent(in) :: pattern2din((ntruncin+1)*(ntruncin+2)) - real(kind_dbl_prec), intent(out) :: pattern2dout((ntruncout+1)*(ntruncout+2)) + real(kind_phys), intent(in) :: pattern2din((ntruncin+1)*(ntruncin+2)) + real(kind_phys), intent(out) :: pattern2dout((ntruncout+1)*(ntruncout+2)) integer, intent(in) :: ntruncin,ntruncout integer :: m,n,nm,ndimsspecin,ndimsspecout integer,allocatable, dimension(:,:):: idxin From 52227ea4ffd340c77988f7f9714fc10e54873fbb Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Sat, 25 May 2024 11:04:09 -0600 Subject: [PATCH 2/8] Comment print statements --- get_stochy_pattern.F90 | 184 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 166 insertions(+), 18 deletions(-) diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index f385856..59d58f4 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -5,22 +5,24 @@ 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, & + nocnsppt,nocnskeb,nepbl,nlndp, & rnlat, rpattern_sfc, rpattern_skeb, & rpattern_shum, rpattern_sppt, rpattern_ocnsppt,& rpattern_epbl1, rpattern_epbl2, skebu_save, & + nspp,rpattern_spp, rpattern_ocnskeb, & skebv_save, skeb_vwts, skeb_vpts, wlon use stochy_patterngenerator_mod, only: random_pattern, ndimspec, & - patterngenerator_advance + patterngenerator_advance, & + patterngenerator_advance_jb use stochy_internal_state_mod, only: stochy_internal_state use mpi_wrapper, only : mp_reduce_sum,is_rootpe use mersenne_twister, only: random_seed 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. @@ -129,8 +131,8 @@ subroutine get_random_pattern_vector(rpattern,npatterns,& if (.not. stochini) call patterngenerator_advance(rpattern(n),k,first_call) ! ke norm (convert streamfunction forcing to vorticity forcing) do nn=1,2 - vrtspec_e(:,nn) = gis_stochy%kenorm_e*rpattern(n)%spec_e(:,nn,k) - vrtspec_o(:,nn) = gis_stochy%kenorm_o*rpattern(n)%spec_o(:,nn,k) + vrtspec_e(:,nn) = gis_stochy%kenorm_e(:)*rpattern(n)%spec_e(:,nn,k) + vrtspec_o(:,nn) = gis_stochy%kenorm_o(:)*rpattern(n)%spec_o(:,nn,k) enddo ! convert to winds call vrtdivspect_to_uvgrid( divspec_e,divspec_o,vrtspec_e,vrtspec_o,& @@ -222,13 +224,14 @@ end subroutine get_random_pattern_vector !>@brief The subroutine 'get_random_pattern_scalar' 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_scalar(rpattern,npatterns,& - gis_stochy,pattern_2d) + gis_stochy,pattern_2d,normalize) ! 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 + logical,intent(in), optional :: normalize integer i,j,lat,n real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d @@ -244,9 +247,21 @@ subroutine get_random_pattern_scalar(rpattern,npatterns,& kmsk0 = 0 glolal = 0. do n=1,npatterns - call patterngenerator_advance(rpattern(n),1,.false.) - call scalarspect_to_gaugrid(rpattern(n),gis_stochy, & - wrk2d,1) + if (present(normalize)) then + if (normalize) then + !print*,'calling patterngenerator_advance norm' + call patterngenerator_advance_jb(rpattern(n)) + call scalarspect_to_gaugrid_norm(rpattern(n),gis_stochy, wrk2d,1) + else + !print*,'calling patterngenerator_advance' + call patterngenerator_advance(rpattern(n),1,.false.) + call scalarspect_to_gaugrid(rpattern(n),gis_stochy, wrk2d,1) + endif + else + !print*,'calling patterngenerator_advance' + call patterngenerator_advance(rpattern(n),1,.false.) + call scalarspect_to_gaugrid(rpattern(n),gis_stochy, wrk2d,1) + endif glolal = glolal + wrk2d(:,:,1) enddo @@ -274,7 +289,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 @@ -315,6 +384,54 @@ subroutine scalarspect_to_gaugrid(rpattern,gis_stochy,datag,n) return end subroutine scalarspect_to_gaugrid +!>@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 +subroutine scalarspect_to_gaugrid_norm(rpattern,gis_stochy,datag,n) +!\callgraph + + implicit none + type(random_pattern), intent(in) :: rpattern + type(stochy_internal_state), intent(in) :: gis_stochy + integer , intent(in) :: n + real(kind=kind_dbl_prec), intent(out) :: datag(lonf,gis_stochy%lats_node_a) +! local vars + real(kind=kind_dbl_prec) :: spec_e(len_trie_ls,2),spec_o(len_trio_ls,2) + real(kind=kind_dbl_prec) for_gr_a_1(gis_stochy%lon_dim_a,1,gis_stochy%lats_dim_a) + real(kind=kind_dbl_prec) for_gr_a_2(lonf,1,gis_stochy%lats_dim_a) + integer i,k + integer lan,lat +! normalize the spectral coefficients + do i=1,2 + !spec_e(:,i) = gis_stochy%kenorm_e(:)*rpattern%spec_e(:,i,n) + !spec_o(:,i) = gis_stochy%kenorm_o(:)*rpattern%spec_o(:,i,n) + spec_e(:,i) = rpattern%spec_e(:,i,n) + spec_o(:,i) = rpattern%spec_o(:,i,n) + enddo + !if (is_rootpe())print*,'spec_e=',gis_stochy%kenorm_e(1:20),rpattern%spec_e(1:20,1,n) + call spec_to_four(spec_e(:,:), spec_o(:,:), & + gis_stochy%plnev_a,gis_stochy%plnod_a,& + gis_stochy%ls_node, & + gis_stochy%lats_dim_a,for_gr_a_1,& + gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,& + gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,& + gis_stochy%lats_node_a,gis_stochy%ipt_lats_node_a,1) + do lan=1,gis_stochy%lats_node_a + lat = gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+lan) + call four_to_grid(for_gr_a_1(:,:,lan),for_gr_a_2(:,:,lan),& + gis_stochy%lon_dim_a,1) + enddo + + datag = 0. + do lan=1,gis_stochy%lats_node_a + lat = gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+lan) + do i=1,lonf + datag(i,lan) = for_gr_a_2(i,1,lan) + enddo + enddo + + return + end subroutine scalarspect_to_gaugrid_norm + !>@brief The subroutine 'write_patterns' writes out a single pattern and the seed associated with the random number sequence in netcdf !>@brief The subroutine 'write_stoch_restart_atm' writes out the speherical harmonics to a file, controlled by restart_interval @@ -322,18 +439,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 @@ -347,6 +466,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 @@ -375,6 +500,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' @@ -406,6 +537,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 @@ -421,14 +557,15 @@ end subroutine write_stoch_restart_atm subroutine write_stoch_restart_ocn(sfile) !\callgraph use netcdf - use stochy_namelist_def, only : do_ocnsppt,pert_epbl + use stochy_namelist_def, only : do_ocnsppt,pert_epbl,do_ocnskeb implicit none character(len=*) :: sfile integer :: stochlun,k,n,isize,ierr - integer :: ncid,varid1a,varid1b,varid2a,varid2b,varid3a,varid3b + integer :: ncid,varid1a,varid1b,varid2a,varid2b,varid3a,varid3b,varid4a,varid4b integer :: seed_dim_id,spec_dim_id,np_dim_id include 'netcdf.inc' - if ( ( .NOT. do_ocnsppt) .AND. (.NOT. pert_epbl) ) return + print*,'in write restart',do_ocnsppt,pert_epbl,do_ocnskeb + if ( ( .NOT. do_ocnsppt) .AND. (.NOT. pert_epbl) .AND. ( .NOT. do_ocnskeb) ) return stochlun=99 if (is_rootpe()) then ierr=nf90_create(trim(sfile),cmode=NF90_CLOBBER,ncid=ncid) @@ -456,6 +593,12 @@ subroutine write_stoch_restart_ocn(sfile) ierr=NF90_DEF_VAR(ncid,"epbl2_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid3b) ierr=NF90_PUT_ATT(ncid,varid3b,"long_name","spectral cofficients EPBL2") endif + if (do_ocnskeb) then + ierr=NF90_DEF_VAR(ncid,"ocnskeb_seed",NF90_DOUBLE,(/seed_dim_id, np_dim_id/), varid4a) + ierr=NF90_PUT_ATT(ncid,varid4a,"long_name","random number seed for SPPT") + ierr=NF90_DEF_VAR(ncid,"ocnskeb_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid4b) + ierr=NF90_PUT_ATT(ncid,varid4b,"long_name","spectral cofficients SPPT") + endif ierr=NF90_ENDDEF(ncid) if (ierr .NE. 0) then write(0,*) 'error creating stochastic restart file' @@ -473,6 +616,11 @@ subroutine write_stoch_restart_ocn(sfile) call write_pattern(rpattern_epbl2(n),ncid,1,n,varid3a,varid3b,.false.,ierr) enddo endif + if (nocnskeb > 0) then + do n=1,nocnskeb + call write_pattern(rpattern_ocnskeb(n),ncid,1,n,varid4a,varid4b,.false.,ierr) + enddo + endif if (is_rootpe() ) then ierr=NF90_CLOSE(ncid) if (ierr .NE. 0) then From a7b4b59abc8c4bd8d5d122adff8b74968d6de778 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Sat, 25 May 2024 11:04:43 -0600 Subject: [PATCH 3/8] change filename --- halo_exchange.fv3.F90 => halo_exchange.F90 | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename halo_exchange.fv3.F90 => halo_exchange.F90 (100%) diff --git a/halo_exchange.fv3.F90 b/halo_exchange.F90 similarity index 100% rename from halo_exchange.fv3.F90 rename to halo_exchange.F90 From 8d5e511479a766a55149a9ee44b69e297ff27b70 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Sat, 25 May 2024 11:05:05 -0600 Subject: [PATCH 4/8] DGEMM Replace the ESMF version of DGEMM with the LAPACK version. --- spectral_transforms.F90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/spectral_transforms.F90 b/spectral_transforms.F90 index 71e14bb..8a89fa6 100644 --- a/spectral_transforms.F90 +++ b/spectral_transforms.F90 @@ -39,7 +39,7 @@ subroutine spec_to_four(flnev,flnod,plnev,plnod, & implicit none ! - external esmf_dgemm + external dgemm ! integer, intent(in) :: nvars real(kind=kind_dbl_prec) flnev(len_trie_ls,2*nvars) @@ -112,14 +112,14 @@ subroutine spec_to_four(flnev,flnod,plnev,plnod, & ! compute the sum of the even real terms for each level ! compute the sum of the even imaginary terms for each level ! - call esmf_dgemm('t', 'n', n2, latg2-lat1+1, (jcap+3-l)/2, & + call dgemm('t', 'n', n2, latg2-lat1+1, (jcap+3-l)/2, & cons1, flnev(indev,1), len_trie_ls, plnev(indev,lat1), & len_trie_ls, cons0, apev(1,lat1), n2 ) ! ! compute the sum of the odd real terms for each level ! compute the sum of the odd imaginary terms for each level ! - call esmf_dgemm('t', 'n', n2, latg2-lat1+1, (jcap+2-l)/2, & + call dgemm('t', 'n', n2, latg2-lat1+1, (jcap+2-l)/2, & cons1, flnod(indod,1), len_trio_ls, plnod(indod,lat1), & len_trio_ls, cons0, apod(1,lat1), n2 ) ! @@ -1539,6 +1539,8 @@ subroutine initialize_spectral(gis_stochy) allocate ( gis_stochy%epsodn(len_trio_ls) ) allocate ( gis_stochy%kenorm_e(len_trie_ls) ) allocate ( gis_stochy%kenorm_o(len_trio_ls) ) + allocate ( gis_stochy%gamma_e(len_trie_ls) ) + allocate ( gis_stochy%gamma_o(len_trio_ls) ) ! allocate ( gis_stochy%snnp1ev(len_trie_ls) ) allocate ( gis_stochy%snnp1od(len_trio_ls) ) @@ -1972,7 +1974,12 @@ subroutine glats_stochy(lgghaf,colrad,wgt,rcs2) cons2 = 2.d0, cons4 = 4.d0, & cons180 = 180.d0, & cons0p25 = 0.25d0 +#ifdef NO_QUAD_PRECISION + real(kind=kind_qdt_prec), parameter :: eps = 1.d-12 +#else real(kind=kind_qdt_prec), parameter :: eps = 1.d-20 +#endif + ! ! for better accuracy to select smaller number ! eps = 1.d-12 From 6065f09797c53f7342278551147cf73cce2b44d5 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Mon, 15 Jul 2024 07:55:14 -0600 Subject: [PATCH 5/8] Fix declarations Added declarations of `gamma_o` and `gamma_e` --- stochy_internal_state_mod.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/stochy_internal_state_mod.F90 b/stochy_internal_state_mod.F90 index cbd8c44..3032b27 100644 --- a/stochy_internal_state_mod.F90 +++ b/stochy_internal_state_mod.F90 @@ -51,6 +51,8 @@ module stochy_internal_state_mod real,allocatable :: epsodn(:) real,allocatable :: kenorm_e(:) real,allocatable :: kenorm_o(:) + real,allocatable :: gamma_e(:) + real,allocatable :: gamma_o(:) real,allocatable :: snnp1ev(:) real,allocatable :: snnp1od(:) From 8f26c01204883cb62c20ebac5e6f731c981a2a14 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Tue, 6 Aug 2024 10:46:01 -0600 Subject: [PATCH 6/8] NOAA Merge This commit merges my changes with the latest noaa-psl version. --- CMakeLists.txt | 13 +- cellular_automata_global.F90 | 33 ++-- cellular_automata_sgs.F90 | 126 +++++++++---- compns_stochy.F90 | 44 ++++- docs/source/references.rst | 2 + docs/source/users_guide.rst | 4 +- get_stochy_pattern.F90 | 11 +- kinddef.F90 | 18 +- lndp_apply_perts.F90 | 310 +++++++++++++++++++------------ main.doc | 2 +- mersenne_twister.F90 | 43 ++--- mpi_wrapper.F90 | 11 +- plumes.F90 | 5 + spectral_transforms.F90 | 99 +++++----- stochastic_physics.F90 | 102 +++++----- stochy_data_mod.F90 | 8 +- stochy_internal_state_mod.F90 | 36 ++-- stochy_namelist_def.F90 | 14 +- stochy_patterngenerator.F90 | 19 +- unit_tests/input.nml.ca_template | 4 +- unit_tests/run_unit_tests_ca.sh | 35 +++- unit_tests/standalone_ca.F90 | 10 +- update_ca.F90 | 238 +++++++++++++++++------- 23 files changed, 755 insertions(+), 432 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f010c9..ae14a65 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,14 @@ -if(32BIT) - message ("Force 64 bits in stochastic_physics") +if(CCPP_32BIT) + message(STATUS "Compile stochastic_physics with 32-bit precision to match CCPP slow physics.") + add_definitions(-DCCPP_32BIT) + if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -real-size 32") + elseif(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fno-default-real-8 -fdefault-double-8") + endif() +else() + message(STATUS "Compile stochastic_physics with 64-bit precision to match CCPP slow physics.") + remove_definitions(-DCCPP_32BIT) if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -real-size 64") elseif(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index a815aae..469db85 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -1,7 +1,7 @@ module cellular_automata_global_mod use update_ca, only : domain_global,iscnx_g,iecnx_g,jscnx_g,jecnx_g,isdnx_g,iednx_g,jsdnx_g,jednx_g, & - nxncells_g,nyncells_g,csum + nxncells_g,nyncells_g,csum,cold_start_ca_global implicit none contains @@ -11,12 +11,13 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp nca,ncells,nlives,nfracseed,nseed,iseed_ca, mytile, & ca_smooth,nspinup,blocksize,nsmooth,ca_amplitude,mpiroot,mpicomm) -use kinddef, only: kind_dbl_prec +use mpi_f08 +use kinddef, only: kind_dbl_prec, kind_phys use update_ca, only: update_cells_global,define_ca_domain use halo_exchange, only: atmosphere_scalar_field_halo use random_numbers, only: random_01_CB use mpp_domains_mod, only: domain2D,mpp_get_global_domain,CENTER, mpp_get_data_domain, mpp_get_compute_domain,mpp_global_sum, & - BITWISE_EFP_SUM, BITWISE_EXACT_SUM + BITWISE_EFP_SUM, BITWISE_EXACT_SUM,mpp_define_io_domain,mpp_get_io_domain_layout use block_control_mod, only: block_control_type, define_blocks_packed use mpi_wrapper, only: mp_reduce_sum,mp_reduce_max,mp_reduce_min, & mpi_wrapper_initialize,mype,is_rootpe @@ -30,13 +31,14 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp !This program evolves a cellular automaton uniform over the globe -integer, intent(in) :: kstep,ncells,nca,nlives,nseed,nspinup,nsmooth,mpiroot,mpicomm +integer, intent(in) :: kstep,ncells,nca,nlives,nseed,nspinup,nsmooth,mpiroot +type(MPI_Comm), intent(in) :: mpicomm integer(kind=kind_dbl_prec), intent(in) :: iseed_ca integer, intent(in) :: mytile -real(kind=kind_dbl_prec), intent(in) :: nfracseed,ca_amplitude +real(kind=kind_phys), intent(in) :: nfracseed,ca_amplitude logical, intent(in) :: ca_smooth,first_time_step, restart integer, intent(in) :: nblks,isc,iec,jsc,jec,npx,npy,nlev,blocksize -real(kind=kind_dbl_prec), intent(out) :: ca1_cpl(:,:),ca2_cpl(:,:),ca3_cpl(:,:) +real(kind=kind_phys), intent(out) :: ca1_cpl(:,:),ca2_cpl(:,:),ca3_cpl(:,:) type(domain2D), intent(inout) :: domain_in type(block_control_type) :: Atm_block integer :: nlon, nlat, isize,jsize,nf,nn @@ -50,8 +52,9 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp integer(8) :: count, count_rate, count_max, count_trunc,nx_full integer(8) :: iscale = 10000000000 integer, allocatable :: iini_g(:,:,:),ilives_g(:,:) -real(kind=kind_dbl_prec), allocatable :: field_out(:,:,:), field_smooth(:,:) -real(kind=kind_dbl_prec), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),CAprime(:,:) +integer, allocatable :: io_layout(:) +real(kind=kind_phys), allocatable :: field_out(:,:,:), field_smooth(:,:) +real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),CAprime(:,:) real*8 , allocatable :: noise(:,:,:) real*8 :: psum,CAmean,sq_diff,CAstdv,inv9 real*8 :: Detmax,Detmin @@ -73,7 +76,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp call mpi_wrapper_initialize(mpiroot,mpicomm) end if -halo=1 +halo=3 k_in=1 !---------------------------------------------------------------------------- @@ -102,8 +105,12 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp !Get CA domain if(first_time_step)then -! if (.not. restart) call define_ca_domain(domain_in,domain_global,ncells,nxncells_g,nyncells_g) - domain_global=domain_in + if (.not. restart) then + allocate(io_layout(2)) + io_layout=mpp_get_io_domain_layout(domain_in) + call define_ca_domain(domain_in,domain_global,halo,ncells,nxncells_g,nyncells_g) + call mpp_define_io_domain(domain_global, io_layout) + endif call mpp_get_data_domain (domain_global,isdnx_g,iednx_g,jsdnx_g,jednx_g) call mpp_get_compute_domain (domain_global,iscnx_g,iecnx_g,jscnx_g,jecnx_g) endif @@ -160,7 +167,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. - count4 = mod(((iseed_ca+7)*mytile)*(i1+nx_full*(j1-1))+ 2147483648, 4294967296) - 2147483648 + count4 = mod(((iseed_ca+7)*mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 endif ct=1 do nf=1,nca @@ -203,7 +210,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp CA(:,:)=0. - call update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,nyc,nxch,nych,nlon,nlat,isc,iec,jsc,jec, & + call update_cells_global(kstep,halo,first_time_step,iseed_ca,restart,nca,nxc,nyc,nxch,nych,nlon,nlat,isc,iec,jsc,jec, & npx,npy,CA,iini_g,ilives_g, & nlives,ncells,nfracseed,nseed,nspinup,nf,mytile) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index ba6f87a..dbfee82 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -1,22 +1,23 @@ module cellular_automata_sgs_mod -use update_ca, only : domain_global,domain_sgs,iscnx,iecnx,jscnx,jecnx,isdnx,iednx,jsdnx,jednx,nxncells,nyncells +use update_ca, only : domain_global,domain_sgs,iscnx,iecnx,jscnx,jecnx,isdnx,iednx,jsdnx,jednx,nxncells,nyncells,cold_start_ca_sgs implicit none contains -subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lake,condition_cpl, & +subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lake,uwind,vwind,height,dx,condition_cpl, & ca_deep_cpl,ca_turb_cpl,ca_shal_cpl,domain_in, & - nblks,isc,iec,jsc,jec,npx,npy,nlev,nthresh,rcell, mytile, & - nca,scells,tlives,nfracseed,nseed,iseed_ca, & + nblks,isc,iec,jsc,jec,npx,npy,nlev,nthresh, mytile, & + nca,ncells,nlives,nfracseed,nseed,iseed_ca,ca_advect, & nspinup,ca_trigger,blocksize,mpiroot,mpicomm) +use mpi_f08 use kinddef, only: kind_phys,kind_dbl_prec use update_ca, only: update_cells_sgs, define_ca_domain use random_numbers, only: random_01_CB use mpp_domains_mod, only: domain2D,mpp_get_global_domain,CENTER, mpp_get_data_domain, mpp_get_compute_domain,& - mpp_define_io_domain,mpp_get_io_domain_layout + mpp_define_io_domain,mpp_get_io_domain_layout,mpp_get_ntile_count use block_control_mod, only: block_control_type, define_blocks_packed use time_manager_mod, only: time_type use mpi_wrapper, only: mype,mp_reduce_max, & @@ -37,17 +38,20 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak ! swtich to new random number generator and improve computational efficiency ! and remove unsued code +!L.Bengtsson, 2023-05 +!Add horizontal advection of CA cells + !This routine produces an output field CA_DEEP for coupling to convection (saSAS). !CA_DEEP can be either number of plumes in a cluster (nca_plumes=true) or updraft !area fraction (nca_plumes=false) -integer,intent(in) :: kstep,scells,nca,tlives,nseed,nspinup,mpiroot,mpicomm,mytile +integer,intent(in) :: kstep,ncells,nca,nlives,nseed,nspinup,mpiroot,mytile +type(MPI_Comm),intent(in) :: mpicomm integer(kind=kind_dbl_prec), intent(in) :: iseed_ca -real(kind=kind_phys), intent(in) :: nfracseed,dtf,rcell -logical,intent(in) :: restart,ca_trigger,first_time_step +real(kind=kind_phys), intent(in) :: nfracseed,dtf,nthresh +logical,intent(in) :: restart,ca_trigger,first_time_step,ca_advect integer, intent(in) :: nblks,isc,iec,jsc,jec,npx,npy,nlev,blocksize -real , intent(out) :: nthresh -real(kind=kind_phys), intent(in) :: sst(:,:),lsmsk(:,:),lake(:,:) +real(kind=kind_phys), intent(in) :: sst(:,:),lsmsk(:,:),lake(:,:),uwind(:,:,:),dx(:,:),vwind(:,:,:),height(:,:,:) real(kind=kind_phys), intent(inout) :: condition_cpl(:,:) real(kind=kind_phys), intent(inout) :: ca_deep_cpl(:,:) real(kind=kind_phys), intent(inout) :: ca_turb_cpl(:,:) @@ -59,29 +63,29 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak integer :: inci, incj, nxc, nyc, nxch, nych, nx, ny integer :: halo, k_in, i, j, k integer :: seed, ierr7,blk, ix, iix, count4,ih,jh -integer :: blocksz,levs -integer :: ncells,nlives +integer :: blocksz,levs,u200,u850 integer, save :: initialize_ca integer(8) :: count, count_rate, count_max, count_trunc,nx_full integer(8) :: iscale = 10000000000 integer, allocatable :: iini(:,:,:),ilives_in(:,:,:),ca_plumes(:,:),io_layout(:) -real(kind=kind_phys), allocatable :: ssti(:,:),lsmski(:,:),lakei(:,:) -real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),conditiongrid(:,:) -real(kind=kind_phys), allocatable :: CA_DEEP(:,:) +real(kind=kind_phys), allocatable :: ssti(:,:),lsmski(:,:),lakei(:,:),uwindi(:,:),vwindi(:,:),dxi(:,:),heighti(:,:,:) +real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),uhigh(:,:),vhigh(:,:),dxhigh(:,:),conditiongrid(:,:) +real(kind=kind_phys), allocatable :: CA_DEEP(:,:),zi(:,:,:),sumx(:,:) real*8 , allocatable :: noise(:,:,:) -real(kind=kind_phys) :: condmax,condmaxinv,livesmax,livesmaxinv,factor,dx,pi,re +real(kind=kind_phys) :: condmax,condmaxinv,livesmax,livesmaxinv,factor,pi,re logical,save :: block_message=.true. logical :: nca_plumes logical,save :: first_flag integer*8 :: i1,j1 -integer :: ct +integer :: ct,ntiles +real :: dz,invgrav !nca :: switch for number of cellular automata to be used. ! :: for the moment only 1 CA can be used !nfracseed :: switch for number of random cells initially seeded -!tlives :: switch for time scale (s) +!nlives :: switch for time scale (s) !nspinup :: switch for number of itterations to spin up the ca -!scells :: switch for CA cell size (m) +!ncells :: switch for CA cell size (m) !nca_plumes :: compute number of CA-cells ("plumes") within a NWP gridbox. if (nca .LT. 1) return @@ -90,8 +94,13 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak call mpi_wrapper_initialize(mpiroot,mpicomm) end if -halo=1 +halo=3 k_in=1 +!Right now these values are experimental: +u200=56 +u850=13 +!Gravitational acceleration: +invgrav=1./9.81 nca_plumes = .true. @@ -117,15 +126,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak !Set time and length scales: call mpp_get_global_domain(domain_in,xsize=nx,ysize=ny,position=CENTER) - pi=3.14159 - re=6371000. - dx=0.5*pi*re/real(nx) - ncells=int(dx/real(scells)) - nlives=int(real(tlives)/dtf) - ncells = MIN(ncells,10) - nlives = MAX(nlives,5) - - nthresh=rcell*real(ncells)*real(ncells) + ntiles = mpp_get_ntile_count(domain_in) if(mype == 1)then write(*,*)'ncells=',ncells @@ -143,15 +144,12 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak if (.not.restart) then allocate(io_layout(2)) io_layout=mpp_get_io_domain_layout(domain_in) - call define_ca_domain(domain_in,domain_sgs,ncells,nxncells,nyncells) + call define_ca_domain(domain_in,domain_sgs,halo,ncells,nxncells,nyncells) call mpp_define_io_domain(domain_sgs, io_layout) endif call mpp_get_data_domain (domain_sgs,isdnx,iednx,jsdnx,jednx) call mpp_get_compute_domain (domain_sgs,iscnx,iecnx,jscnx,jecnx) - !write(1000+mpp_pe(),*) "nxncells,nyncells: ",nxncells,nyncells - !write(1000+mpp_pe(),*) "iscnx,iecnx,jscnx,jecnx: ",iscnx,iecnx,jscnx,jecnx - !write(1000+mpp_pe(),*) "isdnx,iednx,jsdnx,jednx: ",isdnx,iednx,jsdnx,jednx -endif + endif nxc = iecnx-iscnx+1 nyc = jecnx-jscnx+1 @@ -163,15 +161,29 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak allocate(ssti(nlon,nlat)) allocate(lsmski(nlon,nlat)) allocate(lakei(nlon,nlat)) + allocate(uwindi(nlon,nlat)) + allocate(vwindi(nlon,nlat)) + allocate(heighti(nlon,nlat,nlev)) + allocate(zi(nlon,nlat,nlev)) + allocate(sumx(nlon,nlat)) + allocate(dxi(nlon,nlat)) allocate(iini(nxc,nyc,nca)) allocate(ilives_in(nxc,nyc,nca)) allocate(condition(nxc,nyc)) + allocate(uhigh(nxc,nyc)) + allocate(dxhigh(nxc,nyc)) + allocate(vhigh(nxc,nyc)) allocate(conditiongrid(nlon,nlat)) allocate(CA(nlon,nlat)) allocate(ca_plumes(nlon,nlat)) allocate(CA_DEEP(nlon,nlat)) !Initialize: + ilives_in(:,:,:) = 0 + iini(:,:,:) = 0 + sumx(:,:) = 0. + uwindi(:,:) = 0. + vwindi(:,:) =0. !Put the blocks of model fields into a 2d array - can't use nlev and blocksize directly, !because the arguments to define_blocks_packed are intent(inout) and not intent(in). @@ -181,6 +193,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak call define_blocks_packed('cellular_automata', Atm_block, isc, iec, jsc, jec, levs, & blocksz, block_message) + do blk = 1,Atm_block%nblks do ix = 1, Atm_block%blksz(blk) i = Atm_block%index(blk)%ii(ix) - isc + 1 @@ -189,13 +202,39 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak ssti(i,j) = sst(blk,ix) lsmski(i,j) = lsmsk(blk,ix) lakei(i,j) = lake(blk,ix) + dxi(i,j) = dx(blk,ix) + do k = 1,levs + heighti(i,j,k) = height(blk,ix,k)*invgrav + enddo + do k = 2,levs + zi(i,j,k)=0.5*(heighti(i,j,k)+heighti(i,j,k-1)) + enddo + do k = u850,u200 + dz=zi(i,j,k)-zi(i,j,k-1) + uwindi(i,j) = uwindi(i,j) + uwind(blk,ix,k)*dz + vwindi(i,j) = vwindi(i,j) + vwind(blk,ix,k)*dz + sumx(i,j) = sumx(i,j) + dz + enddo enddo enddo + do blk = 1,Atm_block%nblks + do ix = 1, Atm_block%blksz(blk) + i = Atm_block%index(blk)%ii(ix) - isc + 1 + j = Atm_block%index(blk)%jj(ix) - jsc + 1 + uwindi(i,j)=uwindi(i,j)/sumx(i,j) + vwindi(i,j)=vwindi(i,j)/sumx(i,j) + enddo + enddo + + !Initialize the CA when the condition field is populated do j=1,nyc do i=1,nxc condition(i,j)=conditiongrid(inci/ncells,incj/ncells) + uhigh(i,j)=uwindi(inci/ncells,incj/ncells) + vhigh(i,j)=vwindi(inci/ncells,incj/ncells) + dxhigh(i,j)=dxi(inci/ncells,incj/ncells)/ncells !dx on the finer grid if(i.eq.inci)then inci=inci+ncells endif @@ -239,7 +278,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak endif !Generate random number, following stochastic physics code: -if (.not. restart) then +if (cold_start_ca_sgs) then if(kstep == initialize_ca) then nx_full=int(ncells,kind=8)*int(npx-1,kind=8) allocate(noise(nxc,nyc,nca)) @@ -257,7 +296,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. - count4 = mod((iseed_ca+mytile)*(i1+nx_full*(j1-1))+ 2147483648, 4294967296) - 2147483648 + count4 = mod((iseed_ca+mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 endif ct=1 do nf=1,nca @@ -282,13 +321,13 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak deallocate(noise) endif ! -endif ! restart +endif ! cold_start_ca_sgs !Calculate neighbours and update the automata do nf=1,nca - call update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,first_time_step,nca,nxc,nyc, & - nxch,nych,nlon,nlat,isc,iec,jsc,jec, & - npx,npy,CA,ca_plumes,iini,ilives_in, & + call update_cells_sgs(kstep,halo,dtf,initialize_ca,iseed_ca,first_flag,restart,first_time_step,nca,nxc,nyc, & + nxch,nych,nlon,nlat,isc,iec,jsc,jec,ca_advect, & + npx,npy,CA,ca_plumes,iini,ilives_in,uhigh,vhigh,dxhigh, & nlives,nfracseed,nseed,nspinup,nf,nca_plumes,ncells,mytile) if(nca_plumes)then @@ -312,6 +351,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak !Limit CA activity to the Tropical Ocean +if(ntiles==6)then do j=1,nlat do i=1,nlon if(ssti(i,j) < 300. .or. lsmski(i,j) /= 0. .or. lakei(i,j) > 0.0)then @@ -319,6 +359,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak endif enddo enddo +endif !Put back into blocks 1D array to be passed to physics !or diagnostics output @@ -337,6 +378,11 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak deallocate(lakei) deallocate(iini) deallocate(ilives_in) + deallocate(uwindi) + deallocate(vwindi) + deallocate(uhigh) + deallocate(vhigh) + deallocate(dxhigh) deallocate(condition) deallocate(CA) deallocate(ca_plumes) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 06be05c..e11aaf5 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -62,13 +62,13 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2, & shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt, & + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt,pbl_taper, & ocnskeb,ocnskeb_lscale,ocnskeb_tau,iseed_ocnskeb namelist /nam_sfcperts/lndp_type,lndp_model_type, lndp_var_list, lndp_prt_list, & - iseed_lndp, lndp_tau,lndp_lscale + iseed_lndp,lndpint,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 + sppint,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 @@ -83,7 +83,7 @@ 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_var_list = 'XXXXXXXXXX' spp_prt_list = -999. ! logicals do_sppt = .false. @@ -145,6 +145,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) spp_sigtop2 = 0.025 ! reduce amplitude of sppt near surface (lowest 2 levels) sppt_sfclimit = .false. + pbl_taper = (/0.0,0.5,1.0,1.0,1.0,1.0,1.0/) ! gaussian or power law variance spectrum for skeb (0: gaussian, 1: ! power law). If power law, skeb_lscale interpreted as a power not a ! length scale. @@ -156,6 +157,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) 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) + sppint = 0 ! SPP interval in seconds + lndpint = 0 ! lndp interval in seconds #ifdef INTERNAL_FILE_NML read(input_nml_file, nml=nam_stochy) @@ -213,6 +216,9 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) skeb=skeb*1.111e-9*sqrt(deltim) endif ENDIF + IF (spp_prt_list(1) > 0 ) THEN + do_spp=.true. + ENDIF ! compute frequencty to estimate dissipation timescale IF (do_skeb) THEN IF (skebint == 0.) skebint=deltim @@ -241,6 +247,26 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) return ENDIF ENDIF + IF (do_spp) THEN + IF (sppint == 0.) sppint=deltim + nsspp=nint(sppint/deltim) ! sppint in seconds + IF(nsspp<=0 .or. abs(nsspp-sppint/deltim)>tol) THEN + WRITE(0,*) "SPP interval is invalid",sppint + iret=9 + return + ENDIF + ENDIF + + IF (lndp_type > 0) THEN + IF (lndpint == 0.) lndpint=deltim + nslndp=nint(lndpint/deltim) + IF(nslndp<=0 .or. abs(nslndp-lndpint/deltim)>tol) THEN + WRITE(0,*) "lndp interval is invalid",lndpint + iret=9 + return + ENDIF + ENDIF + !calculate ntrunc if not supplied if (ntrunc .LT. 1) then if (me==0) print*,'ntrunc not supplied, calculating' @@ -253,8 +279,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) 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 + ntrunc=1.5*circ/l_min + !ntrunc=circ/l_min if (me==0) print*,'ntrunc calculated from l_min',l_min,ntrunc endif ! ensure lat_s is a mutiple of 4 with a reminader of two @@ -347,7 +373,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! 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 + if ( (spp_var_list(k) .EQ. 'XXXXXXXXXX') .or. (spp_prt_list(k) .LE. 0.) ) then cycle else n_var_spp=n_var_spp+1 @@ -368,7 +394,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) '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') + case('pbl','sfc', 'mp','rad','gwd','cu_deep') 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) @@ -446,7 +472,7 @@ subroutine compns_stochy_ocn (deltim,iret) skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2, & shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt, & + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt,pbl_taper, & ocnskeb,ocnskeb_lscale,ocnskeb_tau,iseed_ocnskeb namelist /nam_sfcperts/lndp_type,lndp_model_type,lndp_var_list, lndp_prt_list, iseed_lndp, & diff --git a/docs/source/references.rst b/docs/source/references.rst index ac2a261..736dd1e 100644 --- a/docs/source/references.rst +++ b/docs/source/references.rst @@ -3,6 +3,8 @@ References Berner, J., G. Shutts, M. Leutbecher, and T. Palmer, 2009: A spectral stochastic kinetic energy backscatter scheme and its impact on flow- dependent predictability in the ECMWF ensemble prediction system. J. Atmos. Sci., 66, 603–626, `doi:10.1175/2008JAS2677.1 `_ +Draper, C, 2021: Accounting for Land Model Uncertainty in Numerical Weather Prediction Ensemble Systems: Toward Ensemble-Based Coupled Land–Atmosphere Data Assimilation, J Hydromet, 22(8), 2089-2104, ``_ + Gehne, M., T. Hamill, G. Bates, P. Pegion, W. Kolczynski 2019: Land-surface parameter and state perturbations in the Global Ensemble Forecast System. Mon. Wea. Rev. 147, 1319–1340 `doi:10.1175/MWR-D-18-0057.1 `_ Palmer, T. N., R. Buizza, F. Doblas-Reyes, T. Jung, M. Leutbecher, G. J. Shutts,M. Steinheimer, and A.Weisheimer, 2009: Stochastic parametrization and model uncertainty. ECMWF Tech. Memo. 598, 42 pp `doi:10.21957/ps8gbwbdv `_ diff --git a/docs/source/users_guide.rst b/docs/source/users_guide.rst index 0d25818..6e03919 100644 --- a/docs/source/users_guide.rst +++ b/docs/source/users_guide.rst @@ -2,7 +2,7 @@ Users Guide ================================================== The stochastic physics currently only works with the UFS-atmosphere model -Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Gehne et al, 2019), and a cellular automata scheme (Bengtsson et al. 2019) which interacts directly with the convective parameterization. +Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Draper, 2021), and a cellular automata scheme (Bengtsson et al. 2019) which interacts directly with the convective parameterization. SKEB adds wind perturbations to model state. Perturbations are random in space/time, but amplitude is determined by a smoothed dissipation estimate provided by the dynamical core. Addresses errors in the dynamics - more active in the mid-latitudes @@ -11,7 +11,7 @@ SPPT multiplies the physics tendencies by a random number O [0,2] before updatin SHUM multiply the low-level specific humidity by a small random number each time-step. It attempts to address missing physics (cold pools, gust fronts), most active in convective regions -Land surface perturbations allow for land surface parameters such as Albedo, Soil Hydraulic Conductivity, LAI, and roughness lengths to vary in space. Addresses error in the land model and land-atmosphere interactions. +Land surface perturbations can be applied to the land model parameters and land model prognostic variables. The scheme is intended to address errors in the land model and land-atmosphere interactions. Due to the model’s numerics, any stochastic perturbation needs to be correlated in space and time in order to have the desired effect of upscale growth of the perturbations. This is achieved by creating a random pattern that has a specified decorrelation length-scale and is a first order auto-regressive process AR(1) in time with a specified decorrelation time-scale. (The CA random pattern generator also satisfies this condition) diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index 59d58f4..231d6dc 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -1,6 +1,6 @@ !>@brief The module 'get_stochy_pattern_mod' contains the subroutines to retrieve the random pattern in the cubed-sphere grid module get_stochy_pattern_mod - use kinddef, only : kind_dbl_prec, kind_evod + use kinddef, only : kind_dbl_prec use spectral_transforms, only : len_trie_ls, & len_trio_ls, ls_dim, stochy_la2ga, & coslat_a, latg, levs, lonf, skeblevs,& @@ -103,8 +103,8 @@ subroutine get_random_pattern_vector(rpattern,npatterns,& type(stochy_internal_state), intent(in) :: gis_stochy type(random_pattern), intent(inout) :: rpattern(npatterns) - real(kind=kind_evod), dimension(len_trie_ls,2) :: vrtspec_e,divspec_e - real(kind=kind_evod), dimension(len_trio_ls,2) :: vrtspec_o,divspec_o + real(kind=kind_dbl_prec), dimension(len_trie_ls,2) :: vrtspec_e,divspec_e + real(kind=kind_dbl_prec), dimension(len_trio_ls,2) :: vrtspec_o,divspec_o integer:: npatterns real(kind=kind_dbl_prec) :: upattern_3d(gis_stochy%nx,gis_stochy%ny,levs) @@ -115,7 +115,7 @@ subroutine get_random_pattern_vector(rpattern,npatterns,& ! logical lprint - real, allocatable, dimension(:,:) :: workgu,workgv + real(kind_dbl_prec), allocatable, dimension(:,:) :: workgu,workgv integer kmsk0(lonf,gis_stochy%lats_node_a) kmsk0 = 0 allocate(workgu(lonf,latg)) @@ -249,16 +249,13 @@ subroutine get_random_pattern_scalar(rpattern,npatterns,& do n=1,npatterns if (present(normalize)) then if (normalize) then - !print*,'calling patterngenerator_advance norm' call patterngenerator_advance_jb(rpattern(n)) call scalarspect_to_gaugrid_norm(rpattern(n),gis_stochy, wrk2d,1) else - !print*,'calling patterngenerator_advance' call patterngenerator_advance(rpattern(n),1,.false.) call scalarspect_to_gaugrid(rpattern(n),gis_stochy, wrk2d,1) endif else - !print*,'calling patterngenerator_advance' call patterngenerator_advance(rpattern(n),1,.false.) call scalarspect_to_gaugrid(rpattern(n),gis_stochy, wrk2d,1) endif diff --git a/kinddef.F90 b/kinddef.F90 index ef369d9..f427925 100644 --- a/kinddef.F90 +++ b/kinddef.F90 @@ -4,19 +4,21 @@ module kinddef private - public :: kind_evod, kind_phys + public :: kind_phys public :: kind_dbl_prec, kind_qdt_prec - public :: kind_io4, kind_io8 + public :: kind_io8 - integer, parameter :: kind_io4 = 4 - - ! DH* TODO - stochastic physics / CA should be using only one of these - integer, parameter :: kind_evod = 8 + ! kind_phys must match CCPP Physics kind_phys +#ifdef CCPP_32BIT + integer, parameter :: kind_phys = 4 +#else integer, parameter :: kind_phys = 8 +#endif + integer, parameter :: kind_dbl_prec = 8 - integer, parameter :: kind_io8 = 8 + integer, parameter :: kind_io8 = kind_dbl_prec -#ifdef __PGI +#ifdef NO_QUAD_PRECISION integer, parameter :: kind_qdt_prec = 8 #else integer, parameter :: kind_qdt_prec = 16 diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index bd3ef9f..835cdf7 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -1,6 +1,7 @@ module lndp_apply_perts_mod - use kinddef, only : kind_dbl_prec + use kinddef, only : kind_dbl_prec, kind_phys + use stochy_namelist_def implicit none @@ -13,95 +14,172 @@ module lndp_apply_perts_mod !==================================================================== ! lndp_apply_perts !==================================================================== -! Driver for applying perturbations to sprecified land states or parameters +! Driver for applying perturbations to specified land states or parameters, +! following the LNDP_TYPE=2 scheme. ! 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, & +! Can select perturbations by specifying lndp_var_list and lndp_pert_list + +! Notes on how the perturbations are added. +! 1. Model prognostic variables. +! If running a long forecast or cycling DA system (as in global UFS), +! perturbing the prognostic variables only at the initial conditions will +! have very limited impact, and they should instead be perturbed at every time step. +! In this case, the pertrubations should be specified as a rate (pert/hr) +! to avoid the ensemble spread being dependent on the model time step. +! +! For a short forecast (~days, as in regional HRRR), can see impact from +! perturbing only the initial conditions. In this case, the perturbation +! is specified as an absolute value (not a rate). +! +! 2. Model parameters: +! The timing of how to perturb the parameters depends on how / whether +! the parameters are updated over time. For the UFS global system, global_cycle +! is periodically called to update the parameters (controlled by FHCYC). +! Each time it's called global_cycle overwrites most of the +! prior parameters (overwriting any perturbations applied to those +! parameters). Hence, the perturbations are applied only immediately after global_cycle +! has been called, and the parameters are not applied as a rate (since they +! don't accumulate). +! For the regional models, FHCYC is 0, and the global_cycle is not called, so +! can perturb parameters every time step. Hence, need to specify the perturbations +! as a rate. +! +! The above cases are controlled by the lndp_model_type variable. +! +! If adding perturbations to new parameters, need to check how/whether +! the parameters are updated by the model. + + subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsm_noahmp, iopt_dveg, & + lsoil, dtf, kdt, 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) + smc, slc, stc, vfrac, alnsf, alnwf, snoalb, semis, zorll, ierr) implicit none ! 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 + integer, intent(in) :: n_var_lndp, lsoil, kdt, iopt_dveg + integer, intent(in) :: lsm, lsm_noah, lsm_ruc, lsm_noahmp 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(:,:) + real(kind=kind_phys), intent(in) :: lndp_prt_list(:) + real(kind=kind_phys), intent(in) :: dtf + real(kind=kind_phys), intent(in) :: sfc_wts(:,:,:) + real(kind=kind_phys), intent(in) :: xlon(:,:) + real(kind=kind_phys), 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) :: smcmax(:) - real(kind=kind_dbl_prec), intent(in) :: smcmin(:) + ! true = parameters have just been updated by global_cycle + integer, intent(in) :: stype(:,:) + real(kind=kind_phys), intent(in) :: smcmax(:) + real(kind=kind_phys), intent(in) :: smcmin(:) ! 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(:,:) - 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(:,:) + real(kind=kind_phys), intent(inout), optional :: smc(:,:,:) + real(kind=kind_phys), intent(inout), optional :: slc(:,:,:) + real(kind=kind_phys), intent(inout), optional :: stc(:,:,:) + real(kind=kind_phys), intent(inout), optional :: vfrac(:,:) + real(kind=kind_phys), intent(inout), optional :: snoalb(:,:) + real(kind=kind_phys), intent(inout), optional :: alnsf(:,:) + real(kind=kind_phys), intent(inout), optional :: alnwf(:,:) + real(kind=kind_phys), intent(inout), optional :: semis(:,:) + real(kind=kind_phys), intent(inout), optional :: 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 + integer :: this_im, v, k + logical :: print_flag, do_pert_state, do_pert_param - 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 + real(kind=kind_phys) :: p, min_bound, max_bound, pert + real(kind=kind_phys) :: tmp_smc + real(kind=kind_phys) :: conv_hr2tstep, tfactor_state, tfactor_param + real(kind=kind_phys), dimension(lsoil) :: zslayer, smc_vertscale, stc_vertscale ! decrease in applied pert with depth !-- 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_phys), dimension(4), parameter :: smc_vertscale_noah = (/1.0,0.5,0.25,0.125/) + real(kind=kind_phys), dimension(4), parameter :: stc_vertscale_noah = (/1.0,0.5,0.25,0.125/) + real(kind=kind_phys), dimension(4), parameter :: zs_noah = (/0.1, 0.3, 0.6, 1.0/) !-- 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/) + real(kind=kind_phys), 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_phys), 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_phys), 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,', & + if (lsm/=lsm_noah .and. lsm/=lsm_ruc .and. lsm/=lsm_noahmp) then + write(6,*) 'ERROR: lndp_apply_pert assumes LSM is Noah, Noah-MP, or RUC,', & ' may need to adapt variable names for a different LSM' ierr=10 return endif - !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. + if (lsm==lsm_noahmp) then + do v = 1,n_var_lndp + select case (trim(lndp_var_list(v))) + case ('alb','sal','emi','zol') + print*, & + 'ERROR: lndp_prt_list option in lndp_apply_pert', trim(lndp_var_list(v)) , & + ' has not been checked for Noah-MP. Please check how the parameter is set/updated ', & + ' before applying. Note: in Noah-MP many variables that have traditionally been', & + ' externally specified parameters are now prognostic. Also, many parameters are', & + ' set at each Noah-MP model call from data tables' + ierr = 10 + return + case ('veg') + if (iopt_dveg .NE. 4 ) then + print*, & + 'ERROR: veg perturbations have not yet been coded for dveg options other than 4' + ierr = 10 + return + endif + end select + enddo endif - if (lsm == lsm_noah) then + ! for perturbations applied as a rate, lndp_prt_list input is per hour. Converts to per timestep + conv_hr2tstep = dtf/3600. ! conversion factor from per hour to per tstep. + + ! determine whether updating state variables and/or parameters + + do_pert_state=.false. + do_pert_param=.false. + + select case (lndp_model_type) + case(1) ! global, perturb states every time step (pert applied as a rate) + ! perturb parameters only when they've been update (pert is not a rate) + do_pert_state=.true. + tfactor_state=conv_hr2tstep + if (param_update_flag) then + do_pert_param=.true. + tfactor_param=1. + endif + case(2) ! regional, perturb states only at first time step (pert is not a rate) + ! perurb parameters at every time step (pert is a rate) + if ( kdt == 2 ) then + do_pert_state=.true. + tfactor_state=1. + endif + do_pert_param = .true. + tfactor_param = conv_hr2tstep + case(3) ! special case to apply perturbations at initial time step only (pert is not a rate) + if ( kdt == 2 ) then + do_pert_state=.true. + tfactor_state=1. + do_pert_param=.true. + tfactor_param=1. + endif + case default + print*, & + 'ERROR: unrecognised lndp_model_type option in lndp_apply_pert, exiting', trim(lndp_var_list(v)) + ierr = 10 + return + end select + + if (lsm == lsm_noah .or. lsm == lsm_noahmp) then do k = 1, lsoil zslayer(k) = zs_noah(k) smc_vertscale(k) = smc_vertscale_noah(k) @@ -122,11 +200,6 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & 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 @@ -141,94 +214,102 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! State updates - performed every cycle !================================================================= case('smc') - p=5. - soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc - 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 - ! 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 - ! (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 - endif + if (do_pert_state) then + p=5. + min_bound = smcmin(stype(nb,i)) + max_bound = smcmax(stype(nb,i)) + + ! with RUC LSM perturb smc only at time step = 2, as in HRRR + do k=1,lsoil + ! apply perts to a copy of smc, retain original smc + ! for later update to liquid soil moisture. + ! note: previously we were saving the ice water content + ! (smc-slc) and subtracting this from the perturbed smc to + ! get the perturbed slc. This was introducing small errors in the slc + ! when passing back to the calling program, I think due to precision issues, + ! as the ice content is typically zero. Clara Draper, March, 2022. + tmp_smc = smc(nb,i,k) + + ! perturb total soil moisture + ! factor of sldepth*1000 converts from mm to m3/m3 + ! NOTE: smc in the pertlist specified in input.nml + ! is in the unit of mm/hour + pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zslayer(k)*1000.) + pert = pert*tfactor_state + + call apply_pert('smc',pert,print_flag, tmp_smc,ierr,p,min_bound, max_bound) + + ! assign all of applied pert to the liquid soil moisture + slc(nb,i,k) = slc(nb,i,k) + tmp_smc - smc(nb,i,k) + smc(nb,i,k) = tmp_smc + + enddo + endif case('stc') + if (do_pert_state) then + do k=1,lsoil + pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) + pert = pert*tfactor_state + call apply_pert('stc',pert,print_flag, stc(nb,i,k),ierr) + enddo + endif - 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 + ! Parameter updates !================================================================= + + ! are all of the params below included in noah? + case('vgf') ! vegetation fraction - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) then p =5. min_bound=0. max_bound=1. pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) endif case('alb') ! albedo - if (param_update_flag .or. lndp_each_step) then + if (do_pert_param) 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) + pert = pert*tfactor_param 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 + if (do_pert_param) then p =5. min_bound=0.3 max_bound=0.85 pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param 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 + if (do_pert_param) then p =5. min_bound=0.8 max_bound=1. pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param 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 + if (do_pert_param) then p =5. min_bound=0. max_bound=300. pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - pert = pert*factor + pert = pert*tfactor_param call apply_pert ('zol',pert,print_flag, zorll(nb,i), ierr,p,min_bound, max_bound) endif case default @@ -251,22 +332,22 @@ 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 + real(kind=kind_phys), 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 + real(kind=kind_phys), 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 + real(kind=kind_phys), optional, intent(in) :: vmin, vmax ! min,max bounds of variable being perturbed ! intent (inout) - real(kind=kind_dbl_prec), intent(inout) :: state + real(kind=kind_phys), intent(inout) :: state ! intent out integer :: ierr !local - real(kind=kind_dbl_prec) :: z + real(kind=kind_phys) :: z if ( print_flag ) then write(*,*) 'LNDP - applying lndp to ',vname, ', initial value', state @@ -282,12 +363,11 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax) z = -1. + 2*(state - vmin)/(vmax - vmin) ! flat-top function state = state + pert*(1-abs(z**p)) else - state = state + pert + 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 @@ -307,8 +387,8 @@ subroutine set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb) ! 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_phys), intent(in) :: xlon(:,:) + real(kind=kind_phys), intent(in) :: xlat(:,:) ! intent (out) diff --git a/main.doc b/main.doc index eeedbef..dfcdb1e 100644 --- a/main.doc +++ b/main.doc @@ -6,7 +6,7 @@ * *The stochastic physics currently only works with the UFS-atmosphere model * - *Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Gehne et al, 2019), and a cellular automata scheme (Bengtsson et al. 20XX) which interacts directly with the convective parameterization. + *Currently, 3 stochastic schemes are used operationally at NCEP/EMC: Stochastic Kinetic Energy Backscatter (SKEB; Berner et al., 2009), Stochastically Perturbed Physics Tendencies (SPPT; Palmer et al., 2009), and Specific Humidity perturbations (SHUM), which is inspired by Tompkins and Berner, 2008. In addition there is the ability to perturb certain land model/surface parameters (Draper, 2021), and a cellular automata scheme (Bengtsson et al. 20XX) which interacts directly with the convective parameterization. * *SKEB adds wind perturbations to model state. Perturbations are random in space/time, but amplitude is determined by a smoothed dissipation estimate provided by the dynamical core. *Addresses errors in the dynamics - more active in the mid-latitudes diff --git a/mersenne_twister.F90 b/mersenne_twister.F90 index 9ab6103..c537681 100644 --- a/mersenne_twister.F90 +++ b/mersenne_twister.F90 @@ -160,6 +160,7 @@ ! !$$$ module mersenne_twister + use kinddef, only: kind_dbl_prec private ! Public declarations public random_stat @@ -188,7 +189,7 @@ module mersenne_twister integer:: mti=n+1 integer:: mt(0:n-1) integer:: iset - real:: gset + real(kind_dbl_prec):: gset end type ! Saved data type(random_stat),save:: sstat @@ -296,8 +297,8 @@ subroutine random_setseed_t(inseed,stat) !> This function generates random numbers in functional mode. function random_number_f() result(harvest) implicit none - real:: harvest - real h(1) + real(kind_dbl_prec):: harvest + real(kind_dbl_prec) :: h(1) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_number_t(h,sstat) harvest=h(1) @@ -306,7 +307,7 @@ function random_number_f() result(harvest) !> This subroutine generates random numbers in interactive mode. subroutine random_number_i(harvest,inseed) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) integer,intent(in):: inseed type(random_stat) stat call random_setseed_t(inseed,stat) @@ -316,7 +317,7 @@ subroutine random_number_i(harvest,inseed) !> This subroutine generates random numbers in saved mode; overloads Fortran 90 standard. subroutine random_number_s(harvest) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_number_t(harvest,sstat) end subroutine @@ -324,7 +325,7 @@ subroutine random_number_s(harvest) !> This subroutine generates random numbers in thread-safe mode. subroutine random_number_t(harvest,stat) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) type(random_stat),intent(inout):: stat integer j,kk,y integer tshftu,tshfts,tshftt,tshftl @@ -352,9 +353,9 @@ subroutine random_number_t(harvest,stat) y=ieor(y,iand(tshftt(y),tmaskc)) y=ieor(y,tshftl(y)) if(y.lt.0) then - harvest(j)=(real(y)+2.0**32)/(2.0**32-1.0) + harvest(j)=(real(y,kind_dbl_prec)+2.0_kind_dbl_prec**32)/(2.0_kind_dbl_prec**32-1.0_kind_dbl_prec) else - harvest(j)=real(y)/(2.0**32-1.0) + harvest(j)=real(y,kind_dbl_prec)/(2.0_kind_dbl_prec**32-1.0_kind_dbl_prec) endif stat%mti=stat%mti+1 enddo @@ -363,8 +364,8 @@ subroutine random_number_t(harvest,stat) !> This subrouitne generates Gaussian random numbers in functional mode. function random_gauss_f() result(harvest) implicit none - real:: harvest - real h(1) + real(kind_dbl_prec):: harvest + real(kind_dbl_prec) :: h(1) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_gauss_t(h,sstat) harvest=h(1) @@ -373,7 +374,7 @@ function random_gauss_f() result(harvest) !> This subrouitne generates Gaussian random numbers in interactive mode. subroutine random_gauss_i(harvest,inseed) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) integer,intent(in):: inseed type(random_stat) stat call random_setseed_t(inseed,stat) @@ -383,7 +384,7 @@ subroutine random_gauss_i(harvest,inseed) !> This subroutine generates Gaussian random numbers in saved mode. subroutine random_gauss_s(harvest) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_gauss_t(harvest,sstat) end subroutine @@ -391,10 +392,10 @@ subroutine random_gauss_s(harvest) !> This subroutine generates Gaussian random numbers in thread-safe mode. subroutine random_gauss_t(harvest,stat) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) type(random_stat),intent(inout):: stat integer mx,my,mz,j - real r2(2),r,g1,g2 + real(kind_dbl_prec) :: r2(2),r,g1,g2 mz=size(harvest) if(mz.le.0) return mx=0 @@ -429,14 +430,14 @@ subroutine random_gauss_t(harvest,stat) contains !> This subroutine contains numerical Recipes algorithm to generate Gaussian random numbers. subroutine rgauss(r1,r2,r,g1,g2) - real,intent(in):: r1,r2 - real,intent(out):: r,g1,g2 - real v1,v2,fac - v1=2.*r1-1. - v2=2.*r2-1. + real(kind_dbl_prec),intent(in):: r1,r2 + real(kind_dbl_prec),intent(out):: r,g1,g2 + real(kind_dbl_prec) :: v1,v2,fac + v1=2.*r1-1._kind_dbl_prec + v2=2.*r2-1._kind_dbl_prec r=v1**2+v2**2 if(r.lt.1.) then - fac=sqrt(-2.*log(r)/r) + fac=sqrt(-2._kind_dbl_prec*log(r)/r) g1=v1*fac g2=v2*fac endif @@ -482,7 +483,7 @@ subroutine random_index_t(imax,iharvest,stat) type(random_stat),intent(inout):: stat integer,parameter:: mh=n integer i1,i2,mz - real h(mh) + real(kind_dbl_prec) :: h(mh) mz=size(iharvest) do i1=1,mz,mh i2=min((i1-1)+mh,mz) diff --git a/mpi_wrapper.F90 b/mpi_wrapper.F90 index e756678..7ccd5a7 100644 --- a/mpi_wrapper.F90 +++ b/mpi_wrapper.F90 @@ -1,5 +1,7 @@ module mpi_wrapper + use mpi_f08 + implicit none private @@ -9,12 +11,10 @@ module mpi_wrapper public :: mp_reduce_min, mp_reduce_max, mp_reduce_sum public :: mp_bcst, mp_alltoall -#include "mpif.h" - integer, save :: mype = -999 integer, save :: npes = -999 integer, save :: root = -999 - integer, save :: comm = -999 + type(MPI_Comm), save :: comm logical, save :: initialized = .false. integer :: ierror @@ -78,7 +78,8 @@ logical function is_rootpe() end function is_rootpe subroutine mpi_wrapper_initialize(mpiroot, mpicomm) - integer, intent(in) :: mpiroot, mpicomm + integer, intent(in) :: mpiroot + type(MPI_Comm), intent(in) :: mpicomm if (initialized) return root = mpiroot comm = mpicomm @@ -92,7 +93,7 @@ subroutine mpi_wrapper_finalize() mype = -999 npes = -999 root = -999 - comm = -999 + comm%mpi_val = -999 initialized = .false. end subroutine mpi_wrapper_finalize diff --git a/plumes.F90 b/plumes.F90 index 80b301f..7f0cbdf 100644 --- a/plumes.F90 +++ b/plumes.F90 @@ -1,3 +1,6 @@ +module plumes_mod + +contains subroutine plumes(V,L,AG,a,row,col,kend) implicit none @@ -165,3 +168,5 @@ subroutine plumes(V,L,AG,a,row,col,kend) end subroutine plumes + +end module plumes_mod diff --git a/spectral_transforms.F90 b/spectral_transforms.F90 index 8a89fa6..5eed1cc 100644 --- a/spectral_transforms.F90 +++ b/spectral_transforms.F90 @@ -20,7 +20,7 @@ module spectral_transforms ! integer, public, allocatable :: lat1s_a(:), lon_dims_a(:) - real, public, allocatable, dimension(:) :: colrad_a, wgt_a, rcs2_a, & + real(kind_dbl_prec), public, allocatable, dimension(:) :: colrad_a, wgt_a, rcs2_a, & sinlat_a, coslat_a @@ -268,7 +268,7 @@ SUBROUTINE dcrft_stochy(init,x,ldx,y,ldy,n,nvars, table,n1,wrk,n2) implicit none integer ,intent(in) :: ldx,ldy,n,nvars integer init,n1,n2,i,j - real x(ldx,nvars),y(ldy,nvars),table(44002),wrk + real(kind_dbl_prec) x(ldx,nvars),y(ldy,nvars),table(44002),wrk IF (init.ne.0) THEN CALL rffti_stochy(n,table) @@ -297,8 +297,8 @@ SUBROUTINE RFFTB_STOCHY (N,R,WSAVE) implicit none - real, intent(inout) :: R(:) - real, intent(inout) :: WSAVE(44002) + real(kind_dbl_prec), intent(inout) :: R(:) + real(kind_dbl_prec), intent(inout) :: WSAVE(44002) integer :: N @@ -311,7 +311,7 @@ SUBROUTINE RFFTI_STOCHY (N,WSAVE) implicit none - REAL, intent(inout) :: WSAVE(44002) + REAL(kind_dbl_prec), intent(inout) :: WSAVE(44002) integer :: N IF (N .EQ. 1) RETURN @@ -325,10 +325,10 @@ SUBROUTINE RFFTB1_STOCHY (N,C,CH,WA,RFAC) implicit none integer, intent(in) :: N - real, intent(inout) :: CH(44002) - real, intent(inout) :: C(:) - real, intent(inout) :: WA(:) - real, intent(inout) :: RFAC(:) + real(kind_dbl_prec), intent(inout) :: CH(44002) + real(kind_dbl_prec), intent(inout) :: C(:) + real(kind_dbl_prec), intent(inout) :: WA(:) + real(kind_dbl_prec), intent(inout) :: RFAC(:) integer :: NF,NA,L1,IW,IP,L2,IDO,IDL1,IX2,IX3,IX4 integer :: K1,I @@ -397,14 +397,14 @@ SUBROUTINE RFFTI1_STOCHY (N,WA,RFAC) implicit none integer, intent(in) :: N - REAL, intent(inout) :: WA(:) - REAL, intent(inout) :: RFAC(:) + REAL(kind_dbl_prec), intent(inout) :: WA(:) + REAL(kind_dbl_prec), intent(inout) :: RFAC(:) integer :: NTRYH(4) integer :: NL,NF, I, J, NQ,NR,LD,FI,IS,ID,L1,L2,IP integer :: NTRY, NFM1, K1,II, IB, IDO, IPM, IC - REAL, parameter :: TPI=6.28318530717959 - real :: ARG,ARGLD,ARGH, TI2,TI4 + REAL(kind_dbl_prec), parameter :: TPI=6.28318530717959 + real(kind_dbl_prec) :: ARG,ARGLD,ARGH, TI2,TI4 DATA NTRYH(:) /4,2,3,5/ @@ -480,12 +480,12 @@ SUBROUTINE RADB2_STOCHY (IDO,L1,CC,CH,WA1) integer, intent(in) :: IDO integer, intent(in) :: L1 - real, intent(inout) :: CC(IDO,2,L1) - real, intent(inout) :: CH(IDO,L1,2) - real, intent(inout) :: WA1(:) + real(kind_dbl_prec), intent(inout) :: CC(IDO,2,L1) + real(kind_dbl_prec), intent(inout) :: CH(IDO,L1,2) + real(kind_dbl_prec), intent(inout) :: WA1(:) integer :: K,I,IC,IDP2 - real :: TR2,TI2 + real(kind_dbl_prec) :: TR2,TI2 DO 101 K=1,L1 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K) CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K) @@ -524,16 +524,16 @@ SUBROUTINE RADB3_STOCHY (IDO,L1,CC,CH,WA1,WA2) implicit none integer, intent(in) :: IDO,L1 - real, intent(inout) :: CC(IDO,3,L1) - real, intent(inout) :: CH(IDO,L1,3) - real, intent(inout) :: WA1(:) - real, intent(inout) :: WA2(:) + real(kind_dbl_prec), intent(inout) :: CC(IDO,3,L1) + real(kind_dbl_prec), intent(inout) :: CH(IDO,L1,3) + real(kind_dbl_prec), intent(inout) :: WA1(:) + real(kind_dbl_prec), intent(inout) :: WA2(:) - REAL, parameter :: TAUR= -.5 - REAL, parameter :: TAUI=.866025403784439 + REAL(kind_dbl_prec), parameter :: TAUR= -.5 + REAL(kind_dbl_prec), parameter :: TAUI=.866025403784439 integer :: I,K,IDP2,IC - real :: TR2,CR2,TI1,CI2,CR3,CI3,DR2,DR3,DI2,DI3 - real :: TI2,TI4 + real(kind_dbl_prec) :: TR2,CR2,TI1,CI2,CR3,CI3,DR2,DR3,DI2,DI3 + real(kind_dbl_prec) :: TI2,TI4 DO 101 K=1,L1 @@ -577,16 +577,16 @@ SUBROUTINE RADB4_STOCHY (IDO,L1,CC,CH,WA1,WA2,WA3) implicit none integer, intent(in) :: IDO,L1 - real, intent(inout) :: CC(IDO,4,L1) - real, intent(inout) :: CH(IDO,L1,4) - real, intent(inout) :: WA1(:) - real, intent(inout) :: WA2(:) - real, intent(inout) :: WA3(:) + real(kind_dbl_prec), intent(inout) :: CC(IDO,4,L1) + real(kind_dbl_prec), intent(inout) :: CH(IDO,L1,4) + real(kind_dbl_prec), intent(inout) :: WA1(:) + real(kind_dbl_prec), intent(inout) :: WA2(:) + real(kind_dbl_prec), intent(inout) :: WA3(:) - REAL, parameter :: SQRT2=1.414213562373095 + REAL(kind_dbl_prec), parameter :: SQRT2=1.414213562373095 integer :: I,K,IDP2,IC - real :: TR1,TR2,TR3,TR4,TI1,TI2,TI3,TI4 - real :: CI2,CI3,CI4,CR2,CR3,CR4 + real(kind_dbl_prec) :: TR1,TR2,TR3,TR4,TI1,TI2,TI3,TI4 + real(kind_dbl_prec) :: CI2,CI3,CI4,CR2,CR3,CR4 DO 101 K=1,L1 TR1 = CC(1,1,K)-CC(IDO,4,K) TR2 = CC(1,1,K)+CC(IDO,4,K) @@ -650,11 +650,19 @@ SUBROUTINE RADB4_STOCHY (IDO,L1,CC,CH,WA1,WA2,WA3) SUBROUTINE RADB5_STOCHY (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - DIMENSION CC(IDO,5,L1), CH(IDO,L1,5), WA1(*), WA2(*), WA3(*), WA4(*) - REAL, parameter :: TR11=0.309016994374947 - REAL, parameter :: TI11= 0.951056516295154 - REAL, parameter :: TR12=-0.809016994374947 - REAL, parameter :: TI12=0.587785252292473 + implicit none + integer, intent(in) :: L1, IDO + REAL(kind_dbl_prec) :: CC(IDO,5,L1), CH(IDO,L1,5), WA1(*), WA2(*), WA3(*), WA4(*) + REAL(kind_dbl_prec), parameter :: TR11=0.309016994374947 + REAL(kind_dbl_prec), parameter :: TI11= 0.951056516295154 + REAL(kind_dbl_prec), parameter :: TR12=-0.809016994374947 + REAL(kind_dbl_prec), parameter :: TI12=0.587785252292473 + integer :: k,IDP2,I,IC + real(kind_dbl_prec) :: & + TI5,TI4,TR2,TR3,CR2,CR3,CI5,CI4, & + TI2,TI3,TR5,TR4, & + CI2,CI3,CR5,CR4, & + DR3,DR4,DI3,DI4,DR5,DR2,DI5,DI2 DO 101 K=1,L1 TI5 = CC(1,3,K)+CC(1,3,K) TI4 = CC(1,5,K)+CC(1,5,K) @@ -714,11 +722,14 @@ SUBROUTINE RADB5_STOCHY (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) RETURN END - SUBROUTINE RADBG_STOCHY (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - DIMENSION CH(IDO,L1,IP), CC(IDO,IP,L1), C1(IDO,L1,IP), C2(IDL1,IP), & + implicit none + INTEGER :: IDO,IP,L1,IDL1 + REAL(kind_dbl_prec) :: CH(IDO,L1,IP), CC(IDO,IP,L1), C1(IDO,L1,IP), C2(IDL1,IP), & CH2(IDL1,IP) , WA(*) - REAL, parameter :: TPI=6.28318530717959 + REAL(kind_dbl_prec), parameter :: TPI=6.28318530717959 + REAL(kind_dbl_prec) :: ARG, DCP, DSP, AI1, AI2, AR1, AR1H, DC2, DS2, AR2, AR2H + INTEGER :: I,J,K,IK, IDP2, IPP2, IPPH, JC, J2, IC, L, IS, IDIJ, NBD, LC ARG = TPI/FLOAT(IP) DCP = COS(ARG) DSP = SIN(ARG) @@ -912,7 +923,7 @@ subroutine dozeuv_stochy(dod,zev,uod,vev,epsedn,epsodn, snnp1ev,snnp1od,ls_node) real(kind_dbl_prec) cons0 !constant integer indlsev,jbasev integer indlsod,jbasod - real(kind_evod) rerth + real(kind_dbl_prec) rerth include 'function2' @@ -1096,7 +1107,7 @@ subroutine dezouv_stochy(dev,zod,uev,vod,epsedn,epsodn,snnp1ev,snnp1od,ls_node) integer indlsev,jbasev integer indlsod,jbasod - real(kind_evod) rerth + real(kind_dbl_prec) rerth include 'function2' !...................................................................... @@ -1804,7 +1815,7 @@ subroutine gozrineo_a_stochy(gis_stochy, num_lat) real(kind=kind_dbl_prec) cons0 !constant real(kind=kind_dbl_prec) cons2 !constant - real rerth + real(kind_dbl_prec) rerth integer indlsev,jbasev integer indlsod,jbasod diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index edabaab..321becb 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -2,7 +2,8 @@ !! the stochastic physics random pattern generators module stochastic_physics -use kinddef, only : kind_dbl_prec +use mpi_f08 +use kinddef, only : kind_phys, kind_dbl_prec implicit none @@ -40,7 +41,8 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in ! Interface variables -integer, intent(in) :: levs, nlunit, nthreads, mpiroot, mpicomm +integer, intent(in) :: levs, nlunit, nthreads, mpiroot +type(MPI_Comm), intent(in) :: mpicomm integer, intent(in) :: blksz(:) real(kind=kind_dbl_prec), intent(in) :: dtp real(kind=kind_dbl_prec), intent(out) :: sppt_amp @@ -54,11 +56,11 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_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(:), intent(out) :: lndp_var_list_out -real(kind=kind_dbl_prec), dimension(:), intent(out) :: lndp_prt_list_out -character(len=3), dimension(:), intent(out) :: spp_var_list_out -real(kind=kind_dbl_prec), dimension(:), intent(out) :: spp_prt_list_out -real(kind=kind_dbl_prec), dimension(:), intent(out) :: spp_stddev_cutoff_out +character(len=3), optional, dimension(:), intent(out) :: lndp_var_list_out +real(kind=kind_phys), optional, dimension(:), intent(out) :: lndp_prt_list_out +character(len=10), optional, dimension(:), intent(out) :: spp_var_list_out +real(kind=kind_phys), optional, dimension(:), intent(out) :: spp_prt_list_out +real(kind=kind_phys), optional, dimension(:), intent(out) :: spp_stddev_cutoff_out ! Local variables @@ -162,8 +164,9 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in endif enddo if (sppt_sfclimit) then - vfact_sppt(2)=vfact_sppt(3)*0.5 - vfact_sppt(1)=0.0 + do k=1,7 + vfact_sppt(k)=pbl_taper(k) + enddo endif if (is_rootpe()) then do k=1,levs @@ -278,8 +281,10 @@ subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_i real :: dx integer :: k,latghf,km +type(MPI_Comm) :: mpicomm_t ! FIXME once MOM6 updates to use mpi_f90 types rad2deg=180.0/con_pi -call mpi_wrapper_initialize(mpiroot,mpicomm) +mpicomm_t%mpi_val = mpicomm +call mpi_wrapper_initialize(mpiroot,mpicomm_t) gis_stochy_ocn%nodes = npes gis_stochy_ocn%mype = mype gis_stochy_ocn%nx=nx @@ -347,7 +352,7 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s 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_spp -use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,sppt_logit, & +use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,nsspp,nslndp,sppt_logit, & lndp_type, n_var_lndp, n_var_spp, do_spp, spp_stddev_cutoff, spp_prt_list use mpi_wrapper, only: is_rootpe implicit none @@ -356,15 +361,15 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s integer, intent(in) :: levs, kdt real(kind=kind_dbl_prec), intent(in) :: fhour integer, intent(in) :: blksz(:) -real(kind=kind_dbl_prec), intent(inout) :: sppt_wts(:,:,:) -real(kind=kind_dbl_prec), intent(inout) :: shum_wts(:,:,:) -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(:,:,:,:) +real(kind=kind_phys), intent(inout), optional :: sppt_wts(:,:,:) +real(kind=kind_phys), intent(inout), optional :: shum_wts(:,:,:) +real(kind=kind_phys), intent(inout), optional :: skebu_wts(:,:,:) +real(kind=kind_phys), intent(inout), optional :: skebv_wts(:,:,:) +real(kind=kind_phys), intent(inout), optional :: sfc_wts(:,:,:) +real(kind=kind_phys), intent(inout), optional :: spp_wts(:,:,:,:) integer, intent(in) :: nthreads -real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:),tmpl_wts(:,:,:),tmp_spp_wts(:,:,:) +real(kind_dbl_prec),allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:),tmpl_wts(:,:,:),tmp_spp_wts(:,:,:) !D-grid integer :: k,v integer j,ierr,i @@ -436,6 +441,7 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s endif if ( lndp_type .EQ. 2 ) then ! add time check? + if (mod(kdt,nslndp) == 1 .or. nslndp == 1) then allocate(tmpl_wts(gis_stochy%nx,gis_stochy%ny,n_var_lndp)) call get_random_pattern_sfc(rpattern_sfc,nlndp,gis_stochy,tmpl_wts) DO blk=1,nblks @@ -446,23 +452,26 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s ENDDO ENDDO deallocate(tmpl_wts) + endif endif if (n_var_spp .GE. 1) then - 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_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 - deallocate(tmp_spp_wts) + if (mod(kdt,nsspp) == 1 .or. nsspp == 1) then + 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_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 + deallocate(tmp_spp_wts) + end if endif deallocate(tmp_wts) deallocate(tmpu_wts) @@ -478,7 +487,7 @@ subroutine run_stochastic_physics_ocn(sppt_wts,skeb_wts,t_rp1,t_rp2) use stochy_namelist_def implicit none real, intent(inout) :: sppt_wts(:,:),t_rp1(:,:),t_rp2(:,:),skeb_wts(:,:) -real, allocatable :: tmp_wts(:,:) +real(kind_dbl_prec), allocatable :: tmp_wts(:,:) if (pert_epbl .OR. do_ocnsppt .OR. do_ocnskeb ) then allocate(tmp_wts(gis_stochy_ocn%nx,gis_stochy_ocn%ny)) if (pert_epbl) then @@ -486,28 +495,27 @@ subroutine run_stochastic_physics_ocn(sppt_wts,skeb_wts,t_rp1,t_rp2) t_rp1(:,:)=2.0/(1+exp(-1*tmp_wts)) call get_random_pattern_scalar(rpattern_epbl2,nepbl,gis_stochy_ocn,tmp_wts) t_rp2(:,:)=2.0/(1+exp(-1*tmp_wts)) -! else -! t_rp1(:,:)=1.0 -! t_rp2(:,:)=1.0 + else + t_rp1(:,:)=1.0 + t_rp2(:,:)=1.0 endif if (do_ocnsppt) then call get_random_pattern_scalar(rpattern_ocnsppt,nocnsppt,gis_stochy_ocn,tmp_wts) sppt_wts=2.0/(1+exp(-1*tmp_wts)) -! else -! sppt_wts=1.0 + else + sppt_wts=1.0 endif if (do_ocnskeb) then call get_random_pattern_scalar(rpattern_ocnskeb,nocnskeb,gis_stochy_ocn,skeb_wts,normalize=.true.) -! else -! skeb_wts=1.0 + else + skeb_wts=1.0 endif - !print*,'after get_random_pattern_scalar',skeb_wts(1,1) deallocate(tmp_wts) -!else -! sppt_wts(:,:)=1.0 -! skeb_wts(:,:)=1.0 -! t_rp1(:,:)=1.0 -! t_rp2(:,:)=1.0 +else + sppt_wts(:,:)=1.0 + skeb_wts(:,:)=1.0 + t_rp1(:,:)=1.0 + t_rp2(:,:)=1.0 endif end subroutine run_stochastic_physics_ocn diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index 3d10b39..95b8518 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -17,6 +17,8 @@ module stochy_data_mod use mersenne_twister, only : random_seed use compns_stochy_mod, only : compns_stochy + use kinddef, only: kind_phys, kind_dbl_prec + implicit none private public :: init_stochdata,init_stochdata_ocn @@ -56,9 +58,9 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) integer, intent(in) :: nlunit,nlevs character(len=*), intent(in) :: input_nml_file(:) character(len=64), intent(in) :: fn_nml - real, intent(in) :: delt + real(kind_phys), intent(in) :: delt integer, intent(out) :: iret - real :: ones(5) + real :: ones(6) real :: rnn1,gamma_sum integer :: nn,k,nm,stochlun,ierr,n @@ -75,7 +77,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) iret=0 ! read in namelist - call compns_stochy (mype,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,delt,iret) + call compns_stochy (mype,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,real(delt,kind=kind_phys),iret) if (iret/=0) return ! need to make sure that non-zero irets are being trapped. if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0) .AND. (.NOT. do_spp)) return diff --git a/stochy_internal_state_mod.F90 b/stochy_internal_state_mod.F90 index 3032b27..093e2c7 100644 --- a/stochy_internal_state_mod.F90 +++ b/stochy_internal_state_mod.F90 @@ -21,7 +21,7 @@ module stochy_internal_state_mod !------ ! use spectral_layout_mod - + use kinddef implicit none private @@ -45,26 +45,26 @@ module stochy_internal_state_mod integer ,allocatable :: global_lats_h (:) integer :: xhalo,yhalo - real,allocatable :: epse (:) - real,allocatable :: epso (:) - real,allocatable :: epsedn(:) - real,allocatable :: epsodn(:) - real,allocatable :: kenorm_e(:) - real,allocatable :: kenorm_o(:) - real,allocatable :: gamma_e(:) - real,allocatable :: gamma_o(:) + real(kind_dbl_prec),allocatable :: epse (:) + real(kind_dbl_prec),allocatable :: epso (:) + real(kind_dbl_prec),allocatable :: epsedn(:) + real(kind_dbl_prec),allocatable :: epsodn(:) + real(kind_dbl_prec),allocatable :: kenorm_e(:) + real(kind_dbl_prec),allocatable :: kenorm_o(:) + real(kind_dbl_prec),allocatable :: gamma_e(:) + real(kind_dbl_prec),allocatable :: gamma_o(:) - real,allocatable :: snnp1ev(:) - real,allocatable :: snnp1od(:) + real(kind_dbl_prec),allocatable :: snnp1ev(:) + real(kind_dbl_prec),allocatable :: snnp1od(:) - real,allocatable :: plnev_a(:,:) - real,allocatable :: plnod_a(:,:) - real,allocatable :: plnew_a(:,:) - real,allocatable :: plnow_a(:,:) + real(kind_dbl_prec),allocatable :: plnev_a(:,:) + real(kind_dbl_prec),allocatable :: plnod_a(:,:) + real(kind_dbl_prec),allocatable :: plnew_a(:,:) + real(kind_dbl_prec),allocatable :: plnow_a(:,:) - real,allocatable :: trie_ls(:,:,:) - real,allocatable :: trio_ls(:,:,:) + real(kind_dbl_prec),allocatable :: trie_ls(:,:,:) + real(kind_dbl_prec),allocatable :: trio_ls(:,:,:) integer lotls integer nlunit @@ -73,7 +73,7 @@ module stochy_internal_state_mod integer lan,lat integer nx,ny,nz integer, allocatable :: len(:) - real, allocatable :: parent_lons(:,:),parent_lats(:,:) + real(kind_dbl_prec), allocatable :: parent_lons(:,:),parent_lats(:,:) ! diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index 922e2a4..39fc0fd 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -12,7 +12,7 @@ module stochy_namelist_def public integer, parameter :: max_n_var_lndp = 20 ! 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,nsocnskeb, nsskeb,lon_s,lat_s,ntrunc + integer nsspp,nssppt,nsshum,nsepbl,nsocnsppt,nsocnskeb,nsskeb,nslndp,lon_s,lat_s,ntrunc ! pjp stochastic phyics integer skeb_varspect_opt,skeb_npass @@ -21,13 +21,15 @@ module stochy_namelist_def real(kind=kind_dbl_prec) :: skeb_sigtop1,skeb_sigtop2, & sppt_sigtop1,sppt_sigtop2,shum_sigefold, & skeb_vdof - real(kind=kind_dbl_prec) skeb_diss_smooth,epblint,ocnspptint,ocnskebint,spptint,shumint,skebint,skebnorm + real(kind=kind_dbl_prec) skeb_diss_smooth,epblint,ocnspptint,ocnskebint,sppint,spptint,shumint,skebint,skebnorm + real(kind=kind_phys) lndpint real(kind=kind_dbl_prec), dimension(5) :: skeb,skeb_lscale,skeb_tau real(kind=kind_dbl_prec), dimension(5) :: sppt,sppt_lscale,sppt_tau real(kind=kind_dbl_prec), dimension(5) :: shum,shum_lscale,shum_tau real(kind=kind_dbl_prec), dimension(5) :: epbl,epbl_lscale,epbl_tau real(kind=kind_dbl_prec), dimension(5) :: ocnsppt,ocnsppt_lscale,ocnsppt_tau real(kind=kind_dbl_prec), dimension(5) :: ocnskeb,ocnskeb_lscale,ocnskeb_tau + real(kind=kind_phys), dimension(7) :: pbl_taper integer,dimension(5) ::skeb_vfilt integer(kind=kind_dbl_prec),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_epbl,iseed_ocnsppt,iseed_ocnskeb,iseed_epbl2 logical stochini,sppt_logit,new_lscale @@ -40,14 +42,14 @@ module stochy_namelist_def integer lndp_type integer lndp_model_type 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_phys), dimension(max_n_var_lndp) :: lndp_prt_list - real(kind=kind_dbl_prec), dimension(max_n_var_spp) :: spp_lscale & + real(kind=kind_phys), 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 + character(len=10), dimension(max_n_var_spp) :: spp_var_list + real(kind=kind_phys), dimension(max_n_var_spp) :: spp_prt_list end module stochy_namelist_def diff --git a/stochy_patterngenerator.F90 b/stochy_patterngenerator.F90 index 0f0bb3b..c2ea3e0 100644 --- a/stochy_patterngenerator.F90 +++ b/stochy_patterngenerator.F90 @@ -1,4 +1,4 @@ -!>nbrief The module 'stochy_patterngenerator_mod' contains the derived type random_pattern +!>@brief The module 'stochy_patterngenerator_mod' contains the derived type random_pattern !! which controls the characteristics of the random pattern !! This is a modified version of the original one stochy_patterngenerator.F90, where the !! the random patterns are not properly normalized @@ -150,7 +150,6 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& rpattern(np)%lengthscale = lscale(np) rpattern(np)%dt = delt rpattern(np)%phi = exp(-delt/tscale(np)) - !rpattern(np)%phi = 0 rpattern(np)%stdev = stdev(np) allocate(rpattern(np)%varspectrum(ndimspec)) allocate(rpattern(np)%varspectrum1d(0:ntrunc)) @@ -170,8 +169,8 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& !count4 = iseed(np) + member_id ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. - !count4 = mod(iseed(np) + member_id + 2147483648, 4294967296) - 2147483648 - count4 = mod(iseed(np) + 2147483648, 4294967296) - 2147483648 + !count4 = mod(iseed(np) + member_id + 2147483648_8, 4294967296_8) - 2147483648_8 + count4 = mod(iseed(np) + 2147483648_8, 4294967296_8) - 2147483648_8 print *,'using seed',count4,iseed(np)!,member_id endif endif @@ -180,7 +179,6 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& rpattern(np)%seed = count4 ! set seed (to be the same) on all tasks. Save random state. call random_setseed(rpattern(np)%seed,rpattern(np)%rstate) - !print*,'calling setvarspect',bn_local if (varspect_opt .lt. 0 .or. varspect_opt .gt. 1) then if (is_rootpe()) then print *,'WARNING: illegal value for varspect_opt (should be 0,1,or 2), using 0 (gaussian spectrum)...' @@ -413,14 +411,12 @@ subroutine setvarspect(rpattern,varspect_opt,new_lscale,berner_normalize) if (new_lscale) then !fix for proper lengthscale print*, 'Proper lengthscale condition' -! print*, 'rpattern lap:',rpattern%lap rpattern%varspectrum = exp((rpattern%lengthscale*0.25)**2*rpattern%lap*inv_rerth_sq) do n=0,ntrunc rpattern%varspectrum1d(n) = exp(-(rpattern%lengthscale*0.25)**2*float(n)*(float(n)+1.)*inv_rerth_sq) enddo else print*, 'Not proper lengthscale condition' - !print*, 'rpattern lap:',rpattern%lap rpattern%varspectrum = exp(rpattern%lengthscale**2*rpattern%lap/(4.*rerth**2)) do n=0,ntrunc rpattern%varspectrum1d(n) = exp(-rpattern%lengthscale**2*(float(n)*(float(n)+1.))/(4.*rerth**2)) @@ -437,15 +433,8 @@ subroutine setvarspect(rpattern,varspect_opt,new_lscale,berner_normalize) ! scaling factors for spectral coeffs of white noise pattern with unit variance rpattern%varspectrum(1) = 0 rpattern%varspectrum(2:) = rpattern%degree(2:)**(rpattern%lengthscale) - - !do n=0,ntrunc - ! rpattern%varspectrum1d(n) = float(n)**(rpattern%lengthscale) - !enddo - ! scaling factors for spectral coeffs of white noise pattern with unit variance - !rpattern%varspectrum = sqrt(ntrunc*(rpattern%degree**(rpattern%lengthscale))) - !rpattern%varspectrum = rpattern%degree**(rpattern%lengthscale)! modified following equation 4 in Berner et al 2009 endif - print*, 'Berner Normalization:',berner_normalize + !print*, 'Berner Normalization:',berner_normalize if (.not. berner_normalize) then ! normalize the GEFS way noise = 0. diff --git a/unit_tests/input.nml.ca_template b/unit_tests/input.nml.ca_template index b4cd043..dcd9745 100644 --- a/unit_tests/input.nml.ca_template +++ b/unit_tests/input.nml.ca_template @@ -51,8 +51,8 @@ ca_sgs = CA_SGS ca_global = CA_GLOBAL nca = 1 - scells = 2600 - tlives = 1800 + ncells = 5 + nlives = 12 nseed = 100 nfracseed = 0.5 rcell = 0.72 diff --git a/unit_tests/run_unit_tests_ca.sh b/unit_tests/run_unit_tests_ca.sh index 77c4f04..0a9553d 100755 --- a/unit_tests/run_unit_tests_ca.sh +++ b/unit_tests/run_unit_tests_ca.sh @@ -10,7 +10,7 @@ RES=96 NPX=`expr $RES + 1` NPY=`expr $RES + 1` -DO_CA_SGS=.false. +DO_CA_SGS=.true. DO_CA_GLOBAL=.true. source ./module-setup.sh @@ -19,7 +19,7 @@ module use $( pwd -P ) module load modules.stoch EXEC=standalone_ca.x # compile codes -#sh compile_standalone_ca.hera.intel +sh compile_standalone_ca.hera.intel if [ ! -f $EXEC ];then echo "compilation errors" exit 1 @@ -78,6 +78,37 @@ mkdir -p RESTART echo "unit test 1 successful" fi set OMP_NUM_THREADS=2 + + cp input.nml.ca_template input.nml + ct=1 + while [ $ct -le 6 ];do + rm INPUT/ca_data.tile${ct}.nc + ct=`expr $ct + 1` + done + sed -i -e "s/LOX/1/g" input.nml + sed -i -e "s/LOY/1/g" input.nml + sed -i -e "s/NPX/$NPX/g" input.nml + sed -i -e "s/NPY/$NPY/g" input.nml + sed -i -e "s/RES/$RES/g" input.nml + sed -i -e "s/CA_SGS/${DO_CA_SGS}/g" input.nml + sed -i -e "s/CA_GLOBAL/${DO_CA_GLOBAL}/g" input.nml + sed -i -e "s/WARM_START/.true./g" input.nml + time srun --label -n 6 $EXEC >& stdout.1x1_restart2 + mkdir ca_layout_1x1_restart2 + mv ca_out* ca_layout_1x1_restart2 + ct=1 + while [ $ct -le 6 ];do + diff RESTART/ca_data.tile${ct}.nc RESTART/run1_end_ca_data.tile${ct}.nc + if [ $? -ne 0 ];then + echo "restart test failed" + exit 1 + fi + ct=`expr $ct + 1` + done + if [ $? -eq 0 ];then + echo "unit test 1 successful" + fi + set OMP_NUM_THREADS=2 cp input.nml.ca_template input.nml sed -i -e "s/LOX/1/g" input.nml sed -i -e "s/LOY/4/g" input.nml diff --git a/unit_tests/standalone_ca.F90 b/unit_tests/standalone_ca.F90 index 498a558..5f56063 100644 --- a/unit_tests/standalone_ca.F90 +++ b/unit_tests/standalone_ca.F90 @@ -38,8 +38,8 @@ program standalone_ca_global type(grid_box_type) :: grid_box !---cellular automata control parameters integer :: nca !< number of independent cellular automata -integer :: tlives !< cellular automata lifetime -integer :: scells !< cellular automata finer grid +integer :: nlives !< cellular automata lifetime +integer :: ncells !< cellular automata finer grid integer :: nca_g !< number of independent cellular automata integer :: nlives_g !< cellular automata lifetime integer :: ncells_g !< cellular automata finer grid @@ -65,7 +65,7 @@ program standalone_ca_global real(kind=kind_phys), dimension(:,:), allocatable :: ca1, ca2, ca3 -NAMELIST /gfs_physics_nml/ do_ca, ca_sgs, ca_global, nca, scells, tlives, nseed, & +NAMELIST /gfs_physics_nml/ do_ca, ca_sgs, ca_global, nca, ncells, nlives, nseed, & nfracseed, rcell, ca_trigger, ca_entr, ca_closure, nca_g, & ncells_g, nlives_g, nseed_g, ca_smooth, nspinup, iseed_ca, & nsmooth, ca_amplitude, warm_start @@ -316,7 +316,7 @@ program standalone_ca_global dump_time=50 if (warm_start) then istart=dump_time+1 - call read_ca_restart(Atm(1)%domain,scells,nca,ncells_g,nca_g) + call read_ca_restart(Atm(1)%domain,ncells,nca,ncells_g,nca_g) else istart=1 endif @@ -328,7 +328,7 @@ program standalone_ca_global sst,lmsk,lake,condition,ca_deep,ca_turb,ca_shal, & Atm(1)%domain_for_coupler,nblks, & isc,iec,jsc,jec,Atm(1)%npx,Atm(1)%npy, levs, & - nthresh,rcell,Atm(1)%tile_of_mosaic,nca,scells,tlives,nfracseed, & ! for new random number + nthresh,Atm(1)%tile_of_mosaic,nca,ncells,nlives,nfracseed, & ! for new random number nseed,iseed_ca ,nspinup,ca_trigger,blksz,root_pe,comm) endif if (ca_global) then diff --git a/update_ca.F90 b/update_ca.F90 index 7205f91..caaa051 100644 --- a/update_ca.F90 +++ b/update_ca.F90 @@ -3,20 +3,20 @@ module update_ca !read and write restart routines, to restart fields !on the ncellsxncells CA grid -use kinddef, only: kind_dbl_prec +use kinddef use halo_exchange, only: atmosphere_scalar_field_halo use random_numbers, only: random_01_CB use mpi_wrapper, only: mype,mp_reduce_min,mp_reduce_max use mpp_domains_mod, only: domain2D,mpp_get_global_domain,CENTER, mpp_get_data_domain, mpp_get_compute_domain,mpp_get_ntile_count,& - mpp_define_mosaic,mpp_get_layout,mpp_define_io_domain,mpp_get_io_domain_layout + mpp_define_mosaic,mpp_get_layout,mpp_define_io_domain,mpp_get_io_domain_layout,mpp_get_domain_extents use mpp_mod, only: mpp_error, mpp_pe, mpp_root_pe, & NOTE, FATAL use fms2_io_mod, only: FmsNetcdfDomainFile_t, unlimited, & - open_file, close_file, & - register_axis, register_restart_field, & - register_variable_attribute, register_field, & - read_restart, write_restart, write_data, & - get_global_io_domain_indices, variable_exists + open_file, close_file, & + register_axis, register_restart_field, & + register_variable_attribute, register_field, & + read_restart, write_restart, write_data, & + get_global_io_domain_indices, variable_exists implicit none @@ -34,20 +34,21 @@ module update_ca integer,public :: iscnx_g,iecnx_g,jscnx_g,jecnx_g integer*8, public :: csum type(domain2D),public :: domain_sgs,domain_global +logical, public :: cold_start_ca_sgs=.true.,cold_start_ca_global=.true. contains !Compute CA domain:-------------------------------------------------------------------------- -subroutine define_ca_domain(domain_in,domain_out,ncells,nxncells_local,nyncells_local) +subroutine define_ca_domain(domain_in,domain_out,halo,ncells,nxncells_local,nyncells_local) implicit none type(domain2D),intent(inout) :: domain_in type(domain2D),intent(inout) :: domain_out -integer,intent(in) :: ncells +integer,intent(in) :: ncells,halo integer,intent(out) :: nxncells_local, nyncells_local -integer :: halo1 = 1 integer :: layout(2) +integer, allocatable :: xextent(:,:), yextent(:,:) integer :: ntiles integer, allocatable :: pe_start(:), pe_end(:) @@ -59,7 +60,12 @@ subroutine define_ca_domain(domain_in,domain_out,ncells,nxncells_local,nyncells_ call mpp_get_global_domain(domain_in,xsize=nx,ysize=ny,position=CENTER) call mpp_get_layout(domain_in,layout) ntiles = mpp_get_ntile_count(domain_in) - !write(1000+mpp_pe(),*) "nx,ny: ",nx,ny + allocate(xextent(layout(1),ntiles)) + allocate(yextent(layout(2),ntiles)) + call mpp_get_domain_extents(domain_in,xextent,yextent) + xextent=xextent*ncells + yextent=yextent*ncells + !write(1000+mpp_pe(),*) "nx,ny: ",nx,ny !write(1000+mpp_pe(),*) "layout: ",layout !--- define mosaic for domain_out refined by 'ncells' from domain_in @@ -72,9 +78,12 @@ subroutine define_ca_domain(domain_in,domain_out,ncells,nxncells_local,nyncells_ pe_start(n) = mpp_root_pe() + (n-1)*layout(1)*layout(2) pe_end(n) = mpp_root_pe() + n*layout(1)*layout(2)-1 enddo - call define_cubic_mosaic(domain_out, nxncells_local-1, nyncells_local-1, layout, pe_start, pe_end, halo1 ) + call define_cubic_mosaic(domain_out, nxncells_local-1, nyncells_local-1, & + layout, pe_start, pe_end, halo, ntiles, xextent, yextent ) deallocate(pe_start) deallocate(pe_end) + deallocate(xextent) + deallocate(yextent) end subroutine define_ca_domain !--------------------------------------------------------------------------------------------- @@ -88,7 +97,7 @@ subroutine write_ca_restart(timestamp) character(len=32) :: fn_ca = 'ca_data.nc' type(FmsNetcdfDomainFile_t) :: CA_restart -integer :: id_restart,ncells,nx,ny,i +integer :: id_restart,nx,ny,i integer :: is,ie,js,je,nca,nca_g integer, allocatable, dimension(:) :: buffer @@ -183,40 +192,31 @@ subroutine write_ca_restart(timestamp) end subroutine write_ca_restart -subroutine read_ca_restart(domain_in,scells,nca,ncells_g,nca_g) +subroutine read_ca_restart(domain_in,halo,ncells,nca,ncells_g,nca_g) !Read restart files implicit none type(FmsNetcdfDomainFile_t) :: CA_restart type(domain2D), intent(inout) :: domain_in -integer,intent(in) :: scells,nca,nca_g,ncells_g +integer,intent(in) :: ncells,nca,nca_g,ncells_g,halo character(len=32) :: fn_ca = 'ca_data.nc' character(len=64) :: fname integer :: id_restart integer :: nxc,nyc,i real :: pi,re,dx -integer :: ncells,nx,ny +integer :: nx,ny character(5) :: indir='INPUT' logical :: amiopen integer, allocatable, dimension(:) :: io_layout(:) - - call mpp_get_global_domain(domain_in,xsize=nx,ysize=ny,position=CENTER) -!Set time and length scales: - pi=3.14159 - re=6371000. - dx=0.5*pi*re/real(nx) - ncells=int(dx/real(scells)) - ncells= MIN(ncells,10) - fname = trim(indir)//'/'//trim(fn_ca) if (nca .gt. 0 ) then allocate(io_layout(2)) io_layout=mpp_get_io_domain_layout(domain_in) - call define_ca_domain(domain_in,domain_sgs,ncells,nxncells,nyncells) + call define_ca_domain(domain_in,domain_sgs,halo,ncells,nxncells,nyncells) call mpp_define_io_domain(domain_sgs, io_layout) call mpp_get_compute_domain (domain_sgs,iscnx,iecnx,jscnx,jecnx) amiopen=open_file(CA_restart, trim(fname), 'read', domain=domain_sgs, is_restart=.true., dont_add_res_to_filename=.true.) @@ -242,19 +242,21 @@ subroutine read_ca_restart(domain_in,scells,nca,ncells_g,nca_g) call mpp_error(NOTE,'reading CA_sgs restart data from INPUT/ca_data.tile*.nc') call read_restart(CA_restart) call close_file(CA_restart) + cold_start_ca_sgs=.false. else call mpp_error(NOTE,'No CA_sgs restarts - cold starting CA') + cold_start_ca_sgs=.true. endif endif if (nca_g .gt. 0 ) then - domain_global=domain_in amiopen=open_file(CA_restart, trim(fname), 'read', domain=domain_global, is_restart=.true., dont_add_res_to_filename=.true.) if( amiopen ) then call register_axis(CA_restart, 'xaxis_2', 'X') call register_axis(CA_restart, 'yaxis_2', 'Y') call register_axis(CA_restart, 'nca_g', nca_g) - !call define_ca_domain(domain_in,domain_global,ncells_g,nxncells_g,nyncells_g) + call define_ca_domain(domain_in,domain_global,halo,ncells_g,nxncells_g,nyncells_g) + call mpp_define_io_domain(domain_global, io_layout) call mpp_get_compute_domain (domain_global,iscnx_g,iecnx_g,jscnx_g,jecnx_g) nxc = iecnx_g-iscnx_g+1 nyc = jecnx_g-jscnx_g+1 @@ -271,44 +273,49 @@ subroutine read_ca_restart(domain_in,scells,nca,ncells_g,nca_g) call mpp_error(NOTE,'reading CA_global restart data from INPUT/ca_data.tile*.nc') call read_restart(CA_restart) call close_file(CA_restart) + cold_start_ca_global=.false. else call mpp_error(NOTE,'No CA_global restarts - cold starting CA') + cold_start_ca_global=.true. endif endif end subroutine read_ca_restart -subroutine update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,first_time_step,nca,nxc,nyc,nxch,nych,nlon,& - nlat,isc,iec,jsc,jec, npx,npy, & - CA,ca_plumes,iini,ilives_in,nlives, & +subroutine update_cells_sgs(kstep,halo,dt,initialize_ca,iseed_ca,first_flag,restart,first_time_step,nca,nxc,nyc,nxch,nych,nlon,& + nlat,isc,iec,jsc,jec,ca_advect, npx,npy, & + CA,ca_plumes,iini,ilives_in,uhigh,vhigh,dxhigh,nlives, & nfracseed,nseed,nspinup,nf,nca_plumes,ncells,mytile) +use plumes_mod + implicit none integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca,isc,iec,jsc,jec,npx,npy -integer(kind=kind_dbl_prec), intent(in) :: iseed_ca +integer(8), intent(in) :: iseed_ca integer, intent(in) :: iini(nxc,nyc,nca),initialize_ca,ilives_in(nxc,nyc,nca) -integer, intent(in) :: mytile -real, intent(out) :: CA(nlon,nlat) +integer, intent(in) :: mytile,halo +real(kind_phys), intent(out) :: CA(nlon,nlat) integer, intent(out) :: ca_plumes(nlon,nlat) integer, intent(in) :: nlives,nseed, nspinup, nf,ncells -real, intent(in) :: nfracseed -logical, intent(in) :: nca_plumes,restart,first_flag,first_time_step +real(kind_phys), intent(in) :: nfracseed,dt,dxhigh(nxc,nyc) +real(kind_phys), intent(inout) :: uhigh(nxc,nyc),vhigh(nxc,nyc) +logical, intent(in) :: nca_plumes,restart,ca_advect,first_flag,first_time_step integer, allocatable :: V(:),L(:),B(:) integer, allocatable :: AG(:,:) -integer :: inci, incj, i, j, k,sub,spinup,it,halo,k_in,isize,jsize -integer :: ih, jh,kend, boardmax,livemax +integer :: inci, incj, i, j, k,sub,spinup,it,k_in,isize,jsize +integer :: ih, jh,kend, boardmax,livemax, Xn,Yn real, allocatable :: board_halo(:,:,:) -integer, dimension(nxc,nyc) :: neighbours, birth, thresh +integer, dimension(nxc,nyc) :: neighbours,birth,thresh,adlives,adgrid integer, dimension(nxc,nyc) :: newcell, temp,newseed integer, dimension(ncells,ncells) :: onegrid integer(8) :: nx_full,ny_full integer(8) :: iscale = 10000000000 logical, save :: start_from_restart - -real, dimension(nxc,nyc) :: noise_b +real, dimension(nxch,nych) :: adlives_halo,adgrid_halo +real, dimension(nxc,nyc) :: noise_b,umax,vmax,umin,vmin,dyhigh integer(8) :: count, count_rate, count_max, count_trunc integer :: count4 integer*8 :: i1,j1 @@ -322,21 +329,21 @@ subroutine update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,firs endif !------------------------------------------------------------------------------------------------- -halo=1 isize=nlon+2*halo jsize=nlat+2*halo k_in=1 if (.not. allocated(board))then allocate(board(nxc,nyc,nca)) + board=0.0 endif if (.not. allocated(lives))then allocate(lives(nxc,nyc,nca)) + lives=0.0 endif if(.not. allocated(board_halo))then allocate(board_halo(nxch,nych,1)) endif - !Step 2: Initialize CA, if restart data exist (board,lives > 0) initialize from restart file, otherwise initialize at time- !step initialize_ca. @@ -390,7 +397,7 @@ subroutine update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,firs count_trunc = iscale*(count/iscale) count4 = count - count_trunc + mytile *( i1+nx_full*(j1-1)) ! no need to multply by 7 since time will be different in sgs else - count4 = mod((iseed_ca*nf+mytile)*(i1+nx_full*(j1-1))+ 2147483648, 4294967296) - 2147483648 + count4 = mod((iseed_ca*nf+mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 endif noise_b(i,j)=real(random_01_CB(kstep,count4),kind=8) enddo @@ -413,7 +420,7 @@ subroutine update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,firs neighbours=0 birth=0 newcell=0 - + board_halo=0 !--- copy board into the halo-augmented board_halo board_halo(1+halo:nxc+halo,1+halo:nyc+halo,1) = real(board(1:nxc,1:nyc,1),kind=8) @@ -491,9 +498,75 @@ subroutine update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,firs enddo enddo +enddo !spinup - enddo !spinup +!ADVECTION OF CA CELLS: +!Let the CFL criteria limit the maximum wind speed, such that the +!advection (distance) of a single CA cell does not move outside the +!Halo region +if(ca_advect)then + do j=1,nyc + do i=1,nxc + dyhigh(i,j)=dxhigh(i,j) + umax(i,j)=((dxhigh(i,j))*halo)/dt + vmax(i,j)=((dyhigh(i,j))*halo)/dt + enddo + enddo + + umin(:,:) = -1.0*umax(:,:) + vmin(:,:) = -1.0*vmax(:,:) + + do j=1,nyc + do i=1,nxc + uhigh(i,j)=MIN(uhigh(i,j),umax(i,j)) + vhigh(i,j)=MIN(vhigh(i,j),vmax(i,j)) + uhigh(i,j)=MAX(uhigh(i,j),umin(i,j)) + vhigh(i,j)=MAX(vhigh(i,j),vmin(i,j)) + enddo + enddo + +!Move CA cells in direction of the wind + + adlives_halo(:,:)=0. + adgrid_halo(:,:)=0. + do j=1,nyc + do i=1,nxc + ih=i+halo + jh=j+halo + Xn=ih+nint((uhigh(i,j)/dxhigh(i,j))*dt) + Yn=jh+nint((vhigh(i,j)/dyhigh(i,j))*dt) + adgrid_halo(Xn,Yn)=adgrid_halo(Xn,Yn)+board(i,j,nf) + adlives_halo(Xn,Yn)=adlives_halo(Xn,Yn)+lives(i,j,nf) + enddo + enddo + + call atmosphere_scalar_field_halo (adgrid_halo, halo, nxch, nych, 1, & + iscnx, iecnx, jscnx, jecnx, & + nxncells, nyncells, domain_sgs) + + call atmosphere_scalar_field_halo (adlives_halo, halo, nxch, nych, 1, & + iscnx, iecnx, jscnx, jecnx, & + nxncells, nyncells, domain_sgs) + + + !--- copy the advected fields from the halo-augmented fields + adgrid(1:nxc,1:nyc) = nint(adgrid_halo(1+halo:nxc+halo,1+halo:nyc+halo)) + adlives(1:nxc,1:nyc) = nint(adlives_halo(1+halo:nxc+halo,1+halo:nyc+halo)) + + + do j=1,nyc + do i=1,nxc + lives(i,j,nf)=0. + board(i,j,nf)=0. + lives(i,j,nf)=adlives(i,j) + if(adgrid(i,j)>=1)then + board(i,j,nf)=1 + endif + enddo + enddo + +endif !advection !COARSE-GRAIN BACK TO NWP GRID @@ -572,7 +645,7 @@ end subroutine update_cells_sgs !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,nyc,nxch,nych,nlon,nlat,isc,iec,jsc,jec, & +subroutine update_cells_global(kstep,halo,first_time_step,iseed_ca,restart,nca,nxc,nyc,nxch,nych,nlon,nlat,isc,iec,jsc,jec, & npx,npy,CA,iini_g,ilives_g, & nlives,ncells,nfracseed,nseed,nspinup,nf,mytile) @@ -580,15 +653,15 @@ subroutine update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,ny integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca,isc,iec,jsc,jec,npx,npy integer, intent(in) :: iini_g(nxc,nyc,nca), ilives_g(nxc,nyc) -integer(kind=kind_dbl_prec), intent(in) :: iseed_ca +integer(8), intent(in) :: iseed_ca real, intent(out) :: CA(nlon,nlat) logical, intent(in) :: first_time_step logical, intent(in) :: restart integer, intent(in) :: nlives, ncells, nseed, nspinup, nf real, intent(in) :: nfracseed -integer, intent(in) :: mytile +integer, intent(in) :: mytile, halo integer,allocatable :: V(:),L(:) -integer :: inci, incj, i, j, k ,sub,spinup,it,halo,k_in,isize,jsize +integer :: inci, incj, i, j, k ,sub,spinup,it,k_in,isize,jsize integer :: ih, jh,kend real, allocatable :: board_halo(:,:,:) integer, dimension(nxc,nyc) :: neighbours, birth, thresh @@ -602,7 +675,6 @@ subroutine update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,ny !------------------------------------------------------------------------------------------------- -halo=1 isize=nlon+2*halo jsize=nlat+2*halo k_in=1 @@ -611,7 +683,7 @@ subroutine update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,ny if (.not. allocated(lives_g)) allocate(lives_g(nxc,nyc,nca)) if (.not. allocated(board_halo)) allocate(board_halo(nxch,nych,1)) - if(first_time_step .and. .not. restart)then + if(first_time_step .and. cold_start_ca_global)then do j=1,nyc do i=1,nxc board_g(i,j,nf) = iini_g(i,j,nf) @@ -636,7 +708,7 @@ subroutine update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,ny count_trunc = iscale*(count/iscale) count4 = count - count_trunc + mytile *( i1+nx_full*(j1-1)) ! no need to multply by 7 since time will be different in sgs else - count4 = mod(iseed_ca*nf+(7*mytile)*(i1+nx_full*(j1-1))+ 2147483648, 4294967296) - 2147483648 + count4 = mod(iseed_ca*nf+(7*mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 endif noise_b(i,j)=real(random_01_CB(kstep,count4),kind=8) enddo @@ -652,7 +724,7 @@ subroutine update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,ny enddo endif - if(first_time_step .and. .not. restart)then + if(first_time_step .and. cold_start_ca_global)then spinup=nspinup else spinup = 1 @@ -767,26 +839,29 @@ end subroutine update_cells_global ! and modified to make it simpler to use. ! domain_decomp in fv_mp_mod.F90 does something similar, but it does a ! few other unnecessary things (and requires more arguments). - subroutine define_cubic_mosaic(domain, ni, nj, layout, pe_start, pe_end, halo) + subroutine define_cubic_mosaic(domain, ni, nj, layout, pe_start, pe_end, halo, ntiles, xextent, yextent) type(domain2d), intent(inout) :: domain - integer, intent(in) :: ni, nj + integer, intent(in) :: ni, nj, ntiles, xextent(:,:), yextent(:,:) integer, intent(in) :: layout(:) integer, intent(in) :: pe_start(:), pe_end(:) integer, intent(in) :: halo !--- local variables - integer :: global_indices(4,6), layout2D(2,6) - integer, dimension(12) :: istart1, iend1, jstart1, jend1, tile1 - integer, dimension(12) :: istart2, iend2, jstart2, jend2, tile2 - integer :: ntiles, num_contact + integer,allocatable :: global_indices(:,:), layout2D(:,:) + integer,allocatable :: istart1(:), iend1(:), jstart1(:), jend1(:), tile1(:) + integer,allocatable :: istart2(:), iend2(:), jstart2(:), jend2(:), tile2(:) + integer :: num_contact integer :: i - ntiles = 6 - num_contact = 12 - if(size(pe_start(:)) .NE. 6 .OR. size(pe_end(:)) .NE. 6 ) call mpp_error(FATAL, & - "define_cubic_mosaic: size of pe_start and pe_end should be 6") if(size(layout) .NE. 2) call mpp_error(FATAL, & "define_cubic_mosaic: size of layout should be 2") - do i = 1, 6 + + if(ntiles==6)then + num_contact = 12 + allocate(global_indices(4,ntiles)) + allocate(layout2D(2,ntiles)) + allocate(istart1(num_contact), iend1(num_contact), jstart1(num_contact), jend1(num_contact), tile1(num_contact) ) + allocate(istart2(num_contact), iend2(num_contact), jstart2(num_contact), jend2(num_contact), tile2(num_contact) ) + do i = 1, ntiles layout2D(:,i) = layout(:) global_indices(1,i) = 1 global_indices(2,i) = ni @@ -853,11 +928,40 @@ subroutine define_cubic_mosaic(domain, ni, nj, layout, pe_start, pe_end, halo) istart1(12) = ni; iend1(12) = ni; jstart1(12) = 1; jend1(12) = nj istart2(12) = 1; iend2(12) = 1; jstart2(12) = 1; jend2(12) = nj + else if(ntiles==1) then !Single tile + + num_contact = 0 + + allocate(global_indices(4,ntiles)) + allocate(layout2D(2,ntiles)) + allocate(istart1(num_contact+1), iend1(num_contact+1), jstart1(num_contact+1), jend1(num_contact+1), tile1(num_contact+1) ) + allocate(istart2(num_contact+1), iend2(num_contact+1), jstart2(num_contact+1), jend2(num_contact+1), tile2(num_contact+1) ) + + do i = 1, ntiles + layout2D(:,i) = layout(:) + global_indices(1,i) = 1 + global_indices(2,i) = ni + global_indices(3,i) = 1 + global_indices(4,i) = nj + enddo + + else + + call mpp_error(FATAL, & + "ntiles should be either 6 or 1 to run cellular automata") + + endif !global or regional domain + call mpp_define_mosaic(global_indices, layout2D, domain, ntiles, & num_contact, tile1, tile2, istart1, iend1, jstart1, jend1, & istart2, iend2, jstart2, jend2, pe_start, pe_end, symmetry=.true., & - whalo=halo, ehalo=halo, shalo=halo, nhalo=halo, & - name='CA cubic mosaic') + whalo=halo, ehalo=halo, shalo=halo, nhalo=halo,name='CA cubic mosaic', & + xextent=xextent, yextent=yextent) + + deallocate(global_indices) + deallocate(layout2D) + deallocate(istart1, iend1, jstart1, jend1, tile1) + deallocate(istart2, iend2, jstart2, jend2, tile2) end subroutine define_cubic_mosaic From 88966c71407ee6c10a3a63f1bc3c79b71fbfad80 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Wed, 7 Aug 2024 18:42:48 -0600 Subject: [PATCH 7/8] Int kinds incorporates @mnlevy1982's changes to integer kinds to get this to run with CESM's DEBUG mode --- cellular_automata_global.F90 | 4 ++-- cellular_automata_sgs.F90 | 4 ++-- mersenne_twister.F90 | 2 +- stochy_patterngenerator.F90 | 4 ++-- update_ca.F90 | 8 ++++---- 5 files changed, 11 insertions(+), 11 deletions(-) diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index 469db85..cf6bae8 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -50,7 +50,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp integer,save :: iscnx,iecnx,jscnx,jecnx integer :: nxncells, nyncells integer(8) :: count, count_rate, count_max, count_trunc,nx_full -integer(8) :: iscale = 10000000000 +integer(8) :: iscale = 10000000000_8 integer, allocatable :: iini_g(:,:,:),ilives_g(:,:) integer, allocatable :: io_layout(:) real(kind=kind_phys), allocatable :: field_out(:,:,:), field_smooth(:,:) @@ -167,7 +167,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. - count4 = mod(((iseed_ca+7)*mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 + count4 = int(mod(int(((iseed_ca+7)*mytile)*(i1+nx_full*(j1-1)), 8) + 2147483648_8, 4294967296_8) - 2147483648_8) endif ct=1 do nf=1,nca diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index dbfee82..75b9f3f 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -66,7 +66,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak integer :: blocksz,levs,u200,u850 integer, save :: initialize_ca integer(8) :: count, count_rate, count_max, count_trunc,nx_full -integer(8) :: iscale = 10000000000 +integer(8) :: iscale = 10000000000_8 integer, allocatable :: iini(:,:,:),ilives_in(:,:,:),ca_plumes(:,:),io_layout(:) real(kind=kind_phys), allocatable :: ssti(:,:),lsmski(:,:),lakei(:,:),uwindi(:,:),vwindi(:,:),dxi(:,:),heighti(:,:,:) real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),uhigh(:,:),vhigh(:,:),dxhigh(:,:),conditiongrid(:,:) @@ -296,7 +296,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. - count4 = mod((iseed_ca+mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 + count4 = int(mod(int((iseed_ca+mytile)*(i1+nx_full*(j1-1)), 8) + 2147483648_8, 4294967296_8) - 2147483648_8) endif ct=1 do nf=1,nca diff --git a/mersenne_twister.F90 b/mersenne_twister.F90 index c537681..2ee2627 100644 --- a/mersenne_twister.F90 +++ b/mersenne_twister.F90 @@ -176,7 +176,7 @@ module mersenne_twister integer,parameter:: n=624 integer,parameter:: m=397 integer,parameter:: mata=-1727483681 !< constant vector a - integer,parameter:: umask=-2147483648 !< most significant w-r bits + integer,parameter:: umask=-2147483647-1 !< most significant w-r bits integer,parameter:: lmask =2147483647 !< least significant r bits integer,parameter:: tmaskb=-1658038656 !< tempering parameter integer,parameter:: tmaskc=-272236544 !< tempering parameter diff --git a/stochy_patterngenerator.F90 b/stochy_patterngenerator.F90 index c2ea3e0..a423427 100644 --- a/stochy_patterngenerator.F90 +++ b/stochy_patterngenerator.F90 @@ -63,7 +63,7 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& integer(8), intent(inout) :: iseed(npatterns) integer m,j,l,n,nm,nn,np,indev1,indev2,indod1,indod2 integer(8) count, count_rate, count_max, count_trunc - integer(8) :: iscale = 10000000000 + integer(8) :: iscale = 10000000000_8 integer count4, ierr logical :: bn_local ! berner normalization ! integer member_id @@ -170,7 +170,7 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. !count4 = mod(iseed(np) + member_id + 2147483648_8, 4294967296_8) - 2147483648_8 - count4 = mod(iseed(np) + 2147483648_8, 4294967296_8) - 2147483648_8 + count4 = int(mod(int(iseed(np), 8) + 2147483648_8, 4294967296_8) - 2147483648_8) print *,'using seed',count4,iseed(np)!,member_id endif endif diff --git a/update_ca.F90 b/update_ca.F90 index caaa051..7c600d7 100644 --- a/update_ca.F90 +++ b/update_ca.F90 @@ -312,7 +312,7 @@ subroutine update_cells_sgs(kstep,halo,dt,initialize_ca,iseed_ca,first_flag,rest integer, dimension(nxc,nyc) :: newcell, temp,newseed integer, dimension(ncells,ncells) :: onegrid integer(8) :: nx_full,ny_full -integer(8) :: iscale = 10000000000 +integer(8) :: iscale = 10000000000_8 logical, save :: start_from_restart real, dimension(nxch,nych) :: adlives_halo,adgrid_halo real, dimension(nxc,nyc) :: noise_b,umax,vmax,umin,vmin,dyhigh @@ -397,7 +397,7 @@ subroutine update_cells_sgs(kstep,halo,dt,initialize_ca,iseed_ca,first_flag,rest count_trunc = iscale*(count/iscale) count4 = count - count_trunc + mytile *( i1+nx_full*(j1-1)) ! no need to multply by 7 since time will be different in sgs else - count4 = mod((iseed_ca*nf+mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 + count4 = int(mod(int((iseed_ca*nf+mytile)*(i1+nx_full*(j1-1)), 8) + 2147483648_8, 4294967296_8) - 2147483648_8) endif noise_b(i,j)=real(random_01_CB(kstep,count4),kind=8) enddo @@ -670,7 +670,7 @@ subroutine update_cells_global(kstep,halo,first_time_step,iseed_ca,restart,nca,n integer(8) :: count, count_rate, count_max, count_trunc integer :: count4 integer(8) :: nx_full,ny_full -integer(8) :: iscale = 10000000000 +integer(8) :: iscale = 10000000000_8 integer*8 :: i1,j1 !------------------------------------------------------------------------------------------------- @@ -708,7 +708,7 @@ subroutine update_cells_global(kstep,halo,first_time_step,iseed_ca,restart,nca,n count_trunc = iscale*(count/iscale) count4 = count - count_trunc + mytile *( i1+nx_full*(j1-1)) ! no need to multply by 7 since time will be different in sgs else - count4 = mod(iseed_ca*nf+(7*mytile)*(i1+nx_full*(j1-1))+ 2147483648_8, 4294967296_8) - 2147483648_8 + count4 = int(mod(int(iseed_ca*nf+(7*mytile)*(i1+nx_full*(j1-1)), 8) + 2147483648_8, 4294967296_8) - 2147483648_8) endif noise_b(i,j)=real(random_01_CB(kstep,count4),kind=8) enddo From 359e236122fd053f7ff89f351eeee3a2879e1338 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Tue, 20 Aug 2024 14:16:31 -0600 Subject: [PATCH 8/8] ePBL/SPPT and SKEB use separate grids This commit updates the code so that the random patterns for stochastic ePBL and SPPT are on the T points while the random pattern for SKEB is on the B points (tracer cell corners). This changes the API, so MOM6 code will need to be updated. --- stochastic_physics.F90 | 227 ++++++++++++++++++++++++----------------- stochy_data_mod.F90 | 14 ++- 2 files changed, 143 insertions(+), 98 deletions(-) diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 321becb..1f945ac 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -16,16 +16,16 @@ module stochastic_physics contains !>@brief The subroutine 'init_stochastic_physics' initializes the stochastic -!!pattern genertors +!!pattern generators !>@details It reads the stochastic physics namelist (nam_stoch and nam_sfcperts) -!allocates and polulates the necessary arrays +!allocates and populates the necessary arrays 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, & n_var_spp_in, spp_var_list_out, spp_prt_list_out, spp_stddev_cutoff_out, do_spp_in, & - ak, bk, nthreads, mpiroot, mpicomm, iret) + ak, bk, nthreads, mpiroot, mpicomm, iret) !\callgraph !use stochy_internal_state_moa use stochy_data_mod, only : init_stochdata,gg_lats,gg_lons,nsppt, & @@ -53,7 +53,7 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_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(:) +real(kind=kind_dbl_prec), intent(in) :: ak(:), bk(:) logical, intent(out) :: use_zmtnblck_out integer, intent(out) :: skeb_npass_out character(len=3), optional, dimension(:), intent(out) :: lndp_var_list_out @@ -101,27 +101,27 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in if (do_sppt_in.neqv.do_sppt) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings do_sppt and sppt' - iret = 20 + iret = 20 return else if (do_shum_in.neqv.do_shum) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings do_shum and shum' - iret = 20 + iret = 20 return else if (do_skeb_in.neqv.do_skeb) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings do_skeb and skeb' - iret = 20 + iret = 20 return else if (lndp_type_in /= lndp_type) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings lndp_type in physics and nam_sfcperts' - iret = 20 + iret = 20 return else if (n_var_lndp_in /= n_var_lndp) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings n_var_lndp in physics nml, and lndp_* in nam_sfcperts' - iret = 20 + iret = 20 return else if (n_var_spp_in .ne. n_var_spp) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & @@ -261,77 +261,106 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in end subroutine init_stochastic_physics !!!!!!!!!!!!!!!!!!!! -subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_in,do_sppt_in, & - do_skeb_in, mpiroot, mpicomm, iret) -use stochy_data_mod, only : init_stochdata_ocn,gg_lats,gg_lons,& - rad2deg,INTTYP,wlon,rnlat,gis_stochy_ocn -use spectral_transforms , only : latg,lonf,colrad_a +!>@brief The subroutine 'init_stochastic_physics_ocn' initializes the stochastic +!!pattern generators for ocean stochastics +!>@details It reads the stochastic physics namelist (nam_stoch and nam_sfcperts) +!allocates and populates the necessary arrays +subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nxT, nyT, nz, & + geoLonB, geoLatB, nxB, nyB, & + pert_epbl_in, do_sppt_in, & + do_skeb_in, mpiroot, mpicomm, iret) + +! Use statements +use stochy_data_mod, only : init_stochdata_ocn, gg_lats, gg_lons, & + rad2deg, INTTYP, wlon, rnlat, & + gis_stochy_ocn, gis_stochy_ocn_skeb +use spectral_transforms , only : latg, lonf, colrad_a use stochy_namelist_def -use mersenne_twister, only: random_gauss -use mpi_wrapper, only : mpi_wrapper_initialize,mype,npes,is_rootpe +use mersenne_twister, only : random_gauss +use mpi_wrapper, only : mpi_wrapper_initialize, mype, npes, is_rootpe +! Arguments implicit none -real,intent(in) :: delt -integer,intent(in) :: nx,ny,nz -real,intent(in) :: geoLonT(nx,ny),geoLatT(nx,ny) -logical,intent(in) :: pert_epbl_in,do_sppt_in,do_skeb_in -integer,intent(in) :: mpiroot, mpicomm +real, intent(in) :: delt +integer, intent(in) :: nxT, nyT, nz +real, intent(in) :: geoLonT(nxT, nyT), geoLatT(nxT, nyT) +integer, intent(in) :: nxB, nyB +real, intent(in) :: geoLonB(nxB, nyB), geoLatB(nxB, nyB) +logical, intent(in) :: pert_epbl_in, do_sppt_in, do_skeb_in +integer, intent(in) :: mpiroot, mpicomm integer, intent(out) :: iret -real(kind=kind_dbl_prec), parameter :: con_pi =4.0d0*atan(1.0d0) -real :: dx -integer :: k,latghf,km -type(MPI_Comm) :: mpicomm_t ! FIXME once MOM6 updates to use mpi_f90 types -rad2deg=180.0/con_pi -mpicomm_t%mpi_val = mpicomm -call mpi_wrapper_initialize(mpiroot,mpicomm_t) -gis_stochy_ocn%nodes = npes -gis_stochy_ocn%mype = mype -gis_stochy_ocn%nx=nx -gis_stochy_ocn%ny=ny -allocate(gis_stochy_ocn%len(ny)) -allocate(gis_stochy_ocn%parent_lons(nx,ny)) -allocate(gis_stochy_ocn%parent_lats(nx,ny)) -gis_stochy_ocn%len(:)=nx -gis_stochy_ocn%parent_lons=geoLonT -gis_stochy_ocn%parent_lats=geoLatT +! Local variables +real(kind=kind_dbl_prec), parameter :: con_pi = 4.0d0 * atan(1.0d0) +real :: dx +integer :: k, latghf, km +type(MPI_Comm) :: mpicomm_t ! FIXME once MOM6 updates to use mpi_f90 types -INTTYP=0 ! bilinear interpolation -km=nz -call init_stochdata_ocn(km,delt,iret) + +rad2deg = 180.0 / con_pi +mpicomm_t%mpi_val = mpicomm +call mpi_wrapper_initialize(mpiroot, mpicomm_t) + +! gis_stochy_ocn is used for ePBL and SPPT +gis_stochy_ocn%nodes = npes +gis_stochy_ocn%mype = mype +gis_stochy_ocn%nx = nxT +gis_stochy_ocn%ny = nyT +allocate(gis_stochy_ocn%len(nyT)) +allocate(gis_stochy_ocn%parent_lons(nxT, nyT)) +allocate(gis_stochy_ocn%parent_lats(nxT, nyT)) +gis_stochy_ocn%len(:) = nxT +gis_stochy_ocn%parent_lons = geoLonT +gis_stochy_ocn%parent_lats = geoLatT + +! gis_stochy_ocn_skeb is used for SKEB +gis_stochy_ocn_skeb%nodes = npes +gis_stochy_ocn_skeb%mype = mype +gis_stochy_ocn_skeb%nx = nxB +gis_stochy_ocn_skeb%ny = nyB +allocate(gis_stochy_ocn_skeb%len(nyB)) +allocate(gis_stochy_ocn_skeb%parent_lons(nxB, nyB)) +allocate(gis_stochy_ocn_skeb%parent_lats(nxB, nyB)) +gis_stochy_ocn_skeb%len(:) = nxB +gis_stochy_ocn_skeb%parent_lons = geoLonB +gis_stochy_ocn_skeb%parent_lats = geoLatB + +INTTYP = 0 ! bilinear interpolation +km = nz +call init_stochdata_ocn(km, delt, iret) if (do_sppt_in.neqv.do_ocnsppt) then - write(0,'(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & - & ' namelist settings do_sppt and sppt' - iret = 20 + write(0, '(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & + & ' namelist settings do_sppt and sppt' + iret = 20 return else if (pert_epbl_in.neqv.pert_epbl) then - write(0,'(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & - & ' namelist settings pert_epbl and epbl' - iret = 20 + write(0, '(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & + & ' namelist settings pert_epbl and epbl' + iret = 20 return else if (do_skeb_in.neqv.do_ocnskeb) then - write(0,'(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & - & ' namelist settings do_skeb and skeb' - iret = 20 + write(0, '(*(a))') 'Logic error in stochastic_physics_ocn_init: incompatible', & + & ' namelist settings do_skeb and skeb' + iret = 20 return end if ! get interpolation weights ! define gaussian grid lats and lons -latghf=latg/2 +latghf = latg / 2 allocate(gg_lats(latg)) allocate(gg_lons(lonf)) -do k=1,latghf - gg_lats(k)=-1.0*colrad_a(latghf-k+1)*rad2deg - gg_lats(latg-k+1)=-1*gg_lats(k) -enddo -dx=360.0/lonf -do k=1,lonf - gg_lons(k)=dx*(k-1) -enddo -WLON=gg_lons(1)-(gg_lons(2)-gg_lons(1)) -RNLAT=gg_lats(1)*2-gg_lats(2) -print*,'finished ocean init' +do k = 1, latghf + gg_lats(k) = -1.0 * colrad_a(latghf - k + 1) * rad2deg + gg_lats(latg - k + 1) = -1 * gg_lats(k) +end do +dx = 360.0 / lonf +do k = 1, lonf + gg_lons(k) = dx * (k - 1) +end do +WLON = gg_lons(1) - (gg_lons(2)-gg_lons(1)) +RNLAT = gg_lats(1)*2 - gg_lats(2) +print*, 'finished ocean init' end subroutine init_stochastic_physics_ocn !!!!!!!!!!!!!!!!!!!! @@ -340,9 +369,9 @@ end subroutine init_stochastic_physics_ocn !>@brief The subroutine 'run_stochastic_physics' updates the random patterns if !!necessary !>@details It updates the AR(1) in spectral space -!allocates and polulates the necessary arrays +!allocates and populates the necessary arrays -subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, skebu_wts, & +subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, skebu_wts, & skebv_wts, sfc_wts, spp_wts, nthreads) !\callgraph @@ -350,9 +379,9 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s 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, & rpattern_spp, nspp, vfact_spp -use get_stochy_pattern_mod,only : get_random_pattern_scalar,get_random_pattern_vector, & +use get_stochy_pattern_mod,only : get_random_pattern_scalar,get_random_pattern_vector, & get_random_pattern_sfc,get_random_pattern_spp -use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,nsspp,nslndp,sppt_logit, & +use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,nssppt,nsshum,nsskeb,nsspp,nslndp,sppt_logit, & lndp_type, n_var_lndp, n_var_spp, do_spp, spp_stddev_cutoff, spp_prt_list use mpi_wrapper, only: is_rootpe implicit none @@ -386,7 +415,7 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s if ( (lndp_type==1) .and. (kdt==0) ) then ! old land pert scheme called once at start - write(0,*) 'calling get_random_pattern_sfc' + write(0,*) 'calling get_random_pattern_sfc' allocate(tmpl_wts(nblks,maxlen,n_var_lndp)) call get_random_pattern_sfc(rpattern_sfc,nlndp,gis_stochy,tmpl_wts) DO blk=1,nblks @@ -439,7 +468,7 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s ENDDO endif endif -if ( lndp_type .EQ. 2 ) then +if ( lndp_type .EQ. 2 ) then ! add time check? if (mod(kdt,nslndp) == 1 .or. nslndp == 1) then allocate(tmpl_wts(gis_stochy%nx,gis_stochy%ny,n_var_lndp)) @@ -480,45 +509,55 @@ subroutine run_stochastic_physics(levs, kdt, fhour, blksz, sppt_wts, shum_wts, s end subroutine run_stochastic_physics -subroutine run_stochastic_physics_ocn(sppt_wts,skeb_wts,t_rp1,t_rp2) +!>@brief The subroutine 'run_stochastic_physics_ocn' updates the random patterns +!!if necessary +!>@details It updates the AR(1) in spectral space +!allocates and populates the necessary arrays +subroutine run_stochastic_physics_ocn(sppt_wts, skeb_wts, t_rp1, t_rp2) use stochy_internal_state_mod -use stochy_data_mod, only : nepbl,nocnsppt,nocnskeb,rpattern_epbl1,rpattern_epbl2,rpattern_ocnsppt, rpattern_ocnskeb, gis_stochy_ocn -use get_stochy_pattern_mod,only : get_random_pattern_scalar +use stochy_data_mod, only : nepbl, nocnsppt, nocnskeb, rpattern_epbl1, & + rpattern_epbl2, rpattern_ocnsppt, & + rpattern_ocnskeb, gis_stochy_ocn, & + gis_stochy_ocn_skeb +use get_stochy_pattern_mod, only : get_random_pattern_scalar use stochy_namelist_def implicit none -real, intent(inout) :: sppt_wts(:,:),t_rp1(:,:),t_rp2(:,:),skeb_wts(:,:) +real, intent(inout) :: sppt_wts(:,:), t_rp1(:,:), t_rp2(:,:), skeb_wts(:,:) real(kind_dbl_prec), allocatable :: tmp_wts(:,:) -if (pert_epbl .OR. do_ocnsppt .OR. do_ocnskeb ) then - allocate(tmp_wts(gis_stochy_ocn%nx,gis_stochy_ocn%ny)) + +if (pert_epbl .OR. do_ocnsppt) then + allocate(tmp_wts(gis_stochy_ocn%nx, gis_stochy_ocn%ny)) if (pert_epbl) then - call get_random_pattern_scalar(rpattern_epbl1,nepbl,gis_stochy_ocn,tmp_wts) - t_rp1(:,:)=2.0/(1+exp(-1*tmp_wts)) - call get_random_pattern_scalar(rpattern_epbl2,nepbl,gis_stochy_ocn,tmp_wts) - t_rp2(:,:)=2.0/(1+exp(-1*tmp_wts)) + call get_random_pattern_scalar(rpattern_epbl1, nepbl, gis_stochy_ocn, tmp_wts) + t_rp1(:,:) = 2.0 / (1.0 + exp(-1 * tmp_wts)) + call get_random_pattern_scalar(rpattern_epbl2, nepbl, gis_stochy_ocn, tmp_wts) + t_rp2(:,:) = 2.0 / (1.0 + exp(-1 * tmp_wts)) else - t_rp1(:,:)=1.0 - t_rp2(:,:)=1.0 + t_rp1(:,:) = 1.0 + t_rp2(:,:) = 1.0 endif if (do_ocnsppt) then - call get_random_pattern_scalar(rpattern_ocnsppt,nocnsppt,gis_stochy_ocn,tmp_wts) - sppt_wts=2.0/(1+exp(-1*tmp_wts)) + call get_random_pattern_scalar(rpattern_ocnsppt, nocnsppt, gis_stochy_ocn, tmp_wts) + sppt_wts = 2.0 / (1.0 + exp(-1 * tmp_wts)) else - sppt_wts=1.0 - endif - if (do_ocnskeb) then - call get_random_pattern_scalar(rpattern_ocnskeb,nocnskeb,gis_stochy_ocn,skeb_wts,normalize=.true.) - else - skeb_wts=1.0 + sppt_wts = 1.0 endif deallocate(tmp_wts) else - sppt_wts(:,:)=1.0 - skeb_wts(:,:)=1.0 - t_rp1(:,:)=1.0 - t_rp2(:,:)=1.0 + sppt_wts(:,:) = 1.0 + t_rp1(:,:) = 1.0 + t_rp2(:,:) = 1.0 +endif + +if (do_ocnskeb) then + call get_random_pattern_scalar(rpattern_ocnskeb, nocnskeb, gis_stochy_ocn_skeb, skeb_wts, normalize=.true.) +else + skeb_wts(:,:) = 1.0 endif 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, & @@ -530,7 +569,7 @@ subroutine finalize_stochastic_physics() if (allocated(gg_lats)) deallocate (gg_lats) if (allocated(gg_lons)) deallocate (gg_lons) if (allocated(sl)) deallocate (sl) - if (nsppt > 0) then + if (nsppt > 0) then if (allocated(rpattern_sppt)) deallocate(rpattern_sppt) if (allocated(vfact_sppt)) deallocate(vfact_sppt) endif @@ -549,7 +588,7 @@ subroutine finalize_stochastic_physics() if (nlndp > 0) then if (allocated(rpattern_sfc)) deallocate(rpattern_sfc) endif - if (nspp > 0) then + if (nspp > 0) then if (allocated(rpattern_spp)) deallocate(rpattern_spp) if (allocated(vfact_spp)) deallocate(vfact_spp) endif diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index 95b8518..aa46eeb 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -43,7 +43,7 @@ module stochy_data_mod real(kind=kind_dbl_prec),public :: wlon,rnlat,rad2deg real(kind=kind_dbl_prec),public, allocatable :: skebu_save(:,:,:),skebv_save(:,:,:) integer,public :: INTTYP - type(stochy_internal_state),public :: gis_stochy,gis_stochy_ocn + type(stochy_internal_state),public :: gis_stochy,gis_stochy_ocn,gis_stochy_ocn_skeb contains !>@brief The subroutine 'init_stochdata' determins which stochastic physics @@ -502,8 +502,14 @@ subroutine init_stochdata_ocn(nlevs,delt,iret) iret=0 call compns_stochy_ocn (delt,iret) if(is_rootpe()) print*,'in init stochdata_ocn' - if ( (.NOT. pert_epbl) .AND. (.NOT. do_ocnsppt) .AND. (.NOT. do_ocnskeb) ) return - call initialize_spectral(gis_stochy_ocn) +! if ( (.NOT. pert_epbl) .AND. (.NOT. do_ocnsppt) .AND. (.NOT. do_ocnskeb) ) return +! call initialize_spectral(gis_stochy_ocn) + if ( pert_epbl .OR. do_ocnsppt .OR. do_ocnskeb ) then + if ( pert_epbl .OR. do_ocnsppt ) call initialize_spectral(gis_stochy_ocn) + if ( do_ocnskeb ) call initialize_spectral(gis_stochy_ocn_skeb) + else + return + endif if (iret/=0) return allocate(noise_e(len_trie_ls,2),noise_o(len_trio_ls,2)) ! determine number of random patterns to be used for each scheme. @@ -723,7 +729,7 @@ subroutine init_stochdata_ocn(nlevs,delt,iret) endif if (is_rootpe()) print *, 'Initialize random pattern for ocnskeb' call patterngenerator_init(ocnskeb_lscale(1:nocnskeb),ocnskebint,ocnskeb_tau(1:nocnskeb),ocnskeb(1:nocnskeb),iseed_ocnskeb,rpattern_ocnskeb, & - lonf,latg,jcap,gis_stochy_ocn%ls_node,nocnskeb,1,0,new_lscale,.true.) + lonf,latg,jcap,gis_stochy_ocn_skeb%ls_node,nocnskeb,1,0,new_lscale,.true.) do n=1,nocnskeb if (stochini) then