Update of GF aerosol treatment and tunings

physics/cu_gf_driver.F90

@@ -7,7 +7,7 @@ module cu_gf_driver
    ! DH* TODO: replace constants with arguments to cu_gf_driver_run
    use physcons  , g => con_g, cp => con_cp, xlv => con_hvap, r_v => con_rv
    use machine   , only: kind_phys
-   use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap,fct1d3
+   use cu_gf_deep, only: cu_gf_deep_run,neg_check,fct1d3
    use cu_gf_sh  , only: cu_gf_sh_run
 
    implicit none
@@ -75,14 +75,14 @@ end subroutine cu_gf_driver_finalize
 !!
 !>\section gen_gf_driver GSD GF Cumulus Scheme General Algorithm
 !> @{
-      subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
-               forcet,forceqv_spechum,phil,raincv,qv_spechum,t,cld1d,           &
-               us,vs,t2di,w,qv2di_spechum,p2di,psuri,                           &
-               hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw,                        &
+      subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,             &
+               cactiv,cactiv_m,forcet,forceqv_spechum,phil,raincv,qv_spechum,   &
+               t,cld1d,us,vs,t2di,w,qv2di_spechum,p2di,psuri,                   &
+               hbot,htop,kcnv,xland,hfx2,qfx2,aod_gf,cliw,clcw,                 &
                pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv,                &
                flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend,           &
                dtend,dtidx,ntqv,ntiw,ntcw,index_of_temperature,index_of_x_wind, &
-               index_of_y_wind,index_of_process_scnv,index_of_process_dcnv,     &
+               index_of_y_wind,index_of_process_scnv,index_of_process_dcnv,     & 
                ldiag3d,qci_conv,errmsg,errflg)
 !-------------------------------------------------------------
       implicit none
@@ -97,22 +97,24 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
      !integer, parameter :: ichoicem=5	! 0 2 5 13
       integer, parameter :: ichoicem=13	! 0 2 5 13
       integer, parameter :: ichoice_s=3	! 0 1 2 3
-      real(kind=kind_phys), parameter :: aodccn=0.1
+
+      real(kind=kind_phys), parameter :: aodc0=0.14
+      real(kind=kind_phys), parameter :: aodreturn=30.   
       real(kind=kind_phys) :: dts,fpi,fp
       integer, parameter :: dicycle=0 ! diurnal cycle flag
       integer, parameter :: dicycle_m=0 !- diurnal cycle flag
       integer            :: ishallow_g3 ! depend on imfshalcnv
 !-------------------------------------------------------------
    integer      :: its,ite, jts,jte, kts,kte 
    integer, intent(in   ) :: im,km,ntracer
