Adding gust and precipitation rates

FV3GFS/FV3GFS_io.F90

@@ -4585,7 +4585,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     Diag(idx)%axes = 2
     Diag(idx)%name = 'rain'
     Diag(idx)%desc = 'total rain at this time step'
-    Diag(idx)%unit = 'XXX'
+    Diag(idx)%unit = 'm'
     Diag(idx)%mod_name = 'gfs_phys'
     allocate (Diag(idx)%data(nblks))
     do nb = 1,nblks
@@ -4596,13 +4596,35 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     Diag(idx)%axes = 2
     Diag(idx)%name = 'rainc'
     Diag(idx)%desc = 'convective rain at this time step'
-    Diag(idx)%unit = 'XXX'
+    Diag(idx)%unit = 'm'
     Diag(idx)%mod_name = 'gfs_phys'
     allocate (Diag(idx)%data(nblks))
     do nb = 1,nblks
       Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rainc(:)
     enddo
 
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = 'precip_rate'
+    Diag(idx)%desc = 'instantaneous precipitation rate at this time step'
+    Diag(idx)%unit = 'kg/m2/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%instant_precip_rate(:)
+    enddo
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = 'precipc_rate'
+    Diag(idx)%desc = 'instantaneous convective precipitation rate at this time step'
+    Diag(idx)%unit = 'kg/m2/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%instant_conv_precip_rate(:)
+    enddo
+
     idx = idx + 1
     Diag(idx)%axes = 2
     Diag(idx)%name = 'ice'
@@ -4742,6 +4764,30 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
       Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v10m(:)
     enddo
 
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = 'u100m'
+    Diag(idx)%desc = '100 meter u wind [m/s]'
+    Diag(idx)%unit = 'm/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'vector_bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%u100m(:)
+    enddo
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = 'v100m'
+    Diag(idx)%desc = '100 meter v wind [m/s]'
+    Diag(idx)%unit = 'm/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'vector_bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v100m(:)
+    enddo
+
     !Sofar added: start: (12/20/23)
     idx = idx + 1
     Diag(idx)%axes = 2
@@ -4766,6 +4812,57 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     do nb = 1,nblks
       Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v10n(:)
     enddo
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = '10mgust'
+    Diag(idx)%desc = '10-meter wind gust [m/s]'
+    Diag(idx)%unit = 'm/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gust10m(:)
+    enddo
+
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = '10mgustmax'
+    Diag(idx)%desc = 'maximum 10-meter wind gust [m/s]'
+    Diag(idx)%unit = 'm/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gustmax10m(:)
+    enddo
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = '100mgust'
+    Diag(idx)%desc = '100-meter wind gust [m/s]'
+    Diag(idx)%unit = 'm/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gust100m(:)
+    enddo
+
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = '100mgustmax'
+    Diag(idx)%desc = 'maximum 100-meter wind gust [m/s]'
+    Diag(idx)%unit = 'm/s'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gustmax100m(:)
+    enddo
+
     !Sofar added: end
 
     idx = idx + 1
@@ -4792,6 +4889,31 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
       Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%zol(:)
     enddo
 
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = 'cd'
+    Diag(idx)%desc = 'Momentum drag coefficient'
+    Diag(idx)%unit = '1'
+    Diag(idx)%mod_name = 'gfs_phys'
+    Diag(idx)%intpl_method = 'bilinear'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%cd(:)
+    enddo
+
+
+    idx = idx + 1
+    Diag(idx)%axes = 2
+    Diag(idx)%name = 'fm10'
+    Diag(idx)%desc = 'fm10 parameter from SL scheme at 10 meters'
+    Diag(idx)%unit = '1'
+    Diag(idx)%mod_name = 'gfs_phys'
+    allocate (Diag(idx)%data(nblks))
+    do nb = 1,nblks
+      Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%fm10(:)
+    enddo
+
     idx = idx + 1
     Diag(idx)%axes = 2
     Diag(idx)%name = 'dpt2m'
@@ -7000,7 +7122,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     Diag(idx)%axes = 2
     Diag(idx)%name = 'f10m'
     Diag(idx)%desc = '10-meter wind speed divided by lowest model wind speed'
-    Diag(idx)%unit = 'N/A'
+    Diag(idx)%unit = '1'
     Diag(idx)%mod_name = 'gfs_sfc'
     allocate (Diag(idx)%data(nblks))
     do nb = 1,nblks
@@ -7033,7 +7155,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     Diag(idx)%axes = 2
     Diag(idx)%name = 'ffhh'
     Diag(idx)%desc = 'fh parameter from PBL scheme'
-    Diag(idx)%unit = 'XXX'
+    Diag(idx)%unit = '1'
     Diag(idx)%mod_name = 'gfs_sfc'
     allocate (Diag(idx)%data(nblks))
     do nb = 1,nblks
@@ -7043,8 +7165,8 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     idx = idx + 1
     Diag(idx)%axes = 2
     Diag(idx)%name = 'ffmm'
