Skip to content

Added code to perform hail size forecasting diagnostic #164

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
laurenchilutti merged 8 commits into from
Apr 29, 2022
Merged

Added code to perform hail size forecasting diagnostic #164

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
9e2653e
Added code to perform hail size forecasting diagnostic
ywangwof Nov 13, 2021
e1bdb6b
Modularized the hailcast module
ywangwof Mar 10, 2022
868a77e
Merge branch 'dev/emc' of https://github.com/NOAA-GFDL/GFDL_atmos_cub…
ywangwof Mar 11, 2022
da64dfd
Tested with Thompson and NSSL 2mom MP schemes
ywangwof Mar 12, 2022
7506f33
Restore to 43f7ed3
ywangwof Mar 18, 2022
e67bb42
Removed hailcast output dependency on other max/min hourly fields in …
ywangwof Mar 21, 2022
018b480
Improved code based PR comments
ywangwof Mar 22, 2022
584aa28
Merge branch 'dev/emc' into feature/HAILCAST
ywangwof Apr 28, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,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)
Expand Down
206 changes: 181 additions & 25 deletions driver/fvGFS/fv_nggps_diag.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
!
!------------------------------------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Comment thread
ywangwof marked this conversation as resolved.
!--- 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)
Expand Down Expand Up @@ -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
Copy link
Copy Markdown
Collaborator

@junwang-noaa junwang-noaa Mar 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove extra spaces?

Copy link
Copy Markdown
Author

@ywangwof ywangwof Mar 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is true. I would suggest that this practice should be established as a code convention. Anyway, almost all editors can be set to prune extra spaces at the end of lines automatically before saving.

Furthermore, we swap the orders of the IF statement and the DO loop statements for code efficiency. Thanks.

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,&
Expand All @@ -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, &
Expand All @@ -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:'
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
Loading