2D advection of cellular automata

.travis.yml

@@ -17,10 +17,7 @@ branches:
   only:
     - ufs_public_release
     - release/public-v1
-    - release/public-v2
-    - release/public-v2
     - release/ufs-v1.1.0
-    - master
 
 # Environment variables
 env:

LICENSE.md

@@ -3,7 +3,7 @@ Copyright 2020 The Regents of the University of Colorado
 The stochastic_physics code incorporated in the Unified Forecast System (UFS) 
 was jointly developed by the National Oceanic and Atmospheric Administration and the 
 University of Colorado, Boulder. The gold standard copy of the Code 
-will be maintained by NOAA at https://github.com/noaa-psl/stochastic_physics.
+will be maintained by NOAA at https://github.com/noaa-psd/stochastic_physics.
 
 In cooperation with the Copyright Holder, the National Oceanic and 
 Atmospheric Administration is releasing this code under the 

cellular_automata_global.F90

@@ -16,7 +16,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp
 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
@@ -50,6 +50,7 @@ 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(:,:)
+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(:,:,:)
@@ -73,7 +74,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 +103,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
@@ -203,7 +208,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)
 

cellular_automata_sgs.F90

@@ -6,10 +6,10 @@ module cellular_automata_sgs_mod
 
 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, mytile, &
-            nca,ncells,nlives,nfracseed,nseed,iseed_ca, &
+            nca,ncells,nlives,nfracseed,nseed,iseed_ca,ca_advect, &
             nspinup,ca_trigger,blocksize,mpiroot,mpicomm)
 
 use kinddef,           only: kind_phys,kind_dbl_prec
@@ -37,16 +37,19 @@ 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,ncells,nca,nlives,nseed,nspinup,mpiroot,mpicomm,mytile
 integer(kind=kind_dbl_prec),           intent(in)    :: iseed_ca
 real(kind=kind_phys), intent(in)    :: nfracseed,dtf,nthresh
-logical,intent(in) :: restart,ca_trigger,first_time_step
+logical,intent(in) :: restart,ca_trigger,first_time_step,ca_advect
 integer, intent(in) :: nblks,isc,iec,jsc,jec,npx,npy,nlev,blocksize
-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(:,:)
@@ -58,21 +61,22 @@ 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 :: 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
+real                 :: dz,invgrav
 
 !nca         :: switch for number of cellular automata to be used.
 !            :: for the moment only 1 CA can be used 
@@ -88,8 +92,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.
 
@@ -132,15 +141,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
@@ -152,9 +158,18 @@ 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))
@@ -163,6 +178,9 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
  !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).
@@ -172,6 +190,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
@@ -180,13 +199,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
@@ -277,9 +322,9 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
 
 !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
@@ -328,6 +373,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)

docs/source/getting_started.rst

@@ -5,3 +5,6 @@ The stochastic physics currently only works with the UFS-atmosphere model
 You should get the full system at https://github.com/ufs-community/ufs-weather-model, which will include the stochastic physics code.
 
 In order to enable stochastic physics in a model run, you will need to turn it on via `namelist options <namelist_options.html>`_
+
+If using the CIME workflow decribed at https://ufs-mrweather-app.readthedocs.io/en/latest/, please add do_sppt=T, etc. to user_nl_ufsatm in the case directory.
+

docs/source/index.rst

@@ -14,5 +14,5 @@ Welcome to UFS stochastic physics's documentation!
 Source Code Documentation
 ==========================
 
-#Doxygen generated (out of date) `content <https://noaa-psd.github.io/stochastic_physics/doxygen/index.html>`_
+Doxygen generated `content <https://noaa-psd.github.io/stochastic_physics/doxygen/index.html>`_
 

docs/source/namelist_options.rst

@@ -8,11 +8,11 @@ General options
    :header: "Option", "Description"
    :widths: 30, 50
 
-   "NEW_LSCALE", "Recommended, set to true if use the correct calculation for decorrelation length scale."
    "NTRUNC", "Optional, Spectral resolution (e.g. T126) of random patterns, default is for model to determine proper truncation"
    "LAT_S", "Optional, number of latitude points for the gaussian grid  (must be even), default is for model to determine gaussian grid"
    "LON_S", "Optional, number of longitude points for the gaussian grid (recommend 2xLAT_S, default is for model to determine gaussian grid"
-   "STOCHINI", "Optional, set to true if wanting to read in a previous random pattern"
+   "FHSTOCH", "Optional, forecast hour to write out random pattern in order to restart the pattern for a different forecast (used in DA), file is stoch_out.F<HHH>"
+   "STOCHINI", "Optional, set to true if wanting to read in a previous random pattern (input file needs to be named stoch_ini)."
 
 SPPT options 
 """"""""""""
@@ -21,31 +21,31 @@ SPPT options
    :header: "Option", "Description"
    :widths: 30, 50
 
-   "SPPT", "Amplitudes of random patterns."
-   "SPPT_TAU", "Decorrelation timescales in seconds."
-   "SPPT_LSCALE", "Decorrelation spatial scales in meters."
+   "DO_SPPT", "logical to tell parent atmospheric model to use SPPT"
+   "SPPT", "Amplitudes of random patterns (0.8,0.4,0.2,0.08,0.04) *"
+   "SPPT_TAU", "Decorrelation timescales in seconds (21600,1.728E5,6.912E5,7.776E6,3.1536E7) *"
+   "SPPT_LSCALE", "Decorrelation spatial scales in meters  (250.E3,1000.E3,2000.E3,2000.E3,2000.E3) *"
    "SPPT_LOGIT", "Should be true to limit the SPPT perturbations between 0 and 2.  Otherwise model will crash."
    "ISEED_SPPT", "Seeds for setting the random number sequence (ignored if stochini is true)"
    "SPPT_SIGTOP1", "lower sigma level to taper perturbations to zero (default is 0.1)"
