Ca develop

atmos_model.F90

@@ -312,14 +312,18 @@ subroutine update_atmos_radiation_physics (Atmos)
     end if
 
     if(IPD_Control%do_ca)then
-       ! DH* The current implementation of cellular_automata assumes that all blocksizes are the
-       ! same, this is tested in the initialization call to cellular_automata, no need to redo *DH
-       call cellular_automata(IPD_Control%kdt, IPD_Data(:)%Statein, IPD_Data(:)%Coupling, IPD_Data(:)%Intdiag, &
-                              Atm_block%nblks, IPD_Control%levs, IPD_Control%nca, IPD_Control%ncells,          &
-                              IPD_Control%nlives, IPD_Control%nfracseed, IPD_Control%nseed,                    &
-                              IPD_Control%nthresh, IPD_Control%ca_global, IPD_Control%ca_sgs,                  &
-                              IPD_Control%iseed_ca, IPD_Control%ca_smooth, IPD_Control%nspinup,                &
-                              Atm_block%blksz(1))
+       if(IPD_Control%ca_sgs == .true.)then
+          call cellular_automata_sgs(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
+            IPD_Control%nca,IPD_Control%ncells,IPD_Control%nlives,IPD_Control%nfracseed,&
+            IPD_Control%nseed,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
+            IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1))
+       endif
+       if(IPD_Control%ca_global == .true.)then
+          call cellular_automata_global(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
+            IPD_Control%nca_g,IPD_Control%ncells_g,IPD_Control%nlives_g,IPD_Control%nfracseed,&
+            IPD_Control%nseed_g,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
+            IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1),IPD_Control%nsmooth,IPD_Control%ca_amplitude)
+      endif
     endif
 
 !--- if coupled, assign coupled fields
@@ -656,12 +660,21 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
          call mpp_error(FATAL, 'Logic errror: cellular_automata not compatible with non-uniform blocksizes')
       end if
       ! *DH
-      call cellular_automata(IPD_Control%kdt, IPD_Data(:)%Statein, IPD_Data(:)%Coupling, IPD_Data(:)%Intdiag, &
-                             Atm_block%nblks, IPD_Control%levs, IPD_Control%nca, IPD_Control%ncells,          &
-                             IPD_Control%nlives, IPD_Control%nfracseed, IPD_Control%nseed,                    &
-                             IPD_Control%nthresh, IPD_Control%ca_global, IPD_Control%ca_sgs,                  &
-                             IPD_Control%iseed_ca, IPD_Control%ca_smooth, IPD_Control%nspinup,                &
-                             Atm_block%blksz(1))
+      if(IPD_Control%do_ca)then
+       if(IPD_Control%ca_sgs == .true.)then
+          call cellular_automata_sgs(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
+            IPD_Control%nca,IPD_Control%ncells,IPD_Control%nlives,IPD_Control%nfracseed,&
+            IPD_Control%nseed,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
+            IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1))
+       endif
+       if(IPD_Control%ca_global == .true.)then
+          call cellular_automata_global(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
+            IPD_Control%nca_g,IPD_Control%ncells_g,IPD_Control%nlives_g,IPD_Control%nfracseed,&
+            IPD_Control%nseed_g,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
+            IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1),IPD_Control%nsmooth,IPD_Control%ca_amplitude)
+       endif
+
+    endif
    endif
 
    Atm(mytile)%flagstruct%do_skeb = IPD_Control%do_skeb

gfsphysics/GFS_layer/GFS_diagnostics.F90

@@ -2028,13 +2028,13 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
 
     idx = idx + 1
     ExtDiag(idx)%axes = 2
-    ExtDiag(idx)%name = 'ca_out'
+    ExtDiag(idx)%name = 'ca1'
     ExtDiag(idx)%desc = 'Cellular Automata'
     ExtDiag(idx)%unit = '%'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     allocate (ExtDiag(idx)%data(nblks))
     do nb = 1,nblks
-      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ca_out(:)
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ca1(:)
     enddo
 
     idx = idx + 1

gfsphysics/GFS_layer/GFS_driver.F90

@@ -165,7 +165,7 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop,     &
     integer :: ntrac
     integer :: ix
 #ifndef CCPP
-    integer :: blocksize
+    integer :: blocksize,k
     real(kind=kind_phys), allocatable :: si(:)
     real(kind=kind_phys), parameter   :: p_ref = 101325.0d0
 #endif
@@ -194,10 +194,9 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop,     &
                      Init_parm%iau_offset, Init_parm%bdat,         &
                      Init_parm%cdat, Init_parm%tracer_names,       &
                      Init_parm%input_nml_file, Init_parm%tile_num, &
