Skip to content

Commit

Permalink
Merge branch 'develop' into feature/post-proc-incr
Browse files Browse the repository at this point in the history
  • Loading branch information
guillaumevernieres committed Aug 3, 2023
2 parents b8f27b2 + 12c3d4f commit 3fe0222
Show file tree
Hide file tree
Showing 6 changed files with 171 additions and 156 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ if(BUILD_GDASBUNDLE)

# External (required) observation operators
option("BUNDLE_SKIP_CRTM" "Don't build CRTM" "OFF") # Don't build crtm unless user passes -DBUNDLE_SKIP_CRTM=OFF
ecbuild_bundle( PROJECT crtm GIT "https://github.com/ADCollard/crtm.git" TAG v2.3-jedi.3_fix )
ecbuild_bundle( PROJECT crtm GIT "https://github.com/JCSDA/crtm.git" TAG v2.4.1-jedi.1 )

# Build GSI-B
option(BUILD_GSIBEC "Build GSI-B" OFF)
Expand Down
5 changes: 3 additions & 2 deletions build.sh
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ usage() {
# Defaults:
INSTALL_PREFIX=""
CMAKE_OPTS=""
BUILD_TARGET="${MACHINE_ID}"
BUILD_TARGET="${MACHINE_ID:-'localhost'}"
BUILD_VERBOSE="NO"
CLONE_JCSDADATA="NO"
CLEAN_BUILD="NO"
Expand Down Expand Up @@ -76,7 +76,6 @@ case ${BUILD_TARGET} in
module use $dir_root/modulefiles
module load GDAS/$BUILD_TARGET
CMAKE_OPTS+=" -DMPIEXEC_EXECUTABLE=$MPIEXEC_EXEC -DMPIEXEC_NUMPROC_FLAG=$MPIEXEC_NPROC -DBUILD_GSIBEC=ON"
CMAKE_OPTS+=" -DCLONE_JCSDADATA=$CLONE_JCSDADATA"
module list
;;
$(hostname))
Expand All @@ -87,6 +86,8 @@ case ${BUILD_TARGET} in
;;
esac

CMAKE_OPTS+=" -DCLONE_JCSDADATA=$CLONE_JCSDADATA"

BUILD_DIR=${BUILD_DIR:-$dir_root/build}
if [[ $CLEAN_BUILD == 'YES' ]]; then
[[ -d ${BUILD_DIR} ]] && rm -rf ${BUILD_DIR}
Expand Down
242 changes: 102 additions & 140 deletions scripts/exgdas_global_marine_analysis_vrfy.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,189 +22,147 @@
import gen_eva_obs_yaml
import marine_eva_post
import diag_statistics
from soca_vrfy import plot_config, plot_horizontal_slice, plot_zonal_slice
from soca_vrfy import statePlotter, plotConfig
import subprocess
from datetime import datetime, timedelta

comout = os.getenv('COM_OCEAN_ANALYSIS')
com_ice_history = os.getenv('COM_ICE_HISTORY_PREV')
com_ocean_history = os.getenv('COM_OCEAN_HISTORY_PREV')
data = os.getenv('DATA')
pdy = os.getenv('PDY')
cyc = os.getenv('cyc')
RUN = os.getenv('RUN')
bcyc = str((int(cyc) - 3) % 24).zfill(2)
gcyc = str((int(cyc) - 6) % 24).zfill(2)
gcdate = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') - timedelta(hours=int(os.getenv('assim_freq')))

bcyc = str((int(cyc) - 3) % 24).zfill(2)
grid_file = os.path.join(comout, f'{RUN}.t'+bcyc+'z.ocngrid.nc')

# for eva
diagdir = os.path.join(comout, 'diags')
HOMEgfs = os.getenv('HOMEgfs')


#######################################
# INCREMENT
# ocean increment
#######################################
incr_cmap = 'RdBu'

