Skip to content

Commit

Permalink
Merge remote-tracking branch 'emc/develop' into CADS_for_Andrew
Browse files Browse the repository at this point in the history
  • Loading branch information
@wx20jjung
wx20jjung committed Oct 3, 2023
2 parents 67014cf + fae4bbf commit c88c1c6
Show file tree
Hide file tree
Showing 12 changed files with 309 additions and 77 deletions.
40 changes: 29 additions & 11 deletions src/gsi/chemmod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,21 +40,23 @@ module chemmod
public :: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3

! fv3smoke
public :: naero_smoke_fv3,aeronames_smoke_fv3,pm2_5_innov_threshold
public :: naero_smoke_fv3,aeronames_smoke_fv3
public :: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold
public :: pm10_innov_threshold,pm10_urban_innov_threshold,pm10_bg_threshold,pm10_obs_threshold

public :: naero_gocart_wrf,aeronames_gocart_wrf

public :: pm2_5_guess,init_pm2_5_guess,&
aerotot_guess,init_aerotot_guess
public :: init_chem
public :: berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,maginnov_chem,magoberr_chem,oneob_type_chem,conconeobs
public :: berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,oneobtest_chem,maginnov_chem,magoberr_chem,oneob_type_chem,conconeobs
public :: oblat_chem,oblon_chem,obpres_chem,diag_incr,oneobschem
public :: site_scale,nsites
public :: tunable_error
public :: in_fname,out_fname,incr_fname,maxstr
public :: code_pm25_ncbufr,code_pm25_anowbufr
public :: code_pm10_ncbufr,code_pm10_anowbufr

public :: anowbufr_ext
public :: l_aoderr_table

public :: laeroana_gocart,laeroana_fv3cmaq,laeroana_fv3smoke,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, &
Expand All @@ -79,7 +81,8 @@ module chemmod
integer(i_kind) :: icvt_cmaq_fv3
real(r_kind) :: raod_radius_mean_scale,raod_radius_std_scale
real(r_kind) :: ppmv_conv = 96.06_r_kind/28.964_r_kind*1.0e+3_r_kind
real(r_kind) :: pm2_5_innov_threshold
real(r_kind) :: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold
real(r_kind) :: pm10_innov_threshold,pm10_urban_innov_threshold,pm10_bg_threshold,pm10_obs_threshold

logical :: wrf_pm2_5

Expand All @@ -90,7 +93,9 @@ module chemmod

logical :: aero_ratios

logical :: oneobtest_chem,diag_incr,berror_chem,berror_fv3_cmaq_regional
logical :: oneobtest_chem,diag_incr,berror_chem
logical :: berror_fv3_cmaq_regional,berror_fv3_sd_regional
logical :: anowbufr_ext
character(len=max_varname_length) :: oneob_type_chem
integer(i_kind), parameter :: maxstr=256
real(r_kind) :: maginnov_chem,magoberr_chem,conconeobs,&
Expand All @@ -103,7 +108,7 @@ module chemmod
real(r_kind),parameter :: pm2_5_teom_max=900.0_r_kind !ug/m3
!some parameters need to be put here since convinfo file won't
!accomodate, stands for maximum realistic value of surface pm2.5
real(r_kind),parameter :: pm10_teom_max=150.0_r_kind !ug/m3
real(r_kind),parameter :: pm10_teom_max=3000.0_r_kind !ug/m3


real(r_kind),parameter :: elev_missing=-9999.0_r_kind
Expand Down Expand Up @@ -157,10 +162,10 @@ module chemmod
'AOLGAJ', 'AISO1J', 'AISO2J', 'AISO3J', 'ATRP1J', 'ATRP2J',&
'ASQTJ', 'AOLGBJ', 'AORGCJ']
! fv3smoke
integer(i_kind), parameter :: naero_smoke_fv3=2
integer(i_kind), parameter :: naero_smoke_fv3=3

character(len=max_varname_length), dimension(naero_smoke_fv3), parameter :: &
aeronames_smoke_fv3=[character(len=max_varname_length) :: 'smoke','dust' ]
aeronames_smoke_fv3=[character(len=max_varname_length) :: 'smoke','dust','coarsepm']

! FV3CMAQ
integer(i_kind), parameter :: naero_cmaq_fv3=70 ! !number of cmaq aerosol species aero6
Expand Down Expand Up @@ -286,8 +291,15 @@ subroutine init_chem
!initialiazes default values to &CHEM namelist parameters

