Merge pull request NCAR#360 from climbfuji/dtc_develop_udpate_from_em…

…c_20191127

dtc/develop: update from EMC 2019/11/27

physics/GFS_surface_composites.F90

@@ -123,7 +123,7 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan
             if (cice(i) < one) then
               wet(i) = .true.
   !           tsfco(i) = tgice
-              tsfco(i) = max(tisfc(i), tgice)
+              if (.not. cplflx) tsfco(i) = max(tisfc(i), tgice)
   !           tsfco(i) = max((tsfc(i) - cice(i)*tisfc(i)) &
   !                                     / (one - cice(i)), tgice)
             endif

physics/drag_suite.F90

@@ -332,11 +332,11 @@ subroutine drag_suite_run(                                           &
      &                                     hpbl(im),             &
      &                                     slmsk(im)
    real(kind=kind_phys), dimension(im)    :: govrth,xland
-   real(kind=kind_phys), dimension(im,km) :: dz2
+   !real(kind=kind_phys), dimension(im,km) :: dz2
    real(kind=kind_phys)                   :: tauwavex0,tauwavey0,  &
      &                                     XNBV,density,tvcon,hpbl2
    integer                          ::     kpbl2,kvar
-   real(kind=kind_phys), dimension(im,km+1)         ::     zq      ! = PHII/g
+   !real(kind=kind_phys), dimension(im,km+1)         ::     zq      ! = PHII/g
    real(kind=kind_phys), dimension(im,km)           ::     zl      ! = PHIL/g
 
 !SPP
@@ -413,10 +413,10 @@ subroutine drag_suite_run(                                           &
 !  local variables
 !
    integer              ::  i,j,k,lcap,lcapp1,nwd,idir,           &
-                            klcap,kp1,ikount,kk
+                            klcap,kp1
 !
-   real(kind=kind_phys) ::  rcs,rclcs,csg,fdir,cleff,cleff_ss,cs, &
-                            rcsks,wdir,ti,rdz,temp,tem2,dw2,shr2, &
+   real(kind=kind_phys) ::  rcs,csg,fdir,cleff,cleff_ss,cs,       &
+                            rcsks,wdir,ti,rdz,tem2,dw2,shr2,      &
                             bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro,  &
                             rim,temc,tem1,efact,temv,dtaux,dtauy, &
                             dtauxb,dtauyb,eng0,eng1
@@ -442,7 +442,6 @@ subroutine drag_suite_run(                                           &
                             coefm(im),coefm_ss(im)
 !
    integer              ::  kbl(im),klowtop(im)
-   logical :: iope
    integer,parameter    ::  mdir=8
    !integer              ::  nwdir(mdir)
    !data nwdir/6,7,5,8,2,3,1,4/
@@ -661,6 +660,17 @@ subroutine drag_suite_run(                                           &
      enddo
    enddo
 !
+!  calculate mid-layer height (zl), interface height (zq), and layer depth (dz2).
+!
+   !zq=0.
+   do k = kts,km
+     do i = its,im
+       !zq(i,k+1) = PHII(i,k+1)*g_inv
+       !dz2(i,k)  = (PHII(i,k+1)-PHII(i,k))*g_inv
+       zl(i,k)   = PHIL(i,k)*g_inv
+     enddo
+   enddo
+!
 !  determine reference level: maximum of 2*var and pbl heights
 !
    do i = its,im
@@ -895,7 +905,6 @@ subroutine drag_suite_run(                                           &
   density=1.2
   utendwave=0.
   vtendwave=0.
-  zq=0.
 !
   IF ( (gwd_opt_ss .EQ. 1).and.(ss_taper.GT.1.E-02) ) THEN
     if (me==master) print *,"in Drag Suite: Running small-scale gravity wave drag"
@@ -914,14 +923,6 @@ subroutine drag_suite_run(                                           &
         thvx(i,k) = thx(i,k)*tvcon
       enddo
     enddo
-    ! Calculate mid-layer height (zl), interface height (zq), and layer depth (dz2).
-    do k = kts,km
-      do i = its,im
-        zq(i,k+1) = PHII(i,k+1)*g_inv
-        dz2(i,k)  = (PHII(i,k+1)-PHII(i,k))*g_inv
-        zl(i,k)   = PHIL(i,k)*g_inv
-      enddo
-    enddo
 
     do i=its,im
        hpbl2 = hpbl(i)+10.
@@ -1027,19 +1028,6 @@ subroutine drag_suite_run(                                           &
 
    utendform=0.
    vtendform=0.
-   zq=0.
-
-   IF ( (gwd_opt_ss .NE. 1).and.(ss_taper.GT.1.E-02) ) THEN
-     ! Defining mid-layer height (zl), interface height (zq), and layer depth (dz2). 
-     ! This is already done above if the small-scale GWD is activated.
-     do k = kts,km
-       do i = its,im
-         zq(i,k+1) = PHII(i,k+1)*g_inv
-         dz2(i,k)  = (PHII(i,k+1)-PHII(i,k))*g_inv
-         zl(i,k)   = PHIL(i,k)*g_inv
-       enddo
-     enddo
-   ENDIF
 
    DO i=its,im
       IF ((xland(i)-1.5) .le. 0.) then

physics/module_gfdl_cloud_microphys.F90

@@ -4729,7 +4729,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
     real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
     real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
     real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
-    real :: qmin = 1.0e-12, beta = 1.22
+    real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6
 
     do k = ks, ke
         do i = is, ie
@@ -4759,7 +4759,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
                 ! cloud ice (Heymsfield and Mcfarquhar, 1996)
             ! -----------------------------------------------------------------------
 
-                if (qmi (i, k) .gt. qmin) then
+                if (qmi (i, k) .gt. qmin1) then
                     qci (i, k) = dpg * qmi (i, k) * 1.0e3
                     rei_fac = log (1.0e3 * qmi (i, k) * den (i, k))
                     if (t (i, k) - tice .lt. - 50) then
@@ -4785,7 +4785,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
                 ! cloud ice (Wyser, 1998)
             ! -----------------------------------------------------------------------
 
-                if (qmi (i, k) .gt. qmin) then
+                if (qmi (i, k) .gt. qmin1) then
                     qci (i, k) = dpg * qmi (i, k) * 1.0e3
                     bw = - 2. + 1.e-3 * log10 (den (i, k) * qmi (i, k) / rho_0) * max (0.0, tice - t (i, k)) ** 1.5
                     rei (i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
@@ -4815,7 +4815,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
             ! snow (Lin et al., 1983)
             ! -----------------------------------------------------------------------
 
-            if (qms (i, k) .gt. qmin) then
+            if (qms (i, k) .gt. qmin1) then
                 qcs (i, k) = dpg * qms (i, k) * 1.0e3
                 lambdas = exp (0.25 * log (pi * rhos * n0s / qms (i, k) / den (i, k)))
                 res (i, k) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6

physics/module_mp_thompson.F90

@@ -1302,7 +1302,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_qc = k
              qc_max = qc1d(k)
             elseif (qc1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative qc ', qc1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qc ', qc1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (qr1d(k) .gt. qr_max) then
@@ -1311,7 +1311,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_qr = k
              qr_max = qr1d(k)
             elseif (qr1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative qr ', qr1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qr ', qr1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (nr1d(k) .gt. nr_max) then
@@ -1320,7 +1320,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_nr = k
              nr_max = nr1d(k)
             elseif (nr1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative nr ', nr1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative nr ', nr1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (qs1d(k) .gt. qs_max) then
@@ -1329,7 +1329,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_qs = k
              qs_max = qs1d(k)
             elseif (qs1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative qs ', qs1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qs ', qs1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (qi1d(k) .gt. qi_max) then
@@ -1338,7 +1338,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_qi = k
              qi_max = qi1d(k)
             elseif (qi1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative qi ', qi1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qi ', qi1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (qg1d(k) .gt. qg_max) then
@@ -1347,7 +1347,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_qg = k
              qg_max = qg1d(k)
             elseif (qg1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative qg ', qg1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qg ', qg1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (ni1d(k) .gt. ni_max) then
@@ -1356,11 +1356,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
              kmax_ni = k
              ni_max = ni1d(k)
             elseif (ni1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative ni ', ni1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative ni ', ni1d(k),        &
                         ' at i,j,k=', i,j,k
             endif
             if (qv1d(k) .lt. 0.0) then
-             write(*,*) 'WARNING, negative qv ', qv1d(k),        &
+             write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qv ', qv1d(k),        &
                         ' at i,j,k=', i,j,k
              if (k.lt.kte-2 .and. k.gt.kts+1) then
                 write(*,*) '   below and above are: ', qv(i,k-1,j), qv(i,k+1,j)

physics/module_mp_thompson_make_number_concentrations.F90

@@ -79,6 +79,11 @@ elemental real function make_IceNumber (Q_ice, temp)
          161.503, 168.262, 175.248, 182.473, 189.952, 197.699,          &
          205.728, 214.055, 222.694, 231.661, 240.971, 250.639 /)
 
+      if (Q_ice == 0) then
+         make_IceNumber = 0
+         return
+      end if
+
 !+---+-----------------------------------------------------------------+ 
 !..From the model 3D temperature field, subtract 179K for which
 !.. index value of retab as a start.  Value of corr is for
@@ -133,6 +138,11 @@ elemental real function make_DropletNumber (Q_cloud, qnwfa)
       real:: q_nwfa, x1, xDc
       integer:: nu_c
 
+      if (Q_cloud == 0) then
+         make_DropletNumber = 0
+         return
+      end if
+
 !+---+
 
       q_nwfa = MAX(99.E6, MIN(qnwfa,5.E10))
@@ -160,6 +170,11 @@ elemental real function make_RainNumber (Q_rain, temp)
       !real, parameter:: PI = 3.1415926536
       real, parameter:: am_r = PI*1000./6.
 
+      if (Q_rain == 0) then
+         make_RainNumber = 0
+         return
+      end if
+
       !+---+-----------------------------------------------------------------+ 
       !.. Not thrilled with it, but set Y-intercept parameter to Marshal-Palmer value
       !.. that basically assumes melting snow becomes typical rain. However, for
@@ -172,7 +187,7 @@ elemental real function make_RainNumber (Q_rain, temp)
       N0 = 8.E6
 
       if (temp .le. 271.15) then
-         N0 = 8.E8      
+         N0 = 8.E8
       elseif (temp .gt. 271.15 .and. temp.lt.273.15) then
          N0 = 8. * 10**(279.15-temp)
       endif

physics/mp_thompson_post.F90

@@ -67,7 +67,7 @@ end subroutine mp_thompson_post_init
 !!
 #endif
    subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, &
-                                        mpicomm, mpirank, mpiroot, errmsg, errflg)
+                                   kdt, mpicomm, mpirank, mpiroot, errmsg, errflg)
 
       implicit none
 
@@ -78,6 +78,7 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, &
       real(kind_phys), dimension(1:ncol,1:nlev), intent(inout) :: tgrs
       real(kind_phys), dimension(1:ncol,1:nlev), intent(in)    :: prslk
       real(kind_phys),                           intent(in)    :: dtp
+      integer,                                   intent(in)    :: kdt
       ! MPI information
       integer,          intent(in   ) :: mpicomm
       integer,          intent(in   ) :: mpirank
@@ -115,8 +116,8 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, &
 
             if (tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) .ne. tgrs(i,k)) then
 #ifdef DEBUG
-              write(0,*) "mp_thompson_post_run mp_tend limiter: i, k, t_old, t_new, t_lim:", &
-                         & i, k, tgrs_save(i,k), tgrs(i,k), tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k)
+              write(0,'(a,3i6,3e16.7)') "mp_thompson_post_run mp_tend limiter: kdt, i, k, t_old, t_new, t_lim:", &
+                                      & kdt, i, k, tgrs_save(i,k), tgrs(i,k), tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k)
 #endif
               events = events + 1
             end if
@@ -125,7 +126,8 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, &
       end do
 
       if (events > 0) then
-        write(0,'(a,i0,a,i0,a)') "mp_thompson_post_run: mp_tend_lim applied ", events, "/", nlev*ncol, " times"
+        write(0,'(a,i0,a,i0,a,i0)') "mp_thompson_post_run: mp_tend_lim applied ", events, "/", nlev*ncol, &
+                                  & " times at timestep ", kdt
       end if
 
    end subroutine mp_thompson_post_run

physics/mp_thompson_post.meta

@@ -92,6 +92,14 @@
   kind = kind_phys
   intent = in
   optional = F
+[kdt]
+  standard_name = index_of_time_step
+  long_name = current forecast iteration
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+  optional = F
 [mpicomm]
   standard_name = mpi_comm
   long_name = MPI communicator