Update to cellular automaton stochastic convective organization scheme

atmos_model.F90

@@ -372,6 +372,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
 #endif
   use fv_mp_mod, only: commglobal
   use mpp_mod, only: mpp_npes
+  use update_ca, only: read_ca_restart
 
   type (atmos_data_type), intent(inout) :: Atmos
   type (time_type), intent(in) :: Time_init, Time, Time_step
@@ -572,7 +573,9 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
                               GFS_data%Coupling, GFS_data%Grid, GFS_data%Tbd, GFS_data%Cldprop, GFS_data%Radtend, &
                               GFS_data%IntDiag, Init_parm, GFS_Diag)
    call FV3GFS_restart_read (GFS_data, GFS_restart_var, Atm_block, GFS_control, Atmos%domain, Atm(mygrid)%flagstruct%warm_start)
-
+   if(GFS_control%ca_sgs)then
+      call read_ca_restart (Atmos%domain,GFS_control%scells)
+   endif
    ! Populate the GFS_data%Statein container with the prognostic state
    ! in Atm_block, which contains the initial conditions/restart data.
    call atmos_phys_driver_statein (GFS_data, Atm_block, flip_vc)
@@ -858,6 +861,7 @@ end subroutine update_atmos_model_state
 
 subroutine atmos_model_end (Atmos)
   use get_stochy_pattern_mod, only: write_stoch_restart_atm
+  use update_ca, only: write_ca_restart
   type (atmos_data_type), intent(inout) :: Atmos
 !---local variables
   integer :: idx, seconds, ierr
@@ -871,6 +875,9 @@ subroutine atmos_model_end (Atmos)
       call FV3GFS_restart_write (GFS_data, GFS_restart_var, Atm_block, &
                                  GFS_control, Atmos%domain)
       call write_stoch_restart_atm('RESTART/atm_stoch.res.nc')
+      if(GFS_control%ca_sgs)then
+         call write_ca_restart(Atmos%domain,GFS_control%scells)
+      endif
     endif
     call stochastic_physics_wrapper_end(GFS_control)
 
@@ -889,13 +896,16 @@ end subroutine atmos_model_end
 !  Write out restart files registered through register_restart_file
 ! </DESCRIPTION>
 subroutine atmos_model_restart(Atmos, timestamp)
+  use update_ca, only: write_ca_restart
   type (atmos_data_type),   intent(inout) :: Atmos
   character(len=*),  intent(in)           :: timestamp
 
     call atmosphere_restart(timestamp)
     call FV3GFS_restart_write (GFS_data, GFS_restart_var, Atm_block, &
                                GFS_control, Atmos%domain, timestamp)
-
+    if(GFS_control%ca_sgs)then
+       call write_ca_restart(Atmos%domain,GFS_control%scells,timestamp)
+    endif
 end subroutine atmos_model_restart
 ! </SUBROUTINE>
 

ccpp/data/GFS_typedefs.F90

@@ -1071,8 +1071,8 @@ module GFS_typedefs
 
 !---cellular automata control parameters
     integer              :: nca             !< number of independent cellular automata
-    integer              :: nlives          !< cellular automata lifetime
-    integer              :: ncells          !< cellular automata finer grid
+    integer              :: tlives          !< cellular automata lifetime
+    integer              :: scells          !< 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
@@ -1085,7 +1085,8 @@ module GFS_typedefs
     logical              :: ca_smooth       !< switch for gaussian spatial filter
     integer              :: iseed_ca        !< seed for random number generation in ca scheme
     integer              :: nspinup         !< number of iterations to spin up the ca
-    real(kind=kind_phys) :: nthresh         !< threshold used for perturbed vertical velocity
+    real(kind=kind_phys) :: rcell           !< threshold used for CA scheme
+    real(kind=kind_phys) :: nthresh         !< threshold used for convection coupling
     real                 :: ca_amplitude    !< amplitude of ca trigger perturbation
     integer              :: nsmooth         !< number of passes through smoother
     logical              :: ca_closure      !< logical switch for ca on closure
