Amend to code review 13424: code cleanup in response to code review

FV3: this commits #refs: 41390, fix to SKEB amplitude to removed time-step dependence and not SPPT peturbations below dividing streamline and other bug fixes

Change-Id: I4a1b3853212f39a301a597eb7a64603c281d0a50

atmos_cubed_sphere/driver/fvGFS/atmosphere.F90

@@ -1012,6 +1012,15 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block)
           enddo
         enddo
       enddo
+      if (.not. Atm(mytile)%flagstruct%hydrostatic) then
+         do k = 1, npz
+            do j = jsc,jec
+               do i = isc,iec
+                  Atm(n)%delz(i,j,k) = Atm(n)%delz(i,j,k)  + IAU_Data%delz_inc(i,j,k)*dt_atmos
+             enddo
+           enddo
+         enddo
+      endif
 !     add analysis increment to tracers to output from physics
       do nb = 1,Atm_block%nblks
          !if (nb.EQ.1) print*,'in block_update',IAU_Data%in_interval,IAU_Data%temp_inc(isc,jsc,30)

atmos_cubed_sphere/tools/fv_iau_mod.F90

@@ -60,13 +60,15 @@ module fv_iau_mod
     real,allocatable :: va_inc(:,:,:)
     real,allocatable :: temp_inc(:,:,:)
     real,allocatable :: delp_inc(:,:,:)
+    real,allocatable :: delz_inc(:,:,:)
     real,allocatable :: tracer_inc(:,:,:,:)
   end type iau_internal_data_type
   type iau_external_data_type
     real,allocatable :: ua_inc(:,:,:)
     real,allocatable :: va_inc(:,:,:)
     real,allocatable :: temp_inc(:,:,:)
     real,allocatable :: delp_inc(:,:,:)
+    real,allocatable :: delz_inc(:,:,:)
     real,allocatable :: tracer_inc(:,:,:,:)
     logical          :: in_interval = .false.
   end type iau_external_data_type
@@ -207,12 +209,14 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
     allocate(IAU_Data%va_inc(is:ie, js:je, km))
     allocate(IAU_Data%temp_inc(is:ie, js:je, km))
     allocate(IAU_Data%delp_inc(is:ie, js:je, km))
+    allocate(IAU_Data%delz_inc(is:ie, js:je, km))
     allocate(IAU_Data%tracer_inc(is:ie, js:je, km,ntracers))
 ! allocate arrays that will hold iau state
     allocate (iau_state%inc1%ua_inc(is:ie, js:je, km))
     allocate (iau_state%inc1%va_inc(is:ie, js:je, km))
     allocate (iau_state%inc1%temp_inc (is:ie, js:je, km))
     allocate (iau_state%inc1%delp_inc (is:ie, js:je, km))
+    allocate (iau_state%inc1%delz_inc (is:ie, js:je, km))
     allocate (iau_state%inc1%tracer_inc(is:ie, js:je, km,ntracers))
     iau_state%hr1=IPD_Control%iaufhrs(1)
     call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1)))
@@ -224,6 +228,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
        allocate (iau_state%inc2%va_inc(is:ie, js:je, km))
        allocate (iau_state%inc2%temp_inc (is:ie, js:je, km))
        allocate (iau_state%inc2%delp_inc (is:ie, js:je, km))
+       allocate (iau_state%inc2%delz_inc (is:ie, js:je, km))
        allocate (iau_state%inc2%tracer_inc(is:ie, js:je, km,ntracers))
        iau_state%hr2=IPD_Control%iaufhrs(2)
        call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
@@ -302,6 +307,7 @@ subroutine updateiauforcing(IPD_Control,IAU_Data)
             IAU_Data%va_inc(i,j,k)    =(delt*IAU_state%inc1%va_inc(i,j,k)    + (1.-delt)* IAU_state%inc2%va_inc(i,j,k))*rdt
             IAU_Data%temp_inc(i,j,k)  =(delt*IAU_state%inc1%temp_inc(i,j,k)  + (1.-delt)* IAU_state%inc2%temp_inc(i,j,k))*rdt
             IAU_Data%delp_inc(i,j,k)  =(delt*IAU_state%inc1%delp_inc(i,j,k)  + (1.-delt)* IAU_state%inc2%delp_inc(i,j,k))*rdt
+            IAU_Data%delz_inc(i,j,k)  =(delt*IAU_state%inc1%delz_inc(i,j,k)  + (1.-delt)* IAU_state%inc2%delz_inc(i,j,k))*rdt
             do l=1,ntracers
                IAU_Data%tracer_inc(i,j,k,l) =(delt*IAU_state%inc1%tracer_inc(i,j,k,l) + (1.-delt)* IAU_state%inc2%tracer_inc(i,j,k,l))*rdt
             enddo
