Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

Fourth reconciliation PR from produciton/RRFS.v1 #363

Open
wants to merge 2 commits into
base: dev/emc
Choose a base branch
from
Open

Fourth reconciliation PR from produciton/RRFS.v1 #363

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a254cfe
"to output the average of smoke and dust"
haiqinli Mar 13, 2024
f28be0e
Merge branch 'dev/emc' into rrfsv1-to-ufs/dev4
grantfirl Nov 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
121 changes: 117 additions & 4 deletions driver/UFS/fv_nggps_diag.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,11 +117,15 @@ module fv_nggps_diags_mod
integer :: kstt_dbz, kend_dbz, kstt_omga, kend_omga
integer :: kstt_windvect, kend_windvect
integer :: kstt_o3_ave, kend_o3_ave, kstt_pm25_ave, kend_pm25_ave
integer :: kstt_smoke_ave, kend_smoke_ave
integer :: kstt_dust_ave, kend_dust_ave
integer :: kstt_coarsepm_ave, kend_coarsepm_ave
integer :: kstt_no_ave, kend_no_ave, kstt_no2_ave, kend_no2_ave
integer :: id_wmaxup,id_wmaxdn,kstt_wup, kend_wup,kstt_wdn,kend_wdn
integer :: id_uhmax03,id_uhmin03,id_uhmax25,id_uhmin25,id_maxvort01
integer :: id_maxvorthy1,kstt_maxvorthy1,kstt_maxvort01,id_ustm
integer :: id_o3_ave,id_pm25_ave,id_no_ave,id_no2_ave
integer :: id_smoke_ave, id_dust_ave, id_coarsepm_ave
integer :: kend_maxvorthy1,kend_maxvort01,id_vstm,id_srh01,id_srh03
integer :: kstt_uhmax03,kstt_uhmin03,kend_uhmax03,kend_uhmin03
integer :: kstt_uhmax25,kstt_uhmin25,kend_uhmax25,kend_uhmin25
Expand Down Expand Up @@ -161,6 +165,7 @@ module fv_nggps_diags_mod
real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01
real, dimension(:,:),allocatable :: maxvorthy1,maxvort02
real, dimension(:,:,:),allocatable :: o3_ave,pm25_ave,no_ave,no2_ave
real, dimension(:,:),allocatable :: smoke_ave,dust_ave,coarsepm_ave

public :: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
#ifdef use_WRTCOMP
Expand Down Expand Up @@ -466,6 +471,27 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
kstt_pm25_ave = nlevs+1; kend_pm25_ave = nlevs+npzo
nlevs = nlevs + npzo
endif
id_smoke_ave = register_diag_field ( trim(file_name), 'smoke_ave',axes(1:2), Time, &
'Hourly averaged surface smoke', 'ug/kg', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_smoke_ave > 0 ) then
allocate ( smoke_ave(isco:ieco,jsco:jeco) )
kstt_smoke_ave = nlevs+1; kend_smoke_ave = nlevs+1
nlevs = nlevs + npzo
endif
id_dust_ave = register_diag_field ( trim(file_name), 'dust_ave',axes(1:2), Time, &
'Hourly averaged surface dust', 'ug/kg', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_dust_ave > 0 ) then
allocate ( dust_ave(isco:ieco,jsco:jeco) )
kstt_dust_ave = nlevs+1; kend_dust_ave = nlevs+1
nlevs = nlevs + npzo
endif
id_coarsepm_ave = register_diag_field ( trim(file_name), 'coarsepm_ave',axes(1:2), Time, &
'Hourly averaged surface coarsepm', 'ug/kg', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_coarsepm_ave > 0 ) then
allocate ( coarsepm_ave(isco:ieco,jsco:jeco) )
kstt_coarsepm_ave = nlevs+1; kend_coarsepm_ave = nlevs+1
nlevs = nlevs + npzo
endif
!
nz = size(atm(1)%ak)
allocate(ak(nz))
Expand Down Expand Up @@ -784,6 +810,18 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
if ( id_pm25_ave > 0) then
call store_data(id_pm25_ave, pm25_ave, Time, kstt_pm25_ave, kend_pm25_ave)
endif
!--- hourly averaged surface smoke
if ( id_smoke_ave > 0) then
call store_data(id_smoke_ave, smoke_ave, Time, kstt_smoke_ave, kend_smoke_ave)
endif
!--- hourly averaged surface dust
if ( id_dust_ave > 0) then
call store_data(id_dust_ave, dust_ave, Time, kstt_dust_ave, kend_dust_ave)
endif
!--- hourly averaged surface coarsepm
if ( id_coarsepm_ave > 0) then
call store_data(id_coarsepm_ave, coarsepm_ave, Time, kstt_coarsepm_ave, kend_coarsepm_ave)
endif
!--- max hourly 0-1km vert. vorticity
if ( id_maxvort01 > 0) then
call store_data(id_maxvort01, maxvort01, Time, kstt_maxvort01, kend_maxvort01)
Expand Down Expand Up @@ -876,6 +914,7 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
real, intent(in):: zvir
integer :: i, j, k, n, ngc, nq, itrac
integer :: o3_idx,no_idx,no2_idx,pm25_idx
integer :: smoke_idx,dust_idx,coarsepm_idx
integer seconds, days, nsteps_per_reset
logical, save :: first_call=.true.
real, save :: first_time = 0.
Expand Down Expand Up @@ -1000,13 +1039,47 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
Atm(n)%q,no2_idx,no2_ave,nsteps_per_reset,ucf1)
call average_tracer_hy1(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
Atm(n)%q,pm25_idx,pm25_ave,nsteps_per_reset,ucf2)