-   "SPPT_SIGTOP2", "upper sigma level to taper perturbations to zero (default is 0.025)"
+   "SPPT_SIGTOP2", "upper sigma level to taper perturbations to zero (0.025)"
    "SPPT_SFCLIMIT", ".T.=tapers the SPPT perturbations to zero at model’s lowest level (helps reduce model crashes)"
    "SPPTINT", "Optional, interval in seconds to update random pattern.  Perturbations still get applied every time-step"
    "USE_ZMTNBLCK", ".T.=do not apply perturbations below the dividing streamline that is diagnosed by the gravity wave drag, mountain blocking scheme"
-   "PERT_MP", ".T.=apply SPPT perturbations to all microphysics specis. If .F. the SPPT is only applied to u,v,t,qv"
-   "PERT_RADTEND", ".T.=apply SPPT perturbations to cloudy sky radiation tendencies. If .F. then do not perturb any radiative tendencies"
-   "PERT_CLDS", ".T.=apply SPPT perturbations to fraction (only works for RRTMG radiation),  if using this option set PERT_RADTEND=.F."
 
+``*``  **SPPT** uses 5 different patterns of varying time/length scales that are added together before being passed to physics
 
 SHUM options 
 """"""""""""
 
 .. csv-table::
    :header: "Option", "Description"
-   :widths: 30, 50
+   :widths: 20, 50
 
-   "SHUM", "Amplitudes of random pattern."
-   "SHUM_TAU", "Decorrelation timescales in seconds."
-   "SHUM_LSCALE", "ecorrelation spatial scales in meters."
+   "DO_SHUM", "logical to tell parent atmospheric model to use SHUM"
+   "SHUM", "Amplitudes of random patterns (0.004)"
+   "SHUM_TAU", "Decorrelation timescales in seconds (21600)"
+   "SHUM_LSCALE", "ecorrelation spatial scales in meters (250000)"
    "SHUM_SIGEFOLD", "e-folding lengthscale (in units of sigma) of specific humidity perturbations, default is 0.2)"
    "SHUMINT", "Optional, interval in seconds to update random pattern.  Perturbations still get applied every time-step"
    "ISEED_SHUM", "Seeds for setting the random number sequence (ignored if stochini is true)."
@@ -57,8 +57,9 @@ SKEB options
    :header: "Option", "Description"
    :widths: 30, 50
 
-   "SKEB", "Amplitudes of random patterns."
-   "SKEB_TAU", "Decorrelation timescales in seconds"
+   "DO_SKEB", "logical to tell parent atmospheric model to use SKEB"
+   "SKEB", "Amplitudes of random patterns (0.5)"
+   "SKEB_TAU", "Decorrelation timescales in seconds (21600)"
    "SKEB_LSCALE", "Decorrelation spatial scales in meters  (250)"
    "ISEED_SKEB", "Seeds for setting the random number sequence (ignored if stochini is true)."
    "SKEBNORM", "0-random pattern is stream function, 1-pattern is K.E. norm, 2-pattern is vorticity (default is 0)"
@@ -69,34 +70,3 @@ SKEB options
    "SKEB_SIGTOP2", "upper sigma level to taper perturbations to zero (0.025)"
    "SKEBINT", "Optional, interval in seconds to update random pattern.  Perturbations still get applied every time-step"
 
-SPP options
-""""""""""""
-
-.. csv-table::
-   :header: "Option", "Description"
-   :widths: 30, 50
-
-   "SPP_VAR_LIST", "List of parameterizations to be perturbed. Check compns_stochy.F90 for options."
-   "SPP_PRT_LIST", "SPP perturbation magnitudes used in each parameterization."
-   "SPP_TAU", "Decorrelation timescales in seconds."
-   "SPP_LSCALE", "Decorrelation spatial scales in meters."
-   "ISEED_SPP", "Seeds for setting the random number sequence (ignored if stochini is true)."
-   "SPP_SIGTOP1", "Lower sigma level to taper perturbations to zero (default is 0.1)"
-   "SPP_SIGTOP2", "Upper sigma level to taper perturbations to zero (0.025)"
-   "SPP_STDDEV_CUTOFF", "Limit for possible perturbation values in standard deviations from the mean."
-
-
-Land perturbation options
-""""""""""""
-
-.. csv-table::
-   :header: "Option", "Description"
-   :widths: 30, 50
-
-   "LNDP_TYPE", "0, no perturbations. 1, old scheme (Gehne et al. 2019); 2, new scheme (Draper 2021)"
-   "LNDP_VAR_LIST", "List of land perturbations parameters. Check compns_stochy.F90 for options"
-   "LNDP_PRT_LIST", "Perturbation magnitudes used for each parameter perturbations."
-   "LNDP_TAU", "Decorrelation timescales in seconds."
-   "LNDP_LSCALE", "Decorrelation spatial scales in meters."
-   "ISEED_LNDP", "Seeds for setting the random number sequence (ignored if stochini is true)."
-