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 3116 for branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Ignore:
Timestamp:
2011-11-15T21:55:40+01:00 (13 years ago)
Author:
cetlod
Message:

dev_NEMO_MERGE_2011: add in changes dev_NOC_UKMO_MERGE developments

Location:
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r2715 r3116  
    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, J. Chanut) 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. 
     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   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp2d_in              ! Damping timescale (days) for 2D solution for inward radiation or FRS  
     82   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp2d_out             ! Damping timescale (days) for 2D solution for outward radiation  
     83   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp3d_in              ! Damping timescale (days) for 3D solution for inward radiation or FRS  
     84   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp3d_out             ! Damping timescale (days) for 3D solution for outward radiation 
    4685 
     86    
    4787   !!---------------------------------------------------------------------- 
    4888   !! Global variables 
    … …  
    5292   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points 
    5393 
     94   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary 
     95 
     96   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !:  
     97   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !:  
     98   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields  
     99   REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !:  
     100   REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !:  
     101 
    54102   !!---------------------------------------------------------------------- 
    55    !! Unstructured open boundary data variables 
     103   !! open boundary data variables 
    56104   !!---------------------------------------------------------------------- 
    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 
    60105 
    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 
     106   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions 
     107                                                                          !: =1 => some data to be read in from data files 
     108   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays 
     109   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
     110   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 
    81111 
    82112   !!---------------------------------------------------------------------- 
    … …  
    94124      !!---------------------------------------------------------------------- 
    95125      ! 
    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 ) 
     126      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                    &   
     127         &      STAT=bdy_oce_alloc ) 
    99128         ! 
    100129      IF( lk_mpp             )   CALL mpp_sum ( bdy_oce_alloc ) 
    … …  
    112141   !!====================================================================== 
    113142END MODULE bdy_oce 
     143 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90

    r2528 r3116  
    1717 
    1818   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 
     19   INTEGER, PUBLIC, PARAMETER ::   jp_bdy  = 10       !: Maximum number of bdy sets 
    2120   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) 
     21   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3        !: Number of horizontal grid types used  (T, U, V) 
     22 
     23   !! Flags for choice of schemes 
     24   INTEGER, PUBLIC, PARAMETER ::   jp_none         = 0        !: Flag for no open boundary condition 
     25   INTEGER, PUBLIC, PARAMETER ::   jp_frs          = 1        !: Flag for Flow Relaxation Scheme 
     26   INTEGER, PUBLIC, PARAMETER ::   jp_flather      = 2        !: Flag for Flather 
    2327#else 
    2428   !!---------------------------------------------------------------------- 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r2977 r3116  
    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  ???????????????? 
    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         
     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 
    1920   !!---------------------------------------------------------------------- 
    2021   USE oce             ! ocean dynamics and tracers 
    2122   USE dom_oce         ! ocean space and time domain 
    2223   USE phycst          ! physical constants 
    23    USE bdy_oce         ! ocean open boundary conditions 
     24   USE bdy_oce         ! ocean open boundary conditions   
    2425   USE bdytides        ! tidal forcing at boundaries 
    25    USE iom 
    26    USE ioipsl 
     26   USE fldread         ! read input fields 
     27   USE iom             ! IOM library 
    2728   USE in_out_manager  ! I/O logical units 
    2829#if defined key_lim2 
    … …  
    3334   PRIVATE 
    3435 
    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  
     36   PUBLIC   bdy_dta          ! routine called by step.F90 and dynspg_ts.F90 
     37   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90 
     38 
     39   INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set. 
     40   INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets. 
     41 
     42   LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions 
     43                                                               ! =F => baroclinic velocities in 3D boundary conditions 
     44 
     45   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read) 
     46 
     47   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap 
     48 
     49#  include "domzgr_substitute.h90" 
    6150   !!---------------------------------------------------------------------- 
    6251   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    6655CONTAINS 
    6756 
    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 ) 
     57      SUBROUTINE bdy_dta( kt, jit, time_offset ) 
    8558      !!---------------------------------------------------------------------- 
    86       !!                   ***  SUBROUTINE bdy_dta_frs  *** 
     59      !!                   ***  SUBROUTINE bdy_dta  *** 
    8760      !!                     
    88       !! ** Purpose :   Read unstructured boundary data for FRS condition. 
     61      !! ** Purpose :   Update external data for open boundary conditions 
    8962      !! 
    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      !! ** Method  :   Use fldread.F90 
     64      !!                 
    9465      !!---------------------------------------------------------------------- 
    95       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index (for timesplitting option, otherwise zero) 
     66      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
     67      USE wrk_nemo, ONLY: wrk_2d_22, wrk_2d_23   ! 2D workspace 
    9668      !! 
    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 
     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      !!--------------------------------------------------------------------------- 
    11883 
    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) 
     84      IF(wrk_in_use(2, 22,23) ) THEN 
     85         CALL ctl_stop('bdy_dta: ERROR: requested workspace arrays are unavailable.')   ;   RETURN 
     86      END IF 
     87 
     88      ! Initialise data arrays once for all from initial conditions where required 
     89      !--------------------------------------------------------------------------- 
     90      IF( kt .eq. nit000 .and. .not. PRESENT(jit) ) THEN 
     91 
     92         ! Calculate depth-mean currents 
     93         !----------------------------- 
     94         pu2d => wrk_2d_22 
     95         pu2d => wrk_2d_23 
     96 
     97         pu2d(:,:) = 0.e0 
     98         pv2d(:,:) = 0.e0 
     99 
     100         DO ik = 1, jpkm1   !! Vertically integrated momentum trends 
     101             pu2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik) 
     102             pv2d(:,:) = pv2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik) 
     103         END DO 
     104         pu2d(:,:) = pu2d(:,:) * hur(:,:) 
     105         pv2d(:,:) = pv2d(:,:) * hvr(:,:) 
     106          
     107         DO ib_bdy = 1, nb_bdy 
     108 
     109            nblen => idx_bdy(ib_bdy)%nblen 
     110            nblenrim => idx_bdy(ib_bdy)%nblenrim 
     111 
     112            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN  
     113               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     114                  ilen1(:) = nblen(:) 
     115               ELSE 
     116                  ilen1(:) = nblenrim(:) 
     117               ENDIF 
     118               igrd = 1 
     119               DO ib = 1, ilen1(igrd) 
     120                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     121                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     122                  dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     123               END DO  
     124               igrd = 2 
     125               DO ib = 1, ilen1(igrd) 
     126                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     127                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     128                  dta_bdy(ib_bdy)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1)          
     129               END DO  
     130               igrd = 3 
     131               DO ib = 1, ilen1(igrd) 
     132                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     133                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     134                  dta_bdy(ib_bdy)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1)          
     135               END DO  
     136            ENDIF 
     137 
     138            IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN  
     139               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     140                  ilen1(:) = nblen(:) 
     141               ELSE 
     142                  ilen1(:) = nblenrim(:) 
     143               ENDIF 
     144               igrd = 2  
     145               DO ib = 1, ilen1(igrd) 
    333146                  DO ik = 1, jpkm1 
    334                      tbdy(ib,ik) = tsn(ii,ij,ik,jp_tem) 
    335                      sbdy(ib,ik) = tsn(ii,ij,ik,jp_sal) 
     147                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     148                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     149                     dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)          
     150                  END DO 
     151               END DO  
     152               igrd = 3  
     153               DO ib = 1, ilen1(igrd) 
     154                  DO ik = 1, jpkm1 
     155                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     156                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     157                     dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik)          
     158                     END DO 
     159               END DO  
     160            ENDIF 
     161 
     162            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN  
     163               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     164                  ilen1(:) = nblen(:) 
     165               ELSE 
     166                  ilen1(:) = nblenrim(:) 
     167               ENDIF 
     168               igrd = 1                       ! Everything is at T-points here 
     169               DO ib = 1, ilen1(igrd) 
     170                  DO ik = 1, jpkm1 
     171                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     172                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     173                     dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
     174                     dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
     175                  END DO 
     176               END DO  
     177            ENDIF 
     178 
     179#if defined key_lim2 
     180            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
     181               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     182                  ilen1(:) = nblen(:) 
     183               ELSE 
     184                  ilen1(:) = nblenrim(:) 
     185               ENDIF 
     186               igrd = 1                       ! Everything is at T-points here 
     187               DO ib = 1, ilen1(igrd) 
     188                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     189                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     190                  dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)          
     191                  dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)          
     192                  dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)          
     193               END DO  
     194            ENDIF 
     195#endif 
     196 
     197         ENDDO ! ib_bdy 
     198 
     199      ENDIF ! kt .eq. nit000 
     200 
     201      ! update external data from files 
     202      !-------------------------------- 
     203      
     204      jstart = 1 
     205      DO ib_bdy = 1, nb_bdy    
     206         IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required 
     207       
     208            IF( PRESENT(jit) ) THEN 
     209               ! Update barotropic boundary conditions only 
     210               ! jit is optional argument for fld_read and tide_update 
     211               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     212                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     213                     dta_bdy(ib_bdy)%ssh(:) = 0.0 
     214                     dta_bdy(ib_bdy)%u2d(:) = 0.0 
     215                     dta_bdy(ib_bdy)%v2d(:) = 0.0 
     216                  ENDIF 
     217                  IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data 
     218                     jend = jstart + 2 
     219                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), 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), jit=jit, time_offset=time_offset ) 
     223                  ENDIF 
     224               ENDIF 
     225            ELSE 
     226               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     227                  dta_bdy(ib_bdy)%ssh(:) = 0.0 
     228                  dta_bdy(ib_bdy)%u2d(:) = 0.0 
     229                  dta_bdy(ib_bdy)%v2d(:) = 0.0 
     230               ENDIF 
     231               IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
     232                  jend = jstart + nb_bdy_fld(ib_bdy) - 1 
     233                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 
     234               ENDIF 
     235               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     236                  CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
     237               ENDIF 
     238            ENDIF 
     239            jstart = jend+1 
     240 
     241            ! If full velocities in boundary data then split into barotropic and baroclinic data 
     242            ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same 
     243            ! time as the dynspg_ts option).  
     244 
     245            IF( ln_full_vel_array(ib_bdy) .and.                                             &  
     246           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN  
     247 
     248               igrd = 2                      ! zonal velocity 
     249               dta_bdy(ib_bdy)%u2d(:) = 0.0 
     250               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     251                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     252                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     253                  DO ik = 1, jpkm1 
     254                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
     255              &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     256                  END DO 
     257                  dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     258                  DO ik = 1, jpkm1 
     259                     dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)  
    336260                  END DO 
    337261               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) 
     262 
     263               igrd = 3                      ! meridional velocity 
     264               dta_bdy(ib_bdy)%v2d(:) = 0.0 
     265               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     266                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     267                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    345268                  DO ik = 1, jpkm1 
    346                      ubdy(ib,ik) = un(ii, ij, ik) 
     269                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
     270              &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     271                  END DO 
     272                  dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     273                  DO ik = 1, jpkm1 
     274                     dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)  
    347275                  END DO 
    348276               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 
     277     
     278            ENDIF 
     279 
     280         END IF ! nn_dta(ib_bdy) = 1 
     281      END DO  ! ib_bdy 
     282 
     283      IF(wrk_not_released(2, 22,23) )    CALL ctl_stop('bdy_dta: ERROR: failed to release workspace arrays.') 
     284 
     285      END SUBROUTINE bdy_dta 
     286 
     287 
     288      SUBROUTINE bdy_dta_init 
     289      !!---------------------------------------------------------------------- 
     290      !!                   ***  SUBROUTINE bdy_dta_init  *** 
     291      !!                     
     292      !! ** Purpose :   Initialise arrays for reading of external data  
     293      !!                for open boundary conditions 
     294      !! 
     295      !! ** Method  :   Use fldread.F90 
     296      !!                 
     297      !!---------------------------------------------------------------------- 
     298      USE dynspg_oce, ONLY: lk_dynspg_ts 
     299      !! 
     300      INTEGER     ::  ib_bdy, jfld, jstart, jend, ierror  ! local indices 
     301      !! 
     302      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files 
     303      CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files 
     304      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data 
     305                                                                ! =F => baroclinic velocities in 3D boundary data 
     306      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays 
     307      INTEGER,              DIMENSION(jpbgrd) ::  ilen0         ! size of local arrays 
     308      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays 
     309      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld 
     310      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V) 
     311      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts 
     312      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures 
     313      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !  
     314      TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
     315#if defined key_lim2 
     316      TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      ! 
     317#endif 
     318      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
     319#if defined key_lim2 
     320      NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 
     321#endif 
     322      NAMELIST/nambdy_dta/ ln_full_vel 
     323      !!--------------------------------------------------------------------------- 
     324 
     325      ! Set nn_dta 
     326      DO ib_bdy = 1, nb_bdy 
     327         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       & 
     328                               ,nn_dyn3d_dta(ib_bdy)       & 
     329                               ,nn_tra_dta(ib_bdy)         & 
     330#if defined key_ice_lim2 
     331                               ,nn_ice_lim2_dta(ib_bdy)    & 
     332#endif 
     333                              ) 
     334         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1 
     335      END DO 
     336 
     337      ! Work out upper bound of how many fields there are to read in and allocate arrays 
     338      ! --------------------------------------------------------------------------- 
     339      ALLOCATE( nb_bdy_fld(nb_bdy) ) 
     340      nb_bdy_fld(:) = 0 
     341      DO ib_bdy = 1, nb_bdy          
     342         IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
     343            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
     344         ENDIF 
     345         IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 
     346            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
     347         ENDIF 
     348         IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN 
     349            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
     350         ENDIF 
     351#if defined key_lim2 
     352         IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN 
     353            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
     354         ENDIF 
     355#endif                
     356      ENDDO             
     357 
     358      nb_bdy_fld_sum = SUM( nb_bdy_fld ) 
     359 
     360      ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror ) 
     361      IF( ierror > 0 ) THEN    
     362         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN   
     363      ENDIF 
     364      ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror ) 
     365      IF( ierror > 0 ) THEN    
     366         CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' )   ;   RETURN   
     367      ENDIF 
     368      ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror ) 
     369      IF( ierror > 0 ) THEN    
     370         CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' )   ;   RETURN   
     371      ENDIF 
     372      ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) )  
     373      ALLOCATE( ibdy(nb_bdy_fld_sum) )  
     374      ALLOCATE( igrid(nb_bdy_fld_sum) )  
     375 
     376      ! Read namelists 
     377      ! -------------- 
     378      REWIND(numnam) 
     379      jfld = 0  
     380      DO ib_bdy = 1, nb_bdy          
     381         IF( nn_dta(ib_bdy) .eq. 1 ) THEN 
     382            ! set file information 
     383            cn_dir = './'        ! directory in which the model is executed 
     384            ln_full_vel = .false. 
     385            ! ... default values (NB: frequency positive => hours, negative => months) 
     386            !                    !  file       ! frequency !  variable        ! time intep !  clim   ! 'yearly' or ! weights  ! rotation  ! 
     387            !                    !  name       !  (hours)  !   name           !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs     ! 
     388            bn_ssh     = FLD_N(  'bdy_ssh'     ,    24     ,  'sossheig'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     389            bn_u2d     = FLD_N(  'bdy_vel2d_u' ,    24     ,  'vobtcrtx'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     390            bn_v2d     = FLD_N(  'bdy_vel2d_v' ,    24     ,  'vobtcrty'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     391            bn_u3d     = FLD_N(  'bdy_vel3d_u' ,    24     ,  'vozocrtx'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     392            bn_v3d     = FLD_N(  'bdy_vel3d_v' ,    24     ,  'vomecrty'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     393            bn_tem     = FLD_N(  'bdy_tem'     ,    24     ,  'votemper'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     394            bn_sal     = FLD_N(  'bdy_sal'     ,    24     ,  'vosaline'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     395#if defined key_lim2 
     396            bn_frld    = FLD_N(  'bdy_frld'    ,    24     ,  'ildsconc'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     397            bn_hicif   = FLD_N(  'bdy_hicif'   ,    24     ,  'iicethic'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     398            bn_hsnif   = FLD_N(  'bdy_hsnif'   ,    24     ,  'isnothic'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        ) 
     399#endif 
     400 
     401            ! Important NOT to rewind here. 
     402            READ( numnam, nambdy_dta ) 
     403 
     404            cn_dir_array(ib_bdy) = cn_dir 
     405            ln_full_vel_array(ib_bdy) = ln_full_vel 
     406 
     407            IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN 
     408               CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',& 
     409            &                  'with dynspg_ts option' )   ;   RETURN   
     410            ENDIF              
     411 
     412            nblen => idx_bdy(ib_bdy)%nblen 
     413            nblenrim => idx_bdy(ib_bdy)%nblenrim 
     414 
     415            ! Only read in necessary fields for this set. 
     416            ! Important that barotropic variables come first. 
     417            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
     418 
     419               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
     420                  jfld = jfld + 1 
     421                  blf_i(jfld) = bn_ssh 
     422                  ibdy(jfld) = ib_bdy 
     423                  igrid(jfld) = 1 
     424                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     425                  ilen3(jfld) = 1 
     426               ENDIF 
     427 
     428               IF( .not. ln_full_vel_array(ib_bdy) ) THEN 
     429 
     430                  jfld = jfld + 1 
     431                  blf_i(jfld) = bn_u2d 
     432                  ibdy(jfld) = ib_bdy 
     433                  igrid(jfld) = 2 
     434                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     435                     ilen1(jfld) = nblen(igrid(jfld)) 
     436                  ELSE 
     437                     ilen1(jfld) = nblenrim(igrid(jfld)) 
     438                  ENDIF 
     439                  ilen3(jfld) = 1 
     440 
     441                  jfld = jfld + 1 
     442                  blf_i(jfld) = bn_v2d 
     443                  ibdy(jfld) = ib_bdy 
     444                  igrid(jfld) = 3 
     445                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     446                     ilen1(jfld) = nblen(igrid(jfld)) 
     447                  ELSE 
     448                     ilen1(jfld) = nblenrim(igrid(jfld)) 
     449                  ENDIF 
     450                  ilen3(jfld) = 1 
     451 
     452               ENDIF 
     453 
     454            ENDIF 
     455 
     456            ! baroclinic velocities 
     457            IF( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) .or. & 
     458           &      ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.  & 
     459           &        ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
     460 
     461               jfld = jfld + 1 
     462               blf_i(jfld) = bn_u3d 
     463               ibdy(jfld) = ib_bdy 
     464               igrid(jfld) = 2 
     465               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     466                  ilen1(jfld) = nblen(igrid(jfld)) 
     467               ELSE 
     468                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     469               ENDIF 
     470               ilen3(jfld) = jpk 
     471 
     472               jfld = jfld + 1 
     473               blf_i(jfld) = bn_v3d 
     474               ibdy(jfld) = ib_bdy 
     475               igrid(jfld) = 3 
     476               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     477                  ilen1(jfld) = nblen(igrid(jfld)) 
     478               ELSE 
     479                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     480               ENDIF 
     481               ilen3(jfld) = jpk 
     482 
     483            ENDIF 
     484 
     485            ! temperature and salinity 
     486            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN 
     487 
     488               jfld = jfld + 1 
     489               blf_i(jfld) = bn_tem 
     490               ibdy(jfld) = ib_bdy 
     491               igrid(jfld) = 1 
     492               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     493                  ilen1(jfld) = nblen(igrid(jfld)) 
     494               ELSE 
     495                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     496               ENDIF 
     497               ilen3(jfld) = jpk 
     498 
     499               jfld = jfld + 1 
     500               blf_i(jfld) = bn_sal 
     501               ibdy(jfld) = ib_bdy 
     502               igrid(jfld) = 1 
     503               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     504                  ilen1(jfld) = nblen(igrid(jfld)) 
     505               ELSE 
     506                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     507               ENDIF 
     508               ilen3(jfld) = jpk 
     509 
     510            ENDIF 
     511 
     512#if defined key_lim2 
     513            ! sea ice 
     514            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
     515 
     516               jfld = jfld + 1 
     517               blf_i(jfld) = bn_frld 
     518               ibdy(jfld) = ib_bdy 
     519               igrid(jfld) = 1 
     520               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     521                  ilen1(jfld) = nblen(igrid(jfld)) 
     522               ELSE 
     523                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     524               ENDIF 
     525               ilen3(jfld) = 1 
     526 
     527               jfld = jfld + 1 
     528               blf_i(jfld) = bn_hicif 
     529               ibdy(jfld) = ib_bdy 
     530               igrid(jfld) = 1 
     531               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     532                  ilen1(jfld) = nblen(igrid(jfld)) 
     533               ELSE 
     534                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     535               ENDIF 
     536               ilen3(jfld) = 1 
     537 
     538               jfld = jfld + 1 
     539               blf_i(jfld) = bn_hsnif 
     540               ibdy(jfld) = ib_bdy 
     541               igrid(jfld) = 1 
     542               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     543                  ilen1(jfld) = nblen(igrid(jfld)) 
     544               ELSE 
     545                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     546               ENDIF 
     547               ilen3(jfld) = 1 
     548 
     549            ENDIF 
     550#endif 
     551            ! Recalculate field counts 
     552            !------------------------- 
     553            nb_bdy_fld_sum = 0 
     554            IF( ib_bdy .eq. 1 ) THEN  
     555               nb_bdy_fld(ib_bdy) = jfld 
     556               nb_bdy_fld_sum     = jfld               
    377557            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. 
     558               nb_bdy_fld(ib_bdy) = jfld - nb_bdy_fld_sum 
     559               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(ib_bdy) 
     560            ENDIF 
     561 
     562         ENDIF ! nn_dta .eq. 1 
     563      ENDDO ! ib_bdy 
     564 
     565 
     566      DO jfld = 1, nb_bdy_fld_sum 
     567         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) 
     568         IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) ) 
     569         nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld)) 
     570      ENDDO 
     571 
     572      ! fill bf with blf_i and control print 
     573      !------------------------------------- 
     574      jstart = 1 
     575      DO ib_bdy = 1, nb_bdy 
     576         jend = jstart + nb_bdy_fld(ib_bdy) - 1 
     577         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta', 'open boundary conditions', 'nambdy_dta' ) 
     578         jstart = jend + 1 
     579      ENDDO 
     580 
     581      ! Initialise local boundary data arrays 
     582      ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later 
     583      ! nn_xxx_dta=1 : point to "fnow" arrays 
     584      !------------------------------------- 
     585 
     586      jfld = 0 
     587      DO ib_bdy=1, nb_bdy 
     588 
     589         nblen => idx_bdy(ib_bdy)%nblen 
     590         nblenrim => idx_bdy(ib_bdy)%nblenrim 
     591 
     592         IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 
     593            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 
     594               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     595                  ilen0(1:3) = nblen(1:3) 
     596               ELSE 
     597                  ilen0(1:3) = nblenrim(1:3) 
     598               ENDIF 
     599               ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) ) 
     600               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
     601               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
    545602            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  
     603               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
     604                  jfld = jfld + 1 
     605                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
     606               ENDIF 
     607               jfld = jfld + 1 
     608               dta_bdy(ib_bdy)%u2d => bf(jfld)%fnow(:,1,1) 
     609               jfld = jfld + 1 
     610               dta_bdy(ib_bdy)%v2d => bf(jfld)%fnow(:,1,1) 
     611            ENDIF 
     612         ENDIF 
     613 
     614         IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
     615            IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     616               ilen0(1:3) = nblen(1:3) 
     617            ELSE 
     618               ilen0(1:3) = nblenrim(1:3) 
     619            ENDIF 
     620            ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 
     621            ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) 
     622         ENDIF 
     623         IF ( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ).or. & 
     624           &  ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.   & 
     625           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
     626            jfld = jfld + 1 
     627            dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 
     628            jfld = jfld + 1 
     629            dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 
     630         ENDIF 
     631 
     632         IF (nn_tra(ib_bdy) .gt. 0) THEN 
     633            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
     634               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     635                  ilen0(1:3) = nblen(1:3) 
     636               ELSE 
     637                  ilen0(1:3) = nblenrim(1:3) 
     638               ENDIF 
     639               ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 
     640               ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) 
     641            ELSE 
     642               jfld = jfld + 1 
     643               dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:) 
     644               jfld = jfld + 1 
     645               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:) 
     646            ENDIF 
     647         ENDIF 
     648 
     649#if defined key_lim2 
     650         IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
     651            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
     652               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     653                  ilen0(1:3) = nblen(1:3) 
     654               ELSE 
     655                  ilen0(1:3) = nblenrim(1:3) 
     656               ENDIF 
     657               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 
     658               ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) 
     659               ALLOCATE( dta_bdy(ib_bdy)%hsnif(ilen0(1)) ) 
     660            ELSE 
     661               jfld = jfld + 1 
     662               dta_bdy(ib_bdy)%frld  => bf(jfld)%fnow(:,1,1) 
     663               jfld = jfld + 1 
     664               dta_bdy(ib_bdy)%hicif => bf(jfld)%fnow(:,1,1) 
     665               jfld = jfld + 1 
     666               dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1) 
     667            ENDIF 
     668         ENDIF 
     669#endif 
     670 
     671      ENDDO ! ib_bdy  
     672 
     673      END SUBROUTINE bdy_dta_init 
    1217674 
    1218675#else 
    1219676   !!---------------------------------------------------------------------- 
    1220    !!   Dummy module                   NO Unstruct Open Boundary Conditions 
     677   !!   Dummy module                   NO Open Boundary Conditions 
    1221678   !!---------------------------------------------------------------------- 
    1222679CONTAINS 
    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 
     680   SUBROUTINE bdy_dta( kt, jit, time_offset ) ! Empty routine 
     681      INTEGER, INTENT( in )           ::   kt     
     682      INTEGER, INTENT( in ), OPTIONAL ::   jit    
     683      INTEGER, INTENT( in ), OPTIONAL ::   time_offset 
     684      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt 
     685   END SUBROUTINE bdy_dta 
     686   SUBROUTINE bdy_dta_init()                  ! Empty routine 
     687      WRITE(*,*) 'bdy_dta_init: You should not have seen this print! error?' 
     688   END SUBROUTINE bdy_dta_init 
    1229689#endif 
    1230690 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r2528 r3116  
    1515   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    1616   !!---------------------------------------------------------------------- 
    17    !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary 
    18    !!   bdy_dyn_fla    : Flather condition for barotropic solution 
     17   !!   bdy_dyn3d        : apply open boundary conditions to baroclinic velocities 
     18   !!   bdy_dyn3d_frs    : apply Flow Relaxation Scheme 
    1919   !!---------------------------------------------------------------------- 
    2020   USE oce             ! ocean dynamics and tracers  
    2121   USE dom_oce         ! ocean space and time domain 
     22   USE dynspg_oce       
    2223   USE bdy_oce         ! ocean open boundary conditions 
    23    USE dynspg_oce      ! for barotropic variables 
    24    USE phycst          ! physical constants 
     24   USE bdydyn2d        ! open boundary conditions for barotropic solution 
     25   USE bdydyn3d        ! open boundary conditions for baroclinic velocities 
    2526   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    26    USE bdytides        ! for tidal harmonic forcing at boundary 
    2727   USE in_out_manager  ! 
    2828 
    … …  
    3030   PRIVATE 
    3131 
    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 
     32   PUBLIC   bdy_dyn     ! routine called in dynspg_flt (if lk_dynspg_flt) or  
     33                        ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp) 
    3634 
     35#  include "domzgr_substitute.h90" 
    3736   !!---------------------------------------------------------------------- 
    3837   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    4241CONTAINS 
    4342 
    44    SUBROUTINE bdy_dyn_frs( kt ) 
     43   SUBROUTINE bdy_dyn( kt, dyn3d_only ) 
    4544      !!---------------------------------------------------------------------- 
    46       !!                  ***  SUBROUTINE bdy_dyn_frs  *** 
     45      !!                  ***  SUBROUTINE bdy_dyn  *** 
    4746      !! 
    48       !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the   
    49       !!                case of unstructured open boundaries. 
     47      !! ** Purpose : - Wrapper routine for bdy_dyn2d and bdy_dyn3d. 
    5048      !! 
    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. 
    5449      !!---------------------------------------------------------------------- 
    55       INTEGER, INTENT( in ) ::   kt   ! Main time step counter 
     50      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
     51      USE wrk_nemo, ONLY: wrk_2d_7, wrk_2d_8      ! 2D workspace 
    5652      !! 
    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 
     53      INTEGER, INTENT( in )           :: kt               ! Main time step counter 
     54      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities 
     55      !! 
     56      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter 
     57      LOGICAL               :: ll_dyn2d, ll_dyn3d   
     58      !! 
    9459 
     60      IF(wrk_in_use(2, 7,8) ) THEN 
     61         CALL ctl_stop('bdy_dyn: ERROR: requested workspace arrays are unavailable.')   ;   RETURN 
     62      END IF 
    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      pu2d => wrk_2d_7 
     79      pv2d => wrk_2d_8 
    13780 
    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 
     81      !------------------------------------------------------- 
     82      ! Split velocities into barotropic and baroclinic parts 
     83      !------------------------------------------------------- 
     84 
     85      pu2d(:,:) = 0.e0 
     86      pv2d(:,:) = 0.e0 
     87      DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     88          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     89          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     90      END DO 
     91      pu2d(:,:) = pu2d(:,:) * phur(:,:) 
     92      pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     93      DO jk = 1 , jpkm1 
     94         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) 
     95         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) 
     96      END DO 
     97 
     98      !------------------------------------------------------- 
     99      ! Apply boundary conditions to barotropic and baroclinic 
     100      ! parts separately 
     101      !------------------------------------------------------- 
     102 
     103      IF( ll_dyn2d ) CALL bdy_dyn2d( kt ) 
     104 
     105      IF( ll_dyn3d ) CALL bdy_dyn3d( kt ) 
     106 
     107      !------------------------------------------------------- 
     108      ! Recombine velocities 
     109      !------------------------------------------------------- 
     110 
     111      DO jk = 1 , jpkm1 
     112         ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk) 
     113         va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk) 
     114      END DO 
     115 
     116      IF(wrk_not_released(2, 7,8) )    CALL ctl_stop('bdy_dyn: ERROR: failed to release workspace arrays.') 
     117 
     118   END SUBROUTINE bdy_dyn 
    181119 
    182120#else 
    … …  
    185123   !!---------------------------------------------------------------------- 
    186124CONTAINS 
    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 
     125   SUBROUTINE bdy_dyn( kt )      ! Empty routine 
     126      WRITE(*,*) 'bdy_dyn: You should not have seen this print! error?', kt 
     127   END SUBROUTINE bdy_dyn 
    194128#endif 
    195129 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r2715 r3116  
    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, J. Chanut) OBC-BDY merge 
     13   !!                 !  --- Renamed bdyini.F90 -> bdyini.F90 --- 
    1214   !!---------------------------------------------------------------------- 
    1315#if defined key_bdy 
    … …  
    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) 
    … …  
    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,                                  & 
     86         &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             & 
     87         &             nn_dmp3d_in, nn_dmp3d_out 
     88      !! 
     89      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
     90                             nbdysegw, jpiwob, jpjwdt, jpjwft,             & 
     91                             nbdysegn, jpjnob, jpindt, jpinft,             & 
     92                             nbdysegs, jpjsob, jpisdt, jpisft 
     93 
    7194      !!---------------------------------------------------------------------- 
    7295 
     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      nn_dmp2d_in(:)    = -1  ! uninitialised flag 
     131      nn_dmp2d_out(:)   = -1  ! uninitialised flag 
     132      nn_dmp3d_in(:)    = -1  ! uninitialised flag 
     133      nn_dmp3d_out(:)   = -1  ! uninitialised flag 
     134 
     135      REWIND( numnam )                     
    89136      READ  ( numnam, nambdy ) 
     137 
     138      ! ----------------------------------------- 
     139      ! Check and write out namelist parameters 
     140      ! ----------------------------------------- 
    90141 
    91142      !                                   ! control prints 
    92143      IF(lwp) WRITE(numout,*) '         nambdy' 
    93144 
    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 ) 
     145      IF( nb_bdy .eq. 0 ) THEN  
     146        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 
     147      ELSE 
     148        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy 
     149      ENDIF 
     150 
     151      DO ib_bdy = 1,nb_bdy 
     152        IF(lwp) WRITE(numout,*) ' '  
     153        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
     154 
     155        IF( ln_coords_file(ib_bdy) ) THEN 
     156           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 
     157        ELSE 
     158           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.' 
     159        ENDIF 
     160        IF(lwp) WRITE(numout,*) 
     161 
     162        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
     163        SELECT CASE( nn_dyn2d(ib_bdy) )                   
     164          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     165          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     166          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
     167          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
     168        END SELECT 
     169        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     170           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
     171              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     172              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     173              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file' 
     174              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
     175              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
     176           END SELECT 
     177        ENDIF 
     178        IF(lwp) WRITE(numout,*) 
     179 
     180        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
     181        SELECT CASE( nn_dyn3d(ib_bdy) )                   
     182          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     183          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     184          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 
     185        END SELECT 
     186        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN 
     187           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
     188              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     189              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     190              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 
     191           END SELECT 
     192        ENDIF 
     193        IF(lwp) WRITE(numout,*) 
     194 
     195        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
     196        SELECT CASE( nn_tra(ib_bdy) )                   
     197          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     198          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     199          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
     200        END SELECT 
     201        IF( nn_tra(ib_bdy) .gt. 0 ) THEN 
     202           SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
     203              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     204              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     205              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
     206           END SELECT 
     207        ENDIF 
     208        IF(lwp) WRITE(numout,*) 
     209 
     210#if defined key_lim2 
     211        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
     212        SELECT CASE( nn_ice_lim2(ib_bdy) )                   
     213          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     214          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     215          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
     216        END SELECT 
     217        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN  
     218           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !  
     219              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     220              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     221              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' ) 
     222           END SELECT 
     223        ENDIF 
     224        IF(lwp) WRITE(numout,*) 
     225#endif 
     226 
     227        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy) 
     228        IF(lwp) WRITE(numout,*) 
     229 
     230      ENDDO 
     231 
     232     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     233       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
     234       IF(lwp) WRITE(numout,*) 
     235       SELECT CASE ( nn_volctl ) 
    111236         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    112237         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    113238         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 
     239       END SELECT 
     240       IF(lwp) WRITE(numout,*) 
     241     ELSE 
     242       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
     243       IF(lwp) WRITE(numout,*) 
     244     ENDIF 
     245 
    150246      ! ------------------------------------------------- 
     247      ! Initialise indices arrays for open boundaries 
     248      ! ------------------------------------------------- 
     249 
     250      ! Work out global dimensions of boundary data 
     251      ! --------------------------------------------- 
     252      REWIND( numnam )                     
     253      DO ib_bdy = 1, nb_bdy 
     254 
     255         jpbdta = 1 
     256         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
     257  
     258            ! No REWIND here because may need to read more than one nambdy_index namelist. 
     259            READ  ( numnam, nambdy_index ) 
     260 
     261            ! Automatic boundary definition: if nbdysegX = -1 
     262            ! set boundary to whole side of model domain. 
     263            IF( nbdysege == -1 ) THEN 
     264               nbdysege = 1 
     265               jpieob(1) = jpiglo - 1 
     266               jpjedt(1) = 2 
     267               jpjeft(1) = jpjglo - 1 
     268            ENDIF 
     269            IF( nbdysegw == -1 ) THEN 
     270               nbdysegw = 1 
     271               jpiwob(1) = 2 
     272               jpjwdt(1) = 2 
     273               jpjwft(1) = jpjglo - 1 
     274            ENDIF 
     275            IF( nbdysegn == -1 ) THEN 
     276               nbdysegn = 1 
     277               jpjnob(1) = jpjglo - 1 
     278               jpindt(1) = 2 
     279               jpinft(1) = jpiglo - 1 
     280            ENDIF 
     281            IF( nbdysegs == -1 ) THEN 
     282               nbdysegs = 1 
     283               jpjsob(1) = 2 
     284               jpisdt(1) = 2 
     285               jpisft(1) = jpiglo - 1 
     286            ENDIF 
     287 
     288            nblendta(:,ib_bdy) = 0 
     289            DO iseg = 1, nbdysege 
     290               igrd = 1 
     291               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
     292               igrd = 2 
     293               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
     294               igrd = 3 
     295               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)                
     296            ENDDO 
     297            DO iseg = 1, nbdysegw 
     298               igrd = 1 
     299               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     300               igrd = 2 
     301               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     302               igrd = 3 
     303               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)                
     304            ENDDO 
     305            DO iseg = 1, nbdysegn 
     306               igrd = 1 
     307               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1                
     308               igrd = 2 
     309               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)                
     310               igrd = 3 
     311               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 
     312            ENDDO 
     313            DO iseg = 1, nbdysegs 
     314               igrd = 1 
     315               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
     316               igrd = 2 
     317               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 
     318               igrd = 3 
     319               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
     320            ENDDO 
     321 
     322            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
     323            jpbdta = MAXVAL(nblendta(:,ib_bdy))                
     324 
     325 
     326         ELSE            ! Read size of arrays in boundary coordinates file. 
     327 
     328 
     329            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
     330            jpbdta = 1 
     331            DO igrd = 1, jpbgrd 
     332               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
     333               nblendta(igrd,ib_bdy) = kdimsz(1) 
     334               jpbdta = MAX(jpbdta, kdimsz(1)) 
     335            ENDDO 
     336 
     337         ENDIF  
     338 
     339      ENDDO ! ib_bdy 
     340 
     341      ! Allocate arrays 
     342      !--------------- 
     343      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
     344         &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     345 
     346      ALLOCATE( dta_global(jpbdta, 1, jpk) ) 
     347 
     348      ! Calculate global boundary index arrays or read in from file 
     349      !------------------------------------------------------------ 
     350      REWIND( numnam )                     
     351      DO ib_bdy = 1, nb_bdy 
     352 
     353         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 
     354 
     355            ! No REWIND here because may need to read more than one nambdy_index namelist. 
     356            READ  ( numnam, nambdy_index ) 
     357 
     358            ! Automatic boundary definition: if nbdysegX = -1 
     359            ! set boundary to whole side of model domain. 
     360            IF( nbdysege == -1 ) THEN 
     361               nbdysege = 1 
     362               jpieob(1) = jpiglo - 1 
     363               jpjedt(1) = 2 
     364               jpjeft(1) = jpjglo - 1 
     365            ENDIF 
     366            IF( nbdysegw == -1 ) THEN 
     367               nbdysegw = 1 
     368               jpiwob(1) = 2 
     369               jpjwdt(1) = 2 
     370               jpjwft(1) = jpjglo - 1 
     371            ENDIF 
     372            IF( nbdysegn == -1 ) THEN 
     373               nbdysegn = 1 
     374               jpjnob(1) = jpjglo - 1 
     375               jpindt(1) = 2 
     376               jpinft(1) = jpiglo - 1 
     377            ENDIF 
     378            IF( nbdysegs == -1 ) THEN 
     379               nbdysegs = 1 
     380               jpjsob(1) = 2 
     381               jpisdt(1) = 2 
     382               jpisft(1) = jpiglo - 1 
     383            ENDIF 
     384 
     385            ! ------------ T points ------------- 
     386            igrd = 1   
     387            icount = 0 
     388            DO ir = 1, nn_rimwidth(ib_bdy) 
     389               ! east 
     390               DO iseg = 1, nbdysege 
     391                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     392                     icount = icount + 1 
     393                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
     394                     nbjdta(icount, igrd, ib_bdy) = ij 
     395                     nbrdta(icount, igrd, ib_bdy) = ir 
     396                  ENDDO 
     397               ENDDO 
     398               ! west 
     399               DO iseg = 1, nbdysegw 
     400                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     401                     icount = icount + 1 
     402                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     403                     nbjdta(icount, igrd, ib_bdy) = ij 
     404                     nbrdta(icount, igrd, ib_bdy) = ir 
     405                  ENDDO 
     406               ENDDO 
     407               ! north 
     408               DO iseg = 1, nbdysegn 
     409                  DO ii = jpindt(iseg), jpinft(iseg) 
     410                     icount = icount + 1 
     411                     nbidta(icount, igrd, ib_bdy) = ii 
     412                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
     413                     nbrdta(icount, igrd, ib_bdy) = ir 
     414                  ENDDO 
     415               ENDDO 
     416               ! south 
     417               DO iseg = 1, nbdysegs 
     418                  DO ii = jpisdt(iseg), jpisft(iseg) 
     419                     icount = icount + 1 
     420                     nbidta(icount, igrd, ib_bdy) = ii 
     421                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 
     422                     nbrdta(icount, igrd, ib_bdy) = ir 
     423                  ENDDO 
     424               ENDDO 
     425            ENDDO 
     426 
     427            ! ------------ U points ------------- 
     428            igrd = 2   
     429            icount = 0 
     430            DO ir = 1, nn_rimwidth(ib_bdy) 
     431               ! east 
     432               DO iseg = 1, nbdysege 
     433                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     434                     icount = icount + 1 
     435                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 
     436                     nbjdta(icount, igrd, ib_bdy) = ij 
     437                     nbrdta(icount, igrd, ib_bdy) = ir 
     438                  ENDDO 
     439               ENDDO 
     440               ! west 
     441               DO iseg = 1, nbdysegw 
     442                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     443                     icount = icount + 1 
     444                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     445                     nbjdta(icount, igrd, ib_bdy) = ij 
     446                     nbrdta(icount, igrd, ib_bdy) = ir 
     447                  ENDDO 
     448               ENDDO 
     449               ! north 
     450               DO iseg = 1, nbdysegn 
     451                  DO ii = jpindt(iseg), jpinft(iseg) - 1 
     452                     icount = icount + 1 
     453                     nbidta(icount, igrd, ib_bdy) = ii 
     454                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
     455                     nbrdta(icount, igrd, ib_bdy) = ir 
     456                  ENDDO 
     457               ENDDO 
     458               ! south 
     459               DO iseg = 1, nbdysegs 
     460                  DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     461                     icount = icount + 1 
     462                     nbidta(icount, igrd, ib_bdy) = ii 
     463                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 
     464                     nbrdta(icount, igrd, ib_bdy) = ir 
     465                  ENDDO 
     466               ENDDO 
     467            ENDDO 
     468 
     469            ! ------------ V points ------------- 
     470            igrd = 3   
     471            icount = 0 
     472            DO ir = 1, nn_rimwidth(ib_bdy) 
     473               ! east 
     474               DO iseg = 1, nbdysege 
     475                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     476                     icount = icount + 1 
     477                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
     478                     nbjdta(icount, igrd, ib_bdy) = ij 
     479                     nbrdta(icount, igrd, ib_bdy) = ir 
     480                  ENDDO 
     481               ENDDO 
     482               ! west 
     483               DO iseg = 1, nbdysegw 
     484                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     485                     icount = icount + 1 
     486                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     487                     nbjdta(icount, igrd, ib_bdy) = ij 
     488                     nbrdta(icount, igrd, ib_bdy) = ir 
     489                  ENDDO 
     490               ENDDO 
     491               ! north 
     492               DO iseg = 1, nbdysegn 
     493                  DO ii = jpindt(iseg), jpinft(iseg) 
     494                     icount = icount + 1 
     495                     nbidta(icount, igrd, ib_bdy) = ii 
     496                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 
     497                     nbrdta(icount, igrd, ib_bdy) = ir 
     498                  ENDDO 
     499               ENDDO 
     500               ! south 
     501               DO iseg = 1, nbdysegs 
     502                  DO ii = jpisdt(iseg), jpisft(iseg) 
     503                     icount = icount + 1 
     504                     nbidta(icount, igrd, ib_bdy) = ii 
     505                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 
     506                     nbrdta(icount, igrd, ib_bdy) = ir 
     507                  ENDDO 
     508               ENDDO 
     509            ENDDO 
     510 
     511         ELSE            ! Read global index arrays from boundary coordinates file. 
     512 
     513            DO igrd = 1, jpbgrd 
     514               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     515               DO ii = 1,nblendta(igrd,ib_bdy) 
     516                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     517               END DO 
     518               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     519               DO ii = 1,nblendta(igrd,ib_bdy) 
     520                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     521               END DO 
     522               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     523               DO ii = 1,nblendta(igrd,ib_bdy) 
     524                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     525               END DO 
     526 
     527               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) ) 
     528               IF(lwp) WRITE(numout,*) 
     529               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 
     530               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 
     531               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
     532                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
     533 
     534            END DO 
     535            CALL iom_close( inum ) 
     536 
     537         ENDIF  
     538 
     539      ENDDO  
     540 
     541      ! Work out dimensions of boundary data on each processor 
     542      ! ------------------------------------------------------ 
     543      
     544      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
     545      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
     546      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
     547      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
     548 
     549      DO ib_bdy = 1, nb_bdy 
     550         DO igrd = 1, jpbgrd 
     551            icount  = 0 
     552            icountr = 0 
     553            idx_bdy(ib_bdy)%nblen(igrd)    = 0 
     554            idx_bdy(ib_bdy)%nblenrim(igrd) = 0 
     555            DO ib = 1, nblendta(igrd,ib_bdy) 
     556               ! check that data is in correct order in file 
     557               ibm1 = MAX(1,ib-1) 
     558               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
     559                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
     560                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', & 
     561                                   'A utility for re-ordering boundary coordinates and data files exists in CDFTOOLS') 
     562                  ENDIF     
     563               ENDIF 
     564               ! check if point is in local domain 
     565               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
     566                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN 
     567                  ! 
     568                  icount = icount  + 1 
     569                  ! 
     570                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1 
     571               ENDIF 
     572            ENDDO 
     573            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 
     574            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc         
     575         ENDDO  ! igrd 
     576 
     577         ! Allocate index arrays for this boundary set 
     578         !-------------------------------------------- 
     579         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 
     580         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 
     581         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
     582         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
     583         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
     584         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
     585         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 
     586         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 
     587 
     588         ! Dispatch mapping indices and discrete distances on each processor 
     589         ! ----------------------------------------------------------------- 
     590 
     591         DO igrd = 1, jpbgrd 
     592            icount  = 0 
     593            ! Loop on rimwidth to ensure outermost points come first in the local arrays. 
     594            DO ir=1, nn_rimwidth(ib_bdy) 
     595               DO ib = 1, nblendta(igrd,ib_bdy) 
     596                  ! check if point is in local domain and equals ir 
     597                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
     598                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   & 
     599                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
     600                     ! 
     601                     icount = icount  + 1 
     602                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
     603                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
     604                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
     605                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
     606                  ENDIF 
     607               ENDDO 
     608            ENDDO 
     609         ENDDO  
     610 
     611         ! Compute rim weights for FRS scheme 
     612         ! ---------------------------------- 
     613         DO igrd = 1, jpbgrd 
     614            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     615               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
     616               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
     617!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic 
     618!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear 
     619            END DO 
     620         END DO  
     621 
     622      ENDDO 
     623 
     624      ! ------------------------------------------------------ 
     625      ! Initialise masks and find normal/tangential directions 
     626      ! ------------------------------------------------------ 
    151627 
    152628      ! Read global 2D mask at T-points: bdytmask 
    153       ! ***************************************** 
     629      ! ----------------------------------------- 
    154630      ! bdytmask = 1  on the computational domain AND on open boundaries 
    155631      !          = 0  elsewhere    
    … …  
    158634         zmask(         :                ,:) = 0.e0 
    159635         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0           
    160       ELSE IF( ln_mask ) THEN 
    161          CALL iom_open( cn_mask, inum ) 
     636      ELSE IF( ln_mask_file ) THEN 
     637         CALL iom_open( cn_mask_file, inum ) 
    162638         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 
    163639         CALL iom_close( inum ) 
    … …  
    184660 
    185661 
    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     
    336662      ! Mask corrections 
    337663      ! ---------------- 
    … …  
    361687      ! bdy masks and bmask are now set to zero on boundary points: 
    362688      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 
     689      DO ib_bdy = 1, nb_bdy 
     690        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)      
     691          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     692        ENDDO 
     693      ENDDO 
    366694      ! 
    367695      igrd = 1 
    368       DO ib = 1, nblenrim(igrd)       
    369         bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    370       END DO 
     696      DO ib_bdy = 1, nb_bdy 
     697        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)       
     698          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     699        ENDDO 
     700      ENDDO 
    371701      ! 
    372702      igrd = 2 
    373       DO ib = 1, nblenrim(igrd) 
    374         bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    375       END DO 
     703      DO ib_bdy = 1, nb_bdy 
     704        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     705          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     706        ENDDO 
     707      ENDDO 
    376708      ! 
    377709      igrd = 3 
    378       DO ib = 1, nblenrim(igrd) 
    379         bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
    380       END DO 
     710      DO ib_bdy = 1, nb_bdy 
     711        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     712          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 
     713        ENDDO 
     714      ENDDO 
    381715 
    382716      ! Lateral boundary conditions 
    … …  
    384718      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 
    385719 
    386       IF( ln_vol .OR. ln_dyn_fla ) THEN      ! Indices and directions of rim velocity components 
    387          ! 
     720      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components 
     721 
     722         idx_bdy(ib_bdy)%flagu(:) = 0.e0 
     723         idx_bdy(ib_bdy)%flagv(:) = 0.e0 
     724         icount = 0  
     725 
    388726         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 
    389727         !flagu =  0 : u is tangential 
    390728         !flagu =  1 : u is normal to the boundary and is direction is inward 
    391          icount = 0  
    392          flagu(:) = 0.e0 
    393   
     729   
    394730         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 
     731         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     732            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     733            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     734            zefl = bdytmask(nbi  ,nbj) 
     735            zwfl = bdytmask(nbi+1,nbj) 
     736            IF( zefl + zwfl == 2 ) THEN 
     737               icount = icount + 1 
    400738            ELSE 
    401                flagu(ib)=-zefl+zwfl 
     739               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl 
    402740            ENDIF 
    403741         END DO 
    … …  
    406744         !flagv =  0 : u is tangential 
    407745         !flagv =  1 : u is normal to the boundary and is direction is inward 
    408          flagv(:) = 0.e0 
    409746 
    410747         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 
     748         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     749            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     750            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     751            znfl = bdytmask(nbi,nbj  ) 
     752            zsfl = bdytmask(nbi,nbj+1) 
     753            IF( znfl + zsfl == 2 ) THEN 
    415754               icount = icount + 1 
    416755            ELSE 
    417                flagv(ib) = -znfl + zsfl 
     756               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl 
    418757            END IF 
    419758         END DO 
    … …  
    422761            IF(lwp) WRITE(numout,*) 
    423762            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   & 
    424                ' are not boundary points. Check nbi, nbj, indices.' 
     763               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy 
    425764            IF(lwp) WRITE(numout,*) ' ========== ' 
    426765            IF(lwp) WRITE(numout,*) 
    … …  
    428767         ENDIF  
    429768     
    430       ENDIF 
     769      ENDDO 
    431770 
    432771      ! Compute total lateral surface for volume correction: 
    … …  
    435774      IF( ln_vol ) THEN   
    436775         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 
     776         DO ib_bdy = 1, nb_bdy 
     777            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     778               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     779               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     780               flagu => idx_bdy(ib_bdy)%flagu(ib) 
     781               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           & 
     782                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) & 
     783                  &                    * tmask_i(nbi  , nbj)                           & 
     784                  &                    * tmask_i(nbi+1, nbj)                    
     785            ENDDO 
     786         ENDDO 
    443787 
    444788         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 
     789         DO ib_bdy = 1, nb_bdy 
     790            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     791               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     792               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     793               flagv => idx_bdy(ib_bdy)%flagv(ib) 
     794               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           & 
     795                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) & 
     796                  &                    * tmask_i(nbi, nbj  )                           & 
     797                  &                    * tmask_i(nbi, nbj+1) 
     798            ENDDO 
     799         ENDDO 
    451800         ! 
    452801         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain 
    453802      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 
    473803      ! 
     804      ! Tidy up 
     805      !-------- 
     806      DEALLOCATE(nbidta, nbjdta, nbrdta) 
     807 
    474808   END SUBROUTINE bdy_init 
    475809 
    476810#else 
    477811   !!--------------------------------------------------------------------------------- 
    478    !!   Dummy module                                   NO unstructured open boundaries 
     812   !!   Dummy module                                   NO open boundaries 
    479813   !!--------------------------------------------------------------------------------- 
    480814CONTAINS 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r2528 r3116  
    1111#if defined key_bdy 
    1212   !!---------------------------------------------------------------------- 
    13    !!   'key_bdy'     Unstructured Open Boundary Condition 
     13   !!   'key_bdy'     Open Boundary Condition 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   PUBLIC 
    16    !!      tide_init     : read of namelist 
    17    !!      tide_data     : read in and initialisation of tidal constituents at boundary 
     16   !!      tide_init     : read of namelist and initialisation of tidal harmonics data 
    1817   !!      tide_update   : calculation of tidal forcing at each timestep 
    1918   !!   PRIVATE 
    … …  
    3332   USE bdy_oce         ! ocean open boundary conditions 
    3433   USE daymod          ! calendar 
     34   USE fldread, ONLY: fld_map 
    3535 
    3636   IMPLICIT NONE 
    3737   PRIVATE 
    3838 
    39    PUBLIC   tide_init     ! routine called in bdyini 
    40    PUBLIC   tide_data     ! routine called in bdyini 
     39   PUBLIC   tide_init     ! routine called in nemo_init 
    4140   PUBLIC   tide_update   ! routine called in bdydyn 
    4241 
    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 
     42   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     43      INTEGER                                ::   ncpt       !: Actual number of tidal components 
     44      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr) 
     45      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH 
     46      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U 
     47      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V 
     48   END TYPE TIDES_DATA 
     49 
     50   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents 
     51 
     52   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data 
    5653    
    5754   !!---------------------------------------------------------------------- 
    … …  
    6663      !!                    ***  SUBROUTINE tide_init  *** 
    6764      !!                      
    68       !! ** Purpose : - Read in namelist for tides 
    69       !! 
    70       !!---------------------------------------------------------------------- 
    71       INTEGER ::   itide                  ! dummy loop index  
     65      !! ** Purpose : - Read in namelist for tides and initialise external 
     66      !!                tidal harmonics data 
     67      !! 
     68      !!---------------------------------------------------------------------- 
     69      !! namelist variables 
     70      !!------------------- 
     71      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date 
     72      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files 
     73      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used. 
     74      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
     75      !! 
     76      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components 
     77      INTEGER                                   ::   ib_bdy, itide, ib  !: dummy loop indices 
     78      INTEGER                                   ::   inum, igrd 
     79      INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays) 
     80      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file  
     81      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t            
     82      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc 
     83      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data 
     84      !! 
     85      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut    
    7286      !! 
    7387      NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed 
    … …  
    7892      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
    7993 
    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       ! 
     94      REWIND(numnam) 
     95      DO ib_bdy = 1, nb_bdy 
     96         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 
     97 
     98            td => tides(ib_bdy) 
     99 
     100            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
     101            ln_tide_date = .false. 
     102            filtide(:) = '' 
     103            tide_speed(:) = 0.0 
     104            tide_cpt(:) = '' 
     105 
     106            ! Don't REWIND here - may need to read more than one of these namelists. 
     107            READ  ( numnam, nambdy_tide ) 
     108            !                                               ! Count number of components specified 
     109            td%ncpt = 0 
     110            DO itide = 1, jptides_max 
     111              IF( tide_cpt(itide) /= '' ) THEN 
     112                 td%ncpt = td%ncpt + 1 
     113              ENDIF 
     114            END DO 
     115 
     116            ! Fill in phase speeds from namelist 
     117            ALLOCATE( td%speed(td%ncpt) ) 
     118            td%speed = tide_speed(1:td%ncpt) 
     119 
     120            ! Find constituents in standard list 
     121            DO itide = 1, td%ncpt 
     122               nindx(itide) = 0 
     123               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
     124               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
     125               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
     126               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
     127               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
     128               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
     129               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
     130               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
     131               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
     132               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
     133               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
     134               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
     135               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
     136               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
     137               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
     138               IF( nindx(itide) == 0  .AND. lwp ) THEN 
     139                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
     140                  CALL ctl_warn( ctmp1 ) 
     141               ENDIF 
     142            END DO 
     143            !                                               ! Parameter control and print 
     144            IF( td%ncpt < 1 ) THEN  
     145               CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' ) 
     146            ELSE 
     147               IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
     148               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt 
     149               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt) 
     150               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
     151               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt) 
     152            ENDIF 
     153 
     154 
     155            ! Allocate space for tidal harmonics data -  
     156            ! get size from OBC data arrays 
     157            ! --------------------------------------- 
     158 
     159            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh )  
     160            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 
     161 
     162            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )  
     163            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 
     164 
     165            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d )  
     166            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 
     167 
     168            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 
     169 
     170 
     171            ! Open files and read in tidal forcing data 
     172            ! ----------------------------------------- 
     173 
     174            DO itide = 1, td%ncpt 
     175               !                                                              ! SSH fields 
     176               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
     177               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     178               CALL iom_open( clfile, inum ) 
     179               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     180               td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) 
     181               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     182               td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     183               CALL iom_close( inum ) 
     184               !                                                              ! U fields 
     185               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
     186               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     187               CALL iom_open( clfile, inum ) 
     188               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     189               td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) 
     190               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     191               td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     192               CALL iom_close( inum ) 
     193               !                                                              ! V fields 
     194               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
     195               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     196               CALL iom_open( clfile, inum ) 
     197               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     198               td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) 
     199               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     200               td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     201               CALL iom_close( inum ) 
     202               ! 
     203            END DO ! end loop on tidal components 
     204 
     205            IF( ln_tide_date ) THEN      ! correct for date factors 
     206 
     207!! used nmonth, nyear and nday from daymod.... 
     208               ! Calculate date corrects for 15 standard consituents 
     209               ! This is the initialisation step, so nday, nmonth, nyear are the  
     210               ! initial date/time of the integration. 
     211                 print *, nday,nmonth,nyear 
     212                 nyear  = int(ndate0 / 10000  )                          ! initial year 
     213                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month 
     214                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100) 
     215 
     216               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
     217 
     218               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
     219 
     220               DO itide = 1, td%ncpt       ! loop on tidal components 
     221                  ! 
     222                  IF( nindx(itide) /= 0 ) THEN 
     223!!gm use rpi  and rad global variable   
     224                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
     225                     z_atde=z_ftc(nindx(itide))*cos(z_arg) 
     226                     z_btde=z_ftc(nindx(itide))*sin(z_arg) 
     227                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   & 
     228                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
     229                  ELSE 
     230                     z_atde = 1.0_wp 
     231                     z_btde = 0.0_wp 
     232                  ENDIF 
     233                  !                                         !  elevation          
     234                  igrd = 1 
     235                  DO ib = 1, ilen0(igrd)                 
     236                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2) 
     237                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1) 
     238                     td%ssh(ib,itide,1) = z1t 
     239                     td%ssh(ib,itide,2) = z2t 
     240                  END DO 
     241                  !                                         !  u        
     242                  igrd = 2 
     243                  DO ib = 1, ilen0(igrd)                 
     244                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2) 
     245                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1) 
     246                     td%u(ib,itide,1) = z1t 
     247                     td%u(ib,itide,2) = z2t 
     248                  END DO 
     249                  !                                         !  v        
     250                  igrd = 3 
     251                  DO ib = 1, ilen0(igrd)                 
     252                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2) 
     253                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1) 
     254                     td%v(ib,itide,1) = z1t 
     255                     td%v(ib,itide,2) = z2t 
     256                  END DO 
     257                  ! 
     258               END DO   ! end loop on tidal components 
     259               ! 
     260            ENDIF ! date correction 
     261            ! 
     262         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
     263         ! 
     264      END DO ! loop on ib_bdy 
     265 
    135266   END SUBROUTINE tide_init 
    136267 
    137268 
    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 ) 
     269   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 
    288270      !!---------------------------------------------------------------------- 
    289271      !!                 ***  SUBROUTINE tide_update  *** 
    290272      !!                 
    291       !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.  
     273      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
    292274      !!                 
    293275      !!---------------------------------------------------------------------- 
    294       INTEGER, INTENT( in ) ::   kt      ! Main timestep counter 
     276      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter 
    295277!!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            !             
     278      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices 
     279      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data 
     280      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data 
     281      INTEGER,INTENT(in),OPTIONAL    ::   jit     ! Barotropic timestep counter (for timesplitting option) 
     282      INTEGER,INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     283                                                       ! is present then units = subcycle timesteps. 
     284                                                       ! time_offset = 0 => get data at "now" time level 
     285                                                       ! time_offset = -1 => get data at "before" time level 
     286                                                       ! time_offset = +1 => get data at "after" time level 
     287                                                       ! etc. 
     288      !! 
     289      INTEGER                          :: itide, igrd, ib      ! dummy loop indices 
     290      INTEGER                          :: time_add             ! time offset in units of timesteps 
     291      REAL(wp)                         :: z_arg, z_sarg       
    300292      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
    301293      !!---------------------------------------------------------------------- 
    302294 
     295      time_add = 0 
     296      IF( PRESENT(time_offset) ) THEN 
     297         time_add = time_offset 
     298      ENDIF 
     299          
    303300      ! Note tide phase speeds are in deg/hour, so we need to convert the 
    304301      ! 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 
     302      IF( PRESENT(jit) ) THEN   
     303         z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
     304      ELSE                               
     305         z_arg = (kt+time_add) * rdt * rad / 3600.0 
    310306      ENDIF 
    311307 
    312       DO itide = 1, ntide 
    313          z_sarg = z_arg * tide_speed(itide) 
     308      DO itide = 1, td%ncpt 
     309         z_sarg = z_arg * td%speed(itide) 
    314310         z_cost(itide) = COS( z_sarg ) 
    315311         z_sist(itide) = SIN( z_sarg ) 
    316312      END DO 
    317313 
    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) 
     314      DO itide = 1, td%ncpt 
     315         igrd=1                              ! SSH on tracer grid. 
     316         DO ib = 1, idx%nblenrim(igrd) 
     317            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 
     318            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2) 
    328319         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) 
     320         igrd=2                              ! U grid 
     321         DO ib=1, idx%nblenrim(igrd) 
     322            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 
     323            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 
    333324         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) 
     325         igrd=3                              ! V grid 
     326         DO ib=1, idx%nblenrim(igrd) 
     327            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 
     328            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 
    338329         END DO 
    339330      END DO 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r2977 r3116  
    1111   !!   'key_bdy'                     Unstructured Open Boundary Conditions 
    1212   !!---------------------------------------------------------------------- 
    13    !!   bdy_tra_frs        : Relaxation of tracers on unstructured open boundaries 
     13   !!   bdy_tra            : Apply open boundary conditions to T and S 
     14   !!   bdy_tra_frs        : Apply Flow Relaxation Scheme 
    1415   !!---------------------------------------------------------------------- 
    1516   USE oce             ! ocean dynamics and tracers variables 
    1617   USE dom_oce         ! ocean space and time domain variables  
    1718   USE bdy_oce         ! ocean open boundary conditions 
     19   USE bdydta, ONLY:   bf 
    1820   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1921   USE in_out_manager  ! I/O manager 
    … …  
    2224   PRIVATE 
    2325 
    24    PUBLIC bdy_tra_frs     ! routine called in tranxt.F90  
     26   PUBLIC bdy_tra      ! routine called in tranxt.F90  
    2527 
    2628   !!---------------------------------------------------------------------- 
    … …  
    3133CONTAINS 
    3234 
    33    SUBROUTINE bdy_tra_frs( kt ) 
     35   SUBROUTINE bdy_tra( kt ) 
     36      !!---------------------------------------------------------------------- 
     37      !!                  ***  SUBROUTINE bdy_dyn3d  *** 
     38      !! 
     39      !! ** Purpose : - Apply open boundary conditions for baroclinic velocities 
     40      !! 
     41      !!---------------------------------------------------------------------- 
     42      INTEGER, INTENT( in ) :: kt     ! Main time step counter 
     43      !! 
     44      INTEGER               :: ib_bdy ! Loop index 
     45 
     46      DO ib_bdy=1, nb_bdy 
     47 
     48         SELECT CASE( nn_tra(ib_bdy) ) 
     49         CASE(jp_none) 
     50            CYCLE 
     51         CASE(jp_frs) 
     52            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     53         CASE DEFAULT 
     54            CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' ) 
     55         END SELECT 
     56      ENDDO 
     57 
     58   END SUBROUTINE bdy_tra 
     59 
     60   SUBROUTINE bdy_tra_frs( idx, dta, kt ) 
    3461      !!---------------------------------------------------------------------- 
    3562      !!                 ***  SUBROUTINE bdy_tra_frs  *** 
    3663      !!                     
    37       !! ** Purpose : Apply the Flow Relaxation Scheme for tracers in the   
    38       !!              case of unstructured open boundaries. 
     64      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries. 
    3965      !!  
    4066      !! Reference : Engedahl H., 1995, Tellus, 365-382. 
    4167      !!---------------------------------------------------------------------- 
    42       INTEGER, INTENT( in ) ::   kt 
     68      INTEGER,         INTENT(in) ::   kt 
     69      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     70      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    4371      !!  
    4472      REAL(wp) ::   zwgt           ! boundary weight 
    … …  
    4775      !!---------------------------------------------------------------------- 
    4876      ! 
    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                tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)          
    64                tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik) 
    65             END DO 
    66          END DO  
    67          !                                              ! Boundary points should be updated 
    68          CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )      
    69          CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )     
    70          ! 
    71       ENDIF ! ln_tra_frs 
    7277      ! 
     78      igrd = 1                       ! Everything is at T-points here 
     79      DO ib = 1, idx%nblen(igrd) 
     80         DO ik = 1, jpkm1 
     81            ii = idx%nbi(ib,igrd) 
     82            ij = idx%nbj(ib,igrd) 
     83            zwgt = idx%nbw(ib,igrd) 
     84            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)          
     85            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) 
     86         END DO 
     87      END DO  
     88      ! 
     89      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
     90      ! 
     91      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     92   ! 
    7393   END SUBROUTINE bdy_tra_frs 
    7494    
    … …  
    7898   !!---------------------------------------------------------------------- 
    7999CONTAINS 
    80    SUBROUTINE bdy_tra_frs(kt)      ! Empty routine 
    81       WRITE(*,*) 'bdy_tra_frs: You should not have seen this print! error?', kt 
    82    END SUBROUTINE bdy_tra_frs 
     100   SUBROUTINE bdy_tra(kt)      ! Empty routine 
     101      WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt 
     102   END SUBROUTINE bdy_tra 
    83103#endif 
    84104 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r2528 r3116  
    7171      !! 
    7272      INTEGER  ::   ji, jj, jk, jb, jgrd 
    73       INTEGER  ::   ii, ij 
     73      INTEGER  ::   ib_bdy, ii, ij 
    7474      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst 
     75      TYPE(OBC_INDEX), POINTER :: idx 
    7576      !!----------------------------------------------------------------------------- 
    7677 
    … …  
    9192      ! ------------------------------------------------ 
    9293      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) 
     94      DO ib_bdy = 1, nb_bdy 
     95         idx => idx_bdy(ib_bdy) 
     96 
     97         jgrd = 2                               ! cumulate u component contribution first  
     98         DO jb = 1, idx%nblenrim(jgrd) 
     99            DO jk = 1, jpkm1 
     100               ii = idx%nbi(jb,jgrd) 
     101               ij = idx%nbj(jb,jgrd) 
     102               zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     103            END DO 
    99104         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)  
     105         jgrd = 3                               ! then add v component contribution 
     106         DO jb = 1, idx%nblenrim(jgrd) 
     107            DO jk = 1, jpkm1 
     108               ii = idx%nbi(jb,jgrd) 
     109               ij = idx%nbj(jb,jgrd) 
     110               zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
     111            END DO 
    107112         END DO 
     113 
    108114      END DO 
    109115      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
    … …  
    118124      ! ------------------------------------------------------------- 
    119125      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) 
     126      DO ib_bdy = 1, nb_bdy 
     127         idx => idx_bdy(ib_bdy) 
     128 
     129         jgrd = 2                               ! correct u component 
     130         DO jb = 1, idx%nblenrim(jgrd) 
     131            DO jk = 1, jpkm1 
     132               ii = idx%nbi(jb,jgrd) 
     133               ij = idx%nbj(jb,jgrd) 
     134               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk) 
     135               ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     136            END DO 
    127137         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) 
     138         jgrd = 3                              ! correct v component 
     139         DO jb = 1, idx%nblenrim(jgrd) 
     140            DO jk = 1, jpkm1 
     141               ii = idx%nbi(jb,jgrd) 
     142               ij = idx%nbj(jb,jgrd) 
     143               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
     144               ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
     145            END DO 
    136146         END DO 
     147 
    137148      END DO 
    138149      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
Note: See TracChangeset for help on using the changeset viewer.