[production/RRFS.v1] Merge the bugfix in the HAILCAST diagnostic code…

… (unit issue) (NOAA-GFDL#320) (#87)

* "include the bugfix in the HAILCAST diagnostic code (unit issue) (NOAA-GFDL#320)"

* "to output the average of smoke and dust"

driver/fvGFS/fv_nggps_diag.F90

@@ -117,11 +117,15 @@ module fv_nggps_diags_mod
  integer :: kstt_dbz, kend_dbz, kstt_omga, kend_omga
  integer :: kstt_windvect, kend_windvect
  integer :: kstt_o3_ave, kend_o3_ave, kstt_pm25_ave, kend_pm25_ave
+ integer :: kstt_smoke_ave, kend_smoke_ave
+ integer :: kstt_dust_ave, kend_dust_ave
+ integer :: kstt_coarsepm_ave, kend_coarsepm_ave
  integer :: kstt_no_ave, kend_no_ave, kstt_no2_ave, kend_no2_ave
  integer :: id_wmaxup,id_wmaxdn,kstt_wup, kend_wup,kstt_wdn,kend_wdn
  integer :: id_uhmax03,id_uhmin03,id_uhmax25,id_uhmin25,id_maxvort01
  integer :: id_maxvorthy1,kstt_maxvorthy1,kstt_maxvort01,id_ustm
  integer :: id_o3_ave,id_pm25_ave,id_no_ave,id_no2_ave
+ integer :: id_smoke_ave, id_dust_ave, id_coarsepm_ave
  integer :: kend_maxvorthy1,kend_maxvort01,id_vstm,id_srh01,id_srh03
  integer :: kstt_uhmax03,kstt_uhmin03,kend_uhmax03,kend_uhmin03
  integer :: kstt_uhmax25,kstt_uhmin25,kend_uhmax25,kend_uhmin25
@@ -161,6 +165,7 @@ module fv_nggps_diags_mod
  real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01
  real, dimension(:,:),allocatable :: maxvorthy1,maxvort02
  real, dimension(:,:,:),allocatable :: o3_ave,pm25_ave,no_ave,no2_ave
+ real, dimension(:,:),allocatable :: smoke_ave,dust_ave,coarsepm_ave
 
  public :: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 #ifdef use_WRTCOMP
@@ -466,6 +471,27 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
            kstt_pm25_ave = nlevs+1; kend_pm25_ave = nlevs+npzo
            nlevs = nlevs + npzo
        endif
+        id_smoke_ave = register_diag_field ( trim(file_name), 'smoke_ave',axes(1:2), Time,      &
+           'Hourly averaged surface smoke', 'ug/kg', missing_value=missing_value )
+       if( .not.Atm(n)%flagstruct%hydrostatic .and. id_smoke_ave > 0 ) then
+           allocate ( smoke_ave(isco:ieco,jsco:jeco) )
+           kstt_smoke_ave = nlevs+1; kend_smoke_ave = nlevs+1
+           nlevs = nlevs + npzo
+       endif
+        id_dust_ave = register_diag_field ( trim(file_name), 'dust_ave',axes(1:2), Time,      &
+           'Hourly averaged surface dust', 'ug/kg', missing_value=missing_value )
+       if( .not.Atm(n)%flagstruct%hydrostatic .and. id_dust_ave > 0 ) then
+           allocate ( dust_ave(isco:ieco,jsco:jeco) )
+           kstt_dust_ave = nlevs+1; kend_dust_ave = nlevs+1
+           nlevs = nlevs + npzo
+       endif
+        id_coarsepm_ave = register_diag_field ( trim(file_name), 'coarsepm_ave',axes(1:2), Time,      &
+           'Hourly averaged surface coarsepm', 'ug/kg', missing_value=missing_value )
+       if( .not.Atm(n)%flagstruct%hydrostatic .and. id_coarsepm_ave > 0 ) then
+           allocate ( coarsepm_ave(isco:ieco,jsco:jeco) )
+           kstt_coarsepm_ave = nlevs+1; kend_coarsepm_ave = nlevs+1
+           nlevs = nlevs + npzo
+       endif
 !
        nz = size(atm(1)%ak)
        allocate(ak(nz))
@@ -784,6 +810,18 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
     if ( id_pm25_ave > 0) then
       call store_data(id_pm25_ave, pm25_ave, Time, kstt_pm25_ave, kend_pm25_ave)
     endif
+    !--- hourly averaged surface smoke
+    if ( id_smoke_ave > 0) then
+      call store_data(id_smoke_ave, smoke_ave, Time, kstt_smoke_ave, kend_smoke_ave)
+    endif
+    !--- hourly averaged surface dust
+    if ( id_dust_ave > 0) then
+      call store_data(id_dust_ave, dust_ave, Time, kstt_dust_ave, kend_dust_ave)
+    endif
+    !--- hourly averaged surface coarsepm
+    if ( id_coarsepm_ave > 0) then
+      call store_data(id_coarsepm_ave, coarsepm_ave, Time, kstt_coarsepm_ave, kend_coarsepm_ave)
+    endif
     !--- max hourly 0-1km vert. vorticity
     if ( id_maxvort01 > 0) then
       call store_data(id_maxvort01, maxvort01, Time, kstt_maxvort01, kend_maxvort01)
@@ -876,6 +914,7 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
     real,                intent(in):: zvir
     integer :: i, j, k, n, ngc, nq, itrac
     integer :: o3_idx,no_idx,no2_idx,pm25_idx
+    integer :: smoke_idx,dust_idx,coarsepm_idx
     integer seconds, days, nsteps_per_reset
     logical, save :: first_call=.true.
     real, save :: first_time = 0.
@@ -1000,13 +1039,47 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
             Atm(n)%q,no2_idx,no2_ave,nsteps_per_reset,ucf1)
     call average_tracer_hy1(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
             Atm(n)%q,pm25_idx,pm25_ave,nsteps_per_reset,ucf2)