-                     Init_parm%blksz                               &
+                     Init_parm%blksz,Init_parm%ak, Init_parm%bk   &
 #ifdef CCPP
-                    ,Init_parm%ak, Init_parm%bk,                   &
-                     Init_parm%restart, Init_parm%hydrostatic,     &
+                     ,Init_parm%restart, Init_parm%hydrostatic,     &
                      communicator, ntasks, nthrds                  &
 #endif
                      )
@@ -469,6 +468,29 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop,     &
     !--- this note is placed here to alert users to study
     !--- the FV3GFS_io.F90 module
 
+#ifndef CCPP
+    if(Model%do_ca .and. Model%ca_global)then
+
+       do nb = 1,nblks
+         do k=1,Model%levs
+          if (Model%si(k) .lt. 0.1 .and. Model%si(k) .gt. 0.025) then
+           Coupling(nb)%vfact_ca(k) = (Model%si(k)-0.025)/(0.1-0.025)
+          else if (Model%si(k) .lt. 0.025) then
+           Coupling(nb)%vfact_ca(k) = 0.0
+          else
+           Coupling(nb)%vfact_ca(k) = 1.0
+          endif
+        enddo
+       enddo
+
+       do nb = 1,nblks
+        Coupling(nb)%vfact_ca(2)=Coupling(nb)%vfact_ca(3)*0.5
+        Coupling(nb)%vfact_ca(1)=0.0
+       enddo
+
+    endif
+#endif
+
   end subroutine GFS_initialize
 
 
@@ -710,8 +732,7 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
     !--- local variables
     integer :: k, i
     real(kind=kind_phys) :: upert, vpert, tpert, qpert, qnew,sppt_vwt
-    !real(kind=kind_phys),dimension(size(Statein%tgrs,1),size(Statein%tgrs,2)) :: tconvtend, &
-    !                     qconvtend,uconvtend,vconvtend
+    real(kind=kind_phys),dimension(size(Statein%tgrs,1),size(Statein%tgrs,2)) :: ca1
 
      if (Model%do_sppt) then
        do k = 1,size(Statein%tgrs,2)
@@ -739,27 +760,11 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
            Diag%sppt_wts(i,Model%levs-k+1)=Coupling%sppt_wts(i,k)
 
 
-          ! if(Model%isppt_deep)then
-
-           ! tconvtend(i,k)=Coupling%tconvtend(i,k)
-           ! qconvtend(i,k)=Coupling%qconvtend(i,k)
-           ! uconvtend(i,k)=Coupling%uconvtend(i,k)
-           ! vconvtend(i,k)=Coupling%vconvtend(i,k)           
-
-
-           ! upert = (Stateout%gu0(i,k)   - Statein%ugrs(i,k) - uconvtend(i,k)) + uconvtend(i,k) * Coupling%sppt_wts(i,k)
-           ! vpert = (Stateout%gv0(i,k)   - Statein%vgrs(i,k) - vconvtend(i,k)) + vconvtend(i,k) * Coupling%sppt_wts(i,k)
-           ! tpert = (Stateout%gt0(i,k)   - Statein%tgrs(i,k) - Tbd%dtdtr(i,k) - tconvtend(i,k)) + tconvtend(i,k) * Coupling%sppt_wts(i,k)
-           ! qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1) - qconvtend(i,k)) + qconvtend(i,k) * Coupling%sppt_wts(i,k)
-
-           !else
-
             upert = (Stateout%gu0(i,k)   - Statein%ugrs(i,k))   * Coupling%sppt_wts(i,k)
             vpert = (Stateout%gv0(i,k)   - Statein%vgrs(i,k))   * Coupling%sppt_wts(i,k)
             tpert = (Stateout%gt0(i,k)   - Statein%tgrs(i,k) - Tbd%dtdtr(i,k)) * Coupling%sppt_wts(i,k)
             qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1)) * Coupling%sppt_wts(i,k)
 
-           !endif
 
            Stateout%gu0(i,k)  = Statein%ugrs(i,k)+upert
            Stateout%gv0(i,k)  = Statein%vgrs(i,k)+vpert
@@ -773,22 +778,6 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
          enddo
        enddo
 