data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocninc.nc')
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap=incr_cmap,
comout=os.path.join(comout, 'vrfy', 'incr'))
config = plotConfig(grid_file=grid_file,
data_file=data_file,
lats=np.arange(-60, 60, 10),
variables_zonal=['Temp', 'Salt'],
variables_horiz=['Temp', 'Salt', 'ave_ssh'],
allbounds={'Temp': [-0.5, 0.5],
'Salt': [-0.1, 0.1],
'ave_ssh': [-0.1, 0.1]},
colormap='RdBu',
comout=os.path.join(comout, 'vrfy', 'incr'))
ocnIncPlotter = statePlotter(config)
ocnIncPlotter.plot()

#######################################
# zonal slices

for lat in np.arange(-60, 60, 10):

for max_depth in [700.0, 5000.0]:
config['lats'] = [lat]
config['max depth'] = max_depth

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-.5, .5]})
plot_zonal_slice(config)

# Salinity
config.update({'variable': 'Salt', 'levels': [1], 'bounds': [-.1, .1]})
plot_zonal_slice(config)

# sea ice increment
#######################################
# Horizontal slices

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-1, 1]})
plot_horizontal_slice(config)

# Salinity
config.update({'variable': 'Salt', 'bounds': [-0.1, 0.1]})
plot_horizontal_slice(config)

# Sea surface height
config.update({'variable': 'ave_ssh', 'bounds': [-0.1, 0.1]})
plot_horizontal_slice(config)

#######################################
# Sea ice
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ice.incr.nc')
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap=incr_cmap,
comout=os.path.join(comout, 'vrfy', 'incr'))

for proj in ['North', 'South']:
# concentration
config.update({'variable': 'aicen', 'bounds': [-0.2, 0.2], 'proj': proj})
plot_horizontal_slice(config)

# thickness
config.update({'variable': 'hicen', 'bounds': [-0.5, 0.5], 'proj': proj})
plot_horizontal_slice(config)
config = plotConfig(grid_file=grid_file,
data_file=data_file,
lats=np.arange(-60, 60, 10),
variables_horiz=['aicen', 'hicen', 'hsnon'],
allbounds={'aicen': [-0.2, 0.2],
'hicen': [-0.5, 0.5],
'hsnon': [-0.1, 0.1]},
colormap='RdBu',
projs=['North', 'South'],
comout=os.path.join(comout, 'vrfy', 'incr'))
iceIncPlotter = statePlotter(config)
iceIncPlotter.plot()

# snow depth
config.update({'variable': 'hsnon', 'bounds': [-0.1, 0.1], 'proj': proj})
plot_horizontal_slice(config)

#######################################
# Analysis/Background
#######################################

# sea ice analysis
#######################################
# Sea ice
data_files = [os.path.join(comout, f'{RUN}.t{cyc}z.iceana.nc'),
os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc')]
dirs_out = ['ana', 'bkg']
ice_vars = {'bkg': ['aice_h', 'hs_h', 'hi_h'], 'ana': ['aicen', 'hicen', 'hsnon']}
for data_file, dir_out in zip(data_files, dirs_out):
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap='jet',
comout=os.path.join(comout, 'vrfy', dir_out))

for proj in ['North', 'South', 'Global']:
# concentration
var = ice_vars[dir_out]
config.update({'variable': var[0], 'bounds': [0.0, 1.0], 'proj': proj})
plot_horizontal_slice(config)

# thickness
config.update({'variable': var[1], 'bounds': [0.0, 4.0], 'proj': proj})
plot_horizontal_slice(config)

# snow depth
config.update({'variable': var[2], 'bounds': [0.0, 0.5], 'proj': proj})
plot_horizontal_slice(config)
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.iceana.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['aicen', 'hicen', 'hsnon'],
allbounds={'aicen': [0.0, 1.0],
'hicen': [0.0, 4.0],
'hsnon': [0.0, 0.5]},
colormap='jet',
projs=['North', 'South', 'Global'],
comout=os.path.join(comout, 'vrfy', 'ana'))
iceAnaPlotter = statePlotter(config)
iceAnaPlotter.plot()

