Version 23 (modified by davestorkey, 11 years ago) (diff) 

New Open Boundary module in NEMO: OBC/BDY merge
Last edited Timestamp?
Introduction
As part of the 2011 workplan, it is planned to merge the OBC and BDY modules to provide a single module to do oneway nesting in NEMO. This page outlines the design for the new module. It is intended as a discussion document. Please feel free to add comments.
There are three main areas of work:
 Create a unified coding framework for the open boundaries.
 Change the data reading routines to use the fldread.F90 module in SBC.
 Implement some of the schemes coded by Jerome Chanut for the Mercator IBIROOS configurations.
I aim to have a branch with a new module with the unified coding framework and the use of fldread.F90 ready by July 2011. This will only have been tested for one or two of the options, but should then be ready for others to work on to implement other options.
Jerome's comments are attached as a Word file.
Dave Storkey
Summary of main features
 The following algorithms will be definitely be included in the first implementation:
 Clamped condition for all variables.
 Zero gradient condition for all variables.
 Flow Relaxation Scheme for all variables. (Maybe explicit and implicit versions).
 Flather radiation condition on barotropic solution.
 Tidal harmonic forcing (included in Flather boundary condition).
 Orlanski radiation condition with normal projection of oblique radiative velocities (NPO) for T and S, 3D and 2D velocities. Plan to code implicit timestepping for greater stability and to avoid OBC restarts (Marchesiello et al 2001).
 sponge layer (enhanced viscosity near the boundary)
 The following algorithms may be included:
 Orlanski radiation condition with normal velocities and fully 2D velocities.
 Other conditions on barotropic solution from Chapman 1985
 There will be two ways of specifying the boundary geometry:
 Using obc_par.h90 to define straightline segments. One will be able to define more than one segment for each of east, west, south, north boundaries.
 Reading in the list of points from a file as is now done in BDY.
 The lowlevel code will use the BDYstyle specification of the boundary geometry (an unstructured list of gridpoints). If obc_par.h90 is used then this will be converted to BDYstyle arrays in the initialisation.
 The read of external boundary data will be handled by fldread.F90 as is done for the SBC fields. This is clean and flexible. In the future this will permit online interpolation of boundary data.
 It will be possible to define more than one open boundary set. This will allow different boundary conditions to be applied on different boundaries.
Questions
 For the timesplitting solution for the free surface (ln_dynspg_ts=.true.) boundary conditions are applied to the barotropic (nn_dyn2d) and baroclinic (nn_dyn3d) solutions separately. For the filtered free surface should we insist that we apply barotropic and baroclinic boundary conditions separately or do we allow boundary conditions to be applied to the full velocity field?
Control and coding structure
The namelist for the new module would look like this:
! \&namobc ! oneway open boundaries ("key_obc") ! nb_obc = 2 ! Number of open boundary sets ln_read = .false.,.true. ! Read in boundary geometry or not cn_file = '','coordinates_bdy.nc' ! Name of file with boundary geometry ln_mask = .false.,.true., ! boundary mask from cn_mask (T), boundaries are on edges of domain (F) cn_mask = '','mask_bdy.nc' ! name of mask file (ln_mask=T) nn_dyn3d = 1, 0, ! Choice of schemes for 3D velocities nn_dyn2d = 1, 0, ! Choice of schemes for barotropic solution (ln_dynspg_ts=.true.) nn_tra = 1, 1, ! Choice of schemes for T and S nn_ice_lim2 = 0, 0, ! Choice of schemes for ice (LIM2) ln_vol = .false.,.false., ! total volume correction (see volbdy parameter) ln_tides = .false.,.false., ! Apply tidal harmonic forcing with Flather condition nn_rimwidth = 9, 1, ! width of the relaxation zone for Flow Relaxation nn_dmp2d_in = 1, 1, ! damping timescale (days) for 2D variables for inward radiation velocity or Flow Relaxation nn_dmp2d_out = 100, 1, ! damping timescale (days) for 2D variables outward radiation velocity nn_dmp3d_in = 1, 1, ! damping timescale (days) for 3D variables inward radiation velocity or Flow Relaxation nn_dmp3d_out = 100, 1, ! damping timescale (days) for 3D variables outward radiation velocity nn_dtactl = 1, 1, ! = 0, bdy data are equal to the initial state ! = 1, bdy data are read in 'bdydata .nc' files nn_volctl = 0, 0, ! = 0, the total water flux across open boundaries is zero ! = 1, the total volume of the system is conserved / ! \&namobc_dta ! oneway open boundaries ("key_obc") ! ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'daily'/ ! weights ! rotation ! ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! bn_tem = 'obcdta_grid_T' , 24 , 'votemper', .true. , .false., 'daily' , '' , '' bn_sal = 'obcdta_grid_T' , 24 , 'vosaline', .true. , .false., 'daily' , '' , '' bn_uvel = 'obcdta_grid_U' , 24 , 'vozocrtx', .true. , .false., 'daily' , '' , '' bn_vvel = 'obcdta_grid_V' , 24 , 'vomecrty', .true. , .false., 'daily' , '' , '' bn_ssh = 'obcdta_bt_grid_T' , 6 , 'sossheig', .true. , .false., 'daily' , '' , '' bn_ubar = 'obcdta_bt_grid_U' , 6 , 'vozocrtx', .true. , .false., 'daily' , '' , '' bn_vbar = 'obcdta_bt_grid_V' , 6 , 'vomecrty', .true. , .false., 'daily' , '' , '' cn_dir = './obcdta' ! root directory for the location of the boundary data files / ! \&namobc_dta ! oneway open boundaries ("key_obc") ! ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'daily'/ ! weights ! rotation ! ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! bn_tem = 'obcdta_clim_grid_T' , 1 , 'votemper', .true. , .true., 'yearly' , '' , '' bn_sal = 'obcdta_clim_grid_T' , 1 , 'vosaline', .true. , .true., 'yearly' , '' , '' cn_dir = './obcdta' ! root directory for the location of the boundary data files /
It is possible to define more than one open boundary set in case you want to use different boundary conditions for different boundaries. For instance one might want to apply different algorithms on the east and west boundaries of a domain, or one might want to apply boundary conditions from an external model over part of the open boundary but climatological boundary conditions over another part of the open boundary (as in the example above).
For each open boundary set you define which conditions to apply to each set of variables using nn_dyn3d, nn_tra etc. For example:
 nn_dyn3d = 0 : Apply no boundary condition (land boundary)
 nn_dyn3d = 1 : Clamped/relaxation boundary condition
 nn_dyn3d = 2 : Orlanski NPO radiation condition
 nn_dyn3d = 3 : Orlanski fully 2d radiation condition
 etc.
