Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
66 commits
Select commit Hold shift + click to select a range
0220e3e
process aerosol analysis via atmanlupp job
Nov 6, 2025
9f29ce9
fix
Nov 7, 2025
92e3cad
fix code check
Nov 7, 2025
a70701b
fix code check
Nov 7, 2025
5aaf628
move aerosol analysis calculation from aeroanlfinal to in analcalc job
Nov 14, 2025
1efe601
link config files from gfs to gcafs
Nov 14, 2025
8a0ae3d
revert changes made to aeroanlinit and aeroanlfinal
Nov 14, 2025
c17a87a
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Nov 14, 2025
0a75cf1
fix failing checks
Nov 14, 2025
b95127e
add gsi info variables to default NO to match the new features in glo…
Nov 17, 2025
78616d0
set resource in anal task for gcdas
Nov 17, 2025
1b91a09
update data dependency in atmanlupp job
Nov 17, 2025
fd1d710
archive master file and grib2 files containing aerosol analysis
Nov 17, 2025
1c76767
set time to 0 in analysis file
Nov 17, 2025
a9a7ce6
discard useless changes
Nov 17, 2025
aca0300
archive if do_aero_da
Nov 18, 2025
7d6222f
link aeroinc if run gcdas
Nov 18, 2025
1b7f10e
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Nov 18, 2025
ddfd4ee
Update jobs/JGLOBAL_ATMOS_ANALYSIS_CALC
ypwang19 Nov 19, 2025
0652121
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Nov 25, 2025
6dc633f
enable atmanlprod task in GCDAS
Nov 25, 2025
d61cda1
resolve conflict
Nov 25, 2025
a79479b
resolve conflict
Nov 25, 2025
768f527
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Nov 25, 2025
1bb6bca
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Nov 25, 2025
c965cb3
fix python code check
Nov 25, 2025
f8d25e6
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 1, 2025
c087ff2
resolve conflict
Dec 1, 2025
ddc9856
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 2, 2025
1e05971
update seas1 in analysis file with background values
Dec 3, 2025
97db7dc
fix conflict
Dec 9, 2025
db9dde1
fix conflict
Dec 9, 2025
d26f7c3
update path of analcalc.sh
Dec 10, 2025
cf5256a
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 10, 2025
d0be549
Keep GDUMP as gdas unless the value is gcdas
ypwang19 Dec 11, 2025
09b07ce
Merge branch 'develop' into feature/aeroanlupp
bbakernoaa Dec 11, 2025
8dc1c0d
set GDUMP to gdas/gcdas if GDUMP is gfs/gcafs respectively
ypwang19 Dec 11, 2025
c184673
Merge remote-tracking branch 'refs/remotes/origin/feature/aeroanlupp'…
ypwang19 Dec 11, 2025
6ea7daa
suggested brief version
ypwang19 Dec 11, 2025
c420754
Merge branch 'develop' into feature/aeroanlupp
bbakernoaa Dec 12, 2025
3720198
discard duplicated changes
Dec 12, 2025
af3cf24
Merge branch 'develop' into feature/aeroanlupp
bbakernoaa Dec 15, 2025
4e2a6ca
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 16, 2025
a8cc512
move towards separate ush script
CoryMartin-NOAA Dec 17, 2025
4658e86
clean up and add a calcanl_gcafs.py
CoryMartin-NOAA Dec 17, 2025
6587813
clean up
CoryMartin-NOAA Dec 17, 2025
760c8e3
Merge branch 'develop' into feature/aeroanlupp
CoryMartin-NOAA Dec 18, 2025
c846c6a
pycodestyle
CoryMartin-NOAA Dec 18, 2025
9bb410e
flip increment and fix filename
CoryMartin-NOAA Dec 18, 2025
61d66cd
remove copy out
CoryMartin-NOAA Dec 19, 2025
7c0607e
Remove icmr increment
CoryMartin-NOAA Dec 19, 2025
e580dcc
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 19, 2025
fd29258
fix python norm
Dec 19, 2025
a9bef87
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 22, 2025
02932fe
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 23, 2025
bc1b8e5
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Dec 29, 2025
f540719
Merge branch 'develop' into feature/aeroanlupp
CoryMartin-NOAA Dec 31, 2025
d8332be
update diagb weight to 1.0
Dec 31, 2025
40aa7e7
temporarily disable pres_b products archiving
Dec 31, 2025
a786205
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Jan 5, 2026
0e5efde
Update dev/jobs/JGLOBAL_ATMOS_ANALYSIS_CALC
ypwang19 Jan 5, 2026
6bc10a7
Update parm/archive/master_gcdas.yaml.j2
ypwang19 Jan 5, 2026
5971619
link aeroinc file if running gcdas
ypwang19 Jan 5, 2026
d09e5e1
declare the chem template when run gcdas
ypwang19 Jan 5, 2026
83f8834
styling change
ypwang19 Jan 6, 2026
061057c
Merge branch 'develop' into feature/aeroanlupp
ypwang19 Jan 8, 2026
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
7 changes: 6 additions & 1 deletion dev/jobs/JGLOBAL_ATMOS_ANALYSIS_CALC
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ export DO_CALC_ANALYSIS=${DO_CALC_ANALYSIS:-"YES"}
GDATE=$(date --utc +%Y%m%d%H -d "${PDY} ${cyc} - ${assim_freq} hours")
export gPDY=${GDATE:0:8}
export gcyc=${GDATE:8:2}
export GDUMP="gdas"
GDUMP="${RUN/gfs/gdas}"
export GDUMP_ENS="enkf${GDUMP}"