@@ -3417,22 +3418,22 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
 
 !---Cellular automaton options
     integer              :: nca            = 1
-    integer              :: ncells         = 5
-    integer              :: nlives         = 30
+    integer              :: scells         = 2600
+    integer              :: tlives         = 1800
     integer              :: nca_g          = 1
     integer              :: ncells_g       = 1
     integer              :: nlives_g       = 100
     real(kind=kind_phys) :: nfracseed      = 0.5
-    integer              :: nseed          = 100000
+    integer              :: nseed          = 1
     integer              :: nseed_g        = 100
-    integer              :: iseed_ca       = 0
+    integer              :: iseed_ca       = 1
     integer              :: nspinup        = 1
     logical              :: do_ca          = .false.
     logical              :: ca_sgs         = .false.
     logical              :: ca_global      = .false.
     logical              :: ca_smooth      = .false.
-    real(kind=kind_phys) :: nthresh        = 0.0
-    real                 :: ca_amplitude   = 500.
+    real(kind=kind_phys) :: rcell          = 0.72
+    real                 :: ca_amplitude   = 0.35
     integer              :: nsmooth        = 100
     logical              :: ca_closure     = .false.
     logical              :: ca_entr        = .false.
@@ -3565,8 +3566,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                           !--- canopy heat storage parameterization
                                z0fac, e0fac,                                                &
                           !--- cellular automata
-                               nca, ncells, nlives, nca_g, ncells_g, nlives_g, nfracseed,   &
-                               nseed, nseed_g, nthresh, do_ca,                              &
+                               nca, scells, tlives, nca_g, ncells_g, nlives_g, nfracseed,   &
+                               nseed, nseed_g, rcell, do_ca,                              &
                                ca_sgs, ca_global,iseed_ca,ca_smooth,                        &
                                nspinup,ca_amplitude,nsmooth,ca_closure,ca_entr,ca_trigger,  &
                           !--- IAU
@@ -4273,8 +4274,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
 
     !--- cellular automata options
     Model%nca              = nca
-    Model%ncells           = ncells
-    Model%nlives           = nlives
+    Model%scells           = scells
+    Model%tlives           = tlives
     Model%nca_g            = nca_g
     Model%ncells_g         = ncells_g
     Model%nlives_g         = nlives_g
@@ -4287,7 +4288,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%iseed_ca         = iseed_ca
     Model%ca_smooth        = ca_smooth
     Model%nspinup          = nspinup
-    Model%nthresh          = nthresh
+    Model%rcell            = rcell
     Model%ca_amplitude     = ca_amplitude
     Model%nsmooth          = nsmooth
     Model%ca_closure       = ca_closure
@@ -5314,8 +5315,8 @@ subroutine control_print(Model)
       print *, ' '
       print *, 'cellular automata'
       print *, ' nca               : ', Model%nca
-      print *, ' ncells            : ', Model%ncells
-      print *, ' nlives            : ', Model%nlives
+      print *, ' scells            : ', Model%scells
+      print *, ' tlives            : ', Model%tlives
       print *, ' nca_g             : ', Model%nca_g
       print *, ' ncells_g          : ', Model%ncells_g
       print *, ' nlives_g          : ', Model%nlives_g
@@ -5328,7 +5329,7 @@ subroutine control_print(Model)
       print *, ' iseed_ca          : ', Model%iseed_ca
       print *, ' ca_smooth         : ', Model%ca_smooth
       print *, ' nspinup           : ', Model%nspinup
-      print *, ' nthresh           : ', Model%nthresh
+      print *, ' rcell             : ', Model%rcell
       print *, ' ca_amplitude      : ', Model%ca_amplitude
       print *, ' nsmooth           : ', Model%nsmooth
       print *, ' ca_closure        : ', Model%ca_closure

ccpp/data/GFS_typedefs.meta

@@ -4316,37 +4316,6 @@
   dimensions = ()
   type = real
   kind = kind_phys
