From 9e2653eaafd09547d62490da495f645731e5971a Mon Sep 17 00:00:00 2001 From: Yunheng Wang Date: Sat, 13 Nov 2021 00:29:18 +0000 Subject: [PATCH 1/6] Added code to perform hail size forecasting diagnostic --- CMakeLists.txt | 1 + driver/fvGFS/fv_nggps_diag.F90 | 351 ++++++- model/fv_arrays.F90 | 4 + tools/module_diag_hailcast.F90 | 1637 ++++++++++++++++++++++++++++++++ 4 files changed, 1986 insertions(+), 7 deletions(-) create mode 100644 tools/module_diag_hailcast.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 888a622d9..9b25ae4af 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index f93df68d2..d16a37e86 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -62,8 +62,9 @@ module fv_nggps_diags_mod ! ! - use mpp_mod, only: mpp_pe, mpp_root_pe,FATAL,mpp_error + use mpp_mod, only: mpp_pe, mpp_root_pe,FATAL,mpp_error,NOTE, stdlog, input_nml_file use constants_mod, only: grav, rdgas + use fms_mod, only: check_nml_error use time_manager_mod, only: time_type, get_time use diag_manager_mod, only: register_diag_field, send_data use diag_axis_mod, only: get_axis_global_length, get_diag_axis, get_diag_axis_name @@ -74,6 +75,7 @@ module fv_nggps_diags_mod use fv_diagnostics_mod, only: range_check, dbzcalc,max_vv,get_vorticity, & max_uh,max_vorticity,bunkers_vector, & helicity_relative_CAPS,max_vorticity_hy1 + use module_diag_hailcast, only: hailcast_diagnostic_driver use fv_arrays_mod, only: fv_atmos_type use mpp_domains_mod, only: domain1d, domainUG #ifdef MULTI_GASES @@ -99,11 +101,18 @@ module fv_nggps_diags_mod 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_hailcast_dhail + integer :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5, id_hailcast_wdur, id_hailcast_wup_mask + integer :: id_hailcast_diam_mean, id_hailcast_diam_std 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 integer :: kstt_ustm,kstt_vstm,kend_ustm,kend_vstm,kstt_srh01 integer :: kstt_srh03,kend_srh01,kend_srh03 + integer :: kstt_hc, kend_hc + integer :: kstt_hc1,kstt_hc2,kstt_hc3,kstt_hc4,kstt_hc5 + integer :: kend_hc1,kend_hc2,kend_hc3,kend_hc4,kend_hc5 + integer :: kstt_hcd, kend_hcd, kstt_hcm, kend_hcm integer :: id_maxvort02,kstt_maxvort02,kend_maxvort02 integer :: isco, ieco, jsco, jeco, npzo, ncnsto integer :: isdo, iedo, jsdo, jedo @@ -136,6 +145,12 @@ module fv_nggps_diags_mod real, dimension(:,:),allocatable :: up2,dn2,uhmax03,uhmin03 real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01 real, dimension(:,:),allocatable :: maxvorthy1,maxvort02 + + real, allocatable :: hailcast_dhail1(:,:), hailcast_dhail2(:,:), hailcast_dhail3(:,:), hailcast_dhail4(:,:), hailcast_dhail5(:,:) !hailstone diameters (mm) + real, allocatable :: hailcast_dhail1_max(:,:), hailcast_dhail2_max(:,:), hailcast_dhail3_max(:,:), hailcast_dhail4_max(:,:), hailcast_dhail5_max(:,:) !hailstone diameters (mm) + real, allocatable :: hailcast_diam_mean(:,:), hailcast_diam_std(:,:) !mean and standard deviation of five hailstones (mm) + real, allocatable :: hailcast_wdur(:,:), hailcast_wup_mask(:,:) !persistent arrays for updraft duration (s) and mask + public :: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg #ifdef use_WRTCOMP public :: fv_dyn_bundle_setup @@ -148,6 +163,32 @@ subroutine fv_nggps_diag_init(Atm, axes, Time) integer, intent(in) :: axes(4) type(time_type), intent(in) :: Time integer :: n, i, j, nz + logical :: do_hailcast = .true. !This controls whether hailcast is used + namelist /fv_diagnostics_nml/ do_hailcast + integer :: ios, ierr + integer :: unit + + !namelist file for hailcast +#ifdef INTERNAL_FILE_NML + ! Read Main namelist + read (input_nml_file,fv_diagnostics_nml,iostat=ios) + ierr = check_nml_error(ios,'fv_diagnostics_nml') +#else + f_unit=open_namelist_file() + rewind (f_unit) + ! Read Main namelist + read (f_unit,fv_diagnostics_nml,iostat=ios) + ierr = check_nml_error(ios,'fv_diagnostics_nml') + call close_file(f_unit) +#endif + + unit = stdlog() + write(unit, nml=fv_diagnostics_nml) + !end hailcast nml + + if (mpp_pe() == mpp_root_pe()) then + print*, 'do_hailcast = ', do_hailcast + end if n = 1 ncnsto = Atm(1)%ncnst @@ -427,7 +468,97 @@ 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 + + !register hailcast arrays + if (do_hailcast) then + id_hailcast_dhail = register_diag_field ( trim(file_name), 'hailcast_dhail_max', axes(1:2), Time, & + 'hourly max hail diameter', 'mm', missing_value=missing_value ) + id_hailcast_dhail1 = register_diag_field ( trim(file_name), 'hailcast_dhail1_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 1)', 'mm', missing_value=missing_value ) + id_hailcast_dhail2 = register_diag_field ( trim(file_name), 'hailcast_dhail2_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 2)', 'mm', missing_value=missing_value ) + id_hailcast_dhail3 = register_diag_field ( trim(file_name), 'hailcast_dhail3_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 3)', 'mm', missing_value=missing_value ) + id_hailcast_dhail4 = register_diag_field ( trim(file_name), 'hailcast_dhail4_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 4)', 'mm', missing_value=missing_value ) + id_hailcast_dhail5 = register_diag_field ( trim(file_name), 'hailcast_dhail5_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 5)', 'mm', missing_value=missing_value ) + id_hailcast_diam_mean = register_diag_field ( trim(file_name), 'hailcast_diam_mean', axes(1:2), Time, & + 'mean hail diameter', 'mm', missing_value=missing_value ) + id_hailcast_diam_std = register_diag_field ( trim(file_name), 'hailcast_diam_std', axes(1:2), Time, & + 'standard deviation of hail diameter', 'mm', missing_value=missing_value ) + id_hailcast_wdur = register_diag_field ( trim(file_name), 'hailcast_wdur', axes(1:2), Time, & + 'updraft duration', 's', missing_value=missing_value ) + id_hailcast_wup_mask = register_diag_field ( trim(file_name), 'hailcast_wup_mask', axes(1:2), Time, & + 'updraft mask', '', missing_value=missing_value ) + + if (id_hailcast_dhail > 0) then + kstt_hc = nlevs+1; kend_hc = nlevs+1 + nlevs = nlevs + 1 + endif + + if (id_hailcast_dhail1 > 0) then + kstt_hc1 = nlevs+1; kend_hc1 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail2 > 0) then + kstt_hc2 = nlevs+1; kend_hc2 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail3 > 0) then + kstt_hc3 = nlevs+1; kend_hc3 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail4 > 0) then + kstt_hc4 = nlevs+1; kend_hc4 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail5 > 0) then + kstt_hc5 = nlevs+1; kend_hc5 = nlevs+1 + nlevs = nlevs + 1 endif + if (id_hailcast_wdur > 0) then + kstt_hcd = nlevs+1; kend_hcd = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_wup_mask > 0) then + kstt_hcm = nlevs+1; kend_hcm = nlevs+1 + nlevs = nlevs + 1 + endif + + if (.not.allocated(hailcast_dhail1)) then + allocate(hailcast_dhail1(isco:ieco,jsco:jeco), & + hailcast_dhail2(isco:ieco,jsco:jeco), & + hailcast_dhail3(isco:ieco,jsco:jeco), & + hailcast_dhail4(isco:ieco,jsco:jeco), & + hailcast_dhail5(isco:ieco,jsco:jeco), & + hailcast_diam_mean(isco:ieco,jsco:jeco), & + hailcast_diam_std(isco:ieco,jsco:jeco)) + allocate ( hailcast_dhail1_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail2_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail3_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail4_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail5_max(isco:ieco,jsco:jeco) ) + + do i=isco,ieco + do j=jsco,jeco + hailcast_dhail1(i,j)=0; hailcast_dhail2(i,j)=0; hailcast_dhail3(i,j)=0; hailcast_dhail4(i,j)=0; hailcast_dhail5(i,j)=0 + hailcast_dhail1_max(i,j)=0; hailcast_dhail2_max(i,j)=0; hailcast_dhail3_max(i,j)=0; hailcast_dhail4_max(i,j)=0; hailcast_dhail5_max(i,j)=0 + enddo + enddo + endif + if (.not.allocated(hailcast_wdur)) then + allocate(hailcast_wdur(isdo:iedo,jsdo:jedo),hailcast_wup_mask(isdo:iedo,jsdo:jedo)) + do i=isdo,iedo + do j=jsdo,jedo + hailcast_wdur(i,j)=0;hailcast_wup_mask(i,j)=0 + enddo + enddo + endif + endif + ! !------------------------------------ ! use wrte grid component for output @@ -447,7 +578,7 @@ subroutine fv_nggps_diag(Atm, zvir, Time) logical :: bad_range real :: ptop, allmax real, allocatable :: wk(:,:,:), wk2(:,:,:) - real, dimension(:,:),allocatable :: ustm,vstm,srh01,srh03 + real, dimension(:,:),allocatable :: ustm,vstm,srh01,srh03,hailcast_dhail_max n = 1 ngc = Atm(n)%ng @@ -734,6 +865,62 @@ subroutine fv_nggps_diag(Atm, zvir, Time) call store_data(id_uhmin25, uhmin25, Time, kstt_uhmin25, kend_uhmin25) endif + !--- max hourly hailcast hail diameter + if ( .not.Atm(n)%flagstruct%hydrostatic .and. id_hailcast_dhail > 0) then + allocate(hailcast_dhail_max(isco:ieco,jsco:jeco)) + + do j=jsco,jeco + do i=isco,ieco + hailcast_dhail_max(i,j) = 0. + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail1_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail2_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail3_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail4_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail5_max(i,j)) + end do + end do + + call store_data(id_hailcast_dhail, hailcast_dhail_max, Time, kstt_hc, kend_hc) + deallocate(hailcast_dhail_max) + endif + + !--- max hourly hailcast hail diameter (embryo 1) + if ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. id_hailcast_dhail5 > 0) then + call store_data(id_hailcast_dhail5, hailcast_dhail5_max, Time, kstt_hc5, kend_hc5) + endif + + !--- hailcast updraft duration + if ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 + + + !call nullify_domain() + end subroutine fv_nggps_diag subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) @@ -747,6 +934,7 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) integer, save :: kdtt = 0 real :: avg_max_length real,dimension(:,:,:),allocatable :: vort + real, allocatable :: rhoair_layer(:,:,:), z(:), z_layer(:,:,:), p_layer(:,:,:),zsfc(:,:) n = 1 ngc = Atm(n)%ng nq = size (Atm(n)%q,4) @@ -770,17 +958,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) + if(mod(kdtt,nsteps_per_reset)==0)then 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 + 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 +986,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) + if(mod(kdtt,nsteps_per_reset)==0)then 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 + endif call max_uh(isco,ieco,jsco,jeco,ngc,npzo,zvir, & sphum,uhmax03,uhmin03,Atm(n)%w,vort,Atm(n)%delz, & @@ -820,7 +1008,6 @@ 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) else print *,'Missing max/min hourly field in diag_table' @@ -830,6 +1017,92 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) stop endif endif + + !allocate hailcast met field arrays + if (id_hailcast_dhail>0 .or. id_hailcast_dhail1>0 .or. id_hailcast_dhail2>0 .or. id_hailcast_dhail3>0 & + .or. id_hailcast_dhail4>0 .or. id_hailcast_dhail5>0 .or. id_hailcast_diam_mean>0 & + .or. id_hailcast_diam_std>0) then + + do j=jsco,jeco + do i=isco,ieco + if(mod(kdtt,nsteps_per_reset)==0)then + hailcast_dhail1_max(i,j) = 0. + hailcast_dhail2_max(i,j) = 0. + hailcast_dhail3_max(i,j) = 0. + hailcast_dhail4_max(i,j) = 0. + hailcast_dhail5_max(i,j) = 0. + endif + + hailcast_dhail1(i,j) = 0 + hailcast_dhail2(i,j) = 0 + hailcast_dhail3(i,j) = 0 + hailcast_dhail4(i,j) = 0 + hailcast_dhail5(i,j) = 0 + enddo + enddo + + + if (.not. allocated(rhoair_layer)) allocate(rhoair_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D air density + if (.not. allocated(z)) allocate(z(1:npzo+1)) !interace levels + if (.not. allocated(z_layer)) allocate(z_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D height ASL + if (.not. allocated(p_layer)) allocate(p_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D pressure + if (.not. allocated(zsfc)) allocate(zsfc(isco:ieco,jsco:jeco)) !terrain height + + do k=npzo,1,-1 + do j=jsco, jeco + if (Atm(n)%flagstruct%hydrostatic) then + call mpp_error(FATAL, 'HAILCAST can only be run with non-hydrostatic FV3') + !do i=is, ie + !rhoair(i,j,k) = Atm(n)%delp(i,j,k)/( ( Atm(n)%peln(i,k+1,j)- Atm(n)%peln(i,k,j)) & + ! * rdgas * Atm(n)%pt(i,j,k) * ( 1. + zvir*q(i,j,k,sphum))) + !enddo + else !non-hydrostatic + do i=isco, ieco + if (k==npzo) then + zsfc(i,j)=Atm(n)%phis(i,j) / grav + z(npzo+1)=zsfc(i,j) + endif + rhoair_layer(i,j,k) = -Atm(n)%delp(i,j,k)/(grav*Atm(n)%delz(i,j,k)) + !height of interfaces and layer height + z(k)=z(k+1)-Atm(n)%delz(i,j,k) + z_layer(i,j,k)=(z(k+1)+z(k))/2 + p_layer(i,j,k)=Atm(n)%delp(i,j,k)*(1.-sum(Atm(n)%q(i,j,k,2:Atm(n)%flagstruct%nwat)))/& + (-Atm(n)%delz(i,j,k)*grav)*rdgas*Atm(n)%pt(i,j,k)*(1.+zvir*Atm(n)%q(i,j,k,sphum)) + enddo + endif + enddo !j loop + enddo !k loop + + !call hailcast diagnostic driver once per time step and provide three-dimensional met fields + call hailcast_diagnostic_driver(Atm(n)%pt(isco:ieco,jsco:jeco,1:npzo), Atm(n)%w(isco:ieco,jsco:jeco,1:npzo), p_layer, z_layer, rhoair_layer, Atm(n)%q(isco:ieco,jsco:jeco,1:npzo,:), & !3D fields + Atm(n)%flagstruct%nwat, sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !number of tracer indices, indices + isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo, & !grid dimensions data array (with halo) and physical grid dimensions + 1,npzo, zsfc, & !vertical dimensions, terrain height + first_time, Atm(n)%domain, & !call hailcast every model step and info for updating haloes + hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5, & !hailcast embryos sizes + hailcast_diam_mean, hailcast_diam_std, & !hailcast embryo mean/std + hailcast_wdur, hailcast_wup_mask) !persistent updraft duration and mask + + do j=jsco,jeco + do i=isco,ieco + hailcast_dhail1_max(i,j) = max(hailcast_dhail1_max(i,j), hailcast_dhail1(i,j)) + hailcast_dhail2_max(i,j) = max(hailcast_dhail2_max(i,j), hailcast_dhail2(i,j)) + hailcast_dhail3_max(i,j) = max(hailcast_dhail3_max(i,j), hailcast_dhail3(i,j)) + hailcast_dhail4_max(i,j) = max(hailcast_dhail4_max(i,j), hailcast_dhail4(i,j)) + hailcast_dhail5_max(i,j) = max(hailcast_dhail5_max(i,j), hailcast_dhail5(i,j)) + enddo + enddo + + + !deallocate hailcast met variables + deallocate(p_layer) + deallocate(z_layer) + deallocate(rhoair_layer) + deallocate(z) + deallocate(zsfc) + endif + + kdtt=kdtt+1 end subroutine fv_nggps_tavg ! subroutine store_data(id, work, Time, nstt, nend) @@ -1335,6 +1608,70 @@ 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. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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 + !jwtest: ! call ESMF_FieldBundleGet(dyn_bundle, fieldCount=fieldCount, rc=rc) ! print *,'in dyn_bundle_setup, fieldCount=',fieldCount diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index a353e0a65..34e5a2fca 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -61,6 +61,10 @@ module fv_arrays_mod real :: efx(max_step), efx_sum, efx_nest(max_step), efx_sum_nest, mtq(max_step), mtq_sum integer :: steps + !for hailcast + integer :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5 + integer :: id_hailcast_diam_mean, id_hailcast_diam_std + end type fv_diag_type diff --git a/tools/module_diag_hailcast.F90 b/tools/module_diag_hailcast.F90 new file mode 100644 index 000000000..8aae38314 --- /dev/null +++ b/tools/module_diag_hailcast.F90 @@ -0,0 +1,1637 @@ +!Hailcase diagnostic driver +!John Henderson AER 20190425 +!time management handling from WRF is preserved but commented out + +!MODULE module_diag_hailcast +!CONTAINS +! SUBROUTINE diag_hailcast_stub +! END SUBROUTINE diag_hailcast_stub +!END MODULE module_diag_hailcast + +MODULE module_diag_hailcast + +CONTAINS + +! SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags & +! , moist & +! , rho & +! , ids, ide, jds, jde, kds, kde & +! , ims, ime, jms, jme, kms, kme & +! , ips, ipe, jps, jpe, kps, kpe & +! , its, ite, jts, jte & +! , k_start, k_end & +! , dt, itimestep & +! , haildt,curr_secs & +! , haildtacttime ) + SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3D temperature, updraft, pressure, height, density and moisture tracers + moist_count,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & ! tracer indices + is,ie,js,je,ishalo,iehalo,jshalo,jehalo, & ! i,j halo data array dimensions and physical grid dimensions (2 smaller on each end) + k_start, k_end, terr_z, & ! bottom and top z indices and terrain height + dt, domainid, & ! call frequency of hailcast = model time step = physics + hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5, hailcast_diam_mean, hailcast_diam_std, hailcast_wdur, hailcast_wup_mask) !hailcast persistent variables + +! USE module_domain, ONLY : domain , domain_clock_get +! USE module_configure, ONLY : grid_config_rec_type, model_config_rec +! USE module_state_description +! USE module_model_constants +! USE module_utility +! USE module_streams, ONLY: history_alarm, auxhist2_alarm +!#ifdef DM_PARALLEL +! USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval +! #endif + + use constants_mod, only: grav, rdgas, rvgas, pi=>pi_8, radius, kappa, WTMAIR, WTMCO2, & + omega, hlv, cp_air, cp_vapor + !use fms_mod, only: write_version_number + use fms_io_mod, only: set_domain, nullify_domain + use time_manager_mod, only: time_type, get_date, get_time + use mpp_domains_mod, only: domain2d, mpp_update_domains, DGRID_NE + use diag_manager_mod, only: diag_axis_init, register_diag_field, & + register_static_field, send_data, diag_grid_init + use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_diag_type, fv_grid_bounds_type, & + R_GRID + use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index + use field_manager_mod, only: MODEL_ATMOS + use fv_mp_mod, only: mp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master + use mpp_mod, only: mpp_error, FATAL, stdlog, mpp_pe, mpp_root_pe, mpp_sum, mpp_max, NOTE + use sat_vapor_pres_mod, only: compute_qs, lookup_es + + use fv_arrays_mod, only: max_step + + + IMPLICIT NONE + +! TYPE ( domain ), INTENT(INOUT) :: grid +! TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags +! type(fv_atmos_type), intent(inout) :: grid +!!!! type(time_type), intent(in) :: Time !will not compile + INTEGER:: ishalo, iehalo, jshalo, jehalo, k_start, k_end !dimensions of halo data array + INTEGER:: is, ie, js, je + +! INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, & +! ims, ime, jms, jme, kms, kme, & +! ips, ipe, jps, jpe, kps, kpe + + logical :: master + INTEGER :: moist_count,sphum,rainwat,snowwat,graupel,liq_wat,ice_wat + + REAL, DIMENSION( is:ie, js:je, k_start: k_end, moist_count), & + INTENT(IN ) :: q !all 3D moisture fields + + REAL, DIMENSION( is:ie, js:je, k_start: k_end), & + INTENT(IN ) :: t, w, p, z, rho !3D non-moisture fields + + REAL, DIMENSION( is:ie, js:je), & + INTENT(IN ) :: terr_z + + REAL, DIMENSION( is:ie, js:je), & + INTENT(INOUT) :: hailcast_dhail1, hailcast_dhail2, & + hailcast_dhail3, hailcast_dhail4, & + hailcast_dhail5, & + hailcast_diam_mean, hailcast_diam_std + + REAL, DIMENSION( ishalo:iehalo, jshalo:jehalo), & + INTENT(INOUT) :: hailcast_wdur, hailcast_wup_mask + +! REAL, INTENT(IN ),OPTIONAL :: haildt +! REAL, INTENT(IN ),OPTIONAL :: curr_secs +! REAL, INTENT(INOUT),OPTIONAL :: haildtacttime +! INTEGER, INTENT(IN ) :: itimestep + REAL, INTENT(IN ) :: dt + TYPE(domain2d), INTENT(INOUT ) :: domainid + + + ! Local + ! ----- + CHARACTER*512 :: message + CHARACTER*256 :: timestr + INTEGER :: i,j,k, numlevs + REAL, DIMENSION( is:ie, js:je, k_start:k_end ) :: qr & + , qs & + , qg & + , qv & + , qc & + , qi & + , ptot + + REAL, DIMENSION( ishalo:iehalo, jshalo:jehalo) :: wdur_prev, wup_mask_prev + + REAL :: dhail1,dhail2,dhail3,dhail4,dhail5 + +! ! Timing +! ! ------ +! TYPE(WRFU_Time) :: hist_time, aux2_time, CurrTime, StartTime +! TYPE(WRFU_TimeInterval) :: dtint, histint, aux2int +! LOGICAL :: is_after_history_dump, is_output_timestep, is_first_timestep +! LOGICAL :: doing_adapt_dt, run_param, decided +! INTEGER :: stephail, itimestep_basezero + + +! ! Chirp the routine name for debugging purposes +! ! --------------------------------------------- +! write ( message, * ) 'inside hailcast_diagnostics_driver' +! CALL wrf_debug( 100 , message ) + +! ! Get timing info +! ! Want to know if when the last history output was +! ! Check history and auxhist2 alarms to check last ring time and how often +! ! they are set to ring +! ! ----------------------------------------------------------------------- +! CALL WRFU_ALARMGET( grid%alarms( HISTORY_ALARM ), prevringtime=hist_time, & +! ringinterval=histint) +! CALL WRFU_ALARMGET( grid%alarms( AUXHIST2_ALARM ), prevringtime=aux2_time, & +! ringinterval=aux2int) + +! ! Get domain clock +! ! ---------------- +! CALL domain_clock_get ( grid, current_time=CurrTime, & +! simulationStartTime=StartTime, & +! current_timestr=timestr, time_step=dtint ) + +! ! Set some booleans for use later +! ! Following uses an overloaded .lt. +! ! --------------------------------- +! is_after_history_dump = ( Currtime .lt. hist_time + dtint ) + +! ! Following uses an overloaded .ge. +! ! --------------------------------- +! is_output_timestep = (Currtime .ge. hist_time + histint - dtint .or. & +! Currtime .ge. aux2_time + aux2int - dtint ) +! write ( message, * ) 'is output timestep? ', is_output_timestep +! CALL wrf_debug( 100 , message ) + +! ! Following uses an overloaded .eq. +! ! --------------------------------- +! is_first_timestep = ( Currtime .eq. StartTime + dtint ) + + master = (mpp_pe()==mpp_root_pe()) .or. is_master() + + !FV3: No need to reset arrays to zero at dump time, thus no need to keep track of dump times + ! + ! After each history dump, reset max/min value arrays + ! ---------------------------------------------------------------------- +! IF ( is_after_history_dump ) THEN + !! DO j = jsc, jec !jts, jte + !! DO i = isc, iec !its, ite + !! hailcast_dhail1(i,j) = 0. + !! hailcast_dhail2(i,j) = 0. + !! hailcast_dhail3(i,j) = 0. + !!! hailcast_dhail4(i,j) = 0. + !! hailcast_dhail5(i,j) = 0. + !! ENDDO + !! ENDDO + !ENDIF ! is_after_history_dump + + ! We need to do some neighboring gridpoint comparisons for the updraft + ! duration calculations; set i,j start and end values so we don't go off + ! the edges of the domain. Updraft duration on domain edges will always be 0. + ! ---------------------------------------------------------------------- + +! !WRF config_flags contains: open_xs, specified, nested,open_xe,open_yx,open_ye,periodic_x +! IF ( config_flags%open_xs .OR. config_flags%specified .OR. & +! config_flags%nested) i_start = MAX( ids+1, its ) +! IF ( config_flags%open_xe .OR. config_flags%specified .OR. & +! config_flags%nested) i_end = MIN( ide-2, ite ) !-1 +! IF ( config_flags%open_ys .OR. config_flags%specified .OR. & +! config_flags%nested) j_start = MAX( jds+1, jts ) +! IF ( config_flags%open_ye .OR. config_flags%specified .OR. & +! config_flags%nested) j_end = MIN( jde-2, jte ) !-1 +! IF ( config_flags%periodic_x ) i_start = its +! IF ( config_flags%periodic_x ) i_end = ite + + + ! Make a copy of the updraft duration, mask variables + ! --------------------------------------------------- + !from WRF: + !!wdur_prev(:,:) = grid%hailcast_wdur(:,:) + !!wup_mask_prev(:,:) = grid%hailcast_wup_mask(:,:) + !DO j = MAX( jts-1 , jds ), MIN( jte+1 , jde-1 ) + ! DO i = MAX( its-1 , ids ), MIN( ite+1 , ide-1) + DO j = jshalo, jehalo + DO i = ishalo, iehalo + wdur_prev(i,j) = hailcast_wdur(i,j) + wup_mask_prev(i,j) = hailcast_wup_mask(i,j) + ENDDO + ENDDO + !update halo: + call mpp_update_domains(wup_mask_prev,domainid) + call mpp_update_domains(wdur_prev,domainid) + + ! Determine updraft mask (where updraft greater than some threshold) + ! --------------------------------------------------- + DO j = js, je + DO i = is, ie + hailcast_wup_mask(i,j) = 0 + hailcast_wdur(i,j) = 0 + + DO k = k_start, k_end + IF ( w(i,j,k) .ge. 10. ) THEN + hailcast_wup_mask(i,j) = 1 + ENDIF + ENDDO + ENDDO + ENDDO + + ! Determine updraft duration; make sure not to call point outside the domain + ! --------------------------------------------------- + DO j = js,je + DO i = is,ie + + ! Determine updraft duration using updraft masks + ! --------------------------------------------------- + IF ( (hailcast_wup_mask(i,j).eq.1) .OR. & + (MAXVAL(wup_mask_prev(i-1:i+1,j-1:j+1)).eq.1) ) THEN + hailcast_wdur(i,j) = MAXVAL(wdur_prev(i-1:i+1,j-1:j+1)) + dt + ENDIF + ENDDO + ENDDO + + !from WRF dealing with calling frequency and adaptive time step; call hailcast every timestep for FV3 + ! + ! Only run the scheme every "haildt" seconds as set in the namelist. + ! Code borrowed from module_pbl_driver, although here haildt is + ! in seconds, not minutes (bldt is in minutes). + ! --------------------------------------------------- + +! ! Are we using adaptive timestepping? +! doing_adapt_dt = .FALSE. +! IF ( (config_flags%use_adaptive_time_step) .and. & +! ( (.not. grid%nested) .or. & +! ( (grid%nested) .and. (abs(grid%dtbc) < 0.0001) ) ) )THEN +! doing_adapt_dt = .TRUE. +! IF ( haildtacttime .eq. 0. ) THEN +! haildtacttime = CURR_SECS + haildt +! END IF +! END IF + +! ! Calculate STEPHAIL +! stephail = NINT(haildt / dt) +! stephail = MAX(stephail,1) + +! ! itimestep starts at 1 - we want it to start at 0 so correctly +! ! divisibile by stephail. +! itimestep_basezero = itimestep -1 + + +! ! Do we run through this scheme or not? +! ! Test 1: If this is the initial model time, then yes. +! ! ITIMESTEP=1 +! ! Test 2: If the user asked for hailcast to be run every time step, then yes. +! ! HAILDT=0 or STEPHAIL=1 +! ! Test 3: If not adaptive dt, and this is on the requested hail frequency, +! ! then yes. +! ! MOD(ITIMESTEP,STEPHAIL)=0 +! ! Test 4: If using adaptive dt and the current time is past the last +! ! requested activate hailcast time, then yes. +! ! CURR_SECS >= HAILDTACTTIME + +! ! If we do run through the scheme, we set the flag run_param to TRUE and we set +! ! the decided flag to TRUE. The decided flag says that one of these tests was +! ! able to say "yes", run the scheme. We only proceed to other tests if the +! ! previous tests all have left decided as FALSE. If we set run_param to TRUE +! ! and this is adaptive time stepping, we set the time to the next hailcast run. + +! run_param = .FALSE. +! decided = .FALSE. +! IF ( ( .NOT. decided ) .AND. & +! ( itimestep_basezero .EQ. 0 ) ) THEN +! run_param = .TRUE. +! decided = .TRUE. +! END IF + +! IF ( PRESENT(haildt) )THEN +! IF ( ( .NOT. decided ) .AND. & +! ( ( haildt .EQ. 0. ) .OR. ( stephail .EQ. 1 ) ) ) THEN +! run_param = .TRUE. +! decided = .TRUE. +! END IF +! ELSE +! IF ( ( .NOT. decided ) .AND. & +! ( stephail .EQ. 1 ) ) THEN +! run_param = .TRUE. +! decided = .TRUE. +! END IF +! END IF + +! IF ( ( .NOT. decided ) .AND. & +! ( .NOT. doing_adapt_dt ) .AND. & +! ( MOD(itimestep_basezero,stephail) .EQ. 0 ) ) THEN +! run_param = .TRUE. +! decided = .TRUE. +! END IF + +! IF ( ( .NOT. decided ) .AND. & +! ( doing_adapt_dt ) .AND. & +! ( curr_secs .GE. haildtacttime ) ) THEN +! run_param = .TRUE. +! decided = .TRUE. +! haildtacttime = curr_secs + haildt +! END IF + +! write ( message, * ) 'timestep to run HAILCAST? ', run_param +! CALL wrf_debug( 100 , message ) + +! IF (run_param) THEN + + ! 3-D arrays for moisture variables + ! --------------------------------- + numlevs=k_end-k_start+1 + DO i=is, ie + DO k=1, numlevs + DO j=js, je + qv(i,j,k)=0;qr(i,j,k)=0;qs(i,j,k)=0;qg(i,j,k)=0;qc(i,j,k)=0;qi(i,j,k)=0 + qv(i,j,k) = q(i,j,k,sphum) + qr(i,j,k) = q(i,j,k,rainwat) + qs(i,j,k) = q(i,j,k,snowwat) + qg(i,j,k) = q(i,j,k,graupel) + qc(i,j,k) = q(i,j,k,liq_wat) + qi(i,j,k) = q(i,j,k,ice_wat) + ENDDO + ENDDO + ENDDO + + ! Hail diameter in millimeters (HAILCAST) + ! --------------------------------------------------- + DO j = js, je + DO i = is, ie + ! Only call hailstone driver if updraft has been + ! around longer than 15 min + ! ---------------------------------------------- + IF (hailcast_wdur(i,j) .gt. 900) THEN + dhail1=0;dhail2=0;dhail3=0;dhail4=0;dhail5=0 + CALL hailstone_driver (i,j, t(i,j,k_end:k_start:-1), & + z(i,j,k_end:k_start:-1), & + terr_z(i,j), & + p(i,j,k_end:k_start:-1), & + rho(i,j,k_end:k_start:-1), & + qv(i,j,k_end:k_start:-1), & + qi(i,j,k_end:k_start:-1), & + qc(i,j,k_end:k_start:-1), & + qr(i,j,k_end:k_start:-1), & + qs(i,j,k_end:k_start:-1), & + qg(i,j,k_end:k_start:-1), & + w(i,j,k_end:k_start:-1), & + hailcast_wdur(i,j), & + numlevs, & + dhail1, dhail2, & + dhail3, dhail4, & + dhail5 ) + + !diag_table should select largest value of these variables during requested period + !thus do not need to accumulate largest value between data dumps + !IF (dhail1 .gt. hailcast_dhail1(i,j)) THEN + hailcast_dhail1(i,j) = dhail1 + !ENDIF + !IF (dhail2 .gt. hailcast_dhail2(i,j)) THEN + hailcast_dhail2(i,j) = dhail2 + !ENDIF + !IF (dhail3 .gt. hailcast_dhail3(i,j)) THEN + hailcast_dhail3(i,j) = dhail3 + !ENDIF + !IF (dhail4 .gt. hailcast_dhail4(i,j)) THEN + hailcast_dhail4(i,j) = dhail4 + !ENDIF + !IF (dhail5 .gt. hailcast_dhail5(i,j)) THEN + hailcast_dhail5(i,j) = dhail5 + !ENDIF + ENDIF + ENDDO + ENDDO + + ! Calculate the mean and standard deviation of the hail diameter + ! distribution over different embryo sizes + ! ---------------------------------------- + DO j = js, je + DO i = is, ie + !mean + hailcast_diam_mean(i,j)=(hailcast_dhail1(i,j)+& + hailcast_dhail2(i,j) +hailcast_dhail3(i,j)+& + hailcast_dhail4(i,j) +hailcast_dhail5(i,j))/5. + !max + !hailcast_diam_max(i,j)=MAX(hailcast_dhail1(i,j),& + ! hailcast_dhail2(i,j), hailcast_dhail3(i,j),& + ! hailcast_dhail4(i,j), hailcast_dhail5(i,j)) + !sample standard deviation + hailcast_diam_std(i,j) = SQRT( ( & + (hailcast_dhail1(i,j)-hailcast_diam_mean(i,j))**2.+& + (hailcast_dhail2(i,j)-hailcast_diam_mean(i,j))**2.+& + (hailcast_dhail3(i,j)-hailcast_diam_mean(i,j))**2.+& + (hailcast_dhail4(i,j)-hailcast_diam_mean(i,j))**2.+& + (hailcast_dhail5(i,j)-hailcast_diam_mean(i,j))**2.)& + / 4.0) + if (hailcast_diam_mean(i,j).ne.0.0 .and. hailcast_diam_std(i,j).ne.0.0) print*,'jmh mean std=',i,j,hailcast_diam_mean(i,j),hailcast_diam_std(i,j) + ENDDO + ENDDO + !END IF + +END SUBROUTINE hailcast_diagnostic_driver + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!! +!!!! Hailstone driver, adapted from hailstone subroutine in HAILCAST +!!!! Inputs: +!!!! 1-d (nz) +!!!! TCA temperature (K) +!!!! h1d height above sea level (m) +!!!! PA total pressure (Pa) +!!!! rho1d density (kg/m3) +!!!! RA vapor mixing ratio (kg/kg) +!!!! qi1d cloud ice mixing ratio (kg/kg) +!!!! qc1d cloud water mixing ratio (kg/kg) +!!!! qr1d rain water mixing ratio (kg/kg) +!!!! qg1d graupel mixing ratio (kg/kg) +!!!! qs1d snow mixing ratio (kg/kg) +!!!! VUU updraft speed at each level (m/s) +!!!! Float +!!!! ht terrain height (m) +!!!! wdur duration of any updraft > 10 m/s within 1 surrounding +!!!! gridpoint +!!!! Integer +!!!! nz number of vertical levels +!!!! +!!!! Output: +!!!! dhail hail diameter in mm +!!!! 1st-5th rank-ordered hail diameters returned +!!!! +!!!! 13 Aug 2013 .................................Becky Adams-Selin AER +!!!! adapted from hailstone subroutine in SPC's HAILCAST +!!!! 18 Mar 2014 .................................Becky Adams-Selin AER +!!!! added variable rime layer density, per Ziegler et al. (1983) +!!!! 4 Jun 2014 ..................................Becky Adams-Selin AER +!!!! removed initial embryo size dependency on microphysics scheme +!!!! 5 Jun 2014 ..................................Becky Adams-Selin AER +!!!! used smaller initial embryo sizes +!!!! 25 Jun 2015..................................Becky Adams-Selin AER +!!!! Significant revamping. Fixed units bug in HEATBUD that caused +!!!! hailstone temperature instabilities. Similar issue fixed in BREAKUP +!!!! subroutine. Removed graupel from ice content. Changed initial +!!!! embryo size and location to better match literature. Added +!!!! enhanced melting when hailstone collides with liquid water +!!!! in regions above freezing. Final diameter returned is ice diameter +!!!! only. Added hailstone temperature averaging over previous timesteps +!!!! to decrease initial temperature instability at small embyro diameters. +!!!! 3 Sep 2015...................................Becky Adams-Selin AER +!!!! Insert embryos at -13C; interpret pressure and other variables to +!!!! that exact temperature level. +!!!! 16 Nov 2015...................................Becky Adams-Selin AER +!!!! Hailstone travels horizontally through updraft instead of being +!!!! locked in the center. +!!!! 9 Jul 2016...................................Becky Adams-Selin AER +!!!! Uses an adiabatic liquid cloud water profile generated from +!!!! the saturation vapor pressure using the model temperature +!!!! profile. +!!!! 04 Feb 2017...................................Becky Adams-Selin AER +!!!! Added adaptive time-stepping option. +!!!! ***Don't set any higher than 60 seconds*** +!!!! 06 May 2017...................................Becky Adams-Selin AER +!!!! Bug fixes for some memory issues in the adiabatic liquid +!!!! water profile code. +!!!! +!!!! See Adams-Selin and Ziegler 2016, MWR for further documentation. +!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& + RA, qi1d,qc1d,qr1d,qs1d,qg1d, & + VUU, wdur, & + nz,dhail1,dhail2,dhail3,dhail4, & + dhail5 ) + IMPLICIT NONE + INTEGER, INTENT(IN ) :: nz, ii,jj + + REAL, DIMENSION( nz ), & + INTENT(IN ) :: TCA & ! temperature (K) + , rho1d & + , h1d & + , PA & ! pressure (Pa) + , RA & ! vapor mixing ratio (kg/kg) + , VUU & ! updraft speed (m/s) + , qi1d,qc1d,qr1d & + , qs1d,qg1d + + REAL, INTENT(IN ) :: ht & + , wdur + + !Output: 1st-5th rank-ordered hail diameters returned + REAL, INTENT(INOUT) :: dhail1 & ! hail diameter (mm); + , dhail2 & + , dhail3 & + , dhail4 & + , dhail5 + !Local variables + REAL ZBAS, TBAS, WBASP ! height, temp, pressure of cloud base + REAL RBAS ! mix ratio of cloud base + REAL cwitot ! total cloud water, ice mix ratio + INTEGER KBAS ! k of cloud base + REAL tk_embryo ! temperature at which initial embryo is inserted + REAL ZFZL, TFZL, WFZLP ! height, temp, pressure of embryo start point + REAL RFZL ! mix ratio of embryo start point + REAL VUFZL, DENSAFZL ! updraft speed, density of embryo start point + INTEGER KFZL ! k of embryo start point + INTEGER nofroze ! keeps track if hailstone has ever been frozen + INTEGER CLOUDON ! use to zero out cloud water, ice once past updraft duration + REAL RTIME ! real updraft duration (sec) + REAL TAU, TAU_1, TAU_2 ! upper time limit of simulation (sec) + REAL delTAU ! difference between TAU_2 and TAU_1 (sec) + REAL g ! gravity (m/s) + REAL r_d ! constant + !hailstone parameters + REAL*8 DD, D, D_ICE ! hail diameter (m) + REAL VT ! terminal velocity (m/s) + REAL V ! actual stone velocity (m/s) + REAL TS ! hailstone temperature (K) + !HAILSTONE temperature differencing + REAL TSm1, TSm2 ! hailstone temperature at previous 3 timesteps + REAL FW ! fraction of stone that is liquid + REAL WATER ! mass of stone that is liquid + REAL CRIT ! critical water mass allowed on stone surface + REAL DENSE ! hailstone density (kg/m3) + INTEGER ITYPE ! wet (2) or dry (1) growth regime + !1-d column arrays of updraft parameters + REAL, DIMENSION( nz ) :: & + RIA, & ! frozen content mix ratio (kg/kg) + RWA, & ! liquid content mix ratio (kg/kg) + VUU_pert ! perturbed updraft profile (m/s) + !1-d column arrays of updraft characteristics for adiabatic water profile + REAL, DIMENSION( nz ) :: & + RWA_adiabat, & ! adiabatic liquid content mixing ratio (kg/kg) + RWA_new, & + ESVA, & ! saturation vapor pressure + RSA ! saturation mixing ratio + !in-cloud updraft parameters at location of hailstone + REAL P ! in-cloud pressure (Pa) + REAL RS ! in-cloud saturation mixing ratio + REAL RI, RW ! ice, liquid water mix. ratio (kg/kg) + REAL XI, XW ! ice, liquid water content (kg/m3 air) + REAL PC ! in-cloud fraction of frozen water + REAL TC ! in-cloud temperature (K) + REAL VU ! in-cloud updraft speed (m/s) + REAL VUMAX ! in-cloud updraft speed read from WRF (max allowed) + REAL VUCORE ! perturbed in-cloud updraft speed + REAL DENSA ! in-cloud updraft density (kg/m3) + REAL Z ! height of hailstone (m) + REAL DELRW ! diff in sat vap. dens. between hail and air (kg/m3) + REAL, DIMENSION(5) :: dhails !hail diameters from 5 different embryo size + !mean sub-cloud layer variables + REAL TLAYER,RLAYER,PLAYER ! mean sub-cloud temp, mix ratio, pres + REAL TSUM,RSUM,PSUM ! sub-cloud layer T, R, P sums + REAL LDEPTH ! layer depth + !internal function variables + REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE + REAL dum, icefactor, adiabatic_factor + + REAL sec, secdel ! time step, increment in seconds + INTEGER i, j, k, IFOUT, ind(1) + CHARACTER*256 :: message + + secdel = 5.0 + g=9.81 + r_d = 287. + + !!!!!!!!!!!!!!!! 1. UPDRAFT PROPERTIES !!!!!!!!!!!!!!!!!!!!!!! + !!! DEFINE UPDRAFT PROPERTIES !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Upper limit of simulation in seconds + TAU = 7200. !simulation ends + + !Set initial updraft strength - you could reduce to simulate the embryo + ! hovering around the edges of the updraft, as in Heymsfield and Musil (1982) + DO i=1,nz + VUU_pert(i) = VUU(i) * 1. + ENDDO + + +! Initialize diameters to 0. + DO i=1,5 + dhails(i) = 0. + ENDDO + +! Cap updraft lifetime at 2000 sec. + IF (wdur .GT. 2000) THEN + RTIME = 2000. + ELSE + RTIME = wdur + ENDIF + + + !!!!!!!!!!!!!!!! 2. INITIAL EMBRYO !!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! FIND PROPERTIES OF INITIAL EMBRYO LOCATION !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !Find the cloud base for end-of-algorithm purposes. + KBAS=nz + !KFZL=nz + DO k=1,nz + cwitot = qi1d(k) + qc1d(k) + !No longer include graupel in in-cloud ice amounts + !RIA(k) = qi1d(k) + qs1d(k) + qg1d(k) + RIA(k) = qi1d(k) + qs1d(k) + !RWA(k) = qc1d(k) + qr1d(k) + RWA(k) = qc1d(k) + !IF ((RIA(k) .ge. 0.0001) .and. (TCA(k).lt.260.155) .and. & + ! (k .lt. KFZL)) THEN + ! KFZL = k + !ENDIF + IF ((cwitot .ge. 1.E-12) .and. (k .lt. KBAS)) THEN + KBAS = k + ENDIF + ENDDO + !QC - our embryo can't start below the cloud base. + !IF (KFZL .lt. KBAS) THEN + ! KFZL = KBAS + !ENDIF + + !Pull heights, etc. of these levels out of 1-d arrays. + ZBAS = h1d(KBAS) + TBAS = TCA(KBAS) + WBASP = PA(KBAS) + RBAS = RA(KBAS) + + !At coarser resolutions WRF underestimates liquid water aloft. + !A fairer estimate is of the adiabatic liquid water profile, or + !the difference between the saturation mixing ratio at each level + !and the cloud base water vapor mixing ratio + DO k=1,nz + RWA_new(k) = RWA(k) + IF (k.LT.KBAS) THEN + RWA_adiabat(k) = 0. + CYCLE + ENDIF + !saturation vapor pressure + !ESVA(k) = 6.112*exp((2.453E6/461)*(1./273. - 1./TCA(k))) !hPa + ESVA(k) = 611.2*exp(17.67*(TCA(k)-273.155)/(TCA(k)-29.655)) !Pa + !saturation vapor mixing ratio + RSA(k) = 0.62197 * ESVA(k) / (PA(k) - ESVA(k)) + !Up above -31, start converting to ice, entirely by -38 + !(Rosenfeld and Woodley 2000) + IF (TCA(k).gt.242.155) THEN + icefactor = 1. + ELSE IF ((TCA(k).LE.242.155).AND.(TCA(k).GT.235.155)) THEN + icefactor = (1-(242.155-TCA(k))/5.) + ELSE + icefactor = 0. + ENDIF + !Don't want any negative liquid water values + IF (RBAS.GT.RSA(k)) THEN + RWA_adiabat(k) = (RBAS - RSA(k))*icefactor + ELSE + RWA_adiabat(k) = RWA(k) + ENDIF + !Remove cloud liquid water outputted at previous lower levels + ! -- bug fix 170506 -- ! + IF (k.eq.KBAS) THEN + RWA_new(k) = RWA_adiabat(k) + ELSE IF ((k.ge.KBAS+1).AND.(RWA_adiabat(k).ge.1.E-12)) THEN + RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1) + IF (RWA_new(k).LT.0) RWA_new(k) = 0. + ENDIF + ! - old code - ! + !IF (k.eq.KBAS+1) THEN + ! RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) + !ELSE IF ((k.ge.KBAS+2).AND.(RWA_adiabat(k).ge.1.E-12)) THEN + ! RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1) + ! IF (RWA_new(k).LT.0) RWA_new(k) = 0. + !ENDIF + ENDDO + !remove the height factor from RWA_new + DO k=KBAS,nz + RWA_new(k) = RWA_new(k) / (h1d(k)-h1d(k-1)) + ENDDO + + + + !!!!!!!!!!!!!!!! 3. INITIAL EMBRYO SIZE !!!!!!!!!!!!!!!!!!!!! + !!! SET CONSTANT RANGE OF INITIAL EMBRYO SIZES !!! + !!! START LOOP !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! See Adams-Selin and Ziegler 2016 MWR for explanation of why + ! these sizes were picked. + !Run each initial embryo size perturbation + DO i=1,5 + SELECT CASE (i) + CASE (1) + !Initial hail embryo diameter in m, at cloud base + DD = 5.E-3 + tk_embryo = 265.155 !-8C + CASE (2) + DD = 7.5E-3 + tk_embryo = 265.155 !-8C + CASE (3) + DD = 5.E-3 + tk_embryo = 260.155 !-13C + CASE (4) + DD = 7.5E-3 + tk_embryo = 260.155 !-13C + CASE (5) + tk_embryo = 260.155 !-13C + DD = 1.E-2 + END SELECT + + + RTIME = 2000. + IF (wdur .LT. RTIME) RTIME = wdur + + TFZL = tk_embryo + CALL INTERPP(PA, WFZLP, TCA, tk_embryo, IFOUT, nz) + CALL INTERP(h1d, ZFZL, WFZLP, IFOUT, PA, nz) + CALL INTERP(RA, RFZL, WFZLP, IFOUT, PA, nz) + CALL INTERP(VUU_pert, VUFZL, WFZLP, IFOUT, PA, nz) + CALL INTERP(rho1d, DENSAFZL, WFZLP, IFOUT, PA, nz) + + !Begin hail simulation time (seconds) + sec = 0. + + + !!!!!!!!!!!!!!!!!! 4. INITIAL PARAMETERS !!!!!!!!!!!!!!!!! + !!! PUT INITIAL PARAMETERS IN VARIABLES !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !Set initial values for parameters at freezing level + P = WFZLP + RS = RFZL + TC = TFZL + VU = VUFZL + Z = ZFZL - ht + LDEPTH = Z + DENSA = DENSAFZL + + !Set initial hailstone parameters + nofroze=1 !Set test for embryo: 0 for never been frozen; 1 frozen + TS = TC + TSm1 = TS + TSm2 = TS + D = DD !hailstone diameter in m + FW = 0.0 + DENSE = 500. !kg/m3 + ITYPE=1. !Assume starts in dry growth. + CLOUDON=1 !we'll eventually turn cloud "off" once updraft past time limit + + !Start time loop. + DO WHILE (sec .lt. TAU) + sec = sec + secdel + + !!!!!!!!!!!!!!!!!! 5. CALCULATE PARAMETERS !!!!!!!!!!!!!!!!! + !!! CALCULATE UPDRAFT PROPERTIES !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !Intepolate vertical velocity to our new pressure level + CALL INTERP(VUU_pert,VUMAX,P,IFOUT,PA,nz) + !print *, 'INTERP VU: ', VU, P + + !Outside pressure levels? If so, exit loop + IF (IFOUT.EQ.1) GOTO 100 + + !Sine wave multiplier on updraft strength + IF (SEC .GT. 0.0 .AND. SEC .LT. RTIME) THEN + VUCORE = VUMAX * (0.5*SIN((3.14159*SEC)/(RTIME))+0.5)*1.2 + VU = VUCORE + ELSEIF (SEC .GE. RTIME) THEN + VU = 0.0 + CLOUDON = 0 + ENDIF + + !Calculate terminal velocity of the hailstone + ! (use previous values) + CALL TERMINL(DENSA,DENSE,D,VT,TC) + + !Actual velocity of hailstone (upwards positive) + V = VU - VT + + !Use hydrostatic eq'n to calc height of next level + P = P - DENSA*g*V*secdel + Z = Z + V*secdel + + !Interpolate cloud temp, qvapor at new p-level + CALL INTERP(TCA,TC,P,IFOUT,PA,nz) + CALL INTERP(RA,RS,P,IFOUT,PA,nz) + + !New density of in-cloud air + DENSA=P/(r_d*(1.+0.609*RS/(1.+RS))*TC) + + !Interpolate liquid, frozen water mix ratio at new level + CALL INTERP(RIA,RI,P,IFOUT,PA,nz) + CALL INTERP(RWA_new,RW,P,IFOUT,PA,nz) + XI = RI * DENSA * CLOUDON + XW = RW * DENSA * CLOUDON + IF( (XW+XI).GT.0) THEN + PC = XI / (XW+XI) + ELSE + PC = 1. + ENDIF + + + ! SATURATION VAPOUR DENSITY DIFFERENCE BETWTEEN STONE AND CLOUD + CALL VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) + + + !!!!!!!!!!!!!!!!!! 6. STONE'S MASS GROWTH !!!!!!!!!!!!!!!!!!!! + CALL MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& + TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,secdel,ITYPE,DELRW) + + + !!!!!!!!!!!!!!!!!! 7. HEAT BUDGET OF HAILSTONE !!!!!!!!!!!!!!! + CALL HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & + DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,secdel,ITYPE,P) + + + !!!!! 8. TEST DIAMETER OF STONE AND HEIGHT ABOVE GROUND !!!!!!! + !!! TEST IF DIAMETER OF STONE IS GREATER THAN CRITICAL LIMIT, IF SO + !!! BREAK UP + WATER=FW*GM !KG + ! CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S SURFACE + CRIT = 2.0E-4 + IF (WATER.GT.CRIT)THEN + CALL BREAKUP(DENSE,D,GM,FW,CRIT) + ENDIF + + !!! Has stone reached below cloud base? + !IF (Z .LE. 0) GOTO 200 + IF (Z .LE. ZBAS) then + GOTO 200 + endif + + !calculate ice-only diameter size + D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 + + !Has the stone entirely melted and it's below the freezing level? + IF ((D_ICE .LT. 1.E-8) .AND. (TC.GT.273.155)) GOTO 300 + + !move values to previous timestep value + TSm1 = TS + TSm2 = TSm1 + + ENDDO !end cloud lifetime loop + +100 CONTINUE !outside pressure levels in model +200 CONTINUE !stone reached surface +300 CONTINUE !stone has entirely melted and is below freezing level + + !!!!!!!!!!!!!!!!!! 9. MELT STONE BELOW CLOUD !!!!!!!!!!!!!!!!!!!! + !Did the stone shoot out the top of the storm? + !Then let's assume it's lost in the murky "outside storm" world. + IF (P.lt.PA(nz)) THEN + D=0.0 + !Is the stone entirely water? Then set D=0 and exit. + ELSE IF(ABS(FW - 1.0).LT.0.001) THEN + D=0.0 + ELSE IF (Z.GT.0) THEN + !If still frozen, then use melt routine to melt below cloud + ! based on mean below-cloud conditions. + + !Calculate mean sub-cloud layer conditions + TSUM = 0. + RSUM = 0. + PSUM = 0. + DO k=1,KBAS + TSUM = TSUM + TCA(k) + PSUM = PSUM + PA(k) + RSUM = RSUM + RA(k) + ENDDO + TLAYER = TSUM / KBAS + PLAYER = PSUM / KBAS + RLAYER = RSUM / KBAS + + !MELT is expecting a hailstone of only ice. At the surface + !we're only interested in the actual ice diameter of the hailstone, + !so let's shed any excess water now. + D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 + D = D_ICE + CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) + + ENDIF !end check for melting call + + !Check to make sure embryo didn’t just immediately fall out + IF (sec .LT. 60) D=0.0 + + !Check to make sure something hasn't gone horribly wrong + 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 + + ENDDO !end embryo size loop + + + dhail1 = dhails(1) + dhail2 = dhails(2) + dhail3 = dhails(3) + dhail4 = dhails(4) + dhail5 = dhails(5) + + !Testing diagnostics that can be commented out, but please keep + if (0.eq.1) then + print*,'jmh KBAS,ZBAS,TBAS',ii,jj,KBAS,ZBAS,TBAS + print*,'jmh WBASP,RBAS=',ii,jj,WBASP,RBAS + print*,'jmh embryo sizes: ',ii,jj,dhail1,dhail2,dhail3,dhail4,dhail5 + if (ii.eq. 847 .and. jj.eq.394) then + print*,'jmh i=847,j=394 TCA:',ii,jj,TCA(1:nz) + print*,'jmh i=847,j=394 h1d:',ii,jj,h1d(1:nz) + print*,'jmh i=847,j=394 ht:',ii,jj,ht + print*,'jmh i=847,j=394 PA:',ii,jj,PA(1:nz) + print*,'jmh i=847,j=394 rho1d:',ii,jj,rho1d(1:nz) + print*,'jmh i=847,j=394 RA:',ii,jj,RA(1:nz) + print*,'jmh i=847,j=394 qi1d:',ii,jj,qi1d(1:nz) + print*,'jmh i=847,j=394 qc1d:',ii,jj,qc1d(1:nz) + print*,'jmh i=847,j=394 qr1d:',ii,jj,qr1d(1:nz) + print*,'jmh i=847,j=394 qs1d:',ii,jj,qs1d(1:nz) + print*,'jmh i=847,j=394 qg1d:',ii,jj,qg1d(1:nz) + print*,'jmh i=847,j=394 VUU:',ii,jj,VUU(1:nz) + print*,'jmh i=847,j=394 wdur:',ii,jj,wdur + print*,'jmh i=847,j=394 nz:',ii,jj,nz + endif + if (KBAS .gt. 30) then + print*,'jmh KBAS>30 TCA:',ii,jj,TCA(1:nz) + print*,'jmh KBAS>30 h1d:',ii,jj,h1d(1:nz) + print*,'jmh KBAS>30 ht:',ii,jj,ht + print*,'jmh KBAS>30 PA:',ii,jj,PA(1:nz) + print*,'jmh KBAS>30 rho1d:',ii,jj,rho1d(1:nz) + print*,'jmh KBAS>30 RA:',ii,jj,RA(1:nz) + print*,'jmh KBAS>30 qi1d:',ii,jj,qi1d(1:nz) + print*,'jmh KBAS>30 qc1d:',ii,jj,qc1d(1:nz) + print*,'jmh KBAS>30 qr1d:',ii,jj,qr1d(1:nz) + print*,'jmh KBAS>30 qs1d:',ii,jj,qs1d(1:nz) + print*,'jmh KBAS>30 qg1d:',ii,jj,qg1d(1:nz) + print*,'jmh KBAS>30 VUU:',ii,jj,VUU(1:nz) + print*,'jmh KBAS>30 wdur:',ii,jj,wdur + print*,'jmh KBAS>30 nz:',ii,jj,nz + endif + print*,'jmh ------end of hailstone driver------- jmh' + endif + END SUBROUTINE hailstone_driver + + + SUBROUTINE INTERPP(PA,PVAL,TA,TVAL,IFOUT,ITEL) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!! + !!!! INTERP: to linearly interpolate value of pval at temperature tval + !!!! between two levels of pressure array pa and temperatures ta + !!!! + !!!! INPUT: PA 1D array of pressure, to be interpolated + !!!! TA 1D array of temperature + !!!! TVAL temperature value at which we want to calculate pressure + !!!! IFOUT set to 0 if TVAL outside range of TA + !!!! ITEL number of vertical levels + !!!! OUTPUT: PVAL interpolated pressure variable at temperature tval + !!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IMPLICIT NONE + + REAL PVAL, TVAL + REAL, DIMENSION( ITEL) :: TA, PA + INTEGER ITEL, IFOUT + !local variables + INTEGER I + REAL FRACT + + IFOUT=1 + + DO I=1,ITEL-1 + IF ( (TVAL .LT. TA(I) .AND. TVAL .GE. TA(I+1)) .or. & ! dT/dz < 0 + (TVAL .GT. TA(I) .AND. TVAL .LE. TA(I+1)) ) THEN ! dT/dz > 0 + + FRACT = (TA(I) - TVAL) / (TA(I) - TA(I+1)) + !.... compute the pressure value pval at temperature tval + PVAL = ((1.0 - FRACT) * PA(I)) + (FRACT * PA(I+1)) + + !End loop + IFOUT=0 + EXIT + ENDIF + ENDDO + + END SUBROUTINE INTERPP + + + + SUBROUTINE INTERP(AA,A,P,IFOUT,PA,ITEL) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!! + !!!! INTERP: to linearly interpolate values of A at level P + !!!! between two levels of AA (at levels PA) + !!!! + !!!! INPUT: AA 1D array of variable + !!!! PA 1D array of pressure + !!!! P new pressure level we want to calculate A at + !!!! IFOUT set to 0 if P outside range of PA + !!!! ITEL number of vertical levels + !!!! OUTPUT: A variable at pressure level P + !!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IMPLICIT NONE + + REAL A, P + REAL, DIMENSION( ITEL) :: AA, PA + INTEGER ITEL, IFOUT + !local variables + INTEGER I + REAL PDIFF, VDIFF, RDIFF, VERH, ADIFF + + IFOUT=1 + + DO I=1,ITEL-1 + IF (P.LE.PA(I) .AND. P.GT.PA(I+1)) THEN + !Calculate ratio between vdiff and pdiff + PDIFF = PA(I)-PA(I+1) + VDIFF = PA(I)-P + VERH = VDIFF/PDIFF + + !Calculate the difference between the 2 A values + RDIFF = AA(I+1) - AA(I) + + !Calculate new value + A = AA(I) + RDIFF*VERH + + !End loop + IFOUT=0 + EXIT + ENDIF + ENDDO + + END SUBROUTINE INTERP + + + SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!! + !!!! INTERP: Calculate terminal velocity of the hailstone + !!!! + !!!! INPUT: DENSA density of updraft air (kg/m3) + !!!! DENSE density of hailstone + !!!! D diameter of hailstone (m) + !!!! TC updraft temperature (K) + !!!! OUTPUT:VT hailstone terminal velocity (m/s) + !!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IMPLICIT NONE + + REAL*8 D + REAL DENSA, DENSE, TC, VT + REAL GMASS, GX, RE, W, Y + REAL, PARAMETER :: PI = 3.141592654, G = 9.78956 + REAL ANU + + !Mass of stone in kg + GMASS = (DENSE * PI * (D**3.)) / 6. + + !Dynamic viscosity + ANU = (0.00001718)*(273.155+120.)/(TC+120.)*(TC/273.155)**(1.5) + + !CALC THE BEST NUMBER, X AND REYNOLDS NUMBER, RE + GX=(8.0*GMASS*G*DENSA)/(PI*(ANU*ANU)) + RE=(GX/0.6)**0.5 + + !SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON + !THE BEST NUMBER + IF (GX.LT.550) THEN + W=LOG10(GX) + Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0) + RE=10**Y + VT=ANU*RE/(D*DENSA) + ELSE IF (GX.GE.550.AND.GX.LT.1800) THEN + W=LOG10(GX) + Y= -1.81391 + 1.34671*W - 0.12427*(W**2.0) + 0.0063*(W**3.0) + RE=10**Y + VT=ANU*RE/(D*DENSA) + ELSE IF (GX.GE.1800.AND.GX.LT.3.45E08) THEN + RE=0.4487*(GX**0.5536) + VT=ANU*RE/(D*DENSA) + ELSE + RE=(GX/0.6)**0.5 + VT=ANU*RE/(D*DENSA) + ENDIF + + END SUBROUTINE TERMINL + + + + SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY + !!! BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD + !!! AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT, + !!! AND IF THE STONE IS IN WET OR DRY GROWTH REGIME + !!! + !!! INPUT: PC fraction of updraft water that is frozen + !!! TS temperature of hailstone (K) + !!! TC temperature of updraft air (K) + !!! ITYPE wet (2) or dry (1) growth regime + !!! OUTPUT: DELRW difference in sat vap. dens. between hail and air + !!! (kg/m3) + !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + IMPLICIT NONE + REAL DELRW, PC, TS, TC + INTEGER ITYPE + !local variables + REAL RV, ALV, ALS, RATIO + DATA RV/461.48/,ALV/2500000./,ALS/2836050./ + REAL ESAT, RHOKOR, ESATW, RHOOMGW, ESATI, RHOOMGI, RHOOMG + + !!! FOR HAILSTONE: FIRST TEST IF STONE IS IN WET OR DRY GROWTH + RATIO = 1./273.155 + IF(ITYPE.EQ.2) THEN !!WET GROWTH + ESAT=611.*EXP(ALV/RV*(RATIO-1./TS)) + ELSE !!DRY GROWTH + ESAT=611.*EXP(ALS/RV*(RATIO-1./TS)) + ENDIF + RHOKOR=ESAT/(RV*TS) + + !!! NOW FOR THE AMBIENT/IN-CLOUD CONDITIONS + ESATW=611.*EXP(ALV/RV*(RATIO-1./TC)) + RHOOMGW=ESATW/(RV*TC) + ESATI=611.*EXP(ALS/RV*(RATIO-1./TC)) + RHOOMGI=ESATI/(RV*TC) + !RHOOMG=PC*(RHOOMGI-RHOOMGW)+RHOOMGW + RHOOMG = RHOOMGI !done as in hailtraj.f + + !!! CALC THE DIFFERENCE(KG/M3): <0 FOR CONDENSATION, + !!! >0 FOR EVAPORATION + DELRW=(RHOKOR-RHOOMG) + + END SUBROUTINE VAPORCLOSE + + + + SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& + TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,ITYPE,DELRW) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! CALC THE STONE'S INCREASE IN MASS + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + IMPLICIT NONE + REAL*8 D + REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DI,ANU,RE,AE, & + TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,DELRW + INTEGER ITYPE + !local variables + REAL PI, D0, GMW2, GMI2, EW, EI,DGMV + REAL DENSEL, DENSELI, DENSELW + REAL DC !MEAN CLOUD DROPLET DIAMETER (MICRONS, 1E-6M) + REAL VOLL, VOLT !VOLUME OF NEW LAYER, TOTAL (M3) + REAL VOL1, DGMW_NOSOAK, SOAK, SOAKM + REAL DENSAC, E, AFACTOR, NS, TSCELSIUS, VIMP, W + PI=3.141592654 + + !!! CALCULATE THE DIFFUSIVITY DI (m2/s) + D0=0.226*1.E-4 ! change to m2/s, not cm2/s + DI=D0*(TC/273.155)**1.81*(100000./P) + + !!! COLLECTION EFFICIENCY FOR WATER AND ICE + !EW=1.0 + !!!! IF TS WARMER THAN -5C THEN ACCRETE ALL THE ICE (EI=1.0) + !!!! OTHERWISE EI=0.21 + !IF(TS.GE.268.15)THEN + ! EI=1.00 + !ELSE + ! EI=0.21 + !ENDIF + + !!! COLLECTION EFFICIENCY FOR WATER AND ICE + EW=1.0 + !!! Linear function for ice accretion efficiency + IF (TC .GE. 273.155) THEN + EI=1.00 + ELSE IF (TC.GE.233.155) THEN + EI= 1.0 - ( (273.155 - TS) / 40. ) + ELSE !cooler than -40C + EI = 0.0 + ENDIF + + !!! CALCULATE THE VENTILATION COEFFICIENT - NEEDED FOR GROWTH FROM VAPOR + !The coefficients in the ventilation coefficient equations have been + !experimentally derived, and are expecting cal-C-g units. Do some conversions. + DENSAC = DENSA * (1.E3) * (1.E-6) + !dynamic viscosity + ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5 + !!! CALCULATE THE REYNOLDS NUMBER - unitless + RE=D*VT*DENSAC/ANU + E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv) + !!! SELECT APPROPRIATE VALUES OF AE ACCORDING TO RE + IF(RE.LT.6000.0)THEN + AE=0.78+0.308*E + ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN + AE=0.76*E + ELSEIF(RE.GE.20000.0) THEN + AE=(0.57+9.0E-6*RE)*E + ENDIF + + + !!! CALC HAILSTONE'S MASS (GM), MASS OF WATER (GMW) AND THE + !!! MASS OF ICE IN THE STONE (GMI) + GM=PI/6.*(D**3.)*DENSE + GMW=FW*GM + GMI=GM-GMW + + !!! STORE THE MASS + GM1=GM + + !!! NEW MASS GROWTH CALCULATIONS WITH VARIABLE RIME + !!! LAYER DENSITY BASED ON ZIEGLER ET AL. (1983) + + !!! CALCULATE INCREASE IN MASS DUE INTERCEPTED CLD WATER, USE + !!! ORIGINAL DIAMETER + GMW2=GMW+SEKDEL*(PI/4.*D**2.*VT*XW*EW) + DGMW=GMW2-GMW + GMW=GMW2 + + !!! CALCULATE THE INCREASE IN MASS DUE INTERCEPTED CLOUD ICE + GMI2=GMI+SEKDEL*(PI/4.*D**2.*VT*XI*EI) + DGMI=GMI2-GMI + GMI=GMI2 + + !!! CALCULATE INCREASE IN MASS DUE TO SUBLIMATION/CONDENSATION OF + !!! WATER VAPOR + DGMV = SEKDEL*2*PI*D*AE*DI*DELRW + IF (DGMV .LT. 0) DGMV=0 + + !!! CALCULATE THE TOTAL MASS CHANGE + DGM=DGMW+DGMI+DGMV + !Dummy arguments in the event of no water, vapor, or ice growth + DENSELW = 900. + DENSELI = 700. + !!! CALCULATE DENSITY OF NEW LAYER, DEPENDS ON FW AND ITYPE + IF (ITYPE.EQ.1) THEN !DRY GROWTH + !If hailstone encountered supercooled water, calculate new layer density + ! using Macklin form + IF ((DGMW.GT.0).OR.(DGMV.GT.0)) THEN + !MEAN CLOUD DROPLET RADIUS, ASSUME CLOUD DROPLET CONC OF 3E8 M-3 (300 CM-3) + DC = (0.74*XW / (3.14159*1000.*3.E8))**0.33333333 * 1.E6 !MICRONS + !!! FIND THE STOKES NUMBER (rasmussen heymsfield 1985) + NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50) !need hail radius in cm + !!! FIND IMPACT VELOCITY (rasmussen heymsfield 1985) + IF (NS.LT.0.1)THEN + W=-1. + ELSE + W = LOG10(NS) + ENDIF + IF (RE.GT.200) THEN + IF (NS.LT.0.1) THEN + VIMP = 0.0 + ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN + VIMP = (0.356 + 0.4738*W - 0.1233*W**2. & + -0.1618*W**3. + 0.0807*W**4.)*VT + ELSEIF (NS.GT.10) THEN + VIMP = 0.63*VT + ENDIF + ELSEIF ((RE.GT.65).AND.(RE.LE.200)) THEN + IF (NS.LT.0.1) THEN + VIMP = 0.0 + ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN + VIMP = (0.3272 + 0.4907*W - 0.09452*W**2. & + -0.1906*W**3. + 0.07105*W**4.)*VT + ELSEIF (NS.GT.10) THEN + VIMP = 0.61*VT + ENDIF + ELSEIF ((RE.GT.20).AND.(RE.LE.65)) THEN + IF (NS.LT.0.1) THEN + VIMP = 0.0 + ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN + VIMP = (0.2927 + 0.5085*W - 0.03453*W**2. & + -0.2184*W**3. + 0.03595*W**4.)*VT + ELSEIF (NS.GT.10) THEN + VIMP = 0.59*VT + ENDIF + ELSEIF (RE.LE.20) THEN + IF (NS.LT.0.4) THEN + VIMP = 0.0 + ELSEIF ((NS.GE.0.4).AND.(NS.LE.10)) THEN + VIMP = (0.1701 + 0.7246*W + 0.2257*W**2. & + -1.13*W**3. + 0.5756*W**4.)*VT + ELSEIF (NS.GT.10) THEN + VIMP = 0.57*VT + ENDIF + ENDIF + + !RIME LAYER DENSITY, HEYMSFIELD AND PFLAUM 1985 FORM + TSCELSIUS = TS - 273.16 + AFACTOR = -DC*VIMP/TSCELSIUS + IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN + DENSELW = 0.30*(AFACTOR)**0.44 + ELSEIF (TSCELSIUS.GT.-5.) THEN + DENSELW = EXP(-0.03115 - 1.7030*AFACTOR + & + 0.9116*AFACTOR**2. - 0.1224*AFACTOR**3.) + ENDIF + + DENSELW = DENSELW * 1000. !KG M-3 + !BOUND POSSIBLE DENSITIES + IF (DENSELW.LT.100) DENSELW=100 + IF (DENSELW.GT.900) DENSELW=900 + !WRITE(12,*) 'MASSAGR, PFLAUM, MACKLIN: ', DENSELW, & + ! MACKLIN_DENSELW + ENDIF + IF (DGMI.GT.0) THEN + !Ice collection main source of growth, so set new density layer + DENSELI = 700. + ENDIF + + !All liquid water contributes to growth, none is soaked into center. + DGMW_NOSOAK = DGMW !All liquid water contributes to growth, + ! none of it is soaked into center. + + ELSE !WET GROWTH + !Collected liquid water can soak into the stone before freezing, + ! increasing mass and density but leaving volume constant. + !Volume of current drop, before growth + VOL1 = GM/DENSE + !Difference b/w mass of stone if density is 900 kg/m3, and + ! current mass + SOAK = 900*VOL1 - GM + !Liquid mass available + SOAKM = DGMW + !Soak up as much liquid as we can, up to a density of 900 kg/m3 + IF (SOAKM.GT.SOAK) SOAKM=SOAK + GM = GM+SOAKM !Mass of current drop, plus soaking + !New density of current drop, including soaking but before growth + DENSE = GM/VOL1 + !Mass increment of liquid water growth that doesn't + ! include the liquid water we just soaked into the stone. + DGMW_NOSOAK = DGMW - SOAKM + + !Whatever growth does occur has high density + DENSELW = 900. !KG M-3 + DENSELI = 900. + + ENDIF + + !!!VOLUME OF NEW LAYER + !VOLL = (DGM) / DENSEL + !VOLL = (DGMI+DGMV+DGMW_NOSOAK) / DENSEL + !VOLL = (DGMI) / DENSELI + (DGMW_NOSOAK+DGMV) / DENSELW + IF (DGMI.LE.0) THEN + VOLL = (DGMW_NOSOAK+DGMV) / DENSELW + ELSE IF (DGMW.LE.0) THEN + VOLL = (DGMI) / DENSELI + ELSE + VOLL = (DGMI) / DENSELI + (DGMW_NOSOAK+DGMV) / DENSELW + ENDIF + + !!!NEW TOTAL VOLUME, DENSITY, DIAMETER + VOLT = VOLL + GM/DENSE + !DENSE = (GM+DGM) / VOLT + DENSE = (GM+DGMI+DGMV+DGMW_NOSOAK) / VOLT + !D=D+SEKDEL*0.5*VT/DENSE*(XW*EW+XI*EI) + GM = GM+DGMI+DGMW_NOSOAK+DGMV + D = ( (6*GM) / (PI*DENSE) )**0.33333333 + + END SUBROUTINE MASSAGR + + + + SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & + DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,ITYPE,P) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! CALCULATE HAILSTONE'S HEAT BUDGET + !!! See Rasmussen and Heymsfield 1987; JAS + !!! Original Hailcast's variable units + !!! TS - Celsius + !!! FW - unitless, between 0 and 1 + !!! TC - Celsius + !!! VT - m/s + !!! D - m + !!! DELRW - g/cm3 (per comment) + !!! DENSA - g/cm3 (per comment) + !!! GM1, DMG, DGMW, DGMV, DGMI, GMW, GMI - should all be kg + !!! DI - cm2 / sec + !!! P - hPa + !!! Original HAILCAST HEATBUD subroutine uses c-g-s units, so do some conversions + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IMPLICIT NONE + REAL*8 D + REAL TS,TSm1,TSm2,FW,TC,VT,DELRW,DENSA,GM1,GM,DGM,DGMW,DGMV, & + DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,P + INTEGER ITYPE + + REAL RV, RD, G, PI, ALF, ALV, ALS, CI, CW, AK + REAL H, AH, TCC, TSC, DELRWC, DENSAC, TDIFF + REAL DMLT + REAL TSCm1, TSCm2 + DATA RV/461.48/,RD/287.04/,G/9.78956/ + DATA PI/3.141592654/,ALF/79.7/,ALV/597.3/ + DATA ALS/677.0/,CI/0.5/,CW/1./ + + !Convert values to non-SI units here + TSC = TS - 273.155 + TSCm1 = TSm1 - 273.155 + TSCm2 = TSm2 - 273.155 + TCC = TC - 273.155 + DELRWC = DELRW * (1.E3) * (1.E-6) + DENSAC = DENSA * (1.E3) * (1.E-6) + !DI still in cm2/sec + + + !!! CALCULATE THE CONSTANTS + AK=(5.8+0.0184*TCC)*1.E-5 !thermal conductivity - cal/(cm*sec*K) + !dynamic viscosity - calculated in MASSAGR + !ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5 + + !!! CALCULATE THE REYNOLDS NUMBER - unitless + !RE=D*VT*DENSAC/ANU - calculated in MASSAGR + + H=(0.71)**(0.333333333)*(RE**0.50) !ventilation coefficient heat (fh) + !E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv) + + !!! SELECT APPROPRIATE VALUES OF AH AND AE ACCORDING TO RE + IF(RE.LT.6000.0)THEN + AH=0.78+0.308*H + !AE=0.78+0.308*E + ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN + AH=0.76*H + !AE=0.76*E + ELSEIF(RE.GE.20000.0) THEN + AH=(0.57+9.0E-6*RE)*H + !AE=(0.57+9.0E-6*RE)*E calculated in MASSAGR + ENDIF + + !!! FOR DRY GROWTH FW=0, CALCULATE NEW TS, ITIPE=1 + !!! FOR WET GROWTH TS=0, CALCULATE NEW FW, ITIPE=2 + + + IF(ITYPE.EQ.1) THEN + !!! DRY GROWTH; CALC NEW TEMP OF THE STONE + !Original Hailcast algorithm (no time differencing) + !TSC=TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* & + ! (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & + ! DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC) + TSC=0.6*(TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* & + (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & + DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)) + & + 0.2*TSCm1 + 0.2*TSCm2 + + TS = TSC+273.155 + IF (TS.GE.273.155) THEN + TS=273.155 + ENDIF + TDIFF = ABS(TS-273.155) + IF (TDIFF.LE.1.E-6) ITYPE=2 !NOW IN WET GROWTH + + ELSE IF (ITYPE.EQ.2) THEN + !!! WET GROWTH; CALC NEW FW + + IF (TCC.LT.0.) THEN + !Original Hailcast algorithm + FW=FW-FW*DGM/GM1+SEKDEL/(GM1*ALF)* & + (2.*PI*D*(AH*AK*TCC-AE*ALV*DI*DELRWC)+ & + DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC) + ELSE + !Calculate decrease in ice mass due to melting + DMLT = (2.*PI*D*AH*AK*TCC + 2.*PI*D*AE*ALV*DI*DELRWC + & + DGMW/SEKDEL*CW*TCC) / ALF + FW = (FW*GM + DMLT) / GM + ENDIF + + IF(FW.GT.1.)FW=1. + IF(FW.LT.0.)FW=0. + + !IF ALL OUR ACCRETED WATER WAS FROZEN, WE ARE BACK IN DRY GROWTH + IF(FW.LE.1.E-6) THEN + ITYPE=1 + ENDIF + + ENDIF + + END SUBROUTINE HEATBUD + + + + SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! TEST IF AMOUNT OF WATER ON SURFACE EXCEEDS CRTICAL LIMIT- + !!! IF SO INVOKE SHEDDING SCHEME + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + IMPLICIT NONE + REAL*8 D + REAL DENSE, GM, FW, CRIT + !local variables + REAL WATER, GMI, WAT, PI + DATA PI/3.141592654/ + + WATER=FW*GM !KG + !GMI=(GM-WATER) !KG + !REMAIN = CRIT*GM + + ! CALC CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S + ! SURFACE + !CRIT=0.268+0.1389*GMI + !CRIT=0.268*1.E-3 + 0.1389*1.E-3*GMI !mass now in kg instead of g + !CRIT = 1.0E-10 + !CRIT - now passed from main subroutine + + WAT=WATER-CRIT + GM=GM-WAT + FW=(CRIT)/GM + + IF(FW.GT.1.0) FW=1.0 + IF(FW.LT.0.0) FW=0.0 + + ! RECALCULATE DIAMETER AFTER SHEDDING + ! Assume density remains the same + D=(6.*GM/(PI*DENSE))**(0.333333333) + END SUBROUTINE BREAKUP + + + SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! This is a spherical hail melting estimate based on the Goyer + !!! et al. (1969) eqn (3). The depth of the warm layer, estimated + !!! terminal velocity, and mean temperature of the warm layer are + !!! used. DRB. 11/17/2003. + !!! + !!! INPUT: TLAYER mean sub-cloud layer temperature (K) + !!! PLAYER mean sub-cloud layer pressure (Pa) + !!! RLAYER mean sub-cloud layer mixing ratio (kg/kg) + !!! VT terminal velocity of stone (m/s) + !!! OUTPUT: D diameter (m) + !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IMPLICIT NONE + + REAL*8 D + REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT + REAL eenv, delta, ewet, de, der, wetold, wetbulb, wetbulbk + REAL tdclayer, tclayer, eps, b, hplayer + REAL*8 a + REAL sd, lt, ka, lf, lv, t0, dv, pi, rv, rhoice, & + tres, re, delt, esenv, rhosenv, essfc, rhosfc, dsig, & + dmdt, mass, massorg, newmass, gamma, r, rho + INTEGER wcnt + + !Convert temp to Celsius, calculate dewpoint in celsius + eps = 0.622 + tclayer = TLAYER - 273.155 + a = 2.53E11 + b = 5.42E3 + tdclayer = b / LOG(a*eps / (rlayer*player)) + hplayer = player / 100. + + !Calculate partial vapor pressure + eenv = (player*rlayer) / (rlayer+eps) + eenv = eenv / 100. !convert to mb + + !Estimate wet bulb temperature (C) + gamma = 6.6E-4*player + delta = (4098.0*eenv)/((tdclayer+237.7)*(tdclayer+237.7)) + wetbulb = ((gamma*tclayer)+(delta*tdclayer))/(gamma+delta) + + !Iterate to get exact wet bulb + wcnt = 0 + DO WHILE (wcnt .lt. 11) + ewet = 6.108*(exp((17.27*wetbulb)/(237.3 + wetbulb))) + de = (0.0006355*hplayer*(tclayer-wetbulb))-(ewet-eenv) + der= (ewet*(.0091379024 - (6106.396/(273.155+wetbulb)**2))) & + - (0.0006355*hplayer) + wetold = wetbulb + wetbulb = wetbulb - de/der + wcnt = wcnt + 1 + IF ((abs(wetbulb-wetold)/wetbulb.gt.0.0001)) THEN + EXIT + ENDIF + ENDDO + + wetbulbk = wetbulb + 273.155 !convert to K + ka = .02 ! thermal conductivity of air + lf = 3.34e5 ! latent heat of melting/fusion + lv = 2.5e6 ! latent heat of vaporization + t0 = 273.155 ! temp of ice/water melting interface + dv = 0.25e-4 ! diffusivity of water vapor (m2/s) + pi = 3.1415927 + rv = 1004. - 287. ! gas constant for water vapor + rhoice = 917.0 ! density of ice (kg/m**3) + r = D/2. ! radius of stone (m) + + !Compute residence time in warm layer + tres = LDEPTH / VT + + !Calculate dmdt based on eqn (3) of Goyer et al. (1969) + !Reynolds number...from pg 317 of Atmo Physics (Salby 1996) + !Just use the density of air at 850 mb...close enough. + rho = 85000./(287.*TLAYER) + re = rho*r*VT*.01/1.7e-5 + + !Temperature difference between environment and hailstone surface + delt = wetbulb !- 0.0 !assume stone surface is at 0C + !wetbulb is in Celsius + + !Difference in vapor density of air stream and equil vapor + !density at the sfc of the hailstone + esenv = 610.8*(exp((17.27*wetbulb)/ & + (237.3 + wetbulb))) ! es environment in Pa + rhosenv = esenv/(rv*wetbulbk) + essfc = 610.8*(exp((17.27*(t0-273.155))/ & + (237.3 + (t0-273.155)))) ! es environment in Pa + rhosfc = essfc/(rv*t0) + dsig = rhosenv - rhosfc + + !Calculate new mass growth + dmdt = (-1.7*pi*r*(re**0.5)/lf)*((ka*delt)+((lv-lf)*dv*dsig)) + IF (dmdt.gt.0.) dmdt = 0 + mass = dmdt*tres + + !Find the new hailstone diameter + massorg = 1.33333333*pi*r*r*r*rhoice + newmass = massorg + mass + if (newmass.lt.0.0) newmass = 0.0 + D = 2.*(0.75*newmass/(pi*rhoice))**0.333333333 + END SUBROUTINE MELT + + +END MODULE module_diag_hailcast + From e1bdb6bc7fdc5a5963bccf6f869dfc4a825e4fb0 Mon Sep 17 00:00:00 2001 From: Yunheng Wang Date: Thu, 10 Mar 2022 23:53:56 +0000 Subject: [PATCH 2/6] Modularized the hailcast module --- driver/fvGFS/fv_nggps_diag.F90 | 443 +++++------------- model/fv_arrays.F90 | 4 +- tools/module_diag_hailcast.F90 | 793 +++++++++++++++++++++------------ 3 files changed, 622 insertions(+), 618 deletions(-) diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index d16a37e86..7d3e09a2f 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -62,9 +62,8 @@ module fv_nggps_diags_mod ! ! - use mpp_mod, only: mpp_pe, mpp_root_pe,FATAL,mpp_error,NOTE, stdlog, input_nml_file + use mpp_mod, only: mpp_pe, mpp_root_pe,FATAL,mpp_error,NOTE use constants_mod, only: grav, rdgas - use fms_mod, only: check_nml_error use time_manager_mod, only: time_type, get_time use diag_manager_mod, only: register_diag_field, send_data use diag_axis_mod, only: get_axis_global_length, get_diag_axis, get_diag_axis_name @@ -75,7 +74,7 @@ module fv_nggps_diags_mod use fv_diagnostics_mod, only: range_check, dbzcalc,max_vv,get_vorticity, & max_uh,max_vorticity,bunkers_vector, & helicity_relative_CAPS,max_vorticity_hy1 - use module_diag_hailcast, only: hailcast_diagnostic_driver + use module_diag_hailcast use fv_arrays_mod, only: fv_atmos_type use mpp_domains_mod, only: domain1d, domainUG #ifdef MULTI_GASES @@ -101,18 +100,11 @@ module fv_nggps_diags_mod 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_hailcast_dhail - integer :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5, id_hailcast_wdur, id_hailcast_wup_mask - integer :: id_hailcast_diam_mean, id_hailcast_diam_std 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 integer :: kstt_ustm,kstt_vstm,kend_ustm,kend_vstm,kstt_srh01 integer :: kstt_srh03,kend_srh01,kend_srh03 - integer :: kstt_hc, kend_hc - integer :: kstt_hc1,kstt_hc2,kstt_hc3,kstt_hc4,kstt_hc5 - integer :: kend_hc1,kend_hc2,kend_hc3,kend_hc4,kend_hc5 - integer :: kstt_hcd, kend_hcd, kstt_hcm, kend_hcm integer :: id_maxvort02,kstt_maxvort02,kend_maxvort02 integer :: isco, ieco, jsco, jeco, npzo, ncnsto integer :: isdo, iedo, jsdo, jedo @@ -132,6 +124,7 @@ module fv_nggps_diags_mod ! file name character(len=64) :: file_name = 'gfs_dyn' + INTEGER :: istatus ! tracers character(len=128) :: tname @@ -146,11 +139,6 @@ module fv_nggps_diags_mod real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01 real, dimension(:,:),allocatable :: maxvorthy1,maxvort02 - real, allocatable :: hailcast_dhail1(:,:), hailcast_dhail2(:,:), hailcast_dhail3(:,:), hailcast_dhail4(:,:), hailcast_dhail5(:,:) !hailstone diameters (mm) - real, allocatable :: hailcast_dhail1_max(:,:), hailcast_dhail2_max(:,:), hailcast_dhail3_max(:,:), hailcast_dhail4_max(:,:), hailcast_dhail5_max(:,:) !hailstone diameters (mm) - real, allocatable :: hailcast_diam_mean(:,:), hailcast_diam_std(:,:) !mean and standard deviation of five hailstones (mm) - real, allocatable :: hailcast_wdur(:,:), hailcast_wup_mask(:,:) !persistent arrays for updraft duration (s) and mask - public :: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg #ifdef use_WRTCOMP public :: fv_dyn_bundle_setup @@ -163,32 +151,6 @@ subroutine fv_nggps_diag_init(Atm, axes, Time) integer, intent(in) :: axes(4) type(time_type), intent(in) :: Time integer :: n, i, j, nz - logical :: do_hailcast = .true. !This controls whether hailcast is used - namelist /fv_diagnostics_nml/ do_hailcast - integer :: ios, ierr - integer :: unit - - !namelist file for hailcast -#ifdef INTERNAL_FILE_NML - ! Read Main namelist - read (input_nml_file,fv_diagnostics_nml,iostat=ios) - ierr = check_nml_error(ios,'fv_diagnostics_nml') -#else - f_unit=open_namelist_file() - rewind (f_unit) - ! Read Main namelist - read (f_unit,fv_diagnostics_nml,iostat=ios) - ierr = check_nml_error(ios,'fv_diagnostics_nml') - call close_file(f_unit) -#endif - - unit = stdlog() - write(unit, nml=fv_diagnostics_nml) - !end hailcast nml - - if (mpp_pe() == mpp_root_pe()) then - print*, 'do_hailcast = ', do_hailcast - end if n = 1 ncnsto = Atm(1)%ncnst @@ -471,93 +433,8 @@ subroutine fv_nggps_diag_init(Atm, axes, Time) endif - !register hailcast arrays - if (do_hailcast) then - id_hailcast_dhail = register_diag_field ( trim(file_name), 'hailcast_dhail_max', axes(1:2), Time, & - 'hourly max hail diameter', 'mm', missing_value=missing_value ) - id_hailcast_dhail1 = register_diag_field ( trim(file_name), 'hailcast_dhail1_max', axes(1:2), Time, & - 'hourly max hail diameter (embryo 1)', 'mm', missing_value=missing_value ) - id_hailcast_dhail2 = register_diag_field ( trim(file_name), 'hailcast_dhail2_max', axes(1:2), Time, & - 'hourly max hail diameter (embryo 2)', 'mm', missing_value=missing_value ) - id_hailcast_dhail3 = register_diag_field ( trim(file_name), 'hailcast_dhail3_max', axes(1:2), Time, & - 'hourly max hail diameter (embryo 3)', 'mm', missing_value=missing_value ) - id_hailcast_dhail4 = register_diag_field ( trim(file_name), 'hailcast_dhail4_max', axes(1:2), Time, & - 'hourly max hail diameter (embryo 4)', 'mm', missing_value=missing_value ) - id_hailcast_dhail5 = register_diag_field ( trim(file_name), 'hailcast_dhail5_max', axes(1:2), Time, & - 'hourly max hail diameter (embryo 5)', 'mm', missing_value=missing_value ) - id_hailcast_diam_mean = register_diag_field ( trim(file_name), 'hailcast_diam_mean', axes(1:2), Time, & - 'mean hail diameter', 'mm', missing_value=missing_value ) - id_hailcast_diam_std = register_diag_field ( trim(file_name), 'hailcast_diam_std', axes(1:2), Time, & - 'standard deviation of hail diameter', 'mm', missing_value=missing_value ) - id_hailcast_wdur = register_diag_field ( trim(file_name), 'hailcast_wdur', axes(1:2), Time, & - 'updraft duration', 's', missing_value=missing_value ) - id_hailcast_wup_mask = register_diag_field ( trim(file_name), 'hailcast_wup_mask', axes(1:2), Time, & - 'updraft mask', '', missing_value=missing_value ) - - if (id_hailcast_dhail > 0) then - kstt_hc = nlevs+1; kend_hc = nlevs+1 - nlevs = nlevs + 1 - endif - - if (id_hailcast_dhail1 > 0) then - kstt_hc1 = nlevs+1; kend_hc1 = nlevs+1 - nlevs = nlevs + 1 - endif - if (id_hailcast_dhail2 > 0) then - kstt_hc2 = nlevs+1; kend_hc2 = nlevs+1 - nlevs = nlevs + 1 - endif - if (id_hailcast_dhail3 > 0) then - kstt_hc3 = nlevs+1; kend_hc3 = nlevs+1 - nlevs = nlevs + 1 - endif - if (id_hailcast_dhail4 > 0) then - kstt_hc4 = nlevs+1; kend_hc4 = nlevs+1 - nlevs = nlevs + 1 - endif - if (id_hailcast_dhail5 > 0) then - kstt_hc5 = nlevs+1; kend_hc5 = nlevs+1 - nlevs = nlevs + 1 - endif - if (id_hailcast_wdur > 0) then - kstt_hcd = nlevs+1; kend_hcd = nlevs+1 - nlevs = nlevs + 1 - endif - if (id_hailcast_wup_mask > 0) then - kstt_hcm = nlevs+1; kend_hcm = nlevs+1 - nlevs = nlevs + 1 - endif - - if (.not.allocated(hailcast_dhail1)) then - allocate(hailcast_dhail1(isco:ieco,jsco:jeco), & - hailcast_dhail2(isco:ieco,jsco:jeco), & - hailcast_dhail3(isco:ieco,jsco:jeco), & - hailcast_dhail4(isco:ieco,jsco:jeco), & - hailcast_dhail5(isco:ieco,jsco:jeco), & - hailcast_diam_mean(isco:ieco,jsco:jeco), & - hailcast_diam_std(isco:ieco,jsco:jeco)) - allocate ( hailcast_dhail1_max(isco:ieco,jsco:jeco) ) - allocate ( hailcast_dhail2_max(isco:ieco,jsco:jeco) ) - allocate ( hailcast_dhail3_max(isco:ieco,jsco:jeco) ) - allocate ( hailcast_dhail4_max(isco:ieco,jsco:jeco) ) - allocate ( hailcast_dhail5_max(isco:ieco,jsco:jeco) ) - - do i=isco,ieco - do j=jsco,jeco - hailcast_dhail1(i,j)=0; hailcast_dhail2(i,j)=0; hailcast_dhail3(i,j)=0; hailcast_dhail4(i,j)=0; hailcast_dhail5(i,j)=0 - hailcast_dhail1_max(i,j)=0; hailcast_dhail2_max(i,j)=0; hailcast_dhail3_max(i,j)=0; hailcast_dhail4_max(i,j)=0; hailcast_dhail5_max(i,j)=0 - enddo - enddo - endif - if (.not.allocated(hailcast_wdur)) then - allocate(hailcast_wdur(isdo:iedo,jsdo:jedo),hailcast_wup_mask(isdo:iedo,jsdo:jedo)) - do i=isdo,iedo - do j=jsdo,jedo - hailcast_wdur(i,j)=0;hailcast_wup_mask(i,j)=0 - enddo - enddo - endif - endif + !initialize hailcast + call hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus) ! !------------------------------------ @@ -578,7 +455,7 @@ subroutine fv_nggps_diag(Atm, zvir, Time) logical :: bad_range real :: ptop, allmax real, allocatable :: wk(:,:,:), wk2(:,:,:) - real, dimension(:,:),allocatable :: ustm,vstm,srh01,srh03,hailcast_dhail_max + real, dimension(:,:),allocatable :: ustm,vstm,srh01,srh03 n = 1 ngc = Atm(n)%ng @@ -865,59 +742,48 @@ subroutine fv_nggps_diag(Atm, zvir, Time) call store_data(id_uhmin25, uhmin25, Time, kstt_uhmin25, kend_uhmin25) endif - !--- max hourly hailcast hail diameter - if ( .not.Atm(n)%flagstruct%hydrostatic .and. id_hailcast_dhail > 0) then - allocate(hailcast_dhail_max(isco:ieco,jsco:jeco)) - - do j=jsco,jeco - do i=isco,ieco - hailcast_dhail_max(i,j) = 0. - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail1_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail2_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail3_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail4_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail5_max(i,j)) - end do - end do - - call store_data(id_hailcast_dhail, hailcast_dhail_max, Time, kstt_hc, kend_hc) - deallocate(hailcast_dhail_max) - endif - - !--- max hourly hailcast hail diameter (embryo 1) - if ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. id_hailcast_dhail5 > 0) then - call store_data(id_hailcast_dhail5, hailcast_dhail5_max, Time, kstt_hc5, kend_hc5) - endif - - !--- hailcast updraft duration - if ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 ( .not.Atm(n)%flagstruct%hydrostatic .and. 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 - + IF ( .not. Atm(n)%flagstruct%hydrostatic .and. do_hailcast ) THEN + !--- max hourly hailcast hail diameter + if (id_hailcast_dhail > 0) then + call hailcast_compute_dhailmax(isco,ieco,jsco,jeco,istatus) + call store_data(id_hailcast_dhail, hailcast_dhail_max, Time, kstt_hc, kend_hc) + if (allocated(hailcast_dhail_max)) deallocate(hailcast_dhail_max) + 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() @@ -934,7 +800,6 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) integer, save :: kdtt = 0 real :: avg_max_length real,dimension(:,:,:),allocatable :: vort - real, allocatable :: rhoair_layer(:,:,:), z(:), z_layer(:,:,:), p_layer(:,:,:),zsfc(:,:) n = 1 ngc = Atm(n)%ng nq = size (Atm(n)%q,4) @@ -1019,87 +884,9 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) endif !allocate hailcast met field arrays - if (id_hailcast_dhail>0 .or. id_hailcast_dhail1>0 .or. id_hailcast_dhail2>0 .or. id_hailcast_dhail3>0 & - .or. id_hailcast_dhail4>0 .or. id_hailcast_dhail5>0 .or. id_hailcast_diam_mean>0 & - .or. id_hailcast_diam_std>0) then - - do j=jsco,jeco - do i=isco,ieco - if(mod(kdtt,nsteps_per_reset)==0)then - hailcast_dhail1_max(i,j) = 0. - hailcast_dhail2_max(i,j) = 0. - hailcast_dhail3_max(i,j) = 0. - hailcast_dhail4_max(i,j) = 0. - hailcast_dhail5_max(i,j) = 0. - endif - - hailcast_dhail1(i,j) = 0 - hailcast_dhail2(i,j) = 0 - hailcast_dhail3(i,j) = 0 - hailcast_dhail4(i,j) = 0 - hailcast_dhail5(i,j) = 0 - enddo - enddo - - - if (.not. allocated(rhoair_layer)) allocate(rhoair_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D air density - if (.not. allocated(z)) allocate(z(1:npzo+1)) !interace levels - if (.not. allocated(z_layer)) allocate(z_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D height ASL - if (.not. allocated(p_layer)) allocate(p_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D pressure - if (.not. allocated(zsfc)) allocate(zsfc(isco:ieco,jsco:jeco)) !terrain height - - do k=npzo,1,-1 - do j=jsco, jeco - if (Atm(n)%flagstruct%hydrostatic) then - call mpp_error(FATAL, 'HAILCAST can only be run with non-hydrostatic FV3') - !do i=is, ie - !rhoair(i,j,k) = Atm(n)%delp(i,j,k)/( ( Atm(n)%peln(i,k+1,j)- Atm(n)%peln(i,k,j)) & - ! * rdgas * Atm(n)%pt(i,j,k) * ( 1. + zvir*q(i,j,k,sphum))) - !enddo - else !non-hydrostatic - do i=isco, ieco - if (k==npzo) then - zsfc(i,j)=Atm(n)%phis(i,j) / grav - z(npzo+1)=zsfc(i,j) - endif - rhoair_layer(i,j,k) = -Atm(n)%delp(i,j,k)/(grav*Atm(n)%delz(i,j,k)) - !height of interfaces and layer height - z(k)=z(k+1)-Atm(n)%delz(i,j,k) - z_layer(i,j,k)=(z(k+1)+z(k))/2 - p_layer(i,j,k)=Atm(n)%delp(i,j,k)*(1.-sum(Atm(n)%q(i,j,k,2:Atm(n)%flagstruct%nwat)))/& - (-Atm(n)%delz(i,j,k)*grav)*rdgas*Atm(n)%pt(i,j,k)*(1.+zvir*Atm(n)%q(i,j,k,sphum)) - enddo - endif - enddo !j loop - enddo !k loop - - !call hailcast diagnostic driver once per time step and provide three-dimensional met fields - call hailcast_diagnostic_driver(Atm(n)%pt(isco:ieco,jsco:jeco,1:npzo), Atm(n)%w(isco:ieco,jsco:jeco,1:npzo), p_layer, z_layer, rhoair_layer, Atm(n)%q(isco:ieco,jsco:jeco,1:npzo,:), & !3D fields - Atm(n)%flagstruct%nwat, sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !number of tracer indices, indices - isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo, & !grid dimensions data array (with halo) and physical grid dimensions - 1,npzo, zsfc, & !vertical dimensions, terrain height - first_time, Atm(n)%domain, & !call hailcast every model step and info for updating haloes - hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5, & !hailcast embryos sizes - hailcast_diam_mean, hailcast_diam_std, & !hailcast embryo mean/std - hailcast_wdur, hailcast_wup_mask) !persistent updraft duration and mask - - do j=jsco,jeco - do i=isco,ieco - hailcast_dhail1_max(i,j) = max(hailcast_dhail1_max(i,j), hailcast_dhail1(i,j)) - hailcast_dhail2_max(i,j) = max(hailcast_dhail2_max(i,j), hailcast_dhail2(i,j)) - hailcast_dhail3_max(i,j) = max(hailcast_dhail3_max(i,j), hailcast_dhail3(i,j)) - hailcast_dhail4_max(i,j) = max(hailcast_dhail4_max(i,j), hailcast_dhail4(i,j)) - hailcast_dhail5_max(i,j) = max(hailcast_dhail5_max(i,j), hailcast_dhail5(i,j)) - enddo - enddo - - - !deallocate hailcast met variables - deallocate(p_layer) - deallocate(z_layer) - deallocate(rhoair_layer) - deallocate(z) - deallocate(zsfc) + if (do_hailcast) then + call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & + isco,jsco,ieco,jeco,npzo,kdtt,nsteps_per_reset) endif kdtt=kdtt+1 @@ -1608,69 +1395,71 @@ 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. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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( .not.hydrostatico .and. 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 + 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) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 34e5a2fca..37a64038f 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -62,8 +62,8 @@ module fv_arrays_mod integer :: steps !for hailcast - integer :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5 - integer :: id_hailcast_diam_mean, id_hailcast_diam_std + !integer :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5 + !integer :: id_hailcast_diam_mean, id_hailcast_diam_std end type fv_diag_type diff --git a/tools/module_diag_hailcast.F90 b/tools/module_diag_hailcast.F90 index 8aae38314..bf9177e79 100644 --- a/tools/module_diag_hailcast.F90 +++ b/tools/module_diag_hailcast.F90 @@ -2,77 +2,316 @@ !John Henderson AER 20190425 !time management handling from WRF is preserved but commented out -!MODULE module_diag_hailcast -!CONTAINS -! SUBROUTINE diag_hailcast_stub -! END SUBROUTINE diag_hailcast_stub -!END MODULE module_diag_hailcast +! Y. Wang NSSL/CIWRO 20220310 +! Remodeled as recommended for PR #164. +! MODULE module_diag_hailcast - + USE time_manager_mod, ONLY: time_type + use mpp_mod, only: mpp_error, FATAL, input_nml_file, stdlog, mpp_pe, mpp_root_pe + use fms_mod, only: check_nml_error + !use fv_mp_mod, only: is_master + + use constants_mod, only: grav + + logical :: do_hailcast = .false. !This controls whether hailcast is used + + INTEGER :: id_hailcast_dhail + INTEGER :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, & + id_hailcast_dhail4, id_hailcast_dhail5, id_hailcast_wdur, & + id_hailcast_wup_mask + INTEGER :: id_hailcast_diam_mean, id_hailcast_diam_std + REAL, ALLOCATABLE :: hailcast_dhail1(:,:), hailcast_dhail2(:,:), & + hailcast_dhail3(:,:), hailcast_dhail4(:,:), & + hailcast_dhail5(:,:) !hailstone diameters (mm) + REAL, ALLOCATABLE :: hailcast_dhail1_max(:,:), hailcast_dhail2_max(:,:), & + hailcast_dhail3_max(:,:), hailcast_dhail4_max(:,:), & + hailcast_dhail5_max(:,:) !hailstone diameters (mm) + REAL, ALLOCATABLE :: hailcast_diam_mean(:,:), hailcast_diam_std(:,:) !mean and standard deviation of five hailstones (mm) + REAL, ALLOCATABLE :: hailcast_wdur(:,:), hailcast_wup_mask(:,:) !persistent arrays for updraft duration (s) and mask + real, allocatable :: hailcast_dhail_max(:,:) + + integer :: kstt_hc, kend_hc + integer :: kstt_hc1,kstt_hc2,kstt_hc3,kstt_hc4,kstt_hc5 + integer :: kend_hc1,kend_hc2,kend_hc3,kend_hc4,kend_hc5 + integer :: kstt_hcd, kend_hcd, kstt_hcm, kend_hcm + + PRIVATE :: INTERPP, INTERP, TERMINL, VAPORCLOSE, MASSAGR, HEATBUD, BREAKUP, MELT + PRIVATE :: hailstone_driver, hailcast_diagnostic_driver CONTAINS -! SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags & -! , moist & -! , rho & -! , ids, ide, jds, jde, kds, kde & -! , ims, ime, jms, jme, kms, kme & -! , ips, ipe, jps, jpe, kps, kpe & -! , its, ite, jts, jte & -! , k_start, k_end & -! , dt, itimestep & -! , haildt,curr_secs & -! , haildtacttime ) - SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3D temperature, updraft, pressure, height, density and moisture tracers + SUBROUTINE hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus) + + use diag_manager_mod, only: register_diag_field + + CHARACTER(LEN=64), INTENT(IN) :: file_name + integer, intent(in) :: axes(4) + type(time_type), intent(in) :: Time + INTEGER, INTENT(IN) :: isdo,iedo,jsdo,jedo + INTEGER, INTENT(INOUT) :: nlevs + real, INTENT(IN) :: missing_value + INTEGER, INTENT(OUT) :: istatus + + namelist /fv_diagnostics_nml/ do_hailcast + integer :: ios, ierr + integer :: unit, f_unit + + !namelist file for hailcast +#ifdef INTERNAL_FILE_NML + ! Read Main namelist + read (input_nml_file,fv_diagnostics_nml,iostat=ios) + ierr = check_nml_error(ios,'fv_diagnostics_nml') +#else + f_unit=open_namelist_file() + rewind (f_unit) + ! Read Main namelist + read (f_unit,fv_diagnostics_nml,iostat=ios) + ierr = check_nml_error(ios,'fv_diagnostics_nml') + call close_file(f_unit) +#endif + + unit = stdlog() + write(unit, nml=fv_diagnostics_nml) + !end hailcast nml + + if (mpp_pe() == mpp_root_pe()) then + print*, 'do_hailcast = ', do_hailcast + end if + + !register hailcast arrays + if (do_hailcast) then + id_hailcast_dhail = register_diag_field ( trim(file_name), 'hailcast_dhail_max', axes(1:2), Time, & + 'hourly max hail diameter', 'mm', missing_value=missing_value ) + id_hailcast_dhail1 = register_diag_field ( trim(file_name), 'hailcast_dhail1_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 1)', 'mm', missing_value=missing_value ) + id_hailcast_dhail2 = register_diag_field ( trim(file_name), 'hailcast_dhail2_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 2)', 'mm', missing_value=missing_value ) + id_hailcast_dhail3 = register_diag_field ( trim(file_name), 'hailcast_dhail3_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 3)', 'mm', missing_value=missing_value ) + id_hailcast_dhail4 = register_diag_field ( trim(file_name), 'hailcast_dhail4_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 4)', 'mm', missing_value=missing_value ) + id_hailcast_dhail5 = register_diag_field ( trim(file_name), 'hailcast_dhail5_max', axes(1:2), Time, & + 'hourly max hail diameter (embryo 5)', 'mm', missing_value=missing_value ) + id_hailcast_diam_mean = register_diag_field ( trim(file_name), 'hailcast_diam_mean', axes(1:2), Time, & + 'mean hail diameter', 'mm', missing_value=missing_value ) + id_hailcast_diam_std = register_diag_field ( trim(file_name), 'hailcast_diam_std', axes(1:2), Time, & + 'standard deviation of hail diameter', 'mm', missing_value=missing_value ) + id_hailcast_wdur = register_diag_field ( trim(file_name), 'hailcast_wdur', axes(1:2), Time, & + 'updraft duration', 's', missing_value=missing_value ) + id_hailcast_wup_mask = register_diag_field ( trim(file_name), 'hailcast_wup_mask', axes(1:2), Time, & + 'updraft mask', '', missing_value=missing_value ) + + if (id_hailcast_dhail > 0) then + kstt_hc = nlevs+1; kend_hc = nlevs+1 + nlevs = nlevs + 1 + endif + + if (id_hailcast_dhail1 > 0) then + kstt_hc1 = nlevs+1; kend_hc1 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail2 > 0) then + kstt_hc2 = nlevs+1; kend_hc2 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail3 > 0) then + kstt_hc3 = nlevs+1; kend_hc3 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail4 > 0) then + kstt_hc4 = nlevs+1; kend_hc4 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_dhail5 > 0) then + kstt_hc5 = nlevs+1; kend_hc5 = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_wdur > 0) then + kstt_hcd = nlevs+1; kend_hcd = nlevs+1 + nlevs = nlevs + 1 + endif + if (id_hailcast_wup_mask > 0) then + kstt_hcm = nlevs+1; kend_hcm = nlevs+1 + nlevs = nlevs + 1 + endif + + if (.not.allocated(hailcast_dhail1)) then + allocate(hailcast_dhail1(isco:ieco,jsco:jeco), & + hailcast_dhail2(isco:ieco,jsco:jeco), & + hailcast_dhail3(isco:ieco,jsco:jeco), & + hailcast_dhail4(isco:ieco,jsco:jeco), & + hailcast_dhail5(isco:ieco,jsco:jeco), & + hailcast_diam_mean(isco:ieco,jsco:jeco), & + hailcast_diam_std(isco:ieco,jsco:jeco)) + allocate ( hailcast_dhail1_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail2_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail3_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail4_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail5_max(isco:ieco,jsco:jeco) ) + + do i=isco,ieco + do j=jsco,jeco + hailcast_dhail1(i,j)=0; hailcast_dhail2(i,j)=0 + hailcast_dhail3(i,j)=0; hailcast_dhail4(i,j)=0 + hailcast_dhail5(i,j)=0 + hailcast_dhail1_max(i,j)=0; hailcast_dhail2_max(i,j)=0 + hailcast_dhail3_max(i,j)=0; hailcast_dhail4_max(i,j)=0 + hailcast_dhail5_max(i,j)=0 + enddo + enddo + endif + if (.not.allocated(hailcast_wdur)) then + allocate(hailcast_wdur(isdo:iedo,jsdo:jedo),hailcast_wup_mask(isdo:iedo,jsdo:jedo)) + do i=isdo,iedo + do j=jsdo,jedo + hailcast_wdur(i,j)=0 + hailcast_wup_mask(i,j)=0 + enddo + enddo + endif + endif + + RETURN + END SUBROUTINE hailcast_init + + SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & + isco,jsco,ieco,jeco,npzo, kdtt, nsteps_per_reset) + + !moved this code from subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) in fv_nggps_diag + use fv_arrays_mod, only: fv_atmos_type + type(fv_atmos_type), intent(inout) :: Atm + integer, INTENT(IN) :: isco, ieco, jsco, jeco, npzo + integer, INTENT(IN) :: sphum, liq_wat, ice_wat !< GFDL physics + integer, INTENT(IN) :: rainwat, snowwat, graupel + integer, INTENT(IN) :: nsteps_per_reset, kdtt + + !declare temporary hailcast variables + real, allocatable :: rhoair_layer(:,:,:), z(:), z_layer(:,:,:), p_layer(:,:,:),zsfc(:,:) + !automatic 3-D arrays for layer air density, height ASL and pressure + !automatic 2-D array for terrain height + !automatic 1-D array for interface levels + + IF (id_hailcast_dhail>0 .or. id_hailcast_dhail1>0 .or. & + id_hailcast_dhail2>0 .or. id_hailcast_dhail3>0 .or. & + id_hailcast_dhail4>0 .or. id_hailcast_dhail5>0 .or. & + id_hailcast_diam_mean>0 .or. id_hailcast_diam_std>0) THEN + do j=jsco,jeco + do i=isco,ieco + if(mod(kdtt,nsteps_per_reset)==0)then + hailcast_dhail1_max(i,j) = 0. + hailcast_dhail2_max(i,j) = 0. + hailcast_dhail3_max(i,j) = 0. + hailcast_dhail4_max(i,j) = 0. + hailcast_dhail5_max(i,j) = 0. + endif + + hailcast_dhail1(i,j) = 0 + hailcast_dhail2(i,j) = 0 + hailcast_dhail3(i,j) = 0 + hailcast_dhail4(i,j) = 0 + hailcast_dhail5(i,j) = 0 + enddo + enddo + + + if (.not. allocated(rhoair_layer)) allocate(rhoair_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D air density + if (.not. allocated(z)) allocate(z(1:npzo+1)) !interace levels + if (.not. allocated(z_layer)) allocate(z_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D height ASL + if (.not. allocated(p_layer)) allocate(p_layer(isco:ieco,jsco:jeco,1:npzo)) !layer 3-D pressure + if (.not. allocated(zsfc)) allocate(zsfc(isco:ieco,jsco:jeco)) !terrain height + + do k=npzo,1,-1 + do j=jsco, jeco + if (Atm%flagstruct%hydrostatic) then + call mpp_error(FATAL, 'HAILCAST can only be run with non-hydrostatic FV3') + !do i=is, ie + !rhoair(i,j,k) = Atm%delp(i,j,k)/( ( Atm%peln(i,k+1,j)- Atm%peln(i,k,j)) & + ! * rdgas * Atm%pt(i,j,k) * ( 1. + zvir*q(i,j,k,sphum))) + !enddo + else !non-hydrostatic + do i=isco, ieco + if (k==npzo) then + zsfc(i,j)=Atm%phis(i,j) / grav + z(npzo+1)=zsfc(i,j) + endif + rhoair_layer(i,j,k) = -Atm%delp(i,j,k)/(grav*Atm%delz(i,j,k)) + !height of interfaces and layer height + z(k)=z(k+1)-Atm%delz(i,j,k) + z_layer(i,j,k)=(z(k+1)+z(k))/2 + p_layer(i,j,k)=Atm%delp(i,j,k)*(1.-sum(Atm%q(i,j,k,2:Atm%flagstruct%nwat)))/& + (-Atm%delz(i,j,k)*grav)*rdgas*Atm%pt(i,j,k)*(1.+zvir*Atm%q(i,j,k,sphum)) + enddo + endif + enddo !j loop + enddo !k loop + + !call hailcast diagnostic driver once per time step and provide three-dimensional met fields + call hailcast_diagnostic_driver(Atm%pt(isco:ieco,jsco:jeco,1:npzo), Atm%w(isco:ieco,jsco:jeco,1:npzo), p_layer, z_layer, rhoair_layer, Atm%q(isco:ieco,jsco:jeco,1:npzo,:), & !3D fields + Atm%flagstruct%nwat, sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !number of tracer indices, indices + isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo, & !grid dimensions data array (with halo) and physical grid dimensions + 1,npzo, zsfc, & !vertical dimensions, terrain height + first_time, Atm%domain) !call hailcast every model step and info for updating haloes + !hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5, & !hailcast embryos sizes + !hailcast_diam_mean, hailcast_diam_std, & !hailcast embryo mean/std + !hailcast_wdur, hailcast_wup_mask) !persistent updraft duration and mask + + do j=jsco,jeco + do i=isco,ieco + hailcast_dhail1_max(i,j) = max(hailcast_dhail1_max(i,j), hailcast_dhail1(i,j)) + hailcast_dhail2_max(i,j) = max(hailcast_dhail2_max(i,j), hailcast_dhail2(i,j)) + hailcast_dhail3_max(i,j) = max(hailcast_dhail3_max(i,j), hailcast_dhail3(i,j)) + hailcast_dhail4_max(i,j) = max(hailcast_dhail4_max(i,j), hailcast_dhail4(i,j)) + hailcast_dhail5_max(i,j) = max(hailcast_dhail5_max(i,j), hailcast_dhail5(i,j)) + enddo + enddo + + + !deallocate hailcast met variables + deallocate(p_layer) + deallocate(z_layer) + deallocate(rhoair_layer) + deallocate(z) + deallocate(zsfc) + END IF + + RETURN + END SUBROUTINE hailcast_compute + + SUBROUTINE hailcast_compute_dhailmax(isco,ieco,jsco,jeco,istatus) + + INTEGER, INTENT(IN) :: isco,ieco,jsco,jeco + INTEGER, INTENT(OUT) :: istatus + + allocate(hailcast_dhail_max(isco:ieco,jsco:jeco)) + + do j=jsco,jeco + do i=isco,ieco + hailcast_dhail_max(i,j) = 0. + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail1_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail2_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail3_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail4_max(i,j)) + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail5_max(i,j)) + end do + end do + + RETURN + END SUBROUTINE hailcast_compute_dhailmax + + SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3D temperature, updraft, pressure, height, density and moisture tracers moist_count,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & ! tracer indices - is,ie,js,je,ishalo,iehalo,jshalo,jehalo, & ! i,j halo data array dimensions and physical grid dimensions (2 smaller on each end) - k_start, k_end, terr_z, & ! bottom and top z indices and terrain height - dt, domainid, & ! call frequency of hailcast = model time step = physics - hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5, hailcast_diam_mean, hailcast_diam_std, hailcast_wdur, hailcast_wup_mask) !hailcast persistent variables - -! USE module_domain, ONLY : domain , domain_clock_get -! USE module_configure, ONLY : grid_config_rec_type, model_config_rec -! USE module_state_description -! USE module_model_constants -! USE module_utility -! USE module_streams, ONLY: history_alarm, auxhist2_alarm -!#ifdef DM_PARALLEL -! USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval -! #endif - - use constants_mod, only: grav, rdgas, rvgas, pi=>pi_8, radius, kappa, WTMAIR, WTMCO2, & - omega, hlv, cp_air, cp_vapor - !use fms_mod, only: write_version_number - use fms_io_mod, only: set_domain, nullify_domain - use time_manager_mod, only: time_type, get_date, get_time - use mpp_domains_mod, only: domain2d, mpp_update_domains, DGRID_NE - use diag_manager_mod, only: diag_axis_init, register_diag_field, & - register_static_field, send_data, diag_grid_init - use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_diag_type, fv_grid_bounds_type, & - R_GRID - use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index - use field_manager_mod, only: MODEL_ATMOS - use fv_mp_mod, only: mp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master - use mpp_mod, only: mpp_error, FATAL, stdlog, mpp_pe, mpp_root_pe, mpp_sum, mpp_max, NOTE - use sat_vapor_pres_mod, only: compute_qs, lookup_es - - use fv_arrays_mod, only: max_step + is,ie,js,je,ishalo,iehalo,jshalo,jehalo, & ! i,j halo data array dimensions and physical grid dimensions (2 smaller on each end) + k_start, k_end, terr_z, & ! bottom and top z indices and terrain height + dt, domainid) ! call frequency of hailcast = model time step = physics + use mpp_domains_mod, only: domain2d, mpp_update_domains IMPLICIT NONE -! TYPE ( domain ), INTENT(INOUT) :: grid -! TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags -! type(fv_atmos_type), intent(inout) :: grid -!!!! type(time_type), intent(in) :: Time !will not compile - INTEGER:: ishalo, iehalo, jshalo, jehalo, k_start, k_end !dimensions of halo data array - INTEGER:: is, ie, js, je - -! INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, & -! ims, ime, jms, jme, kms, kme, & -! ips, ipe, jps, jpe, kps, kpe - - logical :: master + INTEGER:: ishalo, iehalo, jshalo, jehalo, k_start, k_end !dimensions of halo data array + INTEGER:: is, ie, js, je + + !logical :: master INTEGER :: moist_count,sphum,rainwat,snowwat,graupel,liq_wat,ice_wat REAL, DIMENSION( is:ie, js:je, k_start: k_end, moist_count), & @@ -84,15 +323,6 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 REAL, DIMENSION( is:ie, js:je), & INTENT(IN ) :: terr_z - REAL, DIMENSION( is:ie, js:je), & - INTENT(INOUT) :: hailcast_dhail1, hailcast_dhail2, & - hailcast_dhail3, hailcast_dhail4, & - hailcast_dhail5, & - hailcast_diam_mean, hailcast_diam_std - - REAL, DIMENSION( ishalo:iehalo, jshalo:jehalo), & - INTENT(INOUT) :: hailcast_wdur, hailcast_wup_mask - ! REAL, INTENT(IN ),OPTIONAL :: haildt ! REAL, INTENT(IN ),OPTIONAL :: curr_secs ! REAL, INTENT(INOUT),OPTIONAL :: haildtacttime @@ -104,7 +334,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! Local ! ----- CHARACTER*512 :: message - CHARACTER*256 :: timestr + CHARACTER*256 :: timestr INTEGER :: i,j,k, numlevs REAL, DIMENSION( is:ie, js:je, k_start:k_end ) :: qr & , qs & @@ -112,12 +342,12 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 , qv & , qc & , qi & - , ptot + , ptot REAL, DIMENSION( ishalo:iehalo, jshalo:jehalo) :: wdur_prev, wup_mask_prev REAL :: dhail1,dhail2,dhail3,dhail4,dhail5 - + ! ! Timing ! ! ------ ! TYPE(WRFU_Time) :: hist_time, aux2_time, CurrTime, StartTime @@ -132,7 +362,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! write ( message, * ) 'inside hailcast_diagnostics_driver' ! CALL wrf_debug( 100 , message ) -! ! Get timing info +! ! Get timing info ! ! Want to know if when the last history output was ! ! Check history and auxhist2 alarms to check last ring time and how often ! ! they are set to ring @@ -145,7 +375,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! ! Get domain clock ! ! ---------------- ! CALL domain_clock_get ( grid, current_time=CurrTime, & -! simulationStartTime=StartTime, & +! simulationStartTime=StartTime, & ! current_timestr=timestr, time_step=dtint ) ! ! Set some booleans for use later @@ -163,8 +393,8 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! ! Following uses an overloaded .eq. ! ! --------------------------------- ! is_first_timestep = ( Currtime .eq. StartTime + dtint ) - - master = (mpp_pe()==mpp_root_pe()) .or. is_master() + + !master = (mpp_pe()==mpp_root_pe()) .or. is_master() !FV3: No need to reset arrays to zero at dump time, thus no need to keep track of dump times ! @@ -183,7 +413,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 !ENDIF ! is_after_history_dump ! We need to do some neighboring gridpoint comparisons for the updraft - ! duration calculations; set i,j start and end values so we don't go off + ! duration calculations; set i,j start and end values so we don't go off ! the edges of the domain. Updraft duration on domain edges will always be 0. ! ---------------------------------------------------------------------- @@ -198,7 +428,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! config_flags%nested) j_end = MIN( jde-2, jte ) !-1 ! IF ( config_flags%periodic_x ) i_start = its ! IF ( config_flags%periodic_x ) i_end = ite - + ! Make a copy of the updraft duration, mask variables ! --------------------------------------------------- @@ -225,7 +455,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 hailcast_wdur(i,j) = 0 DO k = k_start, k_end - IF ( w(i,j,k) .ge. 10. ) THEN + IF ( w(i,j,k) .ge. 10. ) THEN hailcast_wup_mask(i,j) = 1 ENDIF ENDDO @@ -264,13 +494,13 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! END IF ! END IF -! ! Calculate STEPHAIL -! stephail = NINT(haildt / dt) +! ! Calculate STEPHAIL +! stephail = NINT(haildt / dt) ! stephail = MAX(stephail,1) ! ! itimestep starts at 1 - we want it to start at 0 so correctly ! ! divisibile by stephail. -! itimestep_basezero = itimestep -1 +! itimestep_basezero = itimestep -1 ! ! Do we run through this scheme or not? @@ -286,9 +516,9 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! ! CURR_SECS >= HAILDTACTTIME ! ! If we do run through the scheme, we set the flag run_param to TRUE and we set -! ! the decided flag to TRUE. The decided flag says that one of these tests was -! ! able to say "yes", run the scheme. We only proceed to other tests if the -! ! previous tests all have left decided as FALSE. If we set run_param to TRUE +! ! the decided flag to TRUE. The decided flag says that one of these tests was +! ! able to say "yes", run the scheme. We only proceed to other tests if the +! ! previous tests all have left decided as FALSE. If we set run_param to TRUE ! ! and this is adaptive time stepping, we set the time to the next hailcast run. ! run_param = .FALSE. @@ -332,7 +562,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ! CALL wrf_debug( 100 , message ) ! IF (run_param) THEN - + ! 3-D arrays for moisture variables ! --------------------------------- numlevs=k_end-k_start+1 @@ -340,16 +570,16 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 DO k=1, numlevs DO j=js, je qv(i,j,k)=0;qr(i,j,k)=0;qs(i,j,k)=0;qg(i,j,k)=0;qc(i,j,k)=0;qi(i,j,k)=0 - qv(i,j,k) = q(i,j,k,sphum) - qr(i,j,k) = q(i,j,k,rainwat) - qs(i,j,k) = q(i,j,k,snowwat) - qg(i,j,k) = q(i,j,k,graupel) - qc(i,j,k) = q(i,j,k,liq_wat) - qi(i,j,k) = q(i,j,k,ice_wat) + qv(i,j,k) = q(i,j,k,sphum) + qr(i,j,k) = q(i,j,k,rainwat) + qs(i,j,k) = q(i,j,k,snowwat) + qg(i,j,k) = q(i,j,k,graupel) + qc(i,j,k) = q(i,j,k,liq_wat) + qi(i,j,k) = q(i,j,k,ice_wat) ENDDO ENDDO ENDDO - + ! Hail diameter in millimeters (HAILCAST) ! --------------------------------------------------- DO j = js, je @@ -368,7 +598,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 qi(i,j,k_end:k_start:-1), & qc(i,j,k_end:k_start:-1), & qr(i,j,k_end:k_start:-1), & - qs(i,j,k_end:k_start:-1), & + qs(i,j,k_end:k_start:-1), & qg(i,j,k_end:k_start:-1), & w(i,j,k_end:k_start:-1), & hailcast_wdur(i,j), & @@ -424,7 +654,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3 ENDDO !END IF -END SUBROUTINE hailcast_diagnostic_driver + END SUBROUTINE hailcast_diagnostic_driver !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! @@ -444,13 +674,13 @@ END SUBROUTINE hailcast_diagnostic_driver !!!! VUU updraft speed at each level (m/s) !!!! Float !!!! ht terrain height (m) -!!!! wdur duration of any updraft > 10 m/s within 1 surrounding -!!!! gridpoint +!!!! wdur duration of any updraft > 10 m/s within 1 surrounding +!!!! gridpoint !!!! Integer !!!! nz number of vertical levels !!!! !!!! Output: -!!!! dhail hail diameter in mm +!!!! dhail hail diameter in mm !!!! 1st-5th rank-ordered hail diameters returned !!!! !!!! 13 Aug 2013 .................................Becky Adams-Selin AER @@ -468,8 +698,8 @@ END SUBROUTINE hailcast_diagnostic_driver !!!! embryo size and location to better match literature. Added !!!! enhanced melting when hailstone collides with liquid water !!!! in regions above freezing. Final diameter returned is ice diameter -!!!! only. Added hailstone temperature averaging over previous timesteps -!!!! to decrease initial temperature instability at small embyro diameters. +!!!! only. Added hailstone temperature averaging over previous timesteps +!!!! to decrease initial temperature instability at small embyro diameters. !!!! 3 Sep 2015...................................Becky Adams-Selin AER !!!! Insert embryos at -13C; interpret pressure and other variables to !!!! that exact temperature level. @@ -481,16 +711,16 @@ END SUBROUTINE hailcast_diagnostic_driver !!!! the saturation vapor pressure using the model temperature !!!! profile. !!!! 04 Feb 2017...................................Becky Adams-Selin AER -!!!! Added adaptive time-stepping option. +!!!! Added adaptive time-stepping option. !!!! ***Don't set any higher than 60 seconds*** !!!! 06 May 2017...................................Becky Adams-Selin AER !!!! Bug fixes for some memory issues in the adiabatic liquid !!!! water profile code. -!!!! +!!!! !!!! See Adams-Selin and Ziegler 2016, MWR for further documentation. !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& + SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& RA, qi1d,qc1d,qr1d,qs1d,qg1d, & VUU, wdur, & nz,dhail1,dhail2,dhail3,dhail4, & @@ -510,7 +740,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& REAL, INTENT(IN ) :: ht & , wdur - + !Output: 1st-5th rank-ordered hail diameters returned REAL, INTENT(INOUT) :: dhail1 & ! hail diameter (mm); , dhail2 & @@ -559,7 +789,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& RSA ! saturation mixing ratio !in-cloud updraft parameters at location of hailstone REAL P ! in-cloud pressure (Pa) - REAL RS ! in-cloud saturation mixing ratio + REAL RS ! in-cloud saturation mixing ratio REAL RI, RW ! ice, liquid water mix. ratio (kg/kg) REAL XI, XW ! ice, liquid water content (kg/m3 air) REAL PC ! in-cloud fraction of frozen water @@ -578,7 +808,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& !internal function variables REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE REAL dum, icefactor, adiabatic_factor - + REAL sec, secdel ! time step, increment in seconds INTEGER i, j, k, IFOUT, ind(1) CHARACTER*256 :: message @@ -586,7 +816,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& secdel = 5.0 g=9.81 r_d = 287. - + !!!!!!!!!!!!!!!! 1. UPDRAFT PROPERTIES !!!!!!!!!!!!!!!!!!!!!!! !!! DEFINE UPDRAFT PROPERTIES !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -599,7 +829,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& VUU_pert(i) = VUU(i) * 1. ENDDO - + ! Initialize diameters to 0. DO i=1,5 dhails(i) = 0. @@ -612,7 +842,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& RTIME = wdur ENDIF - + !!!!!!!!!!!!!!!! 2. INITIAL EMBRYO !!!!!!!!!!!!!!!!!!!!!!!!!!! !!! FIND PROPERTIES OF INITIAL EMBRYO LOCATION !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -646,7 +876,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& RBAS = RA(KBAS) !At coarser resolutions WRF underestimates liquid water aloft. - !A fairer estimate is of the adiabatic liquid water profile, or + !A fairer estimate is of the adiabatic liquid water profile, or !the difference between the saturation mixing ratio at each level !and the cloud base water vapor mixing ratio DO k=1,nz @@ -683,7 +913,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1) IF (RWA_new(k).LT.0) RWA_new(k) = 0. ENDIF - ! - old code - ! + ! - old code - ! !IF (k.eq.KBAS+1) THEN ! RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) !ELSE IF ((k.ge.KBAS+2).AND.(RWA_adiabat(k).ge.1.E-12)) THEN @@ -706,7 +936,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& ! these sizes were picked. !Run each initial embryo size perturbation DO i=1,5 - SELECT CASE (i) + SELECT CASE (i) CASE (1) !Initial hail embryo diameter in m, at cloud base DD = 5.E-3 @@ -714,7 +944,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& CASE (2) DD = 7.5E-3 tk_embryo = 265.155 !-8C - CASE (3) + CASE (3) DD = 5.E-3 tk_embryo = 260.155 !-13C CASE (4) @@ -738,7 +968,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& !Begin hail simulation time (seconds) sec = 0. - + !!!!!!!!!!!!!!!!!! 4. INITIAL PARAMETERS !!!!!!!!!!!!!!!!! !!! PUT INITIAL PARAMETERS IN VARIABLES !!! @@ -747,7 +977,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& P = WFZLP RS = RFZL TC = TFZL - VU = VUFZL + VU = VUFZL Z = ZFZL - ht LDEPTH = Z DENSA = DENSAFZL @@ -756,43 +986,43 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& nofroze=1 !Set test for embryo: 0 for never been frozen; 1 frozen TS = TC TSm1 = TS - TSm2 = TS + TSm2 = TS D = DD !hailstone diameter in m FW = 0.0 - DENSE = 500. !kg/m3 + DENSE = 500. !kg/m3 ITYPE=1. !Assume starts in dry growth. CLOUDON=1 !we'll eventually turn cloud "off" once updraft past time limit !Start time loop. DO WHILE (sec .lt. TAU) sec = sec + secdel - + !!!!!!!!!!!!!!!!!! 5. CALCULATE PARAMETERS !!!!!!!!!!!!!!!!! !!! CALCULATE UPDRAFT PROPERTIES !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !Intepolate vertical velocity to our new pressure level CALL INTERP(VUU_pert,VUMAX,P,IFOUT,PA,nz) !print *, 'INTERP VU: ', VU, P - + !Outside pressure levels? If so, exit loop IF (IFOUT.EQ.1) GOTO 100 - + !Sine wave multiplier on updraft strength IF (SEC .GT. 0.0 .AND. SEC .LT. RTIME) THEN VUCORE = VUMAX * (0.5*SIN((3.14159*SEC)/(RTIME))+0.5)*1.2 - VU = VUCORE + VU = VUCORE ELSEIF (SEC .GE. RTIME) THEN VU = 0.0 CLOUDON = 0 ENDIF - - !Calculate terminal velocity of the hailstone + + !Calculate terminal velocity of the hailstone ! (use previous values) CALL TERMINL(DENSA,DENSE,D,VT,TC) - + !Actual velocity of hailstone (upwards positive) V = VU - VT - + !Use hydrostatic eq'n to calc height of next level P = P - DENSA*g*V*secdel Z = Z + V*secdel @@ -800,10 +1030,10 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& !Interpolate cloud temp, qvapor at new p-level CALL INTERP(TCA,TC,P,IFOUT,PA,nz) CALL INTERP(RA,RS,P,IFOUT,PA,nz) - + !New density of in-cloud air DENSA=P/(r_d*(1.+0.609*RS/(1.+RS))*TC) - + !Interpolate liquid, frozen water mix ratio at new level CALL INTERP(RIA,RI,P,IFOUT,PA,nz) CALL INTERP(RWA_new,RW,P,IFOUT,PA,nz) @@ -814,32 +1044,32 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& ELSE PC = 1. ENDIF - - + + ! SATURATION VAPOUR DENSITY DIFFERENCE BETWTEEN STONE AND CLOUD CALL VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) - - + + !!!!!!!!!!!!!!!!!! 6. STONE'S MASS GROWTH !!!!!!!!!!!!!!!!!!!! CALL MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& - TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,secdel,ITYPE,DELRW) + TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,secdel,ITYPE,DELRW) !!!!!!!!!!!!!!!!!! 7. HEAT BUDGET OF HAILSTONE !!!!!!!!!!!!!!! - CALL HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & + CALL HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,secdel,ITYPE,P) - + !!!!! 8. TEST DIAMETER OF STONE AND HEIGHT ABOVE GROUND !!!!!!! - !!! TEST IF DIAMETER OF STONE IS GREATER THAN CRITICAL LIMIT, IF SO - !!! BREAK UP + !!! TEST IF DIAMETER OF STONE IS GREATER THAN CRITICAL LIMIT, IF SO + !!! BREAK UP WATER=FW*GM !KG - ! CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S SURFACE + ! CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S SURFACE CRIT = 2.0E-4 IF (WATER.GT.CRIT)THEN CALL BREAKUP(DENSE,D,GM,FW,CRIT) ENDIF - + !!! Has stone reached below cloud base? !IF (Z .LE. 0) GOTO 200 IF (Z .LE. ZBAS) then @@ -847,15 +1077,15 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& endif !calculate ice-only diameter size - D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 + D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 - !Has the stone entirely melted and it's below the freezing level? + !Has the stone entirely melted and it's below the freezing level? IF ((D_ICE .LT. 1.E-8) .AND. (TC.GT.273.155)) GOTO 300 !move values to previous timestep value TSm1 = TS TSm2 = TSm1 - + ENDDO !end cloud lifetime loop 100 CONTINUE !outside pressure levels in model @@ -863,7 +1093,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& 300 CONTINUE !stone has entirely melted and is below freezing level !!!!!!!!!!!!!!!!!! 9. MELT STONE BELOW CLOUD !!!!!!!!!!!!!!!!!!!! - !Did the stone shoot out the top of the storm? + !Did the stone shoot out the top of the storm? !Then let's assume it's lost in the murky "outside storm" world. IF (P.lt.PA(nz)) THEN D=0.0 @@ -873,7 +1103,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& ELSE IF (Z.GT.0) THEN !If still frozen, then use melt routine to melt below cloud ! based on mean below-cloud conditions. - + !Calculate mean sub-cloud layer conditions TSUM = 0. RSUM = 0. @@ -886,12 +1116,12 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& TLAYER = TSUM / KBAS PLAYER = PSUM / KBAS RLAYER = RSUM / KBAS - + !MELT is expecting a hailstone of only ice. At the surface !we're only interested in the actual ice diameter of the hailstone, !so let's shed any excess water now. - D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 - D = D_ICE + D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 + D = D_ICE CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) ENDIF !end check for melting call @@ -901,14 +1131,14 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& !Check to make sure something hasn't gone horribly wrong 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 ENDDO !end embryo size loop - - + + dhail1 = dhails(1) dhail2 = dhails(2) dhail3 = dhails(3) @@ -920,7 +1150,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& print*,'jmh KBAS,ZBAS,TBAS',ii,jj,KBAS,ZBAS,TBAS print*,'jmh WBASP,RBAS=',ii,jj,WBASP,RBAS print*,'jmh embryo sizes: ',ii,jj,dhail1,dhail2,dhail3,dhail4,dhail5 - if (ii.eq. 847 .and. jj.eq.394) then + if (ii.eq. 847 .and. jj.eq.394) then print*,'jmh i=847,j=394 TCA:',ii,jj,TCA(1:nz) print*,'jmh i=847,j=394 h1d:',ii,jj,h1d(1:nz) print*,'jmh i=847,j=394 ht:',ii,jj,ht @@ -936,7 +1166,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& print*,'jmh i=847,j=394 wdur:',ii,jj,wdur print*,'jmh i=847,j=394 nz:',ii,jj,nz endif - if (KBAS .gt. 30) then + if (KBAS .gt. 30) then print*,'jmh KBAS>30 TCA:',ii,jj,TCA(1:nz) print*,'jmh KBAS>30 h1d:',ii,jj,h1d(1:nz) print*,'jmh KBAS>30 ht:',ii,jj,ht @@ -951,13 +1181,12 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& print*,'jmh KBAS>30 VUU:',ii,jj,VUU(1:nz) print*,'jmh KBAS>30 wdur:',ii,jj,wdur print*,'jmh KBAS>30 nz:',ii,jj,nz - endif + endif print*,'jmh ------end of hailstone driver------- jmh' endif - END SUBROUTINE hailstone_driver - + END SUBROUTINE hailstone_driver - SUBROUTINE INTERPP(PA,PVAL,TA,TVAL,IFOUT,ITEL) + SUBROUTINE INTERPP(PA,PVAL,TA,TVAL,IFOUT,ITEL) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! !!!! INTERP: to linearly interpolate value of pval at temperature tval @@ -972,16 +1201,16 @@ SUBROUTINE INTERPP(PA,PVAL,TA,TVAL,IFOUT,ITEL) !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE - + REAL PVAL, TVAL REAL, DIMENSION( ITEL) :: TA, PA INTEGER ITEL, IFOUT !local variables INTEGER I REAL FRACT - + IFOUT=1 - + DO I=1,ITEL-1 IF ( (TVAL .LT. TA(I) .AND. TVAL .GE. TA(I+1)) .or. & ! dT/dz < 0 (TVAL .GT. TA(I) .AND. TVAL .LE. TA(I+1)) ) THEN ! dT/dz > 0 @@ -989,18 +1218,16 @@ SUBROUTINE INTERPP(PA,PVAL,TA,TVAL,IFOUT,ITEL) FRACT = (TA(I) - TVAL) / (TA(I) - TA(I+1)) !.... compute the pressure value pval at temperature tval PVAL = ((1.0 - FRACT) * PA(I)) + (FRACT * PA(I+1)) - + !End loop IFOUT=0 EXIT ENDIF ENDDO - - END SUBROUTINE INTERPP + END SUBROUTINE INTERPP - - SUBROUTINE INTERP(AA,A,P,IFOUT,PA,ITEL) + SUBROUTINE INTERP(AA,A,P,IFOUT,PA,ITEL) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! !!!! INTERP: to linearly interpolate values of A at level P @@ -1015,39 +1242,38 @@ SUBROUTINE INTERP(AA,A,P,IFOUT,PA,ITEL) !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE - + REAL A, P REAL, DIMENSION( ITEL) :: AA, PA INTEGER ITEL, IFOUT !local variables INTEGER I REAL PDIFF, VDIFF, RDIFF, VERH, ADIFF - + IFOUT=1 - + DO I=1,ITEL-1 IF (P.LE.PA(I) .AND. P.GT.PA(I+1)) THEN !Calculate ratio between vdiff and pdiff PDIFF = PA(I)-PA(I+1) VDIFF = PA(I)-P - VERH = VDIFF/PDIFF - + VERH = VDIFF/PDIFF + !Calculate the difference between the 2 A values RDIFF = AA(I+1) - AA(I) - + !Calculate new value A = AA(I) + RDIFF*VERH - + !End loop IFOUT=0 EXIT ENDIF ENDDO - - END SUBROUTINE INTERP - - SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) + END SUBROUTINE INTERP + + SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! !!!! INTERP: Calculate terminal velocity of the hailstone @@ -1060,28 +1286,28 @@ SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE - + REAL*8 D REAL DENSA, DENSE, TC, VT REAL GMASS, GX, RE, W, Y REAL, PARAMETER :: PI = 3.141592654, G = 9.78956 REAL ANU - + !Mass of stone in kg GMASS = (DENSE * PI * (D**3.)) / 6. - + !Dynamic viscosity ANU = (0.00001718)*(273.155+120.)/(TC+120.)*(TC/273.155)**(1.5) - - !CALC THE BEST NUMBER, X AND REYNOLDS NUMBER, RE + + !CALC THE BEST NUMBER, X AND REYNOLDS NUMBER, RE GX=(8.0*GMASS*G*DENSA)/(PI*(ANU*ANU)) RE=(GX/0.6)**0.5 - !SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON + !SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON !THE BEST NUMBER IF (GX.LT.550) THEN W=LOG10(GX) - Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0) + Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0) RE=10**Y VT=ANU*RE/(D*DENSA) ELSE IF (GX.GE.550.AND.GX.LT.1800) THEN @@ -1092,20 +1318,18 @@ SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) ELSE IF (GX.GE.1800.AND.GX.LT.3.45E08) THEN RE=0.4487*(GX**0.5536) VT=ANU*RE/(D*DENSA) - ELSE + ELSE RE=(GX/0.6)**0.5 VT=ANU*RE/(D*DENSA) ENDIF - - END SUBROUTINE TERMINL - - - SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) + END SUBROUTINE TERMINL + + SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY - !!! BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD - !!! AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT, + !!! VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY + !!! BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD + !!! AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT, !!! AND IF THE STONE IS IN WET OR DRY GROWTH REGIME !!! !!! INPUT: PC fraction of updraft water that is frozen @@ -1122,7 +1346,7 @@ SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) INTEGER ITYPE !local variables REAL RV, ALV, ALS, RATIO - DATA RV/461.48/,ALV/2500000./,ALS/2836050./ + DATA RV/461.48/,ALV/2500000./,ALS/2836050./ REAL ESAT, RHOKOR, ESATW, RHOOMGW, ESATI, RHOOMGI, RHOOMG !!! FOR HAILSTONE: FIRST TEST IF STONE IS IN WET OR DRY GROWTH @@ -1133,8 +1357,8 @@ SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) ESAT=611.*EXP(ALS/RV*(RATIO-1./TS)) ENDIF RHOKOR=ESAT/(RV*TS) - - !!! NOW FOR THE AMBIENT/IN-CLOUD CONDITIONS + + !!! NOW FOR THE AMBIENT/IN-CLOUD CONDITIONS ESATW=611.*EXP(ALV/RV*(RATIO-1./TC)) RHOOMGW=ESATW/(RV*TC) ESATI=611.*EXP(ALS/RV*(RATIO-1./TC)) @@ -1142,28 +1366,26 @@ SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) !RHOOMG=PC*(RHOOMGI-RHOOMGW)+RHOOMGW RHOOMG = RHOOMGI !done as in hailtraj.f - !!! CALC THE DIFFERENCE(KG/M3): <0 FOR CONDENSATION, + !!! CALC THE DIFFERENCE(KG/M3): <0 FOR CONDENSATION, !!! >0 FOR EVAPORATION - DELRW=(RHOKOR-RHOOMG) + DELRW=(RHOKOR-RHOOMG) - END SUBROUTINE VAPORCLOSE - - + END SUBROUTINE VAPORCLOSE - SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& - TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,ITYPE,DELRW) + SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& + TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,ITYPE,DELRW) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! CALC THE STONE'S INCREASE IN MASS + !!! CALC THE STONE'S INCREASE IN MASS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + IMPLICIT NONE REAL*8 D REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DI,ANU,RE,AE, & TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,DELRW - INTEGER ITYPE + INTEGER ITYPE !local variables REAL PI, D0, GMW2, GMI2, EW, EI,DGMV - REAL DENSEL, DENSELI, DENSELW + REAL DENSEL, DENSELI, DENSELW REAL DC !MEAN CLOUD DROPLET DIAMETER (MICRONS, 1E-6M) REAL VOLL, VOLT !VOLUME OF NEW LAYER, TOTAL (M3) REAL VOL1, DGMW_NOSOAK, SOAK, SOAKM @@ -1173,19 +1395,19 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !!! CALCULATE THE DIFFUSIVITY DI (m2/s) D0=0.226*1.E-4 ! change to m2/s, not cm2/s DI=D0*(TC/273.155)**1.81*(100000./P) - - !!! COLLECTION EFFICIENCY FOR WATER AND ICE - !EW=1.0 - !!!! IF TS WARMER THAN -5C THEN ACCRETE ALL THE ICE (EI=1.0) - !!!! OTHERWISE EI=0.21 + + !!! COLLECTION EFFICIENCY FOR WATER AND ICE + !EW=1.0 + !!!! IF TS WARMER THAN -5C THEN ACCRETE ALL THE ICE (EI=1.0) + !!!! OTHERWISE EI=0.21 !IF(TS.GE.268.15)THEN ! EI=1.00 !ELSE ! EI=0.21 !ENDIF - !!! COLLECTION EFFICIENCY FOR WATER AND ICE - EW=1.0 + !!! COLLECTION EFFICIENCY FOR WATER AND ICE + EW=1.0 !!! Linear function for ice accretion efficiency IF (TC .GE. 273.155) THEN EI=1.00 @@ -1196,13 +1418,13 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& ENDIF !!! CALCULATE THE VENTILATION COEFFICIENT - NEEDED FOR GROWTH FROM VAPOR - !The coefficients in the ventilation coefficient equations have been + !The coefficients in the ventilation coefficient equations have been !experimentally derived, and are expecting cal-C-g units. Do some conversions. DENSAC = DENSA * (1.E3) * (1.E-6) - !dynamic viscosity + !dynamic viscosity ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5 !!! CALCULATE THE REYNOLDS NUMBER - unitless - RE=D*VT*DENSAC/ANU + RE=D*VT*DENSAC/ANU E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv) !!! SELECT APPROPRIATE VALUES OF AE ACCORDING TO RE IF(RE.LT.6000.0)THEN @@ -1214,42 +1436,42 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& ENDIF - !!! CALC HAILSTONE'S MASS (GM), MASS OF WATER (GMW) AND THE + !!! CALC HAILSTONE'S MASS (GM), MASS OF WATER (GMW) AND THE !!! MASS OF ICE IN THE STONE (GMI) GM=PI/6.*(D**3.)*DENSE GMW=FW*GM GMI=GM-GMW - + !!! STORE THE MASS GM1=GM - - !!! NEW MASS GROWTH CALCULATIONS WITH VARIABLE RIME + + !!! NEW MASS GROWTH CALCULATIONS WITH VARIABLE RIME !!! LAYER DENSITY BASED ON ZIEGLER ET AL. (1983) - + !!! CALCULATE INCREASE IN MASS DUE INTERCEPTED CLD WATER, USE !!! ORIGINAL DIAMETER GMW2=GMW+SEKDEL*(PI/4.*D**2.*VT*XW*EW) - DGMW=GMW2-GMW + DGMW=GMW2-GMW GMW=GMW2 !!! CALCULATE THE INCREASE IN MASS DUE INTERCEPTED CLOUD ICE GMI2=GMI+SEKDEL*(PI/4.*D**2.*VT*XI*EI) - DGMI=GMI2-GMI + DGMI=GMI2-GMI GMI=GMI2 - - !!! CALCULATE INCREASE IN MASS DUE TO SUBLIMATION/CONDENSATION OF + + !!! CALCULATE INCREASE IN MASS DUE TO SUBLIMATION/CONDENSATION OF !!! WATER VAPOR DGMV = SEKDEL*2*PI*D*AE*DI*DELRW IF (DGMV .LT. 0) DGMV=0 - !!! CALCULATE THE TOTAL MASS CHANGE + !!! CALCULATE THE TOTAL MASS CHANGE DGM=DGMW+DGMI+DGMV !Dummy arguments in the event of no water, vapor, or ice growth DENSELW = 900. DENSELI = 700. !!! CALCULATE DENSITY OF NEW LAYER, DEPENDS ON FW AND ITYPE IF (ITYPE.EQ.1) THEN !DRY GROWTH - !If hailstone encountered supercooled water, calculate new layer density + !If hailstone encountered supercooled water, calculate new layer density ! using Macklin form IF ((DGMW.GT.0).OR.(DGMV.GT.0)) THEN !MEAN CLOUD DROPLET RADIUS, ASSUME CLOUD DROPLET CONC OF 3E8 M-3 (300 CM-3) @@ -1299,9 +1521,9 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& VIMP = 0.57*VT ENDIF ENDIF - + !RIME LAYER DENSITY, HEYMSFIELD AND PFLAUM 1985 FORM - TSCELSIUS = TS - 273.16 + TSCELSIUS = TS - 273.16 AFACTOR = -DC*VIMP/TSCELSIUS IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN DENSELW = 0.30*(AFACTOR)**0.44 @@ -1321,7 +1543,7 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !Ice collection main source of growth, so set new density layer DENSELI = 700. ENDIF - + !All liquid water contributes to growth, none is soaked into center. DGMW_NOSOAK = DGMW !All liquid water contributes to growth, ! none of it is soaked into center. @@ -1329,7 +1551,7 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& ELSE !WET GROWTH !Collected liquid water can soak into the stone before freezing, ! increasing mass and density but leaving volume constant. - !Volume of current drop, before growth + !Volume of current drop, before growth VOL1 = GM/DENSE !Difference b/w mass of stone if density is 900 kg/m3, and ! current mass @@ -1340,15 +1562,15 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& IF (SOAKM.GT.SOAK) SOAKM=SOAK GM = GM+SOAKM !Mass of current drop, plus soaking !New density of current drop, including soaking but before growth - DENSE = GM/VOL1 + DENSE = GM/VOL1 !Mass increment of liquid water growth that doesn't ! include the liquid water we just soaked into the stone. DGMW_NOSOAK = DGMW - SOAKM - + !Whatever growth does occur has high density DENSELW = 900. !KG M-3 DENSELI = 900. - + ENDIF !!!VOLUME OF NEW LAYER @@ -1367,18 +1589,16 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& VOLT = VOLL + GM/DENSE !DENSE = (GM+DGM) / VOLT DENSE = (GM+DGMI+DGMV+DGMW_NOSOAK) / VOLT - !D=D+SEKDEL*0.5*VT/DENSE*(XW*EW+XI*EI) + !D=D+SEKDEL*0.5*VT/DENSE*(XW*EW+XI*EI) GM = GM+DGMI+DGMW_NOSOAK+DGMV - D = ( (6*GM) / (PI*DENSE) )**0.33333333 - - END SUBROUTINE MASSAGR + D = ( (6*GM) / (PI*DENSE) )**0.33333333 + END SUBROUTINE MASSAGR - - SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & + SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,ITYPE,P) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! CALCULATE HAILSTONE'S HEAT BUDGET + !!! CALCULATE HAILSTONE'S HEAT BUDGET !!! See Rasmussen and Heymsfield 1987; JAS !!! Original Hailcast's variable units !!! TS - Celsius @@ -1398,7 +1618,7 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & REAL TS,TSm1,TSm2,FW,TC,VT,DELRW,DENSA,GM1,GM,DGM,DGMW,DGMV, & DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,P INTEGER ITYPE - + REAL RV, RD, G, PI, ALF, ALV, ALS, CI, CW, AK REAL H, AH, TCC, TSC, DELRWC, DENSAC, TDIFF REAL DMLT @@ -1406,7 +1626,7 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & DATA RV/461.48/,RD/287.04/,G/9.78956/ DATA PI/3.141592654/,ALF/79.7/,ALV/597.3/ DATA ALS/677.0/,CI/0.5/,CW/1./ - + !Convert values to non-SI units here TSC = TS - 273.155 TSCm1 = TSm1 - 273.155 @@ -1417,14 +1637,14 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & !DI still in cm2/sec - !!! CALCULATE THE CONSTANTS + !!! CALCULATE THE CONSTANTS AK=(5.8+0.0184*TCC)*1.E-5 !thermal conductivity - cal/(cm*sec*K) !dynamic viscosity - calculated in MASSAGR !ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5 !!! CALCULATE THE REYNOLDS NUMBER - unitless - !RE=D*VT*DENSAC/ANU - calculated in MASSAGR - + !RE=D*VT*DENSAC/ANU - calculated in MASSAGR + H=(0.71)**(0.333333333)*(RE**0.50) !ventilation coefficient heat (fh) !E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv) @@ -1440,12 +1660,12 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & !AE=(0.57+9.0E-6*RE)*E calculated in MASSAGR ENDIF - !!! FOR DRY GROWTH FW=0, CALCULATE NEW TS, ITIPE=1 + !!! FOR DRY GROWTH FW=0, CALCULATE NEW TS, ITIPE=1 !!! FOR WET GROWTH TS=0, CALCULATE NEW FW, ITIPE=2 IF(ITYPE.EQ.1) THEN - !!! DRY GROWTH; CALC NEW TEMP OF THE STONE + !!! DRY GROWTH; CALC NEW TEMP OF THE STONE !Original Hailcast algorithm (no time differencing) !TSC=TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* & ! (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & @@ -1454,17 +1674,17 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ & DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)) + & 0.2*TSCm1 + 0.2*TSCm2 - + TS = TSC+273.155 - IF (TS.GE.273.155) THEN + IF (TS.GE.273.155) THEN TS=273.155 ENDIF - TDIFF = ABS(TS-273.155) + TDIFF = ABS(TS-273.155) IF (TDIFF.LE.1.E-6) ITYPE=2 !NOW IN WET GROWTH - + ELSE IF (ITYPE.EQ.2) THEN - !!! WET GROWTH; CALC NEW FW - + !!! WET GROWTH; CALC NEW FW + IF (TCC.LT.0.) THEN !Original Hailcast algorithm FW=FW-FW*DGM/GM1+SEKDEL/(GM1*ALF)* & @@ -1476,25 +1696,23 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, & DGMW/SEKDEL*CW*TCC) / ALF FW = (FW*GM + DMLT) / GM ENDIF - + IF(FW.GT.1.)FW=1. IF(FW.LT.0.)FW=0. !IF ALL OUR ACCRETED WATER WAS FROZEN, WE ARE BACK IN DRY GROWTH IF(FW.LE.1.E-6) THEN - ITYPE=1 + ITYPE=1 ENDIF - - ENDIF - END SUBROUTINE HEATBUD + ENDIF + END SUBROUTINE HEATBUD - - SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT) + SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! TEST IF AMOUNT OF WATER ON SURFACE EXCEEDS CRTICAL LIMIT- - !!! IF SO INVOKE SHEDDING SCHEME + !!! TEST IF AMOUNT OF WATER ON SURFACE EXCEEDS CRTICAL LIMIT- + !!! IF SO INVOKE SHEDDING SCHEME !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE @@ -1508,9 +1726,9 @@ SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT) !GMI=(GM-WATER) !KG !REMAIN = CRIT*GM - ! CALC CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S - ! SURFACE - !CRIT=0.268+0.1389*GMI + ! CALC CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S + ! SURFACE + !CRIT=0.268+0.1389*GMI !CRIT=0.268*1.E-3 + 0.1389*1.E-3*GMI !mass now in kg instead of g !CRIT = 1.0E-10 !CRIT - now passed from main subroutine @@ -1518,21 +1736,20 @@ SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT) WAT=WATER-CRIT GM=GM-WAT FW=(CRIT)/GM - + IF(FW.GT.1.0) FW=1.0 IF(FW.LT.0.0) FW=0.0 - ! RECALCULATE DIAMETER AFTER SHEDDING + ! RECALCULATE DIAMETER AFTER SHEDDING ! Assume density remains the same D=(6.*GM/(PI*DENSE))**(0.333333333) - END SUBROUTINE BREAKUP - - - SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) + END SUBROUTINE BREAKUP + + SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! This is a spherical hail melting estimate based on the Goyer - !!! et al. (1969) eqn (3). The depth of the warm layer, estimated - !!! terminal velocity, and mean temperature of the warm layer are + !!! This is a spherical hail melting estimate based on the Goyer + !!! et al. (1969) eqn (3). The depth of the warm layer, estimated + !!! terminal velocity, and mean temperature of the warm layer are !!! used. DRB. 11/17/2003. !!! !!! INPUT: TLAYER mean sub-cloud layer temperature (K) @@ -1540,7 +1757,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) !!! RLAYER mean sub-cloud layer mixing ratio (kg/kg) !!! VT terminal velocity of stone (m/s) !!! OUTPUT: D diameter (m) - !!! + !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE @@ -1553,7 +1770,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) tres, re, delt, esenv, rhosenv, essfc, rhosfc, dsig, & dmdt, mass, massorg, newmass, gamma, r, rho INTEGER wcnt - + !Convert temp to Celsius, calculate dewpoint in celsius eps = 0.622 tclayer = TLAYER - 273.155 @@ -1561,20 +1778,20 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) b = 5.42E3 tdclayer = b / LOG(a*eps / (rlayer*player)) hplayer = player / 100. - + !Calculate partial vapor pressure eenv = (player*rlayer) / (rlayer+eps) eenv = eenv / 100. !convert to mb - + !Estimate wet bulb temperature (C) gamma = 6.6E-4*player delta = (4098.0*eenv)/((tdclayer+237.7)*(tdclayer+237.7)) wetbulb = ((gamma*tclayer)+(delta*tdclayer))/(gamma+delta) - + !Iterate to get exact wet bulb wcnt = 0 DO WHILE (wcnt .lt. 11) - ewet = 6.108*(exp((17.27*wetbulb)/(237.3 + wetbulb))) + ewet = 6.108*(exp((17.27*wetbulb)/(237.3 + wetbulb))) de = (0.0006355*hplayer*(tclayer-wetbulb))-(ewet-eenv) der= (ewet*(.0091379024 - (6106.396/(273.155+wetbulb)**2))) & - (0.0006355*hplayer) @@ -1585,7 +1802,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) EXIT ENDIF ENDDO - + wetbulbk = wetbulb + 273.155 !convert to K ka = .02 ! thermal conductivity of air lf = 3.34e5 ! latent heat of melting/fusion @@ -1596,16 +1813,16 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) rv = 1004. - 287. ! gas constant for water vapor rhoice = 917.0 ! density of ice (kg/m**3) r = D/2. ! radius of stone (m) - + !Compute residence time in warm layer tres = LDEPTH / VT - + !Calculate dmdt based on eqn (3) of Goyer et al. (1969) !Reynolds number...from pg 317 of Atmo Physics (Salby 1996) !Just use the density of air at 850 mb...close enough. rho = 85000./(287.*TLAYER) re = rho*r*VT*.01/1.7e-5 - + !Temperature difference between environment and hailstone surface delt = wetbulb !- 0.0 !assume stone surface is at 0C !wetbulb is in Celsius @@ -1624,14 +1841,12 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT) dmdt = (-1.7*pi*r*(re**0.5)/lf)*((ka*delt)+((lv-lf)*dv*dsig)) IF (dmdt.gt.0.) dmdt = 0 mass = dmdt*tres - + !Find the new hailstone diameter massorg = 1.33333333*pi*r*r*r*rhoice newmass = massorg + mass if (newmass.lt.0.0) newmass = 0.0 D = 2.*(0.75*newmass/(pi*rhoice))**0.333333333 - END SUBROUTINE MELT + END SUBROUTINE MELT - END MODULE module_diag_hailcast - From da64dfd5d3f2a934b09dc9ac59c2a42bf7a4f511 Mon Sep 17 00:00:00 2001 From: Yunheng Wang Date: Sat, 12 Mar 2022 15:36:14 +0000 Subject: [PATCH 3/6] Tested with Thompson and NSSL 2mom MP schemes --- driver/fvGFS/fv_nggps_diag.F90 | 59 +++++----- model/fv_arrays.F90 | 8 +- tools/module_diag_hailcast.F90 | 201 ++++++++++++++++----------------- 3 files changed, 132 insertions(+), 136 deletions(-) diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index ee1355686..f70211550 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -62,7 +62,7 @@ module fv_nggps_diags_mod ! ! - use mpp_mod, only: mpp_pe, mpp_root_pe,FATAL,mpp_error,NOTE + use mpp_mod, only: mpp_pe, mpp_root_pe,FATAL,mpp_error 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 @@ -434,7 +434,8 @@ subroutine fv_nggps_diag_init(Atm, axes, Time) endif !initialize hailcast - call hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus) + call hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco, & + isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus) ! !------------------------------------ @@ -615,7 +616,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 @@ -745,9 +746,7 @@ subroutine fv_nggps_diag(Atm, zvir, Time) IF ( .not. Atm(n)%flagstruct%hydrostatic .and. do_hailcast ) THEN !--- max hourly hailcast hail diameter if (id_hailcast_dhail > 0) then - call hailcast_compute_dhailmax(isco,ieco,jsco,jeco,istatus) call store_data(id_hailcast_dhail, hailcast_dhail_max, Time, kstt_hc, kend_hc) - if (allocated(hailcast_dhail_max)) deallocate(hailcast_dhail_max) endif !--- max hourly hailcast hail diameter (embryo 1) @@ -800,6 +799,8 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) integer, save :: kdtt = 0 real :: avg_max_length real,dimension(:,:,:),allocatable :: vort + logical :: reset_step + n = 1 ngc = Atm(n)%ng nq = size (Atm(n)%q,4) @@ -823,17 +824,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) - 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 + 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,& @@ -852,15 +853,15 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) if( .not.Atm(n)%flagstruct%hydrostatic ) then call max_vv(isco,ieco,jsco,jeco,npzo,ngc,up2,dn2,Atm(n)%pe,Atm(n)%w) 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 + 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, & @@ -873,7 +874,7 @@ 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 - deallocate (vort) + 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:' @@ -885,8 +886,10 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) !allocate hailcast met field arrays if (do_hailcast) then - call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & - isco,jsco,ieco,jeco,npzo,kdtt,nsteps_per_reset) + reset_step = (mod(kdtt,nsteps_per_reset)==0) + call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, & + isco,jsco,ieco,jeco,isdo,iedo,jsdo,jedo,npzo, & + reset_step,first_time) endif kdtt=kdtt+1 diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index aa9cb6308..dad1b42ff 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -61,10 +61,6 @@ module fv_arrays_mod real :: efx(max_step), efx_sum, efx_nest(max_step), efx_sum_nest, mtq(max_step), mtq_sum integer :: steps - !for hailcast - !integer :: id_hailcast_dhail1, id_hailcast_dhail2, id_hailcast_dhail3, id_hailcast_dhail4, id_hailcast_dhail5 - !integer :: id_hailcast_diam_mean, id_hailcast_diam_std - end type fv_diag_type @@ -248,7 +244,7 @@ module fv_arrays_mod !< monotonicity constraint of Huynh, which is less diffusive !< but more expensive than other monotonic constraints. For hydrostatic simulation, 8 !< (the L04 monotonicity constraint) or 10 are recommended; for - !< nonhydrostatic simulation, the completely unlimited (“linear” + !< nonhydrostatic simulation, the completely unlimited (linear !< or non-monotone) PPM scheme is recommended. If no monotonicity !< constraint is applied, enabling the flux damping !< (do_vort_damp = .true.) is highly recommended to control grid-scale @@ -643,7 +639,7 @@ module fv_arrays_mod logical :: nudge_dz = .false. !< During the adiabatic initialization (na_init > 0), if set !< to .true., delz is nudged back to the value specified in the initial !< conditions, instead of nudging the temperature back to the initial value. - !< Nudging delz is simpler (faster), doesn’t require consideration of the + !< Nudging delz is simpler (faster), doesnt require consideration of the !< virtual temperature effect, and may be more stable. .false.by default. real :: p_ref = 1.E5 !< Surface pressure used to construct a horizontally-uniform reference !< vertical pressure profile, used in some simple physics packages diff --git a/tools/module_diag_hailcast.F90 b/tools/module_diag_hailcast.F90 index bf9177e79..f15a55d6f 100644 --- a/tools/module_diag_hailcast.F90 +++ b/tools/module_diag_hailcast.F90 @@ -2,7 +2,7 @@ !John Henderson AER 20190425 !time management handling from WRF is preserved but commented out -! Y. Wang NSSL/CIWRO 20220310 +! Yunheng Wang NSSL/CIWRO 20220310 ! Remodeled as recommended for PR #164. ! @@ -12,7 +12,9 @@ MODULE module_diag_hailcast use fms_mod, only: check_nml_error !use fv_mp_mod, only: is_master - use constants_mod, only: grav + use constants_mod, only: grav, rdgas + + implicit none logical :: do_hailcast = .false. !This controls whether hailcast is used @@ -40,7 +42,8 @@ MODULE module_diag_hailcast PRIVATE :: hailstone_driver, hailcast_diagnostic_driver CONTAINS - SUBROUTINE hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus) + SUBROUTINE hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco,& + isdo,iedo,jsdo,jedo,nlevs, missing_value, istatus) use diag_manager_mod, only: register_diag_field @@ -48,27 +51,22 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missi integer, intent(in) :: axes(4) type(time_type), intent(in) :: Time INTEGER, INTENT(IN) :: isdo,iedo,jsdo,jedo + INTEGER, INTENT(IN) :: isco,ieco,jsco,jeco INTEGER, INTENT(INOUT) :: nlevs real, INTENT(IN) :: missing_value INTEGER, INTENT(OUT) :: istatus + INTEGER :: i, j + namelist /fv_diagnostics_nml/ do_hailcast integer :: ios, ierr - integer :: unit, f_unit + integer :: unit !namelist file for hailcast -#ifdef INTERNAL_FILE_NML + ! Read Main namelist read (input_nml_file,fv_diagnostics_nml,iostat=ios) ierr = check_nml_error(ios,'fv_diagnostics_nml') -#else - f_unit=open_namelist_file() - rewind (f_unit) - ! Read Main namelist - read (f_unit,fv_diagnostics_nml,iostat=ios) - ierr = check_nml_error(ios,'fv_diagnostics_nml') - call close_file(f_unit) -#endif unit = stdlog() write(unit, nml=fv_diagnostics_nml) @@ -80,7 +78,7 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missi !register hailcast arrays if (do_hailcast) then - id_hailcast_dhail = register_diag_field ( trim(file_name), 'hailcast_dhail_max', axes(1:2), Time, & + id_hailcast_dhail = register_diag_field ( trim(file_name), 'hailcast_dhail_max', axes(1:2), Time, & 'hourly max hail diameter', 'mm', missing_value=missing_value ) id_hailcast_dhail1 = register_diag_field ( trim(file_name), 'hailcast_dhail1_max', axes(1:2), Time, & 'hourly max hail diameter (embryo 1)', 'mm', missing_value=missing_value ) @@ -136,37 +134,44 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missi endif if (.not.allocated(hailcast_dhail1)) then - allocate(hailcast_dhail1(isco:ieco,jsco:jeco), & - hailcast_dhail2(isco:ieco,jsco:jeco), & - hailcast_dhail3(isco:ieco,jsco:jeco), & - hailcast_dhail4(isco:ieco,jsco:jeco), & - hailcast_dhail5(isco:ieco,jsco:jeco), & - hailcast_diam_mean(isco:ieco,jsco:jeco), & + allocate(hailcast_dhail1(isco:ieco,jsco:jeco), & + hailcast_dhail2(isco:ieco,jsco:jeco), & + hailcast_dhail3(isco:ieco,jsco:jeco), & + hailcast_dhail4(isco:ieco,jsco:jeco), & + hailcast_dhail5(isco:ieco,jsco:jeco), & + hailcast_diam_mean(isco:ieco,jsco:jeco), & hailcast_diam_std(isco:ieco,jsco:jeco)) allocate ( hailcast_dhail1_max(isco:ieco,jsco:jeco) ) allocate ( hailcast_dhail2_max(isco:ieco,jsco:jeco) ) allocate ( hailcast_dhail3_max(isco:ieco,jsco:jeco) ) allocate ( hailcast_dhail4_max(isco:ieco,jsco:jeco) ) allocate ( hailcast_dhail5_max(isco:ieco,jsco:jeco) ) + allocate ( hailcast_dhail_max(isco:ieco,jsco:jeco) ) do i=isco,ieco - do j=jsco,jeco - hailcast_dhail1(i,j)=0; hailcast_dhail2(i,j)=0 - hailcast_dhail3(i,j)=0; hailcast_dhail4(i,j)=0 - hailcast_dhail5(i,j)=0 - hailcast_dhail1_max(i,j)=0; hailcast_dhail2_max(i,j)=0 - hailcast_dhail3_max(i,j)=0; hailcast_dhail4_max(i,j)=0 - hailcast_dhail5_max(i,j)=0 - enddo + do j=jsco,jeco + hailcast_dhail1(i,j)=0 + hailcast_dhail2(i,j)=0 + hailcast_dhail3(i,j)=0 + hailcast_dhail4(i,j)=0 + hailcast_dhail5(i,j)=0 + hailcast_dhail1_max(i,j)=0 + hailcast_dhail2_max(i,j)=0 + hailcast_dhail3_max(i,j)=0 + hailcast_dhail4_max(i,j)=0 + hailcast_dhail5_max(i,j)=0 + hailcast_dhail_max(i,j) =0 + enddo enddo endif if (.not.allocated(hailcast_wdur)) then - allocate(hailcast_wdur(isdo:iedo,jsdo:jedo),hailcast_wup_mask(isdo:iedo,jsdo:jedo)) + allocate(hailcast_wdur(isdo:iedo,jsdo:jedo)) + allocate(hailcast_wup_mask(isdo:iedo,jsdo:jedo)) do i=isdo,iedo - do j=jsdo,jedo - hailcast_wdur(i,j)=0 - hailcast_wup_mask(i,j)=0 - enddo + do j=jsdo,jedo + hailcast_wdur(i,j) = 0 + hailcast_wup_mask(i,j)= 0 + enddo enddo endif endif @@ -174,16 +179,21 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isdo,iedo,jsdo,jedo,nlevs, missi RETURN END SUBROUTINE hailcast_init - SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & - isco,jsco,ieco,jeco,npzo, kdtt, nsteps_per_reset) + SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, & + isco,jsco,ieco,jeco,isdo,iedo,jsdo,jedo,npzo, reset_step, first_time) !moved this code from subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) in fv_nggps_diag use fv_arrays_mod, only: fv_atmos_type type(fv_atmos_type), intent(inout) :: Atm integer, INTENT(IN) :: isco, ieco, jsco, jeco, npzo + integer, INTENT(IN) :: isdo, iedo, jsdo, jedo integer, INTENT(IN) :: sphum, liq_wat, ice_wat !< GFDL physics integer, INTENT(IN) :: rainwat, snowwat, graupel - integer, INTENT(IN) :: nsteps_per_reset, kdtt + real, intent(in) :: zvir + LOGICAL, INTENT(IN) :: reset_step + real, INTENT(IN) :: first_time + + INTEGER :: i, j, k !declare temporary hailcast variables real, allocatable :: rhoair_layer(:,:,:), z(:), z_layer(:,:,:), p_layer(:,:,:),zsfc(:,:) @@ -191,25 +201,26 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !automatic 2-D array for terrain height !automatic 1-D array for interface levels - IF (id_hailcast_dhail>0 .or. id_hailcast_dhail1>0 .or. & + IF (id_hailcast_dhail >0 .or. id_hailcast_dhail1>0 .or. & id_hailcast_dhail2>0 .or. id_hailcast_dhail3>0 .or. & id_hailcast_dhail4>0 .or. id_hailcast_dhail5>0 .or. & id_hailcast_diam_mean>0 .or. id_hailcast_diam_std>0) THEN do j=jsco,jeco do i=isco,ieco - if(mod(kdtt,nsteps_per_reset)==0)then - hailcast_dhail1_max(i,j) = 0. - hailcast_dhail2_max(i,j) = 0. - hailcast_dhail3_max(i,j) = 0. - hailcast_dhail4_max(i,j) = 0. - hailcast_dhail5_max(i,j) = 0. - endif - - hailcast_dhail1(i,j) = 0 - hailcast_dhail2(i,j) = 0 - hailcast_dhail3(i,j) = 0 - hailcast_dhail4(i,j) = 0 - hailcast_dhail5(i,j) = 0 + if (reset_step) then + hailcast_dhail1_max(i,j) = 0. + hailcast_dhail2_max(i,j) = 0. + hailcast_dhail3_max(i,j) = 0. + hailcast_dhail4_max(i,j) = 0. + hailcast_dhail5_max(i,j) = 0. + hailcast_dhail_max(i,j) = 0. + endif + + hailcast_dhail1(i,j) = 0. + hailcast_dhail2(i,j) = 0. + hailcast_dhail3(i,j) = 0. + hailcast_dhail4(i,j) = 0. + hailcast_dhail5(i,j) = 0. enddo enddo @@ -238,22 +249,20 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !height of interfaces and layer height z(k)=z(k+1)-Atm%delz(i,j,k) z_layer(i,j,k)=(z(k+1)+z(k))/2 - p_layer(i,j,k)=Atm%delp(i,j,k)*(1.-sum(Atm%q(i,j,k,2:Atm%flagstruct%nwat)))/& - (-Atm%delz(i,j,k)*grav)*rdgas*Atm%pt(i,j,k)*(1.+zvir*Atm%q(i,j,k,sphum)) + p_layer(i,j,k)=Atm%delp(i,j,k)*(1.-sum(Atm%q(i,j,k,2:Atm%flagstruct%nwat)))/ & + (-Atm%delz(i,j,k)*grav)*rdgas*Atm%pt(i,j,k)*(1.+zvir*Atm%q(i,j,k,sphum)) enddo endif enddo !j loop enddo !k loop !call hailcast diagnostic driver once per time step and provide three-dimensional met fields - call hailcast_diagnostic_driver(Atm%pt(isco:ieco,jsco:jeco,1:npzo), Atm%w(isco:ieco,jsco:jeco,1:npzo), p_layer, z_layer, rhoair_layer, Atm%q(isco:ieco,jsco:jeco,1:npzo,:), & !3D fields - Atm%flagstruct%nwat, sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !number of tracer indices, indices + call hailcast_diagnostic_driver(Atm%pt(isco:ieco,jsco:jeco,1:npzo), Atm%w(isco:ieco,jsco:jeco,1:npzo), & + p_layer, z_layer, rhoair_layer, Atm%q(isco:ieco,jsco:jeco,1:npzo,:), & !3D fields + Atm%flagstruct%nwat, sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !number of tracer indices, indices isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo, & !grid dimensions data array (with halo) and physical grid dimensions 1,npzo, zsfc, & !vertical dimensions, terrain height - first_time, Atm%domain) !call hailcast every model step and info for updating haloes - !hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5, & !hailcast embryos sizes - !hailcast_diam_mean, hailcast_diam_std, & !hailcast embryo mean/std - !hailcast_wdur, hailcast_wup_mask) !persistent updraft duration and mask + first_time, Atm%domain) !call hailcast every model step and info for updating haloes do j=jsco,jeco do i=isco,ieco @@ -262,10 +271,13 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & hailcast_dhail3_max(i,j) = max(hailcast_dhail3_max(i,j), hailcast_dhail3(i,j)) hailcast_dhail4_max(i,j) = max(hailcast_dhail4_max(i,j), hailcast_dhail4(i,j)) hailcast_dhail5_max(i,j) = max(hailcast_dhail5_max(i,j), hailcast_dhail5(i,j)) + + hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail1_max(i,j), & + hailcast_dhail2_max(i,j), hailcast_dhail3_max(i,j), & + hailcast_dhail4_max(i,j), hailcast_dhail5_max(i,j)) enddo enddo - !deallocate hailcast met variables deallocate(p_layer) deallocate(z_layer) @@ -277,27 +289,6 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & RETURN END SUBROUTINE hailcast_compute - SUBROUTINE hailcast_compute_dhailmax(isco,ieco,jsco,jeco,istatus) - - INTEGER, INTENT(IN) :: isco,ieco,jsco,jeco - INTEGER, INTENT(OUT) :: istatus - - allocate(hailcast_dhail_max(isco:ieco,jsco:jeco)) - - do j=jsco,jeco - do i=isco,ieco - hailcast_dhail_max(i,j) = 0. - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail1_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail2_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail3_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail4_max(i,j)) - hailcast_dhail_max(i,j) = max(hailcast_dhail_max(i,j), hailcast_dhail5_max(i,j)) - end do - end do - - RETURN - END SUBROUTINE hailcast_compute_dhailmax - SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! 3D temperature, updraft, pressure, height, density and moisture tracers moist_count,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & ! tracer indices is,ie,js,je,ishalo,iehalo,jshalo,jehalo, & ! i,j halo data array dimensions and physical grid dimensions (2 smaller on each end) @@ -336,13 +327,13 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & CHARACTER*512 :: message CHARACTER*256 :: timestr INTEGER :: i,j,k, numlevs - REAL, DIMENSION( is:ie, js:je, k_start:k_end ) :: qr & - , qs & - , qg & - , qv & - , qc & - , qi & - , ptot + REAL, DIMENSION( is:ie, js:je, k_start:k_end ) :: qr & + , qs & + , qg & + , qv & + , qc & + , qi & + , ptot REAL, DIMENSION( ishalo:iehalo, jshalo:jehalo) :: wdur_prev, wup_mask_prev @@ -439,7 +430,7 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! DO i = MAX( its-1 , ids ), MIN( ite+1 , ide-1) DO j = jshalo, jehalo DO i = ishalo, iehalo - wdur_prev(i,j) = hailcast_wdur(i,j) + wdur_prev(i,j) = hailcast_wdur(i,j) wup_mask_prev(i,j) = hailcast_wup_mask(i,j) ENDDO ENDDO @@ -452,12 +443,12 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & DO j = js, je DO i = is, ie hailcast_wup_mask(i,j) = 0 - hailcast_wdur(i,j) = 0 + hailcast_wdur(i,j) = 0 DO k = k_start, k_end - IF ( w(i,j,k) .ge. 10. ) THEN - hailcast_wup_mask(i,j) = 1 - ENDIF + IF ( w(i,j,k) .ge. 10. ) THEN + hailcast_wup_mask(i,j) = 1 + ENDIF ENDDO ENDDO ENDDO @@ -469,8 +460,8 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! Determine updraft duration using updraft masks ! --------------------------------------------------- - IF ( (hailcast_wup_mask(i,j).eq.1) .OR. & - (MAXVAL(wup_mask_prev(i-1:i+1,j-1:j+1)).eq.1) ) THEN + IF ( (hailcast_wup_mask(i,j).eq.1) .OR. & + (MAXVAL(wup_mask_prev(i-1:i+1,j-1:j+1)) .eq. 1) ) THEN hailcast_wdur(i,j) = MAXVAL(wdur_prev(i-1:i+1,j-1:j+1)) + dt ENDIF ENDDO @@ -569,7 +560,12 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & DO i=is, ie DO k=1, numlevs DO j=js, je - qv(i,j,k)=0;qr(i,j,k)=0;qs(i,j,k)=0;qg(i,j,k)=0;qc(i,j,k)=0;qi(i,j,k)=0 + qv(i,j,k)=0 + qr(i,j,k)=0 + qs(i,j,k)=0 + qg(i,j,k)=0 + qc(i,j,k)=0 + qi(i,j,k)=0 qv(i,j,k) = q(i,j,k,sphum) qr(i,j,k) = q(i,j,k,rainwat) qs(i,j,k) = q(i,j,k,snowwat) @@ -642,14 +638,15 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & ! hailcast_dhail2(i,j), hailcast_dhail3(i,j),& ! hailcast_dhail4(i,j), hailcast_dhail5(i,j)) !sample standard deviation - hailcast_diam_std(i,j) = SQRT( ( & - (hailcast_dhail1(i,j)-hailcast_diam_mean(i,j))**2.+& - (hailcast_dhail2(i,j)-hailcast_diam_mean(i,j))**2.+& - (hailcast_dhail3(i,j)-hailcast_diam_mean(i,j))**2.+& - (hailcast_dhail4(i,j)-hailcast_diam_mean(i,j))**2.+& - (hailcast_dhail5(i,j)-hailcast_diam_mean(i,j))**2.)& + hailcast_diam_std(i,j) = SQRT( ( & + (hailcast_dhail1(i,j)-hailcast_diam_mean(i,j))**2.+ & + (hailcast_dhail2(i,j)-hailcast_diam_mean(i,j))**2.+ & + (hailcast_dhail3(i,j)-hailcast_diam_mean(i,j))**2.+ & + (hailcast_dhail4(i,j)-hailcast_diam_mean(i,j))**2.+ & + (hailcast_dhail5(i,j)-hailcast_diam_mean(i,j))**2.) & / 4.0) - if (hailcast_diam_mean(i,j).ne.0.0 .and. hailcast_diam_std(i,j).ne.0.0) print*,'jmh mean std=',i,j,hailcast_diam_mean(i,j),hailcast_diam_std(i,j) + if (hailcast_diam_mean(i,j).ne.0.0 .and. hailcast_diam_std(i,j).ne.0.0) & + print*,'jmh mean std=',i,j,hailcast_diam_mean(i,j),hailcast_diam_std(i,j) ENDDO ENDDO !END IF From 7506f33e7c0f8987b3dfe88e340da3bcf4949d99 Mon Sep 17 00:00:00 2001 From: Yunheng Wang Date: Fri, 18 Mar 2022 21:46:35 +0000 Subject: [PATCH 4/6] Restore to 43f7ed3 --- model/fv_arrays.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index dad1b42ff..27f44c04b 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -244,7 +244,7 @@ module fv_arrays_mod !< monotonicity constraint of Huynh, which is less diffusive !< but more expensive than other monotonic constraints. For hydrostatic simulation, 8 !< (the L04 monotonicity constraint) or 10 are recommended; for - !< nonhydrostatic simulation, the completely unlimited (linear + !< nonhydrostatic simulation, the completely unlimited (“linear” !< or non-monotone) PPM scheme is recommended. If no monotonicity !< constraint is applied, enabling the flux damping !< (do_vort_damp = .true.) is highly recommended to control grid-scale @@ -639,7 +639,7 @@ module fv_arrays_mod logical :: nudge_dz = .false. !< During the adiabatic initialization (na_init > 0), if set !< to .true., delz is nudged back to the value specified in the initial !< conditions, instead of nudging the temperature back to the initial value. - !< Nudging delz is simpler (faster), doesnt require consideration of the + !< Nudging delz is simpler (faster), doesn’t require consideration of the !< virtual temperature effect, and may be more stable. .false.by default. real :: p_ref = 1.E5 !< Surface pressure used to construct a horizontally-uniform reference !< vertical pressure profile, used in some simple physics packages From e67bb423dd4aa78e4531475d13f06c94fc1ab08d Mon Sep 17 00:00:00 2001 From: Yunheng Wang Date: Mon, 21 Mar 2022 19:17:56 +0000 Subject: [PATCH 5/6] Removed hailcast output dependency on other max/min hourly fields in diag_table --- driver/fvGFS/fv_nggps_diag.F90 | 7 ++----- tools/module_diag_hailcast.F90 | 33 +++++++++++++++++++++++++++------ 2 files changed, 29 insertions(+), 11 deletions(-) diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index f70211550..1d012bb28 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -799,8 +799,6 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) integer, save :: kdtt = 0 real :: avg_max_length real,dimension(:,:,:),allocatable :: vort - logical :: reset_step - n = 1 ngc = Atm(n)%ng nq = size (Atm(n)%q,4) @@ -874,6 +872,7 @@ 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) else print *,'Missing max/min hourly field in diag_table' @@ -886,13 +885,11 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) !allocate hailcast met field arrays if (do_hailcast) then - reset_step = (mod(kdtt,nsteps_per_reset)==0) call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, & isco,jsco,ieco,jeco,isdo,iedo,jsdo,jedo,npzo, & - reset_step,first_time) + Time_step_atmos,avg_max_length) endif - kdtt=kdtt+1 end subroutine fv_nggps_tavg ! subroutine store_data(id, work, Time, nstt, nend) diff --git a/tools/module_diag_hailcast.F90 b/tools/module_diag_hailcast.F90 index f15a55d6f..26b073d90 100644 --- a/tools/module_diag_hailcast.F90 +++ b/tools/module_diag_hailcast.F90 @@ -7,7 +7,7 @@ ! MODULE module_diag_hailcast - USE time_manager_mod, ONLY: time_type + USE time_manager_mod, ONLY: time_type, get_time use mpp_mod, only: mpp_error, FATAL, input_nml_file, stdlog, mpp_pe, mpp_root_pe use fms_mod, only: check_nml_error !use fv_mp_mod, only: is_master @@ -38,8 +38,12 @@ MODULE module_diag_hailcast integer :: kend_hc1,kend_hc2,kend_hc3,kend_hc4,kend_hc5 integer :: kstt_hcd, kend_hcd, kstt_hcm, kend_hcm + real :: time_step + integer :: kdtt = 0 + PRIVATE :: INTERPP, INTERP, TERMINL, VAPORCLOSE, MASSAGR, HEATBUD, BREAKUP, MELT PRIVATE :: hailstone_driver, hailcast_diagnostic_driver + PRIVATE :: time_step, kdtt CONTAINS SUBROUTINE hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco,& @@ -180,20 +184,24 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco,& END SUBROUTINE hailcast_init SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, & - isco,jsco,ieco,jeco,isdo,iedo,jsdo,jedo,npzo, reset_step, first_time) + isco,jsco,ieco,jeco,isdo,iedo,jsdo,jedo,npzo, Time_step_atmos, avg_max_length) !moved this code from subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir) in fv_nggps_diag use fv_arrays_mod, only: fv_atmos_type type(fv_atmos_type), intent(inout) :: Atm + type(time_type), intent(in) :: Time_step_atmos + real, INTENT(IN) :: avg_max_length integer, INTENT(IN) :: isco, ieco, jsco, jeco, npzo integer, INTENT(IN) :: isdo, iedo, jsdo, jedo integer, INTENT(IN) :: sphum, liq_wat, ice_wat !< GFDL physics integer, INTENT(IN) :: rainwat, snowwat, graupel real, intent(in) :: zvir - LOGICAL, INTENT(IN) :: reset_step - real, INTENT(IN) :: first_time INTEGER :: i, j, k + logical, save :: first_call=.true. + integer :: seconds, days + LOGICAL :: reset_step + integer :: nsteps_per_reset !declare temporary hailcast variables real, allocatable :: rhoair_layer(:,:,:), z(:), z_layer(:,:,:), p_layer(:,:,:),zsfc(:,:) @@ -201,6 +209,17 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, z !automatic 2-D array for terrain height !automatic 1-D array for interface levels + if (first_call) then + call get_time (Time_step_atmos, seconds, days) + !write(*,*) "setting hailcast: seconds= ", seconds + time_step = seconds + first_call= .false. + kdtt=0 + endif + + nsteps_per_reset = nint(avg_max_length/time_step) + reset_step = (mod(kdtt,nsteps_per_reset)==0) + IF (id_hailcast_dhail >0 .or. id_hailcast_dhail1>0 .or. & id_hailcast_dhail2>0 .or. id_hailcast_dhail3>0 .or. & id_hailcast_dhail4>0 .or. id_hailcast_dhail5>0 .or. & @@ -262,7 +281,7 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, z Atm%flagstruct%nwat, sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, & !number of tracer indices, indices isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo, & !grid dimensions data array (with halo) and physical grid dimensions 1,npzo, zsfc, & !vertical dimensions, terrain height - first_time, Atm%domain) !call hailcast every model step and info for updating haloes + time_step, Atm%domain) !call hailcast every model step and info for updating haloes do j=jsco,jeco do i=isco,ieco @@ -286,6 +305,8 @@ SUBROUTINE hailcast_compute(Atm,sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, z deallocate(zsfc) END IF + kdtt=kdtt+1 + RETURN END SUBROUTINE hailcast_compute @@ -1123,7 +1144,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& ENDIF !end check for melting call - !Check to make sure embryo didn’t just immediately fall out + !Check to make sure embryo didnt just immediately fall out IF (sec .LT. 60) D=0.0 !Check to make sure something hasn't gone horribly wrong From 018b4809ba270c2eb54b419a34f61dd819dec5e3 Mon Sep 17 00:00:00 2001 From: Yunheng Wang Date: Tue, 22 Mar 2022 19:10:36 +0000 Subject: [PATCH 6/6] Improved code based PR comments --- driver/fvGFS/fv_nggps_diag.F90 | 40 ++++- tools/module_diag_hailcast.F90 | 284 ++++++++++++--------------------- 2 files changed, 141 insertions(+), 183 deletions(-) diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index 1d012bb28..fe5210617 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -62,7 +62,7 @@ module fv_nggps_diags_mod ! ! - 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 @@ -74,8 +74,23 @@ module fv_nggps_diags_mod use fv_diagnostics_mod, only: range_check, dbzcalc,max_vv,get_vorticity, & max_uh,max_vorticity,bunkers_vector, & helicity_relative_CAPS,max_vorticity_hy1 - use module_diag_hailcast 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 @@ -152,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 @@ -433,10 +462,11 @@ subroutine fv_nggps_diag_init(Atm, axes, Time) endif - !initialize hailcast - call hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco, & + 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 ! !------------------------------------ ! use wrte grid component for output diff --git a/tools/module_diag_hailcast.F90 b/tools/module_diag_hailcast.F90 index 26b073d90..e11cec755 100644 --- a/tools/module_diag_hailcast.F90 +++ b/tools/module_diag_hailcast.F90 @@ -8,8 +8,8 @@ MODULE module_diag_hailcast USE time_manager_mod, ONLY: time_type, get_time - use mpp_mod, only: mpp_error, FATAL, input_nml_file, stdlog, mpp_pe, mpp_root_pe - use fms_mod, only: check_nml_error + use mpp_mod, only: mpp_pe, mpp_root_pe, mpp_error, FATAL !, input_nml_file, stdlog + !use fms_mod, only: check_nml_error !use fv_mp_mod, only: is_master use constants_mod, only: grav, rdgas @@ -62,19 +62,21 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco,& INTEGER :: i, j - 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 + !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 if (mpp_pe() == mpp_root_pe()) then print*, 'do_hailcast = ', do_hailcast @@ -335,10 +337,6 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & REAL, DIMENSION( is:ie, js:je), & INTENT(IN ) :: terr_z -! REAL, INTENT(IN ),OPTIONAL :: haildt -! REAL, INTENT(IN ),OPTIONAL :: curr_secs -! REAL, INTENT(INOUT),OPTIONAL :: haildtacttime -! INTEGER, INTENT(IN ) :: itimestep REAL, INTENT(IN ) :: dt TYPE(domain2d), INTENT(INOUT ) :: domainid @@ -360,101 +358,15 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & REAL :: dhail1,dhail2,dhail3,dhail4,dhail5 -! ! Timing -! ! ------ -! TYPE(WRFU_Time) :: hist_time, aux2_time, CurrTime, StartTime -! TYPE(WRFU_TimeInterval) :: dtint, histint, aux2int -! LOGICAL :: is_after_history_dump, is_output_timestep, is_first_timestep -! LOGICAL :: doing_adapt_dt, run_param, decided -! INTEGER :: stephail, itimestep_basezero - - -! ! Chirp the routine name for debugging purposes -! ! --------------------------------------------- -! write ( message, * ) 'inside hailcast_diagnostics_driver' -! CALL wrf_debug( 100 , message ) - -! ! Get timing info -! ! Want to know if when the last history output was -! ! Check history and auxhist2 alarms to check last ring time and how often -! ! they are set to ring -! ! ----------------------------------------------------------------------- -! CALL WRFU_ALARMGET( grid%alarms( HISTORY_ALARM ), prevringtime=hist_time, & -! ringinterval=histint) -! CALL WRFU_ALARMGET( grid%alarms( AUXHIST2_ALARM ), prevringtime=aux2_time, & -! ringinterval=aux2int) - -! ! Get domain clock -! ! ---------------- -! CALL domain_clock_get ( grid, current_time=CurrTime, & -! simulationStartTime=StartTime, & -! current_timestr=timestr, time_step=dtint ) - -! ! Set some booleans for use later -! ! Following uses an overloaded .lt. -! ! --------------------------------- -! is_after_history_dump = ( Currtime .lt. hist_time + dtint ) - -! ! Following uses an overloaded .ge. -! ! --------------------------------- -! is_output_timestep = (Currtime .ge. hist_time + histint - dtint .or. & -! Currtime .ge. aux2_time + aux2int - dtint ) -! write ( message, * ) 'is output timestep? ', is_output_timestep -! CALL wrf_debug( 100 , message ) - -! ! Following uses an overloaded .eq. -! ! --------------------------------- -! is_first_timestep = ( Currtime .eq. StartTime + dtint ) + !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ - !master = (mpp_pe()==mpp_root_pe()) .or. is_master() - - !FV3: No need to reset arrays to zero at dump time, thus no need to keep track of dump times - ! - ! After each history dump, reset max/min value arrays - ! ---------------------------------------------------------------------- -! IF ( is_after_history_dump ) THEN - !! DO j = jsc, jec !jts, jte - !! DO i = isc, iec !its, ite - !! hailcast_dhail1(i,j) = 0. - !! hailcast_dhail2(i,j) = 0. - !! hailcast_dhail3(i,j) = 0. - !!! hailcast_dhail4(i,j) = 0. - !! hailcast_dhail5(i,j) = 0. - !! ENDDO - !! ENDDO - !ENDIF ! is_after_history_dump - - ! We need to do some neighboring gridpoint comparisons for the updraft - ! duration calculations; set i,j start and end values so we don't go off - ! the edges of the domain. Updraft duration on domain edges will always be 0. - ! ---------------------------------------------------------------------- - -! !WRF config_flags contains: open_xs, specified, nested,open_xe,open_yx,open_ye,periodic_x -! IF ( config_flags%open_xs .OR. config_flags%specified .OR. & -! config_flags%nested) i_start = MAX( ids+1, its ) -! IF ( config_flags%open_xe .OR. config_flags%specified .OR. & -! config_flags%nested) i_end = MIN( ide-2, ite ) !-1 -! IF ( config_flags%open_ys .OR. config_flags%specified .OR. & -! config_flags%nested) j_start = MAX( jds+1, jts ) -! IF ( config_flags%open_ye .OR. config_flags%specified .OR. & -! config_flags%nested) j_end = MIN( jde-2, jte ) !-1 -! IF ( config_flags%periodic_x ) i_start = its -! IF ( config_flags%periodic_x ) i_end = ite - - - ! Make a copy of the updraft duration, mask variables - ! --------------------------------------------------- - !from WRF: - !!wdur_prev(:,:) = grid%hailcast_wdur(:,:) - !!wup_mask_prev(:,:) = grid%hailcast_wup_mask(:,:) - !DO j = MAX( jts-1 , jds ), MIN( jte+1 , jde-1 ) - ! DO i = MAX( its-1 , ids ), MIN( ite+1 , ide-1) DO j = jshalo, jehalo DO i = ishalo, iehalo wdur_prev(i,j) = hailcast_wdur(i,j) wup_mask_prev(i,j) = hailcast_wup_mask(i,j) ENDDO ENDDO + !update halo: call mpp_update_domains(wup_mask_prev,domainid) call mpp_update_domains(wdur_prev,domainid) @@ -624,8 +536,6 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & dhail3, dhail4, & dhail5 ) - !diag_table should select largest value of these variables during requested period - !thus do not need to accumulate largest value between data dumps !IF (dhail1 .gt. hailcast_dhail1(i,j)) THEN hailcast_dhail1(i,j) = dhail1 !ENDIF @@ -664,10 +574,9 @@ SUBROUTINE hailcast_diagnostic_driver ( t, w, p, z, rho, q, & (hailcast_dhail2(i,j)-hailcast_diam_mean(i,j))**2.+ & (hailcast_dhail3(i,j)-hailcast_diam_mean(i,j))**2.+ & (hailcast_dhail4(i,j)-hailcast_diam_mean(i,j))**2.+ & - (hailcast_dhail5(i,j)-hailcast_diam_mean(i,j))**2.) & - / 4.0) - if (hailcast_diam_mean(i,j).ne.0.0 .and. hailcast_diam_std(i,j).ne.0.0) & - print*,'jmh mean std=',i,j,hailcast_diam_mean(i,j),hailcast_diam_std(i,j) + (hailcast_dhail5(i,j)-hailcast_diam_mean(i,j))**2.) / 4.0) + !if (hailcast_diam_mean(i,j).ne.0.0 .and. hailcast_diam_std(i,j).ne.0.0) & + ! print*,'jmh mean std=',i,j,hailcast_diam_mean(i,j),hailcast_diam_std(i,j) ENDDO ENDDO !END IF @@ -1023,7 +932,7 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& !print *, 'INTERP VU: ', VU, P !Outside pressure levels? If so, exit loop - IF (IFOUT.EQ.1) GOTO 100 + IF (IFOUT.EQ.1) EXIT !Sine wave multiplier on updraft strength IF (SEC .GT. 0.0 .AND. SEC .LT. RTIME) THEN @@ -1089,16 +998,15 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& ENDIF !!! Has stone reached below cloud base? - !IF (Z .LE. 0) GOTO 200 IF (Z .LE. ZBAS) then - GOTO 200 + EXIT endif !calculate ice-only diameter size D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 !Has the stone entirely melted and it's below the freezing level? - IF ((D_ICE .LT. 1.E-8) .AND. (TC.GT.273.155)) GOTO 300 + IF ((D_ICE .LT. 1.E-8) .AND. (TC.GT.273.155)) EXIT !move values to previous timestep value TSm1 = TS @@ -1106,9 +1014,9 @@ SUBROUTINE hailstone_driver ( ii,jj,TCA, h1d, ht, PA, rho1d,& ENDDO !end cloud lifetime loop -100 CONTINUE !outside pressure levels in model -200 CONTINUE !stone reached surface -300 CONTINUE !stone has entirely melted and is below freezing level + !100 CONTINUE !outside pressure levels in model + !200 CONTINUE !stone reached surface + !300 CONTINUE !stone has entirely melted and is below freezing level !!!!!!!!!!!!!!!!!! 9. MELT STONE BELOW CLOUD !!!!!!!!!!!!!!!!!!!! !Did the stone shoot out the top of the storm? @@ -1292,24 +1200,27 @@ SUBROUTINE INTERP(AA,A,P,IFOUT,PA,ITEL) END SUBROUTINE INTERP SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!! - !!!! INTERP: Calculate terminal velocity of the hailstone - !!!! - !!!! INPUT: DENSA density of updraft air (kg/m3) - !!!! DENSE density of hailstone - !!!! D diameter of hailstone (m) - !!!! TC updraft temperature (K) - !!!! OUTPUT:VT hailstone terminal velocity (m/s) - !!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!! + !!!! INTERP: Calculate terminal velocity of the hailstone + !!!! + !!!! INPUT: DENSA density of updraft air (kg/m3) + !!!! DENSE density of hailstone + !!!! D diameter of hailstone (m) + !!!! TC updraft temperature (K) + !!!! OUTPUT:VT hailstone terminal velocity (m/s) + !!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE - REAL*8 D - REAL DENSA, DENSE, TC, VT - REAL GMASS, GX, RE, W, Y + REAL*8 :: D + REAL :: DENSA, DENSE, TC, VT + REAL :: GMASS, GX, RE, W, Y REAL, PARAMETER :: PI = 3.141592654, G = 9.78956 - REAL ANU + REAL :: ANU + REAL :: W2, W3 + + !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ !Mass of stone in kg GMASS = (DENSE * PI * (D**3.)) / 6. @@ -1325,12 +1236,15 @@ SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) !THE BEST NUMBER IF (GX.LT.550) THEN W=LOG10(GX) - Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0) + W2 = W*W + Y= -1.7095 + 1.33438*W - 0.11591*(W2) RE=10**Y VT=ANU*RE/(D*DENSA) ELSE IF (GX.GE.550.AND.GX.LT.1800) THEN W=LOG10(GX) - Y= -1.81391 + 1.34671*W - 0.12427*(W**2.0) + 0.0063*(W**3.0) + W2 = W*W + W3 = W2*W + Y= -1.81391 + 1.34671*W - 0.12427*(W2) + 0.0063*(W3) RE=10**Y VT=ANU*RE/(D*DENSA) ELSE IF (GX.GE.1800.AND.GX.LT.3.45E08) THEN @@ -1344,20 +1258,20 @@ SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC) END SUBROUTINE TERMINL SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY - !!! BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD - !!! AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT, - !!! AND IF THE STONE IS IN WET OR DRY GROWTH REGIME - !!! - !!! INPUT: PC fraction of updraft water that is frozen - !!! TS temperature of hailstone (K) - !!! TC temperature of updraft air (K) - !!! ITYPE wet (2) or dry (1) growth regime - !!! OUTPUT: DELRW difference in sat vap. dens. between hail and air - !!! (kg/m3) - !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY + !!! BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD + !!! AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT, + !!! AND IF THE STONE IS IN WET OR DRY GROWTH REGIME + !!! + !!! INPUT: PC fraction of updraft water that is frozen + !!! TS temperature of hailstone (K) + !!! TC temperature of updraft air (K) + !!! ITYPE wet (2) or dry (1) growth regime + !!! OUTPUT: DELRW difference in sat vap. dens. between hail and air + !!! (kg/m3) + !!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE REAL DELRW, PC, TS, TC @@ -1392,22 +1306,28 @@ END SUBROUTINE VAPORCLOSE SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,ITYPE,DELRW) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! CALC THE STONE'S INCREASE IN MASS - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!! CALC THE STONE'S INCREASE IN MASS + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE - REAL*8 D - REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DI,ANU,RE,AE, & - TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,DELRW - INTEGER ITYPE + REAL*8 :: D, D2, D3 + REAL :: GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DI,ANU,RE,AE, & + TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,DELRW + INTEGER :: ITYPE + !local variables - REAL PI, D0, GMW2, GMI2, EW, EI,DGMV - REAL DENSEL, DENSELI, DENSELW - REAL DC !MEAN CLOUD DROPLET DIAMETER (MICRONS, 1E-6M) - REAL VOLL, VOLT !VOLUME OF NEW LAYER, TOTAL (M3) - REAL VOL1, DGMW_NOSOAK, SOAK, SOAKM - REAL DENSAC, E, AFACTOR, NS, TSCELSIUS, VIMP, W + REAL :: PI, D0, GMW2, GMI2, EW, EI,DGMV + REAL :: DENSEL, DENSELI, DENSELW + REAL :: DC, DC2 !MEAN CLOUD DROPLET DIAMETER (MICRONS, 1E-6M) + REAL :: VOLL, VOLT !VOLUME OF NEW LAYER, TOTAL (M3) + REAL :: VOL1, DGMW_NOSOAK, SOAK, SOAKM + REAL :: DENSAC, E, AFACTOR, NS, TSCELSIUS, VIMP, W + REAL :: AFACTOR2, AFACTOR3 + REAL :: W2, W3, W4 + + !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ + PI=3.141592654 !!! CALCULATE THE DIFFUSIVITY DI (m2/s) @@ -1453,10 +1373,11 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& AE=(0.57+9.0E-6*RE)*E ENDIF - !!! CALC HAILSTONE'S MASS (GM), MASS OF WATER (GMW) AND THE !!! MASS OF ICE IN THE STONE (GMI) - GM=PI/6.*(D**3.)*DENSE + D2 = D*D + D3 = D2*D + GM=PI/6.*D3*DENSE GMW=FW*GM GMI=GM-GMW @@ -1468,12 +1389,12 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !!! CALCULATE INCREASE IN MASS DUE INTERCEPTED CLD WATER, USE !!! ORIGINAL DIAMETER - GMW2=GMW+SEKDEL*(PI/4.*D**2.*VT*XW*EW) + GMW2=GMW+SEKDEL*(PI/4.*D2*VT*XW*EW) DGMW=GMW2-GMW GMW=GMW2 !!! CALCULATE THE INCREASE IN MASS DUE INTERCEPTED CLOUD ICE - GMI2=GMI+SEKDEL*(PI/4.*D**2.*VT*XI*EI) + GMI2=GMI+SEKDEL*(PI/4.*D2*VT*XI*EI) DGMI=GMI2-GMI GMI=GMI2 @@ -1494,20 +1415,25 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& IF ((DGMW.GT.0).OR.(DGMV.GT.0)) THEN !MEAN CLOUD DROPLET RADIUS, ASSUME CLOUD DROPLET CONC OF 3E8 M-3 (300 CM-3) DC = (0.74*XW / (3.14159*1000.*3.E8))**0.33333333 * 1.E6 !MICRONS + DC2 = DC*DC !!! FIND THE STOKES NUMBER (rasmussen heymsfield 1985) - NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50) !need hail radius in cm + !NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50) !need hail radius in cm + NS = 2*VT*100.*DC2*1.E-8 / (9*ANU*D*50) !need hail radius in cm !!! FIND IMPACT VELOCITY (rasmussen heymsfield 1985) IF (NS.LT.0.1)THEN W=-1. ELSE W = LOG10(NS) ENDIF + W2 = W*W + W3 = W2*W + W4 = W3*W IF (RE.GT.200) THEN IF (NS.LT.0.1) THEN VIMP = 0.0 ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN - VIMP = (0.356 + 0.4738*W - 0.1233*W**2. & - -0.1618*W**3. + 0.0807*W**4.)*VT + VIMP = (0.356 + 0.4738*W - 0.1233*W2 & + -0.1618*W3 + 0.0807*W4)*VT ELSEIF (NS.GT.10) THEN VIMP = 0.63*VT ENDIF @@ -1515,8 +1441,8 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& IF (NS.LT.0.1) THEN VIMP = 0.0 ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN - VIMP = (0.3272 + 0.4907*W - 0.09452*W**2. & - -0.1906*W**3. + 0.07105*W**4.)*VT + VIMP = (0.3272 + 0.4907*W - 0.09452*W2 & + -0.1906*W3 + 0.07105*W4)*VT ELSEIF (NS.GT.10) THEN VIMP = 0.61*VT ENDIF @@ -1524,8 +1450,8 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& IF (NS.LT.0.1) THEN VIMP = 0.0 ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN - VIMP = (0.2927 + 0.5085*W - 0.03453*W**2. & - -0.2184*W**3. + 0.03595*W**4.)*VT + VIMP = (0.2927 + 0.5085*W - 0.03453*W2 & + -0.2184*W3 + 0.03595*W4)*VT ELSEIF (NS.GT.10) THEN VIMP = 0.59*VT ENDIF @@ -1533,8 +1459,8 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& IF (NS.LT.0.4) THEN VIMP = 0.0 ELSEIF ((NS.GE.0.4).AND.(NS.LE.10)) THEN - VIMP = (0.1701 + 0.7246*W + 0.2257*W**2. & - -1.13*W**3. + 0.5756*W**4.)*VT + VIMP = (0.1701 + 0.7246*W + 0.2257*W2 & + -1.13*W3 + 0.5756*W4)*VT ELSEIF (NS.GT.10) THEN VIMP = 0.57*VT ENDIF @@ -1543,11 +1469,13 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& !RIME LAYER DENSITY, HEYMSFIELD AND PFLAUM 1985 FORM TSCELSIUS = TS - 273.16 AFACTOR = -DC*VIMP/TSCELSIUS + AFACTOR2 = AFACTOR*AFACTOR + AFACTOR3 = AFACTOR2*AFACTOR IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN DENSELW = 0.30*(AFACTOR)**0.44 ELSEIF (TSCELSIUS.GT.-5.) THEN DENSELW = EXP(-0.03115 - 1.7030*AFACTOR + & - 0.9116*AFACTOR**2. - 0.1224*AFACTOR**3.) + 0.9116*AFACTOR2 - 0.1224*AFACTOR3) ENDIF DENSELW = DENSELW * 1000. !KG M-3