export OPREFIX="${rCDUMP}.t${cyc}z."
Expand All @@ -31,6 +31,11 @@ YMD=${PDY} HH=${cyc} declare_from_tmpl -rx \
COMIN_ATMOS_ANALYSIS:COM_ATMOS_ANALYSIS_TMPL \
COMIN_ATMOS_RESTART:COM_ATMOS_RESTART_TMPL

if [[ "${RUN}" == "gcdas" ]]; then
YMD=${PDY} HH=${cyc} declare_from_tmpl -rx \
COMIN_CHEM_ANALYSIS:COM_CHEM_ANALYSIS_TMPL
fi

RUN=${GDUMP} YMD=${gPDY} HH=${gcyc} declare_from_tmpl -rx \
COMIN_OBS_PREV:COM_OBS_TMPL \
COMIN_ATMOS_HISTORY_PREV:COM_ATMOS_HISTORY_TMPL
Expand Down
2 changes: 1 addition & 1 deletion dev/parm/config/gcafs/config.aeroanlgenb
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ export aero_diffusion_horiz_len=300e3
export aero_diffusion_fixed_val=20.0
export npx_clim_b=97
export npy_clim_b=97
export aero_diagb_weight=0.9
export aero_diagb_weight=1.0
export aero_staticb_rescaling_factor=2.0
export aero_diagb_n_halo=4
export aero_diagb_n_neighbors=16
Expand Down
1 change: 1 addition & 0 deletions dev/parm/config/gcafs/config.anal
1 change: 1 addition & 0 deletions dev/parm/config/gcafs/config.analcalc
2 changes: 2 additions & 0 deletions dev/parm/config/gcafs/config.base.j2
Original file line number Diff line number Diff line change
Expand Up @@ -381,6 +381,8 @@ export SMOOTH_ENKF="NO"
export l4densvar=".true."
export lwrite4danl=".true."
export DO_CALC_INCREMENT="NO"
export USE_BUILD_GSINFO="NO"
export BUILD_GSINFO_DIR="${PARMgfs}/gsinfo"

# Early-cycle EnKF parameters
export NMEM_ENS_GFS="{{ NMEM_ENS_GFS }}"
Expand Down
2 changes: 1 addition & 1 deletion dev/parm/config/gcafs/config.resources
Original file line number Diff line number Diff line change
Expand Up @@ -680,7 +680,7 @@ case ${step} in