There is a separate &namobc_dta namelist for each boundary set, which defines the input data files using the FLD structure from fldread.F90. (This will allow for online interpolation of boundary conditions in the future).
For each set of variables there is a module which applies the boundary conditions each timestep, eg. for T and S you have obcdyn3d.F90 which looks like this:
MODULE obcdyn3d CONTAINS SUBROUTINE obc_dyn3d( kt ) DO ib_set = 1, nb_set IF ( using Orlanski radiation ) THEN CALL obc_rad( kt, ib_set ) ENDIF SELECT CASE ( nn_dyn3d(ib_set) ) CASE(0) CYCLE CASE(1) CALL obc_dyn3d_frs( kt, ib_set ) CASE(2) CALL obc_dyn3d_orl( kt, ib_set ) CASE(3) CALL obc_dyn3d_orl( kt, ib_set ) END SELECT END DO END SUBROUTINE SUBROUTINE obc_dyn3d_frs( kt, ib_set ) ... SUBROUTINE obc_dyn3d_orl( kt, ib_set ) ...
The module to read in the boundary data from files would be rewritten to use fldread.F90. For each boundary stream a TYPE FLD structure array would be constructed depending on which fields need to be read in. We can merge the obc_dta and obc_dta_bt routines by making the barotropic subloop counter jit an optional argument.
MODULE obcdta USE fldread INTEGER :: nb_dta(nb_obc_max) TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: bf CONTAINS SUBROUTINE obc_dta( kt, jit ) INTEGER :: kt ! time step index INTEGER, OPTIONAL :: jit ! barotropic subloop index TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: blf_i IF( kt == nit000 ) THEN ! First call kt=nit000 REWIND ( numnam ) jfld = 0 jfld_prev = 0 DO ib_set = 1, nb_set IF ( nn_dtactl(ib_set) .eq. 1 ) THEN ! set file information cn_dir = './' ! directory in which the model is executed ! ... default values (NB: frequency positive => hours, negative => months) ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! bn_tem = 'bn_tem' , 24 , 'bn_tem' , .true. , .false. , 'daily' , '' , '' bn_sal = 'bn_sal' , 24 , 'bn_sal' , .true. , .false. , 'daily' , '' , '' bn_uvel = 'bn_uvel' , 24 , 'bn_uvel', .true. , .false. , 'daily' , '' , '' bn_vvel = 'bn_vvel' , 24 , 'bn_vvel', .true. , .false. , 'daily' , '' , '' bn_ssh = 'bn_ssh' , 24 , 'bn_ssh' , .true. , .false. , 'daily' , '' , '' bn_ubar = 'bn_ubar' , 24 , 'bn_ubar', .true. , .false. , 'daily' , '' , '' bn_vbar = 'bn_vbar' , 24 , 'bn_vbar', .true. , .false. , 'daily' , '' , '' READ ( numnam, namobc_dta ) ! Only read in necessary fields for this set. IF( nn_dyn2d(ib_set) .gt. 0 ) THEN jfld = jfld + 1 blf_i(jfd) = bn_ssh jfld = jfld + 1 blf_i(jfld) = bn_ubar jfld = jfld + 1 blf_i(jfld) = bn_vbar ENDIF IF( nn_tra(ib_set) .gt. 0 ) THEN .... nb_dta(ib_set) = jfld  jfld_prev jfld_prev = jfld END IF END DO ALLOCATE( bf(SUM(nb_dta)), STAT=ierror ) ! set sf structure jfld = 0 DO ib_set = 1, nb_set IF ( nn_dtactl(ib_set) .eq. 1 ) THEN DO ji= 1, nb_dta(ibset) jfld = jfld+1 ALLOCATE( sf(jfld)%fnow(nblen_dta(ib_set),1,1) ) IF( slf_i(jfld)%ln_tint ) ALLOCATE( sf(ji)%fdta(nblen_dta(ib_set),1,1,2) ) END DO END IF END DO CALL fld_fill( bf, blf_i, ... ) END IF jstart = 1 jfld = 0 DO ib_set = 1, nb_set IF( nn_dtactl(ib_set) = 1 ) THEN IF( PRESENT(jit) ) THEN ! Update barotropic boundary conditions only ! jit is optional argument for fld_read IF( nn_dyn2d(ib_set) .gt. 0 ) THEN jend = jstart + 2 CALL fld_read( kt, nn_fobc, bf(jstart:jend), jit=jit ) ENDIF ELSE jend = jstart + nb_obc(ib_set)  1 CALL fld_read( kt, nn_fobc, bf(jstart:jend ) ) ENDIF IF( nn_dyn2d .gt. 0 ) THEN jpfld = jpfld + 1 sshbdy(ib_set,:) = bf(jpfld)%fnow jpfld = jpfld + 1 ubar(ib_set,:) = bf(jpfld)%fnow ... ENDIF IF( .NOT. PRESENT(jit) ) THEN IF( nn_tra(ib_set) .gt. 0 ) THEN jpfld = jpfld + 1 tbdy(ib_set,:) = bf(jpfld)%fnow .... END IF ! nn_dtactl(ib_set) = 1 END DO ! ib_set END SUBROUTINE
Boundary geometry and initialisation
The internal code of the new module will use 1D unstructured arrays to specify the boundary as is currently done in BDY. For the Tpoints along the boundary these are nbit, nbjt and nbrt, which specify the (i,j) coordinates of each point and the discrete distance from the boundary. Each boundary set can be initialised using straightline segments in obc_par.F90 (ln_read=.false.) or read in from a coordinates_bdy.nc file (ln_read=.true.). If the boundary coordinates are read in from obc_par.F90 then the nbi, nbj, nbr arrays will be generated internally (and if necessary a coordinates_bdy.nc file could be written out).
For the Flow Relaxation Scheme we require relaxation to external data within the model, not just at the boundary. The BDY module includes these points in the specification of the boundary arrays and uses nbr to tell the discrete distance from the model boundary. The obc_par.F90 specification only specifies the points along the edge of the model domain. If nn_rimwidth > 1 the internal points will be added to the BDYstyle arrays internally.
The obc_par.F90 file will be the easiest way to define rectangular boundaries. More than one segment can be specified for each of the north, east, south and west boundaries as follows. Where there are more than one boundary set, the jpieobtype arrays are specified as 1D arrays and the nobcsegetype arrays are used to split the numbers between the different boundary sets.
MODULE obc_par INTEGER, PARAMETER :: & !: ! Number of open boundary segments nobcsege = 2, 0 INTEGER, PARAMETER, DIMENSION(nobcsege) :: & !: ! First column: Med Sea; Second: Baltic jpieob = (/ 1037 , jpiglo2 /), & !: ilocalization of the East open boundary jpjedt = (/ 473 , 1476 /), & !: jstarting indice of the East open boundary jpjeft = (/ 855 , 1548 /) !: jending indice of the East open boundary !! * WEST open boundary INTEGER, PARAMETER :: & !: ! Number of open boundary segments nobcsegw = 1, 0 INTEGER, PARAMETER, DIMENSION(nobcsegw) :: & jpiwob = (/ 2 /), & !: ilocalization of the West open boundary jpjwdt = (/ 2 /), & !: jstarting indice of the West open boundary jpjwft = (/1875/) !: jending indice of the West open boundary ...
For unstructured boundaries, the easiest method will be to generate them offline and read in the boundary definition from the coordinates_bdy.nc file, which looks like this:
netcdf coordinates_bdy { dimensions: xbT = 8485 ; xbU = 8436 ; xbV = 8455 ; xbF = 8062 ; yb = 1 ; variables: int nbit(yb, xbT) ; int nbiu(yb, xbU) ; int nbiv(yb, xbV) ; int nbif(yb, xbF) ; int nbjt(yb, xbT) ; int nbju(yb, xbU) ; int nbjv(yb, xbV) ; int nbjf(yb, xbF) ; int nbrt(yb, xbT) ; int nbru(yb, xbU) ; int nbrv(yb, xbV) ; int nbrf(yb, xbF) ; float glamt(yb, xbT) ; glamt:units = "degrees_east" ; float glamu(yb, xbU) ; glamu:units = "degrees_east" ; float glamv(yb, xbV) ; glamv:units = "degrees_east" ; float glamf(yb, xbF) ; glamf:units = "degrees_east" ;
Pseudo code for implementation of Orlanski radiation conditions using BDYstyle arrays
... To be done ...
Attachments (1)

NEMO_Bdy_wiki.doc
(55.0 KB) 
added by davestorkey 11 years ago.
Jerome's comments on OBCBDY design
Download all attachments as: .zip