berror_chem=.false.
berror_fv3_cmaq_regional=.false. ! Set .true. to use berror for fv3_cmaq_regional, whose cv has 10 characters
berror_fv3_cmaq_regional=.false. ! .False. : Dont perform aerosal DA for the online RRFS_CMAQ model so dont need to read in B for RRFS_CMAQ.
! .true. : Use berror for fv3_cmaq_regional, whose cv has 10 characters
berror_fv3_sd_regional=.false. ! .False. : Dont perform aerosal DA for the RRFS_SD model so dont need to read in B for RRFS_SD.
! .true. to use berror for rrfs_sd model, whose cv has 10 characters
oneobtest_chem=.false.
anowbufr_ext=.false. ! .False. : use default anowbufr data
! .True. : use the extented bufr data
! that includes PM10, station elevation
! etal in addition to pm2.5.
maginnov_chem=30_r_kind
magoberr_chem=2_r_kind
oneob_type_chem='pm2_5'
Expand All @@ -307,9 +319,15 @@ subroutine init_chem
laeroana_gocart = .false.
laeroana_fv3cmaq = .false. ! .true. for performing aerosol analysis for regional FV3-CMAQ model(Please other parameters requred in gsimod.F90)
laeroana_fv3smoke = .false.
pm2_5_innov_threshold = 20.0_r_kind
pm2_5_innov_threshold = 15.0_r_kind
pm2_5_urban_innov_threshold = 30.0_r_kind
pm2_5_bg_threshold = 2.0_r_kind
pm10_innov_threshold = 15.0_r_kind
pm10_urban_innov_threshold = 30.0_r_kind
pm10_bg_threshold = 2.0_r_kind
pm10_obs_threshold = 140.0_r_kind ! Barry's manuscript
l_aoderr_table = .false.
icvt_cmaq_fv3 = 1 ! 1. Control variable is individual aerosol specie; 2: CV is total mass per I,J,K mode
icvt_cmaq_fv3 = 1 ! 1: Control variable is individual aerosol specie; 2: CV is total mass per I,J,K mode
raod_radius_mean_scale = 1.0_r_kind ! Tune radius of particles when calculating AOD using CRTM
raod_radius_std_scale = 1.0_r_kind ! Tune standard deviation of particles when calculating AOD using CRTM with CMAQ LUTs.
aod_qa_limit = 3
Expand Down
78 changes: 71 additions & 7 deletions src/gsi/gsi_rfv3io_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@ module gsi_rfv3io_mod
! used as background for surface observation operator
! 2022-04-15 Wang - add IO for regional FV3-CMAQ (RRFS-CMAQ) model
! 2022-08-10 Wang - add IO for regional FV3-SMOKE (RRFS-SMOKE) model
! 2023-07-30 Zhao - add IO for the analysis of the significant wave height
! (SWH, aka howv in GSI) in fv3-lam based DA (eg., RRFS-3DRTMA)
!
! subroutines included:
! sub gsi_rfv3io_get_grid_specs
Expand Down Expand Up @@ -56,6 +58,7 @@ module gsi_rfv3io_mod
use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b
use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq
use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke
use rapidrefresh_cldsurf_mod, only: i_howv_3dda

implicit none
public type_fv3regfilenameg
Expand Down Expand Up @@ -133,7 +136,7 @@ module gsi_rfv3io_mod
public :: mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql
public :: mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w
public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc
public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m
public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv
public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g
public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv
public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv
Expand All @@ -144,7 +147,7 @@ module gsi_rfv3io_mod
integer(i_kind) mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w

integer(i_kind) k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc
integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m
integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv
parameter( &
k_f10m =1, & !fact10
k_stype=2, & !soil_type
Expand All @@ -159,7 +162,8 @@ module gsi_rfv3io_mod
k_t2m =11, & ! 2 m T
k_q2m =12, & ! 2 m Q
k_orog =13, & !terrain
n2d=13 )
k_howv =14, & ! significant wave height (aka howv in GSI)
n2d=14 )
logical :: grid_reverse_flag
character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv
! copy of cvars3d excluding uv 3-d fields
Expand Down Expand Up @@ -767,6 +771,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
! 2022-04-01 Y. Wang and X. Wang - add capability to read reflectivity
! for direct radar EnVar DA using reflectivity as state
! variable, poc: [email protected]
! 2023-07-30 Zhao - added code to read significant wave height (howv) field
! from the 2D fv3-lam firstguess file (fv3_sfcdata).
! attributes:
! language: f90
! machine: ibm RS/6000 SP
Expand Down Expand Up @@ -816,6 +822,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL()
real(r_kind),dimension(:,:),pointer::ges_t2m=>NULL()
real(r_kind),dimension(:,:),pointer::ges_q2m=>NULL()
real(r_kind),dimension(:,:),pointer::ges_howv=>NULL()