"anal")
if [[ "${RUN}" = *gdas ]]; then
if [[ "${RUN}" = *gdas || "${RUN}" = gcdas ]]; then
walltime="01:20:00"
case ${CASE} in
"C1152" | "C768")
Expand Down
1 change: 0 additions & 1 deletion dev/parm/config/gcafs/config.upp

This file was deleted.

16 changes: 16 additions & 0 deletions dev/parm/config/gcafs/config.upp
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
#! /usr/bin/env bash

########## config.upp ##########
# UPP specific

echo "BEGIN: config.upp"

# Get task specific resources
. "${EXPDIR}/config.resources" upp

export UPP_CONFIG="${PARMgfs}/post/upp_gcafs.yaml"

# No. of forecast hours to process in a single job
export NFHRS_PER_GROUP=3

echo "END: config.upp"
15 changes: 13 additions & 2 deletions dev/scripts/exglobal_atmos_analysis_calc.sh
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,11 @@ export CHGRESNCEXEC=${CHGRESNCEXEC:-"${EXECgfs}/enkf_chgres_recenter_nc.x"}
export CHGRESINCEXEC=${CHGRESINCEXEC:-"${EXECgfs}/interp_inc.x"}
export NTHREADS_CHGRES=${NTHREADS_CHGRES:-1}
CALCINCPY=${CALCINCPY:-"${USHgfs}/calcinc_gfs.py"}
CALCANLPY=${CALCANLPY:-"${USHgfs}/calcanl_gfs.py"}

if [[ "${RUN}" == "gcdas" ]]; then
CALCANLPY=${CALCANLPY:-"${USHgfs}/calcanl_gcafs.py"}
else
CALCANLPY=${CALCANLPY:-"${USHgfs}/calcanl_gfs.py"}
fi
DOGAUSFCANL=${DOGAUSFCANL-"NO"}
GAUSFCANLSH=${GAUSFCANLSH:-"${USHgfs}/gaussian_sfcanl.sh"}
export GAUSFCANLEXE=${GAUSFCANLEXE:-"${EXECgfs}/gaussian_sfcanl.x"}
Expand All @@ -78,6 +81,9 @@ ATMANL=${ATMANL:-"${COMOUT_ATMOS_ANALYSIS}/${APREFIX}analysis.atm.a006.nc"}
# Increment files
ATMINC=${ATMINC:-"${COMOUT_ATMOS_ANALYSIS}/${APREFIX}increment.atm.i006.nc"}
DTFANL=${DTFANL:-"${COMOUT_ATMOS_ANALYSIS}/${APREFIX}increment.dtf.i006.nc"}
if [[ "${RUN}" == "gcdas" ]]; then
AEROINC=${AEROINC:-"${COMIN_CHEM_ANALYSIS}/${APREFIX}aeroinc.nc"}
fi

# Set script / GSI control parameters
DOHYBVAR=${DOHYBVAR:-"NO"}
Expand Down Expand Up @@ -118,6 +124,11 @@ if [[ "${DO_CALC_ANALYSIS}" == "YES" ]]; then
# link analysis and increment files
${NLN} "${ATMANL}" siganl
${NLN} "${ATMINC}" siginc.nc

if [[ "${RUN}" == "gcdas" ]]; then
${NLN} "${AEROINC}" aeroinc.nc
fi

if [[ "${DOHYBVAR}" == "YES" && "${l4densvar}" == ".true." && "${lwrite4danl}" == ".true." ]]; then
${NLN} "${ATMA03}" siga03
${NLN} "${ATMI03}" sigi03.nc
Expand Down
3 changes: 3 additions & 0 deletions dev/workflow/applications/gcafs_cycled.py
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,7 @@ def _get_app_configs(self, run):
if options['do_aero_anl']:
configs += ['aeroanlgenb', 'aeroanlinit', 'aeroanlvar', 'aeroanlfinal']
configs += ['prep']
configs += ['analcalc']