-[nca]
-  standard_name = number_of_independent_cellular_automata
-  long_name = number of independent cellular automata
-  units = count
-  dimensions = ()
-  type = integer
-[nlives]
-  standard_name = cellular_automata_lifetime
-  long_name = cellular automata lifetime
-  units = count
-  dimensions = ()
-  type = integer
-[ncells]
-  standard_name = cellular_automata_finer_grid
-  long_name = cellular automata finer grid
-  units = count
-  dimensions = ()
-  type = integer
-[nfracseed]
-  standard_name = cellular_automata_seed_probability
-  long_name = cellular automata seed probability
-  units = fraction
-  dimensions = ()
-  type = real
-  kind = kind_phys
-[nseed]
-  standard_name = cellular_automata_seed_frequency
-  long_name = cellular automata seed frequency in units of time steps
-  units = count
-  dimensions = ()
-  type = integer
 [do_ca]
   standard_name = flag_for_cellular_automata
   long_name = cellular automata main switch
@@ -4383,24 +4352,6 @@
   units = flag
   dimensions = ()
   type = logical
-[ca_smooth]
-  standard_name = flag_for_gaussian_spatial_filter
-  long_name = switch for gaussian spatial filter
-  units = flag
-  dimensions = ()
-  type = logical
-[iseed_ca]
-  standard_name = seed_for_random_number_generation_in_cellular_automata_scheme
-  long_name = seed for random number generation in ca scheme
-  units = none
-  dimensions = ()
-  type = integer
-[nspinup]
-  standard_name = number_of_iterations_to_spin_up_cellular_automata
-  long_name = number of iterations to spin up the ca
-  units = count
-  dimensions = ()
-  type = integer
 [nthresh]
   standard_name = threshold_for_perturbed_vertical_velocity
   long_name = threshold used for perturbed vertical velocity

ccpp/driver/GFS_restart.F90

@@ -101,6 +101,10 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
     if (Model%imfdeepcnv == 3) then
       Restart%num2d = Restart%num2d + 1
     endif
+    ! CA
+    if (Model%imfdeepcnv == 2 .and. Model%do_ca) then
+      Restart%num2d = Restart%num2d + 1
+    endif
     ! NoahMP
     if (Model%lsm == Model%lsm_noahmp) then
       Restart%num2d = Restart%num2d + 10
@@ -184,8 +188,17 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
 !      print *,'in restart 2d field, Restart%name2d(',offset+idx,')=',trim(Restart%name2d(offset+idx))
     enddo
 
-    !--- RAP/HRRR-specific variables, 2D
     num = offset + ndiag_rst
+    !--- Celluluar Automaton, 2D
+    !CA
+    if (Model%imfdeepcnv == 2 .and. Model%do_ca) then
+      num = num + 1
+      Restart%name2d(num) = 'ca_condition'
+      do nb = 1,nblks
+        Restart%data(nb,num)%var2p => Coupling(nb)%condition(:)
+      enddo
+    endif
+    !--- RAP/HRRR-specific variables, 2D
     ! GF
     if (Model%imfdeepcnv == Model%imfdeepcnv_gf) then
       num = num + 1

stochastic_physics/stochastic_physics_wrapper.F90

@@ -35,11 +35,9 @@ module stochastic_physics_wrapper_mod
   real(kind=kind_phys), dimension(:,:), allocatable, save :: stype
 
   ! For cellular automata
-  real(kind=kind_phys), dimension(:,:,:), allocatable, save :: ugrs
-  real(kind=kind_phys), dimension(:,:,:), allocatable, save :: qgrs
-  real(kind=kind_phys), dimension(:,:),   allocatable, save :: pgr
-  real(kind=kind_phys), dimension(:,:,:), allocatable, save :: vvl
-  real(kind=kind_phys), dimension(:,:,:), allocatable, save :: prsl
+  real(kind=kind_phys), dimension(:,:),   allocatable, save :: sst
+  real(kind=kind_phys), dimension(:,:),   allocatable, save :: lmsk
+  real(kind=kind_phys), dimension(:,:),   allocatable, save :: lake
   real(kind=kind_phys), dimension(:,:),   allocatable, save :: condition
   real(kind=kind_phys), dimension(:,:),   allocatable, save :: ca_deep_cpl, ca_turb_cpl, ca_shal_cpl
   real(kind=kind_phys), dimension(:,:),   allocatable, save :: ca_deep_diag,ca_turb_diag,ca_shal_diag
