Added code to perform hail size forecasting diagnostic

CMakeLists.txt

@@ -96,6 +96,7 @@ list(APPEND tools_srcs
   tools/fv_surf_map.F90
   tools/fv_timing.F90
   tools/init_hydro.F90
+  tools/module_diag_hailcast.F90
   tools/sim_nc_mod.F90
   tools/sorted_index.F90
   tools/test_cases.F90)

driver/fvGFS/fv_nggps_diag.F90

@@ -62,7 +62,7 @@ module fv_nggps_diags_mod
 !   </tr>
 ! </table>
 
- use mpp_mod,            only: mpp_pe, mpp_root_pe,FATAL,mpp_error
+ use mpp_mod,            only: mpp_pe, mpp_root_pe,FATAL,mpp_error, input_nml_file, stdlog
  use constants_mod,      only: grav, rdgas
  use time_manager_mod,   only: time_type, get_time
  use diag_manager_mod,   only: register_diag_field, send_data
@@ -75,6 +75,22 @@ module fv_nggps_diags_mod
                                max_uh,max_vorticity,bunkers_vector,       &
                                helicity_relative_CAPS,max_vorticity_hy1
  use fv_arrays_mod,      only: fv_atmos_type
+ use module_diag_hailcast, only: do_hailcast, id_hailcast_dhail,         &
+                                 id_hailcast_dhail1, id_hailcast_dhail2, &
+                                 id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5, &
+                                 id_hailcast_wdur, id_hailcast_wup_mask, &
+                                 hailcast_dhail1, hailcast_dhail1_max,   &
+                                 hailcast_dhail2, hailcast_dhail2_max,   &
+                                 hailcast_dhail3, hailcast_dhail3_max,   &
+                                 hailcast_dhail4, hailcast_dhail4_max,   &
+                                 hailcast_dhail5, hailcast_dhail5_max,   &
+                                 hailcast_wdur, hailcast_wup_mask, hailcast_dhail_max, &
+                                 kstt_hc, kend_hc,                       &
+                                 kstt_hc1,kstt_hc2,kstt_hc3,kstt_hc4,kstt_hc5, &
+                                 kend_hc1,kend_hc2,kend_hc3,kend_hc4,kend_hc5, &
+                                 kstt_hcd, kend_hcd, kstt_hcm, kend_hcm,       &
+                                 hailcast_init, hailcast_compute
+ use fms_mod,            only: check_nml_error
  use mpp_domains_mod,    only: domain1d, domainUG
 #ifdef MULTI_GASES
  use multi_gases_mod,  only:  virq
@@ -123,6 +139,7 @@ module fv_nggps_diags_mod
 
 ! file name
  character(len=64) :: file_name = 'gfs_dyn'
+ INTEGER :: istatus
 
 ! tracers
  character(len=128)   :: tname
@@ -136,6 +153,7 @@ module fv_nggps_diags_mod
  real, dimension(:,:),allocatable :: up2,dn2,uhmax03,uhmin03
  real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01
  real, dimension(:,:),allocatable :: maxvorthy1,maxvort02
+
  public :: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 #ifdef use_WRTCOMP
  public :: fv_dyn_bundle_setup
@@ -149,6 +167,20 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
     type(time_type), intent(in) :: Time
     integer :: n, i, j, nz
 
+    namelist /fv_diagnostics_nml/ do_hailcast
+    integer :: ios, ierr
+    integer :: unit
+
+    !namelist file for hailcast
+
+    ! Read Main namelist
+    read (input_nml_file,fv_diagnostics_nml,iostat=ios)
+    ierr = check_nml_error(ios,'fv_diagnostics_nml')
+
+    unit = stdlog()
+    write(unit, nml=fv_diagnostics_nml)
+    !end hailcast nml
+
     n = 1
     ncnsto = Atm(1)%ncnst
     npzo   = Atm(1)%npz
@@ -427,6 +459,13 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
          enddo
 !        if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,lat=',lat(isco,jsco),lat(ieco-2:ieco,jeco-2:jeco)*180./3.14157
        endif
+
+    endif
+
+    if (do_hailcast) then
+        !initialize hailcast
+        call hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco,      &
+                isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus)
     endif
 !
 !------------------------------------
@@ -607,7 +646,7 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
            do i=isco,ieco
              wk(i,j,k) = -Atm(n)%delp(i,j,k)/(Atm(n)%delz(i,j,k)*grav)*rdgas*Atm(n)%pt(i,j,k)
 #ifdef MULTI_GASES
