| Description |
There is a lot of information in this update, please read carefully and understand if you use any of the algorithms mentioned below.
This update is crucial because it includes the capability allows the forward nonlinear solution to be time split into multiple files to avoid creating huge NetCDF files in applications with large grids. Recall that the nonlinear trajectory is used to linearize the tangent linear and adjoint models used in 4D-Var data assimilation, observation impacts, observation sensitivities, adjoint sensitivities, generalized stability analysis, and other adjoint-based algorithms.
- Two subroutines (edit_multifile and edit_file_struct) are included in edit_multifile.F to manipulate the basic state forward trajectory used in various adjoint-based drivers. The TYPE(T_IO) structure can be deallocated and allocate with multiple file capabilities if so desired. For example in edit_multifile, we have:
DO ng=1,Ngrids
!
SELECT CASE (TRIM(ADJUSTL(task)))
!
! Load HIS information into the FWD structure so it can be used to
! process the NLM background trajectory by the ADM and TLM kernels.
!
CASE ('HIS2FWD')
IF (ndefHIS(ng).gt.0) THEN
IF (HIS(ng)%ncid.ne.-1) THEN
CALL netcdf_close (ng, iNLM, HIS(ng)%ncid)
END IF
Nfiles=ntimes(ng)/ndefHIS(ng)
IF (nHIS(ng).eq.ndefHIS(ng)) Nfiles=Nfiles+1
CALL edit_file_struct (ng, Nfiles, FWD)
DO ifile=1,Nfiles
FWD(ng)%files(ifile)=TRIM(HIS(ng)%files(ifile))
END DO
FWD(ng)%name=TRIM(FWD(ng)%files(1))
ELSE
FWD(ng)%ncid=HIS(ng)%ncid
FWD(ng)%name=TRIM(HIS(ng)%name)
FWD(ng)%files(1)=TRIM(HIS(ng)%name)
END IF
CASE ('HIS2BLK')
...
CASE ('FWD2BLK')
...
CASE ('FWD2HIS')
...
CASE ('FCTA2FWD')
...
CASE ('FCTA2BLK')
...
CASE ('FCTB2FWD')
...
CASE ('FCTB2BLK')
...
CASE ('TLM2FWD')
...
END SELECT
END DO
As you can see, there are several possibilities in the manipulations between file structures. In the calling driver, we need to add the specific instructions for such handling. For example in w4dpsas_ocean.h, we have:
!
! Set structure for the nonlinear forward trajectory to be processed
! by the tangent linear and adjoint models. Also, set switches to
! process the FWD structure in routine "check_multifile". Notice that
! it is possible to split solution into multiple NetCDF files to reduce
! their size.
!
CALL edit_multifile ('HIS2FWD')
IF (FoundError(exit_flag, NoError, __LINE__, &
& __FILE__)) RETURN
DO ng=1,Ngrids
LreadFWD(ng)=.TRUE.
END DO
#if defined BULK_FLUXES && defined NL_BULK_FLUXES
!
! Set structure for the nonlinear surface fluxes to be processed by
! by the tangent linear and adjoint models. Also, set switches to
! process the BLK structure in routine "check_multifile". Notice that
! it is possible to split solution into multiple NetCDF files to reduce
! their size.
!
CALL edit_multifile ('HIS2BLK')
IF (FoundError(exit_flag, NoError, __LINE__, &
& __FILE__)) RETURN
DO ng=1,Ngrids
LreadBLK(ng)=.TRUE.
END DO
#endif
- The routine check_multifile was updated to include more TYPE(T_IO) structures to process for information about multiple files or not. Two new generic subroutines multifile_info_s1d and multifile_info_s2d were written to inquire about 1D and 2D TYPE(T_IO) structures, repectively. They are called in a compact way in check_multifile:
#ifdef FORWARD_READ
!
!-----------------------------------------------------------------------
! Nonlinear model forward trajectory data.
!-----------------------------------------------------------------------
!
IF (LreadFWD(ng)) THEN
CALL multifile_info_s1d (ng, model, 'Forward', FWD)
IF (FoundError(exit_flag, NoError, __LINE__, &
& __FILE__)) RETURN
END IF
#endif
...
#ifdef FRC_FILE
!
!-----------------------------------------------------------------------
! Input forcing data.
!-----------------------------------------------------------------------
!
file_type='Forcing'
max_files=MAXVAL(nFfiles)
CALL multifile_info_s2d (ng, model, file_type, nFfiles, &
& max_files, FRC)
IF (FoundError(exit_flag, NoError, __LINE__, &
& __FILE__)) RETURN
#endif
- Several switches (LreadADM, LreadBLK, LreadFWD, and LreadTLM) are introduced in mod_scalar.F to activate the processing of information about file structures at the right time in the algorithms.
- The TYPE(T_IO) file structure definition in mod_param.F now includes the head field, which remains unchanged after initialization and has the I/O heading prefix filename. It is used as follows:
WRITE (HIS(ng)%name,10) TRIM(FWD(ng)%head), outer
lstr=LEN_TRIM(HIS(ng)%name)
HIS(ng)%base=HIS(ng)%name(1:lstr-3)
The HIS(ng)%base is used in output.F to create the multile files. For example, in 4D-Var we can have:
wc13_fwd_outer0_0001.nc
wc13_fwd_outer0_0002.nc
wc13_fwd_outer0_0003.nc
wc13_fwd_outer0_0004.nc
wc13_fwd_outer1_0001.nc
wc13_fwd_outer1_0002.nc
wc13_fwd_outer1_0003.nc
wc13_fwd_outer1_0004.nc
Here, the basic state forward trajectory is split in daily files having data records every two-hours in a 4D-Var 4-days data assimilation cycle for each outer loop.
WARNING: If not using output history multi-files, NDEFHIS=0 in roms.in, we will have
wc13_fwd_outer0.nc instead of wc13_fwd_0000.nc
wc13_fwd_outer1.nc instead of wc13_fwd_0001.nc
as in previous versions of the code. So you need to update your post-processing scripts.
- The multipe output files are activated in roms.in by setting any of their parameters (NDEFHIS, NDEFQCK, NDEFAVG, NDEFDIA) to be greater than zero. ROMS creates new files every specified number of timesteps. The parameters NDEFADJ and NDEFTLM are special and some algorithms do not allow its values to be greater than zero.
- Started to solve the issue with DSTART in the adjoint-based algorithms. The management of DSTART has evolved over the years. It is set in roms.in standard input file. Nowadays, we need to set its value to the actual initialization time (in days) every time that we run any of the 4D-Var drivers or their associated analysis tools. The reason for it is that tl_initial, rp_initial, and ad_initial rarely reads an initial conditions NetCDF file for the TLM, RPM, and ADM, respectively. Usually, those initial conditions are initialized from rest (zeros) or computed by the 4D-Var driver. Therefore, we need to use DSTART to compute the correct and crucial time-stepping parameters: time(ng), tdays(ng), ntstart(ng), ntend(ng), and ntfirst(ng) since we are not reading an initial NetCDF file to set time(ng) for those kernels.
My current strategy is to use the initialization time from the nonlinear model or forward trajectory, which is now saved in INItime(ng). It is already coded but not activated because it requires rigorous testing. I am not fully satisfied yet. This capability is activated with the new GENERIC_DSTART C-preprocessing option. Please do not activate this option! You may get different and wrong results. There is a tricky bug somewhere that I haven't been able to find. In tl_initial.F, we have:
# ifdef GENERIC_DSTART
!
! Rarely, the tangent linear model is initialized from a NetCDF file,
! so we do not know its actual initialization time. Usually, it is
! computed from DSTART, implying that its value is correct in the ROMS
! input script. Therefore, the user needs to check and update its value
! to every time that ROMS is executed. Alternatively, if available, we
! can use the initialization time from the nonlinear model, INItime.
! This variable is assigned when computing or processing the basic
! state trajectory needed to linearize the adjoint model.
!
IF (INItime(ng).lt.0.0_dp) THEN
my_dstart=dstart ! ROMS input script
ELSE
my_dstart=INItime(ng)/86400.0_dp ! NLM IC time is known
END IF
tdays(ng)=my_dstart
# else
tdays(ng)=dstart
# endif
time(ng)=tdays(ng)*day2sec
The reason why I am looking at it is that we want an unrestrained value of DSTART because of the multiple-file option that now is allowed in the adjoint-based algorithms. Recall that DSTART controls the split file enumeration and the starting value for iic(ng). If DSTART is set to the initial conditions time value, the count always starts at zero.
- Corrected few tiny bugs in the 2D applications set-up. Many thanks to Brian Powell for bringing these issues to my attention.
|
| Description |
- The ESM run sequence is now set, by default, from a new input configuration file (Keyword: CONFname) in the coupling standard input script (coupling_esmf.in). It is a more elegant option than the hundreds, if not thousands, lines of code in module esmf_esm.F. It will facilitate complex RunSequence execution with more than one-time loop (time slots) and an increased number of coupled components. The NUOPC layer generates the RunSequence automatically. It can be changed without recompiling. It also facilitates updates to newer versions of the ESMF/NUOPC library.
- Examples of RunSequence templates are provided for coupling ROMS to either COAMPS or WRF explicitly or semi-implicitly:
- ESM/coamps_explicit.runconfig
- ESM/coamps_implicit.runconfig
- ESM/wrf_explicit.runconfig
- ESM/wrf_implicit.runconfig
For the DATA-WRF-ROMS semi-implicit coupling with a single time loop (denoted by the at symbol @), we have:
#
# svn $Id$
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Copyright (c) 2002-2019 The ROMS/TOMS Group :::
# Licensed under a MIT/X style license :::
# See License_ROMS.txt :::
#::::::::::::::::::::::::::::::::::::::::::::::::::::: Hernan G. Arango :::
# :::
# DATA-WRF-ROMS semi-implicit ESMF coupling: :::
# :::
# This ESMF configuration file includes the Run Sequence to couple :::
# DATA, WRF, and ROMS components. All the components interact with :::
# the same coupling time step. The connector from ROMS to WRF is :::
# explicit, whereas the connector from WRF to ROMS is semi-implicit. :::
# Usually, the timestep of the atmosphere kernel is smaller than that :::
# for the ocean. Therefore, WRF export files are time-averaged over :::
# the coupling interval, which is the same as the ROMS timestep. It :::
# is semi-implicit because ROMS right-hand-side terms are forced with :::
# n+1/2 WRF fields because of the time-averaging. :::
# :::
# The time-averaging is activated with WRF_TIMEAVG in the ROMS build :::
# script. The WRF "namelist.input" needs to include the parameters :::
# for the RAMS averaged diagnostics in the &time_control section. :::
# The User needs to check that the single coupling interval specified :::
# here is the same as for the ROMS coupling standard input script :::
# (keyword TimeStep). :::
# :::
# It is not necessary to specify options to the connectors here for :::
# the "remapMethod" since it is specified in the input coupling :::
# metada file (keyword CPLname) for each exchanged field. :::
# :::
# The component label is specific to how is known in the ESMF driver. :::
# It is set in subroutine allocate_esmf_esm, variable MODELS(:)%name. :::
# :::
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Hurricane Irene Application
runSeq::
@60 # 1 min coupling time step, single time loop
DATA -> WRF_01 # DATA to WRF connector
DATA
ROMS_01 -> WRF_01 # ROMS to WRF connector
WRF_01
WRF_01 -> ROMS_01 # WRF to ROMS connector
ROMS_01
@
::
As you can see, the syntax is straightforward and implies thousands of line of code that we can avoid. See NUOPC_refdoc_*.doc for details about the synthax. The User needs to make sure that the coupling time step in seconds (@60) has the same value of TimeStep specified in coupling_esmf.in:
! Set ESM coupling driver clock. If not re-starting, set StartTime and
! RestartTime to the same values. Set coupling time interval (TimeStep)
! to the largest value for the field exchange between activated ESM
! components. Check glossary below for more information.
!
! YYYY MM DD hh mm ss
ReferenceTime = 2006 01 01 00 00 00 ! Simulation reference time
StartTime = 2011 08 27 00 00 00 ! Simulation start time
RestartTime = 2011 08 27 00 00 00 ! Simulation restart time
StopTime = 2011 08 29 00 00 00 ! Simulation stop time
TimeStep = 0000 00 00 00 01 00 ! Coupler driver interval (60 sec ROMS DT)
Here, ROMS DT=60 and larger than the one for WRF. Therefore, the exchange of fields between DATA, WRF, and ROMS components is done every 60 seconds. The execution will stop if both values are not the same.
- A new C-preprocessing option NO_ESMF_CONFIG to preserve the hundreds of line of coding in esmf_esm.F. Such code and option will be removed in the future.
|
| Description |
- The Melding coefficients are used to combine fields from DATA and ESM components. The weight factors are read from the input NetCDF specified in the WeightsFile(atmos) keyword of coupling_esmf.in The user has full control of how the merging is done. It is recommended to provide a gradual transition between the two components.
Recall that the DATA component supplies needed data to a particular ESM component. For example, it may export data to the atmosphere model at locations not covered by the other ESM components because of smaller grid coverage. If the atmosphere and ocean model grids are incongruent, the atmosphere component needs to import sea surface temperature (SST) on those grid points not covered by the ocean component. Thus, the weighting coefficients are used to merge the SST data:
SST_atm(:,:) = Cesm(:,:) * SST_esm(:,:) + Cdat(:,:) * SST_dat(:,:)
where
Cesm(:,:) + Cdat(:,:) = 1.
In coupling_esmf.in, we now have:
! Weighted coefficients for the merging of DATA component fields with the same
! field from other ESM components. Melding coefficients are positive and MUST
! add to unity. They are read from an input NetCDF file ('WeightsFile'). The
! user needs to specify the NetCDF variable names for the weights for the DATA
! ('VnameDATA') and ESM ('VnameESM') components. Currently, the weight values
! are only needed by the atmosphere component.
WeightsFile(atmos) == meld_weights_atm.nc
VnameDATA(atmos) = data_weight
VnameESM(atmos) = ocean_weight
We no longer use the WeightDAT and WeightESM constant parameters from previous versions.
- Several functions to compute melding weights between DATA and ROMS coupling components for exporting fields to the atmosphere component.
- c_weithts.m: Creates a new weight factor NetCDF file for melding DATA and OCEAN components during ESMF coupling.
- coamps_weights.m: Computes and writes the COAMPS melding weights to combine fields from DATA and ROMS components after ESMF regridding to the atmosphere grid.
- smooth_weights.m: Computes smooth melding weights to combine fields from DATA and OCEAN components. The merging factors change gradually in an area next to the OCEAN grid open boundary.
- wrf_weights.m: Computes and writes the WRF melding weights to combine fields from DATA and ROMS components after ESMF regridding to the atmosphere grid.
Many thanks to John Wilkin for his help in coding the generic Matlab script smooth_weights.m to compute the varying weight factors at each grid cell of the atmosphere grid.
- The coamps_weights.m is used to compute the melding weight factors for COAMPS 15km resolution grid that combines SST from the DATA (HyCOM) and ROMS components in the Indian Ocean coupled application.
Here, N=200 is the number of convolutions used to smooth the transition between DATA and ROMS components.
- The wrf_weights.m is used to compute the melding weight factors for WRF 6km resolution grid that combines SST from the DATA (HyCOM) and ROMS components in a U.S. East Coast application to study Hurricanes Irene and Sandy.
- The c_weights.m creates a simple NetCDF file containing both the rigid and smooth weights for the DATA and OCEAN components.
netcdf wrf6km_meld_weights {
dimensions:
lon = 319 ;
lat = 279 ;
variables:
double lon(lat, lon) ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
lon:units = "degrees_east" ;
double lat(lat, lon) ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
lat:units = "degrees_north" ;
double mask(lat, lon) ;
mask:long_name = "land-sea mask" ;
mask:flag_values = 0., 1. ;
mask:flag_meanings = "land sea" ;
mask:coordinates = "lon lat" ;
double data_weight_rigid(lat, lon) ;
data_weight_rigid:long_name = "DATA component rigid melding weights" ;
data_weight_rigid:valid_min = 0. ;
data_weight_rigid:valid_max = 1. ;
data_weight_rigid:coordinates = "lon lat" ;
double ocean_weight_rigid(lat, lon) ;
ocean_weight_rigid:long_name = "OCEAN component rigid melding weights" ;
ocean_weight_rigid:valid_min = 0. ;
ocean_weight_rigid:valid_max = 1. ;
ocean_weight_rigid:coordinates = "lon lat" ;
double data_weight_smooth(lat, lon) ;
data_weight_smooth:long_name = "DATA component smooth melding weights" ;
data_weight_smooth:valid_min = 0. ;
data_weight_smooth:valid_max = 1. ;
data_weight_smooth:coordinates = "lon lat" ;
double ocean_weight_smooth(lat, lon) ;
ocean_weight_smooth:long_name = "OCEAN component smooth melding weights" ;
ocean_weight_smooth:valid_min = 0. ;
ocean_weight_smooth:valid_max = 1. ;
ocean_weight_smooth:coordinates = "lon lat" ;
// global attributes:
:Conventions = "CF-1.0" ;
:title = "WRF coupling import field melding weights between `DATA and ROMS components" ;
:wrf_grid = "/s1/arango/ocean/repository/Projects/irene/Data/WRF/6km/irene_wrf_inp_d01.nc" ;
:roms_grid = "/s1/arango/ocean/repository/Projects/om/grid_doppio_JJA_v12n.nc" ;
:convolutions = "150" ;
:history = "Weights file created using Matlab script /Users/arango/ocean/repository/matlab/coupling/wrf_weights.m 24-Oct-2019 21:39:02" ;
- Added and updated several Matlab scripts.
|