-   logical, intent(in   ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend
+   logical, intent(in   ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend,flag_init
    logical, intent(in   ) :: ldiag3d
 
    real(kind=kind_phys), optional, intent(inout)            :: dtend(:,:,:)
    integer, intent(in)                                      :: dtidx(:,:), &
         index_of_x_wind, index_of_y_wind, index_of_temperature,            &
         index_of_process_scnv, index_of_process_dcnv, ntqv, ntcw, ntiw
-   
+
    real(kind=kind_phys),  dimension( : , : ), intent(in    ) :: forcet,forceqv_spechum,w,phil
    real(kind=kind_phys),  dimension( : , : ), intent(inout ) :: t,us,vs
    real(kind=kind_phys),  dimension( : , : ), intent(inout ) :: qci_conv
@@ -133,14 +135,15 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
    ! Specific humidity from FV3
    real(kind=kind_phys), dimension (:,:), intent(in) :: qv2di_spechum
    real(kind=kind_phys), dimension (:,:), intent(inout) :: qv_spechum
+   real(kind=kind_phys), dimension (:), intent(inout) :: aod_gf
    ! Local water vapor mixing ratios and cloud water mixing ratios
    real(kind=kind_phys), dimension (im,km) :: qv2di, qv, forceqv, cnvw
    !
    real(kind=kind_phys), dimension(:),intent(in) :: garea
    real(kind=kind_phys), intent(in   ) :: dt 
 
    integer, intent(in   ) :: imfshalcnv
-   integer, dimension(:), intent(inout) :: cactiv
+   integer, dimension(:), intent(inout) :: cactiv,cactiv_m
 
    character(len=*), intent(out) :: errmsg
    integer,          intent(out) :: errflg
@@ -151,6 +154,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
    real(kind=kind_phys), dimension (im,4)  :: rand_clos
    real(kind=kind_phys), dimension (im,km,11) :: gdc,gdc2
    real(kind=kind_phys), dimension (im)    :: ht
+   real(kind=kind_phys), dimension (im)    :: ccn_gf,ccn_m
+   real(kind=kind_phys) :: ccnclean
    real(kind=kind_phys), dimension (im)    :: dx
    real(kind=kind_phys), dimension (im,km) :: outt,outq,outqc,phh,subm,cupclw,cupclws
    real(kind=kind_phys), dimension (im,km) :: dhdt,zu,zus,zd,phf,zum,zdm,outum,outvm
@@ -179,9 +184,9 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
 ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off
 ! convection for this call only and at that particular gridpoint
 !
-   real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi
+   real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi,clw_ten
    real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg
-   real(kind=kind_phys), dimension (im)    :: ccn,z1,psur,cuten,cutens,cutenm
+   real(kind=kind_phys), dimension (im)    :: z1,psur,cuten,cutens,cutenm
    real(kind=kind_phys), dimension (im)    :: umean,vmean,pmean
    real(kind=kind_phys), dimension (im)    :: xmbs,xmbs2,xmb,xmbm,xmb_dumm,mconv
 
@@ -190,7 +195,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
    integer :: high_resolution
    real(kind=kind_phys)    :: clwtot,clwtot1,excess,tcrit,tscl_kf,dp,dq,sub_spread,subcenter
    real(kind=kind_phys)    :: dsubclw,dsubclws,dsubclwm,dtime_max,ztm,ztq,hfm,qfm,rkbcon,rktop
-   real(kind=kind_phys), dimension(km)   :: massflx,trcflx_in1,clw_in1,clw_ten1,po_cup
+   real(kind=kind_phys), dimension(km)   :: massflx,trcflx_in1,clw_in1,po_cup
 !  real(kind=kind_phys), dimension(km)   :: trcflx_in2,clw_in2,clw_ten2
    real(kind=kind_phys), dimension (im)  :: flux_tun,tun_rad_mid,tun_rad_shall,tun_rad_deep
    character*50 :: ierrc(im),ierrcm(im)
@@ -280,7 +285,6 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
 ! dx for scale awareness
 !    dx=40075000./float(lonf)
 !    tscl_kf=dx/25000.
-     ccn(its:ite)=150.
 
      if (imfshalcnv == 3) then
       ishallow_g3 = 1
@@ -335,7 +339,24 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
      do i= its,itf
       forcing(i,:)=0.
       forcing2(i,:)=0.
-      ccn(i)=100.
+      ccn_gf(i) = 0.
+      ccn_m(i) = 0.
+
+      ! set aod and ccn
+      if (flag_init) then
+        aod_gf(i)=aodc0
+      else
+        if((cactiv(i).eq.0) .and. (cactiv_m(i).eq.0))then
+          if(aodc0>aod_gf(i)) aod_gf(i)=aod_gf(i)+((aodc0-aod_gf(i))*(dt/(aodreturn*60)))
+          if(aod_gf(i)>aodc0) aod_gf(i)=aodc0
+        endif
+      endif
+
+      ccn_gf(i)=max(5., (aod_gf(i)/0.0027)**(1/0.640))
+      ccn_m(i)=ccn_gf(i)
+
+      ccnclean=max(5., (aodc0/0.0027)**(1/0.640))
+
       hbot(i)  =kte
       htop(i)  =kts
       raincv(i)=0.
@@ -558,7 +579,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
               ,dicycle_m       &
               ,ichoicem       &
               ,ipr           &
-              ,ccn           &
+              ,ccn_m         &
+              ,ccnclean      &
               ,dt            &
               ,imid_gf       &
               ,kpbli         &
@@ -638,7 +660,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
               ,dicycle       &
               ,ichoice       &
               ,ipr           &
-              ,ccn           &
+              ,ccn_gf        &
+              ,ccnclean      &
               ,dt            &
               ,0             &
 
@@ -761,7 +784,9 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
             massflx(:)=0.
             trcflx_in1(:)=0.
             clw_in1(:)=0.
-            clw_ten1(:)=0.
+            do k=kts,ktf 
+              clw_ten(i, k)=0.
+            enddo   
             po_cup(:)=0.
             kstop=kts
             if(ktopm(i).gt.kts .or. ktop(i).gt.kts)kstop=max(ktopm(i),ktop(i))
@@ -851,20 +876,22 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
              massflx   (1)=0.
              trcflx_in1(1)=0.
              call fct1d3 (kstop,kte,dtime_max,po_cup,                  &
-                            clw_in1,massflx,trcflx_in1,clw_ten1,g)
+                            clw_in1,massflx,trcflx_in1,clw_ten(i,:),g)
 
              do k=1,kstop
                tem  = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i)    &
                       +outqcm(i,k)*cutenm(i)                           &