@@ -327,6 +333,7 @@ subroutine setiauforcing(IPD_Control,IAU_Data)
           IAU_Data%va_inc(i,j,k)    =IAU_state%inc1%va_inc(i,j,k)*rdt
           IAU_Data%temp_inc(i,j,k)  =IAU_state%inc1%temp_inc(i,j,k)*rdt
           IAU_Data%delp_inc(i,j,k) =IAU_state%inc1%delp_inc(i,j,k)*rdt
+          IAU_Data%delz_inc(i,j,k) =IAU_state%inc1%delz_inc(i,j,k)*rdt
           do l = 1,ntracers
              IAU_Data%tracer_inc(i,j,k,l) =IAU_state%inc1%tracer_inc(i,j,k,l)*rdt
           enddo
@@ -382,6 +389,7 @@ subroutine read_iau_forcing(IPD_Control,increments,fname)
  ! read in 1 time level
     call interp_inc('T_inc',increments%temp_inc(:,:,:),jbeg,jend)
     call interp_inc('delp_inc',increments%delp_inc(:,:,:),jbeg,jend)
+    call interp_inc('delz_inc',increments%delz_inc(:,:,:),jbeg,jend)
     call interp_inc('u_inc',increments%ua_inc(:,:,:),jbeg,jend)   ! can these be treated as scalars?
     call interp_inc('v_inc',increments%va_inc(:,:,:),jbeg,jend)
     do l=1,ntracers

atmos_cubed_sphere/tools/fv_treat_da_inc.F90

@@ -172,6 +172,9 @@ subroutine read_da_inc(Atm, fv_domain)
 
     call apply_inc_on_3d_scalar('T_inc',Atm(1)%pt)
     call apply_inc_on_3d_scalar('delp_inc',Atm(1)%delp)
+    if (.not. Atm(1)%flagstruct%hydrostatic) then
+        call apply_inc_on_3d_scalar('delz_inc',Atm(1)%delz)
+    endif
     call apply_inc_on_3d_scalar('sphum_inc',Atm(1)%q(:,:,:,sphum))
     call apply_inc_on_3d_scalar('liq_wat_inc',Atm(1)%q(:,:,:,liq_wat))
     call apply_inc_on_3d_scalar('o3mr_inc',Atm(1)%q(:,:,:,o3mr))

gfsphysics/GFS_layer/GFS_diagnostics.F90

@@ -3354,6 +3354,24 @@ subroutine diag_populate (IPD_Diag, Model, Statein, Stateout, Sfcprop, Coupling,
       IPD_Diag(idx)%data(nb)%var3p => Diag(nb)%skebv_wts(:,:)
     enddo
 
+    !---dividing streamline
+    idx = idx + 1
+    IPD_Diag(idx)%name           = 'zmtnblck'
+    IPD_Diag(idx)%output_name    = 'zmtnblck'
+    IPD_Diag(idx)%mod_name       = 'physics'
+    IPD_Diag(idx)%file_name      = 'surface'
+    IPD_Diag(idx)%desc           = 'perturbation velocity'
+    IPD_Diag(idx)%unit           = 'm/s'
+    IPD_Diag(idx)%type_stat_proc = 'inst'
+    IPD_Diag(idx)%level_type     = ''
+    IPD_Diag(idx)%level          = 1
+    IPD_Diag(idx)%cnvfac         = cn_one
+    IPD_Diag(idx)%zhour          = Model%zhour
+    IPD_Diag(idx)%fcst_hour      = Model%fhour
+    do nb = 1,nblks
+      IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%zmtnblck(:)
+    enddo
+
     !---sppt  
     idx = idx + 1
     IPD_Diag(idx)%name           = 'sppt_wts'

gfsphysics/GFS_layer/GFS_driver.F90

@@ -361,13 +361,13 @@ subroutine GFS_time_vary_step (Model, Statein, Stateout, Sfcprop, Coupling, &
           enddo
        enddo
     endif
-    if (Model%do_sppt) then
-       do nb = 1,nblks
-          do k=1,Model%levs
-             Diag(nb)%sppt_wts(:,k)=Coupling(nb)%sppt_wts(:,Model%levs-k+1)
-          enddo
-       enddo
-    endif
+    !if (Model%do_sppt) then
+    !   do nb = 1,nblks
+    !      do k=1,Model%levs
+    !         Diag(nb)%sppt_wts(:,k)=Coupling(nb)%sppt_wts(:,Model%levs-k+1)
+    !      enddo
+    !   enddo
+    !endif
     if (Model%do_shum) then
        do nb = 1,nblks
           do k=1,Model%levs
@@ -407,14 +407,34 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
     type(GFS_diag_type),      intent(inout) :: Diag
     !--- local variables
     integer :: k, i
-    real(kind=kind_phys) :: upert, vpert, tpert, qpert, qnew
+    real(kind=kind_phys) :: upert, vpert, tpert, qpert, qnew,sppt_vwt
      if (Model%do_sppt) 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
+           if (Model%use_zmtnblck)then
+              Coupling%sppt_wts(i,k)=(Coupling%sppt_wts(i,k)-1)*sppt_vwt+1.0
+           endif
+           Diag%sppt_wts(i,Model%levs-k+1)=Coupling%sppt_wts(i,k)
 
            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))   * Coupling%sppt_wts(i,k) - Tbd%dtdtr(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)
 

