New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Ignore:
Timestamp:
2012-01-28T17:44:18+01:00 (12 years ago)
Author:
rblod
Message:

Merge of 3.4beta into the trunk

Location:
trunk/NEMOGCM/NEMO/OPA_SRC/BDY
Files:
1 deleted
8 edited
3 copied

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r2715 r3294  
    77   !!            3.0  !  2008-04  (NEMO team)  add in the reference version      
    88   !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions 
     9   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_bdy  
    … …  
    1920   PUBLIC 
    2021 
     22   TYPE, PUBLIC ::   OBC_INDEX    !: Indices and weights which define the open boundary 
     23      INTEGER,          DIMENSION(jpbgrd) ::  nblen 
     24      INTEGER,          DIMENSION(jpbgrd) ::  nblenrim 
     25      INTEGER, POINTER, DIMENSION(:,:)   ::  nbi 
     26      INTEGER, POINTER, DIMENSION(:,:)   ::  nbj 
     27      INTEGER, POINTER, DIMENSION(:,:)   ::  nbr 
     28      INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap 
     29      REAL   , POINTER, DIMENSION(:,:)   ::  nbw 
     30      REAL   , POINTER, DIMENSION(:)     ::  flagu 
     31      REAL   , POINTER, DIMENSION(:)     ::  flagv 
     32   END TYPE OBC_INDEX 
     33 
     34   TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data 
     35      REAL, POINTER, DIMENSION(:)     ::  ssh 
     36      REAL, POINTER, DIMENSION(:)     ::  u2d 
     37      REAL, POINTER, DIMENSION(:)     ::  v2d 
     38      REAL, POINTER, DIMENSION(:,:)   ::  u3d 
     39      REAL, POINTER, DIMENSION(:,:)   ::  v3d 
     40      REAL, POINTER, DIMENSION(:,:)   ::  tem 
     41      REAL, POINTER, DIMENSION(:,:)   ::  sal 
     42#if defined key_lim2 
     43      REAL, POINTER, DIMENSION(:)     ::  frld 
     44      REAL, POINTER, DIMENSION(:)     ::  hicif 
     45      REAL, POINTER, DIMENSION(:)     ::  hsnif 
     46#endif 
     47   END TYPE OBC_DATA 
     48 
    2149   !!---------------------------------------------------------------------- 
    2250   !! Namelist variables 
    2351   !!---------------------------------------------------------------------- 
    24    CHARACTER(len=80) ::   cn_mask        !: Name of unstruct. bdy mask file 
    25    CHARACTER(len=80) ::   cn_dta_frs_T   !: Name of unstruct. bdy data file at T points for FRS conditions 
    26    CHARACTER(len=80) ::   cn_dta_frs_U   !: Name of unstruct. bdy data file at U points for FRS conditions 
    27    CHARACTER(len=80) ::   cn_dta_frs_V   !: Name of unstruct. bdy data file at V points for FRS conditions 
    28    CHARACTER(len=80) ::   cn_dta_fla_T   !: Name of unstruct. bdy data file at T points for Flather scheme 
    29    CHARACTER(len=80) ::   cn_dta_fla_U   !: Name of unstruct. bdy data file at U points for Flather scheme 
    30    CHARACTER(len=80) ::   cn_dta_fla_V   !: Name of unstruct. bdy data file at V points for Flather scheme 
     52   CHARACTER(len=80), DIMENSION(jp_bdy) ::   cn_coords_file !: Name of bdy coordinates file 
     53   CHARACTER(len=80)                    ::   cn_mask_file   !: Name of bdy mask file 
    3154   ! 
    32    LOGICAL ::   ln_tides = .false.    !: =T apply tidal harmonic forcing along open boundaries 
    33    LOGICAL ::   ln_vol  = .false.     !: =T volume correction              
    34    LOGICAL ::   ln_mask = .false.     !: =T read bdymask from file 
    35    LOGICAL ::   ln_clim = .false.     !: =T bdy data files contain  1 time dump  (-->bdy forcing will be constant)  
    36    !                                  !                         or 12 months     (-->bdy forcing will be cyclic)  
    37    LOGICAL ::   ln_dyn_fla  = .false. !: =T Flather boundary conditions on barotropic velocities 
    38    LOGICAL ::   ln_dyn_frs  = .false. !: =T FRS boundary conditions on velocities 
    39    LOGICAL ::   ln_tra_frs  = .false. !: =T FRS boundary conditions on tracers (T and S) 
    40    LOGICAL ::   ln_ice_frs  = .false. !: =T FRS boundary conditions on seaice (leads fraction, ice depth, snow depth) 
     55   LOGICAL, DIMENSION(jp_bdy) ::   ln_coords_file           !: =T read bdy coordinates from file;  
     56   !                                                        !: =F read bdy coordinates from namelist 
     57   LOGICAL                    ::   ln_mask_file             !: =T read bdymask from file 
     58   LOGICAL                    ::   ln_vol                   !: =T volume correction              
    4159   ! 
    42    INTEGER ::   nn_rimwidth = 7       !: boundary rim width 
    43    INTEGER ::   nn_dtactl   = 1       !: = 0 use the initial state as bdy dta ; = 1 read it in a NetCDF file 
    44    INTEGER ::   nn_volctl   = 1       !: = 0 the total volume will have the variability of the surface Flux E-P  
    45    !                                  !  = 1 the volume will be constant during all the integration. 
    46  
     60   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
     61   INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme 
     62   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
     63   !                                                        !  = 1 the volume will be constant during all the integration. 
     64   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d                 ! Choice of boundary condition for barotropic variables (U,V,SSH) 
     65   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d_dta           !: = 0 use the initial state as bdy dta ;  
     66                                                            !: = 1 read it in a NetCDF file 
     67                                                            !: = 2 read tidal harmonic forcing from a NetCDF file 
     68                                                            !: = 3 read external data AND tidal harmonic forcing from NetCDF files 
     69   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities  
     70   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d_dta           !: = 0 use the initial state as bdy dta ;  
     71                                                            !: = 1 read it in a NetCDF file 
     72   INTEGER, DIMENSION(jp_bdy) ::   nn_tra                   ! Choice of boundary condition for active tracers (T and S) 
     73   INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta             !: = 0 use the initial state as bdy dta ;  
     74                                                            !: = 1 read it in a NetCDF file 
     75#if defined key_lim2 
     76   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables  
     77   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2_dta          !: = 0 use the initial state as bdy dta ;  
     78                                                            !: = 1 read it in a NetCDF file 
     79#endif 
     80   ! 
     81    
    4782   !!---------------------------------------------------------------------- 
    4883   !! Global variables 
    … …  
    5287   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points 
    5388 
     89   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary 
     90 
     91   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !:  
     92   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !:  
     93   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields  
     94   REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !:  
     95   REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !:  
     96 
    5497   !!---------------------------------------------------------------------- 
    55    !! Unstructured open boundary data variables 
     98   !! open boundary data variables 
    5699   !!---------------------------------------------------------------------- 
    57    INTEGER, DIMENSION(jpbgrd) ::   nblen    = 0           !: Size of bdy data on a proc for each grid type 
    58    INTEGER, DIMENSION(jpbgrd) ::   nblenrim = 0           !: Size of bdy data on a proc for first rim ind 
    59    INTEGER, DIMENSION(jpbgrd) ::   nblendta = 0           !: Size of bdy data in file 
    60100 
    61    INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbi, nbj        !: i and j indices of bdy dta 
    62    INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbr             !: Discrete distance from rim points 
    63    INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbmap           !: Indices of data in file for data in memory  
    64      
    65    REAL(wp) ::   bdysurftot                               !: Lateral surface of unstructured open boundary 
    66  
    67    REAL(wp), DIMENSION(jpbdim)        ::   flagu, flagv   !: Flag for normal velocity compnt for velocity components 
    68    REAL(wp), DIMENSION(jpbdim,jpbgrd) ::   nbw            !: Rim weights of bdy data 
    69  
    70    REAL(wp), DIMENSION(jpbdim)     ::   sshbdy            !: Now clim of bdy sea surface height (Flather) 
    71    REAL(wp), DIMENSION(jpbdim)     ::   ubtbdy, vbtbdy    !: Now clim of bdy barotropic velocity components 
    72    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tbdy  , sbdy      !: Now clim of bdy temperature and salinity   
    73    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ubdy  , vbdy    !: Now clim of bdy velocity components 
    74    REAL(wp), DIMENSION(jpbdim) ::   sshtide               !: Tidal boundary array : SSH 
    75    REAL(wp), DIMENSION(jpbdim) ::   utide, vtide          !: Tidal boundary array : U and V 
    76 #if defined key_lim2 
    77    REAL(wp), DIMENSION(jpbdim) ::   frld_bdy    !: now ice leads fraction climatology    
    78    REAL(wp), DIMENSION(jpbdim) ::   hicif_bdy   !: Now ice  thickness climatology 
    79    REAL(wp), DIMENSION(jpbdim) ::   hsnif_bdy   !: now snow thickness 
    80 #endif 
     101   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions 
     102                                                                          !: =1 => some data to be read in from data files 
     103   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays 
     104   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
     105   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 
    81106 
    82107   !!---------------------------------------------------------------------- 
    … …  
    94119      !!---------------------------------------------------------------------- 
    95120      ! 
    96       ALLOCATE( bdytmask(jpi,jpj) , tbdy(jpbdim,jpk) , sbdy(jpbdim,jpk) ,     & 
    97          &      bdyumask(jpi,jpj) , ubdy(jpbdim,jpk) ,                        & 
    98          &      bdyvmask(jpi,jpj) , vbdy(jpbdim,jpk) ,                    STAT=bdy_oce_alloc ) 
     121      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                    &   
     122         &      STAT=bdy_oce_alloc ) 
    99123         ! 
    100124      IF( lk_mpp             )   CALL mpp_sum ( bdy_oce_alloc ) 
    … …  
    112136   !!====================================================================== 
    113137END MODULE bdy_oce 
     138 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90

    r2528 r3294  
    77   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
    88   !!            3.3  !  2010-09  (D. Storkey and E. O'Dea) update for Shelf configurations 
     9   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    910   !!---------------------------------------------------------------------- 
    1011#if defined   key_bdy 
    … …  
    1617   PUBLIC 
    1718 
     19# if ! defined key_agrif 
    1820   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy  = .TRUE.   !: Unstructured Ocean Boundary Condition flag 
    19    INTEGER, PUBLIC, PARAMETER ::   jpbdta  = 20000    !: Max length of bdy field in file 
    20    INTEGER, PUBLIC, PARAMETER ::   jpbdim  = 20000    !: Max length of bdy field on a processor 
     21# else 
     22   LOGICAL, PUBLIC            ::   lk_bdy  = .TRUE.   !: Unstructured Ocean Boundary Condition flag 
     23# endif 
     24   INTEGER, PUBLIC, PARAMETER ::   jp_bdy  = 10       !: Maximum number of bdy sets 
    2125   INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000     !: Max number of time dumps per file 
    22    INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 6        !: Number of horizontal grid types used  (T, u, v, f) 
     26   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3        !: Number of horizontal grid types used  (T, U, V) 
     27 
     28   !! Flags for choice of schemes 
     29   INTEGER, PUBLIC, PARAMETER ::   jp_none         = 0        !: Flag for no open boundary condition 
     30   INTEGER, PUBLIC, PARAMETER ::   jp_frs          = 1        !: Flag for Flow Relaxation Scheme 
     31   INTEGER, PUBLIC, PARAMETER ::   jp_flather      = 2        !: Flag for Flather 
    2332#else 
    2433   !!---------------------------------------------------------------------- 
    2534   !!   Default option :            NO Unstructured open boundary condition 
    2635   !!---------------------------------------------------------------------- 
    27    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .FALSE.   !: Unstructured Ocean Boundary Condition flag 
     36# if ! defined key_agrif 
     37   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy  = .FALSE.   !: Unstructured Ocean Boundary Condition flag 
     38# else 
     39   LOGICAL, PUBLIC            ::   lk_bdy  = .FALSE.   !: Unstructured Ocean Boundary Condition flag 
     40# endif 
    2841#endif 
    2942 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r2715 r3294  
    1010   !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
     12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_bdy 
    1415   !!---------------------------------------------------------------------- 
    15    !!   'key_bdy'                     Unstructured Open Boundary Conditions 
    16    !!---------------------------------------------------------------------- 
    17    !!   bdy_dta_frs    : read u, v, t, s data along open boundaries 
    18    !!   bdy_dta_fla : read depth-mean velocities and elevation along open boundaries         
    19    !!---------------------------------------------------------------------- 
     16   !!   'key_bdy'                     Open Boundary Conditions 
     17   !!---------------------------------------------------------------------- 
     18   !!    bdy_dta        : read external data along open boundaries from file 
     19   !!    bdy_dta_init   : initialise arrays etc for reading of external data 
     20   !!---------------------------------------------------------------------- 
     21   USE wrk_nemo        ! Memory Allocation 
     22   USE timing          ! Timing 
    2023   USE oce             ! ocean dynamics and tracers 
    2124   USE dom_oce         ! ocean space and time domain 
    2225   USE phycst          ! physical constants 
    23    USE bdy_oce         ! ocean open boundary conditions 
     26   USE bdy_oce         ! ocean open boundary conditions   
    2427   USE bdytides        ! tidal forcing at boundaries 
    25    USE iom 
    26    USE ioipsl 
     28   USE fldread         ! read input fields 
     29   USE iom             ! IOM library 
    2730   USE in_out_manager  ! I/O logical units 
    2831#if defined key_lim2 
    … …  
    3336   PRIVATE 
    3437 
    35    PUBLIC   bdy_dta_frs      ! routines called by step.F90 
    36    PUBLIC   bdy_dta_fla  
    37    PUBLIC   bdy_dta_alloc    ! routine called by bdy_init.F90 
    38  
    39    INTEGER ::   numbdyt, numbdyu, numbdyv                      ! logical units for T-, U-, & V-points data file, resp. 
    40    INTEGER ::   ntimes_bdy                                     ! exact number of time dumps in data files 
    41    INTEGER ::   nbdy_b, nbdy_a                                 ! record of bdy data file for before and after time step 
    42    INTEGER ::   numbdyt_bt, numbdyu_bt, numbdyv_bt             ! logical unit for T-, U- & V-points data file, resp. 
    43    INTEGER ::   ntimes_bdy_bt                                  ! exact number of time dumps in data files 
    44    INTEGER ::   nbdy_b_bt, nbdy_a_bt                           ! record of bdy data file for before and after time step 
    45  
    46    INTEGER, DIMENSION (jpbtime) ::   istep, istep_bt           ! time array in seconds in each data file 
    47  
    48    REAL(wp) ::  zoffset                                        ! time offset between time origin in file & start time of model run 
    49  
    50    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tbdydta, sbdydta   ! time interpolated values of T and S bdy data    
    51    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ubdydta, vbdydta   ! time interpolated values of U and V bdy data  
    52    REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta ! Arrays used for time interpolation of bdy data    
    53    REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            ! bdy data of ssh 
    54  
    55 #if defined key_lim2 
    56    REAL(wp), DIMENSION(jpbdim,2)     ::   frld_bdydta          ! } 
    57    REAL(wp), DIMENSION(jpbdim,2)     ::   hicif_bdydta         ! } Arrays used for time interp. of ice bdy data  
    58    REAL(wp), DIMENSION(jpbdim,2)     ::   hsnif_bdydta         ! } 
    59 #endif 
    60  
     38   PUBLIC   bdy_dta          ! routine called by step.F90 and dynspg_ts.F90 
     39   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90 
     40 
     41   INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set. 
     42   INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets. 
     43 
     44   LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions 
     45                                                               ! =F => baroclinic velocities in 3D boundary conditions 
     46 
     47   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read) 
     48 
     49   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap 
     50 
     51#  include "domzgr_substitute.h90" 
    6152   !!---------------------------------------------------------------------- 
    6253   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    6657CONTAINS 
    6758 
    68   FUNCTION bdy_dta_alloc() 
    69      !!---------------------------------------------------------------------- 
    70      USE lib_mpp, ONLY: ctl_warn, mpp_sum 
    71      ! 
    72      INTEGER :: bdy_dta_alloc 
    73      !!---------------------------------------------------------------------- 
    74      ! 
    75      ALLOCATE(tbdydta(jpbdim,jpk,2), sbdydta(jpbdim,jpk,2), & 
    76               ubdydta(jpbdim,jpk,2), vbdydta(jpbdim,jpk,2), Stat=bdy_dta_alloc) 
    77  
    78      IF( lk_mpp           ) CALL mpp_sum ( bdy_dta_alloc ) 
    79      IF(bdy_dta_alloc /= 0) CALL ctl_warn('bdy_dta_alloc: failed to allocate arrays') 
    80  
    81    END FUNCTION bdy_dta_alloc 
    82  
    83  
    84    SUBROUTINE bdy_dta_frs( kt ) 
     59      SUBROUTINE bdy_dta( kt, jit, time_offset ) 
    8560      !!---------------------------------------------------------------------- 
    86       !!                   ***  SUBROUTINE bdy_dta_frs  *** 
     61      !!                   ***  SUBROUTINE bdy_dta  *** 
    8762      !!                     
    88       !! ** Purpose :   Read unstructured boundary data for FRS condition. 
    89       !! 
    90       !! ** Method  :   At the first timestep, read in boundary data for two 
    91       !!                times from the file and time-interpolate. At other  
    92       !!                timesteps, check to see if we need another time from  
    93       !!                the file. If so read it in. Time interpolate. 
     63      !! ** Purpose :   Update external data for open boundary conditions 
     64      !! 
     65      !! ** Method  :   Use fldread.F90 
     66      !!                 
    9467      !!---------------------------------------------------------------------- 
    95       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index (for timesplitting option, otherwise zero) 
    96       !! 
    97       CHARACTER(LEN=80), DIMENSION(3) ::   clfile               ! names of input files 
    98       CHARACTER(LEN=70 )              ::   clunits              ! units attribute of time coordinate 
    99       LOGICAL ::   lect                                         ! flag for reading 
    100       INTEGER ::   it, ib, ik, igrd                             ! dummy loop indices 
    101       INTEGER ::   igrd_start, igrd_end                         ! start and end of loops on igrd 
    102       INTEGER ::   idvar                                        ! netcdf var ID 
    103       INTEGER ::   iman, i15, imois                             ! Time variables for monthly clim forcing 
    104       INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 
    105       INTEGER ::   itimer, totime 
    106       INTEGER ::   ii, ij                                       ! array addresses 
    107       INTEGER ::   ipi, ipj, ipk, inum                          ! local integers (NetCDF read) 
    108       INTEGER ::   iyear0, imonth0, iday0 
    109       INTEGER ::   ihours0, iminutes0, isec0 
    110       INTEGER ::   iyear, imonth, iday, isecs 
    111       INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files 
    112       REAL(wp) ::   dayfrac, zxy, zoffsett 
    113       REAL(wp) ::   zoffsetu, zoffsetv 
    114       REAL(wp) ::   dayjul0, zdayjulini 
    115       REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files 
    116       REAL(wp), DIMENSION(jpbdta,1,jpk) ::   zdta               ! temporary array for data fields 
     68      !! 
     69      INTEGER, INTENT( in )           ::   kt    ! ocean time-step index  
     70      INTEGER, INTENT( in ), OPTIONAL ::   jit   ! subcycle time-step index (for timesplitting option) 
     71      INTEGER, INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     72                                                        ! is present then units = subcycle timesteps. 
     73                                                        ! time_offset = 0 => get data at "now" time level 
     74                                                        ! time_offset = -1 => get data at "before" time level 
     75                                                        ! time_offset = +1 => get data at "after" time level 
     76                                                        ! etc. 
     77      !! 
     78      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices 
     79      INTEGER,          DIMENSION(jpbgrd) ::   ilen1  
     80      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts 
     81      !! 
    11782      !!--------------------------------------------------------------------------- 
    118  
    119  
    120       IF( ln_dyn_frs .OR. ln_tra_frs    & 
    121          &               .OR. ln_ice_frs ) THEN  ! If these are both false then this routine does nothing 
    122  
    123       ! -------------------- ! 
    124       !    Initialization    ! 
    125       ! -------------------- ! 
    126  
    127       lect   = .false.           ! If true, read a time record 
    128  
    129       ! Some time variables for monthly climatological forcing: 
    130       ! ******************************************************* 
    131  
    132 !!gm  here  use directely daymod calendar variables 
    133   
    134       iman = INT( raamo )      ! Number of months in a year 
    135  
    136       i15 = INT( 2*REAL( nday, wp ) / ( REAL( nmonth_len(nmonth), wp ) + 0.5 ) ) 
    137       ! i15=0 if the current day is in the first half of the month, else i15=1 
    138  
    139       imois = nmonth + i15 - 1            ! imois is the first month record 
    140       IF( imois == 0 )   imois = iman 
    141  
    142       ! Time variable for non-climatological forcing: 
    143       ! ********************************************* 
    144       itimer = (kt-nit000+1)*rdt      ! current time in seconds for interpolation  
    145  
    146  
    147       !                                                !-------------------! 
    148       IF( kt == nit000 ) THEN                          !  First call only  ! 
    149          !                                             !-------------------! 
    150          istep(:) = 0 
    151          nbdy_b   = 0 
    152          nbdy_a   = 0 
    153  
    154          ! Get time information from bdy data file 
    155          ! *************************************** 
    156  
    157          IF(lwp) WRITE(numout,*) 
    158          IF(lwp) WRITE(numout,*)    'bdy_dta_frs : Initialize unstructured boundary data' 
    159          IF(lwp) WRITE(numout,*)    '~~~~~~~'  
    160  
    161          IF     ( nn_dtactl == 0 ) THEN 
    162             ! 
    163             IF(lwp) WRITE(numout,*) '          Bdy data are taken from initial conditions' 
    164             ! 
    165          ELSEIF (nn_dtactl == 1) THEN 
    166             ! 
    167             IF(lwp) WRITE(numout,*) '          Bdy data are read in netcdf files' 
    168             ! 
    169             dayfrac = adatrj  - REAL( itimer, wp ) / 86400.   ! day fraction at time step kt-1 
    170             dayfrac = dayfrac - INT ( dayfrac )               ! 
    171             totime  = ( nitend - nit000 + 1 ) * rdt           ! Total time of the run to verify that all the 
    172             !                                                 ! necessary time dumps in file are included 
    173             ! 
    174             clfile(1) = cn_dta_frs_T 
    175             clfile(2) = cn_dta_frs_U 
    176             clfile(3) = cn_dta_frs_V 
    177             !                                                   
    178             ! how many files are we to read in? 
    179             igrd_start = 1 
    180             igrd_end   = 3 
    181             IF(.NOT. ln_tra_frs .AND. .NOT. ln_ice_frs) THEN       ! No T-grid file. 
    182                igrd_start = 2 
    183             ELSEIF ( .NOT. ln_dyn_frs ) THEN                           ! No U-grid or V-grid file. 
    184                igrd_end   = 1          
    185             ENDIF 
    186  
    187             DO igrd = igrd_start, igrd_end                     !  loop over T, U & V grid  ! 
    188                !                                               !---------------------------! 
    189                CALL iom_open( clfile(igrd), inum ) 
    190                CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits )  
    191  
    192                SELECT CASE( igrd ) 
    193                   CASE (1)   ;   numbdyt = inum 
    194                   CASE (2)   ;   numbdyu = inum 
    195                   CASE (3)   ;   numbdyv = inum 
    196                END SELECT 
    197  
    198                ! Calculate time offset  
    199                READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 
    200                ! Convert time origin in file to julian days  
    201                isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 
    202                CALL ymds2ju(iyear0, imonth0, iday0, REAL(isec0, wp), dayjul0) 
    203                ! Compute model initialization time  
    204                iyear  = ndastp / 10000 
    205                imonth = ( ndastp - iyear * 10000 ) / 100 
    206                iday   = ndastp - iyear * 10000 - imonth * 100 
    207                isecs  = dayfrac * 86400 
    208                CALL ymds2ju(iyear, imonth, iday, REAL(isecs, wp) , zdayjulini) 
    209                ! offset from initialization date: 
    210                zoffset = (dayjul0-zdayjulini)*86400 
    211                ! 
    212 7000           FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 
    213  
    214                !! TO BE DONE... Check consistency between calendar from file  
    215                !! (available optionally from iom_gettime) and calendar in model  
    216                !! when calendar in model available outside of IOIPSL. 
    217  
    218                IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy 
    219                IF(lwp) WRITE(numout,*) 'offset: ',zoffset 
    220                IF(lwp) WRITE(numout,*) 'totime: ',totime 
    221                IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr(1:ntimes_bdy) 
    222  
    223                ! Check that there are not too many times in the file.  
    224                IF( ntimes_bdy > jpbtime ) THEN 
    225                   WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy 
    226                   CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 ) 
    227                ENDIF 
    228  
    229                ! Check that time array increases: 
    230                it = 1 
    231                DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 )  
    232                   it = it + 1 
    233                END DO 
    234                ! 
    235                IF( it /= ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN 
    236                      WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
    237                      CALL ctl_stop( 'Time array in unstructured boundary data files',   & 
    238                         &           'does not continuously increase.'               , ctmp1 ) 
    239                ENDIF 
    240                ! 
    241                ! Check that times in file span model run time: 
    242                IF( zstepr(1) + zoffset > 0 ) THEN 
    243                      WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
    244                      CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 ) 
    245                END IF 
    246                IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN 
    247                      WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
    248                      CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 ) 
    249                END IF 
    250                ! 
    251                SELECT CASE( igrd ) 
    252                   CASE (1) 
    253                     ntimes_bdyt = ntimes_bdy 
    254                     zoffsett = zoffset 
    255                     istept(:) = INT( zstepr(:) + zoffset ) 
    256                     numbdyt = inum 
    257                   CASE (2) 
    258                     ntimes_bdyu = ntimes_bdy 
    259                     zoffsetu = zoffset 
    260                     istepu(:) = INT( zstepr(:) + zoffset ) 
    261                     numbdyu = inum 
    262                   CASE (3) 
    263                     ntimes_bdyv = ntimes_bdy 
    264                     zoffsetv = zoffset 
    265                     istepv(:) = INT( zstepr(:) + zoffset ) 
    266                     numbdyv = inum 
    267                END SELECT 
    268                ! 
    269             END DO                                         ! end loop over T, U & V grid  
    270  
    271             IF (igrd_start == 1 .and. igrd_end == 3) THEN 
    272                ! Only test differences if we are reading in 3 files 
    273                ! Verify time consistency between files   
    274                IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 
    275                   CALL ctl_stop( 'Bdy data files must have the same number of time dumps',   & 
    276                   &           'Multiple time frequencies not implemented yet'  ) 
    277                ENDIF 
    278                ntimes_bdy = ntimes_bdyt 
    279                ! 
    280                IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN 
    281                   CALL ctl_stop( 'Bdy data files must have the same time origin',   & 
    282                   &           'Multiple time frequencies not implemented yet' ) 
    283                ENDIF 
    284                zoffset = zoffsett 
    285             ENDIF 
    286  
    287             IF( igrd_start == 1 ) THEN   ;   istep(:) = istept(:) 
    288             ELSE                         ;   istep(:) = istepu(:) 
    289             ENDIF 
    290  
    291             ! Check number of time dumps:               
    292             IF( ntimes_bdy == 1 .AND. .NOT. ln_clim ) THEN 
    293               CALL ctl_stop( 'There is only one time dump in data files',   & 
    294                  &           'Choose ln_clim=.true. in namelist for constant bdy forcing.' ) 
    295             ENDIF 
    296  
    297             IF( ln_clim ) THEN 
    298               IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN 
    299                  CALL ctl_stop( 'For climatological boundary forcing (ln_clim=.true.),',   & 
    300                     &           'bdy data files must contain 1 or 12 time dumps.' ) 
    301               ELSEIF( ntimes_bdy ==  1 ) THEN 
    302                 IF(lwp) WRITE(numout,*) 
    303                 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 
    304               ELSEIF( ntimes_bdy == 12 ) THEN 
    305                 IF(lwp) WRITE(numout,*) 
    306                 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 
    307               ENDIF 
    308             ENDIF 
    309  
    310             ! Find index of first record to read (before first model time).  
    311             it = 1 
    312             DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 ) 
    313                it = it + 1 
    314             END DO 
    315             nbdy_b = it 
    316             ! 
    317             IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset 
    318             IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b 
    319  
    320          ENDIF ! endif (nn_dtactl == 1) 
    321  
    322  
    323          ! 1.2  Read first record in file if necessary (ie if nn_dtactl == 1) 
    324          ! ***************************************************************** 
    325  
    326          IF( nn_dtactl == 0 ) THEN      ! boundary data arrays are filled with initial conditions 
    327             ! 
    328             IF (ln_tra_frs) THEN 
    329                igrd = 1            ! T-points data  
    330                DO ib = 1, nblen(igrd) 
    331                   ii = nbi(ib,igrd) 
    332                   ij = nbj(ib,igrd) 
     83      !! 
     84      IF( nn_timing == 1 ) CALL timing_start('bdy_dta') 
     85 
     86      ! Initialise data arrays once for all from initial conditions where required 
     87      !--------------------------------------------------------------------------- 
     88      IF( kt .eq. nit000 .and. .not. PRESENT(jit) ) THEN 
     89 
     90         ! Calculate depth-mean currents 
     91         !----------------------------- 
     92         CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
     93 
     94         pu2d(:,:) = 0.e0 
     95         pv2d(:,:) = 0.e0 
     96 
     97         DO ik = 1, jpkm1   !! Vertically integrated momentum trends 
     98             pu2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik) 
     99             pv2d(:,:) = pv2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik) 
     100         END DO 
     101         pu2d(:,:) = pu2d(:,:) * hur(:,:) 
     102         pv2d(:,:) = pv2d(:,:) * hvr(:,:) 
     103          
     104         DO ib_bdy = 1, nb_bdy 
     105 
     106            nblen => idx_bdy(ib_bdy)%nblen 
     107            nblenrim => idx_bdy(ib_bdy)%nblenrim 
     108 
     109            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN  
     110               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     111                  ilen1(:) = nblen(:) 
     112               ELSE 
     113                  ilen1(:) = nblenrim(:) 
     114               ENDIF 
     115               igrd = 1 
     116               DO ib = 1, ilen1(igrd) 
     117                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     118                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     119                  dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     120               END DO  
     121               igrd = 2 
     122               DO ib = 1, ilen1(igrd) 
     123                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     124                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     125                  dta_bdy(ib_bdy)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1)          
     126               END DO  
     127               igrd = 3 
     128               DO ib = 1, ilen1(igrd) 
     129                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     130                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     131                  dta_bdy(ib_bdy)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1)          
     132               END DO  
     133            ENDIF 
     134 
     135            IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN  
     136               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     137                  ilen1(:) = nblen(:) 
     138               ELSE 
     139                  ilen1(:) = nblenrim(:) 
     140               ENDIF 
     141               igrd = 2  
     142               DO ib = 1, ilen1(igrd) 
    333143                  DO ik = 1, jpkm1 
    334                      tbdy(ib,ik) = tn(ii,ij,ik) 
    335                      sbdy(ib,ik) = sn(ii,ij,ik) 
     144                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     145                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     146                     dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)          
     147                  END DO 
     148               END DO  
     149               igrd = 3  
     150               DO ib = 1, ilen1(igrd) 
     151                  DO ik = 1, jpkm1 
     152                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     153                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     154                     dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik)          
     155                     END DO 
     156               END DO  
     157            ENDIF 
     158 
     159            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN  
     160               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     161                  ilen1(:) = nblen(:) 
     162               ELSE 
     163                  ilen1(:) = nblenrim(:) 
     164               ENDIF 
     165               igrd = 1                       ! Everything is at T-points here 
     166               DO ib = 1, ilen1(igrd) 
     167                  DO ik = 1, jpkm1 
     168                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     169                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     170                     dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
     171                     dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
     172                  END DO 
     173               END DO  
     174            ENDIF 
     175 
     176#if defined key_lim2 
     177            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
     178               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     179                  ilen1(:) = nblen(:) 
     180               ELSE 
     181                  ilen1(:) = nblenrim(:) 
     182               ENDIF 
     183               igrd = 1                       ! Everything is at T-points here 
     184               DO ib = 1, ilen1(igrd) 
     185                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     186                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     187                  dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)          
     188                  dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)          
     189                  dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)          
     190               END DO  
     191            ENDIF 
     192#endif 
     193 
     194         ENDDO ! ib_bdy 
     195 
     196         CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     197 
     198      ENDIF ! kt .eq. nit000 
     199 
     200      ! update external data from files 
     201      !-------------------------------- 
     202      
     203      jstart = 1 
     204      DO ib_bdy = 1, nb_bdy    
     205         IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required 
     206       
     207            IF( PRESENT(jit) ) THEN 
     208               ! Update barotropic boundary conditions only 
     209               ! jit is optional argument for fld_read and tide_update 
     210               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     211                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     212                     dta_bdy(ib_bdy)%ssh(:) = 0.0 
     213                     dta_bdy(ib_bdy)%u2d(:) = 0.0 
     214                     dta_bdy(ib_bdy)%v2d(:) = 0.0 
     215                  ENDIF 
     216                  IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data 
     217                     jend = jstart + 2 
     218                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),   & 
     219                     &              jit=jit, time_offset=time_offset ) 
     220                  ENDIF 
     221                  IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     222                     CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
     223                     &                 jit=jit, time_offset=time_offset ) 
     224                  ENDIF 
     225               ENDIF 
     226            ELSE 
     227               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     228                  dta_bdy(ib_bdy)%ssh(:) = 0.0 
     229                  dta_bdy(ib_bdy)%u2d(:) = 0.0 
     230                  dta_bdy(ib_bdy)%v2d(:) = 0.0 
     231               ENDIF 
     232               IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
     233                  jend = jstart + nb_bdy_fld(ib_bdy) - 1 
     234                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 
     235               ENDIF 
     236               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     237                  CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
     238               ENDIF 
     239            ENDIF 
     240            jstart = jend+1 
     241 
     242            ! If full velocities in boundary data then split into barotropic and baroclinic data 
     243            ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same 
     244            ! time as the dynspg_ts option).  
     245 
     246            IF( ln_full_vel_array(ib_bdy) .and.                                             &  
     247           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN  
     248 
     249               igrd = 2                      ! zonal velocity 
     250               dta_bdy(ib_bdy)%u2d(:) = 0.0 
     251               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     252                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     253                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     254                  DO ik = 1, jpkm1 
     255                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
     256              &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     257                  END DO 
     258                  dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     259                  DO ik = 1, jpkm1 
     260                     dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)  
    336261                  END DO 
    337262               END DO 
    338             ENDIF 
    339  
    340             IF(ln_dyn_frs) THEN 
    341                igrd = 2            ! U-points data  
    342                DO ib = 1, nblen(igrd) 
    343                   ii = nbi(ib,igrd) 
    344                   ij = nbj(ib,igrd) 
     263 
     264               igrd = 3                      ! meridional velocity 
     265               dta_bdy(ib_bdy)%v2d(:) = 0.0 
     266               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     267                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     268                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    345269                  DO ik = 1, jpkm1 
    346                      ubdy(ib,ik) = un(ii, ij, ik) 
     270                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
     271              &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     272                  END DO 
     273                  dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     274                  DO ik = 1, jpkm1 
     275                     dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)  
    347276                  END DO 
    348277               END DO 
    349                ! 
    350                igrd = 3            ! V-points data  
    351                DO ib = 1, nblen(igrd)             
    352                   ii = nbi(ib,igrd) 
    353                   ij = nbj(ib,igrd) 
    354                   DO ik = 1, jpkm1 
    355                      vbdy(ib,ik) = vn(ii, ij, ik) 
    356                   END DO 
    357                END DO 
    358             ENDIF 
    359             ! 
    360 #if defined key_lim2 
    361             IF( ln_ice_frs ) THEN 
    362                igrd = 1            ! T-points data 
    363                DO ib = 1, nblen(igrd) 
    364                   frld_bdy (ib) =  frld(nbi(ib,igrd), nbj(ib,igrd)) 
    365                   hicif_bdy(ib) = hicif(nbi(ib,igrd), nbj(ib,igrd)) 
    366                   hsnif_bdy(ib) = hsnif(nbi(ib,igrd), nbj(ib,igrd)) 
    367                END DO 
    368             ENDIF 
    369 #endif 
    370          ELSEIF( nn_dtactl == 1 ) THEN    ! Set first record in the climatological case:    
    371             ! 
    372             IF( ln_clim .AND. ntimes_bdy == 1 ) THEN 
    373                nbdy_a = 1 
    374             ELSEIF( ln_clim .AND. ntimes_bdy == iman ) THEN 
    375                nbdy_b = 0 
    376                nbdy_a = imois 
     278     
     279            ENDIF 
     280 
     281         END IF ! nn_dta(ib_bdy) = 1 
     282      END DO  ! ib_bdy 
     283 
     284      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta') 
     285 
     286      END SUBROUTINE bdy_dta 
     287 
     288 
     289      SUBROUTINE bdy_dta_init 
     290      !!---------------------------------------------------------------------- 
     291      !!                   ***  SUBROUTINE bdy_dta_init  *** 
     292      !!                     
     293      !! ** Purpose :   Initialise arrays for reading of external data  
     294      !!                for open boundary conditions 
     295      !! 
     296      !! ** Method  :   Use fldread.F90 
     297      !!                 
     298      !!---------------------------------------------------------------------- 
     299      USE dynspg_oce, ONLY: lk_dynspg_ts 
     300      !! 
     301      INTEGER     ::  ib_bdy, jfld, jstart, jend, ierror  ! local indices 
     302      !! 
     303      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files 
     304      CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files 
     305      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data 
     306                                                                ! =F => baroclinic velocities in 3D boundary data 
     307      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays 
     308      INTEGER,              DIMENSION(jpbgrd) ::  ilen0         ! size of local arrays 
     309      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays 
     310      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld 
     311      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V) 
     312      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts 
     313      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures 
     314      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !  
     315      TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
     316#if defined key_lim2 
     317      TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      ! 
     318#endif 
     319      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
     320#if defined key_lim2 
     321      NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 
     322#endif 
     323      NAMELIST/nambdy_dta/ ln_full_vel 
     324      !!--------------------------------------------------------------------------- 
     325 
     326      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init') 
     327 
     328      ! Set nn_dta 
     329      DO ib_bdy = 1, nb_bdy 
     330         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       & 
     331                               ,nn_dyn3d_dta(ib_bdy)       & 
     332                               ,nn_tra_dta(ib_bdy)         & 
     333#if defined key_lim2 
     334                               ,nn_ice_lim2_dta(ib_bdy)    & 
     335#endif 
     336                              ) 
     337         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1 
     338      END DO 
     339 
     340      ! Work out upper bound of how many fields there are to read in and allocate arrays 
     341      ! --------------------------------------------------------------------------- 
     342      ALLOCATE( nb_bdy_fld(nb_bdy) ) 
     343      nb_bdy_fld(:) = 0 
     344      DO ib_bdy = 1, nb_bdy          
     345         IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
     346            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
     347         ENDIF 
     348         IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 
     349            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
     350         ENDIF 
     351         IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN 
     352            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
     353         ENDIF 
     354#if defined key_lim2 
     355         IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN 
     356            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
     357         ENDIF 
     358#endif                
     359      ENDDO             
     360 
     361      nb_bdy_fld_sum = SUM( nb_bdy_fld ) 
     362 
     363      ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror ) 
     364      IF( ierror > 0 ) THEN    
     365         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN   
     366      ENDIF 
     367      ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror ) 
     368      IF( ierror > 0 ) THEN    
     369         CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' )   ;   RETURN   
     370      ENDIF 
     371      ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror ) 
     372      IF( ierror > 0 ) THEN    
     373         CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' )   ;   RETURN   
     374      ENDIF 
     375      ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) )  
     376      ALLOCATE( ibdy(nb_bdy_fld_sum) )  
     377      ALLOCATE( igrid(nb_bdy_fld_sum) )  
     378 
     379      ! Read namelists 
     380      ! -------------- 
     381      REWIND(numnam) 
     382      jfld = 0  
     383      DO ib_bdy = 1, nb_bdy          
     384         IF( nn_dta(ib_bdy) .eq. 1 ) THEN 
     385            ! set file information 
     386            cn_dir = './'        ! directory in which the model is executed 
     387            ln_full_vel = .false. 
     388            ! ... default values (NB: frequency positive => hours, negative => months) 
     389            !                    !  file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation  ! 
     390            !                    !  name       ! hours !   name     !  (T/F)  !  (T/F)  !  'monthly'  ! filename ! pairs     ! 
     391            bn_ssh     = FLD_N(  'bdy_ssh'     ,  24   , 'sossheig' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     392            bn_u2d     = FLD_N(  'bdy_vel2d_u' ,  24   , 'vobtcrtx' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     393            bn_v2d     = FLD_N(  'bdy_vel2d_v' ,  24   , 'vobtcrty' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     394            bn_u3d     = FLD_N(  'bdy_vel3d_u' ,  24   , 'vozocrtx' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     395            bn_v3d     = FLD_N(  'bdy_vel3d_v' ,  24   , 'vomecrty' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     396            bn_tem     = FLD_N(  'bdy_tem'     ,  24   , 'votemper' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     397            bn_sal     = FLD_N(  'bdy_sal'     ,  24   , 'vosaline' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     398#if defined key_lim2 
     399            bn_frld    = FLD_N(  'bdy_frld'    ,  24   , 'ildsconc' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     400            bn_hicif   = FLD_N(  'bdy_hicif'   ,  24   , 'iicethic' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     401            bn_hsnif   = FLD_N(  'bdy_hsnif'   ,  24   , 'isnothic' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     402#endif 
     403 
     404            ! Important NOT to rewind here. 
     405            READ( numnam, nambdy_dta ) 
     406 
     407            cn_dir_array(ib_bdy) = cn_dir 
     408            ln_full_vel_array(ib_bdy) = ln_full_vel 
     409 
     410            IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN 
     411               CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',& 
     412            &                  'with dynspg_ts option' )   ;   RETURN   
     413            ENDIF              
     414 
     415            nblen => idx_bdy(ib_bdy)%nblen 
     416            nblenrim => idx_bdy(ib_bdy)%nblenrim 
     417 
     418            ! Only read in necessary fields for this set. 
     419            ! Important that barotropic variables come first. 
     420            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
     421 
     422               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
     423                  jfld = jfld + 1 
     424                  blf_i(jfld) = bn_ssh 
     425                  ibdy(jfld) = ib_bdy 
     426                  igrid(jfld) = 1 
     427                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     428                  ilen3(jfld) = 1 
     429               ENDIF 
     430 
     431               IF( .not. ln_full_vel_array(ib_bdy) ) THEN 
     432 
     433                  jfld = jfld + 1 
     434                  blf_i(jfld) = bn_u2d 
     435                  ibdy(jfld) = ib_bdy 
     436                  igrid(jfld) = 2 
     437                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     438                     ilen1(jfld) = nblen(igrid(jfld)) 
     439                  ELSE 
     440                     ilen1(jfld) = nblenrim(igrid(jfld)) 
     441                  ENDIF 
     442                  ilen3(jfld) = 1 
     443 
     444                  jfld = jfld + 1 
     445                  blf_i(jfld) = bn_v2d 
     446                  ibdy(jfld) = ib_bdy 
     447                  igrid(jfld) = 3 
     448                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     449                     ilen1(jfld) = nblen(igrid(jfld)) 
     450                  ELSE 
     451                     ilen1(jfld) = nblenrim(igrid(jfld)) 
     452                  ENDIF 
     453                  ilen3(jfld) = 1 
     454 
     455               ENDIF 
     456 
     457            ENDIF 
     458 
     459            ! baroclinic velocities 
     460            IF( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) .or. & 
     461           &      ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.  & 
     462           &        ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
     463 
     464               jfld = jfld + 1 
     465               blf_i(jfld) = bn_u3d 
     466               ibdy(jfld) = ib_bdy 
     467               igrid(jfld) = 2 
     468               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     469                  ilen1(jfld) = nblen(igrid(jfld)) 
     470               ELSE 
     471                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     472               ENDIF 
     473               ilen3(jfld) = jpk 
     474 
     475               jfld = jfld + 1 
     476               blf_i(jfld) = bn_v3d 
     477               ibdy(jfld) = ib_bdy 
     478               igrid(jfld) = 3 
     479               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     480                  ilen1(jfld) = nblen(igrid(jfld)) 
     481               ELSE 
     482                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     483               ENDIF 
     484               ilen3(jfld) = jpk 
     485 
     486            ENDIF 
     487 
     488            ! temperature and salinity 
     489            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN 
     490 
     491               jfld = jfld + 1 
     492               blf_i(jfld) = bn_tem 
     493               ibdy(jfld) = ib_bdy 
     494               igrid(jfld) = 1 
     495               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     496                  ilen1(jfld) = nblen(igrid(jfld)) 
     497               ELSE 
     498                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     499               ENDIF 
     500               ilen3(jfld) = jpk 
     501 
     502               jfld = jfld + 1 
     503               blf_i(jfld) = bn_sal 
     504               ibdy(jfld) = ib_bdy 
     505               igrid(jfld) = 1 
     506               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     507                  ilen1(jfld) = nblen(igrid(jfld)) 
     508               ELSE 
     509                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     510               ENDIF 
     511               ilen3(jfld) = jpk 
     512 
     513            ENDIF 
     514 
     515#if defined key_lim2 
     516            ! sea ice 
     517            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
     518 
     519               jfld = jfld + 1 
     520               blf_i(jfld) = bn_frld 
     521               ibdy(jfld) = ib_bdy 
     522               igrid(jfld) = 1 
     523               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     524                  ilen1(jfld) = nblen(igrid(jfld)) 
     525               ELSE 
     526                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     527               ENDIF 
     528               ilen3(jfld) = 1 
     529 
     530               jfld = jfld + 1 
     531               blf_i(jfld) = bn_hicif 
     532               ibdy(jfld) = ib_bdy 
     533               igrid(jfld) = 1 
     534               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     535                  ilen1(jfld) = nblen(igrid(jfld)) 
     536               ELSE 
     537                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     538               ENDIF 
     539               ilen3(jfld) = 1 
     540 
     541               jfld = jfld + 1 
     542               blf_i(jfld) = bn_hsnif 
     543               ibdy(jfld) = ib_bdy 
     544               igrid(jfld) = 1 
     545               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     546                  ilen1(jfld) = nblen(igrid(jfld)) 
     547               ELSE 
     548                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     549               ENDIF 
     550               ilen3(jfld) = 1 
     551 
     552            ENDIF 
     553#endif 
     554            ! Recalculate field counts 
     555            !------------------------- 
     556            nb_bdy_fld_sum = 0 
     557            IF( ib_bdy .eq. 1 ) THEN  
     558               nb_bdy_fld(ib_bdy) = jfld 
     559               nb_bdy_fld_sum     = jfld               
    377560            ELSE 
    378                nbdy_a = nbdy_b 
    379             ENDIF 
    380     
    381             ! Read first record: 
    382             ipj  = 1 
    383             ipk  = jpk 
    384             igrd = 1 
    385             ipi  = nblendta(igrd) 
    386  
    387             IF(ln_tra_frs) THEN 
    388                ! 
    389                igrd = 1                                           ! Temperature 
    390                IF( nblendta(igrd) <=  0 ) THEN  
    391                   idvar = iom_varid( numbdyt, 'votemper' ) 
    392                   nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 
    393                ENDIF 
    394                IF(lwp) WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd) 
    395                ipi = nblendta(igrd) 
    396                CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
    397                ! 
    398                DO ib = 1, nblen(igrd) 
    399                   DO ik = 1, jpkm1 
    400                      tbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
    401                   END DO 
    402                END DO 
    403                ! 
    404                igrd = 1                                           ! salinity 
    405                IF( nblendta(igrd) .le. 0 ) THEN  
    406                   idvar = iom_varid( numbdyt, 'vosaline' ) 
    407                   nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 
    408                ENDIF 
    409                IF(lwp) WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd) 
    410                ipi = nblendta(igrd) 
    411                CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
    412                ! 
    413                DO ib = 1, nblen(igrd) 
    414                   DO ik = 1, jpkm1 
    415                      sbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
    416                   END DO 
    417                END DO 
    418             ENDIF  ! ln_tra_frs 
    419   
    420             IF( ln_dyn_frs ) THEN 
    421                ! 
    422                igrd = 2                                           ! u-velocity 
    423                IF ( nblendta(igrd) .le. 0 ) THEN  
    424                  idvar = iom_varid( numbdyu,'vozocrtx' ) 
    425                  nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar) 
    426                ENDIF 
    427                IF(lwp) WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd) 
    428                ipi = nblendta(igrd) 
    429                CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
    430                DO ib = 1, nblen(igrd) 
    431                   DO ik = 1, jpkm1 
    432                      ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
    433                   END DO 
    434                END DO 
    435                ! 
    436                igrd = 3                                           ! v-velocity 
    437                IF ( nblendta(igrd) .le. 0 ) THEN  
    438                  idvar = iom_varid( numbdyv,'vomecrty' ) 
    439                  nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar) 
    440                ENDIF 
    441                IF(lwp) WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd) 
    442                ipi = nblendta(igrd) 
    443                CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
    444                DO ib = 1, nblen(igrd) 
    445                   DO ik = 1, jpkm1 
    446                      vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
    447                   END DO 
    448                END DO 
    449             ENDIF ! ln_dyn_frs 
    450  
    451 #if defined key_lim2 
    452             IF( ln_ice_frs ) THEN 
    453               ! 
    454               igrd=1                                              ! leads fraction 
    455               IF(lwp) WRITE(numout,*) 'Dim size for ildsconc is ',nblendta(igrd) 
    456               ipi=nblendta(igrd) 
    457               CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a ) 
    458               DO ib=1, nblen(igrd) 
    459                 frld_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1) 
    460               END DO 
    461               ! 
    462               igrd=1                                              ! ice thickness 
    463               IF(lwp) WRITE(numout,*) 'Dim size for iicethic is ',nblendta(igrd) 
    464               ipi=nblendta(igrd) 
    465               CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a ) 
    466               DO ib=1, nblen(igrd) 
    467                 hicif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1) 
    468               END DO 
    469               ! 
    470               igrd=1                                              ! snow thickness 
    471               IF(lwp) WRITE(numout,*) 'Dim size for isnowthi is ',nblendta(igrd) 
    472               ipi=nblendta(igrd) 
    473               CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a ) 
    474               DO ib=1, nblen(igrd) 
    475                 hsnif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1) 
    476               END DO 
    477             ENDIF ! just if ln_ice_frs is set 
    478 #endif 
    479  
    480             IF( .NOT.ln_clim .AND. istep(1) > 0 ) THEN     ! First data time is after start of run 
    481                nbdy_b = nbdy_a                                 ! Put first value in both time levels 
    482                IF( ln_tra_frs ) THEN 
    483                  tbdydta(:,:,1) = tbdydta(:,:,2) 
    484                  sbdydta(:,:,1) = sbdydta(:,:,2) 
    485                ENDIF 
    486                IF( ln_dyn_frs ) THEN 
    487                  ubdydta(:,:,1) = ubdydta(:,:,2) 
    488                  vbdydta(:,:,1) = vbdydta(:,:,2) 
    489                ENDIF 
    490 #if defined key_lim2 
    491                IF( ln_ice_frs ) THEN 
    492                   frld_bdydta (:,1) =  frld_bdydta(:,2) 
    493                   hicif_bdydta(:,1) = hicif_bdydta(:,2) 
    494                   hsnif_bdydta(:,1) = hsnif_bdydta(:,2) 
    495                ENDIF 
    496 #endif 
    497             END IF 
    498             ! 
    499          END IF   ! nn_dtactl == 0/1 
    500   
    501          ! In the case of constant boundary forcing fill bdy arrays once for all 
    502          IF( ln_clim .AND. ntimes_bdy == 1 ) THEN 
    503             IF( ln_tra_frs ) THEN 
    504                tbdy  (:,:) = tbdydta  (:,:,2) 
    505                sbdy  (:,:) = sbdydta  (:,:,2) 
    506             ENDIF 
    507             IF( ln_dyn_frs) THEN 
    508                ubdy  (:,:) = ubdydta  (:,:,2) 
    509                vbdy  (:,:) = vbdydta  (:,:,2) 
    510             ENDIF 
    511 #if defined key_lim2 
    512             IF( ln_ice_frs ) THEN 
    513                frld_bdy (:) = frld_bdydta (:,2) 
    514                hicif_bdy(:) = hicif_bdydta(:,2) 
    515                hsnif_bdy(:) = hsnif_bdydta(:,2) 
    516             ENDIF 
    517 #endif 
    518  
    519             IF( ln_tra_frs .OR. ln_ice_frs) CALL iom_close( numbdyt ) 
    520             IF( ln_dyn_frs                    ) CALL iom_close( numbdyu ) 
    521             IF( ln_dyn_frs                    ) CALL iom_close( numbdyv ) 
    522          END IF 
    523          ! 
    524       ENDIF                                            ! End if nit000 
    525  
    526  
    527       !                                                !---------------------! 
    528       IF( nn_dtactl == 1 .AND. ntimes_bdy > 1 ) THEN    !  at each time step  ! 
    529          !                                             !---------------------! 
    530          ! Read one more record if necessary 
    531          !********************************** 
    532  
    533          IF( ln_clim .AND. imois /= nbdy_b ) THEN      ! remember that nbdy_b=0 for kt=nit000 
    534             nbdy_b = imois 
    535             nbdy_a = imois + 1 
    536             nbdy_b = MOD( nbdy_b, iman )   ;   IF( nbdy_b == 0 ) nbdy_b = iman 
    537             nbdy_a = MOD( nbdy_a, iman )   ;   IF( nbdy_a == 0 ) nbdy_a = iman 
    538             lect=.true. 
    539          ELSEIF( .NOT.ln_clim .AND. itimer >= istep(nbdy_a) ) THEN 
    540  
    541             IF( nbdy_a < ntimes_bdy ) THEN 
    542                nbdy_b = nbdy_a 
    543                nbdy_a = nbdy_a + 1 
    544                lect  =.true. 
     561               nb_bdy_fld(ib_bdy) = jfld - nb_bdy_fld_sum 
     562               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(ib_bdy) 
     563            ENDIF 
     564 
     565         ENDIF ! nn_dta .eq. 1 
     566      ENDDO ! ib_bdy 
     567 
     568 
     569      DO jfld = 1, nb_bdy_fld_sum 
     570         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) 
     571         IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) ) 
     572         nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld)) 
     573      ENDDO 
     574 
     575      ! fill bf with blf_i and control print 
     576      !------------------------------------- 
     577      jstart = 1 
     578      DO ib_bdy = 1, nb_bdy 
     579         jend = jstart + nb_bdy_fld(ib_bdy) - 1 
     580         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   & 
     581         &              'open boundary conditions', 'nambdy_dta' ) 
     582         jstart = jend + 1 
     583      ENDDO 
     584 
     585      ! Initialise local boundary data arrays 
     586      ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later 
     587      ! nn_xxx_dta=1 : point to "fnow" arrays 
     588      !------------------------------------- 
     589 
     590      jfld = 0 
     591      DO ib_bdy=1, nb_bdy 
     592 
     593         nblen => idx_bdy(ib_bdy)%nblen 
     594         nblenrim => idx_bdy(ib_bdy)%nblenrim 
     595 
     596         IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 
     597            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 
     598               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     599                  ilen0(1:3) = nblen(1:3) 
     600               ELSE 
     601                  ilen0(1:3) = nblenrim(1:3) 
     602               ENDIF 
     603               ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) ) 
     604               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
     605               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
    545606            ELSE 
    546                ! We have reached the end of the file 
    547                ! put the last data time into both time levels 
    548                nbdy_b = nbdy_a 
    549                IF(ln_tra_frs) THEN 
    550                   tbdydta(:,:,1) =  tbdydta(:,:,2) 
    551                   sbdydta(:,:,1) =  sbdydta(:,:,2) 
    552                ENDIF 
    553                IF(ln_dyn_frs) THEN 
    554                   ubdydta(:,:,1) =  ubdydta(:,:,2) 
    555                   vbdydta(:,:,1) =  vbdydta(:,:,2) 
    556                ENDIF 
    557 #if defined key_lim2 
    558                IF(ln_ice_frs) THEN 
    559                   frld_bdydta (:,1) =  frld_bdydta (:,2) 
    560                   hicif_bdydta(:,1) =  hicif_bdydta(:,2) 
    561                   hsnif_bdydta(:,1) =  hsnif_bdydta(:,2) 
    562                ENDIF 
    563 #endif 
    564             END IF ! nbdy_a < ntimes_bdy 
    565             ! 
    566         END IF 
    567           
    568         IF( lect ) THEN           ! Swap arrays 
    569            IF( ln_tra_frs ) THEN 
    570              tbdydta(:,:,1) =  tbdydta(:,:,2) 
    571              sbdydta(:,:,1) =  sbdydta(:,:,2) 
    572            ENDIF 
    573            IF( ln_dyn_frs ) THEN 
    574              ubdydta(:,:,1) =  ubdydta(:,:,2) 
    575              vbdydta(:,:,1) =  vbdydta(:,:,2) 
    576            ENDIF 
    577 #if defined key_lim2 
    578            IF( ln_ice_frs ) THEN 
    579              frld_bdydta (:,1) =  frld_bdydta (:,2) 
    580              hicif_bdydta(:,1) =  hicif_bdydta(:,2) 
    581              hsnif_bdydta(:,1) =  hsnif_bdydta(:,2) 
    582            ENDIF 
    583 #endif  
    584            ! read another set 
    585            ipj  = 1 
    586            ipk  = jpk 
    587  
    588            IF( ln_tra_frs ) THEN 
    589               !  
    590               igrd = 1                                   ! temperature 
    591               ipi  = nblendta(igrd) 
    592               CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
    593               DO ib = 1, nblen(igrd) 
    594                  DO ik = 1, jpkm1 
    595                     tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 
    596                  END DO 
    597               END DO 
    598               ! 
    599               igrd = 1                                   ! salinity 
    600               ipi  = nblendta(igrd) 
    601               CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
    602               DO ib = 1, nblen(igrd) 
    603                  DO ik = 1, jpkm1 
    604                     sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 
    605                  END DO 
    606               END DO 
    607            ENDIF ! ln_tra_frs 
    608  
    609            IF(ln_dyn_frs) THEN 
    610               ! 
    611               igrd = 2                                   ! u-velocity 
    612               ipi  = nblendta(igrd) 
    613               CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
    614               DO ib = 1, nblen(igrd) 
    615                 DO ik = 1, jpkm1 
    616                   ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
    617                 END DO 
    618               END DO 
    619               ! 
    620               igrd = 3                                   ! v-velocity 
    621               ipi  = nblendta(igrd) 
    622               CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
    623               DO ib = 1, nblen(igrd) 
    624                  DO ik = 1, jpkm1 
    625                     vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
    626                  END DO 
    627               END DO 
    628            ENDIF ! ln_dyn_frs 
    629            ! 
    630 #if defined key_lim2 
    631            IF(ln_ice_frs) THEN 
    632              ! 
    633              igrd = 1                                    ! ice concentration 
    634              ipi=nblendta(igrd) 
    635              CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a ) 
    636              DO ib=1, nblen(igrd) 
    637                frld_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 ) 
    638              END DO 
    639              ! 
    640              igrd=1                                      ! ice thickness 
    641              ipi=nblendta(igrd) 
    642              CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a ) 
    643              DO ib=1, nblen(igrd) 
    644                hicif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 ) 
    645              END DO 
    646              ! 
    647              igrd=1                                      ! snow thickness 
    648              ipi=nblendta(igrd) 
    649              CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a ) 
    650              DO ib=1, nblen(igrd) 
    651                hsnif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 ) 
    652              END DO 
    653            ENDIF ! ln_ice_frs 
    654 #endif 
    655            ! 
    656            IF(lwp) WRITE(numout,*) 'bdy_dta_frs : first record file used nbdy_b ',nbdy_b 
    657            IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a ',nbdy_a 
    658            IF (.NOT.ln_clim) THEN 
    659               IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep(nbdy_b) 
    660               IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer 
    661               IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a) 
    662            ENDIF 
    663            ! 
    664        ENDIF ! end lect=.true. 
    665  
    666  
    667        ! Interpolate linearly 
    668        ! ******************** 
    669        !  
    670        IF( ln_clim ) THEN   ;   zxy = REAL( nday                   ) / REAL( nmonth_len(nbdy_b) ) + 0.5 - i15 
    671        ELSEIF( istep(nbdy_b) == istep(nbdy_a) ) THEN  
    672                                     zxy = 0.0_wp 
    673        ELSE                     ;   zxy = REAL( istep(nbdy_b) - itimer ) / REAL( istep(nbdy_b) - istep(nbdy_a) ) 
    674        END IF 
    675  
    676           IF(ln_tra_frs) THEN 
    677              igrd = 1                                   ! temperature & salinity 
    678              DO ib = 1, nblen(igrd) 
    679                DO ik = 1, jpkm1 
    680                  tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1) 
    681                  sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1) 
    682                END DO 
    683              END DO 
    684           ENDIF 
    685  
    686           IF(ln_dyn_frs) THEN 
    687              igrd = 2                                   ! u-velocity 
    688              DO ib = 1, nblen(igrd) 
    689                DO ik = 1, jpkm1 
    690                  ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1)    
    691                END DO 
    692              END DO 
    693              ! 
    694              igrd = 3                                   ! v-velocity 
    695              DO ib = 1, nblen(igrd) 
    696                DO ik = 1, jpkm1 
    697                  vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1)    
    698                END DO 
    699              END DO 
    700           ENDIF 
    701  
    702 #if defined key_lim2 
    703           IF(ln_ice_frs) THEN 
    704             igrd=1 
    705             DO ib=1, nblen(igrd) 
    706                frld_bdy(ib) = zxy *  frld_bdydta(ib,2) + (1.-zxy) *  frld_bdydta(ib,1) 
    707               hicif_bdy(ib) = zxy * hicif_bdydta(ib,2) + (1.-zxy) * hicif_bdydta(ib,1) 
    708               hsnif_bdy(ib) = zxy * hsnif_bdydta(ib,2) + (1.-zxy) * hsnif_bdydta(ib,1) 
    709             END DO 
    710           ENDIF ! just if ln_ice_frs is true 
    711 #endif 
    712  
    713       END IF                       !end if ((nn_dtactl==1).AND.(ntimes_bdy>1)) 
    714      
    715  
    716       !                                                !---------------------! 
    717       !                                                !     last call       ! 
    718       !                                                !---------------------! 
    719       IF( kt == nitend ) THEN 
    720           IF(ln_tra_frs .or. ln_ice_frs) CALL iom_close( numbdyt )              ! Closing of the 3 files 
    721           IF(ln_dyn_frs) CALL iom_close( numbdyu ) 
    722           IF(ln_dyn_frs) CALL iom_close( numbdyv ) 
    723       ENDIF 
    724       ! 
    725       ENDIF ! ln_dyn_frs .OR. ln_tra_frs 
    726       ! 
    727    END SUBROUTINE bdy_dta_frs 
    728  
    729  
    730    SUBROUTINE bdy_dta_fla( kt, jit, icycl ) 
    731       !!--------------------------------------------------------------------------- 
    732       !!                      ***  SUBROUTINE bdy_dta_fla  *** 
    733       !!                     
    734       !! ** Purpose :   Read unstructured boundary data for Flather condition 
    735       !! 
    736       !! ** Method  :  At the first timestep, read in boundary data for two 
    737       !!               times from the file and time-interpolate. At other  
    738       !!               timesteps, check to see if we need another time from  
    739       !!               the file. If so read it in. Time interpolate. 
    740       !!--------------------------------------------------------------------------- 
    741 !!gm DOCTOR names :   argument integer :  start with "k" 
    742       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    743       INTEGER, INTENT( in ) ::   jit         ! barotropic time step index 
    744       INTEGER, INTENT( in ) ::   icycl       ! number of cycles need for final file close 
    745       !                                      ! (for timesplitting option, otherwise zero) 
    746       !! 
    747       LOGICAL ::   lect                      ! flag for reading 
    748       INTEGER ::   it, ib, igrd              ! dummy loop indices 
    749       INTEGER ::   idvar                     ! netcdf var ID 
    750       INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing 
    751       INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 
    752       INTEGER ::   itimer, totime 
    753       INTEGER ::   ipi, ipj, ipk, inum       ! temporary integers (NetCDF read) 
    754       INTEGER ::   iyear0, imonth0, iday0 
    755       INTEGER ::   ihours0, iminutes0, isec0 
    756       INTEGER ::   iyear, imonth, iday, isecs 
    757       INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files 
    758       REAL(wp) ::   dayfrac, zxy, zoffsett 
    759       REAL(wp) ::   zoffsetu, zoffsetv 
    760       REAL(wp) ::   dayjul0, zdayjulini 
    761       REAL(wp) ::   zinterval_s, zinterval_e                    ! First and last interval in time axis 
    762       REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files 
    763       REAL(wp), DIMENSION(jpbdta,1)     ::   zdta               ! temporary array for data fields 
    764       CHARACTER(LEN=80), DIMENSION(6)   ::   clfile 
    765       CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate 
    766       !!--------------------------------------------------------------------------- 
    767  
    768 !!gm   add here the same style as in bdy_dta_frs 
    769 !!gm      clearly bdy_dta_fla and bdy_dta_frs  can be combined...    
    770 !!gm      too many things duplicated in the read of data...   simplification can be done 
    771  
    772       ! -------------------- ! 
    773       !    Initialization    ! 
    774       ! -------------------- ! 
    775  
    776       lect   = .false.           ! If true, read a time record 
    777  
    778       ! Some time variables for monthly climatological forcing: 
    779       ! ******************************************************* 
    780  !!gm  here  use directely daymod variables 
    781   
    782       iman  = INT( raamo ) ! Number of months in a year 
    783  
    784       i15 = INT( 2*REAL( nday, wp ) / ( REAL( nmonth_len(nmonth), wp ) + 0.5 ) ) 
    785       ! i15=0 if the current day is in the first half of the month, else i15=1 
    786  
    787       imois = nmonth + i15 - 1            ! imois is the first month record 
    788       IF( imois == 0 ) imois = iman 
    789  
    790       ! Time variable for non-climatological forcing: 
    791       ! ********************************************* 
    792  
    793       itimer = ((kt-1)-nit000+1)*rdt                      ! current time in seconds for interpolation  
    794       itimer = itimer + jit*rdt/REAL(nn_baro,wp)      ! in non-climatological case 
    795  
    796       IF ( ln_tides ) THEN 
    797  
    798          ! -------------------------------------! 
    799          ! Update BDY fields with tidal forcing ! 
    800          ! -------------------------------------!   
    801  
    802          CALL tide_update( kt, jit )  
    803    
    804       ENDIF 
    805  
    806       IF ( ln_dyn_fla ) THEN 
    807  
    808          ! -------------------------------------! 
    809          ! Update BDY fields with model data    ! 
    810          ! -------------------------------------!   
    811  
    812       !                                                !-------------------! 
    813       IF( kt == nit000 .and. jit ==2 ) THEN            !  First call only  ! 
    814          !                                             !-------------------! 
    815          istep_bt(:) = 0 
    816          nbdy_b_bt    = 0 
    817          nbdy_a_bt    = 0 
    818  
    819          ! Get time information from bdy data file 
    820          ! *************************************** 
    821  
    822         IF(lwp) WRITE(numout,*) 
    823         IF(lwp) WRITE(numout,*)    'bdy_dta_fla :Initialize unstructured boundary data for barotropic variables.' 
    824         IF(lwp) WRITE(numout,*)    '~~~~~~~'  
    825  
    826         IF( nn_dtactl == 0 ) THEN 
    827           IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions' 
    828  
    829         ELSEIF (nn_dtactl == 1) THEN 
    830           IF(lwp) WRITE(numout,*)  'Bdy data are read in netcdf files' 
    831  
    832           dayfrac = adatrj  - REAL(itimer,wp)/86400. ! day fraction at time step kt-1 
    833           dayfrac = dayfrac - INT (dayfrac)          ! 
    834           totime = (nitend-nit000+1)*rdt             ! Total time of the run to verify that all the 
    835                                                      ! necessary time dumps in file are included 
    836  
    837           clfile(4) = cn_dta_fla_T 
    838           clfile(5) = cn_dta_fla_U 
    839           clfile(6) = cn_dta_fla_V 
    840  
    841           DO igrd = 4,6 
    842  
    843             CALL iom_open( clfile(igrd), inum ) 
    844             CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy_bt, cdunits=clunits )  
    845  
    846             SELECT CASE( igrd ) 
    847                CASE (4)  
    848                   numbdyt_bt = inum 
    849                CASE (5)  
    850                   numbdyu_bt = inum 
    851                CASE (6)  
    852                   numbdyv_bt = inum 
    853             END SELECT 
    854  
    855             ! Calculate time offset  
    856             READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 
    857             ! Convert time origin in file to julian days  
    858             isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 
    859             CALL ymds2ju(iyear0, imonth0, iday0, REAL(isec0, wp), dayjul0) 
    860             ! Compute model initialization time  
    861             iyear  = ndastp / 10000 
    862             imonth = ( ndastp - iyear * 10000 ) / 100 
    863             iday   = ndastp - iyear * 10000 - imonth * 100 
    864             isecs  = dayfrac * 86400 
    865             CALL ymds2ju(iyear, imonth, iday, REAL(isecs, wp) , zdayjulini) 
    866             ! zoffset from initialization date: 
    867             zoffset = (dayjul0-zdayjulini)*86400 
    868             ! 
    869  
    870 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 
    871  
    872             !! TO BE DONE... Check consistency between calendar from file  
    873             !! (available optionally from iom_gettime) and calendar in model  
    874             !! when calendar in model available outside of IOIPSL. 
    875  
    876             ! Check that there are not too many times in the file.  
    877             IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( & 
    878                  'Number of time dumps in bdy file exceed jpbtime parameter', & 
    879                  'Check file:' // TRIM(clfile(igrd))  ) 
    880  
    881             ! Check that time array increases (or interp will fail): 
    882             DO it = 2, ntimes_bdy_bt 
    883                IF ( zstepr(it-1) >= zstepr(it) ) THEN 
    884                   CALL ctl_stop('Time array in unstructured boundary data file', & 
    885                        'does not continuously increase.',               & 
    886                        'Check file:' // TRIM(clfile(igrd))  ) 
    887                   EXIT 
    888                END IF 
    889             END DO 
    890  
    891             IF ( .NOT. ln_clim ) THEN 
    892                ! Check that times in file span model run time: 
    893  
    894                ! Note: the fields may be time means, so we allow nit000 to be before 
    895                ! first time in the file, provided that it falls inside the meaning 
    896                ! period of the first field.  Until we can get the meaning period 
    897                ! from the file, use the interval between fields as a proxy. 
    898                ! If nit000 is before the first time, use the value at first time 
    899                ! instead of extrapolating.  This is done by putting time 1 into 
    900                ! both time levels. 
    901                ! The same applies to the last time level: see setting of lect below. 
    902  
    903                IF ( ntimes_bdy_bt == 1 ) CALL ctl_stop( & 
    904                     'There is only one time dump in data files', & 
    905                     'Set ln_clim=.true. in namelist for constant bdy forcing.' ) 
    906  
    907                zinterval_s = zstepr(2) - zstepr(1) 
    908                zinterval_e = zstepr(ntimes_bdy_bt) - zstepr(ntimes_bdy_bt-1) 
    909  
    910                IF( zstepr(1) + zoffset > 0 ) THEN 
    911                      WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
    912                      CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 ) 
    913                END IF 
    914                IF( zstepr(ntimes_bdy_bt) + zoffset < totime ) THEN 
    915                      WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
    916                      CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 ) 
    917                END IF 
    918             END IF ! .NOT. ln_clim 
    919  
    920             IF ( igrd .EQ. 4) THEN 
    921               ntimes_bdyt = ntimes_bdy_bt 
    922               zoffsett = zoffset 
    923               istept(:) = INT( zstepr(:) + zoffset ) 
    924             ELSE IF (igrd .EQ. 5) THEN 
    925               ntimes_bdyu = ntimes_bdy_bt 
    926               zoffsetu = zoffset 
    927               istepu(:) = INT( zstepr(:) + zoffset ) 
    928             ELSE IF (igrd .EQ. 6) THEN 
    929               ntimes_bdyv = ntimes_bdy_bt 
    930               zoffsetv = zoffset 
    931               istepv(:) = INT( zstepr(:) + zoffset ) 
    932             ENDIF 
    933  
    934           ENDDO 
    935  
    936       ! Verify time consistency between files   
    937  
    938           IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 
    939              CALL ctl_stop( & 
    940              'Time axis lengths differ between bdy data files', & 
    941              'Multiple time frequencies not implemented yet' ) 
    942           ELSE 
    943             ntimes_bdy_bt = ntimes_bdyt 
    944           ENDIF 
    945  
    946           IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN 
    947             CALL ctl_stop( &  
    948             'Bdy data files must have the same time origin', & 
    949             'Multiple time frequencies not implemented yet'  ) 
    950           ENDIF 
    951           zoffset = zoffsett 
    952  
    953       !! Check that times are the same in the three files... HERE. 
    954           istep_bt(:) = istept(:) 
    955  
    956       ! Check number of time dumps:               
    957           IF (ln_clim) THEN 
    958             SELECT CASE ( ntimes_bdy_bt ) 
    959             CASE( 1 ) 
    960               IF(lwp) WRITE(numout,*) 
    961               IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 
    962               IF(lwp) WRITE(numout,*)              
    963             CASE( 12 ) 
    964               IF(lwp) WRITE(numout,*) 
    965               IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 
    966               IF(lwp) WRITE(numout,*)  
    967             CASE DEFAULT 
    968               CALL ctl_stop( & 
    969                 'For climatological boundary forcing (ln_clim=.true.),',& 
    970                 'bdy data files must contain 1 or 12 time dumps.' ) 
    971             END SELECT 
    972           ENDIF 
    973  
    974       ! Find index of first record to read (before first model time).  
    975  
    976           it=1 
    977           DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1))) 
    978             it=it+1 
    979           END DO 
    980           nbdy_b_bt = it 
    981  
    982           IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset 
    983           IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b_bt 
    984  
    985         ENDIF ! endif (nn_dtactl == 1) 
    986  
    987       ! 1.2  Read first record in file if necessary (ie if nn_dtactl == 1) 
    988       ! ***************************************************************** 
    989  
    990         IF ( nn_dtactl == 0) THEN 
    991           ! boundary data arrays are filled with initial conditions 
    992           igrd = 5            ! U-points data  
    993           DO ib = 1, nblen(igrd)               
    994             ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1) 
    995           END DO 
    996  
    997           igrd = 6            ! V-points data  
    998           DO ib = 1, nblen(igrd)               
    999             vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1) 
    1000           END DO 
    1001  
    1002           igrd = 4            ! T-points data  
    1003           DO ib = 1, nblen(igrd)               
    1004             sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd)) 
    1005           END DO 
    1006  
    1007         ELSEIF (nn_dtactl == 1) THEN 
    1008   
    1009         ! Set first record in the climatological case:    
    1010           IF ((ln_clim).AND.(ntimes_bdy_bt==1)) THEN 
    1011             nbdy_a_bt = 1 
    1012           ELSEIF ((ln_clim).AND.(ntimes_bdy_bt==iman)) THEN 
    1013             nbdy_b_bt = 0 
    1014             nbdy_a_bt = imois 
    1015           ELSE 
    1016             nbdy_a_bt = nbdy_b_bt 
    1017           END IF 
    1018   
    1019          ! Open Netcdf files: 
    1020  
    1021           CALL iom_open ( cn_dta_fla_T, numbdyt_bt ) 
    1022           CALL iom_open ( cn_dta_fla_U, numbdyu_bt ) 
    1023           CALL iom_open ( cn_dta_fla_V, numbdyv_bt ) 
    1024  
    1025          ! Read first record: 
    1026           ipj=1 
    1027           igrd=4 
    1028           ipi=nblendta(igrd) 
    1029  
    1030           ! ssh 
    1031           igrd=4 
    1032           IF ( nblendta(igrd) .le. 0 ) THEN  
    1033             idvar = iom_varid( numbdyt_bt,'sossheig' ) 
    1034             nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar) 
    1035           ENDIF 
    1036           WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd) 
    1037           ipi=nblendta(igrd) 
    1038  
    1039           CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
    1040  
    1041           DO ib=1, nblen(igrd) 
    1042             sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    1043           END DO 
    1044   
    1045           ! u-velocity 
    1046           igrd=5 
    1047           IF ( nblendta(igrd) .le. 0 ) THEN  
    1048             idvar = iom_varid( numbdyu_bt,'vobtcrtx' ) 
    1049             nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar) 
    1050           ENDIF 
    1051           WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd) 
    1052           ipi=nblendta(igrd) 
    1053  
    1054           CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
    1055  
    1056           DO ib=1, nblen(igrd) 
    1057             ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    1058           END DO 
    1059  
    1060           ! v-velocity 
    1061           igrd=6 
    1062           IF ( nblendta(igrd) .le. 0 ) THEN  
    1063             idvar = iom_varid( numbdyv_bt,'vobtcrty' ) 
    1064             nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar) 
    1065           ENDIF 
    1066           WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd) 
    1067           ipi=nblendta(igrd) 
    1068  
    1069           CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
    1070  
    1071           DO ib=1, nblen(igrd) 
    1072             vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    1073           END DO 
    1074  
    1075         END IF 
    1076   
    1077         ! In the case of constant boundary forcing fill bdy arrays once for all 
    1078         IF ((ln_clim).AND.(ntimes_bdy_bt==1)) THEN 
    1079  
    1080           ubtbdy  (:) = ubtbdydta  (:,2) 
    1081           vbtbdy  (:) = vbtbdydta  (:,2) 
    1082           sshbdy  (:) = sshbdydta  (:,2) 
    1083  
    1084           CALL iom_close( numbdyt_bt ) 
    1085           CALL iom_close( numbdyu_bt ) 
    1086           CALL iom_close( numbdyv_bt ) 
    1087  
    1088         END IF 
    1089  
    1090       ENDIF ! End if nit000 
    1091  
    1092       ! -------------------- ! 
    1093       ! 2. At each time step ! 
    1094       ! -------------------- ! 
    1095  
    1096       IF ((nn_dtactl==1).AND.(ntimes_bdy_bt>1)) THEN  
    1097  
    1098       ! 2.1 Read one more record if necessary 
    1099       !************************************** 
    1100  
    1101         IF ( (ln_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000 
    1102          nbdy_b_bt = imois 
    1103          nbdy_a_bt = imois+1 
    1104          nbdy_b_bt = MOD( nbdy_b_bt, iman ) 
    1105          IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman 
    1106          nbdy_a_bt = MOD( nbdy_a_bt, iman ) 
    1107          IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman 
    1108          lect=.true. 
    1109  
    1110         ELSEIF ((.NOT.ln_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN 
    1111           nbdy_b_bt=nbdy_a_bt 
    1112           nbdy_a_bt=nbdy_a_bt+1 
    1113           lect=.true. 
    1114         END IF 
    1115           
    1116         IF (lect) THEN 
    1117  
    1118         ! Swap arrays 
    1119           sshbdydta(:,1) =  sshbdydta(:,2) 
    1120           ubtbdydta(:,1) =  ubtbdydta(:,2) 
    1121           vbtbdydta(:,1) =  vbtbdydta(:,2) 
    1122   
    1123         ! read another set 
    1124  
    1125           ipj=1 
    1126           ipk=jpk 
    1127           igrd=4 
    1128           ipi=nblendta(igrd) 
    1129  
    1130            
    1131           ! ssh 
    1132           igrd=4 
    1133           ipi=nblendta(igrd) 
    1134  
    1135           CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
    1136  
    1137           DO ib=1, nblen(igrd) 
    1138             sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    1139           END DO 
    1140  
    1141           ! u-velocity 
    1142           igrd=5 
    1143           ipi=nblendta(igrd) 
    1144  
    1145           CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
    1146  
    1147           DO ib=1, nblen(igrd) 
    1148             ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    1149           END DO 
    1150  
    1151           ! v-velocity 
    1152           igrd=6 
    1153           ipi=nblendta(igrd) 
    1154  
    1155           CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
    1156  
    1157           DO ib=1, nblen(igrd) 
    1158             vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    1159           END DO 
    1160  
    1161  
    1162          IF(lwp) WRITE(numout,*) 'bdy_dta_fla : first record file used nbdy_b_bt ',nbdy_b_bt 
    1163          IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a_bt ',nbdy_a_bt 
    1164          IF (.NOT.ln_clim) THEN 
    1165            IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy_b_bt) 
    1166            IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer 
    1167            IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt) 
    1168          ENDIF 
    1169         END IF ! end lect=.true. 
    1170  
    1171  
    1172       ! 2.2   Interpolate linearly: 
    1173       ! *************************** 
    1174      
    1175         IF (ln_clim) THEN 
    1176           zxy = REAL( nday, wp ) / REAL( nmonth_len(nbdy_b_bt), wp ) + 0.5 - i15 
    1177         ELSE           
    1178           zxy = REAL(istep_bt(nbdy_b_bt)-itimer, wp) / REAL(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt), wp) 
    1179         END IF 
    1180  
    1181           igrd=4 
    1182           DO ib=1, nblen(igrd) 
    1183             sshbdy(ib) = zxy      * sshbdydta(ib,2) + & 
    1184                        (1.-zxy) * sshbdydta(ib,1)    
    1185           END DO 
    1186  
    1187           igrd=5 
    1188           DO ib=1, nblen(igrd) 
    1189             ubtbdy(ib) = zxy      * ubtbdydta(ib,2) + & 
    1190                          (1.-zxy) * ubtbdydta(ib,1)    
    1191           END DO 
    1192  
    1193           igrd=6 
    1194           DO ib=1, nblen(igrd) 
    1195             vbtbdy(ib) = zxy      * vbtbdydta(ib,2) + & 
    1196                          (1.-zxy) * vbtbdydta(ib,1)    
    1197           END DO 
    1198  
    1199  
    1200       END IF !end if ((nn_dtactl==1).AND.(ntimes_bdy_bt>1)) 
    1201      
    1202       ! ------------------- ! 
    1203       ! Last call kt=nitend ! 
    1204       ! ------------------- ! 
    1205  
    1206       ! Closing of the 3 files 
    1207       IF( kt == nitend   .and. jit == icycl ) THEN 
    1208           CALL iom_close( numbdyt_bt ) 
    1209           CALL iom_close( numbdyu_bt ) 
    1210           CALL iom_close( numbdyv_bt ) 
    1211       ENDIF 
    1212  
    1213       ENDIF ! ln_dyn_frs 
    1214  
    1215       END SUBROUTINE bdy_dta_fla 
    1216  
     607               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
     608                  jfld = jfld + 1 
     609                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
     610               ENDIF 
     611               jfld = jfld + 1 
     612               dta_bdy(ib_bdy)%u2d => bf(jfld)%fnow(:,1,1) 
     613               jfld = jfld + 1 
     614               dta_bdy(ib_bdy)%v2d => bf(jfld)%fnow(:,1,1) 
     615            ENDIF 
     616         ENDIF 
     617 
     618         IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
     619            IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     620               ilen0(1:3) = nblen(1:3) 
     621            ELSE 
     622               ilen0(1:3) = nblenrim(1:3) 
     623            ENDIF 
     624            ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 
     625            ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) 
     626         ENDIF 
     627         IF ( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ).or. & 
     628           &  ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.   & 
     629           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
     630            jfld = jfld + 1 
     631            dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 
     632            jfld = jfld + 1 
     633            dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 
     634         ENDIF 
     635 
     636         IF (nn_tra(ib_bdy) .gt. 0) THEN 
     637            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
     638               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     639                  ilen0(1:3) = nblen(1:3) 
     640               ELSE 
     641                  ilen0(1:3) = nblenrim(1:3) 
     642               ENDIF 
     643               ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 
     644               ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) 
     645            ELSE 
     646               jfld = jfld + 1 
     647               dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:) 
     648               jfld = jfld + 1 
     649               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:) 
     650            ENDIF 
     651         ENDIF 
     652 
     653#if defined key_lim2 
     654         IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
     655            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
     656               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     657                  ilen0(1:3) = nblen(1:3) 
     658               ELSE 
     659                  ilen0(1:3) = nblenrim(1:3) 
     660               ENDIF 
     661               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 
     662               ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) 
     663               ALLOCATE( dta_bdy(ib_bdy)%hsnif(ilen0(1)) ) 
     664            ELSE 
     665               jfld = jfld + 1 
     666               dta_bdy(ib_bdy)%frld  => bf(jfld)%fnow(:,1,1) 
     667               jfld = jfld + 1 
     668               dta_bdy(ib_bdy)%hicif => bf(jfld)%fnow(:,1,1) 
     669               jfld = jfld + 1 
     670               dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1) 
     671            ENDIF 
     672         ENDIF 
     673#endif 
     674 
     675      ENDDO ! ib_bdy  
     676 
     677      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_init') 
     678 
     679      END SUBROUTINE bdy_dta_init 
    1217680 
    1218681#else 
    1219682   !!---------------------------------------------------------------------- 
    1220    !!   Dummy module                   NO Unstruct Open Boundary Conditions 
     683   !!   Dummy module                   NO Open Boundary Conditions 
    1221684   !!---------------------------------------------------------------------- 
    1222685CONTAINS 
    1223    SUBROUTINE bdy_dta_frs( kt )              ! Empty routine 
    1224       WRITE(*,*) 'bdy_dta_frs: You should not have seen this print! error?', kt 
    1225    END SUBROUTINE bdy_dta_frs 
    1226    SUBROUTINE bdy_dta_fla( kt, kit, icycle )      ! Empty routine 
    1227       WRITE(*,*) 'bdy_dta_frs: You should not have seen this print! error?', kt, kit 
    1228    END SUBROUTINE bdy_dta_fla 
     686   SUBROUTINE bdy_dta( kt, jit, time_offset ) ! Empty routine 
     687      INTEGER, INTENT( in )           ::   kt     
     688      INTEGER, INTENT( in ), OPTIONAL ::   jit    
     689      INTEGER, INTENT( in ), OPTIONAL ::   time_offset 
     690      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt 
     691   END SUBROUTINE bdy_dta 
     692   SUBROUTINE bdy_dta_init()                  ! Empty routine 
     693      WRITE(*,*) 'bdy_dta_init: You should not have seen this print! error?' 
     694   END SUBROUTINE bdy_dta_init 
    1229695#endif 
    1230696 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r2528 r3294  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  bdydyn  *** 
    4    !! Unstructured Open Boundary Cond. :   Flow relaxation scheme on velocities 
     4   !! Unstructured Open Boundary Cond. :   Apply boundary conditions to velocities 
    55   !!====================================================================== 
    66   !! History :  1.0  !  2005-02  (J. Chanut, A. Sellar)  Original code 
    … …  
    1010   !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
     12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_bdy  
    … …  
    1516   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    1617   !!---------------------------------------------------------------------- 
    17    !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary 
    18    !!   bdy_dyn_fla    : Flather condition for barotropic solution 
     18   !!   bdy_dyn        : split velocities into barotropic and baroclinic parts 
     19   !!                    and call bdy_dyn2d and bdy_dyn3d to apply boundary 
     20   !!                    conditions 
    1921   !!---------------------------------------------------------------------- 
     22   USE wrk_nemo        ! Memory Allocation 
     23   USE timing          ! Timing 
    2024   USE oce             ! ocean dynamics and tracers  
    2125   USE dom_oce         ! ocean space and time domain 
     26   USE dynspg_oce       
    2227   USE bdy_oce         ! ocean open boundary conditions 
    23    USE dynspg_oce      ! for barotropic variables 
    24    USE phycst          ! physical constants 
     28   USE bdydyn2d        ! open boundary conditions for barotropic solution 
     29   USE bdydyn3d        ! open boundary conditions for baroclinic velocities 
    2530   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    26    USE bdytides        ! for tidal harmonic forcing at boundary 
    2731   USE in_out_manager  ! 
    2832 
    … …  
    3034   PRIVATE 
    3135 
    32    PUBLIC   bdy_dyn_frs   ! routine called in dynspg_flt (free surface case ONLY) 
    33 # if defined key_dynspg_exp || defined key_dynspg_ts 
    34    PUBLIC   bdy_dyn_fla   ! routine called in dynspg_flt (free surface case ONLY) 
    35 # endif 
     36   PUBLIC   bdy_dyn     ! routine called in dynspg_flt (if lk_dynspg_flt) or  
     37                        ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp) 
    3638 
     39#  include "domzgr_substitute.h90" 
    3740   !!---------------------------------------------------------------------- 
    3841   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    4245CONTAINS 
    4346 
    44    SUBROUTINE bdy_dyn_frs( kt ) 
     47   SUBROUTINE bdy_dyn( kt, dyn3d_only ) 
    4548      !!---------------------------------------------------------------------- 
    46       !!                  ***  SUBROUTINE bdy_dyn_frs  *** 
     49      !!                  ***  SUBROUTINE bdy_dyn  *** 
    4750      !! 
    48       !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the   
    49       !!                case of unstructured open boundaries. 
     51      !! ** Purpose : - Wrapper routine for bdy_dyn2d and bdy_dyn3d. 
    5052      !! 
    51       !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in  
    52       !!               a three-dimensional baroclinic ocean model with realistic 
    53       !!               topography. Tellus, 365-382. 
    5453      !!---------------------------------------------------------------------- 
    55       INTEGER, INTENT( in ) ::   kt   ! Main time step counter 
    5654      !! 
    57       INTEGER  ::   jb, jk         ! dummy loop indices 
    58       INTEGER  ::   ii, ij, igrd   ! local integers 
    59       REAL(wp) ::   zwgt           ! boundary weight 
    60       !!---------------------------------------------------------------------- 
    61       ! 
    62       IF(ln_dyn_frs) THEN       ! If this is false, then this routine does nothing.  
    63          ! 
    64          IF( kt == nit000 ) THEN 
    65             IF(lwp) WRITE(numout,*) 
    66             IF(lwp) WRITE(numout,*) 'bdy_dyn_frs : Flow Relaxation Scheme on momentum' 
    67             IF(lwp) WRITE(numout,*) '~~~~~~~' 
    68          ENDIF 
    69          ! 
    70          igrd = 2                      ! Relaxation of zonal velocity 
    71          DO jb = 1, nblen(igrd) 
    72             DO jk = 1, jpkm1 
    73                ii   = nbi(jb,igrd) 
    74                ij   = nbj(jb,igrd) 
    75                zwgt = nbw(jb,igrd) 
    76                ua(ii,ij,jk) = ( ua(ii,ij,jk) * ( 1.- zwgt ) + ubdy(jb,jk) * zwgt ) * umask(ii,ij,jk) 
    77             END DO 
    78          END DO 
    79          ! 
    80          igrd = 3                      ! Relaxation of meridional velocity 
    81          DO jb = 1, nblen(igrd) 
    82             DO jk = 1, jpkm1 
    83                ii   = nbi(jb,igrd) 
    84                ij   = nbj(jb,igrd) 
    85                zwgt = nbw(jb,igrd) 
    86                va(ii,ij,jk) = ( va(ii,ij,jk) * ( 1.- zwgt ) + vbdy(jb,jk) * zwgt ) * vmask(ii,ij,jk) 
    87             END DO 
    88          END DO  
    89          CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
    90          ! 
    91       ENDIF ! ln_dyn_frs 
    92       ! 
    93    END SUBROUTINE bdy_dyn_frs 
     55      INTEGER, INTENT( in )           :: kt               ! Main time step counter 
     56      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities 
     57      !! 
     58      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter 
     59      LOGICAL               :: ll_dyn2d, ll_dyn3d   
     60      !! 
    9461 
     62      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn') 
    9563 
    96 # if defined   key_dynspg_exp   ||   defined key_dynspg_ts 
    97    !!---------------------------------------------------------------------- 
    98    !!   'key_dynspg_exp'        OR              explicit sea surface height 
    99    !!   'key_dynspg_ts '                  split-explicit sea surface height 
    100    !!---------------------------------------------------------------------- 
    101     
    102 !! Option to use Flather with dynspg_flt not coded yet... 
     64      ll_dyn2d = .true. 
     65      ll_dyn3d = .true. 
    10366 
    104    SUBROUTINE bdy_dyn_fla( pssh ) 
    105       !!---------------------------------------------------------------------- 
    106       !!                 ***  SUBROUTINE bdy_dyn_fla  *** 
    107       !!              
    108       !!              - Apply Flather boundary conditions on normal barotropic velocities  
    109       !!                (ln_dyn_fla=.true. or ln_tides=.true.) 
    110       !! 
    111       !! ** WARNINGS about FLATHER implementation: 
    112       !!1. According to Palma and Matano, 1998 "after ssh" is used.  
    113       !!   In ROMS and POM implementations, it is "now ssh". In the current  
    114       !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.  
    115       !!   So I use "before ssh" in the following. 
    116       !! 
    117       !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of  
    118       !!   fact, the model ssh just inside the dynamical boundary is used (the outside   
    119       !!   ssh in the code is not updated). 
    120       !! 
    121       !! References:  Flather, R. A., 1976: A tidal model of the northwest European 
    122       !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.      
    123       !!---------------------------------------------------------------------- 
    124       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh 
     67      IF( PRESENT(dyn3d_only) ) THEN 
     68         IF( dyn3d_only ) ll_dyn2d = .false. 
     69      ENDIF 
    12570 
    126       INTEGER  ::   jb, igrd                         ! dummy loop indices 
    127       INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
    128       REAL(wp) ::   zcorr                            ! Flather correction 
    129       REAL(wp) ::   zforc                            ! temporary scalar 
    130       !!---------------------------------------------------------------------- 
     71      !------------------------------------------------------- 
     72      ! Set pointers 
     73      !------------------------------------------------------- 
    13174 
    132       ! ---------------------------------! 
    133       ! Flather boundary conditions     :! 
    134       ! ---------------------------------!  
    135       
    136       IF(ln_dyn_fla .OR. ln_tides) THEN ! If these are both false, then this routine does nothing.  
     75      pssh => sshn 
     76      phur => hur 
     77      phvr => hvr 
     78      CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
    13779 
    138          ! Fill temporary array with ssh data (here spgu): 
    139          igrd = 4 
    140          spgu(:,:) = 0.0 
    141          DO jb = 1, nblenrim(igrd) 
    142             ii = nbi(jb,igrd) 
    143             ij = nbj(jb,igrd) 
    144             IF( ln_dyn_fla ) spgu(ii, ij) = sshbdy(jb) 
    145             IF( ln_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb) 
    146          END DO 
    147          ! 
    148          igrd = 5      ! Flather bc on u-velocity;  
    149          !             ! remember that flagu=-1 if normal velocity direction is outward 
    150          !             ! I think we should rather use after ssh ? 
    151          DO jb = 1, nblenrim(igrd) 
    152             ii  = nbi(jb,igrd) 
    153             ij  = nbj(jb,igrd)  
    154             iim1 = ii + MAX( 0, INT( flagu(jb) ) )   ! T pts i-indice inside the boundary 
    155             iip1 = ii - MIN( 0, INT( flagu(jb) ) )   ! T pts i-indice outside the boundary  
    156             ! 
    157             zcorr = - flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    158             zforc = ubtbdy(jb) + utide(jb) 
    159             ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
    160          END DO 
    161          ! 
    162          igrd = 6      ! Flather bc on v-velocity 
    163          !             ! remember that flagv=-1 if normal velocity direction is outward 
    164          DO jb = 1, nblenrim(igrd) 
    165             ii  = nbi(jb,igrd) 
    166             ij  = nbj(jb,igrd)  
    167             ijm1 = ij + MAX( 0, INT( flagv(jb) ) )   ! T pts j-indice inside the boundary 
    168             ijp1 = ij - MIN( 0, INT( flagv(jb) ) )   ! T pts j-indice outside the boundary  
    169             ! 
    170             zcorr = - flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    171             zforc = vbtbdy(jb) + vtide(jb) 
    172             va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
    173          END DO 
    174          CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated 
    175          CALL lbc_lnk( va_e, 'V', -1. )   ! 
    176          ! 
    177       ENDIF ! ln_dyn_fla .or. ln_tides 
    178       ! 
    179    END SUBROUTINE bdy_dyn_fla 
    180 #endif 
     80      !------------------------------------------------------- 
     81      ! Split velocities into barotropic and baroclinic parts 
     82      !------------------------------------------------------- 
     83 
     84      pu2d(:,:) = 0.e0 
     85      pv2d(:,:) = 0.e0 
     86      DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     87          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     88          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     89      END DO 
     90      pu2d(:,:) = pu2d(:,:) * phur(:,:) 
     91      pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     92      DO jk = 1 , jpkm1 
     93         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) 
     94         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) 
     95      END DO 
     96 
     97      !------------------------------------------------------- 
     98      ! Apply boundary conditions to barotropic and baroclinic 
     99      ! parts separately 
     100      !------------------------------------------------------- 
     101 
     102      IF( ll_dyn2d ) CALL bdy_dyn2d( kt ) 
     103 
     104      IF( ll_dyn3d ) CALL bdy_dyn3d( kt ) 
     105 
     106      !------------------------------------------------------- 
     107      ! Recombine velocities 
     108      !------------------------------------------------------- 
     109 
     110      DO jk = 1 , jpkm1 
     111         ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk) 
     112         va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk) 
     113      END DO 
     114 
     115      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     116 
     117      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn') 
     118 
     119   END SUBROUTINE bdy_dyn 
    181120 
    182121#else 
    … …  
    185124   !!---------------------------------------------------------------------- 
    186125CONTAINS 
    187    SUBROUTINE bdy_dyn_frs( kt )      ! Empty routine 
    188       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
    189    END SUBROUTINE bdy_dyn_frs 
    190    SUBROUTINE bdy_dyn_fla( pssh )    ! Empty routine 
    191       REAL :: pssh(:,:) 
    192       WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?', pssh(1,1) 
    193    END SUBROUTINE bdy_dyn_fla 
     126   SUBROUTINE bdy_dyn( kt )      ! Empty routine 
     127      WRITE(*,*) 'bdy_dyn: You should not have seen this print! error?', kt 
     128   END SUBROUTINE bdy_dyn 
    194129#endif 
    195130 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r2715 r3294  
    1010   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations 
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
     12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_bdy 
    … …  
    1718   !!   bdy_init       : Initialization of unstructured open boundaries 
    1819   !!---------------------------------------------------------------------- 
     20   USE timing          ! Timing 
    1921   USE oce             ! ocean dynamics and tracers variables 
    2022   USE dom_oce         ! ocean space and time domain 
    21    USE obc_par         ! ocean open boundary conditions 
    2223   USE bdy_oce         ! unstructured open boundary conditions 
    23    USE bdydta, ONLY: bdy_dta_alloc ! open boundary data 
    24    USE bdytides        ! tides at open boundaries initialization (tide_init routine) 
    2524   USE in_out_manager  ! I/O units 
    2625   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    3130   PRIVATE 
    3231 
    33    PUBLIC   bdy_init   ! routine called by opa.F90 
     32   PUBLIC   bdy_init   ! routine called in nemo_init 
    3433 
    3534   !!---------------------------------------------------------------------- 
    … …  
    5251      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries 
    5352      !!----------------------------------------------------------------------       
    54       INTEGER  ::   ii, ij, ik, igrd, ib, ir   ! dummy loop indices 
    55       INTEGER  ::   icount, icountr, ib_len, ibr_max   ! local integers 
    56       INTEGER  ::   iw, ie, is, in, inum, id_dummy     !   -       - 
    57       INTEGER  ::   igrd_start, igrd_end               !   -       - 
    58       REAL(wp) ::   zefl, zwfl, znfl, zsfl              ! local scalars 
    59       INTEGER, DIMENSION (2)             ::   kdimsz 
    60       INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    61       INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbrdta           ! Discrete distance from rim points 
    62       REAL(wp), DIMENSION(jpidta,jpjdta) ::   zmask            ! global domain mask 
    63       REAL(wp), DIMENSION(jpbdta,1)      ::   zdta             ! temporary array  
    64       CHARACTER(LEN=80),DIMENSION(6)     ::   clfile 
     53      ! namelist variables 
     54      !------------------- 
     55      INTEGER, PARAMETER          :: jp_nseg = 100 
     56      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
     57      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft 
     58      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft 
     59      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft 
     60      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft 
     61 
     62      ! local variables 
     63      !------------------- 
     64      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 
     65      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers 
     66      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       - 
     67      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       - 
     68      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts 
     69      REAL   , POINTER  ::  flagu, flagv                   !    -   - 
     70      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
     71      INTEGER, DIMENSION (2)                ::   kdimsz 
     72      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
     73      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
     74      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
     75      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask 
     76      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
     77      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid 
    6578      !! 
    66       NAMELIST/nambdy/cn_mask, cn_dta_frs_T, cn_dta_frs_U, cn_dta_frs_V,   & 
    67          &            cn_dta_fla_T, cn_dta_fla_U, cn_dta_fla_V,            & 
    68          &            ln_tides, ln_clim, ln_vol, ln_mask,                  & 
    69          &            ln_dyn_fla, ln_dyn_frs, ln_tra_frs,ln_ice_frs,       & 
    70          &            nn_dtactl, nn_rimwidth, nn_volctl 
     79      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             & 
     80         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, & 
     81         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &   
     82#if defined key_lim2 
     83         &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
     84#endif 
     85         &             ln_vol, nn_volctl, nn_rimwidth 
     86      !! 
     87      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
     88                             nbdysegw, jpiwob, jpjwdt, jpjwft,             & 
     89                             nbdysegn, jpjnob, jpindt, jpinft,             & 
     90                             nbdysegs, jpjsob, jpisdt, jpisft 
     91 
    7192      !!---------------------------------------------------------------------- 
    7293 
     94      IF( nn_timing == 1 ) CALL timing_start('bdy_init') 
     95 
     96      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' ) 
     97 
    7398      IF(lwp) WRITE(numout,*) 
    74       IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries' 
     99      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries' 
    75100      IF(lwp) WRITE(numout,*) '~~~~~~~~' 
    76101      ! 
    77       !                                      ! allocate bdy_oce arrays 
    78       IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' ) 
    79       IF( bdy_dta_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate dta arrays' ) 
    80102 
    81103      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   & 
    82          &                               ' and unstructured open boundary condition are not compatible' ) 
    83  
    84       IF( lk_obc      )   CALL ctl_stop( 'Straight open boundaries,',   & 
    85          &                               ' and unstructured open boundaries are not compatible' ) 
    86  
    87       ! --------------------------- 
    88       REWIND( numnam )                    ! Read namelist parameters 
     104         &                               ' and general open boundary condition are not compatible' ) 
     105 
     106      cgrid= (/'t','u','v'/) 
     107 
     108      ! ----------------------------------------- 
     109      ! Initialise and read namelist parameters 
     110      ! ----------------------------------------- 
     111 
     112      nb_bdy            = 0 
     113      ln_coords_file(:) = .false. 
     114      cn_coords_file(:) = '' 
     115      ln_mask_file      = .false. 
     116      cn_mask_file(:)   = '' 
     117      nn_dyn2d(:)       = 0 
     118      nn_dyn2d_dta(:)   = -1  ! uninitialised flag 
     119      nn_dyn3d(:)       = 0 
     120      nn_dyn3d_dta(:)   = -1  ! uninitialised flag 
     121      nn_tra(:)         = 0 
     122      nn_tra_dta(:)     = -1  ! uninitialised flag 
     123#if defined key_lim2 
     124      nn_ice_lim2(:)    = 0 
     125      nn_ice_lim2_dta(:)= -1  ! uninitialised flag 
     126#endif 
     127      ln_vol            = .false. 
     128      nn_volctl         = -1  ! uninitialised flag 
     129      nn_rimwidth(:)    = -1  ! uninitialised flag 
     130 
     131      REWIND( numnam )                     
    89132      READ  ( numnam, nambdy ) 
     133 
     134      ! ----------------------------------------- 
     135      ! Check and write out namelist parameters 
     136      ! ----------------------------------------- 
    90137 
    91138      !                                   ! control prints 
    92139      IF(lwp) WRITE(numout,*) '         nambdy' 
    93140 
    94       !                                         ! check type of data used (nn_dtactl value) 
    95       IF(lwp) WRITE(numout,*) 'nn_dtactl =', nn_dtactl       
    96       IF(lwp) WRITE(numout,*) 
    97       SELECT CASE( nn_dtactl )                   !  
    98       CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    99       CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    100       CASE DEFAULT   ;   CALL ctl_stop( 'nn_dtactl must be 0 or 1' ) 
    101       END SELECT 
    102  
    103       IF(lwp) WRITE(numout,*) 
    104       IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth 
    105  
    106       IF(lwp) WRITE(numout,*) 
    107       IF(lwp) WRITE(numout,*) '      nn_volctl = ', nn_volctl 
    108  
    109       IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    110          SELECT CASE ( nn_volctl ) 
     141      IF( nb_bdy .eq. 0 ) THEN  
     142        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 
     143      ELSE 
     144        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy 
     145      ENDIF 
     146 
     147      DO ib_bdy = 1,nb_bdy 
     148        IF(lwp) WRITE(numout,*) ' '  
     149        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
     150 
     151        IF( ln_coords_file(ib_bdy) ) THEN 
     152           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 
     153        ELSE 
     154           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.' 
     155        ENDIF 
     156        IF(lwp) WRITE(numout,*) 
     157 
     158        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
     159        SELECT CASE( nn_dyn2d(ib_bdy) )                   
     160          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     161          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     162          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
     163          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
     164        END SELECT 
     165        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     166           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
     167              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     168              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     169              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file' 
     170              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
     171              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
     172           END SELECT 
     173        ENDIF 
     174        IF(lwp) WRITE(numout,*) 
     175 
     176        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
     177        SELECT CASE( nn_dyn3d(ib_bdy) )                   
     178          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     179          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     180          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 
     181        END SELECT 
     182        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN 
     183           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
     184              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     185              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     186              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 
     187           END SELECT 
     188        ENDIF 
     189        IF(lwp) WRITE(numout,*) 
     190 
     191        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
     192        SELECT CASE( nn_tra(ib_bdy) )                   
     193          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     194          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     195          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
     196        END SELECT 
     197        IF( nn_tra(ib_bdy) .gt. 0 ) THEN 
     198           SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
     199              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     200              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     201              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
     202           END SELECT 
     203        ENDIF 
     204        IF(lwp) WRITE(numout,*) 
     205 
     206#if defined key_lim2 
     207        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
     208        SELECT CASE( nn_ice_lim2(ib_bdy) )                   
     209          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     210          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     211          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
     212        END SELECT 
     213        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN  
     214           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !  
     215              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     216              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     217              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' ) 
     218           END SELECT 
     219        ENDIF 
     220        IF(lwp) WRITE(numout,*) 
     221#endif 
     222 
     223        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy) 
     224        IF(lwp) WRITE(numout,*) 
     225 
     226      ENDDO 
     227 
     228     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     229       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
     230       IF(lwp) WRITE(numout,*) 
     231       SELECT CASE ( nn_volctl ) 
    111232         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    112233         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    113234         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    114          END SELECT 
    115          IF(lwp) WRITE(numout,*) 
    116       ELSE 
    117          IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries' 
    118          IF(lwp) WRITE(numout,*) 
    119       ENDIF 
    120  
    121       IF( ln_tides ) THEN 
    122         IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 
    123         IF(lwp) WRITE(numout,*) 
    124       ENDIF 
    125  
    126       IF( ln_dyn_fla ) THEN 
    127         IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries' 
    128         IF(lwp) WRITE(numout,*) 
    129       ENDIF 
    130  
    131       IF( ln_dyn_frs ) THEN 
    132         IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries' 
    133         IF(lwp) WRITE(numout,*) 
    134       ENDIF 
    135  
    136       IF( ln_tra_frs ) THEN 
    137         IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries' 
    138         IF(lwp) WRITE(numout,*) 
    139       ENDIF 
    140  
    141       IF( ln_ice_frs ) THEN 
    142         IF(lwp) WRITE(numout,*) 'FRS condition on ice fields at unstructured open boundaries' 
    143         IF(lwp) WRITE(numout,*) 
    144       ENDIF 
    145  
    146       IF( ln_tides )   CALL tide_init      ! Read tides namelist  
    147  
    148  
    149       ! Read arrays defining unstructured open boundaries 
     235       END SELECT 
     236       IF(lwp) WRITE(numout,*) 
     237     ELSE 
     238       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
     239       IF(lwp) WRITE(numout,*) 
     240     ENDIF 
     241 
    150242      ! ------------------------------------------------- 
     243      ! Initialise indices arrays for open boundaries 
     244      ! ------------------------------------------------- 
     245 
     246      ! Work out global dimensions of boundary data 
     247      ! --------------------------------------------- 
     248      REWIND( numnam )                     
     249      DO ib_bdy = 1, nb_bdy 
     250 
     251         jpbdta = 1 
     252         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
     253  
     254            ! No REWIND here because may need to read more than one nambdy_index namelist. 
     255            READ  ( numnam, nambdy_index ) 
     256 
     257            ! Automatic boundary definition: if nbdysegX = -1 
     258            ! set boundary to whole side of model domain. 
     259            IF( nbdysege == -1 ) THEN 
     260               nbdysege = 1 
     261               jpieob(1) = jpiglo - 1 
     262               jpjedt(1) = 2 
     263               jpjeft(1) = jpjglo - 1 
     264            ENDIF 
     265            IF( nbdysegw == -1 ) THEN 
     266               nbdysegw = 1 
     267               jpiwob(1) = 2 
     268               jpjwdt(1) = 2 
     269               jpjwft(1) = jpjglo - 1 
     270            ENDIF 
     271            IF( nbdysegn == -1 ) THEN 
     272               nbdysegn = 1 
     273               jpjnob(1) = jpjglo - 1 
     274               jpindt(1) = 2 
     275               jpinft(1) = jpiglo - 1 
     276            ENDIF 
     277            IF( nbdysegs == -1 ) THEN 
     278               nbdysegs = 1 
     279               jpjsob(1) = 2 
     280               jpisdt(1) = 2 
     281               jpisft(1) = jpiglo - 1 
     282            ENDIF 
     283 
     284            nblendta(:,ib_bdy) = 0 
     285            DO iseg = 1, nbdysege 
     286               igrd = 1 
     287               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
     288               igrd = 2 
     289               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
     290               igrd = 3 
     291               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)                
     292            ENDDO 
     293            DO iseg = 1, nbdysegw 
     294               igrd = 1 
     295               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     296               igrd = 2 
     297               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     298               igrd = 3 
     299               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)                
     300            ENDDO 
     301            DO iseg = 1, nbdysegn 
     302               igrd = 1 
     303               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1                
     304               igrd = 2 
     305               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)                
     306               igrd = 3 
     307               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 
     308            ENDDO 
     309            DO iseg = 1, nbdysegs 
     310               igrd = 1 
     311               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
     312               igrd = 2 
     313               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 
     314               igrd = 3 
     315               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
     316            ENDDO 
     317 
     318            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
     319            jpbdta = MAXVAL(nblendta(:,ib_bdy))                
     320 
     321 
     322         ELSE            ! Read size of arrays in boundary coordinates file. 
     323 
     324 
     325            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
     326            jpbdta = 1 
     327            DO igrd = 1, jpbgrd 
     328               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
     329               nblendta(igrd,ib_bdy) = kdimsz(1) 
     330               jpbdta = MAX(jpbdta, kdimsz(1)) 
     331            ENDDO 
     332 
     333         ENDIF  
     334 
     335      ENDDO ! ib_bdy 
     336 
     337      ! Allocate arrays 
     338      !--------------- 
     339      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
     340         &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     341 
     342      ALLOCATE( dta_global(jpbdta, 1, jpk) ) 
     343 
     344      ! Calculate global boundary index arrays or read in from file 
     345      !------------------------------------------------------------ 
     346      REWIND( numnam )                     
     347      DO ib_bdy = 1, nb_bdy 
     348 
     349         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 
     350 
     351            ! No REWIND here because may need to read more than one nambdy_index namelist. 
     352            READ  ( numnam, nambdy_index ) 
     353 
     354            ! Automatic boundary definition: if nbdysegX = -1 
     355            ! set boundary to whole side of model domain. 
     356            IF( nbdysege == -1 ) THEN 
     357               nbdysege = 1 
     358               jpieob(1) = jpiglo - 1 
     359               jpjedt(1) = 2 
     360               jpjeft(1) = jpjglo - 1 
     361            ENDIF 
     362            IF( nbdysegw == -1 ) THEN 
     363               nbdysegw = 1 
     364               jpiwob(1) = 2 
     365               jpjwdt(1) = 2 
     366               jpjwft(1) = jpjglo - 1 
     367            ENDIF 
     368            IF( nbdysegn == -1 ) THEN 
     369               nbdysegn = 1 
     370               jpjnob(1) = jpjglo - 1 
     371               jpindt(1) = 2 
     372               jpinft(1) = jpiglo - 1 
     373            ENDIF 
     374            IF( nbdysegs == -1 ) THEN 
     375               nbdysegs = 1 
     376               jpjsob(1) = 2 
     377               jpisdt(1) = 2 
     378               jpisft(1) = jpiglo - 1 
     379            ENDIF 
     380 
     381            ! ------------ T points ------------- 
     382            igrd = 1   
     383            icount = 0 
     384            DO ir = 1, nn_rimwidth(ib_bdy) 
     385               ! east 
     386               DO iseg = 1, nbdysege 
     387                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     388                     icount = icount + 1 
     389                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
     390                     nbjdta(icount, igrd, ib_bdy) = ij 
     391                     nbrdta(icount, igrd, ib_bdy) = ir 
     392                  ENDDO 
     393               ENDDO 
     394               ! west 
     395               DO iseg = 1, nbdysegw 
     396                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     397                     icount = icount + 1 
     398                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     399                     nbjdta(icount, igrd, ib_bdy) = ij 
     400                     nbrdta(icount, igrd, ib_bdy) = ir 
     401                  ENDDO 
     402               ENDDO 
     403               ! north 
     404               DO iseg = 1, nbdysegn 
     405                  DO ii = jpindt(iseg), jpinft(iseg) 
     406                     icount = icount + 1 
     407                     nbidta(icount, igrd, ib_bdy) = ii 
     408                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
     409                     nbrdta(icount, igrd, ib_bdy) = ir 
     410                  ENDDO 
     411               ENDDO 
     412               ! south 
     413               DO iseg = 1, nbdysegs 
     414                  DO ii = jpisdt(iseg), jpisft(iseg) 
     415                     icount = icount + 1 
     416                     nbidta(icount, igrd, ib_bdy) = ii 
     417                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 
     418                     nbrdta(icount, igrd, ib_bdy) = ir 
     419                  ENDDO 
     420               ENDDO 
     421            ENDDO 
     422 
     423            ! ------------ U points ------------- 
     424            igrd = 2   
     425            icount = 0 
     426            DO ir = 1, nn_rimwidth(ib_bdy) 
     427               ! east 
     428               DO iseg = 1, nbdysege 
     429                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     430                     icount = icount + 1 
     431                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 
     432                     nbjdta(icount, igrd, ib_bdy) = ij 
     433                     nbrdta(icount, igrd, ib_bdy) = ir 
     434                  ENDDO 
     435               ENDDO 
     436               ! west 
     437               DO iseg = 1, nbdysegw 
     438                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     439                     icount = icount + 1 
     440                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     441                     nbjdta(icount, igrd, ib_bdy) = ij 
     442                     nbrdta(icount, igrd, ib_bdy) = ir 
     443                  ENDDO 
     444               ENDDO 
     445               ! north 
     446               DO iseg = 1, nbdysegn 
     447                  DO ii = jpindt(iseg), jpinft(iseg) - 1 
     448                     icount = icount + 1 
     449                     nbidta(icount, igrd, ib_bdy) = ii 
     450                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
     451                     nbrdta(icount, igrd, ib_bdy) = ir 
     452                  ENDDO 
     453               ENDDO 
     454               ! south 
     455               DO iseg = 1, nbdysegs 
     456                  DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     457                     icount = icount + 1 
     458                     nbidta(icount, igrd, ib_bdy) = ii 
     459                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 
     460                     nbrdta(icount, igrd, ib_bdy) = ir 
     461                  ENDDO 
     462               ENDDO 
     463            ENDDO 
     464 
     465            ! ------------ V points ------------- 
     466            igrd = 3   
     467            icount = 0 
     468            DO ir = 1, nn_rimwidth(ib_bdy) 
     469               ! east 
     470               DO iseg = 1, nbdysege 
     471                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     472                     icount = icount + 1 
     473                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
     474                     nbjdta(icount, igrd, ib_bdy) = ij 
     475                     nbrdta(icount, igrd, ib_bdy) = ir 
     476                  ENDDO 
     477               ENDDO 
     478               ! west 
     479               DO iseg = 1, nbdysegw 
     480                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     481                     icount = icount + 1 
     482                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     483                     nbjdta(icount, igrd, ib_bdy) = ij 
     484                     nbrdta(icount, igrd, ib_bdy) = ir 
     485                  ENDDO 
     486               ENDDO 
     487               ! north 
     488               DO iseg = 1, nbdysegn 
     489                  DO ii = jpindt(iseg), jpinft(iseg) 
     490                     icount = icount + 1 
     491                     nbidta(icount, igrd, ib_bdy) = ii 
     492                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 
     493                     nbrdta(icount, igrd, ib_bdy) = ir 
     494                  ENDDO 
     495               ENDDO 
     496               ! south 
     497               DO iseg = 1, nbdysegs 
     498                  DO ii = jpisdt(iseg), jpisft(iseg) 
     499                     icount = icount + 1 
     500                     nbidta(icount, igrd, ib_bdy) = ii 
     501                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 
     502                     nbrdta(icount, igrd, ib_bdy) = ir 
     503                  ENDDO 
     504               ENDDO 
     505            ENDDO 
     506 
     507         ELSE            ! Read global index arrays from boundary coordinates file. 
     508 
     509            DO igrd = 1, jpbgrd 
     510               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     511               DO ii = 1,nblendta(igrd,ib_bdy) 
     512                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     513               END DO 
     514               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     515               DO ii = 1,nblendta(igrd,ib_bdy) 
     516                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     517               END DO 
     518               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     519               DO ii = 1,nblendta(igrd,ib_bdy) 
     520                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     521               END DO 
     522 
     523               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) ) 
     524               IF(lwp) WRITE(numout,*) 
     525               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 
     526               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 
     527               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
     528                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
     529 
     530            END DO 
     531            CALL iom_close( inum ) 
     532 
     533         ENDIF  
     534 
     535      ENDDO  
     536 
     537      ! Work out dimensions of boundary data on each processor 
     538      ! ------------------------------------------------------ 
     539      
     540      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
     541      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
     542      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
     543      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
     544 
     545      DO ib_bdy = 1, nb_bdy 
     546         DO igrd = 1, jpbgrd 
     547            icount  = 0 
     548            icountr = 0 
     549            idx_bdy(ib_bdy)%nblen(igrd)    = 0 
     550            idx_bdy(ib_bdy)%nblenrim(igrd) = 0 
     551            DO ib = 1, nblendta(igrd,ib_bdy) 
     552               ! check that data is in correct order in file 
     553               ibm1 = MAX(1,ib-1) 
     554               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
     555                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
     556                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', & 
     557                    'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory') 
     558                  ENDIF     
     559               ENDIF 
     560               ! check if point is in local domain 
     561               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
     562                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN 
     563                  ! 
     564                  icount = icount  + 1 
     565                  ! 
     566                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1 
     567               ENDIF 
     568            ENDDO 
     569            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 
     570            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc         
     571         ENDDO  ! igrd 
     572 
     573         ! Allocate index arrays for this boundary set 
     574         !-------------------------------------------- 
     575         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 
     576         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 
     577         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
     578         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
     579         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
     580         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
     581         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 
     582         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 
     583 
     584         ! Dispatch mapping indices and discrete distances on each processor 
     585         ! ----------------------------------------------------------------- 
     586 
     587         DO igrd = 1, jpbgrd 
     588            icount  = 0 
     589            ! Loop on rimwidth to ensure outermost points come first in the local arrays. 
     590            DO ir=1, nn_rimwidth(ib_bdy) 
     591               DO ib = 1, nblendta(igrd,ib_bdy) 
     592                  ! check if point is in local domain and equals ir 
     593                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
     594                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   & 
     595                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
     596                     ! 
     597                     icount = icount  + 1 
     598                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
     599                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
     600                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
     601                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
     602                  ENDIF 
     603               ENDDO 
     604            ENDDO 
     605         ENDDO  
     606 
     607         ! Compute rim weights for FRS scheme 
     608         ! ---------------------------------- 
     609         DO igrd = 1, jpbgrd 
     610            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     611               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
     612               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
     613!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic 
     614!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear 
     615            END DO 
     616         END DO  
     617 
     618      ENDDO 
     619 
     620      ! ------------------------------------------------------ 
     621      ! Initialise masks and find normal/tangential directions 
     622      ! ------------------------------------------------------ 
    151623 
    152624      ! Read global 2D mask at T-points: bdytmask 
    153       ! ***************************************** 
     625      ! ----------------------------------------- 
    154626      ! bdytmask = 1  on the computational domain AND on open boundaries 
    155627      !          = 0  elsewhere    
    … …  
    158630         zmask(         :                ,:) = 0.e0 
    159631         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0           
    160       ELSE IF( ln_mask ) THEN 
    161          CALL iom_open( cn_mask, inum ) 
     632      ELSE IF( ln_mask_file ) THEN 
     633         CALL iom_open( cn_mask_file, inum ) 
    162634         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 
    163635         CALL iom_close( inum ) 
    … …  
    184656 
    185657 
    186       ! Read discrete distance and mapping indices 
    187       ! ****************************************** 
    188       nbidta(:,:) = 0.e0 
    189       nbjdta(:,:) = 0.e0 
    190       nbrdta(:,:) = 0.e0 
    191  
    192       IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 
    193          icount = 0 
    194          DO ir = 1, nn_rimwidth                  ! Define west boundary (from ii=2 to ii=1+nn_rimwidth): 
    195             DO ij = 3, jpjglo-2 
    196                icount = icount + 1 
    197                nbidta(icount,:) = ir + 1 + (jpizoom-1) 
    198                nbjdta(icount,:) = ij     + (jpjzoom-1)  
    199                nbrdta(icount,:) = ir 
    200             END DO 
    201          END DO 
    202          ! 
    203          DO ir = 1, nn_rimwidth                  ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nn_rimwidth): 
    204             DO ij=3,jpjglo-2 
    205                icount = icount + 1 
    206                nbidta(icount,:) = jpiglo-ir + (jpizoom-1) 
    207                nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points 
    208                nbjdta(icount,:) = ij + (jpjzoom-1) 
    209                nbrdta(icount,:) = ir 
    210             END DO 
    211          END DO 
    212          !        
    213       ELSE            ! Read indices and distances in unstructured boundary data files  
    214          ! 
    215          IF( ln_tides ) THEN             ! Read tides input files for preference in case there are no bdydata files 
    216             clfile(4) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 
    217             clfile(5) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 
    218             clfile(6) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 
    219          ENDIF 
    220          IF( ln_dyn_fla .AND. .NOT. ln_tides ) THEN  
    221             clfile(4) = cn_dta_fla_T 
    222             clfile(5) = cn_dta_fla_U 
    223             clfile(6) = cn_dta_fla_V 
    224          ENDIF 
    225  
    226          IF( ln_tra_frs ) THEN  
    227             clfile(1) = cn_dta_frs_T 
    228             IF( .NOT. ln_dyn_frs ) THEN  
    229                clfile(2) = cn_dta_frs_T     ! Dummy read re read T file for sake of 6 files 
    230                clfile(3) = cn_dta_frs_T     ! 
    231             ENDIF 
    232          ENDIF           
    233          IF( ln_dyn_frs ) THEN  
    234             IF( .NOT. ln_tra_frs )   clfile(1) = cn_dta_frs_U      ! Dummy Read  
    235             clfile(2) = cn_dta_frs_U 
    236             clfile(3) = cn_dta_frs_V  
    237          ENDIF 
    238  
    239          !                                   ! how many files are we to read in? 
    240          IF(ln_tides .OR. ln_dyn_fla)   igrd_start = 4 
    241          ! 
    242          IF(ln_tra_frs    ) THEN   ;   igrd_start = 1 
    243          ELSEIF(ln_dyn_frs) THEN   ;   igrd_start = 2 
    244          ENDIF 
    245          ! 
    246          IF( ln_tra_frs   )   igrd_end = 1 
    247          ! 
    248          IF(ln_dyn_fla .OR. ln_tides) THEN   ;   igrd_end = 6 
    249          ELSEIF( ln_dyn_frs             ) THEN   ;   igrd_end = 3 
    250          ENDIF 
    251  
    252          DO igrd = igrd_start, igrd_end 
    253             CALL iom_open( clfile(igrd), inum ) 
    254             id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz )   
    255             IF(lwp) WRITE(numout,*) 'kdimsz : ',kdimsz 
    256             ib_len = kdimsz(1) 
    257             IF( ib_len > jpbdta)   CALL ctl_stop(  'Boundary data array in file too long.',                  & 
    258                 &                                  'File :', TRIM(clfile(igrd)),'increase parameter jpbdta.' ) 
    259  
    260             CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) ) 
    261             DO ii = 1,ib_len 
    262                nbidta(ii,igrd) = INT( zdta(ii,1) ) 
    263             END DO 
    264             CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) ) 
    265             DO ii = 1,ib_len 
    266                nbjdta(ii,igrd) = INT( zdta(ii,1) ) 
    267             END DO 
    268             CALL iom_get( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) ) 
    269             DO ii = 1,ib_len 
    270                nbrdta(ii,igrd) = INT( zdta(ii,1) ) 
    271             END DO 
    272             CALL iom_close( inum ) 
    273  
    274             IF( igrd < 4) THEN            ! Check that rimwidth in file is big enough for Frs case(barotropic is one): 
    275                ibr_max = MAXVAL( nbrdta(:,igrd) ) 
    276                IF(lwp) WRITE(numout,*) 
    277                IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 
    278                IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth 
    279                IF (ibr_max < nn_rimwidth)   CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file' ) 
    280             ENDIF !Check igrd < 4 
    281             ! 
    282          END DO 
    283          ! 
    284       ENDIF  
    285  
    286       ! Dispatch mapping indices and discrete distances on each processor 
    287       ! ***************************************************************** 
    288       
    289       iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
    290       ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
    291       is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
    292       in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
    293  
    294       DO igrd = igrd_start, igrd_end 
    295          icount  = 0 
    296          icountr = 0 
    297          nblen   (igrd) = 0 
    298          nblenrim(igrd) = 0 
    299          nblendta(igrd) = 0 
    300          DO ir=1, nn_rimwidth 
    301             DO ib = 1, jpbdta 
    302                ! check if point is in local domain and equals ir 
    303                IF(  nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND.   & 
    304                   & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND.   & 
    305                   & nbrdta(ib,igrd) == ir  ) THEN 
    306                   ! 
    307                   icount = icount  + 1 
    308                   ! 
    309                   IF( ir == 1 )   icountr = icountr+1 
    310                   IF (icount > jpbdim) THEN 
    311                      IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small' 
    312                      nstop = nstop + 1 
    313                   ELSE 
    314                      nbi(icount, igrd)  = nbidta(ib,igrd)- mig(1)+1 
    315                      nbj(icount, igrd)  = nbjdta(ib,igrd)- mjg(1)+1 
    316                      nbr(icount, igrd)  = nbrdta(ib,igrd) 
    317                      nbmap(icount,igrd) = ib 
    318                   ENDIF             
    319                ENDIF 
    320             END DO 
    321          END DO 
    322          nblenrim(igrd) = icountr !: length of rim boundary data on each proc 
    323          nblen   (igrd) = icount  !: length of boundary data on each proc         
    324       END DO  
    325  
    326       ! Compute rim weights 
    327       ! ------------------- 
    328       DO igrd = igrd_start, igrd_end 
    329          DO ib = 1, nblen(igrd) 
    330             nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 )                     ! tanh formulation 
    331 !           nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr(ib,igrd))/FLOAT(nn_rimwidth))**2      ! quadratic 
    332 !           nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr(ib,igrd))/FLOAT(nn_rimwidth)          ! linear 
    333          END DO 
    334       END DO  
    335     
    336658      ! Mask corrections 
    337659      ! ---------------- 
    … …  
    361683      ! bdy masks and bmask are now set to zero on boundary points: 
    362684      igrd = 1       ! In the free surface case, bmask is at T-points 
    363       DO ib = 1, nblenrim(igrd)      
    364         bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    365       END DO 
     685      DO ib_bdy = 1, nb_bdy 
     686        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)      
     687          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     688        ENDDO 
     689      ENDDO 
    366690      ! 
    367691      igrd = 1 
    368       DO ib = 1, nblenrim(igrd)       
    369         bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    370       END DO 
     692      DO ib_bdy = 1, nb_bdy 
     693        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)       
     694          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     695        ENDDO 
     696      ENDDO 
    371697      ! 
    372698      igrd = 2 
    373       DO ib = 1, nblenrim(igrd) 
    374         bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    375       END DO 
     699      DO ib_bdy = 1, nb_bdy 
     700        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     701          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     702        ENDDO 
     703      ENDDO 
    376704      ! 
    377705      igrd = 3 
    378       DO ib = 1, nblenrim(igrd) 
    379         bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    380       END DO 
     706      DO ib_bdy = 1, nb_bdy 
     707        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     708          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     709        ENDDO 
     710      ENDDO 
    381711 
    382712      ! Lateral boundary conditions 
    … …  
    384714      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 
    385715 
    386       IF( ln_vol .OR. ln_dyn_fla ) THEN      ! Indices and directions of rim velocity components 
    387          ! 
     716      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components 
     717 
     718         idx_bdy(ib_bdy)%flagu(:) = 0.e0 
     719         idx_bdy(ib_bdy)%flagv(:) = 0.e0 
     720         icount = 0  
     721 
    388722         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 
    389723         !flagu =  0 : u is tangential 
    390724         !flagu =  1 : u is normal to the boundary and is direction is inward 
    391          icount = 0  
    392          flagu(:) = 0.e0 
    393   
     725   
    394726         igrd = 2      ! u-component  
    395          DO ib = 1, nblenrim(igrd)   
    396             zefl=bdytmask(nbi(ib,igrd)  , nbj(ib,igrd)) 
    397             zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd)) 
    398             IF( zefl + zwfl ==2 ) THEN 
    399                icount = icount +1 
     727         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     728            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     729            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     730            zefl = bdytmask(nbi  ,nbj) 
     731            zwfl = bdytmask(nbi+1,nbj) 
     732            IF( zefl + zwfl == 2 ) THEN 
     733               icount = icount + 1 
    400734            ELSE 
    401                flagu(ib)=-zefl+zwfl 
     735               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl 
    402736            ENDIF 
    403737         END DO 
    … …  
    406740         !flagv =  0 : u is tangential 
    407741         !flagv =  1 : u is normal to the boundary and is direction is inward 
    408          flagv(:) = 0.e0 
    409742 
    410743         igrd = 3      ! v-component 
    411          DO ib = 1, nblenrim(igrd)   
    412             znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)) 
    413             zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1) 
    414             IF( znfl + zsfl ==2 ) THEN 
     744         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     745            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     746            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     747            znfl = bdytmask(nbi,nbj  ) 
     748            zsfl = bdytmask(nbi,nbj+1) 
     749            IF( znfl + zsfl == 2 ) THEN 
    415750               icount = icount + 1 
    416751            ELSE 
    417                flagv(ib) = -znfl + zsfl 
     752               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl 
    418753            END IF 
    419754         END DO 
    … …  
    422757            IF(lwp) WRITE(numout,*) 
    423758            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   & 
    424                ' are not boundary points. Check nbi, nbj, indices.' 
     759               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy 
    425760            IF(lwp) WRITE(numout,*) ' ========== ' 
    426761            IF(lwp) WRITE(numout,*) 
    … …  
    428763         ENDIF  
    429764     
    430       ENDIF 
     765      ENDDO 
    431766 
    432767      ! Compute total lateral surface for volume correction: 
    … …  
    435770      IF( ln_vol ) THEN   
    436771         igrd = 2      ! Lateral surface at U-points 
    437          DO ib = 1, nblenrim(igrd) 
    438             bdysurftot = bdysurftot + hu     (nbi(ib,igrd)  ,nbj(ib,igrd))                      & 
    439                &                    * e2u    (nbi(ib,igrd)  ,nbj(ib,igrd)) * ABS( flagu(ib) )   & 
    440                &                    * tmask_i(nbi(ib,igrd)  ,nbj(ib,igrd))                      & 
    441                &                    * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd))                    
    442          END DO 
     772         DO ib_bdy = 1, nb_bdy 
     773            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     774               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     775               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     776               flagu => idx_bdy(ib_bdy)%flagu(ib) 
     777               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           & 
     778                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) & 
     779                  &                    * tmask_i(nbi  , nbj)                           & 
     780                  &                    * tmask_i(nbi+1, nbj)                    
     781            ENDDO 
     782         ENDDO 
    443783 
    444784         igrd=3 ! Add lateral surface at V-points 
    445          DO ib = 1, nblenrim(igrd) 
    446             bdysurftot = bdysurftot + hv     (nbi(ib,igrd),nbj(ib,igrd)  )                      & 
    447                &                    * e1v    (nbi(ib,igrd),nbj(ib,igrd)  ) * ABS( flagv(ib) )   & 
    448                &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)  )                      & 
    449                &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1) 
    450          END DO 
     785         DO ib_bdy = 1, nb_bdy 
     786            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     787               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     788               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     789               flagv => idx_bdy(ib_bdy)%flagv(ib) 
     790               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           & 
     791                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) & 
     792                  &                    * tmask_i(nbi, nbj  )                           & 
     793                  &                    * tmask_i(nbi, nbj+1) 
     794            ENDDO 
     795         ENDDO 
    451796         ! 
    452797         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain 
    453798      END IF    
    454  
    455       ! Initialise bdy data arrays 
    456       ! -------------------------- 
    457       tbdy(:,:) = 0.e0 
    458       sbdy(:,:) = 0.e0 
    459       ubdy(:,:) = 0.e0 
    460       vbdy(:,:) = 0.e0 
    461       sshbdy(:) = 0.e0 
    462       ubtbdy(:) = 0.e0 
    463       vbtbdy(:) = 0.e0 
    464 #if defined key_lim2 
    465       frld_bdy(:) = 0.e0 
    466       hicif_bdy(:) = 0.e0 
    467       hsnif_bdy(:) = 0.e0 
    468 #endif 
    469  
    470       ! Read in tidal constituents and adjust for model start time 
    471       ! ---------------------------------------------------------- 
    472       IF( ln_tides )   CALL tide_data 
    473799      ! 
     800      ! Tidy up 
     801      !-------- 
     802      DEALLOCATE(nbidta, nbjdta, nbrdta) 
     803 
     804      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
     805 
    474806   END SUBROUTINE bdy_init 
    475807 
    476808#else 
    477809   !!--------------------------------------------------------------------------------- 
    478    !!   Dummy module                                   NO unstructured open boundaries 
     810   !!   Dummy module                                   NO open boundaries 
    479811   !!--------------------------------------------------------------------------------- 
    480812CONTAINS 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r2528 r3294  
    88   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
     10   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_bdy 
    1213   !!---------------------------------------------------------------------- 
    13    !!   'key_bdy'     Unstructured Open Boundary Condition 
     14   !!   'key_bdy'     Open Boundary Condition 
    1415   !!---------------------------------------------------------------------- 
    1516   !!   PUBLIC 
    16    !!      tide_init     : read of namelist 
    17    !!      tide_data     : read in and initialisation of tidal constituents at boundary 
     17   !!      tide_init     : read of namelist and initialisation of tidal harmonics data 
    1818   !!      tide_update   : calculation of tidal forcing at each timestep 
    1919   !!   PRIVATE 
    … …  
    2424   !!      vset          :/ 
    2525   !!---------------------------------------------------------------------- 
     26   USE timing          ! Timing 
    2627   USE oce             ! ocean dynamics and tracers  
    2728   USE dom_oce         ! ocean space and time domain 
    … …  
    3334   USE bdy_oce         ! ocean open boundary conditions 
    3435   USE daymod          ! calendar 
     36   USE fldread, ONLY: fld_map 
    3537 
    3638   IMPLICIT NONE 
    3739   PRIVATE 
    3840 
    39    PUBLIC   tide_init     ! routine called in bdyini 
    40    PUBLIC   tide_data     ! routine called in bdyini 
     41   PUBLIC   tide_init     ! routine called in nemo_init 
    4142   PUBLIC   tide_update   ! routine called in bdydyn 
    4243 
    43    LOGICAL, PUBLIC            ::   ln_tide_date          !: =T correct tide phases and amplitude for model start date 
    44    INTEGER, PUBLIC, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents 
    45    INTEGER, PUBLIC            ::   ntide                 !: Actual number of tidal constituents 
    46  
    47    CHARACTER(len=80), PUBLIC                         ::   filtide    !: Filename root for tidal input files 
    48    CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used. 
    49  
    50    INTEGER , PUBLIC, DIMENSION(jptides_max) ::   nindx        !: ??? 
    51    REAL(wp), PUBLIC, DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
    52     
    53    REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   ! Tidal constituents : SSH 
    54    REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     ! Tidal constituents : U 
    55    REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     ! Tidal constituents : V 
     44   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     45      INTEGER                                ::   ncpt       !: Actual number of tidal components 
     46      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr) 
     47      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH 
     48      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U 
     49      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V 
     50   END TYPE TIDES_DATA 
     51 
     52   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents 
     53 
     54   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data 
    5655    
    5756   !!---------------------------------------------------------------------- 
    … …  
    6665      !!                    ***  SUBROUTINE tide_init  *** 
    6766      !!                      
    68       !! ** Purpose : - Read in namelist for tides 
    69       !! 
    70       !!---------------------------------------------------------------------- 
    71       INTEGER ::   itide                  ! dummy loop index  
     67      !! ** Purpose : - Read in namelist for tides and initialise external 
     68      !!                tidal harmonics data 
     69      !! 
     70      !!---------------------------------------------------------------------- 
     71      !! namelist variables 
     72      !!------------------- 
     73      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date 
     74      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files 
     75      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used. 
     76      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
     77      !! 
     78      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components 
     79      INTEGER                                   ::   ib_bdy, itide, ib  !: dummy loop indices 
     80      INTEGER                                   ::   inum, igrd 
     81      INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays) 
     82      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file  
     83      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t            
     84      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc 
     85      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data 
     86      !! 
     87      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut    
    7288      !! 
    7389      NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed 
    7490      !!---------------------------------------------------------------------- 
     91 
     92      IF( nn_timing == 1 ) CALL timing_start('tide_init') 
    7593 
    7694      IF(lwp) WRITE(numout,*) 
    … …  
    7896      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
    7997 
    80       ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
    81       ln_tide_date = .false. 
    82       filtide(:) = '' 
    83       tide_speed(:) = 0.0 
    84       tide_cpt(:) = '' 
    85       REWIND( numnam )                                ! Read namelist parameters 
    86       READ  ( numnam, nambdy_tide ) 
    87       !                                               ! Count number of components specified 
    88       ntide = jptides_max 
    89       DO itide = 1, jptides_max 
    90         IF( tide_cpt(itide) == '' ) THEN 
    91            ntide = itide-1 
    92            exit 
    93         ENDIF 
    94       END DO 
    95  
    96       !                                               ! find constituents in standard list 
    97       DO itide = 1, ntide 
    98          nindx(itide) = 0 
    99          IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
    100          IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
    101          IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
    102          IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
    103          IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
    104          IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
    105          IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
    106          IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
    107          IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
    108          IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
    109          IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
    110          IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
    111          IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
    112          IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
    113          IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
    114          IF( nindx(itide) == 0  .AND. lwp ) THEN 
    115             WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
    116             CALL ctl_warn( ctmp1 ) 
    117          ENDIF 
    118       END DO 
    119       !                                               ! Parameter control and print 
    120       IF( ntide < 1 ) THEN  
    121          CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' ) 
    122       ELSE 
    123          IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
    124          IF(lwp) WRITE(numout,*) '             tidal components specified ', ntide 
    125          IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:ntide) 
    126          IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
    127          IF(lwp) WRITE(numout,*) '                ', tide_speed(1:ntide) 
    128       ENDIF 
    129  
    130       ! Initialisation of tidal harmonics arrays 
    131       sshtide(:) = 0.e0 
    132       utide  (:) = 0.e0 
    133       vtide  (:) = 0.e0 
    134       ! 
     98      REWIND(numnam) 
     99      DO ib_bdy = 1, nb_bdy 
     100         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 
     101 
     102            td => tides(ib_bdy) 
     103 
     104            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
     105            ln_tide_date = .false. 
     106            filtide(:) = '' 
     107            tide_speed(:) = 0.0 
     108            tide_cpt(:) = '' 
     109 
     110            ! Don't REWIND here - may need to read more than one of these namelists. 
     111            READ  ( numnam, nambdy_tide ) 
     112            !                                               ! Count number of components specified 
     113            td%ncpt = 0 
     114            DO itide = 1, jptides_max 
     115              IF( tide_cpt(itide) /= '' ) THEN 
     116                 td%ncpt = td%ncpt + 1 
     117              ENDIF 
     118            END DO 
     119 
     120            ! Fill in phase speeds from namelist 
     121            ALLOCATE( td%speed(td%ncpt) ) 
     122            td%speed = tide_speed(1:td%ncpt) 
     123 
     124            ! Find constituents in standard list 
     125            DO itide = 1, td%ncpt 
     126               nindx(itide) = 0 
     127               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
     128               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
     129               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
     130               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
     131               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
     132               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
     133               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
     134               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
     135               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
     136               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
     137               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
     138               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
     139               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
     140               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
     141               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
     142               IF( nindx(itide) == 0  .AND. lwp ) THEN 
     143                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
     144                  CALL ctl_warn( ctmp1 ) 
     145               ENDIF 
     146            END DO 
     147            !                                               ! Parameter control and print 
     148            IF( td%ncpt < 1 ) THEN  
     149               CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' ) 
     150            ELSE 
     151               IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
     152               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt 
     153               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt) 
     154               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
     155               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt) 
     156            ENDIF 
     157 
     158 
     159            ! Allocate space for tidal harmonics data -  
     160            ! get size from OBC data arrays 
     161            ! --------------------------------------- 
     162 
     163            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh )  
     164            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 
     165 
     166            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )  
     167            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 
     168 
     169            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d )  
     170            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 
     171 
     172            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 
     173 
     174 
     175            ! Open files and read in tidal forcing data 
     176            ! ----------------------------------------- 
     177 
     178            DO itide = 1, td%ncpt 
     179               !                                                              ! SSH fields 
     180               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
     181               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     182               CALL iom_open( clfile, inum ) 
     183               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     184               td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) 
     185               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     186               td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     187               CALL iom_close( inum ) 
     188               !                                                              ! U fields 
     189               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
     190               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     191               CALL iom_open( clfile, inum ) 
     192               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     193               td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) 
     194               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     195               td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     196               CALL iom_close( inum ) 
     197               !                                                              ! V fields 
     198               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
     199               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     200               CALL iom_open( clfile, inum ) 
     201               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     202               td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) 
     203               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     204               td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     205               CALL iom_close( inum ) 
     206               ! 
     207            END DO ! end loop on tidal components 
     208 
     209            IF( ln_tide_date ) THEN      ! correct for date factors 
     210 
     211!! used nmonth, nyear and nday from daymod.... 
     212               ! Calculate date corrects for 15 standard consituents 
     213               ! This is the initialisation step, so nday, nmonth, nyear are the  
     214               ! initial date/time of the integration. 
     215                 print *, nday,nmonth,nyear 
     216                 nyear  = int(ndate0 / 10000  )                          ! initial year 
     217                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month 
     218                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100) 
     219 
     220               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
     221 
     222               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
     223 
     224               DO itide = 1, td%ncpt       ! loop on tidal components 
     225                  ! 
     226                  IF( nindx(itide) /= 0 ) THEN 
     227!!gm use rpi  and rad global variable   
     228                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
     229                     z_atde=z_ftc(nindx(itide))*cos(z_arg) 
     230                     z_btde=z_ftc(nindx(itide))*sin(z_arg) 
     231                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   & 
     232                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
     233                  ELSE 
     234                     z_atde = 1.0_wp 
     235                     z_btde = 0.0_wp 
     236                  ENDIF 
     237                  !                                         !  elevation          
     238                  igrd = 1 
     239                  DO ib = 1, ilen0(igrd)                 
     240                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2) 
     241                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1) 
     242                     td%ssh(ib,itide,1) = z1t 
     243                     td%ssh(ib,itide,2) = z2t 
     244                  END DO 
     245                  !                                         !  u        
     246                  igrd = 2 
     247                  DO ib = 1, ilen0(igrd)                 
     248                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2) 
     249                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1) 
     250                     td%u(ib,itide,1) = z1t 
     251                     td%u(ib,itide,2) = z2t 
     252                  END DO 
     253                  !                                         !  v        
     254                  igrd = 3 
     255                  DO ib = 1, ilen0(igrd)                 
     256                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2) 
     257                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1) 
     258                     td%v(ib,itide,1) = z1t 
     259                     td%v(ib,itide,2) = z2t 
     260                  END DO 
     261                  ! 
     262               END DO   ! end loop on tidal components 
     263               ! 
     264            ENDIF ! date correction 
     265            ! 
     266         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
     267         ! 
     268      END DO ! loop on ib_bdy 
     269 
     270      IF( nn_timing == 1 ) CALL timing_stop('tide_init') 
     271 
    135272   END SUBROUTINE tide_init 
    136273 
    137274 
    138    SUBROUTINE tide_data 
    139       !!---------------------------------------------------------------------- 
    140       !!                    ***  SUBROUTINE tide_data  *** 
    141       !!                     
    142       !! ** Purpose : - Read in tidal harmonics data and adjust for the start  
    143       !!                time of the model run.  
    144       !! 
    145       !!---------------------------------------------------------------------- 
    146       INTEGER ::   itide, igrd, ib        ! dummy loop indices 
    147       CHARACTER(len=80) :: clfile         ! full file name for tidal input file  
    148       INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read) 
    149       INTEGER, DIMENSION(6) :: lendta=0   ! length of data in the file (note may be different from nblendta!) 
    150       REAL(wp) ::  z_arg, z_atde, z_btde, z1t, z2t            
    151       REAL(wp), DIMENSION(jpbdta,1) ::   zdta   ! temporary array for data fields 
    152       REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc 
    153       !!------------------------------------------------------------------------------ 
    154  
    155       ! Open files and read in tidal forcing data 
    156       ! ----------------------------------------- 
    157  
    158       ipj = 1 
    159  
    160       DO itide = 1, ntide 
    161          !                                                              ! SSH fields 
    162          clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
    163          IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    164          CALL iom_open( clfile, inum ) 
    165          igrd = 4 
    166          IF( nblendta(igrd) <= 0 ) THEN  
    167             idvar = iom_varid( inum,'z1' ) 
    168             IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 
    169             nblendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
    170             WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd) 
    171          ENDIF 
    172          ipi = nblendta(igrd) 
    173          CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) ) 
    174          DO ib = 1, nblenrim(igrd) 
    175             ssh1(ib,itide) = zdta(nbmap(ib,igrd),1) 
    176          END DO 
    177          CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) ) 
    178          DO ib = 1, nblenrim(igrd) 
    179             ssh2(ib,itide) = zdta(nbmap(ib,igrd),1) 
    180          END DO 
    181          CALL iom_close( inum ) 
    182          ! 
    183          !                                                              ! U fields 
    184          clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
    185          IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    186          CALL iom_open( clfile, inum ) 
    187          igrd = 5 
    188          IF( lendta(igrd) <= 0 ) THEN  
    189             idvar = iom_varid( inum,'u1' ) 
    190             lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
    191             WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd) 
    192          ENDIF 
    193          ipi = lendta(igrd) 
    194          CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) ) 
    195          DO ib = 1, nblenrim(igrd) 
    196             u1(ib,itide) = zdta(nbmap(ib,igrd),1) 
    197          END DO 
    198          CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) ) 
    199          DO ib = 1, nblenrim(igrd) 
    200             u2(ib,itide) = zdta(nbmap(ib,igrd),1) 
    201          END DO 
    202          CALL iom_close( inum ) 
    203          ! 
    204          !                                                              ! V fields 
    205          clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
    206          if(lwp) write(numout,*) 'Reading data from file ', clfile 
    207          CALL iom_open( clfile, inum ) 
    208          igrd = 6 
    209          IF( lendta(igrd) <= 0 ) THEN  
    210             idvar = iom_varid( inum,'v1' ) 
    211             lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
    212             WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd) 
    213          ENDIF 
    214          ipi = lendta(igrd) 
    215          CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) ) 
    216          DO ib = 1, nblenrim(igrd) 
    217             v1(ib,itide) = zdta(nbmap(ib,igrd),1) 
    218          END DO 
    219          CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) ) 
    220          DO ib=1, nblenrim(igrd) 
    221             v2(ib,itide) = zdta(nbmap(ib,igrd),1) 
    222          END DO 
    223          CALL iom_close( inum ) 
    224          ! 
    225       END DO ! end loop on tidal components 
    226  
    227       IF( ln_tide_date ) THEN      ! correct for date factors 
    228  
    229 !! used nmonth, nyear and nday from daymod.... 
    230          ! Calculate date corrects for 15 standard consituents 
    231          ! This is the initialisation step, so nday, nmonth, nyear are the  
    232          ! initial date/time of the integration. 
    233            print *, nday,nmonth,nyear 
    234            nyear  = int(ndate0 / 10000  )                           ! initial year 
    235            nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month 
    236            nday   = int(ndate0 - nyear * 10000 - nmonth * 100) 
    237  
    238          CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
    239  
    240          IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
    241  
    242          DO itide = 1, ntide       ! loop on tidal components 
    243             ! 
    244             IF( nindx(itide) /= 0 ) THEN 
    245 !!gm use rpi  and rad global variable   
    246                z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
    247                z_atde=z_ftc(nindx(itide))*cos(z_arg) 
    248                z_btde=z_ftc(nindx(itide))*sin(z_arg) 
    249                IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide),   & 
    250                   &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
    251             ELSE 
    252                z_atde = 1.0_wp 
    253                z_btde = 0.0_wp 
    254             ENDIF 
    255             !                                         !  elevation          
    256             igrd = 4 
    257             DO ib = 1, nblenrim(igrd)                 
    258                z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide) 
    259                z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide) 
    260                ssh1(ib,itide) = z1t 
    261                ssh2(ib,itide) = z2t 
    262             END DO 
    263             !                                         !  u        
    264             igrd = 5 
    265             DO ib = 1, nblenrim(igrd)                 
    266                z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide) 
    267                z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide) 
    268                u1(ib,itide) = z1t 
    269                u2(ib,itide) = z2t 
    270             END DO 
    271             !                                         !  v        
    272             igrd = 6 
    273             DO ib = 1, nblenrim(igrd)                 
    274                z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide) 
    275                z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide) 
    276                v1(ib,itide) = z1t 
    277                v2(ib,itide) = z2t 
    278             END DO 
    279             ! 
    280          END DO   ! end loop on tidal components 
    281          ! 
    282       ENDIF ! date correction 
    283       ! 
    284    END SUBROUTINE tide_data 
    285  
    286  
    287    SUBROUTINE tide_update ( kt, jit ) 
     275   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 
    288276      !!---------------------------------------------------------------------- 
    289277      !!                 ***  SUBROUTINE tide_update  *** 
    290278      !!                 
    291       !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.  
     279      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
    292280      !!                 
    293281      !!---------------------------------------------------------------------- 
    294       INTEGER, INTENT( in ) ::   kt      ! Main timestep counter 
     282      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter 
    295283!!gm doctor jit ==> kit 
    296       INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option) 
    297       !! 
    298       INTEGER  ::   itide, igrd, ib      ! dummy loop indices 
    299       REAL(wp) ::   z_arg, z_sarg            !             
     284      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices 
     285      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data 
     286      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data 
     287      INTEGER,INTENT(in),OPTIONAL    ::   jit     ! Barotropic timestep counter (for timesplitting option) 
     288      INTEGER,INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     289                                                       ! is present then units = subcycle timesteps. 
     290                                                       ! time_offset = 0 => get data at "now" time level 
     291                                                       ! time_offset = -1 => get data at "before" time level 
     292                                                       ! time_offset = +1 => get data at "after" time level 
     293                                                       ! etc. 
     294      !! 
     295      INTEGER                          :: itide, igrd, ib      ! dummy loop indices 
     296      INTEGER                          :: time_add             ! time offset in units of timesteps 
     297      REAL(wp)                         :: z_arg, z_sarg       
    300298      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
    301299      !!---------------------------------------------------------------------- 
    302300 
     301      IF( nn_timing == 1 ) CALL timing_start('tide_update') 
     302 
     303      time_add = 0 
     304      IF( PRESENT(time_offset) ) THEN 
     305         time_add = time_offset 
     306      ENDIF 
     307          
    303308      ! Note tide phase speeds are in deg/hour, so we need to convert the 
    304309      ! elapsed time in seconds to hours by dividing by 3600.0 
    305       IF( jit == 0 ) THEN   
    306          z_arg = kt * rdt * rad / 3600.0 
    307       ELSE                              ! we are in a barotropic subcycle (for timesplitting option) 
    308 !         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0 
    309          z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
     310      IF( PRESENT(jit) ) THEN   
     311         z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
     312      ELSE                               
     313         z_arg = (kt+time_add) * rdt * rad / 3600.0 
    310314      ENDIF 
    311315 
    312       DO itide = 1, ntide 
    313          z_sarg = z_arg * tide_speed(itide) 
     316      DO itide = 1, td%ncpt 
     317         z_sarg = z_arg * td%speed(itide) 
    314318         z_cost(itide) = COS( z_sarg ) 
    315319         z_sist(itide) = SIN( z_sarg ) 
    316320      END DO 
    317321 
    318       ! summing of tidal constituents into BDY arrays 
    319       sshtide(:) = 0.0 
    320       utide (:) = 0.0 
    321       vtide (:) = 0.0 
    322       ! 
    323       DO itide = 1, ntide 
    324          igrd=4                              ! SSH on tracer grid. 
    325          DO ib = 1, nblenrim(igrd) 
    326             sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide) 
    327             !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide) 
     322      DO itide = 1, td%ncpt 
     323         igrd=1                              ! SSH on tracer grid. 
     324         DO ib = 1, idx%nblenrim(igrd) 
     325            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 
     326            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2) 
    328327         END DO 
    329          igrd=5                              ! U grid 
    330          DO ib=1, nblenrim(igrd) 
    331             utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide) 
    332             !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide) 
     328         igrd=2                              ! U grid 
     329         DO ib=1, idx%nblenrim(igrd) 
     330            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 
     331            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 
    333332         END DO 
    334          igrd=6                              ! V grid 
    335          DO ib=1, nblenrim(igrd) 
    336             vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide) 
    337             !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide) 
     333         igrd=3                              ! V grid 
     334         DO ib=1, idx%nblenrim(igrd) 
     335            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 
     336            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 
    338337         END DO 
    339338      END DO 
     339      ! 
     340      IF( nn_timing == 1 ) CALL timing_stop('tide_update') 
    340341      ! 
    341342   END SUBROUTINE tide_update 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r2528 r3294  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  bdytra  *** 
    4    !! Ocean tracers:   Flow Relaxation Scheme of tracers on each open boundary 
     4   !! Ocean tracers:   Apply boundary conditions for tracers 
    55   !!====================================================================== 
    66   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
    77   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     8   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_bdy 
    … …  
    1112   !!   'key_bdy'                     Unstructured Open Boundary Conditions 
    1213   !!---------------------------------------------------------------------- 
    13    !!   bdy_tra_frs        : Relaxation of tracers on unstructured open boundaries 
     14   !!   bdy_tra            : Apply open boundary conditions to T and S 
     15   !!   bdy_tra_frs        : Apply Flow Relaxation Scheme 
    1416   !!---------------------------------------------------------------------- 
     17   USE timing          ! Timing 
    1518   USE oce             ! ocean dynamics and tracers variables 
    1619   USE dom_oce         ! ocean space and time domain variables  
    1720   USE bdy_oce         ! ocean open boundary conditions 
     21   USE bdydta, ONLY:   bf 
    1822   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1923   USE in_out_manager  ! I/O manager 
    … …  
    2226   PRIVATE 
    2327 
    24    PUBLIC bdy_tra_frs     ! routine called in tranxt.F90  
     28   PUBLIC bdy_tra      ! routine called in tranxt.F90  
    2529 
    2630   !!---------------------------------------------------------------------- 
    … …  
    3135CONTAINS 
    3236 
    33    SUBROUTINE bdy_tra_frs( kt ) 
     37   SUBROUTINE bdy_tra( kt ) 
     38      !!---------------------------------------------------------------------- 
     39      !!                  ***  SUBROUTINE bdy_tra  *** 
     40      !! 
     41      !! ** Purpose : - Apply open boundary conditions for temperature and salinity 
     42      !! 
     43      !!---------------------------------------------------------------------- 
     44      INTEGER, INTENT( in ) :: kt     ! Main time step counter 
     45      !! 
     46      INTEGER               :: ib_bdy ! Loop index 
     47 
     48      DO ib_bdy=1, nb_bdy 
     49 
     50         SELECT CASE( nn_tra(ib_bdy) ) 
     51         CASE(jp_none) 
     52            CYCLE 
     53         CASE(jp_frs) 
     54            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     55         CASE DEFAULT 
     56            CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' ) 
     57         END SELECT 
     58      ENDDO 
     59 
     60   END SUBROUTINE bdy_tra 
     61 
     62   SUBROUTINE bdy_tra_frs( idx, dta, kt ) 
    3463      !!---------------------------------------------------------------------- 
    3564      !!                 ***  SUBROUTINE bdy_tra_frs  *** 
    3665      !!                     
    37       !! ** Purpose : Apply the Flow Relaxation Scheme for tracers in the   
    38       !!              case of unstructured open boundaries. 
     66      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries. 
    3967      !!  
    4068      !! Reference : Engedahl H., 1995, Tellus, 365-382. 
    4169      !!---------------------------------------------------------------------- 
    42       INTEGER, INTENT( in ) ::   kt 
     70      INTEGER,         INTENT(in) ::   kt 
     71      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     72      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    4373      !!  
    4474      REAL(wp) ::   zwgt           ! boundary weight 
    … …  
    4777      !!---------------------------------------------------------------------- 
    4878      ! 
    49       IF(ln_tra_frs) THEN       ! If this is false, then this routine does nothing.  
    50          ! 
    51          IF( kt == nit000 ) THEN 
    52             IF(lwp) WRITE(numout,*) 
    53             IF(lwp) WRITE(numout,*) 'bdy_tra_frs : Flow Relaxation Scheme for tracers' 
    54             IF(lwp) WRITE(numout,*) '~~~~~~~' 
    55          ENDIF 
    56          ! 
    57          igrd = 1                       ! Everything is at T-points here 
    58          DO ib = 1, nblen(igrd) 
    59             DO ik = 1, jpkm1 
    60                ii = nbi(ib,igrd) 
    61                ij = nbj(ib,igrd) 
    62                zwgt = nbw(ib,igrd) 
    63                ta(ii,ij,ik) = ( ta(ii,ij,ik) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)          
    64                sa(ii,ij,ik) = ( sa(ii,ij,ik) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik) 
    65             END DO 
    66          END DO  
    67          ! 
    68          CALL lbc_lnk( ta, 'T', 1. )   ; CALL lbc_lnk( sa, 'T', 1. )    ! Boundary points should be updated 
    69          ! 
    70       ENDIF ! ln_tra_frs 
     79      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_frs') 
     80      ! 
     81      igrd = 1                       ! Everything is at T-points here 
     82      DO ib = 1, idx%nblen(igrd) 
     83         DO ik = 1, jpkm1 
     84            ii = idx%nbi(ib,igrd) 
     85            ij = idx%nbj(ib,igrd) 
     86            zwgt = idx%nbw(ib,igrd) 
     87            tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta%tem(ib,ik) - tsa(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik)          
     88            tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta%sal(ib,ik) - tsa(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik) 
     89         END DO 
     90      END DO  
     91      ! 
     92      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
     93      ! 
     94      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     95      ! 
     96      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_frs') 
    7197      ! 
    7298   END SUBROUTINE bdy_tra_frs 
    … …  
    77103   !!---------------------------------------------------------------------- 
    78104CONTAINS 
    79    SUBROUTINE bdy_tra_frs(kt)      ! Empty routine 
    80       WRITE(*,*) 'bdy_tra_frs: You should not have seen this print! error?', kt 
    81    END SUBROUTINE bdy_tra_frs 
     105   SUBROUTINE bdy_tra(kt)      ! Empty routine 
     106      WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt 
     107   END SUBROUTINE bdy_tra 
    82108#endif 
    83109 

  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r2528 r3294  
    33   !!                       ***  MODULE  bdyvol  *** 
    44   !! Ocean dynamic :  Volume constraint when unstructured boundary  
    5    !!                  and Free surface are used 
     5   !!                  and filtered free surface are used 
    66   !!====================================================================== 
    77   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
    88   !!             -   !  2006-01  (J. Chanut) Bug correction 
    99   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     10   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    1011   !!---------------------------------------------------------------------- 
    1112#if   defined key_bdy   &&   defined key_dynspg_flt 
    … …  
    1415   !!   'key_dynspg_flt'                              filtered free surface 
    1516   !!---------------------------------------------------------------------- 
     17   USE timing          ! Timing 
    1618   USE oce             ! ocean dynamics and tracers  
    1719   USE dom_oce         ! ocean space and time domain  
    … …  
    7173      !! 
    7274      INTEGER  ::   ji, jj, jk, jb, jgrd 
    73       INTEGER  ::   ii, ij 
     75      INTEGER  ::   ib_bdy, ii, ij 
    7476      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst 
     77      TYPE(OBC_INDEX), POINTER :: idx 
    7578      !!----------------------------------------------------------------------------- 
     79 
     80      IF( nn_timing == 1 ) CALL timing_start('bdy_vol') 
    7681 
    7782      IF( ln_vol ) THEN 
    … …  
    9196      ! ------------------------------------------------ 
    9297      zubtpecor = 0.e0 
    93       jgrd = 2                               ! cumulate u component contribution first  
    94       DO jb = 1, nblenrim(jgrd) 
    95          DO jk = 1, jpkm1 
    96             ii = nbi(jb,jgrd) 
    97             ij = nbj(jb,jgrd) 
    98             zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     98      DO ib_bdy = 1, nb_bdy 
     99         idx => idx_bdy(ib_bdy) 
     100 
     101         jgrd = 2                               ! cumulate u component contribution first  
     102         DO jb = 1, idx%nblenrim(jgrd) 
     103            DO jk = 1, jpkm1 
     104               ii = idx%nbi(jb,jgrd) 
     105               ij = idx%nbj(jb,jgrd) 
     106               zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     107            END DO 
    99108         END DO 
    100       END DO 
    101       jgrd = 3                               ! then add v component contribution 
    102       DO jb = 1, nblenrim(jgrd) 
    103          DO jk = 1, jpkm1 
    104             ii = nbi(jb,jgrd) 
    105             ij = nbj(jb,jgrd) 
    106             zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
     109         jgrd = 3                               ! then add v component contribution 
     110         DO jb = 1, idx%nblenrim(jgrd) 
     111            DO jk = 1, jpkm1 
     112               ii = idx%nbi(jb,jgrd) 
     113               ij = idx%nbj(jb,jgrd) 
     114               zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
     115            END DO 
    107116         END DO 
     117 
    108118      END DO 
    109119      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
    … …  
    118128      ! ------------------------------------------------------------- 
    119129      ztranst = 0.e0 
    120       jgrd = 2                               ! correct u component 
    121       DO jb = 1, nblenrim(jgrd) 
    122          DO jk = 1, jpkm1 
    123             ii = nbi(jb,jgrd) 
    124             ij = nbj(jb,jgrd) 
    125             ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 
    126             ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     130      DO ib_bdy = 1, nb_bdy 
     131         idx => idx_bdy(ib_bdy) 
     132 
     133         jgrd = 2                               ! correct u component 
     134         DO jb = 1, idx%nblenrim(jgrd) 
     135            DO jk = 1, jpkm1 
     136               ii = idx%nbi(jb,jgrd) 
     137               ij = idx%nbj(jb,jgrd) 
     138               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk) 
     139               ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     140            END DO 
    127141         END DO 
    128       END DO 
    129       jgrd = 3                              ! correct v component 
    130       DO jb = 1, nblenrim(jgrd) 
    131          DO jk = 1, jpkm1 
    132             ii = nbi(jb,jgrd) 
    133             ij = nbj(jb,jgrd) 
    134             va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
    135             ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
     142         jgrd = 3                              ! correct v component 
     143         DO jb = 1, idx%nblenrim(jgrd) 
     144            DO jk = 1, jpkm1 
     145               ii = idx%nbi(jb,jgrd) 
     146               ij = idx%nbj(jb,jgrd) 
     147               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
     148               ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
     149            END DO 
    136150         END DO 
     151 
    137152      END DO 
    138153      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
    … …  
    149164         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)' 
    150165      END IF  
     166      ! 
     167      IF( nn_timing == 1 ) CALL timing_stop('bdy_vol') 
    151168      ! 
    152169      END IF ! ln_vol 
Note: See TracChangeset for help on using the changeset viewer.