-                      +clw_ten1(k)                                     &
+                      +clw_ten(i,k)                                    &
                          )
-               tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf))
-               if (clcw(i,k) .gt. -999.0) then
-                cliw(i,k) = max(0.,cliw(i,k) + tem * tem1)            ! ice
-                clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1))       ! water
-               else
-                cliw(i,k) = max(0.,cliw(i,k) + tem)
-               endif
+               !tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf))
+               !if (clcw(i,k) .gt. -999.0) then
+               ! cliw(i,k) = max(0.,cliw(i,k) + tem * tem1)            ! ice
+               ! clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1))       ! water
+               !else
+               ! cliw(i,k) = max(0.,cliw(i,k) + tem)
+               !endif
+               if(t(i,k).le.270.) cliw(i,k) = max(0.,cliw(i,k) + tem) ! HCB 
+               if(t(i,k).gt.270) clcw(i,k) = max(0.,clcw(i,k) + tem) 
 
              enddo
 
@@ -893,6 +920,29 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
                  cactiv(i)=0
                  if(pretm(i).gt.0)raincv(i)=.001*cutenm(i)*pretm(i)*dt
               endif   ! pret > 0
+
+              if(pretm(i).gt.0)then
+                 cactiv_m(i)=1
+              else
+                 cactiv_m(i)=0
+              endif
+
+              ! Unify ccn
+              if(ccn_m(i).lt.ccn_gf(i))then
+                ccn_gf(i)=ccn_m(i)
+              endif
+
+              if(ccn_gf(i)<0) ccn_gf(i)=0
+
+              ! Convert ccn back to aod
+              aod_gf(i)=0.0027*(ccn_gf(i)**0.64)
+              if(aod_gf(i)<0.007)then
+                aod_gf(i)=0.007
+                ccn_gf(i)=(aod_gf(i)/0.0027)**(1/0.640)
+              elseif(aod_gf(i)>aodc0)then
+                aod_gf(i)=aodc0
+                ccn_gf(i)=(aod_gf(i)/0.0027)**(1/0.640)
+              endif
             enddo
  100    continue
 !
@@ -958,7 +1008,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
             if(qidx>=1) then
               do k=kts,ktf
                 do i=its,itf
-                  tem = (cuten(i)*outq(i,k) + cutenm(i)*outqm(i,k))* dt
+                  tem = (cuten(i)*outq(i,k) + cutenm(i)*outqm(i,k))* dt  
                   tem = tem/(1.0_kind_phys+tem)
                   dtend(i,k,qidx) = dtend(i,k,qidx) + tem
                 enddo
@@ -969,14 +1019,14 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv,                &
             do k=kts,ktf
               do i=its,itf
                 tem_shal = dt*(outqcs(i,k)*cutens(i)+outqcm(i,k)*cutenm(i))
-                tem_deep = dt*(outqc(i,k)*cuten(i)+clw_ten1(k))
+                tem_deep = dt*(outqc(i,k)*cuten(i)+clw_ten(i,k))
                 tem  = tem_shal+tem_deep
                 tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf))
                 weight_sum = abs(tem_shal)+abs(tem_deep)
                 if(weight_sum<1e-12) then
                   cycle
                 endif
-                
+
                 if (clcw_save(i,k) .gt. -999.0) then
                   cliw_both = max(0.,cliw_save(i,k) + tem * tem1) - cliw_save(i,k)
                   clcw_both = max(0.,clcw_save(i,k) + tem) - clcw_save(i,k)

physics/cu_gf_driver.meta

@@ -119,6 +119,14 @@
   kind = kind_phys
   intent = in
   optional = F