+
     kdtt1=kdtt1+1
    !else
    ! print *,'calculating hourly-averaegtd o3 or pm25'
    ! call mpp_error(FATAL, 'Missing hourly-averaged o3 or pm25 in diag_table')
    ! stop
    endif
+   if ( id_smoke_ave > 0 .and. id_dust_ave >0 &
+       .and. id_coarsepm_ave > 0 ) then
+       smoke_idx = get_tracer_index(MODEL_ATMOS, 'smoke')
+       dust_idx = get_tracer_index(MODEL_ATMOS, 'dust')
+       coarsepm_idx = get_tracer_index(MODEL_ATMOS, 'coarsepm')
+    if (first_call) then
+       call get_time(Time_step_atmos, seconds,  days)
+       first_time=seconds
+       first_call=.false.
+       kdtt1=0
+    endif
+       nsteps_per_reset = nint(avg_max_length/first_time)
+    if(mod(kdtt1,nsteps_per_reset)==0)then
+        do j=jsco,jeco
+         do i=isco,ieco
+            smoke_ave(i,j)= 0.
+            dust_ave(i,j)= 0.
+            coarsepm_ave(i,j)= 0.
+         enddo
+        enddo
+    endif
+    call average_tracer_hy2(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
+            Atm(n)%q,smoke_idx,smoke_ave,nsteps_per_reset,ucf2)
+    call average_tracer_hy2(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
+            Atm(n)%q,dust_idx,dust_ave,nsteps_per_reset,ucf2)
+    call average_tracer_hy2(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
+            Atm(n)%q,coarsepm_idx,coarsepm_ave,nsteps_per_reset,ucf2)
 
+    kdtt1=kdtt1+1
+   !else
+   ! print *,'calculating hourly-averaegtd o3 or pm25'
+   ! call mpp_error(FATAL, 'Missing hourly-averaged o3 or pm25 in diag_table')
+   ! stop
+   endif
    !allocate hailcast met field arrays
    if (do_hailcast) then
         call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, &
