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

Update fms_mixedmode branch #81

Merged
merged 2 commits into from
Dec 9, 2022
Merged

Update fms_mixedmode branch #81

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f04433c
calculate hourly-averaged O3 and PM2.5 (#226)
JianpingHuang-NOAA Dec 7, 2022
5e1d9cd
Merge branch 'dev/emc' into test/update_fms
binli2337 Dec 7, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
140 changes: 139 additions & 1 deletion driver/fvGFS/fv_nggps_diag.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,9 +116,12 @@ module fv_nggps_diags_mod
integer :: kend_pfnh, kend_w, kend_delz, kend_diss, kend_ps,kend_hs
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_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 :: 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 @@ -157,6 +160,7 @@ module fv_nggps_diags_mod
real, dimension(:,:),allocatable :: up2,dn2,uhmax03,uhmin03
real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01
real, dimension(:,:),allocatable :: maxvorthy1,maxvort02
real, dimension(:,:,:),allocatable :: o3_ave,pm25_ave,no_ave,no2_ave

public :: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
#ifdef use_WRTCOMP
Expand Down Expand Up @@ -434,6 +438,34 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
kstt_uhmin25 = nlevs+1; kend_uhmin25 = nlevs+1
nlevs = nlevs + 1
endif
id_o3_ave = register_diag_field ( trim(file_name), 'o3_ave',axes(1:3), Time, &
'Hourly averaged o3', 'ppbv', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_o3_ave > 0 ) then
allocate ( o3_ave(isco:ieco,jsco:jeco,npzo) )
kstt_o3_ave = nlevs+1; kend_o3_ave = nlevs+npzo
nlevs = nlevs + npzo
endif
id_no_ave = register_diag_field ( trim(file_name), 'no_ave',axes(1:3), Time, &
'Hourly averaged no', 'ppbv', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_no_ave > 0 ) then
allocate ( no_ave(isco:ieco,jsco:jeco,npzo) )
kstt_no_ave = nlevs+1; kend_no_ave = nlevs+npzo
nlevs = nlevs + npzo
endif
id_no2_ave = register_diag_field ( trim(file_name), 'no2_ave',axes(1:3), Time, &
'Hourly averaged no2', 'ppbv', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_no2_ave > 0 ) then
allocate ( no2_ave(isco:ieco,jsco:jeco,npzo) )
kstt_no2_ave = nlevs+1; kend_no2_ave = nlevs+npzo
nlevs = nlevs + npzo
endif
id_pm25_ave = register_diag_field ( trim(file_name), 'pm25_ave',axes(1:3), Time, &
'Hourly averaged pm25', 'ug/m3', missing_value=missing_value )
if( .not.Atm(n)%flagstruct%hydrostatic .and. id_pm25_ave > 0 ) then
allocate ( pm25_ave(isco:ieco,jsco:jeco,npzo) )
kstt_pm25_ave = nlevs+1; kend_pm25_ave = nlevs+npzo
nlevs = nlevs + npzo
endif
!
nz = size(atm(1)%ak)
allocate(ak(nz))
Expand Down Expand Up @@ -736,6 +768,22 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
deallocate ( srh01 )
deallocate ( srh03 )

!--- hourly averaged o3
if ( id_o3_ave > 0) then
call store_data(id_o3_ave, o3_ave, Time, kstt_o3_ave, kend_o3_ave)
endif
!--- hourly averaged no
if ( id_no_ave > 0) then
call store_data(id_no_ave, no_ave, Time, kstt_no_ave, kend_no_ave)
endif
!--- hourly averaged no2
if ( id_no2_ave > 0) then
call store_data(id_no2_ave, no2_ave, Time, kstt_no2_ave, kend_no2_ave)
endif
!--- hourly averaged pm25
if ( id_pm25_ave > 0) then
call store_data(id_pm25_ave, pm25_ave, Time, kstt_pm25_ave, kend_pm25_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 @@ -827,10 +875,12 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
type(time_type), intent(in) :: Time_step_atmos
real, intent(in):: zvir
integer :: i, j, k, n, ngc, nq, itrac
integer :: o3_idx,no_idx,no2_idx,pm25_idx
integer seconds, days, nsteps_per_reset
logical, save :: first_call=.true.
real, save :: first_time = 0.
integer, save :: kdtt = 0
integer, save :: kdtt = 0, kdtt1 = 0
real, save :: ucf1 = 1000., ucf2 = 1.
real :: avg_max_length
real,dimension(:,:,:),allocatable :: vort
n = 1
Expand Down Expand Up @@ -917,6 +967,46 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
endif
endif

if ( id_o3_ave > 0 .and. id_no_ave >0 &
.and. id_no2_ave > 0 .and. id_pm25_ave >0 ) then
o3_idx = get_tracer_index(MODEL_ATMOS, 'O3')
no_idx = get_tracer_index(MODEL_ATMOS, 'NO')
no2_idx = get_tracer_index(MODEL_ATMOS, 'NO2')
pm25_idx = get_tracer_index(MODEL_ATMOS, 'PM25_TOT')
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 k=1,npzo
do j=jsco,jeco
do i=isco,ieco
o3_ave(i,j,k)= 0.
no_ave(i,j,k)= 0.
no2_ave(i,j,k)= 0.
pm25_ave(i,j,k)= 0.
enddo
enddo
enddo
endif
call average_tracer_hy1(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
Atm(n)%q,o3_idx,o3_ave,nsteps_per_reset,ucf1)
call average_tracer_hy1(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
Atm(n)%q,no_idx,no_ave,nsteps_per_reset,ucf1)
call average_tracer_hy1(isco,ieco,jsco,jeco,isdo,iedo,jsdo,jedo,ncnsto,npzo,&
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

!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 @@ -925,6 +1015,26 @@ subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
endif

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)
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,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 store_data(id, work, Time, nstt, nend)
integer, intent(in) :: id
Expand Down Expand Up @@ -1317,6 +1427,34 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,add trac,i=',i,'output_name=',trim(output_name),' rc=',rc
enddo
!
!
if ( id_o3_ave > 0 ) then
call find_outputname(trim(file_name),'o3_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly averaged o3', 'ppbv', "time: point", &
axes(1:3), fcst_grid, kstt_o3_ave,kend_o3_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if ( id_no_ave > 0 ) then
call find_outputname(trim(file_name),'no_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly averaged no', 'ppbv', "time: point", &
axes(1:3), fcst_grid, kstt_no_ave,kend_no_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if ( id_no2_ave > 0 ) then
call find_outputname(trim(file_name),'no2_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly averaged no2', 'ppbv', "time: point", &
axes(1:3), fcst_grid, kstt_no2_ave,kend_no2_ave, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
!
if ( id_pm25_ave > 0 ) then
call find_outputname(trim(file_name),'pm25_ave',output_name)
call add_field_to_bundle(trim(output_name),'hourly averaged pm25', 'ug/m3', "time: point", &
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_ps > 0) then
call find_outputname(trim(file_name),'ps',output_name)
Expand Down