Skip to content

Commit

Permalink
Add support for aerosol analysis fields from Gaussian increment (#50)
Browse files Browse the repository at this point in the history
  • Loading branch information
@CoryMartin-NOAA
CoryMartin-NOAA authored Aug 12, 2024
1 parent 9382fd0 commit 4d4daa2
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 19 deletions.
102 changes: 88 additions & 14 deletions src/netcdf_io/calc_analysis.fd/inc2anl.f90
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! module inc2anl
!! module inc2anl
!! contains subroutines for calculating analysis fields
!! for a given input background and increment
!! for a given input background and increment
!! Original: 2019-09-18 martin - original module
!! 2019-10-24 martin - removed support for NEMSIO background but
!! allows for either NEMSIO or netCDF analysis write
!! 2020-01-21 martin - parallel IO support added
!! 2024-08-12 martin - aerosol support added
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module inc2anl
module inc2anl
implicit none

private
Expand All @@ -17,7 +18,10 @@ module inc2anl
integer, parameter :: nincv=10
character(len=7) :: incvars_nemsio(nincv), incvars_netcdf(nincv), incvars_ncio(nincv)
integer, parameter :: nnciov=20
character(len=7) :: iovars_netcdf(nnciov)
integer, parameter :: naero=14
integer, parameter :: naero_copy=6
character(len=7) :: iovars_netcdf(nnciov), iovars_aero(naero), copyvars_aero(naero_copy)
character(len=50) :: incvars_aero(naero)

data incvars_nemsio / 'ugrd ', 'vgrd ', 'dpres ', 'delz ', 'o3mr ',&
'tmp ', 'spfh ', 'clwmr ', 'icmr ', 'pres '/
Expand All @@ -29,6 +33,20 @@ module inc2anl
'delz ', 'dpres ', 'dzdt ', 'grle ', 'hgtsfc ',&
'icmr ', 'o3mr ', 'pressfc', 'rwmr ', 'snmr ',&
'spfh ', 'tmp ', 'ugrd ', 'vgrd ', 'cld_amt'/
data iovars_aero / 'so4 ', 'bc1 ', 'bc2 ', 'oc1 ', 'oc2 ', &
'dust1 ', 'dust2 ', 'dust3 ', 'dust4 ', 'dust5 ',&
'seas1 ', 'seas2 ', 'seas3 ', 'seas4 '/
data incvars_aero / 'mass_fraction_of_sulfate_in_air',&
'mass_fraction_of_hydrophobic_black_carbon_in_air', &
'mass_fraction_of_hydrophilic_black_carbon_in_air', &
'mass_fraction_of_hydrophobic_organic_carbon_in_air', &
'mass_fraction_of_hydrophilic_organic_carbon_in_air', &
'mass_fraction_of_dust001_in_air', 'mass_fraction_of_dust002_in_air', &
'mass_fraction_of_dust003_in_air', 'mass_fraction_of_dust004_in_air', &
'mass_fraction_of_dust005_in_air', 'mass_fraction_of_sea_salt001_in_air', &
'mass_fraction_of_sea_salt002_in_air', 'mass_fraction_of_sea_salt003_in_air', &
'mass_fraction_of_sea_salt004_in_air' /
data copyvars_aero / 'seas5 ', 'so2 ', 'dms ', 'msa ', 'pm25 ', 'pm10 ' /

contains
subroutine gen_anl
Expand All @@ -38,7 +56,7 @@ subroutine gen_anl
! increment, add the two together, and write out
! the analysis to a new file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use vars_calc_analysis, only: mype
use vars_calc_analysis, only: mype, do_aero
implicit none
! variables local to this subroutine
integer :: i, j, iincvar
Expand Down Expand Up @@ -68,11 +86,24 @@ subroutine gen_anl
end if
else
! otherwise just write out what is in the input to the output
if (mype==0) print *, 'Copying from Background ', iovars_netcdf(i)
if (mype==0) print *, 'Copying from Background ', iovars_netcdf(i)
call copy_ges_to_anl(iovars_netcdf(i))
end if
end do

! if including aerosols, loop and add them
if (do_aero) then
do i=1,naero
if (mype==0) print *, 'Adding Increment to ', iovars_aero(i)
call add_aero_inc(iovars_aero(i), incvars_aero(i))
end do
! need to handle fields that just need copied
do i=1,naero_copy
if (mype==0) print *, 'Copying from Background ', copyvars_aero(i)
call copy_ges_to_anl(copyvars_aero(i))
end do
end if

end subroutine gen_anl

subroutine close_files
Expand Down Expand Up @@ -157,17 +188,60 @@ subroutine copy_ges_to_anl(varname)
end do
end if
end select
else
else
if (mype == 0) write(6,*) varname, 'not in background file, skipping...'
end if

end subroutine copy_ges_to_anl

subroutine add_aero_inc(fcstvar, incvar)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! subroutine add_aero_inc
! generic subroutine for adding increment to background
! and writing out to analysis for aerosol variables
! args:
! fcstvar - input string of netCDF fcst/anal var name
! incvar - input string of netCDF increment var name
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use vars_calc_analysis, only: fcstncfile, anlncfile, aero_file,&
nlat, nlon, nlev, anlfile, levpe, mype
use module_ncio, only: Dataset, read_vardata, write_vardata, &
open_dataset, close_dataset, has_var
use nemsio_module
implicit none
! input vars
character(7), intent(in) :: fcstvar
character(50), intent(in) :: incvar
! local variables
real, allocatable, dimension(:,:,:) :: work3d_bg
real, allocatable, dimension(:,:) :: work3d_inc
real, allocatable, dimension(:) :: work1d
integer :: j,jj,k,krev,iret
type(Dataset) :: incncfile
do k=1,nlev
if (mype == levpe(k)) then
! get first guess
call read_vardata(fcstncfile, trim(fcstvar), work3d_bg, nslice=k, slicedim=3)
! get increment
incncfile = open_dataset(aero_file)
call read_vardata(incncfile, trim(incvar), work3d_inc, nslice=k, slicedim=1)
! add increment to background
work3d_bg(:,:,1) = work3d_bg(:,:,1) + work3d_inc(:,:)
! write out analysis to file
call write_vardata(anlncfile, trim(fcstvar), work3d_bg, nslice=k, slicedim=3)
end if
end do
! clean up and close
deallocate(work3d_bg, work3d_inc)
call close_dataset(incncfile)

end subroutine add_aero_inc

subroutine add_increment(fcstvar, incvar)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! subroutine add_increment
! generic subroutine for adding increment to background
! and writing out to analysis
! and writing out to analysis
! args:
! fcstvar - input string of netCDF fcst/anal var name
! incvar - input string of netCDF increment var name
Expand All @@ -189,7 +263,7 @@ subroutine add_increment(fcstvar, incvar)
real, allocatable, dimension(:) :: work1d
integer :: j,jj,k,krev,iret
type(Dataset) :: incncfile

if (has_var(fcstncfile, fcstvar)) then
do k=1,nlev
if (mype == levpe(k)) then
Expand Down Expand Up @@ -234,17 +308,17 @@ subroutine add_increment(fcstvar, incvar)
end if
deallocate(work3d_bg)
call close_dataset(incncfile)
else
else
write(6,*) fcstvar, ' not in background file, skipping...'
end if

end subroutine add_increment

subroutine add_psfc_increment
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! subroutine add_psfc_increment
! special case of getting surface pressure analysis from
! bk5 and delp increment to get sfc pressure increment
! bk5 and delp increment to get sfc pressure increment
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use vars_calc_analysis, only: fcstncfile, anlncfile, nlat, nlon, incr_file,&
use_nemsio_anl, anlfile, nlev
Expand All @@ -260,7 +334,7 @@ subroutine add_psfc_increment
type(Dataset) :: incncfile

! get bk5 from attributes
call read_attribute(fcstncfile, 'bk', bk5)
call read_attribute(fcstncfile, 'bk', bk5)
! read in delp increment to get ps increment
incncfile = open_dataset(incr_file)
call read_vardata(incncfile, 'delp_inc', work3d_inc)
Expand All @@ -277,7 +351,7 @@ subroutine add_psfc_increment
! write out to file
if (use_nemsio_anl) then
allocate(work1d(nlon*nlat))
! now write out new psfc to NEMSIO analysis file
! now write out new psfc to NEMSIO analysis file
work1d = reshape(work2d,(/size(work1d)/))
call nemsio_writerecv(anlfile, 'pres', 'sfc', 1, work1d, iret=iret)
if (iret /=0) write(6,*) 'Error with NEMSIO write sfc pressure'
Expand Down
12 changes: 10 additions & 2 deletions src/netcdf_io/calc_analysis.fd/init_calc_analysis.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
!! 2019-09-25 martin - update to allow for netCDF I/O
!! 2019-10-24 martin - update to support nemsio output
!! 2020-01-17 martin - parallel IO support added
!! 2024-04-04 martin - aerosol support added
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module init_calc_analysis
implicit none
Expand All @@ -14,20 +15,22 @@ subroutine read_nml
!! read in namelist parameters from
!! calc_analysis.nml file in same directory
!! as executable
use vars_calc_analysis, only: anal_file, fcst_file, incr_file, use_nemsio_anl, fhr, mype, npes, jedi
use vars_calc_analysis, only: anal_file, fcst_file, incr_file, aero_file, use_nemsio_anl, do_aero, fhr, mype, npes, jedi
implicit none
! local variables to this subroutine
character(len=500) :: datapath = './'
character(len=500) :: analysis_filename = 'atmanl.nc'
character(len=500) :: firstguess_filename = 'atmges.nc'
character(len=500) :: increment_filename = 'atminc.nc'
character(len=500) :: aero_inc_filename = 'aeroinc.nc'
character(len=2) :: hrstr
integer, parameter :: lunit = 10
logical :: lexist = .false.
namelist /setup/ datapath, analysis_filename, firstguess_filename, increment_filename, fhr, use_nemsio_anl, jedi
namelist /setup/ datapath, analysis_filename, firstguess_filename, increment_filename, aero_inc_filename, fhr, use_nemsio_anl, do_aero, jedi

fhr = 6 ! default to 6 hour cycle only
use_nemsio_anl = .false. ! default to using netCDF for background and analysis
do_aero = .false. ! default to only do atm, not aero also
jedi = .false. ! default to GSI (not JEDI)

! read in the namelist
Expand All @@ -47,16 +50,21 @@ subroutine read_nml
anal_file = trim(adjustl(datapath)) // '/' // trim(adjustl(analysis_filename)) // '.' // hrstr
fcst_file = trim(adjustl(datapath)) // '/' // trim(adjustl(firstguess_filename)) // '.' // hrstr
incr_file = trim(adjustl(datapath)) // '/' // trim(adjustl(increment_filename)) // '.' // hrstr
aero_file = trim(adjustl(datapath)) // '/' // trim(adjustl(aero_inc_filename)) // '.' // hrstr
else
anal_file = trim(adjustl(datapath)) // '/' // trim(adjustl(analysis_filename))
fcst_file = trim(adjustl(datapath)) // '/' // trim(adjustl(firstguess_filename))
incr_file = trim(adjustl(datapath)) // '/' // trim(adjustl(increment_filename))
aero_file = trim(adjustl(datapath)) // '/' // trim(adjustl(aero_inc_filename))
endif

if (mype == 0) then
write(6,*) 'Analysis File = ', trim(anal_file)
write(6,*) 'First Guess File = ', trim(fcst_file)
write(6,*) 'Increment File = ', trim(incr_file)
if (do_aero) then
write(6,*) 'Aerosol Increment File = ', trim(aero_file)
end if
if (fhr > 0) write(6,*) 'Forecast Hour = ', fhr
write(6,*) 'Number of PEs = ', npes
write(6,*) 'input guess file and increment file should be in netCDF format'
Expand Down
8 changes: 5 additions & 3 deletions src/netcdf_io/calc_analysis.fd/vars_calc_analysis.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,14 +6,15 @@
!! 2019-09-26 martin - add support for netCDF read/write
!! 2019-10-24 martin - support NEMSIO output write
!! 2020-01-17 martin - parallel IO support added
!! 2024-04-04 martin - aerosol support added
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module vars_calc_analysis
use nemsio_module, only: nemsio_gfile
use module_ncio, only: Dataset
implicit none
private

public :: anal_file, fcst_file, incr_file
public :: anal_file, fcst_file, incr_file, aero_file
public :: idate, jdate
public :: idate6, jdate6
public :: nfday, nfhour, nfminute, nfsecondn, nfsecondd
Expand All @@ -23,14 +24,15 @@ module vars_calc_analysis
public :: work1
public :: nhrs_assim
public :: use_nemsio_anl
public :: do_aero
public :: fcstncfile, anlncfile, incncfile
public :: fhrs_pe
public :: fhr
public :: mype, npes
public :: levpe
public :: jedi

character(len=500) :: anal_file, fcst_file, incr_file
character(len=500) :: anal_file, fcst_file, incr_file, aero_file
integer, dimension(7) :: idate, jdate
integer, dimension(6) :: idate6, jdate6
integer :: nfday, nfhour, nfminute, nfsecondn, nfsecondd
Expand All @@ -40,7 +42,7 @@ module vars_calc_analysis
type(nemsio_gfile) :: anlfile
real, allocatable, dimension(:) :: work1
integer :: nhrs_assim, fhr
logical :: use_nemsio_anl
logical :: use_nemsio_anl, do_aero
type(Dataset) :: fcstncfile, anlncfile, incncfile
integer, dimension(7) :: fhrs_pe
integer :: mype, npes
Expand Down

0 comments on commit 4d4daa2

Please sign in to comment.