-             wk(i,j,k) = wk(i,j,k)*virq(Atm(n)%q(i,j,k,:)
+             wk(i,j,k) = wk(i,j,k)*virq(Atm(n)%q(i,j,k,:))
 #else
              wk(i,j,k) = wk(i,j,k)*(1.+zvir*Atm(n)%q(i,j,k,sphum))
 #endif
@@ -734,6 +773,49 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
       call store_data(id_uhmin25, uhmin25, Time, kstt_uhmin25, kend_uhmin25)
     endif
 
+    IF ( .not. Atm(n)%flagstruct%hydrostatic .and. do_hailcast ) THEN
+        !--- max hourly hailcast hail diameter
+        if (id_hailcast_dhail > 0) then
+          call store_data(id_hailcast_dhail, hailcast_dhail_max, Time, kstt_hc, kend_hc)
+        endif
+
+        !--- max hourly hailcast hail diameter (embryo 1)
+        if ( id_hailcast_dhail1 > 0) then
+          call store_data(id_hailcast_dhail1, hailcast_dhail1_max, Time, kstt_hc1, kend_hc1)
+        endif
+
+        !--- max hourly hailcast hail diameter (embryo 2)
+        if ( id_hailcast_dhail2 > 0) then
+          call store_data(id_hailcast_dhail2, hailcast_dhail2_max, Time, kstt_hc2, kend_hc2)
+        endif
+
+        !--- max hourly hailcast hail diameter (embryo 3)
+        if ( id_hailcast_dhail3 > 0) then
+          call store_data(id_hailcast_dhail3, hailcast_dhail3_max, Time, kstt_hc3, kend_hc3)
+        endif
+
+        !--- max hourly hailcast hail diameter (embryo 4)
+        if ( id_hailcast_dhail4 > 0) then
+          call store_data(id_hailcast_dhail4, hailcast_dhail4_max, Time, kstt_hc4, kend_hc4)
+        endif
+
+        !--- max hourly hailcast hail diameter (embryo 5)
+        if ( id_hailcast_dhail5 > 0) then
+          call store_data(id_hailcast_dhail5, hailcast_dhail5_max, Time, kstt_hc5, kend_hc5)
+        endif
+
+        !--- hailcast updraft duration
+        if ( id_hailcast_wdur > 0) then
+          call store_data(id_hailcast_wdur, hailcast_wdur(isco:ieco, jsco:jeco), Time, kstt_hcd, kend_hcd)
+        endif
+        !--- hailcast updraft mask
+        if ( id_hailcast_wup_mask > 0) then
+          call store_data(id_hailcast_wup_mask, hailcast_wup_mask(isco:ieco, jsco:jeco), Time, kstt_hcm, kend_hcm)
+        endif
+    END IF
+
+    !call nullify_domain()
+
  end subroutine fv_nggps_diag
 
  subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
@@ -770,17 +852,17 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
          kdtt=0
        endif
        nsteps_per_reset = nint(avg_max_length/first_time)
-       do j=jsco,jeco
-          do i=isco,ieco
-             if(mod(kdtt,nsteps_per_reset)==0)then
-                up2(i,j) = -999.
-                dn2(i,j) = 999.
-                maxvorthy1(i,j)= 0.
-                maxvort01(i,j)= 0.
-                maxvort02(i,j)= 0.
-             endif
-          enddo
-        enddo
+       if(mod(kdtt,nsteps_per_reset)==0)then
+            do j=jsco,jeco
+              do i=isco,ieco
+                    up2(i,j) = -999.
+                    dn2(i,j) = 999.
+                    maxvorthy1(i,j)= 0.
+                    maxvort01(i,j)= 0.
+                    maxvort02(i,j)= 0.
+              enddo
+            enddo
+       endif
          call get_vorticity(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo, &
                            npzo,Atm(n)%u,Atm(n)%v,vort,    &
                            Atm(n)%gridstruct%dx,Atm(n)%gridstruct%dy,&
@@ -798,16 +880,16 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
                            vort,maxvort02,0., 2.e3)
          if( .not.Atm(n)%flagstruct%hydrostatic ) then
             call max_vv(isco,ieco,jsco,jeco,npzo,ngc,up2,dn2,Atm(n)%pe,Atm(n)%w)
-            do j=jsco,jeco
-               do i=isco,ieco
-                  if(mod(kdtt,nsteps_per_reset)==0)then
-                    uhmax03(i,j)= 0.
-                    uhmin03(i,j)= 0.
-                    uhmax25(i,j)= 0.
-                    uhmin25(i,j)= 0.
-                  endif
-               enddo
-             enddo
+            if(mod(kdtt,nsteps_per_reset)==0)then
+                do j=jsco,jeco
+                  do i=isco,ieco
+                       uhmax03(i,j)= 0.
+                       uhmin03(i,j)= 0.
+                       uhmax25(i,j)= 0.
+                       uhmin25(i,j)= 0.
+                  enddo
+                enddo
+            endif
 
              call max_uh(isco,ieco,jsco,jeco,ngc,npzo,zvir, &
                            sphum,uhmax03,uhmin03,Atm(n)%w,vort,Atm(n)%delz, &
@@ -820,8 +902,8 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
                            Atm(n)%pt,Atm(n)%peln,Atm(n)%phis,grav, &
                            2.e3, 5.e3)
          endif
-    kdtt=kdtt+1
-    deallocate (vort)
+         kdtt=kdtt+1
+         deallocate (vort)
     else
        print *,'Missing max/min hourly field in diag_table'
        print *,'Make sure the following are listed in the diag_table under gfs_dyn:'
@@ -830,6 +912,14 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
        stop
     endif
    endif
+
+   !allocate hailcast met field arrays
+   if (do_hailcast) then
+        call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, &
+                isco,jsco,ieco,jeco,isdo,iedo,jsdo,jedo,npzo,           &
+                Time_step_atmos,avg_max_length)
+   endif
+
  end subroutine fv_nggps_tavg
 !
  subroutine store_data(id, work, Time, nstt, nend)