real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL()
Expand Down Expand Up @@ -1093,6 +1100,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
if(mype == 0) write(6,*)'the metvarname ',trim(vartem),' will be dealt separately'
else if(trim(vartem)=='t2m') then
else if(trim(vartem)=='q2m') then
else if(trim(vartem)=='howv') then
else
write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop'
call stop2(333)
Expand All @@ -1110,7 +1118,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
do i=1,size(name_metvars2d)
vartem=trim(name_metvars2d(i))
if(.not.( (trim(vartem)=='ps'.and.fv3sar_bg_opt==0).or.(trim(vartem)=="z") &
.or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m"))) then !z is treated separately
.or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m") &
.or.(trim(vartem)=="howv"))) then ! z is treated separately
if (ifindstrloc(vardynvars,trim(vartem)) > 0) then
jdynvar=jdynvar+1
fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem)
Expand Down Expand Up @@ -1365,6 +1374,13 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it),'t2m',ges_t2m,istatus );ier=ier+istatus
if (ier/=0) call die(trim(myname),'cannot get pointers for t2m,ier=',ier)
endif

!--- significant wave height (howv)
if ( i_howv_3dda == 1 ) then
call GSI_BundleGetPointer(GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus ); ier=ier+istatus
if (ier/=0) call die(trim(myname),'cannot get pointers for howv, ier=',ier)
endif

if(mype == 0 ) then
call check(nf90_open(fv3filenamegin(it)%dynvars,nf90_nowrite,loc_id))
call check(nf90_inquire(loc_id,formatNum=ncfmt))
Expand Down Expand Up @@ -1546,7 +1562,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
endif


call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m)
call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv)

if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then
! Convert 2m guess mixing ratio to specific humidity
Expand Down Expand Up @@ -1782,7 +1798,7 @@ end subroutine gsi_bundlegetpointer_fv3lam_tracerchem_nouv

end subroutine read_fv3_netcdf_guess

subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
!$$$ subprogram documentation block
! . . . .
! subprogram: gsi_fv3ncdf2d_read
Expand All @@ -1792,6 +1808,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
! Scatter the field to each PE
! program history log:
! 2023-02-14 Hu - Bug fix for read in subdomain surface restart files
! 2023-07-30 Zhao - added IO to read significant wave height (howv) from 2D FV3-LAM
! firstguess file (fv3_sfcdata)
!
! input argument list:
! it - time index for 2d fields
Expand All @@ -1805,23 +1823,28 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype
use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype,mpi_itype
use mpeu_util, only: die
use guess_grids, only: fact10,soil_type,veg_frac,veg_type,sfc_rough, &
sfct,sno,soil_temp,soil_moi,isli
use gridmod, only: lat2,lon2,itotsub,ijn_s
use general_commvars_mod, only: ltosi_s,ltosj_s
use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr
use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension
use netcdf, only: nf90_inquire_variable
use netcdf, only: nf90_inq_varid
use netcdf, only: nf90_noerr
use mod_fv3_lola, only: fv3_h_to_ll,nxa,nya
use constants, only: grav
use constants, only: zero

implicit none

integer(i_kind),intent(in) :: it
real(r_kind),intent(in),dimension(:,:),pointer::ges_z
real(r_kind),intent(in),dimension(:,:),pointer::ges_t2m
real(r_kind),intent(in),dimension(:,:),pointer::ges_q2m
real(r_kind),intent(in),dimension(:,:),pointer::ges_howv
type (type_fv3regfilenameg),intent(in) :: fv3filenamegin
character(len=max_varname_length) :: name
integer(i_kind),allocatable,dimension(:):: dim
Expand All @@ -1835,6 +1858,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
integer(i_kind) kk,n,ns,j,ii,jj,mm1
character(len=:),allocatable :: sfcdata !='fv3_sfcdata'
character(len=:),allocatable :: dynvars !='fv3_dynvars'
! for checking the existence of howv in firstguess file
integer(i_kind) id_howv
integer(i_kind) iret_bcast

! for io_layout > 1
real(r_kind),allocatable,dimension(:,:):: sfc_fulldomain
Expand All @@ -1850,6 +1876,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
allocate(work(itotsub*n2d))
allocate( sfcn2d(lat2,lon2,n2d))

