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 2166 for trunk/NEMO/OPA_SRC/OBC/obcdta.F90 – NEMO

Ignore:
Timestamp:
2010-10-06T15:29:14+02:00 (13 years ago)
Author:
rblod
Message:

Fix various bugs in OBC see ticket #716 #548 #296

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/OBC/obcdta.F90

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