@@ -276,15 +274,31 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
 
     endif initalize_stochastic_physics
 
-    ! Cellular automata code is identical for initialization (kstep=0) and time integration (kstep>0)
     if(GFS_Control%do_ca)then
+       if (GFS_Control%kdt==0) then
+
+          do nb=1,Atm_block%nblks
+             GFS_Data(nb)%Intdiag%ca_deep(:)  = 0.
+             GFS_Data(nb)%Intdiag%ca_turb(:)  = 0.
+             GFS_Data(nb)%Intdiag%ca_shal(:)  = 0.
+             GFS_Data(nb)%Coupling%ca_deep(:) = 0.
+             GFS_Data(nb)%Coupling%ca_turb(:) = 0.
+             GFS_Data(nb)%Coupling%ca_shal(:) = 0.
+             GFS_Data(nb)%Coupling%ca1(:) = 0.
+             GFS_Data(nb)%Coupling%ca2(:) = 0.
+             GFS_Data(nb)%Coupling%ca3(:) = 0.
+             GFS_Data(nb)%Intdiag%ca1(:)  = 0.
+             GFS_Data(nb)%Intdiag%ca2(:)  = 0.
+             GFS_Data(nb)%Intdiag%ca3(:)  = 0.
+          enddo
+
+       else
+
        if(GFS_Control%ca_sgs)then
          ! Allocate contiguous arrays; copy in as needed
