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 1151 – NEMO

Changeset 1151


Ignore:
Timestamp:
2008-06-26T14:54:53+02:00 (16 years ago)
Author:
rblod
Message:

Add new OBC, see ticket #224

Location:
trunk
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/CONFIG/GYRE/EXP00/namelist

    r1149 r1151  
    250250&namobc        !   open boundaries parameters                           ("key_obc") 
    251251!----------------------------------------------------------------------- 
    252     nobc_dta   =    0      !  = 0 the obc data are equal to the initial state 
     252    nobc_dta   =    1      !  = 0 the obc data are equal to the initial state 
    253253                           !  = 1 the obc data are read in 'obc.dta' files 
     254    cffile     = 'annual'  !  set to annual if obc datafile hold 1 year of data 
     255                           !  set to monthly if obc datafile hold 1 month of data 
    254256    rdpein     =    1.     !  ??? 
    255257    rdpwin     =    1.     !  ??? 
    256     rdpnin     =   30.     !  ??? 
     258    rdpnin     =    1.     !  ??? 
    257259    rdpsin     =    1.     !  ??? 
    258     rdpeob     = 1500.     !  time relaxation (days) for the east  open boundary 
     260    rdpeob     = 3000.     !  time relaxation (days) for the east  open boundary 
    259261    rdpwob     =   15.     !    "        "             "     west         " 
    260     rdpnob     =  150.     !    "        "             "     north        " 
     262    rdpnob     = 3000.     !    "        "             "     north        " 
    261263    rdpsob     =   15.     !    "        "             "     south        " 
    262264    zbsic1     =  140.e+6  !   barotropic stream function on first  isolated coastline 
    263265    zbsic2     =    1.e+6  !    "                   "        second       " 
    264266    zbsic3     =    0.     !    "                   "        thrid        " 
    265     ln_obc_clim= .true.    !  climatological obc data files (T) or not (F) 
    266     ln_vol_cst = .false.   !  impose the total volume conservation (T) or not (F) 
     267    ln_obc_clim= .false.   !  climatological obc data files (T) or not (F) 
     268    ln_vol_cst = .true.    !  impose the total volume conservation (T) or not (F) 
    267269/ 
    268270!----------------------------------------------------------------------- 
  • trunk/CONFIG/GYRE_LOBSTER/EXP00/namelist

    r1149 r1151  
    250250&namobc        !   open boundaries parameters                           ("key_obc") 
    251251!----------------------------------------------------------------------- 
    252     nobc_dta   =    0      !  = 0 the obc data are equal to the initial state 
     252    nobc_dta   =    1      !  = 0 the obc data are equal to the initial state 
    253253                           !  = 1 the obc data are read in 'obc.dta' files 
     254    cffile     = 'annual'  !  set to annual if obc datafile hold 1 year of data 
     255                           !  set to monthly if obc datafile hold 1 month of data 
    254256    rdpein     =    1.     !  ??? 
    255257    rdpwin     =    1.     !  ??? 
    256     rdpnin     =   30.     !  ??? 
     258    rdpnin     =    1.     !  ??? 
    257259    rdpsin     =    1.     !  ??? 
    258     rdpeob     = 1500.     !  time relaxation (days) for the east  open boundary 
     260    rdpeob     = 3000.     !  time relaxation (days) for the east  open boundary 
    259261    rdpwob     =   15.     !    "        "             "     west         " 
    260     rdpnob     =  150.     !    "        "             "     north        " 
     262    rdpnob     = 3000.     !    "        "             "     north        " 
    261263    rdpsob     =   15.     !    "        "             "     south        " 
    262264    zbsic1     =  140.e+6  !   barotropic stream function on first  isolated coastline 
    263265    zbsic2     =    1.e+6  !    "                   "        second       " 
    264266    zbsic3     =    0.     !    "                   "        thrid        " 
    265     ln_obc_clim= .true.    !  climatological obc data files (T) or not (F) 
    266     ln_vol_cst = .false.   !  impose the total volume conservation (T) or not (F) 
     267    ln_obc_clim= .false.   !  climatological obc data files (T) or not (F) 
     268    ln_vol_cst = .true.    !  impose the total volume conservation (T) or not (F) 
    267269/ 
    268270!----------------------------------------------------------------------- 
  • trunk/CONFIG/ORCA2_LIM/EXP00/namelist

    r1149 r1151  
    250250&namobc        !   open boundaries parameters                           ("key_obc") 
    251251!----------------------------------------------------------------------- 
    252     nobc_dta   =    0      !  = 0 the obc data are equal to the initial state 
     252    nobc_dta   =    1      !  = 0 the obc data are equal to the initial state 
    253253                           !  = 1 the obc data are read in 'obc.dta' files 
     254    cffile     = 'annual'  !  set to annual if obc datafile hold 1 year of data 
     255                           !  set to monthly if obc datafile hold 1 month of data 
    254256    rdpein     =    1.     !  ??? 
    255257    rdpwin     =    1.     !  ??? 
    256     rdpnin     =   30.     !  ??? 
     258    rdpnin     =    1.     !  ??? 
    257259    rdpsin     =    1.     !  ??? 
    258     rdpeob     = 1500.     !  time relaxation (days) for the east  open boundary 
     260    rdpeob     = 3000.     !  time relaxation (days) for the east  open boundary 
    259261    rdpwob     =   15.     !    "        "             "     west         " 
    260     rdpnob     =  150.     !    "        "             "     north        " 
     262    rdpnob     = 3000.     !    "        "             "     north        " 
    261263    rdpsob     =   15.     !    "        "             "     south        " 
    262264    zbsic1     =  140.e+6  !   barotropic stream function on first  isolated coastline 
    263265    zbsic2     =    1.e+6  !    "                   "        second       " 
    264266    zbsic3     =    0.     !    "                   "        thrid        " 
    265     ln_obc_clim= .true.    !  climatological obc data files (T) or not (F) 
    266     ln_vol_cst = .false.   !  impose the total volume conservation (T) or not (F) 
     267    ln_obc_clim= .false.   !  climatological obc data files (T) or not (F) 
     268    ln_vol_cst = .true.    !  impose the total volume conservation (T) or not (F) 
    267269/ 
    268270!----------------------------------------------------------------------- 
  • trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist

    r1149 r1151  
    112112                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only) 
    113113   ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr) 
    114    ln_rnf      = .false.   !  runoffs (T => fill namsbc_rnf) 
    115    ln_ssr      = .false.   !  Sea Surface Restoring on T and/or S (T => fill namsbc_ssr) 
     114   ln_rnf      = .true.    !  runoffs (T => fill namsbc_rnf) 
     115   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S (T => fill namsbc_ssr) 
    116116   nn_fwb      = 0         !  FreshWater Budget: =0 unchecked                              ,  
    117117                           !                     =1 annual global mean of e-p-r set to zero   , 
     
    146146!              !      file name     ! frequency (hours) !  variable  ! time interpol. !  clim   ! 'yearly' or ! 
    147147!              !                    !  (if <0  months)  !    name    !    (logical)   !  (T/F)  !  'monthly'  ! 
    148    sn_utau     =    'taux_1m'       ,       -1.         , 'sozotaux' ,    .false.     , .true.  ,   'yearly'     
    149    sn_vtau     =    'tauy_1m'       ,       -1.         , 'sometauy' ,    .false.     , .true.  ,   'yearly'     
    150    sn_wndm     =    'flx'           ,       -1.         , 'socliowi' ,    .false.     , .true.  ,   'yearly'       
    151    sn_tair     =    'flx'           ,       -1.         , 'socliot2' ,    .false.     , .true.  ,   'yearly'      
    152    sn_humi     =    'flx'           ,       -1.         , 'socliohu' ,    .false.     , .true.  ,   'yearly'     
     148   sn_utau     =    'taux_1m'       ,       -1.         , 'sozotaux' ,    .true.      , .true.  ,   'yearly'     
     149   sn_vtau     =    'tauy_1m'       ,       -1.         , 'sometauy' ,    .true.      , .true.  ,   'yearly'     
     150   sn_wndm     =    'flx'           ,       -1.         , 'socliowi' ,    .true.      , .true.  ,   'yearly'       
     151   sn_tair     =    'flx'           ,       -1.         , 'socliot2' ,    .true.      , .true.  ,   'yearly'      
     152   sn_humi     =    'flx'           ,       -1.         , 'socliohu' ,    .true.      , .true.  ,   'yearly'     
    153153   sn_ccov     =    'flx'           ,       -1.         , 'socliocl' ,    .false.     , .true.  ,   'yearly'       
    154154   sn_prec     =    'flx'           ,       -1.         , 'socliopl' ,    .false.     , .true.  ,   'yearly'        
     
    206206!              !                    !  (if <0  months)  !   name     !    (logical)   !  (T/F)  !  'monthly'  ! 
    207207   sn_sst      =    'sst_data'     ,        24.        ,  'sst'     ,     .false.    , .false. ,   'yearly' 
    208    sn_sss      =    'sss_data'     ,        -1.        ,  'sss'     ,     .true.     , .false. ,   'yearly' 
     208   sn_sss      =    'sss_data'     ,        -1.        ,  'sss'     ,     .true.     , .true. ,   'yearly' 
    209209!    
    210210   cn_dir      = './'      !  root directory for the location of the runoff files 
    211211   nn_sstr     =     0     !  add a retroaction term in the surface heat       flux (=1) or not (=0) 
    212    nn_sssr     =     0     !  add a damping     term in the surface freshwater flux (=1) or not (=0) 
     212   nn_sssr     =     1     !  add a damping     term in the surface freshwater flux (=1) or not (=0) 
    213213   dqdt        =   -40.    !  magnitude of the retroaction on temperature   [W/m2/K] 
    214214   deds        =   -27.7   !  magnitude of the damping on salinity   [mm/day/psu] 
     
    250250&namobc        !   open boundaries parameters                           ("key_obc") 
    251251!----------------------------------------------------------------------- 
    252     nobc_dta   =    0      !  = 0 the obc data are equal to the initial state 
     252    nobc_dta   =    1      !  = 0 the obc data are equal to the initial state 
    253253                           !  = 1 the obc data are read in 'obc.dta' files 
     254    cffile     = 'annual'  !  set to annual if obc datafile hold 1 year of data 
     255                           !  set to monthly if obc datafile hold 1 month of data 
    254256    rdpein     =    1.     !  ??? 
    255257    rdpwin     =    1.     !  ??? 
    256     rdpnin     =   30.     !  ??? 
     258    rdpnin     =    1.     !  ??? 
    257259    rdpsin     =    1.     !  ??? 
    258     rdpeob     = 1500.     !  time relaxation (days) for the east  open boundary 
     260    rdpeob     = 3000.     !  time relaxation (days) for the east  open boundary 
    259261    rdpwob     =   15.     !    "        "             "     west         " 
    260     rdpnob     =  150.     !    "        "             "     north        " 
     262    rdpnob     = 3000.     !    "        "             "     north        " 
    261263    rdpsob     =   15.     !    "        "             "     south        " 
    262264    zbsic1     =  140.e+6  !   barotropic stream function on first  isolated coastline 
    263265    zbsic2     =    1.e+6  !    "                   "        second       " 
    264266    zbsic3     =    0.     !    "                   "        thrid        " 
    265     ln_obc_clim= .true.    !  climatological obc data files (T) or not (F) 
    266     ln_vol_cst = .false.   !  impose the total volume conservation (T) or not (F) 
     267    ln_obc_clim= .false.   !  climatological obc data files (T) or not (F) 
     268    ln_vol_cst = .true.    !  impose the total volume conservation (T) or not (F) 
    267269/ 
    268270!----------------------------------------------------------------------- 
  • trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r1146 r1151  
    434434      ENDIF 
    435435 
     436#if defined key_agrif 
     437      IF (lk_obc .AND. Agrif_Root() ) THEN 
     438#else 
     439      IF (lk_obc) THEN 
     440#endif 
     441         zssha(:,:)=zssha(:,:)*obctmsk(:,:) 
     442         CALL lbc_lnk(zssha,'T',1.)  ! absolutly compulsory !! (jmm) 
     443      ENDIF 
     444 
    436445      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter 
    437446         ! swap of arrays 
  • trunk/NEMO/OPA_SRC/DYN/wzvmod.F90

    r1146 r1151  
    104104         END DO 
    105105 
    106 #if defined key_obc && ( key_dynspg_exp || key_dynspg_ts ) 
     106#if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 
    107107         ! open boundaries (div must be zero behind the open boundary) 
    108108         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
  • trunk/NEMO/OPA_SRC/OBC/obc_oce.F90

    r1146 r1151  
    55   !!============================================================================== 
    66   !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    7    !! $Header$  
     7   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/OBC/obc_oce.F90,v 1.6 2006/05/10 17:38:46 opalod Exp $  
    88   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    99   !!---------------------------------------------------------------------- 
     
    1616   !!  8.5   06/02   (C. Talandier)  modules 
    1717   !!        06/04   (F. Durand) ORCA2_ZIND config 
    18    !!        06/04   (F. Durand) Dimensions of arrays vsdta, tsdta, ssdta, 
    19    !!                vndta, tndta, sndta, uwdta, twdta, swdta, uedta, tedta, sedta 
    20    !!                are defined to the actual dimensions of the OBs i.e.  
    21    !!         (jpisd:jpisf,jpk,jptobc) for the South OB 
    22    !!            (jpind:jpinf,jpk,jptobc) for the North OB  
    23    !!         (jpjwd:jpjwf,jpk,jptobc) for the West OB 
    24    !!            (jpjed:jpjef,jpk,jptobc) for the East OB 
    2518   !!                  
    2619   !!---------------------------------------------------------------------- 
     
    137130                          !: (if Flather's algoritm applied at open boundary) 
    138131 
    139    REAL(wp), DIMENSION(jpjed:jpjef,0:jptobc+1) ::   &  !: 
    140       sshedta, ubtedta    !: array used for interpolating monthly data on the east boundary 
    141  
    142    REAL(wp), DIMENSION(jpjed:jpjef,jpk,jptobc) ::   &  !: 
    143       uedta, tedta, sedta !: array used for interpolating monthly data on the east boundary 
    144  
    145132   !!------------------------------- 
    146133   !! Arrays for radiative East OBC:  
     
    196183      sshfow_b            !: west boundary ssh correction averaged over the barotropic loop 
    197184                          !: (if Flather's algoritm applied at open boundary) 
    198  
    199    REAL(wp), DIMENSION(jpjwd:jpjwf,0:jptobc+1) ::   &  !: 
    200       sshwdta, ubtwdta    !: array used for interpolating monthly data on the west boundary 
    201  
    202    REAL(wp), DIMENSION(jpjwd:jpjwf,jpk,jptobc) ::   &  !: 
    203       uwdta, twdta, swdta !: array used for interpolating monthly data on the west boundary 
    204185 
    205186   !!------------------------------- 
     
    257238      sshfon_b            !: north boundary ssh correction averaged over the barotropic loop 
    258239                          !: (if Flather's algoritm applied at open boundary) 
    259  
    260    REAL(wp), DIMENSION(jpind:jpinf,0:jptobc+1) ::   &  !: 
    261       sshndta, vbtndta   !: array used for interpolating monthly data on the north boundary 
    262  
    263    REAL(wp), DIMENSION(jpind:jpinf,jpk,jptobc) ::   &  !: 
    264       vndta, tndta, sndta !: array used for interpolating monthly data on the north boundary 
    265240 
    266241   !!-------------------------------- 
     
    318293                          !: (if Flather's algoritm applied at open boundary) 
    319294 
    320    REAL(wp), DIMENSION(jpisd:jpisf,0:jptobc+1) ::    &    !: 
    321       sshsdta, vbtsdta    !: array used for interpolating monthly data on the south boundary 
    322  
    323    REAL(wp), DIMENSION(jpisd:jpisf,jpk,jptobc) ::    &    !: 
    324       vsdta, tsdta, ssdta   !: array used for interpolating monthly data on the south boundary 
    325  
    326295   !!-------------------------------- 
    327296   !! Arrays for radiative South OBC 
     
    344313      usmsk, vsmsk, tsmsk           !: 2D mask for the South OB 
    345314 
    346    ! Note that those arrays are optimized for mpp case  
    347    ! (hence the dimension jpj is the size of one processor subdomain) 
     315   CHARACTER ( len=20 ) :: cffile 
    348316 
    349317#else 
  • trunk/NEMO/OPA_SRC/OBC/obcdta.F90

    r1146 r1151  
    11MODULE obcdta 
    2    !!============================================================================== 
    3    !!                            ***  MODULE obcdta  *** 
    4    !! Open boundary data : read the data for the open boundaries. 
    5    !!============================================================================== 
    6    !! History : 
    7    !!        !  98-05 (J.M. Molines) Original code 
    8    !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90 
    9    !!   9.0  !  04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input 
    10    !!    "   ! 07-07  (D. Storkey) Change to arguments for iom_gettime 
    11    !!------------------------------------------------------------------------------ 
     2  !!============================================================================== 
     3  !!                            ***  MODULE obcdta  *** 
     4  !! Open boundary data : read the data for the open boundaries. 
     5  !!============================================================================== 
    126#if defined key_obc 
    13    !!------------------------------------------------------------------------------ 
    14    !!   'key_obc'         :                                Open Boundary Conditions 
    15    !!------------------------------------------------------------------------------ 
    16    !!   obc_dta           : read u, v, t, s data along each open boundary 
    17    !!   obc_dta_psi       : read psi data along each open boundary (rigid lid only) 
    18    !!------------------------------------------------------------------------------ 
    19    !! * Modules used 
    20    USE oce             ! ocean dynamics and tracers  
    21    USE dom_oce         ! ocean space and time domain 
    22    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    23    USE phycst          ! physical constants 
    24    USE obc_oce         ! ocean open boundary conditions 
    25    USE daymod          ! calendar 
    26    USE in_out_manager  ! I/O logical units 
    27    USE lib_mpp         ! distributed memory computing 
    28    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient 
    29    USE iom 
    30 #  if defined key_dynspg_rl 
    31    USE obccli          ! climatological obc, use only in rigid-lid case 
    32 #  endif 
    33  
    34    IMPLICIT NONE 
    35    PRIVATE 
    36  
    37    !! * Accessibility 
    38    PUBLIC obc_dta        ! routines called by step.F90 
    39    PUBLIC obc_dta_bt     ! routines called by dynspg_ts.F90 
    40  
    41    !! * Shared module variables 
    42    INTEGER ::   & 
    43       nlecto,   &  ! switch for the first read 
    44       ntobc1,   &  ! first record used 
    45       ntobc2,   &  ! second record used 
    46       ntobc        ! number of time steps in OBC files  
    47  
    48    REAL(wp), DIMENSION(:), ALLOCATABLE :: tcobc      ! time_counter variable of BCs 
    49  
    50    !! * Substitutions 
    51 #  include "domzgr_substitute.h90"                                             
     7  !!------------------------------------------------------------------------------ 
     8  !!   'key_obc'         :                                Open Boundary Conditions 
     9  !!------------------------------------------------------------------------------ 
     10  !!   obc_dta           : read u, v, t, s data along each open boundary 
     11  !!   obc_dta_psi       : read psi data along each open boundary (rigid lid only) 
     12  !!------------------------------------------------------------------------------ 
     13  !! * Modules used 
     14  USE oce             ! ocean dynamics and tracers  
     15  USE dom_oce         ! ocean space and time domain 
     16  USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     17  USE phycst          ! physical constants 
     18  USE obc_oce         ! ocean open boundary conditions 
     19  USE daymod          ! calendar 
     20  USE in_out_manager  ! I/O logical units 
     21  USE lib_mpp         ! distributed memory computing 
     22  USE dynspg_oce 
     23  USE ioipsl          ! now only for  ymds2ju function  
     24  USE iom             !  
     25  USE obccli          ! used when dynspg_rl 
     26 
     27  IMPLICIT NONE 
     28  PRIVATE 
     29 
     30  !! * Accessibility 
     31  PUBLIC obc_dta      ! routines called by step.F90 
     32  PUBLIC obc_dta_bt   ! routines called by dynspg_ts.F90 
     33 
     34  !! * Shared module variables 
     35  REAL(wp),  DIMENSION(2)              ::  zjcnes_obc   !  
     36  REAL(wp),  DIMENSION(:), ALLOCATABLE :: ztcobc 
     37  REAL(wp) :: rdt_obc 
     38  REAL(wp) :: zjcnes 
     39  INTEGER :: imm0, iyy0, idd0, iyy, imm, idd 
     40  INTEGER :: nt_a=2, nt_b=1, itobc, ndate0_cnes, nday_year0 
     41  INTEGER ::  itobce, itobcw, itobcs, itobcn, itobc_b  ! number of time steps in OBC files 
     42  REAL(wp) :: zraaj=365. 
     43 
     44  INTEGER ::   & 
     45       ntobc      , &     !:  where we are in the obc file 
     46       ntobc_b    , &     !:  first record used 
     47       ntobc_a            !:  second record used 
     48 
     49  CHARACTER (len=40) :: &    ! name of data files 
     50       cl_obc_eTS   , cl_obc_eU,  & 
     51       cl_obc_wTS   , cl_obc_wU,  & 
     52       cl_obc_nTS   , cl_obc_nV,  & 
     53       cl_obc_sTS   , cl_obc_sV 
     54 
     55  ! arrays used for interpolating time dependent data on the boundaries 
     56  REAL(wp), DIMENSION(jpj,jpk,jptobc) :: uedta, vedta, tedta, sedta    ! East 
     57  REAL(wp), DIMENSION(jpj,jpk,jptobc) :: uwdta, vwdta, twdta, swdta    ! West 
     58  REAL(wp), DIMENSION(jpi,jpk,jptobc) :: undta, vndta, tndta, sndta    ! North 
     59  REAL(wp), DIMENSION(jpi,jpk,jptobc) :: usdta, vsdta, tsdta, ssdta    ! South 
     60 
     61  LOGICAL, DIMENSION (jpj,jpk ) :: ltemsk=.TRUE., luemsk=.TRUE., lvemsk=.TRUE.  ! boolean msks 
     62  LOGICAL, DIMENSION (jpj,jpk ) :: ltwmsk=.TRUE., luwmsk=.TRUE., lvwmsk=.TRUE.  ! used for outliers 
     63  LOGICAL, DIMENSION (jpi,jpk ) :: ltnmsk=.TRUE., lunmsk=.TRUE., lvnmsk=.TRUE.  ! checks 
     64  LOGICAL, DIMENSION (jpi,jpk ) :: ltsmsk=.TRUE., lusmsk=.TRUE., lvsmsk=.TRUE. 
     65 
     66  !! * Substitutions 
    5267#  include "obc_vectopt_loop_substitute.h90" 
    53    !!--------------------------------------------------------------------------------- 
    54    !!   OPA 9.0 , LODYC-IPSL  (2003) 
    55    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/OBC/obcdta.F90,v 1.14 2007/06/29 17:01:51 opalod Exp $ 
    56    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    57    !!--------------------------------------------------------------------------------- 
     68  !!--------------------------------------------------------------------------------- 
     69  !!   OPA 9.0 , LODYC-IPSL  (2003) 
     70  !!--------------------------------------------------------------------------------- 
    5871 
    5972CONTAINS 
    6073 
    61    SUBROUTINE obc_dta( kt ) 
    62       !!-------------------------------------------------------------------- 
    63       !!              ***  SUBROUTINE obc_dta  *** 
    64       !!                    
    65       !! ** Purpose :   Find the climatological boundary arrays for the specified date, 
    66       !!   The boundary arrays are netcdf files. Three possible cases: 
    67       !!   - one time frame only in the file (time dimension = 1). 
    68       !!     in that case the boundary data does not change in time. 
    69       !!   - many time frames. In that case,  if we have 12 frames 
    70       !!     we assume monthly fields. 
    71       !!     Else, we assume that time_counter is in seconds 
    72       !!     since the beginning of either the current year or a reference 
    73       !!     year given in the namelist. 
    74       !!     (no check is done so far but one would have to check the "unit" 
    75       !!     attribute of variable time_counter). 
    76       !! 
    77       !! History : 
    78       !!        !  98-05 (J.M. Molines) Original code 
    79       !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90 
    80       !!   9.0  !  04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input 
    81       !!-------------------------------------------------------------------- 
    82       !! * Arguments 
    83       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    84        
    85       !! * Local declarations 
    86       INTEGER ::   ji, jj, jk, ii, ij   ! dummy loop indices 
    87       INTEGER ::   itimo, iman, imois 
    88       INTEGER ::   i15 
    89       REAL(wp) ::   zxy 
    90       !! * Ajouts FD 
    91       INTEGER ::  isrel              ! number of seconds since 1/1/1992 
    92       INTEGER, DIMENSION(1) ::  itobce, itobcw,  & ! number of time steps in OBC files 
    93                                 itobcs, itobcn     !    "       "       "       " 
    94       INTEGER ::  istop         
    95       INTEGER ::  iprint                              ! frequency for printouts. 
    96       INTEGER ::  idvar, id_e, id_w, id_n, id_s       ! file identifiers 
    97       LOGICAL :: llnot_done 
    98       CHARACTER(LEN=25) :: cl_vname 
    99       !!-------------------------------------------------------------------- 
    100  
    101       IF( lk_dynspg_rl )  THEN 
    102          CALL obc_dta_psi (kt)     ! update bsf data at open boundaries 
    103          IF ( nobc_dta == 1 .AND. kt == nit000 ) CALL ctl_stop( 'obcdta : time-variable psi boundary data not allowed yet' ) 
    104       ENDIF 
    105             
    106       ! 1.   First call: check time frames available in files. 
    107       ! ------------------------------------------------------- 
    108  
    109       IF ( kt == nit000 ) THEN 
    110        
    111          nlecto = 0 
    112          itobce(1) = 0   ;    itobcw(1) = 0 
    113          itobcn(1) = 0   ;    itobcs(1) = 0 
    114  
    115          IF (lwp) WRITE(numout,*) 
    116          IF (lwp) WRITE(numout,*)     'obc_dta : find boundary data' 
    117          IF (lwp) WRITE(numout,*)     '~~~~~~~' 
    118               
    119          IF ( nobc_dta == 0 ) THEN 
    120             IF(lwp) WRITE(numout,*)  '  OBC data taken from initial conditions.' 
    121             ntobc1 = 1 
    122             ntobc2 = 1 
    123          ELSE     
    124             IF (lwp) WRITE(numout,*)  '  OBC data taken from netcdf files.' 
    125             IF (lwp) WRITE(numout,*)  '  climatology (T/F):',ln_obc_clim 
    126             ! check the number of time steps in the files. 
    127             cl_vname = 'time_counter' 
    128             IF ( lp_obc_east ) THEN 
    129                CALL iom_open ( 'obceast_TS.nc' , id_e ) 
    130                idvar = iom_varid( id_e, cl_vname, kdimsz = itobce ) 
    131             ENDIF 
    132             IF ( lp_obc_west ) THEN 
    133                CALL iom_open ( 'obcwest_TS.nc' , id_w ) 
    134                idvar = iom_varid( id_w, cl_vname, kdimsz = itobcw ) 
    135             ENDIF 
    136             IF ( lp_obc_north ) THEN 
    137                CALL iom_open ( 'obcnorth_TS.nc', id_n ) 
    138                idvar = iom_varid( id_n, cl_vname, kdimsz = itobcn ) 
    139             ENDIF 
    140             IF ( lp_obc_south ) THEN 
    141                CALL iom_open ( 'obcsouth_TS.nc', id_s ) 
    142                idvar = iom_varid( id_s, cl_vname, kdimsz = itobcs ) 
    143             ENDIF 
    144  
    145             ntobc = MAX(itobce(1),itobcw(1),itobcn(1),itobcs(1)) 
    146             istop = 0 
    147             IF ( lp_obc_east  .AND. itobce(1) /= ntobc ) istop = 1  
    148             IF ( lp_obc_west  .AND. itobcw(1) /= ntobc ) istop = 1       
    149             IF ( lp_obc_north .AND. itobcn(1) /= ntobc ) istop = 1 
    150             IF ( lp_obc_south .AND. itobcs(1) /= ntobc ) istop = 1  
    151             IF ( istop /= 0 )  THEN 
    152                WRITE(ctmp1,*) ' east, west, north, south: ', itobce(1), itobcw(1), itobcn(1), itobcs(1) 
    153                CALL ctl_stop( 'obcdta : all files must have the same number of time steps', ctmp1 ) 
    154             ENDIF 
    155             IF ( ntobc == 1 ) THEN 
    156                IF ( lwp ) WRITE(numout,*) ' obcdta found one time step only in the OBC files' 
    157             ELSE 
    158                ALLOCATE (tcobc(ntobc)) 
    159                llnot_done = .TRUE. 
    160                IF ( lp_obc_east ) THEN 
    161                   IF ( llnot_done ) THEN 
    162                      CALL iom_gettime ( id_e, tcobc, TRIM(cl_vname) ) 
    163                      llnot_done = .FALSE. 
    164                   ENDIF 
    165                   CALL iom_close (id_e) 
    166                ENDIF 
    167                IF ( lp_obc_west ) THEN 
    168                   IF ( llnot_done ) THEN 
    169                      CALL iom_gettime ( id_w, tcobc, TRIM(cl_vname) ) 
    170                      llnot_done = .FALSE. 
    171                  ENDIF 
    172                  CALL iom_close (id_w) 
    173                ENDIF 
    174                IF ( lp_obc_north ) THEN 
    175                   IF ( llnot_done ) THEN 
    176                      CALL iom_gettime ( id_n, tcobc, TRIM(cl_vname) ) 
    177                      llnot_done = .FALSE. 
    178                   ENDIF 
    179                   CALL iom_close (id_n) 
    180                ENDIF 
    181                IF ( lp_obc_south ) THEN 
    182                   IF ( llnot_done ) THEN 
    183                      CALL iom_gettime ( id_s, tcobc, TRIM(cl_vname) ) 
    184                      llnot_done = .FALSE. 
    185                   ENDIF 
    186                   CALL iom_close (id_s) 
    187                ENDIF 
    188                IF ( lwp ) WRITE(numout,*) ' obcdta found', ntobc,' time steps in the OBC files' 
    189                IF ( .NOT. ln_obc_clim .AND. ntobc == 12 ) THEN 
    190                   IF ( lwp ) WRITE(numout,*) '  WARNING: With monthly data we assume climatology' 
    191                   ln_obc_clim = .true. 
    192                ENDIF 
    193             ENDIF 
    194          ENDIF 
    195        
    196        ! 1.1  Tangential velocities set to zero 
    197        ! -------------------------------------- 
    198          IF( lp_obc_east  ) vfoe = 0.e0 
    199          IF( lp_obc_west  ) vfow = 0.e0 
    200          IF( lp_obc_south ) ufos = 0.e0 
    201          IF( lp_obc_north ) ufon = 0.e0 
    202        
    203        ! 1.2  Data temperature, salinity, normal velocities set to zero 
    204        !                        or initial conditions if nobc_dta == 0 
    205        ! -------------------------------------------------------------- 
    206  
    207          IF( lp_obc_east )   THEN 
    208             ! initialisation to zero 
    209             sedta(:,:,:) = 0.e0 
    210             tedta(:,:,:) = 0.e0 
    211             uedta(:,:,:) = 0.e0 
    212             !                                    ! ================== ! 
    213             IF( nobc_dta == 0 )   THEN           ! initial state used 
    214                !                                 ! ================== ! 
    215                !  Fills sedta, tedta, uedta (global arrays) 
    216                !  Remark: this works for njzoom = 1. 
    217                !          Should the definition of ij include njzoom? 
    218                DO ji = nie0, nie1 
    219                   DO jk = 1, jpkm1 
    220                      DO jj = nje0p1, nje1m1 
    221                         ij = jj -1 + njmpp 
    222                         sedta(ij,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk) 
    223                         tedta(ij,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk) 
    224                         uedta(ij,jk,1) = un(ji,jj,jk)*umask(ji,jj,jk) 
    225                      END DO 
    226                   END DO 
    227                END DO 
    228             ENDIF 
    229          ENDIF 
    230  
    231          IF( lp_obc_west )   THEN 
    232             ! initialisation to zero 
    233             swdta(:,:,:) = 0.e0 
    234             twdta(:,:,:) = 0.e0 
    235             uwdta(:,:,:) = 0.e0 
    236             !                                    ! ================== ! 
    237             IF( nobc_dta == 0 )   THEN           ! initial state used ! 
    238                !                                 ! ================== ! 
    239                !  Fills swdta, twdta, uwdta (global arrays) 
    240                !  Remark: this works for njzoom = 1.  
    241                !          Should the definition of ij include njzoom? 
    242                DO ji = niw0, niw1 
    243                   DO jk = 1, jpkm1 
    244                      DO jj = njw0p1, njw1m1 
    245                         ij = jj -1 + njmpp 
    246                         swdta(ij,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk) 
    247                         twdta(ij,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk) 
    248                         uwdta(ij,jk,1) = un(ji,jj,jk)*umask(ji,jj,jk) 
    249                      END DO 
    250                   END DO 
    251                END DO 
    252             ENDIF 
    253          ENDIF 
    254  
    255          IF( lp_obc_north)   THEN 
    256             ! initialisation to zero 
    257             sndta(:,:,:) = 0.e0 
    258             tndta(:,:,:) = 0.e0 
    259             vndta(:,:,:) = 0.e0 
    260             !                                    ! ================== ! 
    261             IF( nobc_dta == 0 )   THEN           ! initial state used 
    262                !                                 ! ================== ! 
    263                !  Fills sndta, tndta, vndta (global arrays) 
    264                !  Remark: this works for njzoom = 1.  
    265                !          Should the definition of ij include njzoom? 
    266                DO jj = njn0, njn1 
    267                   DO jk = 1, jpkm1 
    268                      DO ji = nin0p1, nin1m1 
    269                         ii = ji -1 + nimpp 
    270                         sndta(ii,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk) 
    271                         tndta(ii,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk) 
    272                         vndta(ii,jk,1) = vn(ji,jj,jk)*vmask(ji,jj,jk) 
    273                      END DO 
    274                   END DO 
    275                END DO 
    276             ENDIF 
    277          ENDIF 
    278  
    279          IF( lp_obc_south )   THEN 
    280             ! initialisation to zero 
    281             ssdta(:,:,:) = 0.e0 
    282             tsdta(:,:,:) = 0.e0 
    283             vsdta(:,:,:) = 0.e0 
    284             !                                    ! ================== ! 
    285             IF( nobc_dta == 0 )   THEN           ! initial state used 
    286                !                                 ! ================== ! 
    287                !  Fills ssdta, tsdta, vsdta (global arrays) 
    288                !  Remark: this works for njzoom = 1.  
    289                !          Should the definition of ij include njzoom? 
    290                DO jj = njs0, njs1 
    291                   DO jk = 1, jpkm1 
    292                      DO ji = nis0p1, nis1m1 
    293                         ii = ji -1 + nimpp 
    294                         ssdta(ii,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk) 
    295                         tsdta(ii,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk) 
    296                         vsdta(ii,jk,1) = vn(ji,jj,jk)*vmask(ji,jj,jk) 
    297                      END DO 
    298                   END DO 
    299                END DO 
    300             ENDIF 
    301          ENDIF 
    302  
    303       ENDIF        !       end if kt == nit000 
    304        
    305       ! 2.  Initialize the time we are at. 
    306       !     Does this every time the routine is called, 
    307       !     excepted when nobc_dta = 0 
    308       !--------------------------------------------------------------------- 
    309       IF( nobc_dta == 0 )   THEN 
    310          itimo = 1 
    311          zxy   = 0 
    312       ELSE 
    313          IF( ntobc == 1 )   THEN 
    314             itimo = 1 
    315          ELSE IF( ntobc == 12 )   THEN      !   BC are monthly    
    316             ! we assume we have climatology in that case 
    317             iman  = 12 
    318             i15   = nday / 16 
    319             imois = nmonth + i15 - 1 
    320             IF( imois == 0 )   imois = iman 
    321             itimo = imois    
    322          ELSE 
    323             IF(lwp) WRITE(numout,*) 'data other than constant or monthly', kt 
    324             iman  = ntobc 
    325             itimo = FLOOR( kt*rdt / (tcobc(2)-tcobc(1)) ) 
    326             isrel = kt*rdt 
    327          ENDIF 
    328       ENDIF 
    329        
    330       ! 2.1 Read two records in the file if necessary 
    331       ! --------------------------------------------- 
    332       IF( ( nobc_dta == 1 ) .AND. ( ( kt == nit000 .AND. nlecto == 0 ) .OR. itimo  /= ntobc1 ) )   THEN 
    333          nlecto = 1 
    334        
    335          ! Calendar computation 
    336          IF( ntobc == 1 )   THEN            !  BC are constant in time 
    337             ntobc1 = 1 
    338             ntobc2 = 1   
    339          ELSE IF( ntobc == 12 )   THEN      !   BC are monthly    
    340             ntobc1 = itimo         ! first file record used 
    341             ntobc2 = ntobc1 + 1    ! last  file record used 
    342             ntobc1 = MOD( ntobc1, iman ) 
    343             IF( ntobc1 == 0 )   ntobc1 = iman 
    344             ntobc2 = MOD( ntobc2, iman ) 
    345             IF( ntobc2 == 0 )   ntobc2 = iman 
    346             IF( lwp )   THEN 
    347                WRITE(numout,*) ' read monthly obc first record file used ntobc1 ', ntobc1 
    348                WRITE(numout,*) ' read monthly obc last  record file used ntobc2 ', ntobc2 
    349             ENDIF 
    350          ELSE 
    351             isrel=kt*rdt 
    352             ntobc1 = itimo         ! first file record used 
    353             ntobc2 = ntobc1 + 1    ! last  file record used 
    354             ntobc1 = MOD( ntobc1, iman ) 
    355             IF( ntobc1 == 0 )   ntobc1 = iman 
    356             ntobc2 = MOD( ntobc2, iman ) 
    357             IF( ntobc2 == 0 )   ntobc2 = iman 
    358             IF(lwp) WRITE(numout,*) ' read obc first record file used ntobc1 ', ntobc1 
    359             IF(lwp) WRITE(numout,*) ' read obc last  record file used ntobc2 ', ntobc2 
    360          ENDIF 
    361                                ! ======================= ! 
    362                                !  BCs read               ! 
    363          !                     ! ======================= ! 
    364          IF( lp_obc_east )   THEN 
    365             ! ... Read datafile and set temperature, salinity and normal velocity 
    366             ! ... initialise the sedta, tedta, uedta arrays 
    367             CALL iom_open ( 'obceast_TS.nc' , id_e ) 
    368             CALL iom_get ( id_e, jpdom_unknown, 'votemper', tedta(:,:,1), ktime=ntobc1 ) 
    369             CALL iom_get ( id_e, jpdom_unknown, 'votemper', tedta(:,:,2), ktime=ntobc2 ) 
    370             CALL iom_get ( id_e, jpdom_unknown, 'vosaline', sedta(:,:,1), ktime=ntobc1 ) 
    371             CALL iom_get ( id_e, jpdom_unknown, 'vosaline', sedta(:,:,2), ktime=ntobc2 ) 
    372             CALL iom_close (id_e) 
    373             ! 
    374             CALL iom_open ( 'obceast_U.nc' , id_e ) 
    375             CALL iom_get ( id_e, jpdom_unknown, 'vozocrtx', uedta(:,:,1), ktime=ntobc1 ) 
    376             CALL iom_get ( id_e, jpdom_unknown, 'vozocrtx', uedta(:,:,2), ktime=ntobc2 ) 
    377             CALL iom_close ( id_e ) 
    378             !  Usually printout is done only once at kt = nit000, 
    379             !  unless nprint (namelist) > 1       
    380             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN 
    381                WRITE(numout,*) 
    382                WRITE(numout,*) ' Read East OBC data records ', ntobc1, ntobc2 
    383                iprint = (jpjef-jpjed+1)/20 +1 
    384                WRITE(numout,*) ' Temperature record 1 - printout every 3 level' 
    385                CALL prihre( tedta(:,:,1),jpjef-jpjed+1,jpk,1,jpjef-jpjed+1,iprint, & 
    386                   &         jpk, 1, -3, 1., numout ) 
    387                WRITE(numout,*) 
    388                WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
    389                CALL prihre( sedta(:,:,1), jpjef-jpjed+1, jpk, 1, jpjef-jpjed+1, iprint, & 
    390                   &        jpk, 1, -3, 1., numout ) 
    391                WRITE(numout,*) 
    392                WRITE(numout,*) ' Normal velocity U  record 1  - printout every 3 level' 
    393                CALL prihre( uedta(:,:,1), jpjef-jpjed+1, jpk, 1, jpjef-jpjed+1, iprint, & 
    394                   &         jpk, 1, -3, 1., numout ) 
    395             ENDIF 
    396          ENDIF 
    397  
    398          IF( lp_obc_west )   THEN 
    399             ! ... Read datafile and set temperature, salinity and normal velocity 
    400             ! ... initialise the swdta, twdta, uwdta arrays 
    401             CALL iom_open ( 'obcwest_TS.nc' , id_w ) 
    402             CALL iom_get ( id_w, jpdom_unknown, 'votemper', twdta(:,:,1), ktime=ntobc1 ) 
    403             CALL iom_get ( id_w, jpdom_unknown, 'votemper', twdta(:,:,2), ktime=ntobc2 ) 
    404             CALL iom_get ( id_w, jpdom_unknown, 'vosaline', swdta(:,:,1), ktime=ntobc1 ) 
    405             CALL iom_get ( id_w, jpdom_unknown, 'vosaline', swdta(:,:,2), ktime=ntobc2 ) 
    406             CALL iom_close (id_w) 
    407             ! 
    408             CALL iom_open ( 'obcwest_U.nc' , id_w ) 
    409             CALL iom_get ( id_w, jpdom_unknown, 'vozocrtx', uwdta(:,:,1), ktime=ntobc1 ) 
    410             CALL iom_get ( id_w, jpdom_unknown, 'vozocrtx', uwdta(:,:,2), ktime=ntobc2 ) 
    411             CALL iom_close ( id_w ) 
    412             ! 
    413             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN 
    414                WRITE(numout,*) 
    415                WRITE(numout,*) ' Read West OBC data records ', ntobc1, ntobc2 
    416                iprint = (jpjwf-jpjwd+1)/20 +1 
    417                WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
    418                CALL prihre( twdta(:,:,1), jpjwf-jpjwd+1, jpk, 1, jpjwf-jpjwd+1, iprint, & 
    419                   &         jpk, 1, -3, 1., numout ) 
    420                WRITE(numout,*) 
    421                WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
    422                CALL prihre( swdta(:,:,1), jpjwf-jpjwd+1, jpk, 1, jpjwf-jpjwd+1, iprint, & 
    423                   &         jpk, 1, -3, 1., numout ) 
    424                WRITE(numout,*) 
    425                WRITE(numout,*) ' Normal velocity U  record 1  - printout every 3 level' 
    426                CALL prihre( uwdta(:,:,1), jpjwf-jpjwd+1, jpk, 1, jpjwf-jpjwd+1, iprint, & 
    427                   &         jpk, 1, -3, 1., numout ) 
    428             ENDIF 
    429          ENDIF 
    430  
    431          IF( lp_obc_north )   THEN       
    432             CALL iom_open ( 'obcnorth_TS.nc', id_n ) 
    433             CALL iom_get ( id_n, jpdom_unknown, 'votemper', tndta(:,:,1), ktime=ntobc1 ) 
    434             CALL iom_get ( id_n, jpdom_unknown, 'votemper', tndta(:,:,2), ktime=ntobc2 ) 
    435             CALL iom_get ( id_n, jpdom_unknown, 'vosaline', sndta(:,:,1), ktime=ntobc1 ) 
    436             CALL iom_get ( id_n, jpdom_unknown, 'vosaline', sndta(:,:,2), ktime=ntobc2 ) 
    437             CALL iom_close ( id_n )                                                            
    438             ! 
    439             CALL iom_open ( 'obcnorth_V.nc', id_n )                                           
    440             CALL iom_get ( id_n, jpdom_unknown, 'vomecrty', vndta(:,:,1), ktime=ntobc1 ) 
    441             CALL iom_get ( id_n, jpdom_unknown ,'vomecrty', vndta(:,:,2), ktime=ntobc2 ) 
    442             CALL iom_close ( id_n ) 
    443             ! 
    444             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN 
    445                WRITE(numout,*) 
    446                WRITE(numout,*) ' Read North OBC data records ', ntobc1, ntobc2 
    447                iprint = (jpinf-jpind+1)/20 +1 
    448                WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
    449                CALL prihre( tndta(:,:,1), jpinf-jpind+1, jpk, 1, jpinf-jpind+1, iprint, & 
    450                   &         jpk, 1, -3, 1., numout ) 
    451                WRITE(numout,*) 
    452                WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
    453                CALL prihre( sndta(:,:,1), jpinf-jpind+1, jpk, 1, jpinf-jpind+1, iprint, & 
    454                   &         jpk, 1, -3, 1., numout ) 
    455                WRITE(numout,*) 
    456                WRITE(numout,*) ' Normal velocity V  record 1  - printout every 3 level' 
    457                CALL prihre( vndta(:,:,1), jpinf-jpind+1, jpk, 1, jpinf-jpind+1, iprint, & 
    458                   &         jpk, 1, -3, 1., numout ) 
    459             ENDIF 
    460          ENDIF 
    461  
    462          IF( lp_obc_south )   THEN       
    463             CALL iom_open ( 'obcsouth_TS.nc', id_s ) 
    464             CALL iom_get ( id_s, jpdom_unknown, 'votemper', tsdta(:,:,1), ktime=ntobc1 ) 
    465             CALL iom_get ( id_s, jpdom_unknown, 'votemper', tsdta(:,:,2), ktime=ntobc2 ) 
    466             CALL iom_get ( id_s, jpdom_unknown, 'vosaline', ssdta(:,:,1), ktime=ntobc1 ) 
    467             CALL iom_get ( id_s, jpdom_unknown, 'vosaline', ssdta(:,:,2), ktime=ntobc2 ) 
    468             CALL iom_close ( id_s )                                                            
    469             ! 
    470             CALL iom_open ( 'obcsouth_V.nc', id_s )                                           
    471             CALL iom_get ( id_s, jpdom_unknown, 'vomecrty', vsdta(:,:,1), ktime=ntobc1 ) 
    472             CALL iom_get ( id_s, jpdom_unknown ,'vomecrty', vsdta(:,:,2), ktime=ntobc2 ) 
    473             CALL iom_close ( id_s ) 
    474             ! 
    475             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN 
    476                WRITE(numout,*) 
    477                WRITE(numout,*) ' Read South OBC data records ', ntobc1, ntobc2 
    478                iprint = (jpisf-jpisd+1)/20 +1 
    479                WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
    480                CALL prihre( tsdta(:,:,1), jpisf-jpisd+1, jpk, 1, jpisf-jpisd+1, iprint, & 
    481                   &         jpk, 1, -3, 1., numout ) 
    482                WRITE(numout,*) 
    483                WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
    484                CALL prihre( ssdta(:,:,1), jpisf-jpisd+1, jpk, 1, jpisf-jpisd+1, iprint, & 
    485                   &         jpk, 1, -3, 1., numout ) 
    486                WRITE(numout,*) 
    487                WRITE(numout,*) ' Normal velocity V  record 1  - printout every 3 level' 
    488                CALL prihre( vsdta(:,:,1), jpisf-jpisd+1, jpk, 1, jpisf-jpisd+1, iprint, & 
    489                   &         jpk, 1, -3, 1., numout ) 
    490             ENDIF 
    491          ENDIF 
    492  
    493       ELSE 
    494           
    495          nlecto = 0        !      no reading of OBC barotropic data                          
    496  
    497       ENDIF                !      end of the test on the condition to read or not the files  
    498        
    499       ! 3.  Call at every time step : 
    500       !     Linear interpolation of BCs to current time step  
    501       ! ---------------------------------------------------- 
    502  
    503       IF( ntobc == 1 .OR. nobc_dta == 0 )   THEN  
    504          zxy = 0. 
    505       ELSE IF( ntobc == 12 )   THEN          
    506          zxy = FLOAT( nday + 15 - 30 * i15 ) / 30. 
    507       ELSE 
    508          zxy = (tcobc(ntobc1)-FLOAT(isrel))/(tcobc(ntobc1)-tcobc(ntobc2)) 
    509       ENDIF 
    510        
    511       IF( lp_obc_east )   THEN 
    512          !  fills sfoe, tfoe, ufoe (local to each processor) 
    513          DO jk = 1, jpkm1 
    514             DO jj = nje0p1, nje1m1 
    515                ij = jj -1 + njmpp 
    516                sfoe(jj,jk) =  ( zxy * sedta(ij,jk,2) + & 
    517                   &           (1.-zxy) * sedta(ij,jk,1) ) * temsk(jj,jk) 
    518                tfoe(jj,jk) =  ( zxy * tedta(ij,jk,2) + & 
    519                   &           (1.-zxy) * tedta(ij,jk,1) ) * temsk(jj,jk) 
    520                ufoe(jj,jk) =  ( zxy * uedta(ij,jk,2) + & 
    521                   &           (1.-zxy) * uedta(ij,jk,1) ) * uemsk(jj,jk) 
    522             END DO 
    523          END DO 
    524       ENDIF 
    525        
    526       IF( lp_obc_west )   THEN 
    527          !  fills sfow, tfow, ufow (local to each processor) 
    528          DO jk = 1, jpkm1 
    529             DO jj = njw0p1, njw1m1 
    530                ij = jj -1 + njmpp 
    531                sfow(jj,jk) =  ( zxy * swdta(ij,jk,2) + & 
    532                   &           (1.-zxy) * swdta(ij,jk,1) ) * twmsk(jj,jk) 
    533                tfow(jj,jk) =  ( zxy * twdta(ij,jk,2) + & 
    534                   &           (1.-zxy) * twdta(ij,jk,1) ) * twmsk(jj,jk) 
    535                ufow(jj,jk) =  ( zxy * uwdta(ij,jk,2) + & 
    536                   &           (1.-zxy) * uwdta(ij,jk,1) ) * uwmsk(jj,jk) 
    537             END DO 
    538          END DO 
    539       ENDIF 
    540        
    541       IF( lp_obc_north )   THEN 
    542          !  fills sfon, tfon, vfon (local to each processor) 
    543          DO jk = 1, jpkm1 
    544             DO ji = nin0p1, nin1m1 
    545                ii = ji -1 + nimpp 
    546                sfon(ji,jk) =  ( zxy * sndta(ii,jk,2) + & 
    547                   &           (1.-zxy) * sndta(ii,jk,1) ) * tnmsk(ji,jk) 
    548                tfon(ji,jk) =  ( zxy * tndta(ii,jk,2) + & 
    549                   &           (1.-zxy) * tndta(ii,jk,1) ) * tnmsk(ji,jk) 
    550                vfon(ji,jk) =  ( zxy * vndta(ii,jk,2) + & 
    551                   &           (1.-zxy) * vndta(ii,jk,1) ) * vnmsk(ji,jk) 
    552             END DO 
    553          END DO 
    554       ENDIF 
    555        
    556       IF( lp_obc_south )   THEN 
    557          !  fills sfos, tfos, vfos (local to each processor) 
    558          DO jk = 1, jpkm1 
    559            DO ji = nis0p1, nis1m1 
    560               ii = ji -1 + nimpp 
    561               sfos(ji,jk) = ( zxy * ssdta(ii,jk,2) + & 
    562                  &          (1.-zxy) * ssdta(ii,jk,1) ) * tsmsk(ji,jk) 
    563               tfos(ji,jk) = ( zxy * tsdta(ii,jk,2) + & 
    564                  &          (1.-zxy) * tsdta(ii,jk,1) ) * tsmsk(ji,jk) 
    565               vfos(ji,jk) = ( zxy * vsdta(ii,jk,2) + & 
    566                  &          (1.-zxy) * vsdta(ii,jk,1) ) * vsmsk(ji,jk) 
    567            END DO 
    568          END DO              
    569       ENDIF 
    570        
    571    END SUBROUTINE obc_dta 
    572        
     74  SUBROUTINE obc_dta( kt ) 
     75    !!--------------------------------------------------------------------------- 
     76    !!                      ***  SUBROUTINE obc_dta  *** 
     77    !!                     
     78    !! ** Purpose :   Find the climatological  boundary arrays for the specified date,  
     79    !!                The boundary arrays are netcdf files. Three possible cases:  
     80    !!                - one time frame only in the file (time dimension = 1). 
     81    !!                in that case the boundary data does not change in time. 
     82    !!                - many time frames. In that case,  if we have 12 frames 
     83    !!                we assume monthly fields.  
     84    !!                Else, we assume that time_counter is in seconds  
     85    !!                since the beginning of either the current year or a reference 
     86    !!                year given in the namelist. 
     87    !!                (no check is done so far but one would have to check the "unit" 
     88    !!                 attribute of variable time_counter). 
     89    !! 
     90    !! 
     91    !! History : 
     92    !!        !  98-05 (J.M. Molines) Original code 
     93    !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90 
     94    !! 
     95    !!   9.0  !  04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input 
     96    !!        !  2007-2008 (C. Langlais, P. Mathiot, J.M. Molines) high frequency boundaries data 
     97    !!--------------------------------------------------------------------------- 
     98    !! * Arguments 
     99    INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
     100 
     101    !! * Local declarations 
     102    INTEGER ::   ji, jj, jk, ii, ij ,it  ! dummy loop indices 
     103    INTEGER ::  ikprint                                ! frequency for printouts. 
     104    INTEGER, SAVE :: immfile, iyyfile                     ! 
     105    INTEGER :: nt              !  record indices (incrementation) 
     106    INTEGER :: istop           ! local error check 
     107 
     108    REAL(wp) ::   zxy, znum, zden ! time interpolation weight 
     109 
     110    ! variables for the julian day calculation 
     111    INTEGER :: iyear, imonth, iday 
     112    REAL(wp) :: zsec , zjulian, zjuliancnes    
     113 
     114    ! IOM STUFF 
     115    INTEGER ::  idvar, id_e, id_w, id_n, id_s, id_x       ! file identifiers 
     116    INTEGER, DIMENSION(1)  :: itmp 
     117    CHARACTER(LEN=25) :: cl_vname 
     118 
     119    !!--------------------------------------------------------------------------- 
     120    IF( lk_dynspg_rl )  THEN 
     121       CALL obc_dta_psi( kt )     ! update bsf data at open boundaries 
     122       IF ( nobc_dta == 1 .AND. kt == nit000 ) CALL ctl_stop( 'obcdta : time-variable psi boundary data not allowed yet' ) 
     123    ENDIF 
     124 
     125    ! 0.  initialisation : 
     126    ! -------------------- 
     127    IF ( kt == nit000  )  CALL obc_dta_ini ( kt ) 
     128    IF ( nobc_dta == 0 )  RETURN   ! already done in obc_dta_ini 
     129    IF ( itobc == 1    )  RETURN   ! case of only one time frame in file done in obc_dta_ini 
     130 
     131    ! in the following code, we assume that obc data are read from files, with more than 1 time frame in it 
     132 
     133    iyyfile=iyy ; immfile = 00  ! set component of the current file name 
     134    IF ( cffile /= 'annual') immfile = imm   !  
     135    IF ( ln_obc_clim       ) iyyfile = 0000  ! assume that climatological files are labeled y0000 
     136 
     137    ! 1. Synchronize time of run with time of data files 
     138    !--------------------------------------------------- 
     139    ! nday_year is the day number in the current year ( 1 for 01/01 ) 
     140    zsec=MOD( (kt-nit000)*rdt - (nday_year - nday_year0 )*rday, rday ) ! number of seconds in the current day 
     141    IF (ln_obc_clim)  THEN  
     142      zjcnes = nday_year - 1  + zsec/rday 
     143    ELSE 
     144      zjcnes = zjcnes + rdt/rday 
     145    ENDIF 
     146 
     147    ! look for 'before' record number in the current file 
     148    ntobc = nrecbef ()  ! this function return the record number for 'before', relative to zjcnes 
     149 
     150    IF (MOD(kt-1,10)==0) THEN 
     151       IF (lwp) WRITE(numout,*) 'kt= ',kt,' zjcnes =', zjcnes,' ndastp =',ndastp, 'mm =',imm  
     152    END IF 
     153 
     154    ! 2. read a new data if necessary  
     155    !-------------------------------- 
     156    IF ( ntobc /= ntobc_b ) THEN 
     157    ! we need to read the 'after' record 
     158    ! swap working index: 
     159    nt=nt_b ; nt_b=nt_a ; nt_a=nt 
     160    ntobc_b = ntobc 
     161 
     162    ! new record number : 
     163    ntobc_a = ntobc_a + 1  
     164 
     165    ! all tricky things related to record number, changing files etc... are managed by obc_read 
     166      
     167    CALL obc_read (kt, nt_a, ntobc_a, iyyfile, immfile ) 
     168 
     169    ! update zjcnes_obc 
     170    zjcnes_obc(nt_b)= ztcobc(ntobc_b) 
     171    zjcnes_obc(nt_a)= ztcobc(ntobc_a) 
     172    ENDIF 
     173 
     174    ! 3.   interpolation at each time step 
     175    ! ------------------------------------ 
     176    IF ( ln_obc_clim) THEN 
     177      znum= MOD(zjcnes           - zjcnes_obc(nt_b), zraaj ) ; IF ( znum < 0 ) znum = znum + zraaj 
     178      zden= MOD(zjcnes_obc(nt_a) - zjcnes_obc(nt_b), zraaj ) ; IF ( zden < 0 ) zden = zden + zraaj 
     179    ELSE 
     180      znum= zjcnes           - zjcnes_obc(nt_b) 
     181      zden= zjcnes_obc(nt_a) - zjcnes_obc(nt_b) 
     182    ENDIF 
     183    zxy = znum / zden 
     184 
     185    IF( lp_obc_east ) THEN 
     186       !  fills sfoe, tfoe, ufoe ,vfoe 
     187       sfoe(:,:) = zxy * sedta (:,:,nt_a) + (1. - zxy)*sedta(:,:,nt_b) 
     188       tfoe(:,:) = zxy * tedta (:,:,nt_a) + (1. - zxy)*tedta(:,:,nt_b) 
     189       ufoe(:,:) = zxy * uedta (:,:,nt_a) + (1. - zxy)*uedta(:,:,nt_b) 
     190       vfoe(:,:) = zxy * vedta (:,:,nt_a) + (1. - zxy)*vedta(:,:,nt_b) 
     191    ENDIF 
     192 
     193    IF( lp_obc_west) THEN 
     194       !  fills sfow, tfow, ufow ,vfow 
     195       sfow(:,:) = zxy * swdta (:,:,nt_a) + (1. - zxy)*swdta(:,:,nt_b) 
     196       tfow(:,:) = zxy * twdta (:,:,nt_a) + (1. - zxy)*twdta(:,:,nt_b) 
     197       ufow(:,:) = zxy * uwdta (:,:,nt_a) + (1. - zxy)*uwdta(:,:,nt_b) 
     198       vfow(:,:) = zxy * vwdta (:,:,nt_a) + (1. - zxy)*vwdta(:,:,nt_b) 
     199    ENDIF 
     200 
     201    IF( lp_obc_north) THEN 
     202       !  fills sfon, tfon, ufon ,vfon 
     203       sfon(:,:) = zxy * sndta (:,:,nt_a) + (1. - zxy)*sndta(:,:,nt_b) 
     204       tfon(:,:) = zxy * tndta (:,:,nt_a) + (1. - zxy)*tndta(:,:,nt_b) 
     205       ufon(:,:) = zxy * undta (:,:,nt_a) + (1. - zxy)*undta(:,:,nt_b) 
     206       vfon(:,:) = zxy * vndta (:,:,nt_a) + (1. - zxy)*vndta(:,:,nt_b) 
     207    ENDIF 
     208 
     209    IF( lp_obc_south) THEN 
     210       !  fills sfos, tfos, ufos ,vfos 
     211       sfos(:,:) = zxy * ssdta (:,:,nt_a) + (1. - zxy)*ssdta(:,:,nt_b) 
     212       tfos(:,:) = zxy * tsdta (:,:,nt_a) + (1. - zxy)*tsdta(:,:,nt_b) 
     213       ufos(:,:) = zxy * usdta (:,:,nt_a) + (1. - zxy)*usdta(:,:,nt_b) 
     214       vfos(:,:) = zxy * vsdta (:,:,nt_a) + (1. - zxy)*vsdta(:,:,nt_b) 
     215    ENDIF 
     216  END SUBROUTINE obc_dta 
     217 
     218 
     219  SUBROUTINE obc_dta_ini (kt) 
     220    !!----------------------------------------------------------------------------- 
     221    !!                       ***  SUBROUTINE obc_dta_ini  *** 
     222    !! 
     223    !! ** Purpose : 
     224    !!      When obc_dta first call, realize some data initialization 
     225    !! 
     226    !! ** Method : 
     227    !! 
     228    !! History : 
     229    !!   9.0  ! 07-10 (J.M. Molines ) 
     230    !!---------------------------------------------------------------------------- 
     231    !! * Argument 
     232    INTEGER, INTENT(in)  :: kt      ! ocean time-step index 
     233 
     234    !! * Local declarations 
     235    INTEGER ::   ji,jj, it   ! dummy loop indices 
     236 
     237    REAL(wp) ::   zxy                                    ! time interpolation weight 
     238 
     239    INTEGER ::  ikprint                                ! frequency for printouts. 
     240 
     241    INTEGER, SAVE :: immfile, iyyfile                     ! 
     242    INTEGER :: nt              !  record indices (incrementation) 
     243    INTEGER :: istop           ! local error check 
     244 
     245    ! variables for the julian day calculation 
     246    INTEGER :: iyear, imonth, iday 
     247    REAL(wp) :: zsec , zjulian, zjuliancnes    
     248 
     249    ! IOM STUFF 
     250    INTEGER ::  idvar, id_e, id_w, id_n, id_s, id_x       ! file identifiers 
     251    INTEGER, DIMENSION(1)  :: itmp 
     252    CHARACTER(LEN=25) :: cl_vname 
     253 
     254    IF(lwp) WRITE(numout,*) 
     255    IF(lwp) WRITE(numout,*)  'obc_dta : find boundary data' 
     256    IF(lwp) WRITE(numout,*)  '~~~~~~~' 
     257    IF (lwp) THEN 
     258       IF ( nobc_dta == 0 ) THEN  
     259          WRITE(numout,*)  '          OBC data taken from initial conditions.' 
     260       ELSE       
     261          WRITE(numout,*)  '          OBC data taken from netcdf files.' 
     262       ENDIF 
     263    ENDIF 
     264    nday_year0 = nday_year  ! to remember the day when kt=nit000 
     265    SELECT CASE (nleapy) 
     266    CASE ( 0 )  ; zraaj = 365.e0  
     267    CASE ( 1 )  ; zraaj = 366.e0   !! ERROR TO CORRECT:  NOT EVERY YEAR IS LEAP YEAR !! JMM 
     268    CASE DEFAULT; zraaj = 12.* nleapy 
     269    END SELECT 
     270 
     271    sedta(:,:,:) = 0.e0 ; tedta(:,:,:) = 0.e0 ; uedta(:,:,:) = 0.e0 ; vedta(:,:,:) = 0.e0 ! East 
     272    swdta(:,:,:) = 0.e0 ; twdta(:,:,:) = 0.e0 ; uwdta(:,:,:) = 0.e0 ; vwdta(:,:,:) = 0.e0 ! West 
     273    sndta(:,:,:) = 0.e0 ; tndta(:,:,:) = 0.e0 ; undta(:,:,:) = 0.e0 ; vndta(:,:,:) = 0.e0 ! North 
     274    ssdta(:,:,:) = 0.e0 ; tsdta(:,:,:) = 0.e0 ; usdta(:,:,:) = 0.e0 ; vsdta(:,:,:) = 0.e0 ! South 
     275 
     276    sfoe(:,:) = 0.e0  ; tfoe(:,:) = 0.e0 ; ufoe(:,:) = 0.e0 ; vfoe(:,:) = 0.e0   ! East 
     277    sfow(:,:) = 0.e0  ; tfow(:,:) = 0.e0 ; ufow(:,:) = 0.e0 ; vfow(:,:) = 0.e0   ! West 
     278    sfon(:,:) = 0.e0  ; tfon(:,:) = 0.e0 ; ufon(:,:) = 0.e0 ; vfon(:,:) = 0.e0   ! North 
     279    sfos(:,:) = 0.e0  ; tfos(:,:) = 0.e0 ; ufos(:,:) = 0.e0 ; vfos(:,:) = 0.e0   ! South 
     280 
     281    IF (nobc_dta == 0 ) THEN   ! boundary data are the initial data of this run (set only at nit000) 
     282       IF (lp_obc_east) THEN  ! East 
     283          DO ji = nie0 , nie1     
     284             sfoe(nje0p1:nje1m1,:) = temsk(nje0p1:nje1m1,:) * sn (ji+1 , nje0p1:nje1m1 , :) 
     285             tfoe(nje0p1:nje1m1,:) = temsk(nje0p1:nje1m1,:) * tn (ji+1 , nje0p1:nje1m1 , :) 
     286             ufoe(nje0p1:nje1m1,:) = uemsk(nje0p1:nje1m1,:) * un (ji   , nje0p1:nje1m1 , :) 
     287             vfoe(nje0p1:nje1m1,:) = vemsk(nje0p1:nje1m1,:) * vn (ji+1 , nje0p1:nje1m1 , :) 
     288          END DO 
     289       ENDIF 
     290 
     291       IF (lp_obc_west) THEN  ! West 
     292          DO ji = niw0 , niw1     
     293             sfow(njw0p1:njw1m1,:) = twmsk(njw0p1:njw1m1,:) * sn (ji , njw0p1:njw1m1 , :) 
     294             tfow(njw0p1:njw1m1,:) = twmsk(njw0p1:njw1m1,:) * tn (ji , njw0p1:njw1m1 , :) 
     295             ufow(njw0p1:njw1m1,:) = uwmsk(njw0p1:njw1m1,:) * un (ji , njw0p1:njw1m1 , :) 
     296             vfow(njw0p1:njw1m1,:) = vwmsk(njw0p1:njw1m1,:) * vn (ji , njw0p1:njw1m1 , :) 
     297          END DO 
     298       ENDIF 
     299 
     300       IF (lp_obc_north) THEN ! North 
     301          DO jj = njn0 , njn1 
     302             sfon(nin0p1:nin1m1,:) = tnmsk(nin0p1:nin1m1,:) * sn (nin0p1:nin1m1 , jj+1 , :) 
     303             tfon(nin0p1:nin1m1,:) = tnmsk(nin0p1:nin1m1,:) * tn (nin0p1:nin1m1 , jj+1 , :) 
     304             ufon(nin0p1:nin1m1,:) = unmsk(nin0p1:nin1m1,:) * un (nin0p1:nin1m1 , jj+1 , :) 
     305             vfon(nin0p1:nin1m1,:) = vnmsk(nin0p1:nin1m1,:) * vn (nin0p1:nin1m1 , jj   , :) 
     306          END DO 
     307       ENDIF 
     308 
     309       IF (lp_obc_south) THEN ! South 
     310          DO jj = njs0 , njs1 
     311             sfos(nis0p1:nis1m1,:) = tsmsk(nis0p1:nis1m1,:) * sn (nis0p1:nis1m1 , jj , :) 
     312             tfos(nis0p1:nis1m1,:) = tsmsk(nis0p1:nis1m1,:) * tn (nis0p1:nis1m1 , jj , :) 
     313             ufos(nis0p1:nis1m1,:) = usmsk(nis0p1:nis1m1,:) * un (nis0p1:nis1m1 , jj , :) 
     314             vfos(nis0p1:nis1m1,:) = vsmsk(nis0p1:nis1m1,:) * vn (nis0p1:nis1m1 , jj , :) 
     315          END DO 
     316       ENDIF 
     317       RETURN  ! exit the routine all is done 
     318    ENDIF  ! nobc_dta = 0  
     319 
     320!!!! In the following OBC data are read from files. 
     321    ! all logical-mask are initialzed to true when declared 
     322    WHERE ( temsk == 0 ) ltemsk=.FALSE.  
     323    WHERE ( uemsk == 0 ) luemsk=.FALSE.  
     324    WHERE ( vemsk == 0 ) lvemsk=.FALSE.  
     325 
     326    WHERE ( twmsk == 0 ) ltwmsk=.FALSE.  
     327    WHERE ( uwmsk == 0 ) luwmsk=.FALSE.  
     328    WHERE ( vwmsk == 0 ) lvwmsk=.FALSE.  
     329 
     330    WHERE ( tnmsk == 0 ) ltnmsk=.FALSE.  
     331    WHERE ( unmsk == 0 ) lunmsk=.FALSE.  
     332    WHERE ( vnmsk == 0 ) lvnmsk=.FALSE.  
     333 
     334    WHERE ( tsmsk == 0 ) ltsmsk=.FALSE.  
     335    WHERE ( usmsk == 0 ) lusmsk=.FALSE.  
     336    WHERE ( vsmsk == 0 ) lvsmsk=.FALSE.  
     337 
     338    iyear=1950;  imonth=01; iday=01;  zsec=0.  
     339    ! zjuliancnes : julian day corresonding  to  01/01/1950 
     340    CALL ymds2ju(iyear, imonth, iday,zsec , zjuliancnes) 
     341 
     342    !current year and curent month  
     343    iyy=INT(ndastp/10000) ; imm=INT((ndastp -iyy*10000)/100) ; idd=(ndastp-iyy*10000-imm*100) 
     344    IF (iyy <  1900)  iyy = iyy+1900  ! always assume that years are on 4 digits. 
     345    CALL ymds2ju(iyy, imm, idd ,zsec , zjulian) 
     346    ndate0_cnes = zjulian - zjuliancnes   ! jcnes day when call to obc_dta_ini 
     347 
     348    iyyfile=iyy ; immfile=0  ! set component of the current file name 
     349    IF ( cffile /= 'annual') immfile=imm 
     350    IF ( ln_obc_clim) iyyfile = 0  ! assume that climatological files are labeled y0000 
     351 
     352    CALL obc_dta_chktime ( iyyfile, immfile ) 
     353 
     354    IF ( itobc == 1 ) THEN  
     355       ! in this case we will provide boundary data only once. 
     356       nt_a=1 ; ntobc_a=1 
     357       CALL obc_read (nit000, nt_a, ntobc_a, iyyfile, immfile)  
     358       IF( lp_obc_east ) THEN 
     359          !  fills sfoe, tfoe, ufoe ,vfoe 
     360          sfoe(:,:) =  sedta (:,:,1) ; tfoe(:,:) =  tedta (:,:,1) 
     361          ufoe(:,:) =  uedta (:,:,1) ; vfoe(:,:) =  vedta (:,:,1) 
     362       ENDIF 
     363 
     364       IF( lp_obc_west) THEN 
     365          !  fills sfow, tfow, ufow ,vfow 
     366          sfow(:,:) =  swdta (:,:,1) ; tfow(:,:) =  twdta (:,:,1) 
     367          ufow(:,:) =  uwdta (:,:,1) ; vfow(:,:) =  vwdta (:,:,1) 
     368       ENDIF 
     369 
     370       IF( lp_obc_north) THEN 
     371          !  fills sfon, tfon, ufon ,vfon 
     372          sfon(:,:) =  sndta (:,:,1) ; tfon(:,:) =  tndta (:,:,1) 
     373          ufon(:,:) =  undta (:,:,1) ; vfon(:,:) =  vndta (:,:,1) 
     374       ENDIF 
     375 
     376       IF( lp_obc_south) THEN 
     377          !  fills sfos, tfos, ufos ,vfos 
     378          sfos(:,:) =  ssdta (:,:,1) ; tfos(:,:) =  tsdta (:,:,1) 
     379          ufos(:,:) =  usdta (:,:,1) ; vfos(:,:) =  vsdta (:,:,1) 
     380       ENDIF 
     381       RETURN  ! we go out of obc_dta_ini -------------------------------------->>>>> 
     382    ENDIF 
     383 
     384    ! nday_year is the day number in the current year ( 1 for 01/01 ) 
     385    ! we suppose that we always start from the begining of a day 
     386    !    zsec=MOD( (kt-nit000)*rdt - (nday_year - nday_year0 )*rday, rday ) ! number of seconds in the current day 
     387    zsec=0.e0  ! here, kt=nit000, nday_year = ndat_year0  
     388 
     389    IF (ln_obc_clim)  THEN  
     390      zjcnes = nday_year - 1  + zsec/rday  ! for clim file time is in days in a year 
     391    ELSE 
     392      zjcnes = ndate0_cnes + (nday_year - nday_year0 ) + zsec/rday 
     393    ENDIF 
     394 
     395    ! look for 'before' record number in the current file 
     396    ntobc = nrecbef () 
     397 
     398    IF (lwp) WRITE(numout,*) 'obc files frequency :',cffile 
     399    IF (lwp) WRITE(numout,*) ' zjcnes0 =',zjcnes,' ndastp0 =',ndastp 
     400    IF (lwp) WRITE(numout,*) ' annee0 ',iyy,' month0 ', imm,' day0 ', idd 
     401    IF (lwp) WRITE(numout,*) 'first file open :',cl_obc_nTS 
     402 
     403    ! record initialisation 
     404    !-------------------- 
     405    nt_b = 1 ; nt_a = 2 
     406 
     407    ntobc_a = ntobc + 1 
     408    ntobc_b = ntobc 
     409 
     410    CALL obc_read (kt, nt_b, ntobc_b, iyyfile, immfile)  ! read 'before' fields 
     411    CALL obc_read (kt, nt_a, ntobc_a, iyyfile, immfile)  ! read 'after' fields 
     412 
     413    zjcnes_obc(nt_b)= ztcobc(ntobc_b) 
     414    zjcnes_obc(nt_a)= ztcobc(ntobc_a) 
     415    !  
     416  END SUBROUTINE obc_dta_ini 
     417 
     418 
     419  SUBROUTINE obc_dta_chktime (kyyfile, kmmfile) 
     420   ! 
     421   ! check the number of time steps in the files and read ztcobc  
     422   ! 
     423   ! * Arguments 
     424   INTEGER, INTENT(in) :: kyyfile, kmmfile 
     425   ! * local variables 
     426   INTEGER :: istop       ! error control 
     427   INTEGER :: ji          ! dummy loop index 
     428 
     429    INTEGER ::  idvar, id_e, id_w, id_n, id_s, id_x       ! file identifiers 
     430    INTEGER, DIMENSION(1)  :: itmp 
     431    CHARACTER(LEN=25) :: cl_vname 
     432 
     433    ntobc_a = 0; itobce =0 ; itobcw = 0; itobcn = 0; itobcs = 0 
     434    ! build file name 
     435    WRITE(cl_obc_eTS ,'("obc_east_TS_y",i4.4,"m",i2.2,".nc")'  ) kyyfile,kmmfile 
     436    WRITE(cl_obc_wTS ,'("obc_west_TS_y",i4.4,"m",i2.2,".nc")'  ) kyyfile,kmmfile 
     437    WRITE(cl_obc_nTS ,'("obc_north_TS_y",i4.4,"m",i2.2,".nc")' ) kyyfile,kmmfile 
     438    WRITE(cl_obc_sTS ,'("obc_south_TS_y",i4.4,"m",i2.2,".nc")' ) kyyfile,kmmfile 
     439 
     440    cl_vname = 'time_counter' 
     441    IF ( lp_obc_east ) THEN 
     442       CALL iom_open ( cl_obc_eTS , id_e ) 
     443       idvar = iom_varid( id_e, cl_vname, kdimsz = itmp ); itobce=itmp(1) 
     444    ENDIF 
     445    IF ( lp_obc_west ) THEN 
     446       CALL iom_open ( cl_obc_wTS , id_w ) 
     447       idvar = iom_varid( id_w, cl_vname, kdimsz = itmp ) ; itobcw=itmp(1) 
     448    ENDIF 
     449    IF ( lp_obc_north ) THEN 
     450       CALL iom_open ( cl_obc_nTS , id_n ) 
     451       idvar = iom_varid( id_n, cl_vname, kdimsz = itmp ) ; itobcn=itmp(1) 
     452    ENDIF 
     453    IF ( lp_obc_south ) THEN 
     454       CALL iom_open ( cl_obc_sTS , id_s ) 
     455       idvar = iom_varid( id_s, cl_vname, kdimsz = itmp ) ; itobcs=itmp(1) 
     456    ENDIF 
     457 
     458    itobc = MAX( itobce, itobcw, itobcn, itobcs ) 
     459    istop = 0 
     460    IF ( lp_obc_east  .AND. itobce /= itobc ) istop = istop+1  
     461    IF ( lp_obc_west  .AND. itobcw /= itobc ) istop = istop+1       
     462    IF ( lp_obc_north .AND. itobcn /= itobc ) istop = istop+1 
     463    IF ( lp_obc_south .AND. itobcs /= itobc ) istop = istop+1  
     464    nstop = nstop + istop 
     465 
     466    IF ( istop /=  0 )  THEN 
     467       WRITE(ctmp1,*) ' east, west, north, south: ', itobce, itobcw, itobcn, itobcs 
     468       CALL ctl_stop( 'obcdta : all files must have the same number of time steps', ctmp1 ) 
     469    ENDIF 
     470 
     471    IF ( itobc == 1 ) THEN  
     472       IF (lwp) THEN 
     473          WRITE(numout,*) ' obcdta found one time step only in the OBC files' 
     474          IF (ln_obc_clim) THEN 
     475             ! OK no problem 
     476          ELSE 
     477             ln_obc_clim=.true. 
     478             WRITE(numout,*) ' we force ln_obc_clim to T' 
     479          ENDIF 
     480       ENDIF 
     481    ELSE 
     482       IF ( ALLOCATED(ztcobc) ) DEALLOCATE ( ztcobc ) 
     483       ALLOCATE (ztcobc(itobc)) 
     484       DO ji=1,1   ! use a dummy loop to read ztcobc only once 
     485          IF ( lp_obc_east ) THEN 
     486             CALL iom_gettime ( id_e, ztcobc, cl_vname ) ; CALL iom_close (id_e) ; EXIT 
     487          ENDIF 
     488          IF ( lp_obc_west ) THEN 
     489             CALL iom_gettime ( id_w, ztcobc, cl_vname ) ; CALL iom_close (id_w) ; EXIT 
     490          ENDIF 
     491          IF ( lp_obc_north ) THEN 
     492             CALL iom_gettime ( id_n, ztcobc, cl_vname ) ; CALL iom_close (id_n) ; EXIT 
     493          ENDIF 
     494          IF ( lp_obc_south ) THEN 
     495             CALL iom_gettime ( id_s, ztcobc, cl_vname ) ; CALL iom_close (id_s) ; EXIT 
     496          ENDIF 
     497       END DO 
     498       rdt_obc = ztcobc(2)-ztcobc(1)  !  just an information, not used for any computation 
     499       IF (lwp) WRITE(numout,*) ' obcdta found', itobc,' time steps in the OBC files' 
     500       IF (lwp) WRITE(numout,*) ' time step of obc data :', rdt_obc,' days'             
     501     ENDIF 
     502     zjcnes = zjcnes - rdt/rday  ! trick : zcnes is always incremented by rdt/rday in obc_dta! 
     503  END SUBROUTINE obc_dta_chktime 
     504 
    573505# if defined key_dynspg_rl 
    574    !!----------------------------------------------------------------------------- 
    575    !!   Rigid-lid 
    576    !!----------------------------------------------------------------------------- 
    577  
    578    SUBROUTINE obc_dta_psi ( kt ) 
    579       !!----------------------------------------------------------------------------- 
    580       !!                       ***  SUBROUTINE obc_dta_psi  *** 
    581       !! 
    582       !! ** Purpose : 
    583       !!      Update the climatological streamfunction OBC at each time step. 
    584       !!      Depends on the user's configuration.  Here data are read only once 
    585       !!      at the beginning of the run. 
    586       !! 
    587       !! ** Method : 
    588       !!      1. initialization 
    589       !!         kbsfstart: number of time steps over which increase bsf 
    590       !!         during initialization. This is provided for a smooth start 
    591       !!         in cases where the transport is large (like on the Antarctic 
    592       !!         continent). also note that when kbfstart=1, the transport 
    593       !!         increases a lot in one time step and the precision usually 
    594       !!         required for the solver may not be enough. 
    595       !!      2. set the time evolution of the climatological barotropic streamfunction 
    596       !!         along the isolated coastlines ( gcbic(jnic) ).  
    597       !!      3. set the climatological barotropic streamfunction at the boundary. 
    598       !! 
    599       !!      The last two steps are done only at first step (nit000) or if kt <= kbfstart 
    600       !! 
    601       !! History : 
    602       !!        ! 97-08 (G. Madec, J.M. Molines) 
    603       !!   8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 
    604       !!   9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization 
    605       !!---------------------------------------------------------------------------- 
    606       !! * Arguments 
    607       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    608  
    609       !! * Local declarations 
    610       INTEGER ::   ji, jj, jnic, jip         ! dummy loop indices 
    611       INTEGER ::   inum                      ! temporary logical unit 
    612       INTEGER ::   ip, ii, ij, iii, ijj 
    613       INTEGER ::   kbsfstart 
    614       REAL(wp) ::   zsver1, zsver2, zsver3, z2dtr, zcoef 
    615       !!---------------------------------------------------------------------------- 
    616  
    617       ! 1. initialisation 
    618       ! ----------------- 
    619  
    620       kbsfstart =  1 
    621       zsver1 =  bsfic0(1) 
    622       zsver2 =  zsver1 
    623       IF( kt <= kbsfstart )   THEN 
    624          zcoef = float(kt)/float(kbsfstart) 
    625       ELSE 
    626          zcoef = 1. 
    627       END IF 
    628       bsfic(1) = zsver1*zcoef 
    629       IF( lwp .AND. ( kt <= kbsfstart ) )   THEN 
    630          IF(lwp) WRITE(numout,*)'            ' 
    631          IF(lwp) WRITE(numout,*)'obcdta: spinup phase in obc_dta_psi routine' 
    632          IF(lwp) WRITE(numout,*)'~~~~~~  it=',kt,'  OBC: spinup coef: ', & 
    633                                            zcoef, ' and transport: ',bsfic(1) 
    634       END IF 
    635  
    636       zsver2 =  bsfic(1)-bsfic(2) 
    637       zsver3 =  bsfic(2) 
    638  
    639       ! 2. Right hand side of the barotropic elliptic equation (isolated coastlines) 
    640       ! ---------------------------------------------------------------------------- 
    641  
    642       IF( ( neuler == 0 ) .AND. ( kt == nit000 ) )   THEN 
    643          z2dtr = 1./rdt 
    644       ELSE 
    645          z2dtr = 1./2./rdt 
    646       END IF 
    647       ! ... bsfb(ii,ij) should be constant but due to the Asselin filter it 
    648       ! ... converges asymptotically towards bsfic(jnic) 
    649       ! ... However, bsfb(ii,ij) is constant along the same coastlines 
    650       ! ... ---> can be improved using an extra array for storing bsficb (before) 
    651       IF( nbobc > 1 )   THEN 
    652          DO jnic = 1,nbobc - 1 
    653             gcbic(jnic) = 0.e0 
    654             ip=mnic(0,jnic) 
    655             DO jip = 1,ip 
    656                ii = miic(jip,0,jnic)  
    657                ij = mjic(jip,0,jnic)  
    658                IF( ii >= nldi+ nimpp - 1 .AND. ii <= nlei+ nimpp - 1 .AND. & 
    659                    ij >= nldj+ njmpp - 1 .AND. ij <= nlej+ njmpp - 1 )   THEN 
    660                   iii=ii-nimpp+1 
    661                   ijj=ij-njmpp+1 
    662                   gcbic(jnic) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr 
    663                END IF 
    664             END DO 
    665          END DO 
    666       END IF 
    667  
    668       IF( lk_mpp )   CALL mpp_isl( gcbic, 3 ) 
    669  
    670       ! 3. Update the climatological barotropic function at the boundary  
    671       ! ---------------------------------------------------------------- 
    672  
    673       IF( lpeastobc )   THEN 
    674  
    675          IF( kt == nit000 .OR. kt <= kbsfstart )   THEN 
    676             CALL ctlopn( inum, 'obceastbsf.dta', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    677                &         1, numout, lwp, 1 ) 
    678             READ(inum,*) 
    679             READ(inum,*) 
    680             READ(inum,*) 
    681             READ(inum,*) 
    682             READ(inum,*) 
    683             READ(inum,*) (bfoe(jj),jj=jpjed, jpjef) 
    684             CLOSE(inum) 
    685          END IF 
    686          DO jj=jpjed, jpjefm1 
    687             bfoe(jj)=bfoe(jj)*zcoef 
    688          END DO 
    689  
    690       END IF 
    691  
    692       IF( lpwestobc)   THEN 
    693  
    694          IF( kt == nit000 .OR. kt <= kbsfstart ) then 
    695             CALL ctlopn( inum, 'obcwestbsf.dta', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    696                &         1, numout, lwp, 1 ) 
    697             READ(inum,*) 
    698             READ(inum,*) 
    699             READ(inum,*) 
    700             READ(inum,*) 
    701             READ(inum,*) 
    702             READ(inum,*) (bfow(jj),jj=jpjwd, jpjwf) 
    703             CLOSE(inum) 
    704          END IF 
    705          DO jj=jpjwd, jpjwfm1 
    706             bfow(jj)=bfow(jj)*zcoef 
    707          END DO 
    708  
    709       END IF 
    710  
    711       IF( lpsouthobc)   THEN 
    712  
    713          IF( kt == nit000 .OR. kt <= kbsfstart )   THEN 
    714             CALL ctlopn( inum, 'obcsouthbsf.dta', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    715                &         1, numout, lwp, 1 ) 
    716             READ(inum,*) 
    717             READ(inum,*) 
    718             READ(inum,*) 
    719             READ(inum,*) 
    720             READ(inum,*) 
    721             READ(inum,*) (bfos(jj),jj=jpisd, jpisf) 
    722             CLOSE(inum) 
    723          END IF 
    724          DO ji=jpisd, jpisfm1 
    725             bfos(ji)=bfos(ji)*zcoef 
    726          END DO 
    727  
    728       END IF 
    729  
    730       IF( lpnorthobc)   THEN 
    731          IF( kt == nit000 .OR. kt <= kbsfstart )   THEN 
    732             CALL ctlopn( inum, 'obcnorthbsf.dta', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    733                &         1, numout, lwp, 1 ) 
    734             READ(inum,*) 
    735             READ(inum,*) 
    736             READ(inum,*) 
    737             READ(inum,*) 
    738             READ(inum,*) 
    739             READ(inum,*) (bfon(jj),jj=jpind, jpinf) 
    740             CLOSE(inum) 
    741          END IF 
    742          DO ji=jpind, jpinfm1 
    743             bfon(ji)=bfon(ji)*zcoef 
    744          END DO 
    745  
    746       END IF 
    747  
    748    END SUBROUTINE obc_dta_psi 
     506  !!----------------------------------------------------------------------------- 
     507  !!   Rigid-lid 
     508  !!----------------------------------------------------------------------------- 
     509 
     510  SUBROUTINE obc_dta_psi ( kt ) 
     511    !!----------------------------------------------------------------------------- 
     512    !!                       ***  SUBROUTINE obc_dta_psi  *** 
     513    !! 
     514    !! ** Purpose : 
     515    !!      Update the climatological streamfunction OBC at each time step. 
     516    !!      Depends on the user's configuration.  Here data are read only once 
     517    !!      at the beginning of the run. 
     518    !! 
     519    !! ** Method : 
     520    !!      1. initialization 
     521    !!         kbsfstart: number of time steps over which increase bsf 
     522    !!         during initialization. This is provided for a smooth start 
     523    !!         in cases where the transport is large (like on the Antarctic 
     524    !!         continent). also note that when kbfstart=1, the transport 
     525    !!         increases a lot in one time step and the precision usually 
     526    !!         required for the solver may not be enough. 
     527    !!      2. set the time evolution of the climatological barotropic streamfunction 
     528    !!         along the isolated coastlines ( gcbic(jnic) ).  
     529    !!      3. set the climatological barotropic streamfunction at the boundary. 
     530    !! 
     531    !!      The last two steps are done only at first step (nit000) or if kt <= kbfstart 
     532    !! 
     533    !! History : 
     534    !!        ! 97-08 (G. Madec, J.M. Molines) 
     535    !!   8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 
     536    !!---------------------------------------------------------------------------- 
     537    !! * Arguments 
     538    INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
     539 
     540    !! * Local declarations 
     541    INTEGER ::   ji, jj, jnic, jip         ! dummy loop indices 
     542    INTEGER ::   inum = 11                 ! temporary logical unit 
     543    INTEGER ::   ip, ii, ij, iii, ijj 
     544    INTEGER ::   kbsfstart 
     545    REAL(wp) ::   zsver1, zsver2, zsver3, z2dtr, zcoef 
     546    !!---------------------------------------------------------------------------- 
     547 
     548    ! 1. initialisation 
     549    ! ----------------- 
     550    kbsfstart =  1 
     551    zsver1 =  bsfic0(1) 
     552    zsver2 =  zsver1 
     553    IF( kt <= kbsfstart ) THEN 
     554       zcoef = float(kt)/float(kbsfstart) 
     555    ELSE 
     556       zcoef = 1. 
     557    END IF 
     558    bsfic(1) = zsver1*zcoef 
     559    IF( lwp .AND. ( kt <= kbsfstart ) ) THEN 
     560       IF(lwp) WRITE(numout,*)'            ' 
     561       IF(lwp) WRITE(numout,*)'obcdta: spinup phase in obc_dta_psi routine' 
     562       IF(lwp) WRITE(numout,*)'~~~~~~  it=',kt,'  OBC: spinup coef: ', & 
     563            zcoef, ' and transport: ',bsfic(1) 
     564    END IF 
     565 
     566    zsver2 =  bsfic(1)-bsfic(2) 
     567    zsver3 =  bsfic(2) 
     568 
     569    ! 2. Right hand side of the barotropic elliptic equation (isolated coastlines) 
     570    ! ---------------------------------------------------------------------------- 
     571 
     572    IF( ( neuler == 0 ) .AND. ( kt == nit000 ) ) THEN 
     573       z2dtr = 1./rdt 
     574    ELSE 
     575       z2dtr = 1./2./rdt 
     576    END IF 
     577    ! ... bsfb(ii,ij) should be constant but due to the Asselin filter it 
     578    ! ... converges asymptotically towards bsfic(jnic) 
     579    ! ... However, bsfb(ii,ij) is constant along the same coastlines 
     580    ! ... ---> can be improved using an extra array for storing bsficb (before) 
     581    IF( nbobc > 1 ) THEN 
     582       DO jnic = 1,nbobc - 1 
     583          gcbic(jnic) = 0. 
     584          ip=mnic(0,jnic) 
     585          DO jip = 1,ip 
     586             ii = miic(jip,0,jnic)  
     587             ij = mjic(jip,0,jnic)  
     588             IF( ii >= nldi+ nimpp - 1 .AND. ii <= nlei+ nimpp - 1 .AND. & 
     589                  ij >= nldj+ njmpp - 1 .AND. ij <= nlej+ njmpp - 1 ) THEN 
     590                iii=ii-nimpp+1 
     591                ijj=ij-njmpp+1 
     592                gcbic(jnic) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr 
     593             END IF 
     594          END DO 
     595       END DO 
     596    END IF 
     597 
     598    IF( lk_mpp )   CALL mpp_isl( gcbic, 3 ) 
     599 
     600    ! 3. Update the climatological barotropic function at the boundary  
     601    ! ---------------------------------------------------------------- 
     602 
     603    IF( lpeastobc ) THEN 
     604 
     605       IF( kt == nit000 .OR. kt <= kbsfstart ) THEN 
     606          OPEN(inum,file='obceastbsf.dta') 
     607          READ(inum,*) 
     608          READ(inum,*) 
     609          READ(inum,*) 
     610          READ(inum,*) 
     611          READ(inum,*) 
     612          READ(inum,*) (bfoe(jj),jj=jpjed, jpjef) 
     613          CLOSE(inum) 
     614       END IF 
     615       DO jj=jpjed, jpjefm1 
     616          bfoe(jj)=bfoe(jj)*zcoef 
     617       END DO 
     618 
     619    END IF 
     620 
     621    IF( lpwestobc) THEN 
     622 
     623       IF( kt == nit000 .OR. kt <= kbsfstart ) THEN 
     624          OPEN(inum,file='obcwestbsf.dta') 
     625          READ(inum,*) 
     626          READ(inum,*) 
     627          READ(inum,*) 
     628          READ(inum,*) 
     629          READ(inum,*) 
     630          READ(inum,*) (bfow(jj),jj=jpjwd, jpjwf) 
     631          CLOSE(inum) 
     632       END IF 
     633       DO jj=jpjwd, jpjwfm1 
     634          bfow(jj)=bfow(jj)*zcoef 
     635       END DO 
     636 
     637    END IF 
     638 
     639    IF( lpsouthobc) THEN 
     640 
     641       IF( kt == nit000.OR.kt <= kbsfstart ) THEN 
     642          OPEN(inum,file='obcsouthbsf.dta') 
     643          READ(inum,*) 
     644          READ(inum,*) 
     645          READ(inum,*) 
     646          READ(inum,*) 
     647          READ(inum,*) 
     648          READ(inum,*) (bfos(jj),jj=jpisd, jpisf) 
     649          CLOSE(inum) 
     650       END IF 
     651       DO ji=jpisd, jpisfm1 
     652          bfos(ji)=bfos(ji)*zcoef 
     653       END DO 
     654 
     655    END IF 
     656 
     657    IF( lpnorthobc) THEN 
     658       IF( kt == nit000.OR.kt <= kbsfstart ) THEN 
     659          OPEN(inum,file='obcnorthbsf.dta') 
     660          READ(inum,*) 
     661          READ(inum,*) 
     662          READ(inum,*) 
     663          READ(inum,*) 
     664          READ(inum,*) 
     665          READ(inum,*) (bfon(jj),jj=jpind, jpinf) 
     666          CLOSE(inum) 
     667       END IF 
     668       DO ji=jpind, jpinfm1 
     669          bfon(ji)=bfon(ji)*zcoef 
     670       END DO 
     671 
     672    END IF 
     673  END SUBROUTINE obc_dta_psi 
     674 
    749675#else 
    750    !!----------------------------------------------------------------------------- 
    751    !!   Default option                     
    752    !!----------------------------------------------------------------------------- 
    753    SUBROUTINE obc_dta_psi ( kt )       ! Empty routine 
    754       !! * Arguments 
    755       INTEGER,INTENT(in) :: kt  
    756       WRITE(*,*) 'obc_dta_psi: You should not have seen this print! error?', kt 
    757    END SUBROUTINE obc_dta_psi 
     676  !!----------------------------------------------------------------------------- 
     677  !!   'key_dynspg_rl not defined ( dynspg_flt or dynspg_exp or dynspg_ts  ) 
     678  !!----------------------------------------------------------------------------- 
     679  SUBROUTINE obc_dta_psi ( kt )       ! Empty routine 
     680    !! * Arguments 
     681    INTEGER,INTENT(in) :: kt  
     682    WRITE(*,*) 'obc_dta_psi: You should not have seen this print! error?', kt 
     683    STOP 
     684  END SUBROUTINE obc_dta_psi 
    758685# endif 
    759686 
    760  
    761687#if defined key_dynspg_ts || defined key_dynspg_exp 
    762    SUBROUTINE obc_dta_bt( kt, kbt ) 
    763       !!--------------------------------------------------------------------------- 
    764       !!                      ***  SUBROUTINE obc_dta  *** 
    765       !! 
    766       !! ** Purpose :   time interpolation of barotropic data for time-splitting scheme 
    767       !!                Data at the boundary must be in m2/s  
    768       !! 
    769       !! History : 
    770       !!   9.0  !  05-11 (V. garnier) Original code 
    771       !!--------------------------------------------------------------------------- 
    772       !! * Arguments 
    773       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    774       INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index 
    775  
    776       !! * Local declarations 
    777       INTEGER ::   ji, jj, jk, ii, ij   ! dummy loop indices 
    778       INTEGER ::   id_e, id_w, id_n, id_s, fid  ! file identifiers 
    779       INTEGER ::   itimo, iman, imois, i15 
    780       INTEGER ::   itobcm, itobcp, itimom, itimop 
    781       REAL(wp) ::  zxy 
    782       INTEGER ::   isrel, ikt           ! number of seconds since 1/1/1992 
    783       INTEGER ::   iprint              ! frequency for printouts. 
    784  
    785       !!--------------------------------------------------------------------------- 
    786  
    787       ! 1.   First call: check time frames available in files. 
    788       ! ------------------------------------------------------- 
    789  
    790       IF( kt == nit000 ) THEN 
    791  
    792          ! 1.1  Barotropic tangential velocities set to zero 
    793          ! ------------------------------------------------- 
    794          IF( lp_obc_east  ) vbtfoe(:) = 0.e0 
    795          IF( lp_obc_west  ) vbtfow(:) = 0.e0 
    796          IF( lp_obc_south ) ubtfos(:) = 0.e0 
    797          IF( lp_obc_north ) ubtfon(:) = 0.e0 
    798  
    799          ! 1.2  Sea surface height and normal barotropic velocities set to zero 
    800          !                               or initial conditions if nobc_dta == 0 
    801          ! -------------------------------------------------------------------- 
    802  
    803           IF( lp_obc_east ) THEN 
    804              ! initialisation to zero 
    805              sshedta(:,:) = 0.e0 
    806              ubtedta(:,:) = 0.e0 
    807              !                                        ! ================== ! 
    808              IF( nobc_dta == 0 )   THEN               ! initial state used ! 
    809                 !                                     ! ================== ! 
    810                 !  Fills sedta, tedta, uedta (global arrays) 
    811                 !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
    812                 DO ji = nie0, nie1 
    813                    DO jj = nje0p1, nje1m1 
    814                       ij = jj -1 + njmpp 
    815                       sshedta(ij,1) = sshn(ji+1,jj) * tmask(ji+1,jj,1) 
     688  SUBROUTINE obc_dta_bt( kt, kbt ) 
     689    !!--------------------------------------------------------------------------- 
     690    !!                      ***  SUBROUTINE obc_dta  *** 
     691    !! 
     692    !! ** Purpose :   time interpolation of barotropic data for time-splitting scheme 
     693    !!                Data at the boundary must be in m2/s  
     694    !! 
     695    !! History : 
     696    !!   9.0  !  05-11 (V. garnier) Original code 
     697    !!--------------------------------------------------------------------------- 
     698    !! * Arguments 
     699    INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
     700    INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index 
     701 
     702    !! * Local declarations 
     703    INTEGER ::   ji, jj, jk, ii, ij   ! dummy loop indices 
     704    INTEGER ::   id_e, id_w, id_n, id_s, fid  ! file identifiers 
     705    INTEGER ::   itimo, iman, imois, i15 
     706    INTEGER ::   itobcm, itobcp, itimom, itimop 
     707    REAL(wp) ::  zxy 
     708    INTEGER ::   isrel, ikt           ! number of seconds since 1/1/1992 
     709    INTEGER ::   iprint              ! frequency for printouts. 
     710 
     711    !!--------------------------------------------------------------------------- 
     712 
     713    ! 1.   First call: check time frames available in files. 
     714    ! ------------------------------------------------------- 
     715 
     716    IF( kt == nit000 ) THEN 
     717 
     718       ! 1.1  Barotropic tangential velocities set to zero 
     719       ! ------------------------------------------------- 
     720       IF( lp_obc_east  ) vbtfoe(:) = 0.e0 
     721       IF( lp_obc_west  ) vbtfow(:) = 0.e0 
     722       IF( lp_obc_south ) ubtfos(:) = 0.e0 
     723       IF( lp_obc_north ) ubtfon(:) = 0.e0 
     724 
     725       ! 1.2  Sea surface height and normal barotropic velocities set to zero 
     726       !                               or initial conditions if nobc_dta == 0 
     727       ! -------------------------------------------------------------------- 
     728 
     729       IF( lp_obc_east ) THEN 
     730          ! initialisation to zero 
     731          sshedta(:,:) = 0.e0 
     732          ubtedta(:,:) = 0.e0 
     733          !                                        ! ================== ! 
     734          IF( nobc_dta == 0 )   THEN               ! initial state used ! 
     735             !                                     ! ================== ! 
     736             !  Fills sedta, tedta, uedta (global arrays) 
     737             !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
     738             DO ji = nie0, nie1 
     739                DO jj = nje0p1, nje1m1 
     740                   ij = jj -1 + njmpp 
     741                   sshedta(ij,1) = sshn(ji+1,jj) * tmask(ji+1,jj,1) 
     742                END DO 
     743             END DO 
     744          ENDIF 
     745       ENDIF 
     746 
     747       IF( lp_obc_west) THEN 
     748          ! initialisation to zero 
     749          sshwdta(:,:) = 0.e0 
     750          ubtwdta(:,:) = 0.e0 
     751          !                                        ! ================== ! 
     752          IF( nobc_dta == 0 )   THEN               ! initial state used ! 
     753             !                                     ! ================== ! 
     754             !  Fills swdta, twdta, uwdta (global arrays) 
     755             !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
     756             DO ji = niw0, niw1 
     757                DO jj = njw0p1, njw1m1 
     758                   ij = jj -1 + njmpp 
     759                   sshwdta(ij,1) = sshn(ji,jj) * tmask(ji,jj,1) 
     760                END DO 
     761             END DO 
     762          ENDIF 
     763       ENDIF 
     764 
     765       IF( lp_obc_north) THEN 
     766          ! initialisation to zero 
     767          sshndta(:,:) = 0.e0 
     768          vbtndta(:,:) = 0.e0 
     769          !                                        ! ================== ! 
     770          IF( nobc_dta == 0 )   THEN               ! initial state used ! 
     771             !                                     ! ================== ! 
     772             !  Fills sndta, tndta, vndta (global arrays) 
     773             !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
     774             DO jj = njn0, njn1 
     775                DO ji = nin0p1, nin1m1 
     776                   DO jk = 1, jpkm1 
     777                      ii = ji -1 + nimpp 
     778                      vbtndta(ii,1) = vbtndta(ii,1) + vndta(ii,jk,1)*fse3v(ji,jj,jk) 
    816779                   END DO 
     780                   sshndta(ii,1) = sshn(ii,jj+1) * tmask(ji,jj+1,1) 
    817781                END DO 
    818              ENDIF 
    819           ENDIF 
    820  
    821           IF( lp_obc_west) THEN 
    822              ! initialisation to zero 
    823              sshwdta(:,:) = 0.e0 
    824              ubtwdta(:,:) = 0.e0 
    825              !                                        ! ================== ! 
    826              IF( nobc_dta == 0 )   THEN               ! initial state used ! 
    827                 !                                     ! ================== ! 
    828                 !  Fills swdta, twdta, uwdta (global arrays) 
    829                 !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
    830                 DO ji = niw0, niw1 
    831                    DO jj = njw0p1, njw1m1 
    832                       ij = jj -1 + njmpp 
    833                       sshwdta(ij,1) = sshn(ji,jj) * tmask(ji,jj,1) 
     782             END DO 
     783          ENDIF 
     784       ENDIF 
     785 
     786       IF( lp_obc_south) THEN 
     787          ! initialisation to zero 
     788          ssdta(:,:,:) = 0.e0 
     789          tsdta(:,:,:) = 0.e0 
     790          vsdta(:,:,:) = 0.e0 
     791          sshsdta(:,:) = 0.e0 
     792          vbtsdta(:,:) = 0.e0 
     793          !                                        ! ================== ! 
     794          IF( nobc_dta == 0 )   THEN               ! initial state used ! 
     795             !                                     ! ================== ! 
     796             !  Fills ssdta, tsdta, vsdta (global arrays) 
     797             !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
     798             DO jj = njs0, njs1 
     799                DO ji = nis0p1, nis1m1 
     800                   DO jk = 1, jpkm1 
     801                      ii = ji -1 + nimpp 
     802                      vbtsdta(ii,1) = vbtsdta(ii,1) + vsdta(ii,jk,1)*fse3v(ji,jj,jk) 
    834803                   END DO 
     804                   sshsdta(ii,1) = sshn(ji,jj) * tmask(ii,jj,1) 
    835805                END DO 
    836              ENDIF 
    837           ENDIF 
    838  
    839           IF( lp_obc_north) THEN 
    840              ! initialisation to zero 
    841              sshndta(:,:) = 0.e0 
    842              vbtndta(:,:) = 0.e0 
    843              !                                        ! ================== ! 
    844              IF( nobc_dta == 0 )   THEN               ! initial state used ! 
    845                 !                                     ! ================== ! 
    846                 !  Fills sndta, tndta, vndta (global arrays) 
    847                 !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
    848                 DO jj = njn0, njn1 
    849                    DO ji = nin0p1, nin1m1 
    850                       DO jk = 1, jpkm1 
    851                          ii = ji -1 + nimpp 
    852                          vbtndta(ii,1) = vbtndta(ii,1) + vndta(ii,jk,1)*fse3v(ji,jj,jk) 
    853                       END DO 
    854                       sshndta(ii,1) = sshn(ii,jj+1) * tmask(ji,jj+1,1) 
    855                    END DO 
    856                 END DO 
    857              ENDIF 
    858           ENDIF 
    859  
    860           IF( lp_obc_south) THEN 
    861              ! initialisation to zero 
    862              ssdta(:,:,:) = 0.e0 
    863              tsdta(:,:,:) = 0.e0 
    864              vsdta(:,:,:) = 0.e0 
    865              sshsdta(:,:) = 0.e0 
    866              vbtsdta(:,:) = 0.e0 
    867              !                                        ! ================== ! 
    868              IF( nobc_dta == 0 )   THEN               ! initial state used ! 
    869                 !                                     ! ================== ! 
    870                 !  Fills ssdta, tsdta, vsdta (global arrays) 
    871                 !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom? 
    872                 DO jj = njs0, njs1 
    873                    DO ji = nis0p1, nis1m1 
    874                       DO jk = 1, jpkm1 
    875                          ii = ji -1 + nimpp 
    876                          vbtsdta(ii,1) = vbtsdta(ii,1) + vsdta(ii,jk,1)*fse3v(ji,jj,jk) 
    877                       END DO 
    878                       sshsdta(ii,1) = sshn(ji,jj) * tmask(ii,jj,1) 
    879                    END DO 
    880                 END DO 
    881              ENDIF 
    882           ENDIF 
    883  
    884        ENDIF        !       END IF kt == nit000 
    885  
    886       !!------------------------------------------------------------------------------------ 
    887       ! 2.      Initialize the time we are at. Does this every time the routine is called, 
    888       !         excepted when nobc_dta = 0 
    889       ! 
    890       IF( nobc_dta == 0) THEN 
    891          itimo = 1 
    892          zxy   = 0 
    893       ELSE 
    894          IF(ntobc == 1) THEN 
    895             itimo = 1 
    896          ELSE IF (ntobc == 12) THEN      !   BC are monthly 
    897             ! we assume we have climatology in that case 
    898             iman  = 12 
    899             i15   = nday / 16 
    900             imois = nmonth + i15 - 1 
    901             IF( imois == 0 )   imois = iman 
    902             itimo = imois 
    903          ELSE 
    904             IF(lwp) WRITE(numout,*) 'data other than constant or monthly',kt 
    905             iman  = ntobc 
    906             itimo = FLOOR( kt*rdt / tcobc(1)) 
    907             isrel=kt*rdt 
    908          ENDIF 
    909       ENDIF 
    910  
    911       ! 2. Read two records in the file if necessary 
    912       ! --------------------------------------------- 
    913  
    914       IF( nobc_dta == 1 .AND. nlecto == 1 ) THEN 
    915  
    916          IF( lp_obc_east ) THEN 
    917             ! ... Read datafile and set sea surface height and barotropic velocity 
    918             ! ... initialise the sshedta, ubtedta arrays 
    919             sshedta(:,0) = sshedta(:,1) 
    920             ubtedta(:,0) = ubtedta(:,1) 
    921             CALL iom_open ( 'obceast_TS.nc', id_e ) 
    922             CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,1), ktime=ntobc1 ) 
    923             CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,2), ktime=ntobc2 ) 
    924             IF( lk_dynspg_ts ) THEN 
    925                CALL iom_get (id_e, jpdom_unknown, 'vossurfh', sshedta(:,3), ktime=ntobc2+1 ) 
    926             ENDIF 
    927             CALL iom_close ( id_e ) 
    928             ! 
    929             CALL iom_open ( 'obceast_U.nc', id_e ) 
    930             CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,1), ktime=ntobc1 ) 
    931             CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,2), ktime=ntobc2 ) 
    932             IF( lk_dynspg_ts ) THEN 
    933                CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,3), ktime=ntobc2+1 ) 
    934             ENDIF 
    935             CALL iom_close ( id_e ) 
    936             ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
    937             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
    938                WRITE(numout,*) 
    939                WRITE(numout,*) ' Read East OBC barotropic data records ', ntobc1, ntobc2 
    940                iprint = (jpjef-jpjed+1)/20 +1 
    941                WRITE(numout,*) 
    942                WRITE(numout,*) ' Sea surface height record 1' 
    943                CALL prihre( sshedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout ) 
    944                WRITE(numout,*) 
    945                WRITE(numout,*) ' Normal transport (m2/s) record 1',iprint 
    946                CALL prihre( ubtedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout ) 
    947             ENDIF 
    948          ENDIF 
    949  
    950          IF( lp_obc_west ) THEN 
    951             ! ... Read datafile and set temperature, salinity and normal velocity 
    952             ! ... initialise the swdta, twdta, uwdta arrays 
    953             sshwdta(:,0) = sshwdta(:,1) 
    954             ubtwdta(:,0) = ubtwdta(:,1) 
    955             CALL iom_open ( 'obcwest_TS.nc', id_w ) 
    956             CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,1), ktime=ntobc1 ) 
    957             CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,2), ktime=ntobc2 ) 
    958             IF( lk_dynspg_ts ) THEN 
    959                CALL  ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,3), ktime=ntobc2+1 ) 
    960             ENDIF 
    961             CALL iom_close ( id_w ) 
    962             ! 
    963             CALL iom_open ( 'obcwest_U.nc', id_w ) 
    964             CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,1), ktime=ntobc1 ) 
    965             CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,2), ktime=ntobc2 ) 
    966             IF( lk_dynspg_ts ) THEN 
    967                CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,3), ktime=ntobc2+1 ) 
    968             ENDIF 
    969             CALL iom_close ( id_w ) 
    970             ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
    971             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
    972                WRITE(numout,*) 
    973                WRITE(numout,*) ' Read West OBC barotropic data records ', ntobc1, ntobc2 
    974                iprint = (jpjwf-jpjwd+1)/20 +1 
    975                WRITE(numout,*) 
    976                WRITE(numout,*) ' Sea surface height record 1  - printout surface level' 
    977                CALL prihre( sshwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout ) 
    978                WRITE(numout,*) 
    979                WRITE(numout,*) ' Normal transport (m2/s) record 1' 
    980                CALL prihre( ubtwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout ) 
    981             ENDIF 
    982          ENDIF 
    983  
    984          IF( lp_obc_north) THEN 
    985             ! ... Read datafile and set sea surface height and barotropic velocity 
    986             ! ... initialise the sshndta, ubtndta arrays 
    987             sshndta(:,0) = sshndta(:,1) 
    988             vbtndta(:,0) = vbtndta(:,1) 
    989             CALL iom_open ( 'obcnorth_TS.nc', id_n ) 
    990             CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,1), ktime=ntobc1 ) 
    991             CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,2), ktime=ntobc2 ) 
    992             IF( lk_dynspg_ts ) THEN 
    993                CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,3), ktime=ntobc2+1 ) 
    994             ENDIF 
    995             CALL iom_close ( id_n ) 
    996  
    997             CALL iom_open ( 'obcnorth_V.nc', id_n ) 
    998             CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,1), ktime=ntobc1 ) 
    999             CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,2), ktime=ntobc2 ) 
    1000             IF( lk_dynspg_ts ) THEN 
    1001                CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,3), ktime=ntobc2+1 ) 
    1002             ENDIF 
    1003             CALL iom_close ( id_n ) 
    1004  
    1005             ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
    1006             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
    1007                WRITE(numout,*) 
    1008                WRITE(numout,*) ' Read North OBC barotropic data records ', ntobc1, ntobc2 
    1009                iprint = (jpinf-jpind+1)/20 +1 
    1010                WRITE(numout,*) 
    1011                WRITE(numout,*) ' Sea surface height record 1  - printout surface level' 
    1012                CALL prihre( sshndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout ) 
    1013                WRITE(numout,*) 
    1014                WRITE(numout,*) ' Normal transport (m2/s) record 1' 
    1015                CALL prihre( vbtndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout ) 
    1016             ENDIF 
    1017          ENDIF 
    1018  
    1019          IF( lp_obc_south) THEN 
    1020             ! ... Read datafile and set sea surface height and barotropic velocity 
    1021             ! ... initialise the sshsdta, ubtsdta arrays 
    1022             sshsdta(:,0) = sshsdta(:,1) 
    1023             vbtsdta(:,0) = vbtsdta(:,1) 
    1024             CALL iom_open ( 'obcsouth_TS.nc', id_s ) 
    1025             CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,1), ktime=ntobc1 ) 
    1026             CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,2), ktime=ntobc2 ) 
    1027             IF( lk_dynspg_ts ) THEN 
    1028                CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,3), ktime=ntobc2+1 ) 
    1029             ENDIF 
    1030             CALL iom_close ( id_s ) 
    1031  
    1032             CALL iom_open ( 'obcsouth_V.nc', id_s ) 
    1033             CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,1), ktime=ntobc1 ) 
    1034             CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,2), ktime=ntobc2 ) 
    1035             IF( lk_dynspg_ts ) THEN 
    1036                CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,3), ktime=ntobc2+1 ) 
    1037             ENDIF 
    1038             CALL iom_close ( id_s ) 
    1039                  
    1040             ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
    1041             IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
    1042                WRITE(numout,*) 
    1043                WRITE(numout,*) ' Read South OBC barotropic data records ', ntobc1, ntobc2 
    1044                iprint = (jpisf-jpisd+1)/20 +1 
    1045                WRITE(numout,*) 
    1046                WRITE(numout,*) ' Sea surface height record 1  - printout surface level' 
    1047                CALL prihre( sshsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout ) 
    1048                WRITE(numout,*) 
    1049                WRITE(numout,*) ' Normal transport (m2/s) record 1' 
    1050                CALL prihre( vbtsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout ) 
    1051             ENDIF 
    1052          ENDIF 
    1053  
    1054        ENDIF        !      end of the test on the condition to read or not the files 
    1055  
    1056       ! 3.  Call at every time step : Linear interpolation of BCs to current time step 
    1057       ! ---------------------------------------------------------------------- 
    1058  
    1059        IF( lk_dynspg_ts ) THEN 
     806             END DO 
     807          ENDIF 
     808       ENDIF 
     809 
     810    ENDIF        !       END IF kt == nit000 
     811 
     812    !!------------------------------------------------------------------------------------ 
     813    ! 2.      Initialize the time we are at. Does this every time the routine is called, 
     814    !         excepted when nobc_dta = 0 
     815    ! 
     816    IF( nobc_dta == 0) THEN 
     817       itimo = 1 
     818       zxy   = 0 
     819    ELSE 
     820       IF(ntobc == 1) THEN 
     821          itimo = 1 
     822       ELSE IF (ntobc == 12) THEN      !   BC are monthly 
     823          ! we assume we have climatology in that case 
     824          iman  = 12 
     825          i15   = nday / 16 
     826          imois = nmonth + i15 - 1 
     827          IF( imois == 0 )   imois = iman 
     828          itimo = imois 
     829       ELSE 
     830          IF(lwp) WRITE(numout,*) 'data other than constant or monthly',kt 
     831          iman  = ntobc 
     832          itimo = FLOOR( kt*rdt / tcobc(1)) 
     833          isrel=kt*rdt 
     834       ENDIF 
     835    ENDIF 
     836 
     837    ! 2. Read two records in the file if necessary 
     838    ! --------------------------------------------- 
     839 
     840    IF( nobc_dta == 1 .AND. nlecto == 1 ) THEN 
     841 
     842       IF( lp_obc_east ) THEN 
     843          ! ... Read datafile and set sea surface height and barotropic velocity 
     844          ! ... initialise the sshedta, ubtedta arrays 
     845          sshedta(:,0) = sshedta(:,1) 
     846          ubtedta(:,0) = ubtedta(:,1) 
     847          CALL iom_open ( 'obceast_TS.nc', id_e ) 
     848          CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,1), ktime=ntobc1 ) 
     849          CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,2), ktime=ntobc2 ) 
     850          IF( lk_dynspg_ts ) THEN 
     851             CALL iom_get (id_e, jpdom_unknown, 'vossurfh', sshedta(:,3), ktime=ntobc2+1 ) 
     852          ENDIF 
     853          CALL iom_close ( id_e ) 
     854          ! 
     855          CALL iom_open ( 'obceast_U.nc', id_e ) 
     856          CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,1), ktime=ntobc1 ) 
     857          CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,2), ktime=ntobc2 ) 
     858          IF( lk_dynspg_ts ) THEN 
     859             CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,3), ktime=ntobc2+1 ) 
     860          ENDIF 
     861          CALL iom_close ( id_e ) 
     862          ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
     863          IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     864             WRITE(numout,*) 
     865             WRITE(numout,*) ' Read East OBC barotropic data records ', ntobc1, ntobc2 
     866             iprint = (jpjef-jpjed+1)/20 +1 
     867             WRITE(numout,*) 
     868             WRITE(numout,*) ' Sea surface height record 1' 
     869             CALL prihre( sshedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout ) 
     870             WRITE(numout,*) 
     871             WRITE(numout,*) ' Normal transport (m2/s) record 1',iprint 
     872             CALL prihre( ubtedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout ) 
     873          ENDIF 
     874       ENDIF 
     875 
     876       IF( lp_obc_west ) THEN 
     877          ! ... Read datafile and set temperature, salinity and normal velocity 
     878          ! ... initialise the swdta, twdta, uwdta arrays 
     879          sshwdta(:,0) = sshwdta(:,1) 
     880          ubtwdta(:,0) = ubtwdta(:,1) 
     881          CALL iom_open ( 'obcwest_TS.nc', id_w ) 
     882          CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,1), ktime=ntobc1 ) 
     883          CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,2), ktime=ntobc2 ) 
     884          IF( lk_dynspg_ts ) THEN 
     885             CALL  ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,3), ktime=ntobc2+1 ) 
     886          ENDIF 
     887          CALL iom_close ( id_w ) 
     888          ! 
     889          CALL iom_open ( 'obcwest_U.nc', id_w ) 
     890          CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,1), ktime=ntobc1 ) 
     891          CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,2), ktime=ntobc2 ) 
     892          IF( lk_dynspg_ts ) THEN 
     893             CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,3), ktime=ntobc2+1 ) 
     894          ENDIF 
     895          CALL iom_close ( id_w ) 
     896          ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
     897          IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     898             WRITE(numout,*) 
     899             WRITE(numout,*) ' Read West OBC barotropic data records ', ntobc1, ntobc2 
     900             iprint = (jpjwf-jpjwd+1)/20 +1 
     901             WRITE(numout,*) 
     902             WRITE(numout,*) ' Sea surface height record 1  - printout surface level' 
     903             CALL prihre( sshwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout ) 
     904             WRITE(numout,*) 
     905             WRITE(numout,*) ' Normal transport (m2/s) record 1' 
     906             CALL prihre( ubtwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout ) 
     907          ENDIF 
     908       ENDIF 
     909 
     910       IF( lp_obc_north) THEN 
     911          ! ... Read datafile and set sea surface height and barotropic velocity 
     912          ! ... initialise the sshndta, ubtndta arrays 
     913          sshndta(:,0) = sshndta(:,1) 
     914          vbtndta(:,0) = vbtndta(:,1) 
     915          CALL iom_open ( 'obcnorth_TS.nc', id_n ) 
     916          CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,1), ktime=ntobc1 ) 
     917          CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,2), ktime=ntobc2 ) 
     918          IF( lk_dynspg_ts ) THEN 
     919             CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,3), ktime=ntobc2+1 ) 
     920          ENDIF 
     921          CALL iom_close ( id_n ) 
     922 
     923          CALL iom_open ( 'obcnorth_V.nc', id_n ) 
     924          CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,1), ktime=ntobc1 ) 
     925          CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,2), ktime=ntobc2 ) 
     926          IF( lk_dynspg_ts ) THEN 
     927             CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,3), ktime=ntobc2+1 ) 
     928          ENDIF 
     929          CALL iom_close ( id_n ) 
     930 
     931          ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
     932          IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     933             WRITE(numout,*) 
     934             WRITE(numout,*) ' Read North OBC barotropic data records ', ntobc1, ntobc2 
     935             iprint = (jpinf-jpind+1)/20 +1 
     936             WRITE(numout,*) 
     937             WRITE(numout,*) ' Sea surface height record 1  - printout surface level' 
     938             CALL prihre( sshndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout ) 
     939             WRITE(numout,*) 
     940             WRITE(numout,*) ' Normal transport (m2/s) record 1' 
     941             CALL prihre( vbtndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout ) 
     942          ENDIF 
     943       ENDIF 
     944 
     945       IF( lp_obc_south) THEN 
     946          ! ... Read datafile and set sea surface height and barotropic velocity 
     947          ! ... initialise the sshsdta, ubtsdta arrays 
     948          sshsdta(:,0) = sshsdta(:,1) 
     949          vbtsdta(:,0) = vbtsdta(:,1) 
     950          CALL iom_open ( 'obcsouth_TS.nc', id_s ) 
     951          CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,1), ktime=ntobc1 ) 
     952          CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,2), ktime=ntobc2 ) 
     953          IF( lk_dynspg_ts ) THEN 
     954             CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,3), ktime=ntobc2+1 ) 
     955          ENDIF 
     956          CALL iom_close ( id_s ) 
     957 
     958          CALL iom_open ( 'obcsouth_V.nc', id_s ) 
     959          CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,1), ktime=ntobc1 ) 
     960          CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,2), ktime=ntobc2 ) 
     961          IF( lk_dynspg_ts ) THEN 
     962             CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,3), ktime=ntobc2+1 ) 
     963          ENDIF 
     964          CALL iom_close ( id_s ) 
     965 
     966          ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1 
     967          IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     968             WRITE(numout,*) 
     969             WRITE(numout,*) ' Read South OBC barotropic data records ', ntobc1, ntobc2 
     970             iprint = (jpisf-jpisd+1)/20 +1 
     971             WRITE(numout,*) 
     972             WRITE(numout,*) ' Sea surface height record 1  - printout surface level' 
     973             CALL prihre( sshsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout ) 
     974             WRITE(numout,*) 
     975             WRITE(numout,*) ' Normal transport (m2/s) record 1' 
     976             CALL prihre( vbtsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout ) 
     977          ENDIF 
     978       ENDIF 
     979 
     980    ENDIF        !      end of the test on the condition to read or not the files 
     981 
     982    ! 3.  Call at every time step : Linear interpolation of BCs to current time step 
     983    ! ---------------------------------------------------------------------- 
     984 
     985    IF( lk_dynspg_ts ) THEN 
     986       isrel = (kt-1)*rdt + kbt*rdtbt 
     987 
     988       IF( nobc_dta == 1 ) THEN 
    1060989          isrel = (kt-1)*rdt + kbt*rdtbt 
    1061  
    1062           IF( nobc_dta == 1 ) THEN 
    1063              isrel = (kt-1)*rdt + kbt*rdtbt 
    1064              itimo  = FLOOR(  kt*rdt    / (tcobc(2)-tcobc(1)) ) 
    1065              itimom = FLOOR( (kt-1)*rdt / (tcobc(2)-tcobc(1)) ) 
    1066              itimop = FLOOR( (kt+1)*rdt / (tcobc(2)-tcobc(1)) ) 
    1067              IF( itimom == itimo .AND. itimop == itimo ) THEN 
     990          itimo  = FLOOR(  kt*rdt    / (tcobc(2)-tcobc(1)) ) 
     991          itimom = FLOOR( (kt-1)*rdt / (tcobc(2)-tcobc(1)) ) 
     992          itimop = FLOOR( (kt+1)*rdt / (tcobc(2)-tcobc(1)) ) 
     993          IF( itimom == itimo .AND. itimop == itimo ) THEN 
     994             itobcm = ntobc1 
     995             itobcp = ntobc2 
     996 
     997          ELSEIF ( itimom <= itimo .AND. itimop == itimo ) THEN 
     998             IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN 
     999                itobcm = ntobc1-1 
     1000                itobcp = ntobc2-1 
     1001             ELSE 
    10681002                itobcm = ntobc1 
    10691003                itobcp = ntobc2 
    1070  
    1071              ELSEIF ( itimom <= itimo .AND. itimop == itimo ) THEN 
    1072                 IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN 
    1073                    itobcm = ntobc1-1 
    1074                    itobcp = ntobc2-1 
    1075                 ELSE 
    1076                    itobcm = ntobc1 
    1077                    itobcp = ntobc2 
    1078                 ENDIF 
    1079  
    1080              ELSEIF ( itimom == itimo .AND. itimop >= itimo ) THEN 
    1081                 IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN 
    1082                    itobcm = ntobc1 
    1083                    itobcp = ntobc2 
    1084                 ELSE 
    1085                    itobcm = ntobc1+1 
    1086                    itobcp = ntobc2+1 
    1087                 ENDIF 
    1088  
    1089              ELSEIF ( itimom == itimo-1 .AND. itimop == itimo+1 ) THEN 
    1090                 IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN 
    1091                    itobcm = ntobc1-1 
    1092                    itobcp = ntobc2-1 
    1093                 ELSEIF (  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN 
    1094                    itobcm = ntobc1 
    1095                    itobcp = ntobc2 
    1096                 ELSEIF (  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) == itimop ) THEN 
    1097                    itobcm = ntobc1+1 
    1098                    itobcp = ntobc2+2 
    1099                 ELSE 
    1100                    IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 1?' 
    1101                 ENDIF 
     1004             ENDIF 
     1005 
     1006          ELSEIF ( itimom == itimo .AND. itimop >= itimo ) THEN 
     1007             IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN 
     1008                itobcm = ntobc1 
     1009                itobcp = ntobc2 
    11021010             ELSE 
    1103                 IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 2?' 
     1011                itobcm = ntobc1+1 
     1012                itobcp = ntobc2+1 
    11041013             ENDIF 
    11051014 
    1106           ENDIF 
    1107  
    1108        ELSE IF( lk_dynspg_exp ) THEN 
    1109           isrel=kt*rdt 
    1110           itobcm = ntobc1 
    1111           itobcp = ntobc2 
    1112        ENDIF 
    1113  
    1114        IF( ntobc == 1 .OR. nobc_dta == 0 ) THEN 
    1115           zxy = 0.e0 
    1116        ELSE IF( ntobc == 12 ) THEN 
    1117           zxy = FLOAT( nday + 15 - 30 * i15 ) / 30. 
     1015          ELSEIF ( itimom == itimo-1 .AND. itimop == itimo+1 ) THEN 
     1016             IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN 
     1017                itobcm = ntobc1-1 
     1018                itobcp = ntobc2-1 
     1019             ELSEIF (  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN 
     1020                itobcm = ntobc1 
     1021                itobcp = ntobc2 
     1022             ELSEIF (  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) == itimop ) THEN 
     1023                itobcm = ntobc1+1 
     1024                itobcp = ntobc2+2 
     1025             ELSE 
     1026                IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 1?' 
     1027             ENDIF 
     1028          ELSE 
     1029             IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 2?' 
     1030          ENDIF 
     1031 
     1032       ENDIF 
     1033 
     1034    ELSE IF( lk_dynspg_exp ) THEN 
     1035       isrel=kt*rdt 
     1036       itobcm = ntobc1 
     1037       itobcp = ntobc2 
     1038    ENDIF 
     1039 
     1040    IF( ntobc == 1 .OR. nobc_dta == 0 ) THEN 
     1041       zxy = 0.e0 
     1042    ELSE IF( ntobc == 12 ) THEN 
     1043       zxy = FLOAT( nday + 15 - 30 * i15 ) / 30. 
     1044    ELSE 
     1045       zxy = (tcobc(itobcm)-FLOAT(isrel)) / (tcobc(itobcm)-tcobc(itobcp)) 
     1046    ENDIF 
     1047 
     1048    IF( lp_obc_east ) THEN           !  fills sshfoe, ubtfoe (local to each processor) 
     1049       DO jj = nje0p1, nje1m1 
     1050          ij = jj -1 + njmpp 
     1051          sshfoe(jj) = ( zxy * sshedta(ij,2) + (1.-zxy) * sshedta(ij,1) ) * temsk(jj,1) 
     1052          ubtfoe(jj) = ( zxy * ubtedta(ij,2) + (1.-zxy) * ubtedta(ij,1) ) * uemsk(jj,1) 
     1053       END DO 
     1054    ENDIF 
     1055 
     1056    IF( lp_obc_west) THEN            !  fills sshfow, ubtfow (local to each processor) 
     1057       DO jj = njw0p1, njw1m1 
     1058          ij = jj -1 + njmpp 
     1059          sshfow(jj) = ( zxy * sshwdta(ij,2) + (1.-zxy) * sshwdta(ij,1) ) * twmsk(jj,1) 
     1060          ubtfow(jj) = ( zxy * ubtwdta(ij,2) + (1.-zxy) * ubtwdta(ij,1) ) * uwmsk(jj,1) 
     1061       END DO 
     1062    ENDIF 
     1063 
     1064    IF( lp_obc_north) THEN           !  fills sshfon, vbtfon (local to each processor) 
     1065       DO ji = nin0p1, nin1m1 
     1066          ii = ji -1 + nimpp 
     1067          sshfon(ji) = ( zxy * sshndta(ii,2) + (1.-zxy) * sshndta(ii,1) ) * tnmsk(ji,1) 
     1068          vbtfon(ji) = ( zxy * vbtndta(ii,2) + (1.-zxy) * vbtndta(ii,1) ) * vnmsk(ji,1) 
     1069       END DO 
     1070    ENDIF 
     1071 
     1072    IF( lp_obc_south) THEN           !  fills sshfos, vbtfos (local to each processor) 
     1073       DO ji = nis0p1, nis1m1 
     1074          ii = ji -1 + nimpp 
     1075          sshfos(ji) = ( zxy * sshsdta(ii,2) + (1.-zxy) * sshsdta(ii,1) ) * tsmsk(ji,1) 
     1076          vbtfos(ji) = ( zxy * vbtsdta(ii,2) + (1.-zxy) * vbtsdta(ii,1) ) * vsmsk(ji,1) 
     1077       END DO 
     1078    ENDIF 
     1079 
     1080  END SUBROUTINE obc_dta_bt 
     1081 
     1082#else 
     1083  !!----------------------------------------------------------------------------- 
     1084  !!   Default option 
     1085  !!----------------------------------------------------------------------------- 
     1086  SUBROUTINE obc_dta_bt ( kt, kbt )       ! Empty routine 
     1087    !! * Arguments 
     1088    INTEGER,INTENT(in) :: kt 
     1089    INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index 
     1090    WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt 
     1091    WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kbt 
     1092  END SUBROUTINE obc_dta_bt 
     1093#endif 
     1094 
     1095 
     1096  !!============================================================================== 
     1097  SUBROUTINE obc_read (kt, nt_x, ntobc_x, iyy, imm) 
     1098     !!------------------------------------------------------------------------- 
     1099     !!                      ***  ROUTINE obc_read  *** 
     1100     !! 
     1101     !! ** Purpose :  Read the boundary data in files identified by iyy and imm 
     1102     !!               According to the validated open boundaries, return the  
     1103     !!               following arrays : 
     1104     !!                sedta, tedta : East OBC salinity and temperature 
     1105     !!                uedta, vedta :   "   "  u and v velocity component       
     1106     !! 
     1107     !!                swdta, twdta : West OBC salinity and temperature 
     1108     !!                uwdta, vwdta :   "   "  u and v velocity component       
     1109     !! 
     1110     !!                sndta, tndta : North OBC salinity and temperature 
     1111     !!                undta, vndta :   "   "  u and v velocity component       
     1112     !! 
     1113     !!                ssdta, tsdta : South OBC salinity and temperature 
     1114     !!                usdta, vsdta :   "   "  u and v velocity component       
     1115     !! 
     1116     !! ** Method  :  These fields are read in the record ntobc_x of the files. 
     1117     !!               The number of records is already known. If  ntobc_x is greater 
     1118     !!               than the number of record, this routine will look for next file, 
     1119     !!               updating the indices (case of inter-annual obcs) or loop at the 
     1120     !!               begining in case of climatological file (ln_obc_clim = true ). 
     1121     !! ------------------------------------------------------------------------- 
     1122     !! History:     !  2005  ( P. Mathiot, C. Langlais ) Original code 
     1123     !!              !  2008  ( J,M, Molines ) Use IOM and cleaning 
     1124     !!-------------------------------------------------------------------------- 
     1125 
     1126    ! * Arguments 
     1127    INTEGER, INTENT( in ) :: kt, nt_x 
     1128    INTEGER, INTENT( inout ) :: ntobc_x , iyy, imm      ! yes ! inout ! 
     1129 
     1130    ! * Local variables 
     1131    CHARACTER (len=40) :: &    ! file names 
     1132         cl_obc_eTS   , cl_obc_eU,  cl_obc_eV,& 
     1133         cl_obc_wTS   , cl_obc_wU,  cl_obc_wV,& 
     1134         cl_obc_nTS   , cl_obc_nU,  cl_obc_nV,& 
     1135         cl_obc_sTS   , cl_obc_sU,  cl_obc_sV 
     1136 
     1137    INTEGER :: ikprint 
     1138    REAL(wp) :: zmin, zmax   ! control of boundary values 
     1139 
     1140    !IOM stuff 
     1141    INTEGER :: id_e, id_w, id_n, id_s, ji, jj 
     1142    INTEGER, DIMENSION(2) :: istart, icount 
     1143    
     1144    !-------------------------------------------------------------------------- 
     1145    IF ( ntobc_x > itobc ) THEN 
     1146      IF (ln_obc_clim) THEN  ! just loop on the same file 
     1147        ntobc_x = 1  
     1148      ELSE 
     1149        ! need to change file : it is always for an 'after' data 
     1150        IF ( cffile == 'annual' ) THEN ! go to next year file 
     1151          iyy = iyy + 1 
     1152        ELSE IF ( cffile =='monthly' ) THEN  ! go to next month file 
     1153          imm = imm + 1  
     1154          IF ( imm == 13 ) THEN  
     1155            imm = 1 ; iyy = iyy + 1 
     1156          ENDIF 
     1157        ELSE 
     1158         ctmp1='obcread : this type of obc file is not supported :( ' 
     1159         ctmp2=TRIM(cffile) 
     1160         CALL ctl_stop (ctmp1, ctmp2) 
     1161         ! cffile should be either annual or monthly ... 
     1162        ENDIF 
     1163       ! as the file is changed, need to update itobc etc ... 
     1164        CALL obc_dta_chktime (iyy,imm) 
     1165        ntobc_x = nrecbef() + 1 ! remember : this case occur for an after data 
     1166      ENDIF 
     1167    ENDIF 
     1168 
     1169    IF ( lp_obc_east ) THEN  
     1170       ! ... Read datafile and set temperature, salinity and normal velocity 
     1171       ! ... initialise the sedta, tedta, uedta arrays 
     1172       WRITE(cl_obc_eTS ,'("obc_east_TS_y"  ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1173       WRITE(cl_obc_eU  ,'("obc_east_U_y"   ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1174       WRITE(cl_obc_eV  ,'("obc_east_V_y"   ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1175       ! JMM this may change depending on the obc data format ... 
     1176       istart(:)=(/nje0+njmpp-1,1/) ; icount(:)=(/nje1-nje0 +1,jpk/) 
     1177       IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_eTS) 
     1178       IF (nje1 >= nje0 ) THEN 
     1179          CALL iom_open ( cl_obc_eTS , id_e ) 
     1180          CALL iom_get ( id_e, jpdom_unknown, 'votemper', tedta(nje0:nje1,:,nt_x), & 
     1181               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1182          CALL iom_get ( id_e, jpdom_unknown, 'vosaline', sedta(nje0:nje1,:,nt_x), & 
     1183               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1184          CALL iom_close (id_e) 
     1185          ! 
     1186          CALL iom_open ( cl_obc_eU , id_e ) 
     1187          CALL iom_get  ( id_e, jpdom_unknown, 'vozocrtx', uedta(nje0:nje1,:,nt_x), & 
     1188               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1189          CALL iom_close ( id_e ) 
     1190          ! 
     1191          CALL iom_open ( cl_obc_eV , id_e ) 
     1192          CALL iom_get ( id_e, jpdom_unknown, 'vomecrty', vedta(nje0:nje1,:,nt_x), & 
     1193               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1194          CALL iom_close ( id_e ) 
     1195 
     1196          ! mask the boundary values 
     1197          tedta(:,:,nt_x) = tedta(:,:,nt_x)*temsk(:,:) ;  sedta(:,:,nt_x) = sedta(:,:,nt_x)*temsk(:,:) 
     1198          uedta(:,:,nt_x) = uedta(:,:,nt_x)*uemsk(:,:) ;  vedta(:,:,nt_x) = vedta(:,:,nt_x)*vemsk(:,:) 
     1199 
     1200          ! check any outliers  
     1201          zmin=MINVAL( sedta(:,:,nt_x), mask=ltemsk ) ; zmax=MAXVAL(sedta(:,:,nt_x), mask=ltemsk) 
     1202          IF (  zmin < 5 .OR. zmax > 50)   THEN 
     1203             CALL ctl_stop('Error in sedta',' routine obcdta') 
     1204          ENDIF 
     1205          zmin=MINVAL( tedta(:,:,nt_x), mask=ltemsk ) ; zmax=MAXVAL(tedta(:,:,nt_x), mask=ltemsk) 
     1206          IF (  zmin < -10. .OR. zmax > 40)   THEN 
     1207             CALL ctl_stop('Error in tedta',' routine obcdta') 
     1208          ENDIF 
     1209          zmin=MINVAL( uedta(:,:,nt_x), mask=luemsk ) ; zmax=MAXVAL(uedta(:,:,nt_x), mask=luemsk) 
     1210          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1211             CALL ctl_stop('Error in uedta',' routine obcdta') 
     1212          ENDIF 
     1213          zmin=MINVAL( vedta(:,:,nt_x), mask=lvemsk ) ; zmax=MAXVAL(vedta(:,:,nt_x), mask=lvemsk) 
     1214          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1215             CALL ctl_stop('Error in vedta',' routine obcdta') 
     1216          ENDIF 
     1217 
     1218          !               Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1       
     1219          IF ( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     1220             WRITE(numout,*) 
     1221             WRITE(numout,*) ' Read East OBC data records ', ntobc_x 
     1222             ikprint = jpj/20 +1 
     1223             WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
     1224             CALL prihre( tedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1225             WRITE(numout,*) 
     1226             WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
     1227             CALL prihre( sedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1228             WRITE(numout,*) 
     1229             WRITE(numout,*) ' Normal velocity U  record 1  - printout every 3 level' 
     1230             CALL prihre( uedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1231             WRITE(numout,*) 
     1232             WRITE(numout,*) ' Tangential velocity V  record 1  - printout every 3 level' 
     1233             CALL prihre( vedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1234          ENDIF 
     1235       ENDIF 
     1236    ENDIF 
     1237!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1238    IF ( lp_obc_west ) THEN 
     1239       ! ... Read datafile and set temperature, salinity and normal velocity 
     1240       ! ... initialise the swdta, twdta, uwdta arrays 
     1241       WRITE(cl_obc_wTS ,'("obc_west_TS_y"  ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1242       WRITE(cl_obc_wU  ,'("obc_west_U_y"   ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1243       WRITE(cl_obc_wV  ,'("obc_west_V_y"   ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1244       istart(:)=(/njw0+njmpp-1,1/) ; icount(:)=(/njw1-njw0 +1,jpk/) 
     1245       IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_wTS) 
     1246 
     1247       IF ( njw1 >= njw0 ) THEN 
     1248          CALL iom_open ( cl_obc_wTS , id_w ) 
     1249          CALL iom_get ( id_w, jpdom_unknown, 'votemper', twdta(njw0:njw1,:,nt_x), &  
     1250               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1251 
     1252          CALL iom_get ( id_w, jpdom_unknown, 'vosaline', swdta(njw0:njw1,:,nt_x), & 
     1253               &               ktime=ntobc_x , kstart=istart, kcount= icount) 
     1254          CALL iom_close (id_w) 
     1255          ! 
     1256          CALL iom_open ( cl_obc_wU , id_w ) 
     1257          CALL iom_get  ( id_w, jpdom_unknown, 'vozocrtx', uwdta(njw0:njw1,:,nt_x),& 
     1258               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1259          CALL iom_close ( id_w ) 
     1260          ! 
     1261          CALL iom_open ( cl_obc_wV , id_w ) 
     1262          CALL iom_get ( id_w, jpdom_unknown, 'vomecrty', vwdta(njw0:njw1,:,nt_x), & 
     1263               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1264          CALL iom_close ( id_w ) 
     1265 
     1266          ! mask the boundary values 
     1267          twdta(:,:,nt_x) = twdta(:,:,nt_x)*twmsk(:,:) ;  swdta(:,:,nt_x) = swdta(:,:,nt_x)*twmsk(:,:) 
     1268          uwdta(:,:,nt_x) = uwdta(:,:,nt_x)*uwmsk(:,:) ;  vwdta(:,:,nt_x) = vwdta(:,:,nt_x)*vwmsk(:,:) 
     1269 
     1270          ! check any outliers 
     1271          zmin=MINVAL( swdta(:,:,nt_x), mask=ltwmsk ) ; zmax=MAXVAL(swdta(:,:,nt_x), mask=ltwmsk) 
     1272          IF (  zmin < 5 .OR. zmax > 50)   THEN 
     1273             CALL ctl_stop('Error in swdta',' routine obcdta') 
     1274          ENDIF 
     1275          zmin=MINVAL( twdta(:,:,nt_x), mask=ltwmsk ) ; zmax=MAXVAL(twdta(:,:,nt_x), mask=ltwmsk) 
     1276          IF (  zmin < -10. .OR. zmax > 40)   THEN 
     1277             CALL ctl_stop('Error in twdta',' routine obcdta') 
     1278          ENDIF 
     1279          zmin=MINVAL( uwdta(:,:,nt_x), mask=luwmsk ) ; zmax=MAXVAL(uwdta(:,:,nt_x), mask=luwmsk) 
     1280          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1281             CALL ctl_stop('Error in uwdta',' routine obcdta') 
     1282          ENDIF 
     1283          zmin=MINVAL( vwdta(:,:,nt_x), mask=lvwmsk ) ; zmax=MAXVAL(vwdta(:,:,nt_x), mask=lvwmsk) 
     1284          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1285             CALL ctl_stop('Error in vwdta',' routine obcdta') 
     1286          ENDIF 
     1287 
     1288 
     1289          IF ( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     1290             WRITE(numout,*) 
     1291             WRITE(numout,*) ' Read West OBC data records ', ntobc_x 
     1292             ikprint = jpj/20 +1 
     1293             WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
     1294             CALL prihre( twdta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1295             WRITE(numout,*) 
     1296             WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
     1297             CALL prihre( swdta(:,:,nt_x),jpj,jpk, 1, jpj, ikprint,   jpk, 1, -3, 1., numout ) 
     1298             WRITE(numout,*) 
     1299             WRITE(numout,*) ' Normal velocity U  record 1  - printout every 3 level' 
     1300             CALL prihre( uwdta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1301             WRITE(numout,*) 
     1302             WRITE(numout,*) ' Tangential velocity V  record 1  - printout every 3 level' 
     1303             CALL prihre( vwdta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout ) 
     1304          ENDIF 
     1305       END IF 
     1306    ENDIF 
     1307!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1308    IF( lp_obc_north) THEN 
     1309       WRITE(cl_obc_nTS ,'("obc_north_TS_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1310       WRITE(cl_obc_nV  ,'("obc_north_V_y"  ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1311       WRITE(cl_obc_nU  ,'("obc_north_U_y"  ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1312       istart(:)=(/nin0+nimpp-1,1/) ; icount(:)=(/nin1-nin0 +1,jpk/) 
     1313       IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_nTS) 
     1314       IF ( nin1 >= nin0 ) THEN 
     1315          CALL iom_open ( cl_obc_nTS , id_n ) 
     1316          CALL iom_get ( id_n, jpdom_unknown, 'votemper', tndta(nin0:nin1,:,nt_x), & 
     1317               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1318          CALL iom_get ( id_n, jpdom_unknown, 'vosaline', sndta(nin0:nin1,:,nt_x), & 
     1319               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1320          CALL iom_close (id_n) 
     1321          ! 
     1322          CALL iom_open ( cl_obc_nU , id_n ) 
     1323          CALL iom_get  ( id_n, jpdom_unknown, 'vozocrtx', undta(nin0:nin1,:,nt_x), & 
     1324               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1325          CALL iom_close ( id_n ) 
     1326          ! 
     1327          CALL iom_open ( cl_obc_nV , id_n ) 
     1328          CALL iom_get  ( id_n, jpdom_unknown, 'vomecrty', vndta(nin0:nin1,:,nt_x), & 
     1329               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1330          CALL iom_close ( id_n ) 
     1331 
     1332          ! mask the boundary values 
     1333          tndta(:,:,nt_x) = tndta(:,:,nt_x)*tnmsk(:,:) ;  sndta(:,:,nt_x) = sndta(:,:,nt_x)*tnmsk(:,:) 
     1334          undta(:,:,nt_x) = undta(:,:,nt_x)*unmsk(:,:) ;  vndta(:,:,nt_x) = vndta(:,:,nt_x)*vnmsk(:,:) 
     1335 
     1336          ! check any outliers 
     1337          zmin=MINVAL( sndta(:,:,nt_x), mask=ltnmsk ) ; zmax=MAXVAL(sndta(:,:,nt_x), mask=ltnmsk) 
     1338          IF (  zmin < 5 .OR. zmax > 50)   THEN 
     1339             CALL ctl_stop('Error in sndta',' routine obcdta') 
     1340          ENDIF 
     1341          zmin=MINVAL( tndta(:,:,nt_x), mask=ltnmsk ) ; zmax=MAXVAL(tndta(:,:,nt_x), mask=ltnmsk) 
     1342          IF (  zmin < -10. .OR. zmax > 40)   THEN 
     1343             CALL ctl_stop('Error in tndta',' routine obcdta') 
     1344          ENDIF 
     1345          zmin=MINVAL( undta(:,:,nt_x), mask=lunmsk ) ; zmax=MAXVAL(undta(:,:,nt_x), mask=lunmsk) 
     1346          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1347             CALL ctl_stop('Error in undta',' routine obcdta') 
     1348          ENDIF 
     1349          zmin=MINVAL( vndta(:,:,nt_x), mask=lvnmsk ) ; zmax=MAXVAL(vndta(:,:,nt_x), mask=lvnmsk) 
     1350          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1351             CALL ctl_stop('Error in vndta',' routine obcdta') 
     1352          ENDIF 
     1353 
     1354          IF ( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     1355             WRITE(numout,*) 
     1356             WRITE(numout,*) ' Read North OBC data records ', ntobc_x 
     1357             ikprint = jpi/20 +1 
     1358             WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
     1359             CALL prihre( tndta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1360             WRITE(numout,*) 
     1361             WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
     1362             CALL prihre( sndta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1363             WRITE(numout,*) 
     1364             WRITE(numout,*) ' Normal velocity V  record 1  - printout every 3 level' 
     1365             CALL prihre( vndta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1366             WRITE(numout,*) 
     1367             WRITE(numout,*) ' Tangential  velocity U  record 1  - printout every 3 level' 
     1368             CALL prihre( undta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1369          ENDIF 
     1370       ENDIF 
     1371    ENDIF 
     1372!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1373    IF( lp_obc_south) THEN  
     1374       WRITE(cl_obc_sTS ,'("obc_south_TS_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1375       WRITE(cl_obc_sV  ,'("obc_south_V_y"  ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1376       WRITE(cl_obc_sU  ,'("obc_south_U_y"  ,i4.4,"m",i2.2,".nc")' ) iyy,imm 
     1377       istart(:)=(/nis0+nimpp-1,1/) ; icount(:)=(/nis1-nis0 +1,jpk/) 
     1378       IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_sTS) 
     1379       IF ( nis1 >= nis0 ) THEN  
     1380          CALL iom_open ( cl_obc_sTS , id_s ) 
     1381          CALL iom_get ( id_s, jpdom_unknown, 'votemper', tsdta(nis0:nis1,:,nt_x), & 
     1382               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1383          CALL iom_get ( id_s, jpdom_unknown, 'vosaline', ssdta(nis0:nis1,:,nt_x), & 
     1384               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1385          CALL iom_close (id_s) 
     1386          ! 
     1387          CALL iom_open ( cl_obc_sU , id_s ) 
     1388          CALL iom_get  ( id_s, jpdom_unknown, 'vozocrtx', usdta(nis0:nis1,:,nt_x), & 
     1389               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1390          CALL iom_close ( id_s ) 
     1391          ! 
     1392          CALL iom_open ( cl_obc_sV , id_s ) 
     1393          CALL iom_get  ( id_s, jpdom_unknown, 'vomecrty', vsdta(nis0:nis1,:,nt_x), & 
     1394               &               ktime=ntobc_x , kstart=istart, kcount= icount ) 
     1395          CALL iom_close ( id_s ) 
     1396 
     1397          ! mask the boundary values 
     1398          tsdta(:,:,nt_x) = tsdta(:,:,nt_x)*tsmsk(:,:) ;  ssdta(:,:,nt_x) = ssdta(:,:,nt_x)*tsmsk(:,:) 
     1399          usdta(:,:,nt_x) = usdta(:,:,nt_x)*usmsk(:,:) ;  vsdta(:,:,nt_x) = vsdta(:,:,nt_x)*vsmsk(:,:) 
     1400 
     1401          ! check any outliers 
     1402          zmin=MINVAL( ssdta(:,:,nt_x), mask=ltsmsk ) ; zmax=MAXVAL(ssdta(:,:,nt_x), mask=ltsmsk) 
     1403          IF (  zmin < 5 .OR. zmax > 50)   THEN 
     1404             CALL ctl_stop('Error in ssdta',' routine obcdta') 
     1405          ENDIF 
     1406          zmin=MINVAL( tsdta(:,:,nt_x), mask=ltsmsk ) ; zmax=MAXVAL(tsdta(:,:,nt_x), mask=ltsmsk) 
     1407          IF (  zmin < -10. .OR. zmax > 40)   THEN 
     1408             CALL ctl_stop('Error in tsdta',' routine obcdta') 
     1409          ENDIF 
     1410          zmin=MINVAL( usdta(:,:,nt_x), mask=lusmsk ) ; zmax=MAXVAL(usdta(:,:,nt_x), mask=lusmsk) 
     1411          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1412             CALL ctl_stop('Error in usdta',' routine obcdta') 
     1413          ENDIF 
     1414          zmin=MINVAL( vsdta(:,:,nt_x), mask=lvsmsk ) ; zmax=MAXVAL(vsdta(:,:,nt_x), mask=lvsmsk) 
     1415          IF (  zmin < -5. .OR. zmax > 5.)   THEN 
     1416             CALL ctl_stop('Error in vsdta',' routine obcdta') 
     1417          ENDIF 
     1418 
     1419          IF ( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN 
     1420             WRITE(numout,*) 
     1421             WRITE(numout,*) ' Read South OBC data records ', ntobc_x 
     1422             ikprint = jpi/20 +1 
     1423             WRITE(numout,*) ' Temperature  record 1 - printout every 3 level' 
     1424             CALL prihre( tsdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1425             WRITE(numout,*) 
     1426             WRITE(numout,*) ' Salinity  record 1 - printout every 3 level' 
     1427             CALL prihre( ssdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1428             WRITE(numout,*) 
     1429             WRITE(numout,*) ' Normal velocity V  record 1  - printout every 3 level' 
     1430             CALL prihre( vsdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1431             WRITE(numout,*) 
     1432             WRITE(numout,*) ' Tangential velocity U  record 1  - printout every 3 level' 
     1433             CALL prihre( usdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout ) 
     1434          ENDIF 
     1435       ENDIF 
     1436    ENDIF 
     1437!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1438  END SUBROUTINE obc_read  
     1439 
     1440  INTEGER FUNCTION nrecbef() 
     1441      !!----------------------------------------------------------------------- 
     1442      !!                     ***    FUNCTION nrecbef   *** 
     1443      !! 
     1444      !!  Purpose : - provide the before record number in files, with respect to zjcnes 
     1445      !! 
     1446      !!    History : 2008-04 : ( J.M. Molines ) Original code 
     1447      !!----------------------------------------------------------------------- 
     1448 
     1449      INTEGER :: it , idum 
     1450 
     1451    idum = itobc 
     1452    DO it =1, itobc 
     1453       IF ( ztcobc(it) > zjcnes ) THEN ;  idum = it - 1 ; EXIT ;  ENDIF 
     1454    ENDDO 
     1455    ! idum can be 0 (climato, before first record) 
     1456    IF ( idum == 0 ) THEN 
     1457       IF ( ln_obc_clim ) THEN 
     1458         idum = itobc 
    11181459       ELSE 
    1119           zxy = (tcobc(itobcm)-FLOAT(isrel)) / (tcobc(itobcm)-tcobc(itobcp)) 
    1120        ENDIF 
    1121  
    1122        IF( lp_obc_east ) THEN           !  fills sshfoe, ubtfoe (local to each processor) 
    1123           DO jj = nje0p1, nje1m1 
    1124              ij = jj -1 + njmpp 
    1125              sshfoe(jj) = ( zxy * sshedta(ij,2) + (1.-zxy) * sshedta(ij,1) ) * temsk(jj,1) 
    1126              ubtfoe(jj) = ( zxy * ubtedta(ij,2) + (1.-zxy) * ubtedta(ij,1) ) * uemsk(jj,1) 
    1127           END DO 
    1128        ENDIF 
    1129  
    1130        IF( lp_obc_west) THEN            !  fills sshfow, ubtfow (local to each processor) 
    1131           DO jj = njw0p1, njw1m1 
    1132              ij = jj -1 + njmpp 
    1133              sshfow(jj) = ( zxy * sshwdta(ij,2) + (1.-zxy) * sshwdta(ij,1) ) * twmsk(jj,1) 
    1134              ubtfow(jj) = ( zxy * ubtwdta(ij,2) + (1.-zxy) * ubtwdta(ij,1) ) * uwmsk(jj,1) 
    1135           END DO 
    1136        ENDIF 
    1137  
    1138        IF( lp_obc_north) THEN           !  fills sshfon, vbtfon (local to each processor) 
    1139           DO ji = nin0p1, nin1m1 
    1140              ii = ji -1 + nimpp 
    1141              sshfon(ji) = ( zxy * sshndta(ii,2) + (1.-zxy) * sshndta(ii,1) ) * tnmsk(ji,1) 
    1142              vbtfon(ji) = ( zxy * vbtndta(ii,2) + (1.-zxy) * vbtndta(ii,1) ) * vnmsk(ji,1) 
    1143           END DO 
    1144        ENDIF 
    1145  
    1146        IF( lp_obc_south) THEN           !  fills sshfos, vbtfos (local to each processor) 
    1147           DO ji = nis0p1, nis1m1 
    1148              ii = ji -1 + nimpp 
    1149              sshfos(ji) = ( zxy * sshsdta(ii,2) + (1.-zxy) * sshsdta(ii,1) ) * tsmsk(ji,1) 
    1150              vbtfos(ji) = ( zxy * vbtsdta(ii,2) + (1.-zxy) * vbtsdta(ii,1) ) * vsmsk(ji,1) 
    1151           END DO 
    1152        ENDIF 
    1153  
    1154    END SUBROUTINE obc_dta_bt 
     1460         ctmp1='obc_dta: find ntobc == 0 for  non climatological file ' 
     1461         ctmp2='consider adding a first record in your data file ' 
     1462         CALL ctl_stop(ctmp1, ctmp2) 
     1463       ENDIF 
     1464    ENDIF 
     1465    ! idum can be itobc ( zjcnes > ztcobc (itobc) ) 
     1466    !  This is not a problem ... 
     1467    nrecbef = idum 
     1468 
     1469  END FUNCTION nrecbef 
    11551470 
    11561471#else 
    1157    !!----------------------------------------------------------------------------- 
    1158    !!   Default option 
    1159    !!----------------------------------------------------------------------------- 
    1160    SUBROUTINE obc_dta_bt ( kt, kbt )       ! Empty routine 
    1161       !! * Arguments 
    1162       INTEGER,INTENT(in) :: kt 
    1163       INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index 
    1164       WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt 
    1165       WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kbt 
    1166    END SUBROUTINE obc_dta_bt 
     1472  !!------------------------------------------------------------------------------ 
     1473  !!   default option:           Dummy module          NO Open Boundary Conditions 
     1474  !!------------------------------------------------------------------------------ 
     1475CONTAINS 
     1476  SUBROUTINE obc_dta( kt )             ! Dummy routine 
     1477    INTEGER, INTENT (in) :: kt 
     1478    WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt 
     1479  END SUBROUTINE obc_dta 
    11671480#endif 
    1168  
    1169 #else 
    1170    !!-------------------------------------------------------------------- 
    1171    !!  default option  :  Dummy module    NO Open Boundary Conditions 
    1172    !!-------------------------------------------------------------------- 
    1173 CONTAINS 
    1174    SUBROUTINE obc_dta( kt )             ! Dummy routine 
    1175      INTEGER, INTENT (in) :: kt 
    1176      WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt 
    1177    END SUBROUTINE obc_dta 
    1178    SUBROUTINE obc_dta_bt( kt, jn)             ! Dummy routine 
    1179      INTEGER, INTENT (in) :: kt, jn 
    1180      WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt 
    1181      WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', jn 
    1182    END SUBROUTINE obc_dta_bt 
    1183 #endif 
    1184  
    1185    !!===================================================================== 
    11861481END MODULE obcdta 
  • trunk/NEMO/OPA_SRC/OBC/obcdyn_bt.F90

    r1146 r1151  
    11MODULE obcdyn_bt 
    2 #if ( defined key_dynspg_ts || defined key_dynspg_exp ) & defined key_obc 
     2#if ( defined key_dynspg_ts || defined key_dynspg_exp ) && defined key_obc 
    33   !!================================================================================= 
    44   !!                       ***  MODULE  obcdyn_bt  *** 
  • trunk/NEMO/OPA_SRC/OBC/obcini.F90

    r1146 r1151  
    1818   USE lib_mpp         ! for mpp_sum 
    1919   USE in_out_manager  ! I/O units 
     20   USE dynspg_oce      ! flag lk_dynspg_flt 
    2021 
    2122   IMPLICIT NONE 
     
    2930   !!--------------------------------------------------------------------------------- 
    3031   !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    31    !! $Header$  
     32   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/OBC/obcini.F90,v 1.9 2006/05/11 15:15:31 opalod Exp $  
    3233   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3334   !!--------------------------------------------------------------------------------- 
     
    5859      !!---------------------------------------------------------------------- 
    5960      !! * Modules used 
    60       USE obcrst,   ONLY :   obc_rst_lec   ! Make obc_rst_lec routine available 
     61      USE obcrst,   ONLY :   obc_rst_read   ! Make obc_rst_read routine available 
    6162      USE obcdom,   ONLY :   obc_dom       ! Make obc_dom routine available 
    6263 
     
    7071         &             rdpeob, rdpwob, rdpnob, rdpsob,   & 
    7172         &             zbsic1, zbsic2, zbsic3,           & 
    72          &             nbic, volemp, nobc_dta,           & 
     73         &             nbic, volemp, nobc_dta, cffile,   & 
    7374         &             ln_obc_clim, ln_vol_cst, ln_obc_fla 
    7475      !!---------------------------------------------------------------------- 
     
    9697      ! By security we set rdpxin and rdpxob respectively 
    9798      ! to 1. and 15. if the corresponding OBC is not activated 
    98       IF( .NOT.lp_obc_east ) THEN 
    99          rdpein = 1. 
    100          rdpeob = 15. 
    101       END IF 
    102       IF( .NOT.lp_obc_west ) THEN 
    103          rdpwin = 1. 
    104          rdpwob = 15. 
    105       END IF 
    106       IF( .NOT.lp_obc_north ) THEN 
    107          rdpnin = 1. 
    108          rdpnob = 15. 
    109       END IF 
    110       IF( .NOT.lp_obc_south ) THEN 
    111          rdpsin = 1. 
    112          rdpsob = 15. 
    113       END IF 
     99      IF( .NOT.lp_obc_east  ) THEN ; rdpein = 1. ; rdpeob = 15. ; END IF 
     100      IF( .NOT.lp_obc_west  ) THEN ; rdpwin = 1. ; rdpwob = 15. ; END IF 
     101      IF( .NOT.lp_obc_north ) THEN ; rdpnin = 1. ; rdpnob = 15. ; END IF 
     102      IF( .NOT.lp_obc_south ) THEN ; rdpsin = 1. ; rdpsob = 15. ; END IF 
    114103 
    115104      ! number of open boudaries and open boundary indicators 
     
    131120      IF(lwp) WRITE(numout,*) '         initial state used (=0)             ' 
    132121      IF(lwp) WRITE(numout,*) '         climatology (true) or not:', ln_obc_clim 
     122      IF(lwp) WRITE(numout,*) '         vol_cst (true) or not:', ln_vol_cst 
     123      IF(lwp) THEN 
     124          IF ( lk_dynspg_flt ) WRITE(numout,*) '                   dynspg_flt T ==> vol_cst forced to T' 
     125      ENDIF 
    133126      IF(lwp) WRITE(numout,*) ' ' 
    134127      IF(lwp) WRITE(numout,*) '                                 WARNING                                                  ' 
     
    141134      IF(lwp) WRITE(numout,*) '         For the rigid-lid case or the filtered free surface case,                        ' 
    142135      IF(lwp) WRITE(numout,*) '         radiation, relaxation or presciption of data can be applied                      ' 
     136 
    143137      IF( lwp.AND.lp_obc_east ) THEN 
    144138         WRITE(numout,*) '         East open boundary :' 
     
    173167      ! ------------------------------ 
    174168 
    175       ! ... convert rdp$ob in seconds 
    176       rdpein = rdpein * rday 
    177       rdpwin = rdpwin * rday 
    178       rdpnin = rdpnin * rday 
    179       rdpsin = rdpsin * rday 
    180       rdpeob = rdpeob * rday 
    181       rdpwob = rdpwob * rday 
    182       rdpnob = rdpnob * rday 
    183       rdpsob = rdpsob * rday 
    184       lfbceast  = .FALSE. 
    185       lfbcwest  = .FALSE. 
    186       lfbcnorth = .FALSE. 
    187       lfbcsouth = .FALSE. 
     169      ! ...                          convert rdp$ob in seconds 
     170      ! Fixed Bdy flag              inbound                outbound 
     171      lfbceast  = .FALSE. ; rdpein = rdpein * rday  ; rdpeob = rdpeob * rday 
     172      lfbcwest  = .FALSE. ; rdpwin = rdpwin * rday  ; rdpwob = rdpwob * rday 
     173      lfbcnorth = .FALSE. ; rdpnin = rdpnin * rday  ; rdpnob = rdpnob * rday 
     174      lfbcsouth = .FALSE. ; rdpsin = rdpsin * rday  ; rdpsob = rdpsob * rday 
    188175      inumfbc = 0 
    189176      ! ... look for Fixed Boundaries (rdp = 0 ) 
     
    193180      IF( lp_obc_east )  THEN 
    194181         IF( (rdpein+rdpeob) == 0 )  THEN 
    195             lfbceast = .TRUE. 
    196             rdpein = 1e-3 
    197             rdpeob = 1e-3 
     182            lfbceast = .TRUE. ; rdpein = 1e-3 ; rdpeob = 1e-3 
    198183            inumfbc = inumfbc+1 
    199184         ELSEIF ( (rdpein*rdpeob) == 0 )  THEN 
     
    201186         END IF 
    202187      END IF 
     188 
    203189      IF( lp_obc_west )  THEN 
    204190         IF( (rdpwin + rdpwob) == 0 )  THEN 
    205             lfbcwest = .TRUE. 
    206             rdpwin = 1e-3 
    207             rdpwob = 1e-3 
     191            lfbcwest = .TRUE. ; rdpwin = 1e-3 ; rdpwob = 1e-3 
    208192            inumfbc = inumfbc+1 
    209193         ELSEIF ( (rdpwin*rdpwob) == 0 )  THEN 
     
    213197      IF( lp_obc_north )  THEN 
    214198         IF( (rdpnin + rdpnob) == 0 )  THEN 
    215             lfbcnorth = .TRUE. 
    216             rdpnin = 1e-3 
    217             rdpnob = 1e-3 
     199            lfbcnorth = .TRUE. ; rdpnin = 1e-3 ; rdpnob = 1e-3 
    218200            inumfbc = inumfbc+1 
    219201         ELSEIF ( (rdpnin*rdpnob) == 0 )  THEN 
     
    223205      IF( lp_obc_south )  THEN 
    224206         IF( (rdpsin + rdpsob) == 0 )  THEN 
    225             lfbcsouth = .TRUE. 
    226             rdpsin = 1e-3 
    227             rdpsob = 1e-3 
     207            lfbcsouth = .TRUE. ; rdpsin = 1e-3 ; rdpsob = 1e-3 
    228208            inumfbc = inumfbc+1 
    229209         ELSEIF ( (rdpsin*rdpsob) == 0 )  THEN 
     
    338318      IF( lp_obc_east ) THEN 
    339319         !... (jpjed,jpjefm1),jpieob 
    340          DO jj = nje0, nje1m1 
    341 # if defined key_dynspg_rl 
    342             DO ji = nie0, nie1 
    343 # else 
    344             DO ji = nie0p1, nie1p1 
    345 # endif 
    346                bmask(ji,jj) = 0.e0 
    347             END DO 
    348          END DO 
     320         IF (lk_dynspg_rl ) THEN ; bmask(nie0:nie1,nje0:nje1m1) = 0.e0 ; ELSE ; bmask(nie0p1:nie1p1,nje0:nje1m1) = 0.e0 ; ENDIF 
    349321 
    350322         ! ... initilization to zero 
    351          uemsk(:,:) = 0.e0 
    352          vemsk(:,:) = 0.e0 
    353          temsk(:,:) = 0.e0 
     323         uemsk(:,:) = 0.e0 ; vemsk(:,:) = 0.e0 ; temsk(:,:) = 0.e0 
    354324 
    355325         ! ... set 2D mask on East OBC,  Vopt 
    356326         DO ji = fs_nie0, fs_nie1 
    357327            DO jj = nje0, nje1 
    358                uemsk(jj,:) = umask(ji,  jj,:) 
    359                vemsk(jj,:) = vmask(ji+1,jj,:) 
    360                temsk(jj,:) = tmask(ji+1,jj,:) 
     328               uemsk(jj,:) = umask(ji,  jj,:) * tmask_i(ji,jj)   * tmask_i(ji+1,jj) 
     329               vemsk(jj,:) = vmask(ji+1,jj,:) * tmask_i(ji+1,jj)  
     330               temsk(jj,:) = tmask(ji+1,jj,:) * tmask_i(ji+1,jj)  
    361331            END DO 
    362332         END DO 
     
    366336      IF( lp_obc_west ) THEN 
    367337         ! ... (jpjwd,jpjwfm1),jpiwob 
    368          DO jj = njw0, njw1m1 
    369             DO ji = niw0, niw1 
    370                bmask(ji,jj) = 0.e0 
    371             END DO 
    372          END DO 
     338         bmask(niw0:niw1,njw0:njw1m1) = 0.e0 
    373339 
    374340         ! ... initilization to zero 
    375          uwmsk(:,:) = 0.e0 
    376          vwmsk(:,:) = 0.e0 
    377          twmsk(:,:) = 0.e0 
     341         uwmsk(:,:) = 0.e0 ; vwmsk(:,:) = 0.e0 ; twmsk(:,:) = 0.e0   
    378342 
    379343         ! ... set 2D mask on West OBC,  Vopt 
    380344         DO ji = fs_niw0, fs_niw1 
    381345            DO jj = njw0, njw1 
    382                uwmsk(jj,:) = umask(ji,jj,:) 
    383                vwmsk(jj,:) = vmask(ji,jj,:) 
    384                twmsk(jj,:) = tmask(ji,jj,:) 
     346               uwmsk(jj,:) = umask(ji,jj,:) * tmask_i(ji,jj)   * tmask_i(ji+1,jj) 
     347               vwmsk(jj,:) = vmask(ji,jj,:) * tmask_i(ji,jj)   
     348               twmsk(jj,:) = tmask(ji,jj,:) * tmask_i(ji,jj) 
    385349            END DO 
    386350         END DO 
     
    390354      IF( lp_obc_north ) THEN 
    391355         ! ... jpjnob,(jpind,jpisfm1) 
    392 # if defined key_dynspg_rl 
    393          DO jj = njn0, njn1 
    394 # else 
    395          DO jj = njn0p1, njn1p1 
    396 # endif 
    397             DO ji = nin0, nin1m1 
    398                bmask(ji,jj) = 0.e0 
    399             END DO 
    400          END DO 
     356         IF (lk_dynspg_rl ) THEN ; bmask(nin0:nin1m1,njn0:njn1) = 0.e0 ; ELSE ; bmask(nin0:nin1m1,njn0p1:njn1p1) = 0.e0 ; ENDIF 
    401357 
    402358         ! ... initilization to zero 
    403          unmsk(:,:) = 0.e0 
    404          vnmsk(:,:) = 0.e0 
    405          tnmsk(:,:) = 0.e0 
     359         unmsk(:,:) = 0.e0 ; vnmsk(:,:) = 0.e0 ; tnmsk(:,:) = 0.e0 
    406360 
    407361         ! ... set 2D mask on North OBC,  Vopt 
    408362         DO jj = fs_njn0, fs_njn1 
    409363            DO ji = nin0, nin1 
    410                unmsk(ji,:) = umask(ji,jj+1,:) 
    411                vnmsk(ji,:) = vmask(ji,jj  ,:) 
    412                tnmsk(ji,:) = tmask(ji,jj+1,:) 
     364               unmsk(ji,:) = umask(ji,jj+1,:) * tmask_i(ji,jj+1)  
     365               vnmsk(ji,:) = vmask(ji,jj  ,:) * tmask_i(ji,jj)   * tmask_i(ji,jj+1) 
     366               tnmsk(ji,:) = tmask(ji,jj+1,:) * tmask_i(ji,jj+1) 
    413367            END DO 
    414368         END DO 
     
    418372      IF( lp_obc_south ) THEN  
    419373         ! ... jpjsob,(jpisd,jpisfm1) 
    420          DO jj = njs0, njs1 
    421             DO ji = nis0, nis1m1 
    422                bmask(ji,jj) = 0.e0 
     374         bmask(nis0:nis1m1,njs0:njs1) = 0.e0 
     375 
     376         ! ... initilization to zero 
     377         usmsk(:,:) = 0.e0 ; vsmsk(:,:) = 0.e0 ; tsmsk(:,:) = 0.e0 
     378 
     379         ! ... set 2D mask on South OBC,  Vopt 
     380         DO jj = fs_njs0, fs_njs1  
     381            DO ji = nis0, nis1 
     382               usmsk(ji,:) = umask(ji,jj,:) * tmask_i(ji,jj)  
     383               vsmsk(ji,:) = vmask(ji,jj,:) * tmask_i(ji,jj) * tmask_i(ji,jj+1) 
     384               tsmsk(ji,:) = tmask(ji,jj,:) * tmask_i(ji,jj) 
    423385            END DO 
    424386         END DO 
    425387 
    426          ! ... initilization to zero 
    427          usmsk(:,:) = 0.e0 
    428          vsmsk(:,:) = 0.e0 
    429          tsmsk(:,:) = 0.e0 
    430  
    431          ! ... set 2D mask on South OBC,  Vopt 
    432          DO jj = njs0, njs1 
    433             DO ji = nis0, nis1 
    434                usmsk(ji,:) = umask(ji,jj,:) 
    435                vsmsk(ji,:) = vmask(ji,jj,:) 
    436                tsmsk(ji,:) = tmask(ji,jj,:) 
    437             END DO 
    438          END DO 
    439  
    440       END IF 
    441  
    442 # if defined key_dynspg_flt 
    443  
    444       ! ... Initialize obcumask and obcvmask for the Force filtering  
    445       !     boundary condition in dynspg_flt 
    446       obcumask(:,:) = umask(:,:,1) 
    447       obcvmask(:,:) = vmask(:,:,1) 
    448  
    449       ! ... Initialize obctmsk on overlap region and obcs. This mask 
    450       !     is used in obcvol.F90 to calculate cumulate flux E-P.  
    451       !     - no flux E-P on obcs and overlap region (jpereci = jprecj = 1) 
    452       obctmsk(:,:) = tmask(:,:,1)      
    453       obctmsk(1  ,:) = 0.e0            
    454       obctmsk(jpi,:) = 0.e0            
    455       obctmsk(:  ,1) = 0.e0            
    456       obctmsk(:,jpj) = 0.e0            
    457  
    458       IF( lp_obc_east ) THEN 
    459          ! ... East obc Force filtering mask for the grad D 
    460          DO ji = nie0, nie1 
    461             DO jj = nje0p1, nje1m1 
    462                obcumask(ji  ,jj)=0.e0 
    463                obcvmask(ji+1,jj)=0.e0 
    464             END DO 
    465          END DO 
    466  
    467          ! ... set to 0 on East OBC 
    468          DO jj = nje0p1, nje1m1 
    469             DO ji = nie0p1, nie1p1 
    470                obctmsk(ji,jj) = 0.e0 
    471             END DO 
    472          END DO 
    473       END IF 
    474  
    475       IF( lp_obc_west ) THEN 
    476          ! ... West obc Force filtering mask for the grad D 
    477          DO ji = niw0, niw1 
    478             DO jj = njw0p1, njw1m1 
    479                obcumask(ji,jj)=0.e0 
    480                obcvmask(ji,jj)=0.e0 
    481             END DO 
    482          END DO 
    483  
    484          ! ... set to 0 on West OBC 
    485          DO jj = njw0p1, njw1m1 
    486             DO ji = niw0, niw1 
    487                obctmsk(ji,jj) = 0.e0 
    488             END DO 
    489          END DO 
    490       END IF 
    491  
    492       IF( lp_obc_north ) THEN 
    493          ! ... North obc Force filtering mask for the grad D 
    494          DO jj = njn0, njn1 
    495             DO ji = nin0p1, nin1m1 
    496                obcvmask(ji,jj  )=0.e0 
    497                obcumask(ji,jj+1)=0.e0 
    498             END DO 
    499          END DO 
    500  
    501          ! ... set to 0 on North OBC 
    502          DO jj = njn0p1, njn1p1 
    503             DO ji = nin0p1, nin1m1 
    504                obctmsk(ji,jj) = 0.e0 
    505             END DO 
    506          END DO 
    507       END IF 
    508  
    509       IF( lp_obc_south ) THEN 
    510          ! ... South obc Force filtering mask for the grad D 
    511          DO jj = njs0, njs1  
    512             DO ji = nis0p1, nis1m1 
    513                obcumask(ji,jj)=0.e0 
    514                obcvmask(ji,jj)=0.e0 
    515             END DO 
    516          END DO 
    517  
    518          ! ... set to 0 on South OBC 
    519          DO jj = njs0, njs1 
    520             DO ji = nis0p1, nis1m1 
    521                obctmsk(ji,jj) = 0.e0 
    522             END DO 
    523          END DO 
    524       END IF 
    525  
    526 # endif 
    527  
    528 # if ! defined key_dynspg_rl 
    529  
    530       IF ( ln_vol_cst ) THEN 
    531  
    532          ! 3.1 Total lateral surface for each open boundary 
    533          ! ------------------------------------------------ 
    534  
    535          ! ... West open boundary surface 
    536          IF( lp_obc_west ) THEN 
    537             DO ji = niw0, niw1 
    538                DO jj = 1, jpj  
    539                   obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1) 
    540                END DO 
    541             END DO 
    542          END IF 
    543  
    544          ! ... East open boundary surface 
    545          IF( lp_obc_east ) THEN 
    546             DO ji = nie0, nie1 
    547                DO jj = 1, jpj  
    548                   obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1) 
    549                END DO 
    550             END DO 
    551          END IF 
    552  
    553          ! ... North open boundary vertical surface 
    554          IF( lp_obc_north ) THEN 
    555             DO jj = njn0, njn1 
    556                DO ji = 1, jpi 
    557                   obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vnmsk(ji,1) 
    558                END DO 
    559             END DO 
    560          END IF 
    561  
    562          ! ... South open boundary vertical surface 
    563          IF( lp_obc_south ) THEN 
    564             DO jj = njs0, njs1 
    565                DO ji = 1, jpi 
    566                   obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vsmsk(ji,1) 
    567                END DO 
    568             END DO 
    569          END IF 
    570          IF( lk_mpp )   CALL mpp_sum( obcsurftot )   ! sum over the global domain 
    571       ENDIF 
    572 # endif 
     388      END IF 
     389 
     390      IF ( ln_vol_cst .OR. lk_dynspg_flt ) THEN 
     391        ! ... Initialize obcumask and obcvmask for the Force filtering  
     392        !     boundary condition in dynspg_flt 
     393        obcumask(:,:) = umask(:,:,1) 
     394        obcvmask(:,:) = vmask(:,:,1) 
     395 
     396        ! ... Initialize obctmsk on overlap region and obcs. This mask 
     397        !     is used in obcvol.F90 to calculate cumulate flux E-P.  
     398        !     obc Tracer point are outside the domain ( U/V obc points) ==> masked by obctmsk 
     399        !     - no flux E-P on obcs and overlap region (jpreci = jprecj = 1) 
     400        obctmsk(:,:) = tmask_i(:,:)      
     401 
     402        IF( lp_obc_east ) THEN 
     403           ! ... East obc Force filtering mask for the grad D 
     404           obcumask(nie0  :nie1  ,nje0p1:nje1m1) = 0.e0 
     405           obcvmask(nie0p1:nie1p1,nje0p1:nje1m1) = 0.e0 
     406           ! ... set to 0 on East OBC 
     407           obctmsk(nie0p1:nie1p1,nje0:nje1) = 0.e0 
     408        END IF 
     409   
     410        IF( lp_obc_west ) THEN 
     411           ! ... West obc Force filtering mask for the grad D 
     412           obcumask(niw0:niw1,njw0:njw1) = 0.e0 
     413           obcvmask(niw0:niw1,njw0:njw1) = 0.e0 
     414           ! ... set to 0 on West OBC 
     415           obctmsk(niw0:niw1,njw0:njw1) = 0.e0 
     416        END IF 
     417   
     418        IF( lp_obc_north ) THEN 
     419           ! ... North obc Force filtering mask for the grad D 
     420           obcumask(nin0p1:nin1m1,njn0p1:njn1p1) = 0.e0 
     421           obcvmask(nin0p1:nin1m1,njn0  :njn1  ) = 0.e0 
     422           ! ... set to 0 on North OBC 
     423           obctmsk(nin0:nin1,njn0p1:njn1p1) = 0.e0 
     424        END IF 
     425   
     426        IF( lp_obc_south ) THEN 
     427           ! ... South obc Force filtering mask for the grad D 
     428           obcumask(nis0p1:nis1m1,njs0:njs1) = 0.e0 
     429           obcvmask(nis0p1:nis1m1,njs0:njs1) = 0.e0 
     430           ! ... set to 0 on South OBC 
     431           obctmsk(nis0:nis1,njs0:njs1) = 0.e0 
     432        END IF 
     433      ENDIF 
     434 
     435      IF (lk_dynspg_rl ) THEN  
     436        ! do nothing particular 
     437      ElSE  !  one of the available free surfaces 
     438 
     439        IF ( ln_vol_cst .OR. lk_dynspg_flt  ) THEN 
     440 
     441           ! 3.1 Total lateral surface  
     442           ! ------------------------- 
     443           obcsurftot = 0.e0 
     444   
     445           IF( lp_obc_east ) THEN ! ... East open boundary lateral surface 
     446              DO ji = nie0, nie1 
     447                 DO jj = 1, jpj  
     448                    obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
     449                 END DO 
     450              END DO 
     451           END IF 
     452   
     453           IF( lp_obc_west ) THEN ! ... West open boundary lateral surface 
     454              DO ji = niw0, niw1 
     455                 DO jj = 1, jpj  
     456                    obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
     457                 END DO 
     458              END DO 
     459           END IF 
     460   
     461           IF( lp_obc_north ) THEN ! ... North open boundary lateral surface 
     462              DO jj = njn0, njn1 
     463                 DO ji = 1, jpi 
     464                    obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vnmsk(ji,1) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 
     465                 END DO 
     466              END DO 
     467           END IF 
     468   
     469           IF( lp_obc_south ) THEN ! ... South open boundary lateral surface 
     470              DO jj = njs0, njs1 
     471                 DO ji = 1, jpi 
     472                    obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vsmsk(ji,1) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 
     473                 END DO 
     474              END DO 
     475           END IF 
     476   
     477           IF( lk_mpp )   CALL mpp_sum( obcsurftot )   ! sum over the global domain 
     478        ENDIF 
     479      ENDIF  ! rigid lid 
    573480 
    574481      ! 5. Control print on mask  
     
    696603      ! -------------------------------------------------------------- 
    697604      !   only if at least one boundary is  radiative  
    698  
    699       ! ... Restart from restart.obc 
    700605      IF ( inumfbc < nbobc .AND.  ln_rstart ) THEN 
    701          CALL obc_rst_lec 
     606         !  Restart from restart.obc 
     607         CALL obc_rst_read 
    702608      ELSE 
    703609 
    704           ! ... Initialization to zero of radiation arrays. 
    705           !     Those have dimensions of local subdomains 
     610!         ! ... Initialization to zero of radiation arrays. 
     611!         !     Those have dimensions of local subdomains 
    706612 
    707613          bebnd(:,:,:)   = 0.e0   ;   bnbnd(:,:,:)   = 0.e0 
    708614          uebnd(:,:,:,:) = 0.e0   ;   unbnd(:,:,:,:) = 0.e0 
    709615          vebnd(:,:,:,:) = 0.e0   ;   vnbnd(:,:,:,:) = 0.e0 
    710           tebnd(:,:,:,:) = 0.e0   ;   tnbnd(:,:,:,:) = 0.e0  
     616          tebnd(:,:,:,:) = 0.e0   ;   tnbnd(:,:,:,:) = 0.e0 
    711617          sebnd(:,:,:,:) = 0.e0   ;   snbnd(:,:,:,:) = 0.e0 
    712618 
     
    714620          uwbnd(:,:,:,:) = 0.e0   ;   usbnd(:,:,:,:) = 0.e0 
    715621          vwbnd(:,:,:,:) = 0.e0   ;   vsbnd(:,:,:,:) = 0.e0 
    716           twbnd(:,:,:,:) = 0.e0   ;   tsbnd(:,:,:,:) = 0.e0  
     622          twbnd(:,:,:,:) = 0.e0   ;   tsbnd(:,:,:,:) = 0.e0 
    717623          swbnd(:,:,:,:) = 0.e0   ;   ssbnd(:,:,:,:) = 0.e0 
    718624 
    719625      END IF 
    720626 
    721 # if defined key_dynspg_rl 
    722627      ! 7. Isolated coastline arrays initialization (rigid lid case only) 
    723628      ! ----------------------------------------------------------------- 
    724       CALL obc_dom 
    725 # endif 
     629      IF (lk_dynspg_rl )      CALL obc_dom 
    726630 
    727631      ! 8. Control print 
  • trunk/NEMO/OPA_SRC/OBC/obcrst.F90

    r1146 r1151  
    2121 
    2222   !! * Accessibility 
    23    PUBLIC obc_rst_lec        ! routine called by iniobc.F90 
    24    PUBLIC obc_rst_wri        ! routine called by step.F90 
     23   PUBLIC obc_rst_read       ! routine called by obc_ini 
     24   PUBLIC obc_rst_write      ! routine called by step 
    2525 
    2626   !!--------------------------------------------------------------------------------- 
     
    3232CONTAINS 
    3333 
    34    SUBROUTINE obc_rst_wri ( kt ) 
     34   SUBROUTINE obc_rst_write ( kt ) 
    3535      !!-------------------------------------------------------------------------------- 
    36       !!                  ***  SUBROUTINE obc_rst_wri  *** 
     36      !!                  ***  SUBROUTINE obc_rst_write  *** 
    3737      !!                 
    3838      !! ** Purpose :   Write open boundary restart fields in restart.obc.output file  
     
    4141      !!      Each nstock time step , save fields which are necessary for restart. 
    4242      !!      - This routine is called if at least the key_obc is defined. It is called 
    43       !!        at the same time step than rstwri. 
     43      !!        at the same time step than rstwrite. 
    4444      !!      - First record holds OBC parameters nbobc,jpieob,jpiwob,jpjnob,jpjsob and  
    4545      !!        the OBC layout jpjed, jpjef ... for checking purposes. 
     
    8282         IF(lwp) THEN 
    8383              WRITE(numout,*) ' ' 
    84               WRITE(numout,*) 'obcrst: OBC output for restart with obc_rst_wri routine' 
     84              WRITE(numout,*) 'obcrst: OBC output for restart with obc_rst_write routine' 
    8585              WRITE(numout,*) '~~~~~~' 
    8686              WRITE(numout,*) '        output done in restart.obc.output file at it= ', kt, ' date= ', ndastp 
     
    283283            END IF 
    284284         END IF 
    285       END IF 
    286285      CLOSE(inum) 
    287  
    288    END SUBROUTINE obc_rst_wri 
    289  
    290  
    291    SUBROUTINE obc_rst_lec 
     286      END IF 
     287 
     288   END SUBROUTINE obc_rst_write 
     289 
     290 
     291   SUBROUTINE obc_rst_read 
    292292      !!---------------------------------------------------------------------------- 
    293       !!                   ***  SUBROUTINE obc_rst_lec  *** 
     293      !!                   ***  SUBROUTINE obc_rst_read  *** 
    294294      !!                    
    295295      !! ** Purpose :   Read files for restart at open boundaries 
     
    337337  
    338338      IF(lwp) THEN 
    339          WRITE(numout,*) 'obcrst: beginning of restart with obc_rst_lec routine' 
     339         WRITE(numout,*) 'obcrst: beginning of restart with obc_rst_read routine' 
    340340         WRITE(numout,*) '~~~~~~' 
    341341         WRITE(numout,*) ' ' 
     
    348348      ! 0.1 Open files 
    349349      ! --------------- 
    350       CALL ctlopn( inum, 'restart.obc.output', 'UNKNOWN', 'UNFORMATTED', 'DIRECT',   & 
     350      CALL ctlopn( inum, 'restart.obc', 'UNKNOWN', 'UNFORMATTED', 'DIRECT',   & 
    351351         &         nreclo, numout, lwp, 1 ) 
    352352 
     
    370370          CALL ctl_stop( '        ===>>>> : problem with nit000 for the restart',   & 
    371371               &         '        ==============',   & 
    372                &         '        we stop in obc_rst_lec routine. Verify the file or rerun with the value',   & 
     372               &         '        we stop in obc_rst_read routine. Verify the file or rerun with the value',   & 
    373373               &         '        0 for the control of time parameter nrstdt' ) 
    374374              
     
    391391            WRITE(numout,*) '         ' 
    392392            WRITE(numout,*) '        East open boundary' 
    393             IF( jpieob0 /= jpieob1 ) CALL ctl_stop( '         ==>>>> : Problem in obc_rst_lec, jpieob have changed' ) 
     393            IF( jpieob0 /= jpieob1 ) CALL ctl_stop( '         ==>>>> : Problem in obc_rst_read, jpieob have changed' ) 
    394394         END IF 
    395395      END IF 
     
    399399            WRITE(numout,*) '         ' 
    400400            WRITE(numout,*) '        West open boundary' 
    401             IF( jpiwob0 /= jpiwob1 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpiwob has changed' ) 
     401            IF( jpiwob0 /= jpiwob1 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpiwob has changed' ) 
    402402         END IF 
    403403      END IF 
     
    407407            WRITE(numout,*) '         ' 
    408408            WRITE(numout,*) '        North open boundary' 
    409             IF( jpjnob0 /= jpjnob1 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpjnob has changed' ) 
     409            IF( jpjnob0 /= jpjnob1 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpjnob has changed' ) 
    410410         END IF 
    411411      END IF 
     
    415415            WRITE(numout,*) '         ' 
    416416            WRITE(numout,*) '        South open boundary' 
    417             IF( jpjsob0 /= jpjsob1) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpjsob has changed' ) 
     417            IF( jpjsob0 /= jpjsob1) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpjsob has changed' ) 
    418418         END IF 
    419419      END IF 
     
    423423      ! ------------------------------------------ 
    424424      IF( lp_obc_east .AND. ( jpieob1 /= 0 ) ) THEN 
    425          IF( ied1 /= ied0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpjed has changed' ) 
    426          IF( ief1 /= ief0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpjef has changed' ) 
     425         IF( ied1 /= ied0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpjed has changed' ) 
     426         IF( ief1 /= ief0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpjef has changed' ) 
    427427      END IF 
    428428 
    429429      IF( lp_obc_west .AND. ( jpiwob1 /= 0 ) ) THEN 
    430          IF( iwd1 /= iwd0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpjwd has changed' ) 
    431          IF( iwf1 /= iwf0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpjwf has changed' ) 
     430         IF( iwd1 /= iwd0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpjwd has changed' ) 
     431         IF( iwf1 /= iwf0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpjwf has changed' ) 
    432432      END IF 
    433433  
    434434      IF( lp_obc_north .AND. ( jpjnob1 /= 0 ) ) THEN 
    435          IF( ind1 /= ind0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpind has changed' ) 
    436          IF( inf1 /= inf0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpinf has changed' ) 
     435         IF( ind1 /= ind0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpind has changed' ) 
     436         IF( inf1 /= inf0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpinf has changed' ) 
    437437      END IF 
    438438  
    439439      IF( lp_obc_south .AND. ( jpjsob1 /= 0 ) ) THEN 
    440          IF( isd1 /= isd0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpisd has changed' ) 
    441          IF( isf1 /= isf0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_lec, jpisf has changed' ) 
     440         IF( isd1 /= isd0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpisd has changed' ) 
     441         IF( isf1 /= isf0 ) CALL ctl_stop( '        ==>>>> : Problem in obc_rst_read, jpisf has changed' ) 
    442442      END IF 
    443443  
     
    653653      ENDIF 
    654654  
    655    END SUBROUTINE obc_rst_lec 
     655   END SUBROUTINE obc_rst_read 
    656656#else 
    657657   !!================================================================================= 
     
    660660   !!================================================================================= 
    661661CONTAINS 
    662    SUBROUTINE obc_rst_wri( kt )           !  No Open boundary ==> empty routine 
     662   SUBROUTINE obc_rst_write( kt )           !  No Open boundary ==> empty routine 
    663663      INTEGER,INTENT(in) :: kt 
    664       WRITE(*,*) 'obc_rst_wri: You should not have seen this print! error?', kt 
    665    END SUBROUTINE obc_rst_wri 
    666    SUBROUTINE obc_rst_lec                 !  No Open boundary ==> empty routine 
    667    END SUBROUTINE obc_rst_lec 
     664      WRITE(*,*) 'obc_rst_write: You should not have seen this print! error?', kt 
     665   END SUBROUTINE obc_rst_write 
     666   SUBROUTINE obc_rst_read                 !  No Open boundary ==> empty routine 
     667   END SUBROUTINE obc_rst_read 
    668668#endif 
    669669 
  • trunk/NEMO/OPA_SRC/OBC/obcvol.F90

    r1146 r1151  
    1212   USE oce             ! ocean dynamics and tracers  
    1313   USE dom_oce         ! ocean space and time domain  
    14    USE sbc_oce         ! surface boundary condition: ocean 
     14   USE sbc_oce         ! ocean surface boundary conditions 
    1515   USE phycst          ! physical constants 
    1616   USE obc_oce         ! ocean open boundary conditions 
     
    2929   !!--------------------------------------------------------------------------------- 
    3030   !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    31    !! $Id$ 
     31   !! $Header$  
    3232   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3333   !!--------------------------------------------------------------------------------- 
     
    9696      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain. 
    9797      ! --------------------------------------------------------------------------- 
    98   
    99       zCflxemp = 0.e0 
    100  
    101       DO jj = 1, jpj 
    102          DO ji = 1, jpi 
    103             zCflxemp = zCflxemp + ( (emp(ji,jj)*obctmsk(ji,jj) )/rauw)*e1v(ji,jj)*e2u(ji,jj) 
    104          END DO 
    105       END DO 
     98 
     99      zCflxemp = SUM ( emp(:,:)*obctmsk(:,:)* e1t(:,:) * e2t(:,:)  / rauw ) 
     100 
    106101      IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain 
    107102 
     
    110105 
    111106      zubtpecor = 0.e0 
     107 
     108      ! ... East open boundary 
     109      IF( lp_obc_east ) THEN                      ! ... Total transport through the East OBC 
     110         DO ji = fs_nie0, fs_nie1 ! Vector opt. 
     111            DO jk = 1, jpkm1 
     112               DO jj = 1, jpj 
     113                  zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * & 
     114             &     uemsk(jj,jk)*MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
     115               END DO 
     116            END DO 
     117         END DO 
     118      END IF  
    112119 
    113120      ! ... West open boundary 
     
    116123            DO jk = 1, jpkm1 
    117124               DO jj = 1, jpj 
    118                   zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * uwmsk(jj,jk) 
    119                END DO 
    120             END DO 
    121          END DO 
    122       END IF  
    123  
    124       ! ... East open boundary 
    125       IF( lp_obc_east ) THEN                      ! ... Total transport through the East OBC 
    126          DO ji = fs_nie0, fs_nie1 ! Vector opt. 
    127             DO jk = 1, jpkm1 
    128                DO jj = 1, jpj 
    129                   zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * uemsk(jj,jk) 
    130                END DO 
    131             END DO 
    132          END DO 
    133       END IF  
     125                  zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * & 
     126             &    uwmsk(jj,jk) *MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
     127               END DO 
     128            END DO 
     129         END DO 
     130       ENDIF 
    134131 
    135132      ! ... North open boundary 
     
    138135            DO jk = 1, jpkm1 
    139136               DO ji = 1, jpi 
    140                   zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * vnmsk(ji,jk) 
    141                END DO 
    142             END DO 
    143          END DO 
    144       END IF  
     137                  zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * & 
     138             &    vnmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 
     139               END DO 
     140            END DO 
     141         END DO 
     142       ENDIF 
    145143 
    146144      ! ... South open boundary 
     
    149147            DO jk = 1, jpkm1 
    150148               DO ji = 1, jpi 
    151                   zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * vsmsk(ji,jk) 
    152                END DO 
    153             END DO 
    154          END DO 
    155       END IF  
     149                  zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * & 
     150             &    vsmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 
     151               END DO 
     152            END DO 
     153         END DO 
     154       ENDIF 
    156155 
    157156      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
     
    160159      ! 3. The normal velocity correction 
    161160      ! --------------------------------- 
    162  
    163       zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot) 
    164  
    165161      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
    166162         IF(lwp) WRITE(numout,*)'        ' 
    167163         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 
    168164         IF(lwp) WRITE(numout,*)'~~~~~~~ ' 
    169          IF(lwp) WRITE(numout,*)'        ' 
    170165         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)' 
     166         IF(lwp) WRITE(numout,*)'          lateral transport :',zubtpecor,'(m3/s)' 
     167         IF(lwp) WRITE(numout,*)'          net inflow        :',zubtpecor-zCflxemp,'(m3/s)' 
     168      ENDIF 
     169 
     170      zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot) 
     171 
     172      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
    171173         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',obcsurftot,'(m2)' 
    172174         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)' 
     
    175177 
    176178      ! 4. Correction of the total velocity on each open  
    177       !    boundary torespect the mass flux conservation 
     179      !    boundary to respect the mass flux conservation 
    178180      ! ------------------------------------------------- 
    179181 
    180       ztransw = 0.e0 
    181       ztranse = 0.e0 
    182       ztransn = 0.e0 
    183       ztranss = 0.e0 
    184       ztranst = 0.e0 
     182      ztranse = 0.e0   ; ztransw = 0.e0 ; ztransn = 0.e0 ; ztranss = 0.e0 
     183      ztranst = 0.e0  ! total 
    185184 
    186185      IF( lp_obc_west ) THEN 
    187  
    188186         ! ... correction of the west velocity 
    189187         DO ji = fs_niw0, fs_niw1 ! Vector opt. 
     
    191189               DO jj = 1, jpj 
    192190                  ua(ji,jj,jk) = ua(ji,jj,jk) - zubtpecor*uwmsk(jj,jk) 
    193                   ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk) 
     191                  ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk) * & 
     192             &    MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
    194193               END DO 
    195194            END DO 
     
    198197         IF( lk_mpp )   CALL mpp_sum( ztransw )   ! sum over the global domain 
    199198 
    200          IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
    201             IF(lwp) WRITE(numout,*)'          West OB transport ztransw :', ztransw,'(m3/s)' 
    202          END IF  
    203  
     199         IF( lwp .AND. MOD( kt, nwrite ) == 0)  WRITE(numout,*)'          West OB transport ztransw :', ztransw,'(m3/s)' 
    204200      END IF  
    205201 
     
    211207               DO jj = 1, jpj 
    212208                  ua(ji,jj,jk) = ua(ji,jj,jk) + zubtpecor*uemsk(jj,jk) 
    213                   ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk) 
     209                  ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk) * & 
     210            &     MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
    214211               END DO 
    215212            END DO 
     
    231228               DO ji =  1, jpi 
    232229                  va(ji,jj,jk) = va(ji,jj,jk) + zubtpecor*vnmsk(ji,jk) 
    233                   ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk) 
     230                  ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk) * & 
     231           &      MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 
    234232               END DO 
    235233            END DO 
     
    250248               DO ji =  1, jpi 
    251249                  va(ji,jj,jk) = va(ji,jj,jk) - zubtpecor*vsmsk(ji,jk) 
    252                   ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk) 
     250                  ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk) * & 
     251            &     MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 
    253252               END DO 
    254253            END DO 
     
    266265      ! ------------------------------------------- 
    267266 
    268       ztranst = ztransw - ztranse + ztranss - ztransn 
    269267 
    270268      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
     269         ztranst = ztransw - ztranse + ztranss - ztransn 
    271270         IF(lwp) WRITE(numout,*)'        ' 
    272271         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)' 
     272         IF(lwp) WRITE(numout,*)'          Balance  =', ztranst - zCflxemp ,'(m3/s)' 
    273273         IF(lwp) WRITE(numout,*)'        ' 
    274274      END IF  
  • trunk/NEMO/OPA_SRC/step.F90

    r1146 r1151  
    332332      IF( indic < 0          )   CALL ctl_stop( 'step: indic < 0' ) 
    333333      IF( kstp == nit000     )   CALL iom_close( numror )             ! close input  ocean restart file 
    334       IF( lrst_oce           )   CALL rst_write  ( kstp )             ! write output ocean restart file 
    335       IF( lk_obc             )   CALL obc_rst_wri( kstp )             ! write open boundary restart file 
     334      IF( lrst_oce           )   CALL rst_write    ( kstp )           ! write output ocean restart file 
     335      IF( lk_obc             )   CALL obc_rst_write( kstp )           ! write open boundary restart file 
    336336 
    337337      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.