if options['do_anlstat']:
configs += ['anlstat']
Expand Down Expand Up @@ -204,6 +205,8 @@ def get_task_names(self):
task_names[run] += ['aeroanlgenb']
task_names[run] += ['aeroanlinit', 'aeroanlvar', 'aeroanlfinal']
task_names[run] += ['prep']
task_names[run] += ['analcalc']
task_names[run] += ['atmanlupp', 'atmanlprod']

if options['do_anlstat']:
task_names[run] += ['anlstat']
Expand Down
40 changes: 31 additions & 9 deletions dev/workflow/rocoto/gcafs_tasks.py
Original file line number Diff line number Diff line change
Expand Up @@ -534,6 +534,34 @@ def aeroanlfinal(self):

return task

def analcalc(self):

deps = []
dep_dict = {'type': 'task', 'name': f'{self.run}_offlineanl'}
deps.append(rocoto.add_dependency(dep_dict))
dep_dict = {'type': 'task', 'name': f'{self.run}_sfcanl'}
deps.append(rocoto.add_dependency(dep_dict))
dep_dict = {'type': 'task', 'name': f'{self.run}_aeroanlfinal'}
deps.append(rocoto.add_dependency(dep_dict))
dependencies = rocoto.create_dependency(dep_condition='and', dep=deps)

resources = self.get_resource('analcalc')
task_name = f'{self.run}_analcalc'
task_dict = {'task_name': task_name,
'resources': resources,
'dependency': dependencies,
'envars': self.envars,
'cycledef': self.run.replace('enkf', ''),
'command': f'{self.HOMEgfs}/dev/job_cards/rocoto/analcalc.sh',
'job_name': f'{self.pslot}_{task_name}_@H',
'log': f'{self.rotdir}/logs/@Y@m@d@H/{task_name}.log',
'maxtries': '&MAXTRIES;'
}

task = rocoto.create_task(task_dict)

return task