@@ -1336,6 +1426,72 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
      if(rc==0)  num_field_dyn=num_field_dyn+1
    endif
 
+   if( .not.hydrostatico .and. do_hailcast) then
+        if( id_hailcast_dhail > 0 ) then
+            call find_outputname(trim(file_name),'hailcast_dhail_max',output_name)
+            if(mpp_pe()==mpp_root_pe())print *,'max hourly hailcast hail diameter, output name=',trim(output_name)
+            call add_field_to_bundle(trim(output_name),'Max hourly hailcast hail diameter', 'mm', "time: point",   &
+                 axes(1:2), fcst_grid, kstt_hc,kend_hc, dyn_bundle, output_file, rcd=rc)
+            if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_dhail1 > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_dhail1_max',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'max hourly hailcast hail diameter (embryo 1), output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Max hourly hailcast hail diameter (embryo 1)', 'mm', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hc1,kend_hc1, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_dhail2 > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_dhail2_max',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'max hourly hailcast hail diameter (embryo 2), output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Max hourly hailcast hail diameter (embryo 2)', 'mm', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hc2,kend_hc2, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_dhail3 > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_dhail3_max',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'max hourly hailcast hail diameter (embryo 3), output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Max hourly hailcast hail diameter (embryo 3)', 'mm', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hc3,kend_hc3, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_dhail4 > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_dhail4_max',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'max hourly hailcast hail diameter (embryo 4), output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Max hourly hailcast hail diameter (embryo 4)', 'mm', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hc4,kend_hc4, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_dhail5 > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_dhail5_max',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'max hourly hailcast hail diameter (embryo 5), output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Max hourly hailcast hail diameter (embryo 5)', 'mm', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hc5,kend_hc5, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_wdur > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_wdur',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'hailcast updraft duration, output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Hailcast updraft duration', 's', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hcd,kend_hcd, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+
+        if( id_hailcast_wup_mask > 0 ) then
+          call find_outputname(trim(file_name),'hailcast_wup_mask',output_name)
+          if(mpp_pe()==mpp_root_pe())print *,'hailcast updraft mask, output name=',trim(output_name)
+          call add_field_to_bundle(trim(output_name),'Hailcast updraft mask', '', "time: point",   &
+               axes(1:2), fcst_grid, kstt_hcm,kend_hcm, dyn_bundle, output_file, rcd=rc)
+          if(rc==0)  num_field_dyn=num_field_dyn+1
+        endif
+    endif
+
 !jwtest:
 !   call ESMF_FieldBundleGet(dyn_bundle, fieldCount=fieldCount, rc=rc)
 !   print *,'in dyn_bundle_setup, fieldCount=',fieldCount