-    Diag(idx)%desc = 'fm parameter from PBL scheme'
-    Diag(idx)%unit = 'XXX'
+    Diag(idx)%desc = 'fm parameter from SL scheme'
+    Diag(idx)%unit = '1'
     Diag(idx)%mod_name = 'gfs_sfc'
     allocate (Diag(idx)%data(nblks))
     do nb = 1,nblks
@@ -7055,7 +7177,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block
     Diag(idx)%axes = 2
     Diag(idx)%name = 'uustar'
     Diag(idx)%desc = 'uustar surface frictional wind'
-    Diag(idx)%unit = 'XXX'
+    Diag(idx)%unit = 'm/s'
     Diag(idx)%mod_name = 'gfs_sfc'
     allocate (Diag(idx)%data(nblks))
     do nb = 1,nblks

GFS_layer/GFS_diagnostics.F90

@@ -1480,6 +1480,43 @@ subroutine diag_populate (IPD_Diag, Model, Statein, Stateout, Sfcprop, Coupling,
       IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%wind10mmax
     enddo
 
+
+    !---wind10mmax
+    idx = idx + 1
+    IPD_Diag(idx)%name           = 'gust10mmaxn'
+    IPD_Diag(idx)%output_name    = 'gust10mmax'
+    IPD_Diag(idx)%mod_name       = 'physics'
+    IPD_Diag(idx)%file_name      = 'flx'
+    IPD_Diag(idx)%desc           = 'maximum wind gust (m/s) at 10 m above ground'
+    IPD_Diag(idx)%unit           = 'm/s'
+    IPD_Diag(idx)%type_stat_proc = 'max'
+    IPD_Diag(idx)%level_type     = '10 m above grnd'
+    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)%gustmax10m
+    enddo
+
+    !---wind100mmax
+    idx = idx + 1
+    IPD_Diag(idx)%name           = 'gust100mmaxn'
+    IPD_Diag(idx)%output_name    = 'gust100mmax'
+    IPD_Diag(idx)%mod_name       = 'physics'
+    IPD_Diag(idx)%file_name      = 'flx'
+    IPD_Diag(idx)%desc           = 'maximum wind gust (m/s) at 100 m above ground'
+    IPD_Diag(idx)%unit           = 'm/s'
+    IPD_Diag(idx)%type_stat_proc = 'max'
+    IPD_Diag(idx)%level_type     = '100 m above grnd'
+    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)%gustmax100m
+    enddo
+
     !---rain
     idx = idx + 1
     IPD_Diag(idx)%name           = 'rain'

GFS_layer/GFS_physics_driver.F90

@@ -1196,28 +1196,29 @@ subroutine GFS_physics_driver                         &
 ! a new and more flexible version of sfc_diff by kgao
             call sfc_diff_gfdl(im,Statein%pgr, Statein%ugrs, Statein%vgrs,&
                  Statein%tgrs, Statein%qgrs, Diag%zlvl, Sfcprop%snowd, &
-                 Sfcprop%tsfc, Sfcprop%zorl, Sfcprop%ztrl, cd,      &
+                 Sfcprop%tsfc, Sfcprop%qsfc, Sfcprop%zorl, Sfcprop%ztrl, cd,      &
                  cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, &
                  Sfcprop%ffmm,  Sfcprop%ffhh,                       &
                  Sfcprop%charnock,                                  &              
                  fm10_neutral,                                      &  
                  Sfcprop%uustar,                                    &
-                 wind,  Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2,    &
+                 wind,  Tbd%phy_f2d(1,Model%num_p2d), Diag%fm10, fh2,    &
                  sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc,    &
                  tsurf, flag_iter, Model%redrag, Model%z0s_max,     &
                  Model%do_z0_moon, Model%do_z0_hwrf15,              &
                  Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly,     &
                  Model%wind_th_hwrf,                                &
                  Model%alpha_stable, Model%alpha_unstable,          &
-                 Model%tune_ocean_surface_layer, Diag%zol)
+                 Model%tune_ocean_surface_layer, Diag%zol,          & 
+                 Model%add_w_freeconv)
             else
 ! GFS original sfc_diff modified by kgao 
             call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,&
                  Statein%tgrs, Statein%qgrs, Diag%zlvl,             &
                  Sfcprop%snowd, Sfcprop%tsfc,  Sfcprop%zorl, cd,    &
                  cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, &
                  Sfcprop%ffmm,  Sfcprop%ffhh, Sfcprop%uustar,       &
-                 wind,  Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2,    &
+                 wind,  Tbd%phy_f2d(1,Model%num_p2d), Diag%fm10, fh2,    &
                  sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc,    &
                  tsurf, flag_iter, Model%redrag, Model%czil_sfc,    &
                  Model%z0s_max,     &
@@ -1226,6 +1227,7 @@ subroutine GFS_physics_driver                         &
                  Model%wind_th_hwrf, Diag%zol)
             endif
             Diag%rb = rb