!-- initialisation of the array for howv
sfcn2d(:,:,k_howv) = zero

if(mype==mype_2d ) then
allocate(sfc_fulldomain(nx,ny))

Expand Down Expand Up @@ -1877,6 +1906,24 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
iret=nf90_inquire_dimension(gfile_loc,k,name,len)
dim(k)=len
enddo

!--- check the existence of significant wave height (howv) in 2D FV3-LAM firstguess file
! if howv is set in anavinfo (as i_howv_3dda=1), then check its existence in firstguess,
! but if it is not found in firstguess, then stop GSI run and set i_howv_3dda = 0.
if ( i_howv_3dda == 1 ) then
iret = nf90_inq_varid(gfile_loc,'howv',id_howv)
if ( iret /= nf90_noerr ) then
iret = nf90_inq_varid(gfile_loc,'HOWV',id_howv) ! double check with name in uppercase
end if
if ( iret /= nf90_noerr ) then
i_howv_3dda = 0 ! howv does not exist in firstguess, then stop GSI run.
call die('gsi_fv3ncdf2d_read','Warning: CANNOT find howv in firstguess, aborting..., iret = ', iret)
else
write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,I4,1x,A,1x,I4.4,A)') 'gsi_fv3ncdf2d_read:: Found howv in firstguess ', &
trim(sfcdata), ', iret, varid = ',iret, id_howv,' (on pe: ', mype,').'
end if
end if

!!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!!
do i=ndimensions+1,nvariables
iret=nf90_inquire_variable(gfile_loc,i,name,len)
Expand Down Expand Up @@ -1904,6 +1951,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
k=k_t2m
else if( trim(name)=='Q2M'.or.trim(name)=='q2m' ) then
k=k_q2m
else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then
k=k_howv
else
cycle
endif
Expand Down Expand Up @@ -2036,6 +2085,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain)
endif ! mype

!-- broadcast the updated i_howv_3dda to all tasks (!!!!)
call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)

!!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call mpi_scatterv(work,ijns2d,displss2d,mpi_rtype,&
Expand All @@ -2058,6 +2109,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m)
ges_t2m(:,:)=sfcn2d(:,:,k_t2m)
ges_q2m(:,:)=sfcn2d(:,:,k_q2m)
endif
if ( i_howv_3dda == 1 ) then
ges_howv(:,:)=sfcn2d(:,:,k_howv)
endif
deallocate (sfcn2d,a)
return
end subroutine gsi_fv3ncdf2d_read
Expand Down Expand Up @@ -3192,6 +3246,8 @@ subroutine wrfv3_netcdf(fv3filenamegin)
! 2019-11-22 CAPS(C. Tong) - modify "add_saved" to properly output analyses
! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da
! 2022-04-01 Y. Wang and X. Wang - add code for updating reflectivity
! 2023-07-30 Zhao - added code for the output of the analysis of
! significant wave height (howv)
!
! input argument list:
!
Expand Down Expand Up @@ -3234,6 +3290,7 @@ subroutine wrfv3_netcdf(fv3filenamegin)
real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL()
real(r_kind),pointer,dimension(:,: ):: ges_t2m =>NULL()
real(r_kind),pointer,dimension(:,: ):: ges_q2m =>NULL()
real(r_kind),pointer,dimension(:,: ):: ges_howv =>NULL()

integer(i_kind) i,k

Expand Down Expand Up @@ -3350,6 +3407,9 @@ subroutine wrfv3_netcdf(fv3filenamegin)
call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q2m',ges_q2m,istatus); ier=ier+istatus
call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'t2m',ges_t2m,istatus );ier=ier+istatus
endif
if ( i_howv_3dda == 1 ) then
call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus); ier=ier+istatus
endif
if (ier/=0) call die('wrfv3_netcdf','cannot get pointers for fv3 met-fields, ier =',ier)

if (laeroana_fv3cmaq) then
Expand Down Expand Up @@ -3559,6 +3619,10 @@ subroutine wrfv3_netcdf(fv3filenamegin)
call gsi_fv3ncdf_write_sfc(fv3filenamegin,'t2m',ges_t2m,add_saved)
call gsi_fv3ncdf_write_sfc(fv3filenamegin,'q2m',ges_q2m,add_saved)
endif
!-- output analysis of howv
if ( i_howv_3dda == 1 ) then
call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved)
endif

if(allocated(g_prsi)) deallocate(g_prsi)

Expand Down
Loading

0 comments on commit c88c1c6

Please sign in to comment.