-       !if(Model%isppt_deep == .true.)then
-       !  Sfcprop%tprcp(:) = Sfcprop%tprcp(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
-       !  Diag%totprcp(:)  = Diag%totprcp(:)  + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:) 
-       !  Diag%cnvprcp(:)  = Diag%cnvprcp(:)  + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
-       !!  ! bucket precipitation adjustment due to sppt                                                                                                                                                                     
-       !  Diag%totprcpb(:) = Diag%totprcpb(:)  + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
-       !  Diag%cnvprcpb(:) = Diag%cnvprcpb(:)  + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
-
-
-       !  if (Model%cplflx) then !Need to make proper adjustments for deep convection only perturbations
-       !    Coupling%rain_cpl(:) = Coupling%rain_cpl(:) + (Coupling%sppt_wts(:,15) - 1.0)*Tbd%drain_cpl(:)
-       !    Coupling%snow_cpl(:) = Coupling%snow_cpl(:) + (Coupling%sppt_wts(:,15) - 1.0)*Tbd%dsnow_cpl(:)
-       !  endif
-
-       !else
-
         ! instantaneous precip rate going into land model at the next time step
         Sfcprop%tprcp(:) = Coupling%sppt_wts(:,15)*Sfcprop%tprcp(:)
         Diag%totprcp(:)      = Diag%totprcp(:)      + (Coupling%sppt_wts(:,15) - 1 )*Diag%rain(:)
@@ -804,10 +793,71 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
            Coupling%snow_cpl(:) = Coupling%snow_cpl(:) + (Coupling%sppt_wts(:,15) - 1.0)*Tbd%dsnow_cpl(:)
         endif
 
-       !endif 
+     endif
+
+
+     if (Model%do_ca .and. Model%ca_global) then
+       do k = 1,size(Statein%tgrs,2)
+         do i = 1,size(Statein%tgrs,1)
+           sppt_vwt=1.0
+           if (Diag%zmtnblck(i).EQ.0.0) then
+              sppt_vwt=1.0
+           else
+              if (k.GT.Diag%zmtnblck(i)+2) then
+                 sppt_vwt=1.0
+              endif
+              if (k.LE.Diag%zmtnblck(i)) then
+                 sppt_vwt=0.0
+              endif
+              if (k.EQ.Diag%zmtnblck(i)+1) then
+                 sppt_vwt=0.333333
+              endif
+              if (k.EQ.Diag%zmtnblck(i)+2) then
+                 sppt_vwt=0.666667
+              endif
+           endif
+
+           ca1(i,k)=((Coupling%ca1(i)-1.)*sppt_vwt*Coupling%vfact_ca(k))+1.0
+
+           upert = (Stateout%gu0(i,k)   - Statein%ugrs(i,k))   * ca1(i,k)
+           vpert = (Stateout%gv0(i,k)   - Statein%vgrs(i,k))   * ca1(i,k)
+           tpert = (Stateout%gt0(i,k)   - Statein%tgrs(i,k) - Tbd%dtdtr(i,k)) * ca1(i,k)
+           qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1)) * ca1(i,k)
+
+           Stateout%gu0(i,k)  = Statein%ugrs(i,k)+upert
+           Stateout%gv0(i,k)  = Statein%vgrs(i,k)+vpert
+
+           !negative humidity check                                                                                                                                                                                                                     
+           qnew = Statein%qgrs(i,k,1)+qpert
+           if (qnew >= 1.0e-10) then
+              Stateout%gq0(i,k,1) = qnew
+              Stateout%gt0(i,k)   = Statein%tgrs(i,k) + tpert + Tbd%dtdtr(i,k)
+           endif
+
+         enddo
+       enddo
+
+
+
+        ! instantaneous precip rate going into land model at the next time step                                                                                                                                                                         
+        Sfcprop%tprcp(:) = ca1(:,15)*Sfcprop%tprcp(:)
+        Diag%totprcp(:)      = Diag%totprcp(:)      + (ca1(:,15) - 1 )*Diag%rain(:)
+        ! acccumulated total and convective preciptiation                                                                                                                                                                                               
+        Diag%cnvprcp(:)      = Diag%cnvprcp(:)      + (ca1(:,15) - 1 )*Diag%rainc(:)
+        ! bucket precipitation adjustment due to sppt                                                                                                                                                                                                   
+        Diag%totprcpb(:)      = Diag%totprcpb(:)      + (ca1(:,15) - 1 )*Diag%rain(:)
+        Diag%cnvprcpb(:)      = Diag%cnvprcpb(:)      + (ca1(:,15) - 1 )*Diag%rainc(:)
+
+        if (Model%cplflx) then
+           Coupling%rain_cpl(:) = Coupling%rain_cpl(:) + (ca1(:,15) - 1.0)*Tbd%drain_cpl(:)
+           Coupling%snow_cpl(:) = Coupling%snow_cpl(:) + (ca1(:,15) - 1.0)*Tbd%dsnow_cpl(:)
+        endif
 
      endif
 
+
+
+
      if (Model%do_shum) then
        Stateout%gq0(:,:,1) = Stateout%gq0(:,:,1)*(1.0 + Coupling%shum_wts(:,:))
      endif