diff --git a/scripts/da_blending_fv3.py b/scripts/da_blending_fv3.py new file mode 100755 index 000000000..79df7fbc5 --- /dev/null +++ b/scripts/da_blending_fv3.py @@ -0,0 +1,173 @@ +import numpy as np +from netCDF4 import Dataset +import raymond +import sys + +print("Starting blending code") + +Lx = float(sys.argv[1]) # BLENDING_LENGTHSCALE +pi = np.pi +nbdy = 40 # 20 on each side + +blend = str(sys.argv[5]) # TRUE: Blend RRFS and GDAS EnKF + # FALSE: Don't blend, activate cold2warm start only, and use either GDAS or RRFS +if blend == "TRUE": + blend = True + print("Blending is activated") +elif blend == "FALSE": + blend = False + print("Blending is **NOT** activated! Will perform cold2warm start conversion only.") +else: + print("variable 'blend' not set correctly") + exit() + +if not blend: + use_host_EnKF = str(sys.argv[6]) # TRUE: Final EnKF will be GDAS (no blending) + # FALSE: Final EnKF will be RRFS (no blending) + if use_host_EnKF == "TRUE": + use_host_EnKF = True + print("Using GDAS EnKF (no blending)") + elif use_host_EnKF == "FALSE": + use_host_EnKF = False + print("Using RRFS EnKF (no blending)") + else: + print("variable 'use_host_EnKF' not set correctly") + exit() + + +# List of variables from the regional (fg) and global (bg) to blend respectively. +vars_fg = ["u", "v", "T", "sphum", "delp"] +vars_bg = ["u_cold2fv3", "v_cold2fv3", "t_cold2fv3", "sphum_cold2fv3", "delp_cold2fv3"] + +# GDAS EnKF file chgres_cube-ed from gaussian grid to ESG grid. +# There is one more step to make sure the winds are on the same +# grid staggering and have the same orientation as the RRFS winds. +glb_fg = str(sys.argv[2]) +glb_fg_nc = Dataset(glb_fg) +glb_nlon = glb_fg_nc.dimensions["lon"].size # 1820 (lonp=1821) +glb_nlat = glb_fg_nc.dimensions["lat"].size # 1092 (latp=1093) +glb_nlev = glb_fg_nc.dimensions["lev"].size # 66 (levp=67) +glb_Dx = 3.0 + +# RRFS EnKF restart file fv_core.res.tile1 on ESG grid. +reg_fg = str(sys.argv[3]) +# Open the blended file for updating the required vars (use a copy of the regional file) +reg_fg_nc = Dataset(reg_fg, mode="a") +nlon = reg_fg_nc.dimensions["xaxis_1"].size # 1820 (xaxis_2=1821) +nlat = reg_fg_nc.dimensions["yaxis_2"].size # 1092 (yaxis_1=1093) +nlev = reg_fg_nc.dimensions["zaxis_1"].size # 65 +Dx = 3.0 + +# RRFS EnKF restart file fv_tracer.res.tile1 on ESG grid. +reg_fg_t = str(sys.argv[4]) +# Open the blended file for updating the required vars (use a copy of the regional file) +reg_fg_t_nc = Dataset(reg_fg_t, mode="a") + +# Check matching grids +# Note: global_hyblev_fcst_rrfsL65.txt has 0.000 0.0000000 as the 66th row, so +# don't compare glb_nlev and nlev because glb_nlev will be 66 and nlev will be 65. +# As a work around for now, we will just slice the top (or bottom?) 65 levels of +# the global file and blend those with the regional file. +if (glb_nlon != nlon or glb_nlat != nlat or glb_Dx != Dx): + print(f"glb_nlon:{glb_nlon} vs nlon:{nlon}") + print(f"glb_nlat:{glb_nlat} vs nlat:{nlat}") + print(f"glb_Dx:{glb_Dx} vs Dx:{Dx}") + print("grids don't match") + exit() + +eps = (np.tan(pi*Dx/Lx))**-6 # 131319732.431162 + +print(f"Input") +print(f" RRFS restart (core) : {reg_fg}") +print(f" RRFS restart (tracer) : {reg_fg_t}") +print(f" GDAS coldstart from chgres : {glb_fg}") +print(f" Lx : {Lx}") +print(f" Dx : {Dx}") +print(f" NLON : {nlon}") +print(f" NLAT : {nlat}") +print(f" NLEV : {nlev}") +print(f" eps : {eps}") +print(f"Output") +print(f" Blended background file : {reg_fg}, {reg_fg_t}") + +# Step 1. blend. +for (var_fg, var_bg) in zip(vars_fg, vars_bg): + i = vars_fg.index(var_fg) + + if var_fg == "sphum": + reg_nc = reg_fg_t_nc + else: + reg_nc = reg_fg_nc + glb_nc = glb_fg_nc + + dim = len(np.shape(reg_nc[var_fg]))-1 + if dim == 2: # 2D vars + glb = np.float64(glb_nc[var_bg][:, :]) # (1093 1820) + reg = np.float64(reg_nc[var_fg][:, :, :]) # (1, 1093, 1820) + ntim = np.shape(reg)[0] + nlat = np.shape(reg)[1] + nlon = np.shape(reg)[2] + nlev = 1 + glb = np.reshape(glb, [ntim, nlat, nlon]) # add time dim bc missing from chgres + var_out = np.zeros(shape=(nlon, nlat, 1), dtype=np.float64) + field = np.zeros(shape=(nlon*nlat), dtype=np.float64) + var_work = np.zeros(shape=((nlon+nbdy), (nlat+nbdy), 1), dtype=np.float64) + field_work = np.zeros(shape=((nlon+nbdy)*(nlat+nbdy)), dtype=np.float64) + if dim == 3: # 3D vars + glb = np.float64(glb_nc[var_bg][:, :, :]) + reg = np.float64(reg_nc[var_fg][:, :, :, :]) # (1, 65, 1093, 1820) + ntim = np.shape(reg)[0] + nlev = np.shape(reg)[1] + nlat = np.shape(reg)[2] + nlon = np.shape(reg)[3] + glb = np.reshape(glb, [ntim, nlev, nlat, nlon]) # add time dim bc missing from chgres + var_out = np.zeros(shape=(nlon, nlat, nlev, 1), dtype=np.float64) + field = np.zeros(shape=(nlon*nlat, nlev), dtype=np.float64) + var_work = np.zeros(shape=((nlon+nbdy), (nlat+nbdy), nlev, 1), dtype=np.float64) + field_work = np.zeros(shape=((nlon+nbdy)*(nlat+nbdy), nlev), dtype=np.float64) + glbT = np.transpose(glb) # (1820, 1093, 65) + regT = np.transpose(reg) # (1820, 1093, 65, 1) + + nlon_start = int(nbdy/2) + nlon_end = int(nlon+nbdy/2) + nlat_start = int(nbdy/2) + nlat_end = int(nlat+nbdy/2) + + if blend: + print(f"Blending backgrounds for {var_fg}/{var_bg}") + var_work[nlon_start:nlon_end, nlat_start:nlat_end, :] = glbT - regT + field_work = var_work.reshape((nlon+nbdy)*(nlat+nbdy), nlev, order="F") # order="F" (FORTRAN) + field_work = raymond.raymond(field_work, nlon+nbdy, nlat+nbdy, eps, nlev) + var_work = field_work.reshape(nlon+nbdy, nlat+nbdy, nlev, order="F") + var_out = var_work[nlon_start:nlon_end, nlat_start:nlat_end, :] + if dim == 2: # 2D vars + var_out = var_out[:, :, 0] + regT[:, :, 0] + var_out = np.reshape(var_out, [nlon, nlat, 1]) # add the time ("1") dimension back + if dim == 3: # 3D vars + var_out = var_out + regT[:, :, :, 0] + var_out = np.reshape(var_out, [nlon, nlat, nlev, 1]) # add the time ("1") dimension back + else: # skip blending and use either host enkf (GDAS) or the RRFS enkf + print(f"Blending code is NOT executing blending!") + print(f" This is used for finishing converting the cold start files into warm start format.") + if use_host_EnKF: + print(f"---> Use the GDAS EnKF") + var_out = glbT + else: + print(f"---> Use the RRFS EnKF") + var_out = regT + + var_out = np.transpose(var_out) # (1, 50, 834, 954) + + # Overwrite blended fields to blended file. + if dim == 2: # 2D vars + reg_nc.variables[var_fg][:, :, :] = var_out + if dim == 3: # 3D vars + reg_nc.variables[var_fg][:, :, :, :] = var_out + +# Close nc files +reg_nc.close() # blended file +glb_nc.close() + +print("Blending finished successfully.") + +exit(0) diff --git a/scripts/da_chgres_cold2fv3.py b/scripts/da_chgres_cold2fv3.py new file mode 100755 index 000000000..a4e194c4a --- /dev/null +++ b/scripts/da_chgres_cold2fv3.py @@ -0,0 +1,193 @@ +import numpy as np +from netCDF4 import Dataset +import remap_dwinds +import remap_scalar +import chgres_winds # might need to rename in the future +import sys + +print("Reading in NETCDF4 Files... ", end="\r") +warm = str(sys.argv[1]) +cold = str(sys.argv[2]) +grid = str(sys.argv[3]) +akbk = str(sys.argv[4]) +akbkcold = str(sys.argv[5]) +orog = str(sys.argv[6]) + +warmnc = Dataset(warm) +coldnc = Dataset(cold, mode="a") +akbknc = Dataset(akbk) +gridnc = Dataset(grid) +akbkcoldnc = Dataset(akbkcold) +orognc = Dataset(orog) +print("Reading in NETCDF4 Files... Done.") + +ColdStartWinds = True +VertRemapScalar = True +VertRemapWinds = True +WriteData = True + +# STEP 1. ROTATE THE WINDS FROM CHGRES +if ColdStartWinds: + print("Starting ColdStartWinds.... ", end="\r") + # Data from warm restarts + u = np.float64(warmnc["u"][0, :, :, :]) + v = np.float64(warmnc["v"][0, :, :, :]) + nlev = np.shape(u)[0] # 127,z + nlat = np.shape(u)[1] # 769,y + nlon = np.shape(u)[2] # 768,x + + km = nlev + nlev = coldnc.createDimension("nlev", nlev) # 127 + + # Data from cold chgres + u_s = np.float64(coldnc["u_s"][:, :, :]) # (128, 769, 768) + v_s = np.float64(coldnc["v_s"][:, :, :]) # (128, 769, 768) + u_w = np.float64(coldnc["u_w"][:, :, :]) # (128, 768, 769) + v_w = np.float64(coldnc["v_w"][:, :, :]) # (128, 768, 769) + + # grid data - usually has 2x as many grid points, so we need every other value. + gridx = np.float64(gridnc["x"][0:-1:2, 0:-1:2]) + gridy = np.float64(gridnc["y"][0:-1:2, 0:-1:2]) + + # Fortran wants everything transposed and in fortran array type + gridx = np.asfortranarray(gridx.transpose()) + gridy = np.asfortranarray(gridy.transpose()) + u_s = np.asfortranarray(u_s.transpose()) + v_s = np.asfortranarray(v_s.transpose()) + u_w = np.asfortranarray(u_w.transpose()) + v_w = np.asfortranarray(v_w.transpose()) + + # Initialize some computed fields to zero + ud = np.float64(0.0*u_s) # initialize to zero + vd = np.float64(0.0*u_w) # initialize to zero + + # rotate winds to model d-grid (~30s; nodes=2; cpus=128) + chgres_winds.main(gridx, gridy, u_s, v_s, u_w, v_w, ud, vd) + + print("Starting ColdStartWinds.... Done.") + +# STEP 2. VERTICAL REMAPPING OF SCALARS +if VertRemapScalar: + print("Starting VertRemapScalar... ", end="\r") + + # Data from cold restarts + ak0 = np.float64(akbkcoldnc["vcoord"][0, :]) # ( lev, ) == (128, ) + bk0 = np.float64(akbkcoldnc["vcoord"][1, :]) # ( lev, ) == (128, ) + ak = ak0[1:] + bk = bk0[1:] + ps = np.float64(coldnc["ps"][:, :]) # ( lat, lon) == ( 768, 768) + zh = np.float64(coldnc["zh"][:, :, :]) # (levp, lat, lon) == (129, 768, 768) + omga = np.float64(coldnc["w"][:, :, :]) # ( lev, lat, lon) == (128, 768, 768) + delp_cold = np.float64(coldnc["delp"][:, :, :]) # (127, 768, 768) + t_cold = np.float64(coldnc["t"][:, :, :]) # (128, 768, 768) + Atm_phis = np.float64(orognc["orog_filt"][:, :])*9.80665 + + ak0[0] = 1.000000000000000E-009 + bk0[0] = 1.000000000000000E-009 + + sphum = np.float64(coldnc["sphum"][:, :, :]) # ( lev, lat, lon) == (128, 768, 768) + liq_wat = np.float64(coldnc["liq_wat"][:, :, :]) + o3mr = np.float64(coldnc["o3mr"][:, :, :]) + ice_wat = np.float64(coldnc["ice_wat"][:, :, :]) + rainwat = np.float64(coldnc["rainwat"][:, :, :]) + snowwat = np.float64(coldnc["snowwat"][:, :, :]) + graupel = np.float64(coldnc["graupel"][:, :, :]) + ntracers = 7 + qa = np.array([sphum, liq_wat, o3mr, ice_wat, rainwat, snowwat, graupel]) + + # Fortran wants everything transposed and in fortran array type + ak = np.asfortranarray(ak) # Don't transpose 1D array + bk = np.asfortranarray(bk) # Don't transpose 1D array + ak0 = np.asfortranarray(ak0) + bk0 = np.asfortranarray(bk0) + Atm_phis = np.asfortranarray(np.transpose(Atm_phis)) + ps = np.asfortranarray(ps.transpose()) + zh = np.asfortranarray(zh.transpose()) + omga = np.asfortranarray(omga.transpose()) + qa = np.asfortranarray(qa.transpose()) + delp_cold = np.asfortranarray(delp_cold.transpose()) + t_cold = np.asfortranarray(t_cold.transpose()) + + isrt = 1 + jsrt = 1 + iend = np.shape(t_cold)[0] + jend = np.shape(t_cold)[1] + npz = np.shape(t_cold)[2]-1 + levp = npz + 1 # (km) + + # Initialize some computed fields to zero + Atm_delp = 1.0*delp_cold[:, :, 1:] # initialize to zero (delp for sfcp) + Atm_q = 1.0*qa[:, :, 1:, :] # initialize to zero (tracers... sphum=1) + Atm_pt = 1.0*t_cold[:, :, 1:] # initialize to zero (temperature) + Atm_ps = 1.0*ps[:, :] # initialize to zero (need for remap_dwinds) + + remap_scalar.main(levp, npz, ntracers, ak0, bk0, ak, bk, ps, qa, zh, omga, t_cold, + isrt, iend, jsrt, jend, Atm_pt, Atm_q, Atm_delp, Atm_phis, Atm_ps) + + print("Starting VertRemapScalar... Done.") + + +# STEP 3. VERTICAL REMAPPING OF WINDS +if VertRemapWinds: + print("Starting VertRemapWinds....", end="\r") + + Atm_u = 1.0*ud[:, :, 1:] # Atm_u has levs=npz, ud has levs km + Atm_v = 1.0*vd[:, :, 1:] + + # vertically remap the dwinds + remap_dwinds.main(levp, npz, ak0, bk0, ak, bk, ps, ud, vd, + isrt, iend, jsrt, jend, Atm_u, Atm_v, Atm_ps) + + print("Starting VertRemapWinds.... Done.") + +else: + Atm_u = ud[:, :, 1:] + Atm_v = vd[:, :, 1:] + + +# STEP 4. WRITE OUT DATA +if WriteData: + # tranpose ud, vd back to original shape, cutoff one of levels (there is an extra level), + # add a new variable to the nc file by duplicating the corresponding u/v variable and + # redefining the shape of the array, finally, assign ud/vd into the u/v variable in nc file. + + if ColdStartWinds: + # For ud + new_var = "u_cold2fv3" + ud = np.transpose(Atm_u) + var_to_duplicate = coldnc.variables["u_s"] + coldnc.createVariable(new_var, var_to_duplicate.datatype, ('nlev', 'latp', 'lon')) + coldnc.variables[new_var][:, :, :] = ud + + # For vd + new_var = "v_cold2fv3" + vd = np.transpose(Atm_v) + var_to_duplicate = coldnc.variables["v_w"] + coldnc.createVariable(new_var, var_to_duplicate.datatype, ('nlev', 'lat', 'lonp')) + coldnc.variables[new_var][:, :, :] = vd + + if VertRemapScalar: + # For Temperature + new_var = "t_cold2fv3" + t_cold = np.transpose(Atm_pt) + var_to_duplicate = coldnc.variables["t"] + coldnc.createVariable(new_var, var_to_duplicate.datatype, ('nlev', 'lat', 'lon')) + coldnc.variables[new_var][:, :, :] = t_cold + + # For delp + new_var = "delp_cold2fv3" + delp = Atm_delp.T + var_to_duplicate = coldnc.variables["delp"] + coldnc.createVariable(new_var, var_to_duplicate.datatype, ('nlev', 'lat', 'lon')) + coldnc.variables[new_var][:, :, :] = delp + + # For sphum + new_var = "sphum_cold2fv3" + sphum = Atm_q[:, :, :, 0].T + var_to_duplicate = coldnc.variables["sphum"] + coldnc.createVariable(new_var, var_to_duplicate.datatype, ('nlev', 'lat', 'lon')) + coldnc.variables[new_var][:, :, :] = sphum + +# close the nc files +warmnc.close() +coldnc.close() diff --git a/scripts/exregional_make_ics.sh b/scripts/exregional_make_ics.sh index 3a2a5f163..cf1562921 100755 --- a/scripts/exregional_make_ics.sh +++ b/scripts/exregional_make_ics.sh @@ -681,20 +681,140 @@ located in the following directory: # #----------------------------------------------------------------------- # -# Move initial condition, surface, control, and 0-th hour lateral bound- -# ary files to ICs_BCs directory. +# Run large-scale blending +# +#----------------------------------------------------------------------- +# NOTES: +# * The large-scale blending is broken down into 4 major parts +# 1) chgres_winds: This part rotates the coldstart winds from chgres to the model D-grid. +# It is based on atmos_cubed_sphere/tools/external_ic.F90#L433, and it +# is equivalent to the fv3jedi tool called ColdStartWinds. +# 2) remap_dwinds: This part vertically remaps the D-grid winds. +# It is based on atmos_cubed_sphere/tools/external_ic.F90#L3485, and it +# is part of the fv3jedi tool called VertRemap. +# 3) remap_scalar: This part vertically remaps all other variables. +# It is based on atmos_cubed_sphere/tools/external_ic.F90#L2942, and it +# is the other part of the fv3jedi tool called VertRemap. +# 4) raymond: This is the actual blending code which uses the raymond filter. The +# raymond filter is a sixth-order tangent low-pass implicit filter +# and can be controlled via the cutoff length scale (Lx). +# +# * Currently blended fields: u, v, t, dpres, and sphum +# -) Blending only works with GDASENKF (netcdf) +# +# * Two RRFS EnKF member files are needed: fv_core and fv_tracer. +# -) fv_core contains u, v, t, and dpres +# -) fv_tracer contains sphum +# +# * Before we can do any blending, the coldstart files from chgres need to be +# processed. This includes rotating the winds and vertically remapping all the +# variables. The cold start file has u_w, v_w, u_s, and v_s which correspond +# to the D-grid staggering. +# -) u_s is the D-grid south tangential wind component (m/s) +# -) v_s is the D-grid south face normal wind component (m/s) +# -) u_w is the D-grid west face tangential wind component (m/s) +# -) v_w is the D-grid west face normal wind component (m/s) +# -) https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/blob/bdeee64e860c5091da2d169b1f4307ad466eca2c/tools/external_ic.F90 +# -) https://dtcenter.org/sites/default/files/events/2020/20201105-1300p-fv3-gfdl-1.pdf +# +cdate_crnt_fhr_m1=$( date --utc --date "$yyyymmdd $hh UTC - 1 hours" "+%Y%m%d%H" ) +if [[ $DO_ENS_BLENDING == "TRUE" && + $cdate_crnt_fhr -ge ${FIRST_BLENDED_CYCLE_DATE} && + $EXTRN_MDL_NAME_ICS == "GDASENKF" ]]; then + + echo "Blending Starting." + ulimit -s unlimited + export OMP_STACKSIZE=2G + ncores=$(( NNODES_MAKE_ICS*PPN_MAKE_ICS )) # OMP_NUM_THREADS=ntasks*cpus-per-task + export OMP_NUM_THREADS=$ncores #WCOSS2:"96", Hera/Orion:"80" + export FI_OFI_RXM_SAR_LIMIT=3145728 + export FI_MR_CACHE_MAX_COUNT=0 + export MPICH_OFI_STARTUP_CONNECT=1 + + # Python/F2Py scripts + cp_vrfy $SCRIPTSDIR/da_blending_fv3.py . + cp_vrfy $SCRIPTSDIR/da_chgres_cold2fv3.py . + + # F2Py shared object files + ln_vrfy -sf $LIB64DIR/raymond.so . + ln_vrfy -sf $LIB64DIR/chgres_winds.so . + ln_vrfy -sf $LIB64DIR/remap_scalar.so . + ln_vrfy -sf $LIB64DIR/remap_dwinds.so . + + # Required NETCDF files - HOST MODEL (e.g., GDAS; these files should already be present) + #cp_vrfy out.atm.tile${TILE_RGNL}.nc . + #cp_vrfy out.sfc.tile${TILE_RGNL}.nc . + #cp_vrfy gfs_ctrl.nc . + + # Required NETCDF files - RRFS + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.fv_core.res.tile1.nc ./fv_core.res.tile1.nc + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.fv_tracer.res.tile1.nc ./fv_tracer.res.tile1.nc + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.fv_core.res.nc ./fv_core.res.nc + + # Required FIX files + cp_vrfy $FIXLAM/C3359_grid.tile7.nc . + cp_vrfy $FIXLAM/C3359_oro_data.tile7.halo0.nc . + + # Shortcut the file names + warm=./fv_core.res.tile1.nc + cold=./out.atm.tile7.nc + grid=./C3359_grid.tile7.nc + akbk=./fv_core.res.nc + akbkcold=./gfs_ctrl.nc + orog=./C3359_oro_data.tile7.halo0.nc + bndy=./gfs.bndy.nc + + # Run convert coldstart files to fv3 restart (rotate winds and remap). + ${BLENDINGPYTHON} da_chgres_cold2fv3.py $warm $cold $grid $akbk $akbkcold $orog + + # Shortcut the file names/arguments. + Lx=$ENS_BLENDING_LENGTHSCALE + glb=./out.atm.tile${TILE_RGNL}.nc + reg=./fv_core.res.tile1.nc + trcr=./fv_tracer.res.tile1.nc + + # Blend OR finish convert cold2warm start without blending. + blend=${BLEND} # TRUE: Blend RRFS and GDAS EnKF + # FALSE: Don't blend, activate cold2warm start only, and use either GDAS or RRFS + use_host_enkf=${USE_HOST_ENKF} # ignored if blend="TRUE". + # TRUE: Final EnKF will be GDAS (no blending) + # FALSE: Final EnKF will be RRFS (no blending) + ${BLENDINGPYTHON} da_blending_fv3.py $Lx $glb $reg $trcr $blend $use_host_enkf + cp_vrfy ./fv_core.res.tile1.nc ${ics_dir}/. + cp_vrfy ./fv_tracer.res.tile1.nc ${ics_dir}/. + + # Move the remaining RESTART files to INPUT + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.coupler.res ${ics_dir}/coupler.res + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.fv_core.res.nc ${ics_dir}/fv_core.res.nc + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.fv_srf_wnd.res.tile1.nc ${ics_dir}/fv_srf_wnd.res.tile1.nc + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.phy_data.nc ${ics_dir}/phy_data.nc + cp_vrfy ${NWGES_BASEDIR}/${cdate_crnt_fhr_m1}${SLASH_ENSMEM_SUBDIR}/fcst_fv3lam/RESTART/${yyyymmdd}.${hh}0000.sfc_data.nc ${ics_dir}/sfc_data.nc + cp_vrfy gfs_ctrl.nc ${ics_dir} + cp_vrfy gfs.bndy.nc ${ics_dir}/gfs_bndy.tile${TILE_RGNL}.000.nc + + +fi + # #----------------------------------------------------------------------- # -mv_vrfy out.atm.tile${TILE_RGNL}.nc \ +# Move initial condition, surface, control, and 0-th hour lateral bound- +# ary files to ICs_BCs directory. Only do this if blending is off or on- +# ly for the first DA cycle if blending is on inorder to coldstart the +# system. +#----------------------------------------------------------------------- +# +if [[ ($cdate_crnt_fhr -lt ${FIRST_BLENDED_CYCLE_DATE} && $DO_ENS_BLENDING == "TRUE") || $DO_ENS_BLENDING == "FALSE" ]]; then + mv_vrfy out.atm.tile${TILE_RGNL}.nc \ ${ics_dir}/gfs_data.tile${TILE_RGNL}.halo${NH0}.nc -mv_vrfy out.sfc.tile${TILE_RGNL}.nc \ + mv_vrfy out.sfc.tile${TILE_RGNL}.nc \ ${ics_dir}/sfc_data.tile${TILE_RGNL}.halo${NH0}.nc -mv_vrfy gfs_ctrl.nc ${ics_dir} + mv_vrfy gfs_ctrl.nc ${ics_dir} -mv_vrfy gfs.bndy.nc ${ics_dir}/gfs_bndy.tile${TILE_RGNL}.000.nc + mv_vrfy gfs.bndy.nc ${ics_dir}/gfs_bndy.tile${TILE_RGNL}.000.nc +fi # #----------------------------------------------------------------------- # @@ -704,6 +824,9 @@ mv_vrfy gfs.bndy.nc ${ics_dir}/gfs_bndy.tile${TILE_RGNL}.000.nc # cp_vrfy ${ics_dir}/*.nc ${ics_nwges_dir}/. +if [[ $DO_ENS_BLENDING == "TRUE" && $cdate_crnt_fhr -ge ${FIRST_BLENDED_CYCLE_DATE} ]]; then + cp_vrfy ${ics_dir}/coupler.res ${ics_nwges_dir}/. +fi # #----------------------------------------------------------------------- diff --git a/scripts/exregional_run_prepstart.sh b/scripts/exregional_run_prepstart.sh index 89cff5bf0..c6406bd41 100755 --- a/scripts/exregional_run_prepstart.sh +++ b/scripts/exregional_run_prepstart.sh @@ -142,6 +142,7 @@ YYYYMMDD=${YYYYMMDDHH:0:8} YYYYJJJHH=${YYYY}${JJJ}${HH} current_time=$(date "+%T") +cdate_crnt_fhr=$( date --utc --date "${YYYYMMDD} ${HH} UTC" "+%Y%m%d%H" ) YYYYMMDDm1=$(date +%Y%m%d -d "${START_DATE} 1 days ago") YYYYMMDDm2=$(date +%Y%m%d -d "${START_DATE} 2 days ago") @@ -210,6 +211,10 @@ if [ ${cycle_type} == "spinup" ]; then for cyc_start in "${CYCL_HRS_SPINSTART[@]}"; do if [ ${HH} -eq ${cyc_start} ]; then BKTYPE=1 + if [[ ${DO_ENS_BLENDING} == "TRUE" && $cdate_crnt_fhr -ge ${FIRST_BLENDED_CYCLE_DATE} ]]; then + echo "do blending" + BKTYPE=3 # warm start from blended ics + fi fi done if [ ${cycle_subtype} == "spinup" ]; then @@ -281,6 +286,45 @@ if [ ${BKTYPE} -eq 1 ] ; then # cold start, use prepare cold strat initial file print_err_msg_exit "Error: cannot find cold start initial condition from : ${bkpath}" fi +elif [[ $BKTYPE == 3 ]]; then + bkpath=${lbcs_root}/$YYYYMMDD$HH${SLASH_ENSMEM_SUBDIR}/ics + if [ -r "${bkpath}/coupler.res" ]; then + cp_vrfy ${bkpath}/fv_core.res.nc fv_core.res.nc + cp_vrfy ${bkpath}/fv_core.res.tile1.nc fv_core.res.tile1.nc + cp_vrfy ${bkpath}/fv_srf_wnd.res.tile1.nc fv_srf_wnd.res.tile1.nc + cp_vrfy ${bkpath}/fv_tracer.res.tile1.nc fv_tracer.res.tile1.nc + cp_vrfy ${bkpath}/phy_data.nc phy_data.nc + cp_vrfy ${bkpath}/sfc_data.nc sfc_data.nc + cp_vrfy ${bkpath}/gfs_ctrl.nc gfs_ctrl.nc + + ln_vrfy -s ${bkpath}/coupler.res bk_coupler.res + ln_vrfy -s ${bkpath}/fv_core.res.nc bk_fv_core.res.nc + ln_vrfy -s ${bkpath}/fv_core.res.tile1.nc bk_fv_core.res.tile1.nc + ln_vrfy -s ${bkpath}/fv_srf_wnd.res.tile1.nc bk_fv_srf_wnd.res.tile1.nc + ln_vrfy -s ${bkpath}/fv_tracer.res.tile1.nc bk_fv_tracer.res.tile1.nc + ln_vrfy -s ${bkpath}/phy_data.nc bk_phy_data.nc + ln_vrfy -s ${bkpath}/sfc_data.nc bk_sfc_data.nc + + if [ ${SAVE_CYCLE_LOG} == "TRUE" ] ; then + echo "${YYYYMMDDHH}(${cycle_type}): blended warm start at ${current_time} from $bkpath " >> ${EXPTDIR}/log.cycles + fi + else + print_err_msg_exit "Error: cannot find blended warm start initial condition from : ${bkpath}" + fi +# generate coupler.res with right date + head -1 bk_coupler.res > coupler.res + tail -1 bk_coupler.res >> coupler.res + tail -1 bk_coupler.res >> coupler.res + +# +# remove checksum from restart files. Checksum will cause trouble if model initializes from blended ics +# + filelistn="fv_core.res.nc fv_core.res.tile1.nc fv_srf_wnd.res.tile1.nc fv_tracer.res.tile1.nc phy_data.nc sfc_data.nc" + + for file in ${filelistn}; do + ncatted -a checksum,,d,, ${file} + done + ncatted -O -a source,global,c,c,'FV3GFS GAUSSIAN NETCDF FILE' fv_core.res.tile1.nc else # Setup the INPUT directory for warm start cycles, which can be spin-up cycle or product cycle. diff --git a/ush/config.sh_rrfs_a_enkf_c3_retro b/ush/config.sh_rrfs_a_enkf_c3_retro index 4335f18eb..1eb245ad5 100644 --- a/ush/config.sh_rrfs_a_enkf_c3_retro +++ b/ush/config.sh_rrfs_a_enkf_c3_retro @@ -13,6 +13,7 @@ PREDEF_GRID_NAME=RRFS_CONUS_3km DO_ENSEMBLE="TRUE" #DO_ENSFCST="TRUE" +DO_ENS_BLENDING="TRUE" # Use this for blending or, instead of do_ensinit, for cold2warm conversion. #DO_DACYCLE="TRUE" #DO_SURFACE_CYCLE="TRUE" DO_SPINUP="TRUE" @@ -34,6 +35,13 @@ DO_SMOKE_DUST="FALSE" DO_PARALLEL_PRDGEN="FALSE" DO_GSIDIAG_OFFLINE="FALSE" +# Options for blending or cold2warm start conversion. +ENS_BLENDING_LENGTHSCALE=960 +BLEND="TRUE" # TRUE: Blend RRFS and GDAS EnKF + # FALSE: Don't blend, activate cold2warm start only, and use either GDAS or RRFS +USE_HOST_ENKF="TRUE" # TRUE: Final EnKF (u,v,t,delp,sphum) will be GDAS (no blending) + # FALSE: Final EnKF (u,v,t,delp,sphum) will be RRFS (no blending) + if [[ ${DO_ENSFCST} == "TRUE" ]] ; then EXPT_SUBDIR="rrfs_conus_enfcst" DO_SPINUP="FALSE" @@ -42,7 +50,7 @@ if [[ ${DO_ENSFCST} == "TRUE" ]] ; then DO_POST_PROD="TRUE" fi -EXTRN_MDL_ICS_OFFSET_HRS="30" +EXTRN_MDL_ICS_OFFSET_HRS="6" LBC_SPEC_INTVL_HRS="1" EXTRN_MDL_LBCS_OFFSET_HRS="6" BOUNDARY_LEN_HRS="12" @@ -122,9 +130,9 @@ WTIME_RUN_FCST="00:30:00" WTIME_RUN_FCST_LONG="03:45:00" NNODES_RUN_ANAL="1" -EXTRN_MDL_NAME_ICS="GEFS" +EXTRN_MDL_NAME_ICS="GDASENKF" EXTRN_MDL_NAME_LBCS="GEFS" -FV3GFS_FILE_FMT_ICS="grib2" +FV3GFS_FILE_FMT_ICS="netcdf" FV3GFS_FILE_FMT_LBCS="grib2" EXTRN_MDL_SAVETYPE="GSL" @@ -143,7 +151,7 @@ if [[ ${DO_ENSEMBLE} == "TRUE" ]]; then DO_ENS_GRAPHICS="TRUE" DO_ENKF_RADAR_REF="TRUE" DO_ENSPOST="FALSE" - DO_ENSINIT="TRUE" + DO_ENSINIT="FALSE" # No longer used for 1tstep to get warm starts RADAR_REF_THINNING="2" ARCHIVEDIR="/5year/BMC/wrfruc/rrfs_ens" diff --git a/ush/config_defaults.sh b/ush/config_defaults.sh index 24153655e..617e60b08 100644 --- a/ush/config_defaults.sh +++ b/ush/config_defaults.sh @@ -598,6 +598,8 @@ DA_CYCLE_INTERV="1" RESTART_INTERVAL="1 2" RESTART_INTERVAL_LONG="1 2" CYCL_HRS_HYB_FV3LAM_ENS=( "99" ) +FIRST_BLENDED_CYCLE="18" +FIRST_BLENDED_CYCLE_DATE="YYYYMMDDHH" #----------------------------------------------------------------------- # @@ -2099,6 +2101,28 @@ TILE_SETS="full" # 'DO_ENS_RADDA="TRUE"', the radiance DA must be true, i.e., 'DO_RADDA="TRUE"'. This # is because the radiance DA in EnKF relies the radiance procedures in the GSI-observer, # which is mainly controled by DO_RADDA. +# +# DO_ENS_BLENDING: +# Flag that can enable two things: +# 1) large-scale blending during initialization. +# 2) activate cold2warm start only (replaces ensinit step). +# When this is activated there are two other flags that are relevant: +# 1) BLEND +# 2) USE_HOST_ENKF +# +# BLEND: Only relevant when DO_ENS_BLENDING=TRUE. Flag to perform large scale +# blending during initialization. If this is set to "TRUE", then the RRFS +# EnKF will be blended with the external model ICS using the Raymond filter +# (a low-pass, sixth-order implicit tangent filter). +# TRUE: Blend RRFS and GDAS EnKF +# FALSE: Don't blend, activate cold2warm start only, and use either GDAS or +# RRFS; default +# +# USE_HOST_ENKF: Only relevant when DO_ENS_BLENDING=TRUE and BLEND=FALSE. +# Flag for which EnKF to use during cold2warm start conversion. +# TRUE: Final EnKF will be GDAS (no blending); default +# FALSE: Final EnKF will be RRFS (no blending) +# #----------------------------------------------------------------------- # DO_ENSEMBLE="FALSE" @@ -2118,6 +2142,10 @@ DO_ENSINIT="FALSE" DO_SAVE_DA_OUTPUT="FALSE" DO_GSIDIAG_OFFLINE="FALSE" DO_ENS_RADDA="FALSE" +DO_ENS_BLENDING="FALSE" +ENS_BLENDING_LENGTHSCALE="960" # (Lx) in kilometers +BLEND="FALSE" +USE_HOST_ENKF="TRUE" # #----------------------------------------------------------------------- # diff --git a/ush/generate_FV3LAM_wflow.sh b/ush/generate_FV3LAM_wflow.sh index 6ee814cb4..e4a336a5b 100755 --- a/ush/generate_FV3LAM_wflow.sh +++ b/ush/generate_FV3LAM_wflow.sh @@ -560,6 +560,10 @@ settings="\ # retrospective experiments # 'do_retro': ${DO_RETRO} +# +# large-scale blending EnKF initialization +# + 'do_ens_blending': ${DO_ENS_BLENDING} " # End of "settings" variable. print_info_msg $VERBOSE " diff --git a/ush/setup.sh b/ush/setup.sh index 1e0917bf2..790fb916d 100644 --- a/ush/setup.sh +++ b/ush/setup.sh @@ -808,6 +808,9 @@ NUM_CYCLES="${#ALL_CDATES[@]}" # EXECDIR: # Directory containing various executable files. # +# LIB64DIR: +# Directory containing various library files. +# # TEMPLATE_DIR: # Directory in which templates of various FV3-LAM input files are locat- # ed. @@ -837,6 +840,7 @@ SRC_DIR="${SR_WX_APP_TOP_DIR}/src" PARMDIR="$HOMErrfs/parm" MODULES_DIR="$HOMErrfs/modulefiles" EXECDIR="${SR_WX_APP_TOP_DIR}/bin" +LIB64DIR="${SR_WX_APP_TOP_DIR}/lib64" TEMPLATE_DIR="$USHDIR/templates" if [ "${RUN_ENVIR}" = "nco" ]; then FIXgsm=${FIXgsm:-"$HOMErrfs/fix/am"} @@ -2661,6 +2665,7 @@ SRC_DIR="$SRC_DIR" PARMDIR="$PARMDIR" MODULES_DIR="${MODULES_DIR}" EXECDIR="$EXECDIR" +LIB64DIR="$LIB64DIR" FIXam="$FIXam" FIXLAM="$FIXLAM" FIXgsm="$FIXgsm" diff --git a/ush/templates/FV3LAM_wflow.xml b/ush/templates/FV3LAM_wflow.xml index 8739cdaeb..2365e1369 100644 --- a/ush/templates/FV3LAM_wflow.xml +++ b/ush/templates/FV3LAM_wflow.xml @@ -993,11 +993,15 @@ MODULES_RUN_TASK_FP script. CDATE@Y@m@d@H CYCLE_DIR&CYCLE_BASEDIR;/@Y@m@d@H NWGES_DIR&NWGES_BASEDIR;/@Y@m@d@H + FG_ROOT&FG_ROOT; SLASH_ENSMEM_SUBDIR{{ slash_ensmem_subdir }} + {% if do_ens_blending %} + &FG_ROOT;/@Y@m@d@H/mem#mem#/fcst_fv3lam/RESTART/@Y@m@d.@H0000.coupler.res + {% endif %} &LOGDIR;/&MAKE_GRID_TN;_task_complete.txt @@ -1273,6 +1277,8 @@ MODULES_RUN_TASK_FP script. +{%- endif %} +{%- if do_ens_blending or do_ensinit %}