-         allocate(ugrs        (1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
-         allocate(qgrs        (1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
-         allocate(pgr         (1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
-         allocate(vvl         (1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
-         allocate(prsl        (1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
+         allocate(sst         (1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
+         allocate(lmsk        (1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
+         allocate(lake        (1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
          allocate(ca_deep_diag(1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
          allocate(ca_turb_diag(1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
          allocate(ca_shal_diag(1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
@@ -293,22 +307,21 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
          allocate(ca_turb_cpl (1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
          allocate(ca_shal_cpl (1:Atm_block%nblks,maxval(GFS_Control%blksz)                   ))
          do nb=1,Atm_block%nblks
-             ugrs       (nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Statein%ugrs(:,:)
-             qgrs       (nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Statein%qgrs(:,:,1)
-             pgr        (nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Statein%pgr(:)
-             vvl        (nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Statein%vvl(:,:)
-             prsl       (nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Statein%prsl(:,:)
+             sst        (nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Sfcprop%tsfco(:)
+             lmsk       (nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Sfcprop%slmsk(:)
+             lake       (nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Sfcprop%lakefrac(:)
              condition  (nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Coupling%condition(:)
              ca_deep_cpl(nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Coupling%ca_deep(:)
              ca_turb_cpl(nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Coupling%ca_turb(:)
              ca_shal_cpl(nb,1:GFS_Control%blksz(nb))   = GFS_Data(nb)%Coupling%ca_shal(:)
          enddo
-         call cellular_automata_sgs(GFS_Control%kdt,ugrs,qgrs,pgr,vvl,prsl,condition,ca_deep_cpl,ca_turb_cpl,ca_shal_cpl, &
-            ca_deep_diag,ca_turb_diag,ca_shal_diag,Atm(mygrid)%domain_for_coupler,Atm_block%nblks,                        &
+         call cellular_automata_sgs(GFS_Control%kdt,GFS_control%dtf,GFS_control%restart,GFS_Control%first_time_step,       &
+            sst,lmsk,lake,condition,ca_deep_cpl,ca_turb_cpl,ca_shal_cpl,ca_deep_diag,ca_turb_diag,                        &
+            ca_shal_diag,Atm(mygrid)%domain_for_coupler,Atm_block%nblks,                                                  &
             Atm_block%isc,Atm_block%iec,Atm_block%jsc,Atm_block%jec,Atm(mygrid)%npx,Atm(mygrid)%npy, GFS_Control%levs,    &
-            GFS_Control%nca,GFS_Control%ncells,GFS_Control%nlives,GFS_Control%nfracseed,                                  &
-            GFS_Control%nseed,GFS_Control%nthresh,GFS_Control%ca_global,GFS_Control%ca_sgs,GFS_Control%iseed_ca,          &
-            GFS_Control%ca_smooth,GFS_Control%nspinup,Atm_block%blksz(1),GFS_Control%master,GFS_Control%communicator)
+            GFS_Control%nthresh,GFS_Control%rcell,GFS_Control%nca,GFS_Control%scells,GFS_Control%tlives,GFS_Control%nfracseed,    &
+            GFS_Control%nseed,GFS_Control%ca_global,GFS_Control%ca_sgs,GFS_Control%iseed_ca,          &
+            GFS_Control%ca_smooth,GFS_Control%nspinup,GFS_Control%ca_trigger,Atm_block%blksz(1),GFS_Control%master,GFS_Control%communicator)
          ! Copy contiguous data back as needed
          do nb=1,Atm_block%nblks
              GFS_Data(nb)%Intdiag%ca_deep(:)  = ca_deep_diag(nb,1:GFS_Control%blksz(nb))
@@ -318,11 +331,9 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
              GFS_Data(nb)%Coupling%ca_turb(:) = ca_turb_cpl (nb,1:GFS_Control%blksz(nb))
              GFS_Data(nb)%Coupling%ca_shal(:) = ca_shal_cpl (nb,1:GFS_Control%blksz(nb))
          enddo
-         deallocate(ugrs        )
-         deallocate(qgrs        )
-         deallocate(pgr         )
-         deallocate(vvl         )
-         deallocate(prsl        )
+         deallocate(sst         )
+         deallocate(lmsk        )
+         deallocate(lake        )
          deallocate(condition   )
          deallocate(ca_deep_cpl )
          deallocate(ca_turb_cpl )
@@ -339,9 +350,10 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
           allocate(ca1_diag(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
           allocate(ca2_diag(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
           allocate(ca3_diag(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
-          call cellular_automata_global(GFS_Control%kdt,ca1_cpl,ca2_cpl,ca3_cpl,ca1_diag,ca2_diag,ca3_diag,Atm(mygrid)%domain_for_coupler, &
+          call cellular_automata_global(GFS_Control%kdt,GFS_Control%first_time_step,ca1_cpl,ca2_cpl,ca3_cpl,ca1_diag,ca2_diag,ca3_diag,     &
+            Atm(mygrid)%domain_for_coupler, &
             Atm_block%nblks,Atm_block%isc,Atm_block%iec,Atm_block%jsc,Atm_block%jec,Atm(mygrid)%npx,Atm(mygrid)%npy,GFS_Control%levs,      &
-            GFS_Control%nca_g,GFS_Control%ncells_g,GFS_Control%nlives_g,GFS_Control%nfracseed,GFS_Control%nseed_g,GFS_Control%nthresh,     &
+            GFS_Control%nca_g,GFS_Control%ncells_g,GFS_Control%nlives_g,GFS_Control%nfracseed,GFS_Control%nseed_g,                         &
             GFS_Control%ca_global,GFS_Control%ca_sgs,GFS_Control%iseed_ca,GFS_Control%ca_smooth,GFS_Control%nspinup,Atm_block%blksz(1),    &
             GFS_Control%nsmooth,GFS_Control%ca_amplitude,GFS_Control%master,GFS_Control%communicator)
           ! Copy contiguous data back
@@ -360,7 +372,10 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
           deallocate(ca2_diag)
           deallocate(ca3_diag)
        endif
-    endif
+
+    endif !kdt =0
+
+    endif !do_ca
 
   end subroutine stochastic_physics_wrapper