#######################################
# Ocean surface
data_files = [os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc'),
os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc')]
dirs_out = ['ana', 'bkg']
ocn_vars = ['ave_ssh', 'Temp', 'Salt']
for data_file, dir_out in zip(data_files, dirs_out):
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap='jet',
comout=os.path.join(comout, 'vrfy', dir_out))

# ssh
config.update({'variable': 'ave_ssh', 'bounds': [-1.8, 1.3], 'proj': proj, 'levels': [1]})
plot_horizontal_slice(config)

# sst
config.update({'variable': 'Temp', 'bounds': [-1.8, 34.0], 'proj': proj, 'levels': [1]})
plot_horizontal_slice(config)
# sea ice background
#######################################

# sss
config.update({'variable': 'Salt', 'bounds': [30, 38], 'proj': proj, 'levels': [1]})
plot_horizontal_slice(config)
data_file = os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['aice_h', 'hs_h', 'hi_h'],
allbounds={'aice_h': [0.0, 1.0],
'hs_h': [0.0, 4.0],
'hi_h': [0.0, 0.5]},
colormap='jet',
projs=['North', 'South', 'Global'],
comout=os.path.join(comout, 'vrfy', 'bkg'))
iceBkgPlotter = statePlotter(config)
iceBkgPlotter.plot()

#######################################
# Std Bkg. Error
# ocean surface analysis
#######################################
bmat_cmap = 'jet'
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc')
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap=bmat_cmap,
comout=os.path.join(comout, 'vrfy', 'bkgerr'))

#######################################
# zonal slices

for lat in np.arange(-60, 60, 10):

for max_depth in [700.0, 5000.0]:
config['lats'] = [lat]
config['max depth'] = max_depth

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 1.5]})
plot_zonal_slice(config)

# Salinity
config.update({'variable': 'Salt', 'levels': [1], 'bounds': [0, .2]})
plot_zonal_slice(config)
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['ave_ssh', 'Temp', 'Salt'],
allbounds={'ave_ssh': [-1.8, 1.3],
'Temp': [-1.8, 34.0],
'Salt': [30, 38]},
colormap='jet',
comout=os.path.join(comout, 'vrfy', 'ana'))
ocnAnaPlotter = statePlotter(config)
ocnAnaPlotter.plot()

#######################################
# Horizontal slices
# ocean surface background
#######################################

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 2]})
plot_horizontal_slice(config)
data_file = os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['ave_ssh', 'Temp', 'Salt'],
allbounds={'ave_ssh': [-1.8, 1.3],
'Temp': [-1.8, 34.0],
'Salt': [30, 38]},
colormap='jet',
comout=os.path.join(comout, 'vrfy', 'bkg'))
ocnBkgPlotter = statePlotter(config)
ocnBkgPlotter.plot()

# Salinity
config.update({'variable': 'Salt', 'bounds': [0, 0.2]})
plot_horizontal_slice(config)
#######################################
# background error
#######################################

# Sea surface height
config.update({'variable': 'ave_ssh', 'bounds': [0, 0.1]})
plot_horizontal_slice(config)
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
lats=np.arange(-60, 60, 10),
variables_zonal=['Temp', 'Salt'],
variables_horiz=['Temp', 'Salt', 'ave_ssh'],
allbounds={'Temp': [0, 2],
'Salt': [0, 0.2],
'ave_ssh': [0, 0.1]},
colormap='jet',
comout=os.path.join(comout, 'vrfy', 'bkgerr'))
bkgErrPlotter = statePlotter(config)
bkgErrPlotter.plot()

#######################################
# eva plots
#######################################

evadir = os.path.join(HOMEgfs, 'sorc', f'{RUN}.cd', 'ush', 'eva')
marinetemplate = os.path.join(evadir, 'marine_gdas_plots.yaml')
Expand All @@ -231,4 +189,8 @@
print('running eva on', infile)
subprocess.run(['eva', infile], check=True)

#######################################
# calculate diag statistics
#######################################

diag_statistics.get_diag_stats()
File renamed without changes.
Loading

0 comments on commit 3fe0222

Please sign in to comment.