def aerosol_init(self):
"""
Create a task for aerosol initialization.
Expand Down Expand Up @@ -848,18 +876,12 @@ def atmanlupp(self):
for key, value in postenvar_dict.items():
postenvars.append(rocoto.create_envar(name=key, value=str(value)))

atm_anl_path = self._template_to_rocoto_cycstring(self._base["COM_ATMOS_ANALYSIS_TMPL"])
deps = []
data = f'{atm_anl_path}/{self.run}[email protected]'
dep_dict = {'type': 'data', 'data': data, 'age': 120}
deps.append(rocoto.add_dependency(dep_dict))
data = f'{atm_anl_path}/{self.run}[email protected]'
dep_dict = {'type': 'data', 'data': data, 'age': 120}
deps.append(rocoto.add_dependency(dep_dict))
data = f'{atm_anl_path}/{self.run}[email protected]'
dep_dict = {'type': 'data', 'data': data, 'age': 60}
dep_dict = {'type': 'task', 'name': f'{self.run}_analcalc'}
deps.append(rocoto.add_dependency(dep_dict))

dependencies = rocoto.create_dependency(dep=deps, dep_condition='and')

resources = self.get_resource('upp')
task_name = f'{self.run}_atmanlupp'
task_dict = {'task_name': task_name,
Expand Down
23 changes: 23 additions & 0 deletions parm/archive/master_gcdas.yaml.j2
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,28 @@ datasets:
{% for itile in range(1,7) %}
- "{{ COMIN_CHEM_ANALYSIS | relpath(ROTDIR) }}/aeroinc.{{ cycle_YMD }}.{{ cycle_HH }}0000.fv_tracer.res.tile{{ itile }}.nc"
{% endfor %}
# Analysis Master GRIB2 data
- "{{ COMIN_ATMOS_MASTER | relpath(ROTDIR) }}/{{ head }}master.analysis.grib2"
- "{{ COMIN_ATMOS_MASTER | relpath(ROTDIR) }}/{{ head }}master.analysis.grib2.idx"

# Analysis GRIB2 (sub-sampled) data
- "{{ COMIN_ATMOS_GRIB_0p25 | relpath(ROTDIR) }}/{{ head }}pres_a.0p25.analysis.grib2"
- "{{ COMIN_ATMOS_GRIB_0p25 | relpath(ROTDIR) }}/{{ head }}pres_a.0p25.analysis.grib2.idx"
- "{{ COMIN_ATMOS_GRIB_0p50 | relpath(ROTDIR) }}/{{ head }}pres_a.0p50.analysis.grib2"
- "{{ COMIN_ATMOS_GRIB_0p50 | relpath(ROTDIR) }}/{{ head }}pres_a.0p50.analysis.grib2.idx"
- "{{ COMIN_ATMOS_GRIB_1p00 | relpath(ROTDIR) }}/{{ head }}pres_a.1p00.analysis.grib2"
- "{{ COMIN_ATMOS_GRIB_1p00 | relpath(ROTDIR) }}/{{ head }}pres_a.1p00.analysis.grib2.idx"
# temporarily disable pres_b products archiving.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bbakernoaa note this for when we get your previous PR un-reverted

#- "{{ COMIN_ATMOS_GRIB_0p25 | relpath(ROTDIR) }}/{{ head }}pres_b.0p25.analysis.grib2"
#- "{{ COMIN_ATMOS_GRIB_0p25 | relpath(ROTDIR) }}/{{ head }}pres_b.0p25.analysis.grib2.idx"
#- "{{ COMIN_ATMOS_GRIB_0p25 | relpath(ROTDIR) }}/{{ head }}pres_b.0p50.analysis.grib2"
#- "{{ COMIN_ATMOS_GRIB_0p25 | relpath(ROTDIR) }}/{{ head }}pres_b.0p50.analysis.grib2.idx"
#- "{{ COMIN_ATMOS_GRIB_1p00 | relpath(ROTDIR) }}/{{ head }}pres_b.1p00.analysis.grib2"
#- "{{ COMIN_ATMOS_GRIB_1p00 | relpath(ROTDIR) }}/{{ head }}pres_b.1p00.analysis.grib2.idx"

# Analysis netCDF (raw) data
- "{{ COMIN_ATMOS_ANALYSIS | relpath(ROTDIR) }}/{{ head }}analysis.atm.a006.nc"
- "{{ COMIN_ATMOS_ANALYSIS | relpath(ROTDIR) }}/{{ head }}analysis.sfc.a006.nc"
{% endif %}
{% if DO_PREP_OBS_AERO %}
- "{{ COMIN_OBS | relpath(ROTDIR) }}/{{ head }}aeroobs"
Expand Down Expand Up @@ -65,6 +87,7 @@ datasets:
- "{{ COMIN_ATMOS_RESTART | relpath(ROTDIR) }}/{{ r_prefix }}.coupler.res"
- "{{ COMIN_ATMOS_RESTART | relpath(ROTDIR) }}/{{ r_prefix }}.fv_core.res.nc"
{% endfor %}

# Archive the EXPDIR if requested
{% if archive_expdir %}
{% filter indent(width=4) %}
Expand Down
1 change: 1 addition & 0 deletions parm/post/itag.jinja
Original file line number Diff line number Diff line change
Expand Up @@ -12,5 +12,6 @@
kpo = {{ po | length }},
po = {{ po | join(', ') }},
rdaod = {{ rdaod | to_f90bool }}
nasa_on = {{ nasa_on | to_f90bool }}
/

37 changes: 37 additions & 0 deletions parm/post/upp_gcafs.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
upp:
config:
grib_version: "grib2"
ioform: "netcdfpara"
po: [1000.,975.,950.,925.,900.,875.,850.,825.,800.,775.,750.,725.,700.,675.,650.,625.,600.,575.,550.,525.,500.,475.,450.,425.,400.,375.,350.,325.,300.,275.,250.,225.,200.,175.,150.,125.,100.,70.,50.,40.,30.,20.,15.,10.,7.,5.,3.,2.,1.,0.7,0.4,0.2,0.1,0.07,0.04,0.02,0.01]
rdaod: False
nasa_on: False
fix_data:
mkdir:
- "{{ DATA }}"
copy:
- ["{{ 'g2tmpl_ROOT' | getenv }}/share/params_grib2_tbl_new", "{{ DATA }}/params_grib2_tbl_new"]
- ["{{ PARMgfs }}/post/nam_micro_lookup.dat", "{{ DATA }}/eta_micro_lookup.dat"]
- ["{{ EXECgfs }}/upp.x", "{{ DATA }}/"]
- ["{{ PARMgfs }}/post/itag.jinja", "{{ DATA }}/"]
- ["{{ PARMgfs }}/post/optics_luts_DUST_nasa.dat", "{{ DATA }}/optics_luts_DUST_nasa.dat"]
- ["{{ PARMgfs }}/post/optics_luts_NITR_nasa.dat", "{{ DATA }}/optics_luts_NITR_nasa.dat"]
- ["{{ PARMgfs }}/post/optics_luts_SALT_nasa.dat", "{{ DATA }}/optics_luts_SALT_nasa.dat"]
- ["{{ PARMgfs }}/post/optics_luts_SOOT_nasa.dat", "{{ DATA }}/optics_luts_SOOT_nasa.dat"]
- ["{{ PARMgfs }}/post/optics_luts_SUSO_nasa.dat", "{{ DATA }}/optics_luts_SUSO_nasa.dat"]
- ["{{ PARMgfs }}/post/optics_luts_WASO_nasa.dat", "{{ DATA }}/optics_luts_WASO_nasa.dat"]

analysis:
config:
rdaod: False
nasa_on: True
NET: GFS # upp doesn't work with GCAFS,set to GFS instead
data_in:
copy:
- ["{{ PARMgfs }}/post/gcafs/postxconfig-NT-gcafs.txt", "{{ DATA }}/postxconfig-NT.txt"]
- ["{{ COMIN_ATMOS_ANALYSIS }}/{{ RUN }}.t{{ current_cycle | strftime('%H') }}z.analysis.atm.a006.nc", "{{ DATA }}/{{ atmos_filename }}"]
- ["{{ COMIN_ATMOS_ANALYSIS }}/{{ RUN }}.t{{ current_cycle | strftime('%H') }}z.analysis.sfc.a006.nc", "{{ DATA }}/{{ flux_filename }}"]
data_out:
copy:
- ["{{ DATA }}/GFSPRS.GrbF00", "{{ COMOUT_ATMOS_MASTER }}/{{ RUN }}.t{{ current_cycle | strftime('%H') }}z.master.analysis.grib2"]
- ["{{ DATA }}/GFSPRS.GrbF00.idx", "{{ COMOUT_ATMOS_MASTER }}/{{ RUN }}.t{{ current_cycle | strftime('%H') }}z.master.analysis.grib2.idx"]

121 changes: 121 additions & 0 deletions ush/calcanl_gcafs.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@
#!/usr/bin/env python
# calcanl_gcafs.py
# script to run executables to produce netCDF analysis
# on GCAFS gaussian grid for downstream users
# based on calcanl_gfs.py
import os
import datetime
from wxflow import FileHandler
from netCDF4 import Dataset, num2date
import numpy as np

python2fortran_bool = {True: '.true.', False: '.false.'}


# function to calculate analysis from two increment files and background
def calcanl_gcafs(RunDir, ComOut, APrefix):
print('calcanl_gcafs beginning at: ', datetime.datetime.utcnow())
# add aerosol increments to background aerosol fields
# add meteorological increments to background meteorological fields

# analysis increment is already assumed to be
# at the GCAFS C384 equivalent Gaussian resolution

# define path variables
inc_file = os.path.join(RunDir, 'siginc.nc')
anl_file = os.path.join(RunDir, 'siganl')
ges_file = os.path.join(RunDir, 'sigf06')

# add meteorological increments to background meteorological fields
metvars = [['spfh', 'sphum'],
['tmp', 'T'],
['ugrd', 'u'],
['vgrd', 'v'],
['dpres', 'delp'],
['delz', 'delz'],
['o3mr', 'o3mr'],
['clwmr', 'liq_wat']]

with Dataset(inc_file, mode='r') as incfile, Dataset(ges_file, mode='r') as gesfile, Dataset(anl_file, mode='a') as anlfile:
# loop over meteorological variables and add increments to background
for ioname, incname in metvars:
print(f"Adding increment to background for variable: {ioname}")
bkg = gesfile.variables[ioname][:]
increment = incfile.variables[incname + '_inc'][:]
anl = bkg + np.flip(increment, axis=1)

anlfile.variables[ioname][:] = anl[:]

# handle pressfc as a special case
print("Adding increment to background for variable: pressfc")
# read bk attribute and compute ps_inc from delp_inc
bk = gesfile.ncattrs()
if 'bk' in bk:
bk = gesfile.getncattr('bk')
else:
# try to find bk as a variable if not an attribute
bk = gesfile.variables['bk'][:]

pressfc = gesfile.variables['pressfc'][:]
delp_inc = incfile.variables['delp_inc'][:]

# compute surface pressure increment
ps_inc = delp_inc[-1] / (bk[-1] - bk[-2])

# add increment to background surface pressure
pressfc_anl = pressfc + np.flip(ps_inc, axis=0)
anlfile.variables['pressfc'][:] = pressfc_anl[:]

# add aerosol increments to background aerosol fields
aerovars = [['so4', 'mass_fraction_of_sulfate_in_air'],
['bc1', 'mass_fraction_of_hydrophobic_black_carbon_in_air'],
['bc2', 'mass_fraction_of_hydrophilic_black_carbon_in_air'],
['oc1', 'mass_fraction_of_hydrophobic_organic_carbon_in_air'],
['oc2', 'mass_fraction_of_hydrophilic_organic_carbon_in_air'],
['dust1', 'mass_fraction_of_dust001_in_air'],
['dust2', 'mass_fraction_of_dust002_in_air'],
['dust3', 'mass_fraction_of_dust003_in_air'],
['dust4', 'mass_fraction_of_dust004_in_air'],
['dust5', 'mass_fraction_of_dust005_in_air'],
['seas1', ''], # no seas1 increment
['seas2', 'mass_fraction_of_sea_salt001_in_air'],
['seas3', 'mass_fraction_of_sea_salt002_in_air'],
['seas4', 'mass_fraction_of_sea_salt003_in_air'],
['seas5', 'mass_fraction_of_sea_salt004_in_air']
]

inc_file = os.path.join(RunDir, 'aeroinc.nc')
with Dataset(inc_file, mode='r') as incfile, Dataset(ges_file, mode='r') as gesfile, Dataset(anl_file, mode='a') as anlfile:
for ioname, incname in aerovars:
print(f"Adding increment to background for variable: {ioname}")
bkg = gesfile.variables[ioname][:]
# no seas1 increment
if ioname == 'seas1':
anl = bkg
else:
increment = incfile.variables[incname][:]
# reordering the dimensions of increment (latitude, longitude, levels) to macth background (time, levs, lat, lon)
increment_reshape = np.transpose(increment, (2, 0, 1))
anl = bkg + increment_reshape[np.newaxis, :, :, :]

anlfile.variables[ioname][:] = anl[:]
# update time (from 6 to 0) and time units in anlfile so UPP can create anl variables
time = gesfile.variables['time']
time_val = time[:]
time_units = time.units
time_calendar = getattr(time, "calendar", "standard")
cycle_time = num2date(time_val, units=time_units, calendar=time_calendar)
time_units_new = f"hours since {cycle_time[0]}"
anlfile.variables['time'][:] = 0.0
anlfile.variables['time'].setncattr("units", time_units_new)

print('calcanl_gcafs successfully completed at: ', datetime.datetime.utcnow())


# run the function if this script is called from the command line
if __name__ == '__main__':
ComOut = os.getenv('COMOUT_ATMOS_ANALYSIS', './')
APrefix = os.getenv('APREFIX', '')
RunDir = os.getenv('DATA', './')

calcanl_gcafs(RunDir, ComOut, APrefix)
Loading