gfsphysics/GFS_layer/GFS_physics_driver.F90

@@ -1510,7 +1510,8 @@ subroutine GFS_physics_driver                         &
                  Sfcprop%hprime(1,1), oc, oa4, clx, theta,  &
                  sigma, gamma, elvmax, dusfcg, dvsfcg,      &
                  con_g, con_cp, con_rd, con_rv, Model%lonr, &
-                 Model%nmtvr, Model%cdmbgwd, me, lprnt,ipr)
+                 Model%nmtvr, Model%cdmbgwd, me, lprnt,ipr, &
+                 Diag%zmtnblck)
 
 !     if (lprnt)  print *,' dudtg=',dudt(ipr,:)
 

gfsphysics/GFS_layer/GFS_typedefs.F90

@@ -555,6 +555,7 @@ module GFS_typedefs
 
     !--- stochastic physics control parameters
     logical              :: do_sppt
+    logical              :: use_zmtnblck
     logical              :: do_shum
     logical              :: do_skeb
     integer              :: skeb_npass 
@@ -842,10 +843,11 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: wet1   (:)    => null()   !< normalized soil wetness
     real (kind=kind_phys), pointer :: sr     (:)    => null()   !< snow ratio : ratio of snow to total precipitation
 
-    real (kind=kind_phys), pointer :: skebu_wts(:,:)    => null()   !< 10 meater u/v wind speed
-    real (kind=kind_phys), pointer :: skebv_wts(:,:)    => null()   !< 10 meater u/v wind speed
+    real (kind=kind_phys), pointer :: skebu_wts(:,:)   => null()   !< 10 meater u/v wind speed
+    real (kind=kind_phys), pointer :: skebv_wts(:,:)   => null()   !< 10 meater u/v wind speed
     real (kind=kind_phys), pointer :: sppt_wts(:,:)    => null()   !< 10 meater u/v wind speed
     real (kind=kind_phys), pointer :: shum_wts(:,:)    => null()   !< 10 meater u/v wind speed
+    real (kind=kind_phys), pointer :: zmtnblck(:)      => null()   !<mountain blocking evel
     !--- accumulated quantities for 3D diagnostics
     real (kind=kind_phys), pointer :: du3dt (:,:,:) => null()   !< u momentum change due to physics
     real (kind=kind_phys), pointer :: dv3dt (:,:,:) => null()   !< v momentum change due to physics
@@ -1558,9 +1560,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     logical              :: debug          = .false.
     logical              :: pre_rad        = .false.         !< flag for testing purpose
     !--- stochastic physics control parameters
-    logical :: do_sppt   = .false.
-    logical :: do_shum   = .false.
-    logical :: do_skeb   = .false.
+    logical :: do_sppt      = .false.
+    logical :: use_zmtnblck = .false.
+    logical :: do_shum      = .false.
+    logical :: do_skeb      = .false.
     integer :: skeb_npass = 11
     !--- END NAMELIST VARIABLES
 
@@ -1813,6 +1816,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
 
     !--- stochastic physics options
     Model%do_sppt          = do_sppt
+    Model%use_zmtnblck     = use_zmtnblck
     Model%do_shum          = do_shum
     Model%do_skeb          = do_skeb
 
@@ -2304,11 +2308,6 @@ subroutine control_print(Model)
       print *, ' xkzminv           : ', Model%xkzminv
       print *, ' moninq_fac        : ', Model%moninq_fac
       print *, ' '
-      print *, 'stochastic physics'
-      print *, ' do_sppt           : ', Model%do_sppt
-      print *, ' do_shum           : ', Model%do_shum
-      print *, ' do_skeb           : ', Model%do_skeb
-      print *, ' '
       print *, 'tracers'
       print *, ' tracer_names      : ', Model%tracer_names
       print *, ' ntrac             : ', Model%ntrac
@@ -2633,6 +2632,7 @@ subroutine diag_create (Diag, IM, Model)
     allocate (Diag%skebv_wts(IM,Model%levs))
     allocate (Diag%sppt_wts(IM,Model%levs))
     allocate (Diag%shum_wts(IM,Model%levs))