+            Diag%cd = cd
          !endif
 
          !endif
@@ -1484,11 +1486,31 @@ subroutine GFS_physics_driver                         &
               Statein%tgrs, Statein%qgrs, Sfcprop%tsfc, qss,   &
               Sfcprop%f10m, Diag%u10m, Diag%v10m,        &
               Sfcprop%t2m, Sfcprop%q2m, work3, evap,           &
-              Sfcprop%ffmm, Sfcprop%ffhh, fm10, fh2, &
+              Sfcprop%ffmm, Sfcprop%ffhh, Diag%fm10, fh2, &
               Sfcprop%u10m, Sfcprop%v10m,&
               fm10_neutral, Diag%u10n, Diag%v10n, &
               Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa)  
       !endif
+      ! get 100-m wind components using linear interpolation same as in FV3, then get wind gust at 10 and 100 m
+      call interpolate_z(im, levs, 100., Statein%phii, Statein%ugrs, Diag%u100m)
+      call interpolate_z(im, levs, 100., Statein%phii, Statein%vgrs, Diag%v100m)
+      call compute_gust(im, Diag%u10m, Diag%v10m, Sfcprop%uustar, &
+            Diag%zol, Diag%zlvl, Model%gust_parameter, &
+            Diag%gust10m)
+      call compute_gust(im, Diag%u100m, Diag%v100m, Sfcprop%uustar, &
+            Diag%zol, Diag%zlvl, Model%gust_parameter, &
+            Diag%gust100m)
+      do i=1, im
+           !find max wind gusts
+           tem = Diag%gust10m(i)
+           if (tem > Diag%gustmax10m(i)) then
+              Diag%gustmax10m(i) = tem
+           endif
+           tem = Diag%gust100m(i)
+           if (tem > Diag%gustmax100m(i)) then
+              Diag%gustmax100m(i) = tem
+           endif
+      end do
 
       Tbd%phy_f2d(:,Model%num_p2d) = 0.0
 
@@ -3850,7 +3872,7 @@ subroutine GFS_physics_driver                         &
                        Stateout%gt0, Stateout%gq0, Sfcprop%tsfc, qss,   &
                        Sfcprop%f10m, Diag%u10m, Diag%v10m, Sfcprop%t2m, &
                        Sfcprop%q2m, work3, evap, Sfcprop%ffmm,          &
-                       Sfcprop%ffhh, fm10, fh2, &
+                       Sfcprop%ffhh, Diag%fm10, fh2, &
                        Sfcprop%u10m, Sfcprop%v10m,&
                        fm10_neutral, Diag%u10n, Diag%v10n, &
                        Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa)  
@@ -3949,6 +3971,11 @@ subroutine GFS_physics_driver                         &
         enddo
       endif
 
+      do i = 1, im
+        Diag%instant_precip_rate(i) = Diag%rain(i) / (dtf*con_p001)
+        Diag%instant_conv_precip_rate(i) = Diag%rainc(i) / (dtf*con_p001)
+      end do
+
       deallocate (clw)
       deallocate (clw_trac_idx)
       if (Model%do_shoc) then
@@ -4322,6 +4349,46 @@ subroutine compute_diagnostics_with_scaled_co2(Model, Statein, Sfcprop, Coupling
       enddo
 
   end subroutine compute_diagnostics_with_scaled_co2
+
+  subroutine interpolate_z(im, km, zl, phi, a3_in, a2)
+    use physcons, only: con_g
+    implicit none
+    integer,  intent(in):: im, km
+    real, intent(in):: phi(im,km+1) ! phi(k+1) > phi(k)
+    real, intent(in):: a3_in(im,km) 
+    real, intent(in):: zl
+    real, intent(out):: a2(im)
+    ! local:
+    real hght(im,km+1), a3(im, km) ! hght(k) > hght(k+1)
+    real zm(km)
+    integer i,k
+
+    ! flip z and a3
+    do k=1,km+1
+      do i=1,im
+        hght(i,k) = phi(i,km+2-k) / con_g
+        if (k <= km) a3(i,k) = a3_in(i,km+1-k)
+      enddo
+    enddo
+
+    do i=1,im
+      do k=1,km
+        zm(k) = 0.5*(hght(i,k)+hght(i,k+1))
+      enddo
+      if( zl >= zm(1) ) then
+        a2(i) = a3(i,1)
+      elseif ( zl <= zm(km) ) then
+        a2(i) = a3(i,km)
+      else
+        do k=1,km-1
+          if( zl <= zm(k) .and. zl >= zm(k+1) ) then
+            a2(i) = a3(i,k) + (a3(i,k+1)-a3(i,k))*(zm(k)-zl)/(zm(k)-zm(k+1))
+            exit
+          endif
+        enddo
+      endif
+    end do
+ end subroutine interpolate_z
 !> @}
 
 end module module_physics_driver