+[flag_init]
+  standard_name = flag_for_first_time_step
+  long_name = flag signaling first time step for time integration loop
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
 [cactiv]
   standard_name = conv_activity_counter
   long_name = convective activity memory
@@ -127,6 +135,14 @@
   type = integer
   intent = inout
   optional = F
+[cactiv_m]
+  standard_name = mid_conv_activity_counter
+  long_name = mid-level cloud convective activity memory
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent = inout
+  optional = F
 [forcet]
   standard_name = temperature_tendency_due_to_dynamics
   long_name = temperature tendency due to dynamics only
@@ -303,6 +319,15 @@
   kind = kind_phys
   intent = in
   optional = F
+[aod_gf]
+  standard_name = aod_gf_deep
+  long_name = aerosol optical depth used in Grell-Freitas Convective Parameterization
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+  optional = F
 [cliw]
   standard_name = ice_water_mixing_ratio_convective_transport_tracer
   long_name = ratio of mass of ice water to mass of dry air plus vapor (without condensates) in the convectively transported tracer array

physics/cu_gf_driver_post.F90

@@ -20,7 +20,7 @@ end subroutine cu_gf_driver_post_finalize
 !> \section arg_table_cu_gf_driver_post_run Argument Table
 !! \htmlinclude cu_gf_driver_post_run.html
 !!
-   subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, conv_act, errmsg, errflg)
+   subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, cactiv_m, conv_act, conv_act_m, errmsg, errflg)
 
       use machine, only: kind_phys
 
@@ -33,7 +33,9 @@ subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, conv_act, er
       real(kind_phys),  intent(out) :: prevst(:,:)
       real(kind_phys),  intent(out) :: prevsq(:,:)
       integer,          intent(in)  :: cactiv(:)
+      integer,          intent(in)  :: cactiv_m(:)
       real(kind_phys),  intent(out) :: conv_act(:)
+      real(kind_phys),  intent(out) :: conv_act_m(:)
       character(len=*), intent(out) :: errmsg
       integer, intent(out)          :: errflg
 
@@ -53,6 +55,11 @@ subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, conv_act, er
         else
           conv_act(i)=0.0
         endif
+        if (cactiv_m(i).gt.0) then
+          conv_act_m(i) = conv_act_m(i)+1.0
+        else
+          conv_act_m(i)=0.0
+        endif
       enddo
 
    end subroutine cu_gf_driver_post_run

physics/cu_gf_driver_post.meta

@@ -59,6 +59,14 @@
   type = integer
   intent = in
   optional = F
+[cactiv_m]
+  standard_name = conv_mid_activity_counter
+  long_name = midlevel convective activity memory
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent = in
+  optional = F
 [conv_act]
   standard_name = gf_memory_counter
   long_name = Memory counter for GF
@@ -68,6 +76,15 @@
   kind = kind_phys
   intent = out
   optional = F
+[conv_act_m]
+  standard_name = gf_mid_memory_counter
+  long_name = Memory counter for GF midlevel
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = out
+  optional = F
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/cu_gf_driver_pre.F90

@@ -21,7 +21,8 @@ end subroutine cu_gf_driver_pre_finalize
 !! \htmlinclude cu_gf_driver_pre_run.html
 !!
    subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, prevst, prevsq, &
-                                    forcet, forceq, cactiv, conv_act, errmsg, errflg)
+                                    forcet, forceq, cactiv, cactiv_m, conv_act, conv_act_m,         &
+                                    errmsg, errflg)
 
       use machine, only: kind_phys
 
@@ -39,7 +40,9 @@ subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q,
       real(kind_phys),  intent(out) :: forcet(:,:)
       real(kind_phys),  intent(out) :: forceq(:,:)
       integer,          intent(out) :: cactiv(:)
+      integer,          intent(out) :: cactiv_m(:)
       real(kind_phys),  intent(in)  :: conv_act(:)
+      real(kind_phys),  intent(in)  :: conv_act_m(:)
       character(len=*), intent(out) :: errmsg
       integer,          intent(out) :: errflg
 
@@ -68,6 +71,7 @@ subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q,
       endif
 
       cactiv(:)=nint(conv_act(:))
+      cactiv_m(:)=nint(conv_act_m(:))
 
    end subroutine cu_gf_driver_pre_run