+    allocate (Diag%zmtnblck(IM))
     !--- 3D diagnostics
     if (Model%ldiag3d) then
       allocate (Diag%du3dt  (IM,Model%levs,4))

gfsphysics/physics/gwdps.f

@@ -184,7 +184,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL,              &
      &               PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DELTIM,KDT,         &
      &               HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX,       &
      &               DUSFC,DVSFC,G, CP, RD, RV, IMX,                    &
-     &               nmtvr, cdmbgwd, me, lprnt, ipr)
+     &               nmtvr, cdmbgwd, me, lprnt, ipr,  rdxzb)
 !
 !   ********************************************************************
 ! ----->  I M P L E M E N T A T I O N    V E R S I O N   <----------
@@ -286,7 +286,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL,              &
      &                     U1(IX,KM),   V1(IX,KM),     T1(IX,KM),       &
      &                     Q1(IX,KM),   PRSI(IX,KM+1), DEL(IX,KM),      &
      &                     PRSL(IX,KM), PRSLK(IX,KM),  PHIL(IX,KM),     &
-     &                     PHII(IX,KM+1)
+     &                     PHII(IX,KM+1),RDXZB(IY)
       real(kind=kind_phys) OC(IM),     OA4(IY,4), CLX4(IY,4)            &
      &,                    HPRIME(IM)
 ! for lm mtn blocking
@@ -408,6 +408,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL,              &
 !
       IF ( NMTVR .eq. 14) then 
 ! ----  for lm and gwd calculation points
+        RDXZB(:)  = 0 
         ipt = 0
         npt = 0
         DO I = 1,IM
@@ -585,7 +586,10 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL,              &
               EK(I)  = 0.5 *  UP(I) * UP(I) 
 
 ! --- Dividing Stream lime  is found when PE =exceeds EK.
-              IF ( PE(I) .ge.  EK(I) ) IDXZB(I) = K
+              IF ( PE(I) .ge.  EK(I) ) THEN
+                 IDXZB(I) = K
+                 RDXZB(J) = real(K,kind=kind_phys)
+              ENDIF
 ! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface
 !
 !> - The dividing streamline height (idxzb), of a subgrid scale 
@@ -709,6 +713,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL,              &
 !
         do i=1,npt
           IDXZB(i) = 0
+          RDXZB(i) = 0.
         enddo
       ENDIF
 !

gfsphysics/physics/h2ointerp.f90

@@ -126,7 +126,7 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy)
       integer  idat(8),jdat(8)
 !
       real(kind=kind_phys) ddy(npts)
-      real(kind=kind_phys) h2oplout(levh2o,npts,h2o_coeff)
+      real(kind=kind_phys) h2oplout(npts,levh2o,h2o_coeff)
       real(kind=kind_phys) rinc(5), rjday
       integer              jdow, jdoy, jday
       real(4)              rinc4(5)
@@ -152,35 +152,32 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy)
       jday = 0
       call w3doxdat(jdat,jdow,jdoy,jday)
       rjday = jdoy + jdat(5) / 24.
-      if (rjday < h2o_time(1)) rjday = rjday+365.
+      if (rjday < h2o_time(1)) rjday = rjday + 365.
 !
       n2 = timeh2o + 1
-      do j=1,timeh2o
+      do j=2,timeh2o
         if (rjday < h2o_time(j)) then
           n2 = j
           exit
         endif
       enddo
       n1 = n2 - 1
-      if (n1 <= 0)      n1 = n1 + timeh2o
-      if (n2 > timeh2o) n2 = n2 - timeh2o
-
 !
 !     if (me .eq. 0) print *,' n1=',n1,' n2=',n2,' rjday=',rjday
 !    &,'h2o_time=',h2o_time(n1),h2o_time(n2)
 !
-
       tx1 = (h2o_time(n2) - rjday) / (h2o_time(n2) - h2o_time(n1))
       tx2 = 1.0 - tx1
+      if (n2 > timeh2o) n2 = n2 - timeh2o
 !
       do nc=1,h2o_coeff
         do l=1,levh2o
           do j=1,npts
             j1  = jindx1(j)
             j2  = jindx2(j)
             tem = 1.0 - ddy(j)
-            h2oplout(j,l,nc) = & 
-              tx1*(tem*h2oplin(j1,l,nc,n1)+ddy(j)*h2oplin(j2,l,nc,n1)) & 
+            h2oplout(j,l,nc) =                                         &
+              tx1*(tem*h2oplin(j1,l,nc,n1)+ddy(j)*h2oplin(j2,l,nc,n1)) &
             + tx2*(tem*h2oplin(j1,l,nc,n2)+ddy(j)*h2oplin(j2,l,nc,n2))
           enddo
         enddo