kdtt1=kdtt1+1
!else
! print *,'calculating hourly-averaegtd o3 or pm25'
! call mpp_error(FATAL, 'Missing hourly-averaged o3 or pm25 in diag_table')
! stop
endif
if ( id_smoke_ave > 0 .and. id_dust_ave >0 &
.and. id_coarsepm_ave > 0 ) then
smoke_idx = get_tracer_index(MODEL_ATMOS, 'smoke')
dust_idx = get_tracer_index(MODEL_ATMOS, 'dust')
coarsepm_idx = get_tracer_index(MODEL_ATMOS, 'coarsepm')
if (first_call) then
call get_time(Time_step_atmos, seconds, days)
first_time=seconds
first_call=.false.
kdtt1=0
endif
nsteps_per_reset = nint(avg_max_length/first_time)
if(mod(kdtt1,nsteps_per_reset)==0)then
do j=jsco,jeco
do i=isco,ieco
smoke_ave(i,j)= 0.
dust_ave(i,j)= 0.
coarsepm_ave(i,j)= 0.
enddo
enddo
endif
call average_tracer_hy2(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
Atm(n)%q,smoke_idx,smoke_ave,nsteps_per_reset,ucf2)
call average_tracer_hy2(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
Atm(n)%q,dust_idx,dust_ave,nsteps_per_reset,ucf2)
call average_tracer_hy2(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
Atm(n)%q,coarsepm_idx,coarsepm_ave,nsteps_per_reset,ucf2)

kdtt1=kdtt1+1
!else
! print *,'calculating hourly-averaegtd o3 or pm25'
! call mpp_error(FATAL, 'Missing hourly-averaged o3 or pm25 in diag_table')
! stop
endif
!allocate hailcast met field arrays
if (do_hailcast) then
call hailcast_compute(Atm(n),sphum,liq_wat,ice_wat,rainwat,snowwat,graupel, zvir, &
Expand All @@ -1017,23 +1090,42 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
end subroutine fv_nggps_tavg
!
subroutine average_tracer_hy1(is,ie,js,je,isd,ied,jsd,jed, &
ncns,npz,tracer,tr_idx,tracer_ave,nstp,unitcf)
ncns,npz,tracer,tr_idx,tracer_ave,nstp,unitcf)
integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
integer, intent(in):: ncns, npz, nstp, tr_idx
real, intent(in) :: unitcf
real, intent(in), dimension(isd:ied,jsd:jed,npz,ncns):: tracer
real, intent(in), dimension(isd:ied,jsd:jed,npz,ncns):: tracer
real, intent(inout), dimension(is:ie,js:je,npz):: tracer_ave
integer i, j, k
!
do k=1,npz
do j=js,je
do i=is,ie
tracer_ave(i,j,k)=tracer_ave(i,j,k)+tracer(i,j,k,tr_idx)/nstp*unitcf
enddo
enddo
enddo

end subroutine average_tracer_hy1

!
subroutine average_tracer_hy2(is,ie,js,je,isd,ied,jsd,jed, &
ncns,npz,tracer,tr_idx,tracer_ave,nstp,unitcf)
integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
integer, intent(in):: ncns, npz, nstp, tr_idx
real, intent(in) :: unitcf
real, intent(in), dimension(isd:ied,jsd:jed,npz,ncns):: tracer
real, intent(inout), dimension(is:ie,js:je):: tracer_ave
integer i, j
!

do j=js,je
do i=is,ie
tracer_ave(i,j)=tracer_ave(i,j)+tracer(i,j,npz,tr_idx)/nstp*unitcf
enddo
enddo
enddo

end subroutine average_tracer_hy1
end subroutine average_tracer_hy2

!
subroutine store_data(id, work, Time, nstt, nend)
Expand Down Expand Up @@ -1455,6 +1547,27 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
axes(1:3), fcst_grid, kstt_pm25_ave,kend_pm25_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if ( id_smoke_ave > 0 ) then
call find_outputname(trim(file_name),'smoke_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly averaged surface smoke', 'ug/kg', "time: point", &
axes(1:2), fcst_grid, kstt_smoke_ave,kend_smoke_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if ( id_dust_ave > 0 ) then
call find_outputname(trim(file_name),'dust_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly surface dust', 'ug/kg', "time: point", &
axes(1:2), fcst_grid, kstt_dust_ave,kend_dust_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if ( id_coarsepm_ave > 0 ) then
call find_outputname(trim(file_name),'coarsepm_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly averaged surface coarsepm', 'ug/kg', "time: point", &
axes(1:2), fcst_grid, kstt_coarsepm_ave,kend_coarsepm_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if( id_ps > 0) then
call find_outputname(trim(file_name),'ps',output_name)
Expand Down