@@ -1017,23 +1090,42 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
  end subroutine fv_nggps_tavg
 !
  subroutine average_tracer_hy1(is,ie,js,je,isd,ied,jsd,jed, &
-                 ncns,npz,tracer,tr_idx,tracer_ave,nstp,unitcf) 
+                 ncns,npz,tracer,tr_idx,tracer_ave,nstp,unitcf)
    integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
    integer, intent(in):: ncns, npz, nstp, tr_idx
    real, intent(in) :: unitcf
-   real, intent(in), dimension(isd:ied,jsd:jed,npz,ncns):: tracer 
+   real, intent(in), dimension(isd:ied,jsd:jed,npz,ncns):: tracer
    real, intent(inout), dimension(is:ie,js:je,npz):: tracer_ave
    integer i, j, k
 !
    do k=1,npz
     do j=js,je
      do i=is,ie
         tracer_ave(i,j,k)=tracer_ave(i,j,k)+tracer(i,j,k,tr_idx)/nstp*unitcf
+     enddo
+    enddo
+   enddo
+
+ end subroutine average_tracer_hy1
+
+!
+ subroutine average_tracer_hy2(is,ie,js,je,isd,ied,jsd,jed, &
+                 ncns,npz,tracer,tr_idx,tracer_ave,nstp,unitcf)
+   integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
+   integer, intent(in):: ncns, npz, nstp, tr_idx
+   real, intent(in) :: unitcf
+   real, intent(in), dimension(isd:ied,jsd:jed,npz,ncns):: tracer
+   real, intent(inout), dimension(is:ie,js:je):: tracer_ave
+   integer i, j
+!
+
+    do j=js,je
+     do i=is,ie
+        tracer_ave(i,j)=tracer_ave(i,j)+tracer(i,j,npz,tr_idx)/nstp*unitcf
      enddo 
     enddo 
-   enddo  
 
- end subroutine average_tracer_hy1
+ end subroutine average_tracer_hy2
 
 !
  subroutine store_data(id, work, Time, nstt, nend)
@@ -1455,6 +1547,27 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
           axes(1:3), fcst_grid, kstt_pm25_ave,kend_pm25_ave, dyn_bundle, output_file, rcd=rc)
      if(rc==0)  num_field_dyn=num_field_dyn+1
    endif
+!
+   if ( id_smoke_ave > 0 ) then
+     call find_outputname(trim(file_name),'smoke_ave',output_name)
+     call add_field_to_bundle(trim(output_name),'hourly averaged surface smoke', 'ug/kg', "time: point",   &
+          axes(1:2), fcst_grid, kstt_smoke_ave,kend_smoke_ave, dyn_bundle, output_file, rcd=rc)
+     if(rc==0)  num_field_dyn=num_field_dyn+1
+   endif
+!
+   if ( id_dust_ave > 0 ) then
+     call find_outputname(trim(file_name),'dust_ave',output_name)
+     call add_field_to_bundle(trim(output_name),'hourly surface dust', 'ug/kg', "time: point",   &
+          axes(1:2), fcst_grid, kstt_dust_ave,kend_dust_ave, dyn_bundle, output_file, rcd=rc)
+     if(rc==0)  num_field_dyn=num_field_dyn+1
+   endif
+!
+   if ( id_coarsepm_ave > 0 ) then
+     call find_outputname(trim(file_name),'coarsepm_ave',output_name)
+     call add_field_to_bundle(trim(output_name),'hourly averaged surface coarsepm', 'ug/kg', "time: point",   &
+          axes(1:2), fcst_grid, kstt_coarsepm_ave,kend_coarsepm_ave, dyn_bundle, output_file, rcd=rc)
+     if(rc==0)  num_field_dyn=num_field_dyn+1
+   endif
 !
    if( id_ps > 0) then
      call find_outputname(trim(file_name),'ps',output_name)

tools/module_diag_hailcast.F90

@@ -1092,8 +1092,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,&
       IF (D.GT.0.254) D = 0.  !just consider missing for now if > 10 in
 
       !assign hail size in mm for output
-      !dhails(i) = D * 1000
-      dhails(i) = D
+      dhails(i) = D * 1000
 
     ENDDO  !end embryo size loop