Changeset 1125


Ignore:
Timestamp:
2008-06-23T11:05:02+02:00 (12 years ago)
Author:
ctlod
Message:

trunk: BDY package code review (coding rules), see ticket: #214

Location:
trunk/NEMO/OPA_SRC
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r1058 r1125  
    44   !! Unstructured Open Boundary Cond. :   define related variables 
    55   !!====================================================================== 
    6 #if defined key_bdy || defined key_bdy_tides 
     6   !! History :  1.0  !  2001-05  (J. Chanut, A. Sellar)  Original code 
     7   !!            3.0  !  2008-04  (NEMO team)  add in the reference version      
    78   !!---------------------------------------------------------------------- 
    8    !!   'key_bdy' :                    Unstructured Open Boundary Condition 
     9#if defined key_bdy  
    910   !!---------------------------------------------------------------------- 
    10    !! History : 
    11    !!  9.0   01/05   (J. Chanut, A. Sellar)  Original code 
     11   !!   'key_bdy'                      Unstructured Open Boundary Condition 
    1212   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1413   USE par_oce         ! ocean parameters 
    1514   USE bdy_par         ! Unstructured boundary parameters 
     
    2120   !! Namelist variables 
    2221   !!---------------------------------------------------------------------- 
    23  
    24    CHARACTER(len=80) :: & 
    25       filbdy_mask, &         !: Name of unstruct. bdy mask file 
    26       filbdy_data_T, &       !: Name of unstruct. bdy data file at T points 
    27       filbdy_data_U, &       !: Name of unstruct. bdy data file at U points 
    28       filbdy_data_V, &       !: Name of unstruct. bdy data file at V points 
    29       filbdy_data_bt_T, &    !: Name of unstruct. bdy data file at T points for barotropic variables 
    30       filbdy_data_bt_U, &    !: Name of unstruct. bdy data file at U points for barotropic variables 
    31       filbdy_data_bt_V       !: Name of unstruct. bdy data file at V points for barotropic variables 
    32  
    33    LOGICAL ::  & 
    34       ln_bdy_clim =.false.!: if true, we assume that bdy data files contain  
    35                           !  1 time dump (-->bdy forcing will be constant)  
    36                           !  or 12 months  (-->bdy forcing will be cyclic)  
    37    LOGICAL ::  & 
    38       ln_bdy_fla =.false. !: if true, flather boundary conditions 
    39  
    40    LOGICAL ::  & 
    41       ln_bdy_vol =.false. !: if true, volume correction              
    42  
    43    LOGICAL ::  & 
    44       ln_bdy_mask = .false. !: if true, read bdymask from file 
    45  
    46    INTEGER ::  & 
    47       nb_rimwidth = 7,  & !: boundary rim width 
    48       nbdy_dta = 1     !: = 0 use the initial state as bdy dta 
    49                           !: = 1 read bdy data in netcdf file 
    50    INTEGER ::  &  
    51       volbdy = 1          !: = 0 the total volume will have the variability  
    52                           !       of the surface Flux E-P else (volbdy = 1)  
    53                           !       the volume will be constant 
    54                           !  = 1 the volume will be constant during all the  
    55                           !       integration. 
     22   CHARACTER(len=80) ::   filbdy_mask        !: Name of unstruct. bdy mask file 
     23   CHARACTER(len=80) ::   filbdy_data_T      !: Name of unstruct. bdy data file at T points 
     24   CHARACTER(len=80) ::   filbdy_data_U      !: Name of unstruct. bdy data file at U points 
     25   CHARACTER(len=80) ::   filbdy_data_V      !: Name of unstruct. bdy data file at V points 
     26   CHARACTER(len=80) ::   filbdy_data_bt_T   !: Name of unstruct. bdy data file at T points for barotropic variables 
     27   CHARACTER(len=80) ::   filbdy_data_bt_U   !: Name of unstruct. bdy data file at U points for barotropic variables 
     28   CHARACTER(len=80) ::   filbdy_data_bt_V   !: Name of unstruct. bdy data file at V points for barotropic variables 
     29   ! 
     30   LOGICAL ::   ln_bdy_tides = .false.  !: =T apply tidal harmonic forcing along open boundaries 
     31   LOGICAL ::   ln_bdy_vol  = .false.   !: =T volume correction              
     32   LOGICAL ::   ln_bdy_mask = .false.   !: =T read bdymask from file 
     33   LOGICAL ::   ln_bdy_clim = .false.   !: if true, we assume that bdy data files contain  
     34   !                                    !  1 time dump  (-->bdy forcing will be constant)  
     35   !                                    !  or 12 months (-->bdy forcing will be cyclic)  
     36   LOGICAL ::   ln_bdy_dyn_fla  = .false. !: =T Flather boundary conditions on barotropic velocities 
     37   LOGICAL ::   ln_bdy_dyn_frs  = .false. !: =T FRS boundary conditions on velocities 
     38   LOGICAL ::   ln_bdy_tra_frs  = .false. !: =T FRS boundary conditions on tracers (T and S) 
     39   LOGICAL ::   ln_bdy_ice_frs  = .false. !: =T FRS boundary conditions on seaice (leads fraction, ice depth, snow depth) 
     40   ! 
     41   INTEGER ::   nb_rimwidth = 7         !: boundary rim width 
     42   INTEGER ::   nbdy_dta    = 1          !: = 0 use the initial state as bdy dta or = 1 read it in a NetCDF file 
     43   INTEGER ::   volbdy      = 1         !: = 0 the total volume will have the variability of the surface Flux E-P  
     44   !                                    !  = 1 the volume will be constant during all the integration. 
    5645 
    5746   !!---------------------------------------------------------------------- 
    5847   !! Global variables 
    5948   !!---------------------------------------------------------------------- 
    60  
    61    REAL(wp), DIMENSION(jpi, jpj) ::  & !: 
    62       bdytmask, &  !: Mask defining computational domain at T-points 
    63       bdyumask, &  !: Mask defining computational domain at U-points 
    64       bdyvmask     !: Mask defining computational domain at V-points 
     49   REAL(wp), DIMENSION(jpi,jpj) ::   bdytmask   !: Mask defining computational domain at T-points 
     50   REAL(wp), DIMENSION(jpi,jpj) ::   bdyumask   !: Mask defining computational domain at U-points 
     51   REAL(wp), DIMENSION(jpi,jpj) ::   bdyvmask   !: Mask defining computational domain at V-points 
    6552 
    6653   !!---------------------------------------------------------------------- 
    6754   !! Unstructured open boundary data variables 
    6855   !!---------------------------------------------------------------------- 
     56   INTEGER, DIMENSION(jpbgrd) ::   nblen                  !: Size of bdy data on a proc for each grid type 
     57   INTEGER, DIMENSION(jpbgrd) ::   nblenrim               !: Size of bdy data on a proc for first rim ind 
     58   INTEGER, DIMENSION(jpbgrd) ::   nblendta               !: Size of bdy data in file 
    6959 
    70    INTEGER, DIMENSION(jpbgrd) ::  & 
    71       nblen  = 0, & !: Size of bdy data on a proc for each grid type 
    72       nblenrim = 0,  & !: Size of bdy data on a proc for first rim ind 
    73       nblendta = 0        !: Size of bdy data in file 
     60   INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbi, nbj        !: i and j indices of bdy dta 
     61   INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbr             !: Discrete distance from rim points 
     62   INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbmap           !: Indices of data in file for data in memory  
     63     
     64   REAL(wp) ::   bdysurftot                             !: Lateral surface of unstructured open boundary 
    7465 
    75    INTEGER, DIMENSION(jpbdim, jpbgrd) ::  &  
    76       nbi, nbj,      & !: i and j indices of bdy dta 
    77       nbr,     & !: Discrete distance from rim points 
    78       nbmap      !: Indices of data in file for data in memory  
    79      
    80    REAL(wp) :: & 
    81       bdysurftot    !: Lateral surface of unstructured open boundary 
     66   REAL(wp), DIMENSION(jpbdim)        ::   flagu, flagv   !: Flag for normal velocity compnt for velocity components 
     67   REAL(wp), DIMENSION(jpbdim,jpbgrd) ::   nbw            !: Rim weights of bdy data 
    8268 
    83    REAL(wp), DIMENSION(jpbdim) ::  & 
    84       flagu,      & !: Flag for normal velocity compnt for u-velocity 
    85       flagv      !: Flag for normal velocity compnt for v-velocity 
     69   REAL(wp), DIMENSION(jpbdim)     ::   sshbdy            !: Now clim of bdy sea surface height (Flather) 
     70   REAL(wp), DIMENSION(jpbdim)     ::   ubtbdy, vbtbdy    !: Now clim of bdy barotropic velocity components 
     71   REAL(wp), DIMENSION(jpbdim,jpk) ::   tbdy  , sbdy      !: Now clim of bdy temperature and salinity   
     72   REAL(wp), DIMENSION(jpbdim,jpk) ::   ubdy  , vbdy    !: Now clim of bdy velocity components 
     73   REAL(wp), DIMENSION(jpbdim) ::   sshtide               !: Tidal boundary array : SSH 
     74   REAL(wp), DIMENSION(jpbdim) ::   utide, vtide          !: Tidal boundary array : U and V 
    8675 
    87    REAL(wp), DIMENSION(jpbdim, jpbgrd) ::  & 
    88       nbw        !: Rim weights of bdy data 
    89  
    90    REAL(wp), DIMENSION(jpbdim) ::  & 
    91       sshbdy,     & !: Now clim of bdy sea surface height (Flather) 
    92       ubtbdy, vbtbdy   !: Now clim of bdy barotropic velocity components 
    93  
    94    REAL(wp), DIMENSION(jpbdim,jpk) ::  & 
    95       tbdy, sbdy, & !: Now clim of bdy temperature and salinity   
    96       ubdy, vbdy    !: Now clim of bdy velocity components 
    97  
    98 #if defined key_bdy_tides 
    99    REAL(wp), DIMENSION(jpbdim) ::  & 
    100       sshtide,          & !: Tidal boundary array : SSH 
    101       utide,            & !: Tidal boundary array : U 
    102       vtide           !: Tidal boundary array : V 
    103 #endif 
    104            
    10576#else 
    10677   !!---------------------------------------------------------------------- 
    107    !!   Default option :                                       Empty module 
     78   !!   Dummy module                NO Unstructured Open Boundary Condition 
    10879   !!---------------------------------------------------------------------- 
    10980#endif 
    11081 
     82   !!---------------------------------------------------------------------- 
     83   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     84   !! $Id: $  
     85   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    11186   !!====================================================================== 
    11287END MODULE bdy_oce 
  • trunk/NEMO/OPA_SRC/BDY/bdy_par.F90

    r1058 r1125  
    11MODULE bdy_par 
    22   !!====================================================================== 
    3    !!                  ***  MODULE bdy_par   *** 
     3   !!                      ***  MODULE bdy_par   *** 
    44   !! Unstructured Open Boundary Cond. :   define related parameters 
    55   !!====================================================================== 
    6 #if defined key_bdy || defined key_bdy_tides 
     6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
     7   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     8   !!---------------------------------------------------------------------- 
     9#if defined key_bdy 
    710   !!---------------------------------------------------------------------- 
    811   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    912   !!---------------------------------------------------------------------- 
    10    !! history : 
    11    !!  9.0   01/05   (J. Chanut, A. Sellar)  Original code 
    12    !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    14    USE par_oce         ! ocean parameters 
    1513 
    1614   IMPLICIT NONE 
    1715   PUBLIC 
    1816 
    19 #if defined key_bdy 
    20    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .TRUE. !: Unstructured Ocean  
    21                                                    !Boundary Condition flag 
     17 
     18   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy  = .TRUE.  !: Unstructured Ocean Boundary Condition flag 
     19   INTEGER, PUBLIC, PARAMETER ::   jpbdta  = 5000    !: Max length of bdy field in file 
     20   INTEGER, PUBLIC, PARAMETER ::   jpbdim  = 5000    !: Max length of bdy field on a processor 
     21   INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000    !: Max number of time dumps per file 
     22   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3       !: Number of horizontal grid types used  (T, u, v, f) 
    2223#else 
    23    ! tides only at boundaries 
    24    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .FALSE. !: Unstructured Ocean  
    25                                                     !Boundary Condition flag 
     24   !!---------------------------------------------------------------------- 
     25   !!   Default option :            NO Unstructured open boundary condition 
     26   !!---------------------------------------------------------------------- 
     27   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .FALSE.   !: Unstructured Ocean Boundary Condition flag 
    2628#endif 
    2729 
    2830   !!---------------------------------------------------------------------- 
    29    !! Unstructured open boundary parameters 
    30    !!---------------------------------------------------------------------- 
    31  
    32    INTEGER, PARAMETER ::   &  
    33       jpbdta     = 5000,  & !: Max length of bdy field in file 
    34       jpbdim     = 5000,  & !: Max length of bdy field on a processor 
    35       jpbtime    = 1000,  & !: Max number of time dumps per file 
    36       jpbgrd     = 3     !: Number of horizontal grid types used  
    37                             !  (T, u, v, f) 
    38 #else 
    39    !!---------------------------------------------------------------------- 
    40    !!   Default option :                         NO open boundary condition 
    41    !!---------------------------------------------------------------------- 
    42    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .FALSE.!: Unstructured Ocean  
    43                                                    !Boundary Condition flag 
    44 #endif 
    45  
     31   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     32   !! $Id: $  
     33   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4634   !!====================================================================== 
    4735END MODULE bdy_par 
  • trunk/NEMO/OPA_SRC/BDY/bdydta.F90

    r1058 r1125  
    11MODULE bdydta 
    2    !!============================================================================== 
    3    !!                            ***  MODULE bdydta  *** 
     2   !!====================================================================== 
     3   !!                       ***  MODULE bdydta  *** 
    44   !! Open boundary data : read the data for the unstructured open boundaries. 
    5    !!============================================================================== 
    6 #if defined key_bdy || defined key_bdy_tides 
    7    !!------------------------------------------------------------------------------ 
    8    !!   'key_bdy'         :                   Unstructured Open Boundary Conditions 
    9    !!------------------------------------------------------------------------------ 
    10    !!   bdy_dta           : read u, v, t, s data along each open boundary 
    11    !!------------------------------------------------------------------------------ 
    12    !! * Modules used 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
     7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module 
     8   !!             -   !  2007-07  (D. Storkey) add bdy_dta_bt 
     9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     10   !!---------------------------------------------------------------------- 
     11#if defined key_bdy 
     12   !!---------------------------------------------------------------------- 
     13   !!   'key_bdy'                     Unstructured Open Boundary Conditions 
     14   !!---------------------------------------------------------------------- 
     15   !!   bdy_dta    : read u, v, t, s data along open boundaries 
     16   !!   bdy_dta_bt : read depth-mean velocities and elevation along open 
     17   !!                boundaries         
     18   !!---------------------------------------------------------------------- 
    1319   USE oce             ! ocean dynamics and tracers 
    1420   USE dom_oce         ! ocean space and time domain 
     
    2430   PRIVATE 
    2531 
    26    !! * Accessibility 
    27    PUBLIC bdy_dta                         ! routines called by step.F90 
    28    PUBLIC bdy_dta_bt  
     32   PUBLIC   bdy_dta      ! routines called by step.F90 
     33   PUBLIC   bdy_dta_bt  
     34 
     35   INTEGER ::   numbdyt, numbdyu, numbdyv                      !: logical units for T-, U-, & V-points data file, resp. 
     36   INTEGER ::   ntimes_bdy                                     !: exact number of time dumps in data files 
     37   INTEGER ::   nbdy_b, nbdy_a                                 !: record of bdy data file for before and after model time step 
     38   INTEGER ::   numbdyt_bt, numbdyu_bt, numbdyv_bt             !: logical unit for T-, U- & V-points data file, resp. 
     39   INTEGER ::   ntimes_bdy_bt                                  !: exact number of time dumps in data files 
     40   INTEGER ::   nbdy_b_bt, nbdy_a_bt                           !: record of bdy data file for before and after model time step 
     41 
     42   INTEGER, DIMENSION (jpbtime) ::   istep, istep_bt           !: time array in seconds in each data file 
     43 
     44   REAL(wp) ::  zoffset                                        !: time offset between time origin in file & start time of model run 
     45 
     46   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   tbdydta, sbdydta     !: time interpolated values of T and S bdy data    
     47   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   ubdydta, vbdydta     !: time interpolated values of U and V bdy data  
     48   REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta !: Arrays used for time interpolation of bdy data    
     49   REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            !: bdy data of ssh 
     50 
     51   !!---------------------------------------------------------------------- 
     52   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     53   !! $Id: $  
     54   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     55   !!---------------------------------------------------------------------- 
     56 
     57CONTAINS 
     58 
     59   SUBROUTINE bdy_dta( kt ) 
     60      !!---------------------------------------------------------------------- 
     61      !!                   ***  SUBROUTINE bdy_dta  *** 
     62      !!                     
     63      !! ** Purpose :   Read unstructured boundary data for FRS condition. 
     64      !! 
     65      !! ** Method  :   At the first timestep, read in boundary data for two 
     66      !!                times from the file and time-interpolate. At other  
     67      !!                timesteps, check to see if we need another time from  
     68      !!                the file. If so read it in. Time interpolate. 
     69      !!---------------------------------------------------------------------- 
     70      INTEGER, INTENT( in ) ::   kt                             ! ocean time-step index (for timesplitting option, otherwise zero) 
     71      !! 
     72      CHARACTER(LEN=80), DIMENSION(3) ::   clfile               ! names of input files 
     73      CHARACTER(LEN=70 )              ::   clunits              ! units attribute of time coordinate 
     74      LOGICAL ::   lect                                         ! flag for reading 
     75      INTEGER ::   it, ib, ik, igrd                             ! dummy loop indices 
     76      INTEGER ::   igrd_start, igrd_end                         ! start and end of loops on igrd 
     77      INTEGER ::   idvar                                        ! netcdf var ID 
     78      INTEGER ::   iman, i15, imois                             ! Time variables for monthly clim forcing 
     79      INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 
     80      INTEGER ::   itimer, totime 
     81      INTEGER ::   ii, ij                                       ! array addresses 
     82      INTEGER ::   ipi, ipj, ipk, inum                          ! temporary integers (NetCDF read) 
     83      INTEGER ::   iyear0, imonth0, iday0 
     84      INTEGER ::   ihours0, iminutes0, isec0 
     85      INTEGER ::   iyear, imonth, iday, isecs 
     86      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files 
     87      REAL(wp) ::   dayfrac, zxy, zoffsett 
     88      REAL(wp) ::   zoffsetu, zoffsetv 
     89      REAL(wp) ::   dayjul0, zdayjulini 
     90      REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files 
     91      REAL(wp), DIMENSION(jpbdta,1,jpk) ::   zdta               ! temporary array for data fields 
     92      !!--------------------------------------------------------------------------- 
     93 
     94      IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs ) THEN  ! If these are both false then this routine  
     95                                                     ! does nothing. 
     96 
     97      ! -------------------- ! 
     98      !    Initialization    ! 
     99      ! -------------------- ! 
     100 
     101      lect   = .false.           ! If true, read a time record 
     102 
     103      ! Some time variables for monthly climatological forcing: 
     104      ! ******************************************************* 
     105 !!gm  here  use directely daymod variables 
     106  
     107      iman = INT( raamo )      ! Number of months in a year 
     108 
     109      i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 
     110      ! i15=0 if the current day is in the first half of the month, else i15=1 
     111 
     112      imois = nmonth + i15 - 1            ! imois is the first month record 
     113      IF( imois == 0 )   imois = iman 
     114 
     115      ! Time variable for non-climatological forcing: 
     116      ! ********************************************* 
     117      itimer = (kt-nit000+1)*rdt      ! current time in seconds for interpolation  
     118 
     119 
     120      !                                                !-------------------! 
     121      IF( kt == nit000 ) THEN                          !  First call only  ! 
     122         !                                             !-------------------! 
     123         istep(:) = 0 
     124         nbdy_b    = 0 
     125         nbdy_a    = 0 
     126 
     127         ! Get time information from bdy data file 
     128         ! *************************************** 
     129 
     130         IF(lwp) WRITE(numout,*) 
     131         IF(lwp) WRITE(numout,*)    'bdy_dta : Initialize unstructured boundary data' 
     132         IF(lwp) WRITE(numout,*)    '~~~~~~~'  
     133 
     134         IF     ( nbdy_dta == 0 ) THEN 
     135            ! 
     136            IF(lwp) WRITE(numout,*) '          Bdy data are taken from initial conditions' 
     137            ! 
     138         ELSEIF (nbdy_dta == 1) THEN 
     139            ! 
     140            IF(lwp) WRITE(numout,*) '          Bdy data are read in netcdf files' 
     141            ! 
     142            dayfrac = adatrj  - FLOAT( itimer ) / 86400.   ! day fraction at time step kt-1 
     143            dayfrac = dayfrac - INT  ( dayfrac )           ! 
     144            totime  = ( nitend - nit000 + 1 ) * rdt        ! Total time of the run to verify that all the 
     145            !                                              ! necessary time dumps in file are included 
     146            ! 
     147            clfile(1) = filbdy_data_T 
     148            clfile(2) = filbdy_data_U 
     149            clfile(3) = filbdy_data_V 
     150            !                                                   
     151            ! how many files are we to read in? 
     152            igrd_start = 1 
     153            igrd_end   = 3 
     154            IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN 
     155               ! No T-grid file. 
     156               igrd_start = 2 
     157            ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN 
     158               ! No U-grid or V-grid file. 
     159               igrd_end   = 1          
     160            ENDIF 
     161 
     162            DO igrd = igrd_start, igrd_end                     !  loop over T, U & V grid  ! 
     163               !                                               !---------------------------! 
     164               CALL iom_open( clfile(igrd), inum ) 
     165               CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits )  
     166 
     167               SELECT CASE( igrd ) 
     168                  CASE (1)  
     169                     numbdyt = inum 
     170                  CASE (2)  
     171                     numbdyu = inum 
     172                  CASE (3)  
     173                     numbdyv = inum 
     174               END SELECT 
     175 
     176               ! Calculate time offset  
     177               READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 
     178               ! Convert time origin in file to julian days  
     179               isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 
     180               CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 
     181               ! Compute model initialization time  
     182               iyear  = ndastp / 10000 
     183               imonth = ( ndastp - iyear * 10000 ) / 100 
     184               iday   = ndastp - iyear * 10000 - imonth * 100 
     185               isecs  = dayfrac * 86400 
     186               CALL ymds2ju(iyear, imonth, iday, real(isecs) , zdayjulini) 
     187               ! offset from initialization date: 
     188               zoffset = (dayjul0-zdayjulini)*86400 
     189               ! 
     1907000           FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 
     191 
     192               !! TO BE DONE... Check consistency between calendar from file  
     193               !! (available optionally from iom_gettime) and calendar in model  
     194               !! when calendar in model available outside of IOIPSL. 
     195 
     196               IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy 
     197               IF(lwp) WRITE(numout,*) 'offset: ',zoffset 
     198               IF(lwp) WRITE(numout,*) 'totime: ',totime 
     199               IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr 
     200 
     201               ! Check that there are not too many times in the file.  
     202               IF( ntimes_bdy > jpbtime ) THEN 
     203                  WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy 
     204                  CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 ) 
     205               ENDIF 
     206 
     207               ! Check that time array increases: 
     208         
     209               it = 1 
     210               DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 ) 
     211                 it = it + 1 
     212               END DO 
     213 
     214               IF( it.NE.ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN 
     215                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
     216                     CALL ctl_stop( 'Time array in unstructured boundary data files',   & 
     217                        &           'does not continuously increase.'               , ctmp1 ) 
     218               ENDIF 
     219               ! 
     220               ! Check that times in file span model run time: 
     221               IF( zstepr(1) + zoffset > 0 ) THEN 
     222                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
     223                     CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 ) 
     224               END IF 
     225               IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN 
     226                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 
     227                     CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 ) 
     228               END IF 
     229               ! 
     230               IF    ( igrd == 1 ) THEN 
     231                 ntimes_bdyt = ntimes_bdy 
     232                 zoffsett = zoffset 
     233                 istept(:) = INT( zstepr(:) + zoffset ) 
     234               ELSEIF(igrd == 2 ) THEN 
     235                 ntimes_bdyu = ntimes_bdy 
     236                 zoffsetu = zoffset 
     237                 istepu(:) = INT( zstepr(:) + zoffset ) 
     238               ELSEIF(igrd == 3 ) THEN 
     239                 ntimes_bdyv = ntimes_bdy 
     240                 zoffsetv = zoffset 
     241                 istepv(:) = INT( zstepr(:) + zoffset ) 
     242               ENDIF 
     243               ! 
     244            END DO                                         ! end loop over T, U & V grid  
     245 
     246            IF (igrd_start == 1 .and. igrd_end == 3) THEN 
     247               ! Only test differences if we are reading in 3 files 
     248               ! Verify time consistency between files   
     249               IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 
     250                  CALL ctl_stop( 'Bdy data files must have the same number of time dumps',   & 
     251                  &           'Multiple time frequencies not implemented yet'  ) 
     252               ENDIF 
     253               ntimes_bdy = ntimes_bdyt 
     254               ! 
     255               IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN 
     256                  CALL ctl_stop( 'Bdy data files must have the same time origin',   & 
     257                  &           'Multiple time frequencies not implemented yet' ) 
     258               ENDIF 
     259               zoffset = zoffsett 
     260            ENDIF 
     261 
     262            IF( igrd_start == 1 ) THEN 
     263               istep(:) = istept(:) 
     264            ELSE 
     265               istep(:) = istepu(:) 
     266            ENDIF 
     267 
     268            ! Check number of time dumps:               
     269            IF( ntimes_bdy == 1 .AND. .NOT. ln_bdy_clim ) THEN 
     270              CALL ctl_stop( 'There is only one time dump in data files',   & 
     271                 &           'Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' ) 
     272            ENDIF 
     273 
     274            IF( ln_bdy_clim ) THEN 
     275              IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN 
     276                 CALL ctl_stop( 'For climatological boundary forcing (ln_bdy_clim=.true.),',   & 
     277                    &           'bdy data files must contain 1 or 12 time dumps.' ) 
     278              ELSEIF( ntimes_bdy ==  1 ) THEN 
     279                IF(lwp) WRITE(numout,*) 
     280                IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 
     281              ELSEIF( ntimes_bdy == 12 ) THEN 
     282                IF(lwp) WRITE(numout,*) 
     283                IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 
     284              ENDIF 
     285            ENDIF 
     286 
     287            ! Find index of first record to read (before first model time).  
     288            it = 1 
     289            DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 ) 
     290              it = it + 1 
     291            END DO 
     292            nbdy_b = it 
     293            ! 
     294            WRITE(numout,*) 'Time offset is ',zoffset 
     295            WRITE(numout,*) 'First record to read is ',nbdy_b 
     296 
     297         ENDIF ! endif (nbdy_dta == 1) 
     298 
     299 
     300         ! 1.2  Read first record in file if necessary (ie if nbdy_dta == 1) 
     301         ! ***************************************************************** 
     302 
     303         IF( nbdy_dta == 0) THEN      ! boundary data arrays are filled with initial conditions 
     304            ! 
     305            IF (ln_bdy_tra_frs) THEN 
     306              igrd = 1            ! T-points data  
     307              DO ib = 1, nblen(igrd) 
     308                ii = nbi(ib,igrd) 
     309                ij = nbj(ib,igrd) 
     310                DO ik = 1, jpkm1 
     311                  tbdy(ib,ik) = tn(ii, ij, ik) 
     312                  sbdy(ib,ik) = sn(ii, ij, ik) 
     313                ENDDO 
     314              END DO 
     315            ENDIF 
     316 
     317            IF(ln_bdy_dyn_frs) THEN 
     318              igrd = 2            ! U-points data  
     319              DO ib = 1, nblen(igrd) 
     320                ii = nbi(ib,igrd) 
     321                ij = nbj(ib,igrd) 
     322                DO ik = 1, jpkm1 
     323                  ubdy(ib,ik) = un(ii, ij, ik) 
     324                ENDDO 
     325              END DO 
     326 
     327              igrd = 3            ! V-points data  
     328              DO ib = 1, nblen(igrd)             
     329                ii = nbi(ib,igrd) 
     330                ij = nbj(ib,igrd) 
     331                DO ik = 1, jpkm1 
     332                  vbdy(ib,ik) = vn(ii, ij, ik) 
     333                ENDDO 
     334              END DO 
     335            ENDIF 
     336            ! 
     337         ELSEIF( nbdy_dta == 1 ) THEN    ! Set first record in the climatological case:    
     338            ! 
     339            IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN 
     340               nbdy_a = 1 
     341            ELSEIF( ln_bdy_clim .AND. ntimes_bdy == iman ) THEN 
     342               nbdy_b = 0 
     343               nbdy_a = imois 
     344            ELSE 
     345               nbdy_a = nbdy_b 
     346            ENDIF 
    29347    
    30    !! * local modules variables 
    31    INTEGER ::   & 
    32       bdynumt,         &  ! logical unit for T-points data file 
    33       bdynumu,         &  ! logical unit for U-points data file 
    34       bdynumv,         &  ! logical unit for V-points data file 
    35       kbdy,            &  ! Exact number of time dumps in data files 
    36       nbdy1,           &  ! Time record in bdy data file BEFORE the model current time 
    37       nbdy2               ! Time record in bdy data file AFTER the model current time 
    38  
    39    !! * local modules variables for barotropic variables 
    40    INTEGER ::   & 
    41       bdynumt_bt,         &  ! logical unit for T-points data file 
    42       bdynumu_bt,         &  ! logical unit for U-points data file 
    43       bdynumv_bt,         &  ! logical unit for V-points data file 
    44       kbdy_bt,            &  ! Exact number of time dumps in data files 
    45       nbdy1_bt,           &  ! Time record in bdy data file BEFORE the model current time 
    46       nbdy2_bt               ! Time record in bdy data file AFTER the model current time 
    47  
    48    INTEGER, DIMENSION (jpbtime) ::   &  
    49       istep, istep_bt        ! time array in seconds in each data file 
    50  
    51    REAL(wp) ::  & 
    52       offset              ! time offset between time origin in file and start time of model run 
    53  
    54    REAL(wp), DIMENSION(jpbdim,jpk,2) ::   & 
    55       tbdydta, sbdydta, & !: Arrays used for time interpolation of bdy data    
    56       ubdydta, vbdydta 
    57  
    58    REAL(wp), DIMENSION(jpbdim,2) ::   & 
    59       ubtbdydta, vbtbdydta,   & !: Arrays used for time interpolation of bdy data    
    60       sshbdydta                   !: bdy data of ssh 
    61  
    62    !!--------------------------------------------------------------------------------- 
    63    !!   OPA 9.0 , LODYC-IPSL  (2003) 
    64    !!--------------------------------------------------------------------------------- 
    65  
    66 CONTAINS 
    67  
    68    SUBROUTINE bdy_dta( kt ) 
     348            ! Read first record: 
     349            ipj  = 1 
     350            ipk  = jpk 
     351            igrd = 1 
     352            ipi  = nblendta(igrd) 
     353 
     354            IF(ln_bdy_tra_frs) THEN 
     355               igrd = 1                                           ! Temperature 
     356               IF( nblendta(igrd) <=  0 ) THEN  
     357                  idvar = iom_varid( numbdyt, 'votemper' ) 
     358                  nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 
     359               ENDIF 
     360               WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd) 
     361               ipi = nblendta(igrd) 
     362               CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
     363 
     364               DO ib = 1, nblen(igrd) 
     365                  DO ik = 1, jpkm1 
     366                     tbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
     367                  END DO 
     368               END DO 
     369               ! 
     370               igrd = 1                                           ! salinity 
     371               IF( nblendta(igrd) .le. 0 ) THEN  
     372                  idvar = iom_varid( numbdyt, 'vosaline' ) 
     373                  nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 
     374               ENDIF 
     375               WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd) 
     376               ipi = nblendta(igrd) 
     377               CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
     378 
     379               DO ib = 1, nblen(igrd) 
     380                  DO ik = 1, jpkm1 
     381                     sbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
     382                  END DO 
     383               END DO 
     384            ENDIF  ! ln_bdy_tra_frs 
     385  
     386            IF(ln_bdy_dyn_frs) THEN 
     387 
     388               igrd = 2                                           ! u-velocity 
     389               IF ( nblendta(igrd) .le. 0 ) THEN  
     390                 idvar = iom_varid( numbdyu,'vozocrtx' ) 
     391                 nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar) 
     392               ENDIF 
     393               WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd) 
     394               ipi = nblendta(igrd) 
     395               CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
     396               DO ib = 1, nblen(igrd) 
     397                  DO ik = 1, jpkm1 
     398                     ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
     399                  END DO 
     400               END DO 
     401               ! 
     402               igrd = 3                                           ! v-velocity 
     403               IF ( nblendta(igrd) .le. 0 ) THEN  
     404                 idvar = iom_varid( numbdyv,'vomecrty' ) 
     405                 nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar) 
     406               ENDIF 
     407               WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd) 
     408               ipi = nblendta(igrd) 
     409               CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
     410               DO ib = 1, nblen(igrd) 
     411                  DO ik = 1, jpkm1 
     412                     vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
     413                  END DO 
     414               END DO 
     415            ENDIF ! ln_bdy_dyn_frs 
     416 
     417 
     418            IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN 
     419               ! First data time is after start of run 
     420               ! Put first value in both time levels 
     421               nbdy_b = nbdy_a 
     422               IF(ln_bdy_tra_frs) THEN 
     423                 tbdydta(:,:,1) = tbdydta(:,:,2) 
     424                 sbdydta(:,:,1) = sbdydta(:,:,2) 
     425               ENDIF 
     426               IF(ln_bdy_dyn_frs) THEN 
     427                 ubdydta(:,:,1) = ubdydta(:,:,2) 
     428                 vbdydta(:,:,1) = vbdydta(:,:,2) 
     429               ENDIF 
     430            END IF 
     431 
     432         END IF ! nbdy_dta == 0/1 
     433  
     434         ! In the case of constant boundary forcing fill bdy arrays once for all 
     435         IF ((ln_bdy_clim).AND.(ntimes_bdy==1)) THEN 
     436            IF(ln_bdy_tra_frs) THEN 
     437              tbdy  (:,:) = tbdydta  (:,:,2) 
     438              sbdy  (:,:) = sbdydta  (:,:,2) 
     439            ENDIF 
     440            IF(ln_bdy_dyn_frs) THEN 
     441              ubdy  (:,:) = ubdydta  (:,:,2) 
     442              vbdy  (:,:) = vbdydta  (:,:,2) 
     443            ENDIF 
     444 
     445            IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt ) 
     446            IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu ) 
     447            IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv ) 
     448         END IF 
     449 
     450      ENDIF                                            ! End if nit000 
     451 
     452 
     453      !                                                !---------------------! 
     454      !                                                !  at each time step  ! 
     455      !                                                !---------------------! 
     456 
     457      IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN  
     458         ! 
     459         ! Read one more record if necessary 
     460         !********************************** 
     461 
     462        IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN      ! remember that nbdy_b=0 for kt=nit000 
     463           nbdy_b = imois 
     464           nbdy_a = imois + 1 
     465           nbdy_b = MOD( nbdy_b, iman )   ;   IF( nbdy_b == 0 ) nbdy_b = iman 
     466           nbdy_a = MOD( nbdy_a, iman )   ;   IF( nbdy_a == 0 ) nbdy_a = iman 
     467           lect=.true. 
     468        ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN 
     469 
     470           IF ( nbdy_a < ntimes_bdy ) THEN 
     471              nbdy_b = nbdy_a 
     472              nbdy_a = nbdy_a + 1 
     473              lect  =.true. 
     474           ELSE 
     475              ! We have reached the end of the file 
     476              ! put the last data time into both time levels 
     477              nbdy_b = nbdy_a 
     478              IF(ln_bdy_tra_frs) THEN 
     479                tbdydta(:,:,1) =  tbdydta(:,:,2) 
     480                sbdydta(:,:,1) =  sbdydta(:,:,2) 
     481              ENDIF 
     482              IF(ln_bdy_dyn_frs) THEN 
     483                ubdydta(:,:,1) =  ubdydta(:,:,2) 
     484                vbdydta(:,:,1) =  vbdydta(:,:,2) 
     485              ENDIF 
     486            END IF ! nbdy_a < ntimes_bdy 
     487 
     488        END IF 
     489          
     490        IF( lect ) THEN 
     491           ! Swap arrays 
     492           IF(ln_bdy_tra_frs) THEN 
     493             tbdydta(:,:,1) =  tbdydta(:,:,2) 
     494             sbdydta(:,:,1) =  sbdydta(:,:,2) 
     495           ENDIF 
     496           IF(ln_bdy_dyn_frs) THEN 
     497             ubdydta(:,:,1) =  ubdydta(:,:,2) 
     498             vbdydta(:,:,1) =  vbdydta(:,:,2) 
     499           ENDIF 
     500  
     501           ! read another set 
     502           ipj  = 1 
     503           ipk  = jpk 
     504 
     505           IF(ln_bdy_tra_frs) THEN 
     506              !  
     507              igrd = 1                                   ! temperature 
     508              ipi  = nblendta(igrd) 
     509              CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
     510              DO ib = 1, nblen(igrd) 
     511                 DO ik = 1, jpkm1 
     512                    tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 
     513                 END DO 
     514              END DO 
     515              ! 
     516              igrd = 1                                   ! salinity 
     517              ipi  = nblendta(igrd) 
     518              CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 
     519              DO ib = 1, nblen(igrd) 
     520                 DO ik = 1, jpkm1 
     521                    sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 
     522                 END DO 
     523              END DO 
     524           ENDIF ! ln_bdy_tra_frs 
     525 
     526           IF(ln_bdy_dyn_frs) THEN 
     527              ! 
     528              igrd = 2                                   ! u-velocity 
     529              ipi  = nblendta(igrd) 
     530              CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
     531              DO ib = 1, nblen(igrd) 
     532                DO ik = 1, jpkm1 
     533                  ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
     534                END DO 
     535              END DO 
     536              ! 
     537              igrd = 3                                   ! v-velocity 
     538              ipi  = nblendta(igrd) 
     539              CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 
     540              DO ib = 1, nblen(igrd) 
     541                 DO ik = 1, jpkm1 
     542                    vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik) 
     543                 END DO 
     544              END DO 
     545           ENDIF ! ln_bdy_dyn_frs 
     546 
     547           ! 
     548           IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b 
     549           IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a ',nbdy_a 
     550           IF (.NOT.ln_bdy_clim) THEN 
     551              IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep(nbdy_b) 
     552              IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer 
     553              IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a) 
     554           ENDIF 
     555           ! 
     556       ENDIF ! end lect=.true. 
     557 
     558 
     559       ! Interpolate linearly 
     560       ! ******************** 
     561       !  
     562       IF( ln_bdy_clim ) THEN   ;   zxy = FLOAT( nday                  ) / FLOAT( nobis(nbdy_b) ) + 0.5 - i15 
     563       ELSE                     ;   zxy = FLOAT( istep(nbdy_b) - itimer ) / FLOAT( istep(nbdy_b) - istep(nbdy_a) ) 
     564       END IF 
     565 
     566          IF(ln_bdy_tra_frs) THEN 
     567             igrd = 1                                   ! temperature & salinity 
     568             DO ib = 1, nblen(igrd) 
     569               DO ik = 1, jpkm1 
     570                 tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1) 
     571                 sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1) 
     572               END DO 
     573             END DO 
     574          ENDIF 
     575 
     576          IF(ln_bdy_dyn_frs) THEN 
     577             igrd = 2                                   ! u-velocity 
     578             DO ib = 1, nblen(igrd) 
     579               DO ik = 1, jpkm1 
     580                 ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1)    
     581               END DO 
     582             END DO 
     583             ! 
     584             igrd = 3                                   ! v-velocity 
     585             DO ib = 1, nblen(igrd) 
     586               DO ik = 1, jpkm1 
     587                 vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1)    
     588               END DO 
     589             END DO 
     590          ENDIF 
     591 
     592      END IF                       !end if ((nbdy_dta==1).AND.(ntimes_bdy>1)) 
     593     
     594 
     595      !                                                !---------------------! 
     596      !                                                !     last call       ! 
     597      !                                                !---------------------! 
     598      IF( kt == nitend ) THEN 
     599          IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )              ! Closing of the 3 files 
     600          IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu ) 
     601          IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv ) 
     602      ENDIF 
     603      ! 
     604      ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs 
     605 
     606   END SUBROUTINE bdy_dta 
     607 
     608 
     609   SUBROUTINE bdy_dta_bt( kt, jit ) 
    69610      !!--------------------------------------------------------------------------- 
    70       !!                      ***  SUBROUTINE bdy_dta  *** 
     611      !!                      ***  SUBROUTINE bdy_dta_bt  *** 
    71612      !!                     
    72       !! ** Purpose :   Read unstructured boundary data 
     613      !! ** Purpose :   Read unstructured boundary data for Flather condition 
    73614      !! 
    74       !! ** Method  :    
     615      !! ** Method  :  At the first timestep, read in boundary data for two 
     616      !!               times from the file and time-interpolate. At other  
     617      !!               timesteps, check to see if we need another time from  
     618      !!               the file. If so read it in. Time interpolate. 
     619      !!--------------------------------------------------------------------------- 
     620!!gm DOCTOR names :   argument integer :  start with "k" 
     621      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
     622      INTEGER, INTENT( in ) ::   jit         ! barotropic time step index 
     623      !                                      ! (for timesplitting option, otherwise zero) 
    75624      !! 
    76       !! History :  OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original 
    77       !!                 "  ! 07-01 (D. Storkey) Update to use IOM module. 
    78       !!--------------------------------------------------------------------------- 
    79       !! * Arguments 
    80       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    81                                              ! (for timesplitting option, otherwise zero) 
    82  
    83       !! * Local declarations 
    84625      LOGICAL ::   lect                      ! flag for reading 
    85       INTEGER ::   jt, jb, jk, jgrd          ! dummy loop indices 
     626      INTEGER ::   it, ib, igrd              ! dummy loop indices 
    86627      INTEGER ::   idvar                     ! netcdf var ID 
    87628      INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing 
    88       INTEGER ::   kbdyt, kbdyu, kbdyv 
     629      INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 
    89630      INTEGER ::   itimer, totime 
    90631      INTEGER ::   ipi, ipj, ipk, inum       ! temporary integers (NetCDF read) 
    91632      INTEGER ::   iyear0, imonth0, iday0 
    92633      INTEGER ::   ihours0, iminutes0, isec0 
    93       INTEGER ::   kyear, kmonth, kday, ksecs 
     634      INTEGER ::   iyear, imonth, iday, isecs 
    94635      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files 
    95       REAL(wp) ::   dayfrac, zxy, offsett 
    96       REAL(wp) ::   offsetu, offsetv 
    97       REAL(wp) ::   dayjul0, zdayjulini 
    98       REAL(wp), DIMENSION(jpbtime)      ::   istepr             ! REAL time array from data files 
    99       REAL(wp), DIMENSION(jpbdta,1,jpk) ::   pdta3              ! temporary array for data fields 
    100       CHARACTER(LEN=80), DIMENSION(3)   ::   bdyfile 
     636      REAL(wp) ::   dayfrac, zxy, zoffsett 
     637      REAL(wp) ::   zoffsetu, zoffsetv 
     638      REAL(wp) ::   dayjul0, zdayjulini, zntotime 
     639      REAL(wp) ::   zinterval_s, zinterval_e                    ! First and last interval in time axis 
     640      REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files 
     641      REAL(wp), DIMENSION(jpbdta,1)     ::   zdta               ! temporary array for data fields 
     642      CHARACTER(LEN=80), DIMENSION(3)   ::   clfile 
    101643      CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate 
    102  
    103644      !!--------------------------------------------------------------------------- 
     645 
     646!!gm   add here the same style as in bdy_dta 
     647!!gm      clearly bdy_dta_bt and bdy_dta  can be combined...    
     648!!gm      too many things duplicated in the read of data...   simplification can be done 
    104649 
    105650      ! -------------------- ! 
     
    111656      ! Some time variables for monthly climatological forcing: 
    112657      ! ******************************************************* 
     658 !!gm  here  use directely daymod variables 
    113659  
    114660      iman  = INT( raamo ) ! Number of months in a year 
     
    123669      ! ********************************************* 
    124670 
    125       itimer = (kt-nit000+1)*rdt      ! current time in seconds for interpolation  
    126  
    127       ! -------------------- ! 
    128       ! 1.   First call only ! 
    129       ! -------------------- ! 
    130  
    131       IF ( kt == nit000 ) THEN 
    132  
    133         istep(:) = 0 
    134  
    135         nbdy1=0 
    136         nbdy2=0 
    137  
    138       ! 1.1  Get time information from bdy data file 
    139       ! ******************************************** 
     671      itimer = ((kt-1)-nit000+1)*rdt      ! current time in seconds for interpolation  
     672      itimer = itimer + jit*rdtbt         ! in non-climatological case 
     673 
     674      IF ( ln_bdy_tides ) THEN 
     675 
     676         ! -------------------------------------! 
     677         ! Update BDY fields with tidal forcing ! 
     678         ! -------------------------------------!   
     679 
     680         CALL tide_update( kt, jit )  
     681   
     682      ENDIF 
     683 
     684      IF ( ln_bdy_dyn_fla ) THEN 
     685 
     686         ! -------------------------------------! 
     687         ! Update BDY fields with model data    ! 
     688         ! -------------------------------------!   
     689 
     690      !                                                !-------------------! 
     691      IF( kt == nit000 ) THEN                          !  First call only  ! 
     692         !                                             !-------------------! 
     693         istep_bt(:) = 0 
     694         nbdy_b_bt    = 0 
     695         nbdy_a_bt    = 0 
     696 
     697         ! Get time information from bdy data file 
     698         ! *************************************** 
    140699 
    141700        IF(lwp) WRITE(numout,*) 
    142         IF(lwp) WRITE(numout,*)    'bdy_dta :Initialize unstructured boundary data' 
     701        IF(lwp) WRITE(numout,*)    'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.' 
    143702        IF(lwp) WRITE(numout,*)    '~~~~~~~'  
    144703 
    145         IF (nbdy_dta == 0) THEN 
     704        IF( nbdy_dta == 0 ) THEN 
    146705          IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions' 
    147706 
     
    154713                                           ! necessary time dumps in file are included 
    155714 
    156           bdyfile(1) = filbdy_data_T 
    157           bdyfile(2) = filbdy_data_U 
    158           bdyfile(3) = filbdy_data_V 
    159  
    160           DO jgrd = 1,3 
    161  
    162             CALL iom_open( bdyfile(jgrd), inum ) 
    163             CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits )  
    164             CALL iom_close( inum ) 
     715          clfile(1) = filbdy_data_bt_T 
     716          clfile(2) = filbdy_data_bt_U 
     717          clfile(3) = filbdy_data_bt_V 
     718 
     719          DO igrd = 1,3 
     720 
     721            CALL iom_open( clfile(igrd), inum ) 
     722            CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits )  
     723 
     724            SELECT CASE( igrd ) 
     725               CASE (1)  
     726                  numbdyt = inum 
     727               CASE (2)  
     728                  numbdyu = inum 
     729               CASE (3)  
     730                  numbdyv = inum 
     731            END SELECT 
    165732 
    166733            ! Calculate time offset  
     
    170737            CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 
    171738            ! Compute model initialization time  
    172             kyear  = ndastp / 10000 
    173             kmonth = ( ndastp - kyear * 10000 ) / 100 
    174             kday   = ndastp - kyear * 10000 - kmonth * 100 
    175             ksecs  = dayfrac * 86400 
    176             CALL ymds2ju(kyear, kmonth, kday, real(ksecs) , zdayjulini) 
    177             ! offset from initialization date: 
    178             offset = (dayjul0-zdayjulini)*86400 
     739            iyear  = ndastp / 10000 
     740            imonth = ( ndastp - iyear * 10000 ) / 100 
     741            iday   = ndastp - iyear * 10000 - imonth * 100 
     742            isecs  = dayfrac * 86400 
     743            CALL ymds2ju(iyear, imonth, iday, real(isecs) , zdayjulini) 
     744            ! zoffset from initialization date: 
     745            zoffset = (dayjul0-zdayjulini)*86400 
    179746            ! 
    180747 
    1817487000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 
    182  
    183749 
    184750            !! TO BE DONE... Check consistency between calendar from file  
     
    186752            !! when calendar in model available outside of IOIPSL. 
    187753 
    188             IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy 
    189             IF(lwp) WRITE(numout,*) 'offset: ',offset 
    190             IF(lwp) WRITE(numout,*) 'totime: ',totime 
    191             IF(lwp) WRITE(numout,*) 'istepr: ',istepr 
    192  
    193754            ! Check that there are not too many times in the file.  
    194             IF (kbdy > jpbtime) THEN 
    195               IF(lwp) WRITE(numout,*) 
    196               IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter' 
    197               IF(lwp) WRITE(numout,*) ' ==========  jpbtime= ',jpbtime,' kbdy= ', kbdy              
    198               IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd) 
    199               nstop = nstop + 1 
     755            IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( & 
     756                 'Number of time dumps in bdy file exceed jpbtime parameter', & 
     757                 'Check file:' // TRIM(clfile(igrd))  ) 
     758 
     759            ! Check that time array increases (or interp will fail): 
     760            DO it = 2, ntimes_bdy 
     761               IF ( zstepr(it-1) >= zstepr(it) ) THEN 
     762                  CALL ctl_stop('Time array in unstructured boundary data file', & 
     763                       'does not continuously increase.',               & 
     764                       'Check file:' // TRIM(clfile(igrd))  ) 
     765                  EXIT 
     766               END IF 
     767            END DO 
     768 
     769            IF ( .NOT. ln_bdy_clim ) THEN 
     770               ! Check that times in file span model run time: 
     771 
     772               ! Note: the fields may be time means, so we allow nit000 to be before 
     773               ! first time in the file, provided that it falls inside the meaning 
     774               ! period of the first field.  Until we can get the meaning period 
     775               ! from the file, use the interval between fields as a proxy. 
     776               ! If nit000 is before the first time, use the value at first time 
     777               ! instead of extrapolating.  This is done by putting time 1 into 
     778               ! both time levels. 
     779               ! The same applies to the last time level: see setting of lect below. 
     780 
     781               IF ( ntimes_bdy == 1 ) CALL ctl_stop( & 
     782                    'There is only one time dump in data files', & 
     783                    'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' ) 
     784 
     785               zinterval_s = zstepr(2) - zstepr(1) 
     786               zinterval_e = zstepr(ntimes_bdy) - zstepr(ntimes_bdy-1) 
     787 
     788               IF ( zstepr(1) - zinterval_s / 2.0 > 0 ) THEN              
     789                  IF(lwp) WRITE(numout,*) 'First bdy time relative to nit000:', zstepr(1) 
     790                  IF(lwp) WRITE(numout,*) 'Interval between first two times: ', zinterval_s 
     791                  CALL ctl_stop( 'First data time is after start of run', &  
     792                       'by more than half a meaning period', & 
     793                       'Check file: ' // TRIM(clfile(igrd)) ) 
     794               END IF 
     795 
     796               IF ( zstepr(ntimes_bdy) + zinterval_e / 2.0 < zntotime ) THEN 
     797                  IF(lwp) WRITE(numout,*) 'Last bdy time relative to nit000:', zstepr(ntimes_bdy) 
     798                  IF(lwp) WRITE(numout,*) 'Interval between last two times: ', zinterval_e 
     799                  CALL ctl_stop( 'Last data time is before end of run', &  
     800                       'by more than half a meaning period', & 
     801                       'Check file: ' // TRIM(clfile(igrd))  ) 
     802               END IF 
     803 
     804            END IF ! .NOT. ln_bdy_clim 
     805 
     806            IF ( igrd .EQ. 1) THEN 
     807              ntimes_bdyt = ntimes_bdy_bt 
     808              zoffsett = zoffset 
     809              istept(:) = INT( zstepr(:) + zoffset ) 
     810            ELSE IF (igrd .EQ. 2) THEN 
     811              ntimes_bdyu = ntimes_bdy_bt 
     812              zoffsetu = zoffset 
     813              istepu(:) = INT( zstepr(:) + zoffset ) 
     814            ELSE IF (igrd .EQ. 3) THEN 
     815              ntimes_bdyv = ntimes_bdy_bt 
     816              zoffsetv = zoffset 
     817              istepv(:) = INT( zstepr(:) + zoffset ) 
    200818            ENDIF 
    201819 
    202             ! Check that time array increases: 
    203          
    204             jt=1 
    205             DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy-1))) 
    206               jt=jt+1 
    207             END DO 
    208  
    209             IF ((jt.NE.kbdy).AND.(kbdy.GT.1)) THEN 
    210               IF(lwp) WRITE(numout,*) 
    211               IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files' 
    212               IF(lwp) WRITE(numout,*) ' ==========  does not continuously increase. Check file: '          
    213               IF(lwp) WRITE(numout,*) bdyfile(jgrd) 
    214               IF(lwp) WRITE(numout,*) 
    215               nstop = nstop + 1 
    216             ENDIF 
    217  
    218             ! Check that times in file span model run time: 
    219   
    220             IF (istepr(1)+offset > 0) THEN 
    221               IF(lwp) WRITE(numout,*) 
    222               IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time' 
    223               IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd) 
    224               IF(lwp) WRITE(numout,*) 
    225               nstop = nstop + 1 
    226             END IF 
    227  
    228             IF (istepr(kbdy)+offset < totime) THEN 
    229               IF(lwp) WRITE(numout,*) 
    230               IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time' 
    231               IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd)  
    232               IF(lwp) WRITE(numout,*) 
    233               nstop = nstop + 1 
    234             END IF 
    235  
    236             IF ( jgrd .EQ. 1) THEN 
    237               kbdyt = kbdy 
    238               offsett = offset 
    239               istept(:) = INT( istepr(:) + offset ) 
    240             ELSE IF (jgrd .EQ. 2) THEN 
    241               kbdyu = kbdy 
    242               offsetu = offset 
    243               istepu(:) = INT( istepr(:) + offset ) 
    244             ELSE IF (jgrd .EQ. 3) THEN 
    245               kbdyv = kbdy 
    246               offsetv = offset 
    247               istepv(:) = INT( istepr(:) + offset ) 
    248             ENDIF 
    249  
    250820          ENDDO 
    251821 
    252822      ! Verify time consistency between files   
    253823 
    254           IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN 
    255             IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps' 
    256             IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 
    257             nstop = nstop + 1 
     824          IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 
     825             CALL ctl_stop( & 
     826             'Time axis lengths differ between bdy data files', & 
     827             'Multiple time frequencies not implemented yet' ) 
     828          ELSE 
     829            ntimes_bdy_bt = ntimes_bdyt 
    258830          ENDIF 
    259           kbdy = kbdyt 
    260  
    261           IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN 
    262             IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin' 
    263             IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 
    264             nstop = nstop + 1 
     831 
     832          IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN 
     833            CALL ctl_stop( &  
     834            'Bdy data files must have the same time origin', & 
     835            'Multiple time frequencies not implemented yet'  ) 
    265836          ENDIF 
    266           offset = offsett 
     837          zoffset = zoffsett 
    267838 
    268839      !! Check that times are the same in the three files... HERE. 
    269           istep(:) = istept(:) 
     840          istep_bt(:) = istept(:) 
    270841 
    271842      ! Check number of time dumps:               
    272           IF ((kbdy.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN 
    273             IF(lwp) WRITE(numout,*) 
    274             IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files' 
    275             IF(lwp) WRITE(numout,*) ' ==========  Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.'             
    276             IF(lwp) WRITE(numout,*) 
    277             nstop = nstop + 1 
    278           ENDIF 
    279  
    280843          IF (ln_bdy_clim) THEN 
    281             IF ((kbdy.NE.1).AND.(kbdy.NE.12)) THEN 
    282               IF(lwp) WRITE(numout,*) 
    283               IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),' 
    284               IF(lwp) WRITE(numout,*) ' ==========  bdy data files must contain 1 or 12 time dumps.'          
    285               IF(lwp) WRITE(numout,*) 
    286               nstop = nstop + 1 
    287             ELSEIF (kbdy.EQ.1 ) THEN 
     844            SELECT CASE ( ntimes_bdy_bt ) 
     845            CASE( 1 ) 
    288846              IF(lwp) WRITE(numout,*) 
    289847              IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 
    290848              IF(lwp) WRITE(numout,*)              
    291             ELSEIF (kbdy.EQ.12) THEN 
     849            CASE( 12 ) 
    292850              IF(lwp) WRITE(numout,*) 
    293851              IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 
    294852              IF(lwp) WRITE(numout,*)  
    295             ENDIF 
     853            CASE DEFAULT 
     854              CALL ctl_stop( & 
     855                'For climatological boundary forcing (ln_bdy_clim=.true.),',& 
     856                'bdy data files must contain 1 or 12 time dumps.' ) 
     857            END SELECT 
    296858          ENDIF 
    297859 
    298860      ! Find index of first record to read (before first model time).  
    299861 
    300           jt=1 
    301           DO WHILE ( ((istep(jt+1)) <= 0 ).AND.(jt.LE.(kbdy-1))) 
    302             jt=jt+1 
    303           END DO 
    304           nbdy1 = jt 
    305  
    306           WRITE(numout,*) 'Time offset is ',offset 
    307           WRITE(numout,*) 'First record to read is ',nbdy1 
     862          it=1 
     863          DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1))) 
     864            it=it+1 
     865          END DO 
     866          nbdy_b_bt = it 
     867 
     868          WRITE(numout,*) 'Time offset is ',zoffset 
     869          WRITE(numout,*) 'First record to read is ',nbdy_b_bt 
    308870 
    309871        ENDIF ! endif (nbdy_dta == 1) 
     
    314876        IF ( nbdy_dta == 0) THEN 
    315877          ! boundary data arrays are filled with initial conditions 
    316           DO jk = 1, jpkm1  
    317  
    318 #if ! defined key_barotropic 
    319             jgrd = 1            ! T-points data  
    320             DO jb = 1, nblen(jgrd)               
    321               tbdy(jb,jk) = tn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    322               sbdy(jb,jk) = sn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    323             END DO 
    324 #endif 
    325  
    326             jgrd = 2            ! U-points data  
    327             DO jb = 1, nblen(jgrd)               
    328               ubdy(jb,jk) = un(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    329             END DO 
    330  
    331             jgrd = 3            ! V-points data  
    332             DO jb = 1, nblen(jgrd)               
    333               vbdy(jb,jk) = vn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    334             END DO 
    335           END DO  
     878          igrd = 2            ! U-points data  
     879          DO ib = 1, nblen(igrd)               
     880            ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1) 
     881          END DO 
     882 
     883          igrd = 3            ! V-points data  
     884          DO ib = 1, nblen(igrd)               
     885            vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1) 
     886          END DO 
     887 
     888          igrd = 1            ! T-points data  
     889          DO ib = 1, nblen(igrd)               
     890            sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd)) 
     891          END DO 
    336892 
    337893        ELSEIF (nbdy_dta == 1) THEN 
    338894  
    339895        ! Set first record in the climatological case:    
    340           IF ((ln_bdy_clim).AND.(kbdy==1)) THEN 
    341             nbdy2 = 1 
    342           ELSEIF ((ln_bdy_clim).AND.(kbdy==iman)) THEN 
    343             nbdy1 = 0 
    344             nbdy2 = imois 
     896          IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN 
     897            nbdy_a_bt = 1 
     898          ELSEIF ((ln_bdy_clim).AND.(ntimes_bdy_bt==iman)) THEN 
     899            nbdy_b_bt = 0 
     900            nbdy_a_bt = imois 
    345901          ELSE 
    346             nbdy2 = nbdy1 
     902            nbdy_a_bt = nbdy_b_bt 
    347903          END IF 
    348904  
    349905         ! Open Netcdf files: 
    350906 
    351           CALL iom_open ( filbdy_data_T, bdynumt ) 
    352           CALL iom_open ( filbdy_data_U, bdynumu ) 
    353           CALL iom_open ( filbdy_data_V, bdynumv ) 
     907          CALL iom_open ( filbdy_data_bt_T, numbdyt_bt ) 
     908          CALL iom_open ( filbdy_data_bt_U, numbdyu_bt ) 
     909          CALL iom_open ( filbdy_data_bt_V, numbdyv_bt ) 
    354910 
    355911         ! Read first record: 
    356912          ipj=1 
    357           ipk=jpk 
    358           jgrd=1 
    359           ipi=nblendta(jgrd) 
    360  
    361 #if ! defined key_barotropic 
    362           !temperature 
    363  
    364           jgrd=1 
    365           IF ( nblendta(jgrd) .le. 0 ) THEN  
    366             idvar = iom_varid( bdynumt,'votemper' ) 
    367             nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar) 
     913          igrd=1 
     914          ipi=nblendta(igrd) 
     915 
     916          ! ssh 
     917          igrd=1 
     918          IF ( nblendta(igrd) .le. 0 ) THEN  
     919            idvar = iom_varid( numbdyt_bt,'sossheig' ) 
     920            nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar) 
    368921          ENDIF 
    369           WRITE(numout,*) 'Dim size for votemper is ',nblendta(jgrd) 
    370           ipi=nblendta(jgrd) 
    371  
    372           CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    373  
    374           DO jb=1, nblen(jgrd) 
    375             DO jk=1,jpkm1 
    376               tbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    377             END DO 
    378           END DO 
    379  
    380           ! salinity 
    381  
    382           jgrd=1 
    383           IF ( nblendta(jgrd) .le. 0 ) THEN  
    384             idvar = iom_varid( bdynumt,'vosaline' ) 
    385             nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar) 
     922          WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd) 
     923          ipi=nblendta(igrd) 
     924 
     925          CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
     926 
     927          DO ib=1, nblen(igrd) 
     928            sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
     929          END DO 
     930  
     931          ! u-velocity 
     932          igrd=2 
     933          IF ( nblendta(igrd) .le. 0 ) THEN  
     934            idvar = iom_varid( numbdyu_bt,'vobtcrtx' ) 
     935            nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar) 
    386936          ENDIF 
    387           WRITE(numout,*) 'Dim size for vosaline is ',nblendta(jgrd) 
    388           ipi=nblendta(jgrd) 
    389  
    390           CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    391  
    392           DO jb=1, nblen(jgrd) 
    393             DO jk=1,jpkm1 
    394               sbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    395             END DO 
    396           END DO 
    397 #endif 
    398            
    399           ! u-velocity 
    400  
    401           jgrd=2 
    402           IF ( nblendta(jgrd) .le. 0 ) THEN  
    403             idvar = iom_varid( bdynumu,'vozocrtx' ) 
    404             nblendta(jgrd) = iom_file(bdynumu)%dimsz(1,idvar) 
     937          WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd) 
     938          ipi=nblendta(igrd) 
     939 
     940          CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
     941 
     942          DO ib=1, nblen(igrd) 
     943            ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
     944          END DO 
     945 
     946          ! v-velocity 
     947          igrd=3 
     948          IF ( nblendta(igrd) .le. 0 ) THEN  
     949            idvar = iom_varid( numbdyv_bt,'vobtcrty' ) 
     950            nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar) 
    405951          ENDIF 
    406           WRITE(numout,*) 'Dim size for vozocrtx is ',nblendta(jgrd) 
    407           ipi=nblendta(jgrd) 
    408  
    409           CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    410  
    411           DO jb=1, nblen(jgrd) 
    412             DO jk=1,jpkm1 
    413               ubdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    414             END DO 
    415           END DO 
    416  
    417           ! v-velocity 
    418           jgrd=3 
    419           IF ( nblendta(jgrd) .le. 0 ) THEN  
    420             idvar = iom_varid( bdynumv,'vomecrty' ) 
    421             nblendta(jgrd) = iom_file(bdynumv)%dimsz(1,idvar) 
    422           ENDIF 
    423           WRITE(numout,*) 'Dim size for vomecrty is ',nblendta(jgrd) 
    424           ipi=nblendta(jgrd) 
    425  
    426           CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    427  
    428           DO jb=1, nblen(jgrd) 
    429             DO jk=1,jpkm1 
    430               vbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    431             END DO 
     952          WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd) 
     953          ipi=nblendta(igrd) 
     954 
     955          CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
     956 
     957          DO ib=1, nblen(igrd) 
     958            vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    432959          END DO 
    433960 
     
    435962  
    436963        ! In the case of constant boundary forcing fill bdy arrays once for all 
    437         IF ((ln_bdy_clim).AND.(kbdy==1)) THEN 
    438 #if ! defined key_barotropic 
    439           tbdy  (:,:) = tbdydta  (:,:,2) 
    440           sbdy  (:,:) = sbdydta  (:,:,2) 
    441 #endif 
    442           ubdy  (:,:) = ubdydta  (:,:,2) 
    443           vbdy  (:,:) = vbdydta  (:,:,2) 
    444  
    445           CALL iom_close( bdynumt ) 
    446           CALL iom_close( bdynumu ) 
    447           CALL iom_close( bdynumv ) 
     964        IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN 
     965 
     966          ubtbdy  (:) = ubtbdydta  (:,2) 
     967          vbtbdy  (:) = vbtbdydta  (:,2) 
     968          sshbdy  (:) = sshbdydta  (:,2) 
     969 
     970          CALL iom_close( numbdyt_bt ) 
     971          CALL iom_close( numbdyu_bt ) 
     972          CALL iom_close( numbdyv_bt ) 
     973 
    448974        END IF 
    449975 
     
    454980      ! -------------------- ! 
    455981 
    456       IF ((nbdy_dta==1).AND.(kbdy>1)) THEN  
     982      IF ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) THEN  
    457983 
    458984      ! 2.1 Read one more record if necessary 
    459985      !************************************** 
    460986 
    461         IF ( (ln_bdy_clim).AND.(imois/=nbdy1) ) THEN ! remember that nbdy1=0 for kt=nit000 
    462          nbdy1 = imois 
    463          nbdy2 = imois+1 
    464          nbdy1 = MOD( nbdy1, iman ) 
    465          IF( nbdy1 == 0 ) nbdy1 = iman 
    466          nbdy2 = MOD( nbdy2, iman ) 
    467          IF( nbdy2 == 0 ) nbdy2 = iman 
     987        IF ( (ln_bdy_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000 
     988         nbdy_b_bt = imois 
     989         nbdy_a_bt = imois+1 
     990         nbdy_b_bt = MOD( nbdy_b_bt, iman ) 
     991         IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman 
     992         nbdy_a_bt = MOD( nbdy_a_bt, iman ) 
     993         IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman 
    468994         lect=.true. 
    469995 
    470         ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep(nbdy2))) THEN 
    471           nbdy1=nbdy2 
    472           nbdy2=nbdy2+1 
    473           lect=.true. 
    474         END IF 
    475           
    476         IF (lect) THEN 
    477  
    478         ! Swap arrays 
    479 #if ! defined key_barotropic 
    480           tbdydta(:,:,1) =  tbdydta(:,:,2) 
    481           sbdydta(:,:,1) =  sbdydta(:,:,2) 
    482 #endif 
    483           ubdydta(:,:,1) =  ubdydta(:,:,2) 
    484           vbdydta(:,:,1) =  vbdydta(:,:,2) 
    485   
    486         ! read another set 
    487  
    488           ipj=1 
    489           ipk=jpk 
    490           jgrd=1 
    491           ipi=nblendta(jgrd) 
    492  
    493 #if ! defined key_barotropic 
    494           !temperature 
    495   
    496           jgrd=1 
    497           ipi=nblendta(jgrd) 
    498  
    499           CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    500  
    501           DO jb=1, nblen(jgrd) 
    502             DO jk=1,jpkm1 
    503               tbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    504             END DO 
    505           END DO 
    506  
    507           ! salinity 
    508  
    509           jgrd=1 
    510           ipi=nblendta(jgrd) 
    511  
    512           CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    513  
    514           DO jb=1, nblen(jgrd) 
    515             DO jk=1,jpkm1 
    516               sbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    517             END DO 
    518           END DO 
    519 #endif 
    520            
    521           ! u-velocity 
    522  
    523           jgrd=2 
    524           ipi=nblendta(jgrd) 
    525  
    526           CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    527  
    528           DO jb=1, nblen(jgrd) 
    529             DO jk=1,jpkm1 
    530               ubdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    531             END DO 
    532           END DO 
    533  
    534           ! v-velocity 
    535           jgrd=3 
    536           ipi=nblendta(jgrd) 
    537  
    538           CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 
    539  
    540           DO jb=1, nblen(jgrd) 
    541             DO jk=1,jpkm1 
    542               vbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk) 
    543             END DO 
    544           END DO 
    545  
    546          IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy1 ',nbdy1 
    547          IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy2 ',nbdy2 
    548          IF (.NOT.ln_bdy_clim) THEN 
    549            IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep(nbdy1) 
    550            IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer 
    551            IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy2) 
    552          ENDIF 
    553         END IF ! end lect=.true. 
    554  
    555  
    556       ! 2.2   Interpolate linearly: 
    557       ! *************************** 
    558      
    559         IF (ln_bdy_clim) THEN 
    560           zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1) ) + 0.5 - i15 
    561         ELSE           
    562           zxy = FLOAT(istep(nbdy1)-itimer) / FLOAT(istep(nbdy1)-istep(nbdy2)) 
    563         END IF 
    564  
    565 #if ! defined key_barotropic 
    566           jgrd=1 
    567           DO jb=1, nblen(jgrd) 
    568             DO jk=1, jpkm1 
    569               tbdy(jb,jk) = zxy      * tbdydta(jb,jk,2) + & 
    570                             (1.-zxy) * tbdydta(jb,jk,1) 
    571               sbdy(jb,jk) = zxy      * sbdydta(jb,jk,2) + & 
    572                             (1.-zxy) * sbdydta(jb,jk,1) 
    573             END DO 
    574           END DO 
    575 #endif 
    576  
    577           jgrd=2 
    578           DO jb=1, nblen(jgrd) 
    579             DO jk=1, jpkm1 
    580               ubdy(jb,jk) = zxy      * ubdydta(jb,jk,2) + & 
    581                             (1.-zxy) * ubdydta(jb,jk,1)    
    582             END DO 
    583           END DO 
    584  
    585           jgrd=3 
    586           DO jb=1, nblen(jgrd) 
    587             DO jk=1, jpkm1 
    588               vbdy(jb,jk) = zxy      * vbdydta(jb,jk,2) + & 
    589                             (1.-zxy) * vbdydta(jb,jk,1)    
    590             END DO 
    591           END DO 
    592  
    593       END IF !end if ((nbdy_dta==1).AND.(kbdy>1)) 
    594      
    595       ! ------------------- ! 
    596       ! Last call kt=nitend ! 
    597       ! ------------------- ! 
    598  
    599       ! Closing of the 3 files 
    600       IF( kt == nitend ) THEN 
    601           CALL iom_close( bdynumt ) 
    602           CALL iom_close( bdynumu ) 
    603           CALL iom_close( bdynumv ) 
    604       ENDIF 
    605  
    606    END SUBROUTINE bdy_dta 
    607  
    608    SUBROUTINE bdy_dta_bt( kt, jit ) 
    609       !!--------------------------------------------------------------------------- 
    610       !!                      ***  SUBROUTINE bdy_dta_bt  *** 
    611       !!                     
    612       !! ** Purpose :   Read unstructured boundary data for barotropic variables 
    613       !! 
    614       !! ** Method  :    
    615       !! 
    616       !! History : OPA 9.0  ! 07-07 (D. Storkey) Original 
    617       !!--------------------------------------------------------------------------- 
    618       !! * Arguments 
    619       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    620       INTEGER, INTENT( in ) ::   jit         ! barotropic time step index 
    621                                              ! (for timesplitting option, otherwise zero) 
    622  
    623       !! * Local declarations 
    624       LOGICAL ::   lect                      ! flag for reading 
    625       INTEGER ::   jt, jb, jgrd              ! dummy loop indices 
    626       INTEGER ::   idvar                     ! netcdf var ID 
    627       INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing 
    628       INTEGER ::   kbdyt, kbdyu, kbdyv 
    629       INTEGER ::   itimer, totime 
    630       INTEGER ::   ipi, ipj, ipk, inum       ! temporary integers (NetCDF read) 
    631       INTEGER ::   iyear0, imonth0, iday0 
    632       INTEGER ::   ihours0, iminutes0, isec0 
    633       INTEGER ::   kyear, kmonth, kday, ksecs 
    634       INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files 
    635       REAL(wp) ::   dayfrac, zxy, offsett 
    636       REAL(wp) ::   offsetu, offsetv 
    637       REAL(wp) ::   dayjul0, zdayjulini 
    638       REAL(wp), DIMENSION(jpbtime)      ::   istepr             ! REAL time array from data files 
    639       REAL(wp), DIMENSION(jpbdta,1)     ::   pdta2              ! temporary array for data fields 
    640       REAL(wp), DIMENSION(jpbdta,1,jpk) ::   pdta3              ! temporary array for data fields 
    641       CHARACTER(LEN=80), DIMENSION(3)   ::   bdyfile 
    642       CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate 
    643  
    644       !!--------------------------------------------------------------------------- 
    645  
    646       ! -------------------- ! 
    647       !    Initialization    ! 
    648       ! -------------------- ! 
    649  
    650       lect   = .false.           ! If true, read a time record 
    651  
    652       ! Some time variables for monthly climatological forcing: 
    653       ! ******************************************************* 
    654   
    655       iman  = INT( raamo ) ! Number of months in a year 
    656  
    657       i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 
    658       ! i15=0 if the current day is in the first half of the month, else i15=1 
    659  
    660       imois = nmonth + i15 - 1            ! imois is the first month record 
    661       IF( imois == 0 ) imois = iman 
    662  
    663       ! Time variable for non-climatological forcing: 
    664       ! ********************************************* 
    665  
    666       itimer = ((kt-1)-nit000+1)*rdt      ! current time in seconds for interpolation  
    667       itimer = itimer + jit*rdtbt         ! in non-climatological case 
    668  
    669       ! -------------------------------------! 
    670       ! Update BDY fields with tidal forcing ! 
    671       ! -------------------------------------!   
    672  
    673       IF ( lk_bdy_tides ) CALL tide_update( kt, jit )  
    674    
    675       IF ( lk_bdy ) THEN  
    676  
    677       ! -------------------- ! 
    678       ! 1.   First call only ! 
    679       ! -------------------- ! 
    680  
    681       IF ( kt == nit000 ) THEN 
    682  
    683         istep_bt(:) = 0 
    684  
    685         nbdy1_bt=0 
    686         nbdy2_bt=0 
    687  
    688       ! 1.1  Get time information from bdy data file 
    689       ! ******************************************** 
    690  
    691         IF(lwp) WRITE(numout,*) 
    692         IF(lwp) WRITE(numout,*)    'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.' 
    693         IF(lwp) WRITE(numout,*)    '~~~~~~~'  
    694  
    695         IF (nbdy_dta == 0) THEN 
    696           IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions' 
    697  
    698         ELSEIF (nbdy_dta == 1) THEN 
    699           IF(lwp) WRITE(numout,*)  'Bdy data are read in netcdf files' 
    700  
    701           dayfrac = adatrj-float(itimer)/86400. ! day fraction at time step kt-1 
    702           dayfrac = dayfrac - INT(dayfrac)      ! 
    703           totime = (nitend-nit000+1)*rdt ! Total time of the run to verify that all the 
    704                                            ! necessary time dumps in file are included 
    705  
    706           bdyfile(1) = filbdy_data_bt_T 
    707           bdyfile(2) = filbdy_data_bt_U 
    708           bdyfile(3) = filbdy_data_bt_V 
    709  
    710           DO jgrd = 1,3 
    711  
    712             CALL iom_open( bdyfile(jgrd), inum ) 
    713             CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits )  
    714             CALL iom_close( inum ) 
    715  
    716             ! Calculate time offset  
    717             READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 
    718             ! Convert time origin in file to julian days  
    719             isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 
    720             CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 
    721             ! Compute model initialization time  
    722             kyear  = ndastp / 10000 
    723             kmonth = ( ndastp - kyear * 10000 ) / 100 
    724             kday   = ndastp - kyear * 10000 - kmonth * 100 
    725             ksecs  = dayfrac * 86400 
    726             CALL ymds2ju(kyear, kmonth, kday, real(ksecs) , zdayjulini) 
    727             ! offset from initialization date: 
    728             offset = (dayjul0-zdayjulini)*86400 
    729             ! 
    730  
    731 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 
    732  
    733             !! TO BE DONE... Check consistency between calendar from file  
    734             !! (available optionally from iom_gettime) and calendar in model  
    735             !! when calendar in model available outside of IOIPSL. 
    736  
    737             IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy_bt 
    738             IF(lwp) WRITE(numout,*) 'offset: ',offset 
    739             IF(lwp) WRITE(numout,*) 'totime: ',totime 
    740             IF(lwp) WRITE(numout,*) 'istepr: ',istepr 
    741  
    742             ! Check that there are not too many times in the file.  
    743             IF (kbdy_bt > jpbtime) THEN 
    744               IF(lwp) WRITE(numout,*) 
    745               IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter' 
    746               IF(lwp) WRITE(numout,*) ' ==========  jpbtime= ',jpbtime,' kbdy_bt= ', kbdy_bt              
    747               IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd) 
    748               nstop = nstop + 1 
    749             ENDIF 
    750  
    751             ! Check that time array increases: 
    752          
    753             jt=1 
    754             DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy_bt-1))) 
    755               jt=jt+1 
    756             END DO 
    757  
    758             IF ((jt.NE.kbdy_bt).AND.(kbdy_bt.GT.1)) THEN 
    759               IF(lwp) WRITE(numout,*) 
    760               IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files' 
    761               IF(lwp) WRITE(numout,*) ' ==========  does not continuously increase. Check file: '          
    762               IF(lwp) WRITE(numout,*) bdyfile(jgrd) 
    763               IF(lwp) WRITE(numout,*) 
    764               nstop = nstop + 1 
    765             ENDIF 
    766  
    767             ! Check that times in file span model run time: 
    768   
    769             IF (istepr(1)+offset > 0) THEN 
    770               IF(lwp) WRITE(numout,*) 
    771               IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time' 
    772               IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd) 
    773               IF(lwp) WRITE(numout,*) 
    774               nstop = nstop + 1 
    775             END IF 
    776  
    777             IF (istepr(kbdy_bt)+offset < totime) THEN 
    778               IF(lwp) WRITE(numout,*) 
    779               IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time' 
    780               IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd)  
    781               IF(lwp) WRITE(numout,*) 
    782               nstop = nstop + 1 
    783             END IF 
    784  
    785             IF ( jgrd .EQ. 1) THEN 
    786               kbdyt = kbdy_bt 
    787               offsett = offset 
    788               istept(:) = INT( istepr(:) + offset ) 
    789             ELSE IF (jgrd .EQ. 2) THEN 
    790               kbdyu = kbdy_bt 
    791               offsetu = offset 
    792               istepu(:) = INT( istepr(:) + offset ) 
    793             ELSE IF (jgrd .EQ. 3) THEN 
    794               kbdyv = kbdy_bt 
    795               offsetv = offset 
    796               istepv(:) = INT( istepr(:) + offset ) 
    797             ENDIF 
    798  
    799           ENDDO 
    800  
    801       ! Verify time consistency between files   
    802  
    803           IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN 
    804             IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps' 
    805             IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 
    806             nstop = nstop + 1 
    807           ENDIF 
    808           kbdy_bt = kbdyt 
    809  
    810           IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN 
    811             IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin' 
    812             IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 
    813             nstop = nstop + 1 
    814           ENDIF 
    815           offset = offsett 
    816  
    817       !! Check that times are the same in the three files... HERE. 
    818           istep_bt(:) = istept(:) 
    819  
    820       ! Check number of time dumps:               
    821           IF ((kbdy_bt.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN 
    822             IF(lwp) WRITE(numout,*) 
    823             IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files' 
    824             IF(lwp) WRITE(numout,*) ' ==========  Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.'             
    825             IF(lwp) WRITE(numout,*) 
    826             nstop = nstop + 1 
    827           ENDIF 
    828  
    829           IF (ln_bdy_clim) THEN 
    830             IF ((kbdy_bt.NE.1).AND.(kbdy_bt.NE.12)) THEN 
    831               IF(lwp) WRITE(numout,*) 
    832               IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),' 
    833               IF(lwp) WRITE(numout,*) ' ==========  bdy data files must contain 1 or 12 time dumps.'          
    834               IF(lwp) WRITE(numout,*) 
    835               nstop = nstop + 1 
    836             ELSEIF (kbdy_bt.EQ.1 ) THEN 
    837               IF(lwp) WRITE(numout,*) 
    838               IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 
    839               IF(lwp) WRITE(numout,*)              
    840             ELSEIF (kbdy_bt.EQ.12) THEN 
    841               IF(lwp) WRITE(numout,*) 
    842               IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 
    843               IF(lwp) WRITE(numout,*)  
    844             ENDIF 
    845           ENDIF 
    846  
    847       ! Find index of first record to read (before first model time).  
    848  
    849           jt=1 
    850           DO WHILE ( ((istep_bt(jt+1)) <= 0 ).AND.(jt.LE.(kbdy_bt-1))) 
    851             jt=jt+1 
    852           END DO 
    853           nbdy1_bt = jt 
    854  
    855           WRITE(numout,*) 'Time offset is ',offset 
    856           WRITE(numout,*) 'First record to read is ',nbdy1_bt 
    857  
    858         ENDIF ! endif (nbdy_dta == 1) 
    859  
    860       ! 1.2  Read first record in file if necessary (ie if nbdy_dta == 1) 
    861       ! ***************************************************************** 
    862  
    863         IF ( nbdy_dta == 0) THEN 
    864           ! boundary data arrays are filled with initial conditions 
    865           jgrd = 2            ! U-points data  
    866           DO jb = 1, nblen(jgrd)               
    867             ubtbdy(jb) = un(nbi(jb,jgrd), nbj(jb,jgrd), 1) 
    868           END DO 
    869  
    870           jgrd = 3            ! V-points data  
    871           DO jb = 1, nblen(jgrd)               
    872             vbtbdy(jb) = vn(nbi(jb,jgrd), nbj(jb,jgrd), 1) 
    873           END DO 
    874  
    875           jgrd = 1            ! T-points data  
    876           DO jb = 1, nblen(jgrd)               
    877             sshbdy(jb) = sshn(nbi(jb,jgrd), nbj(jb,jgrd)) 
    878           END DO 
    879  
    880         ELSEIF (nbdy_dta == 1) THEN 
    881   
    882         ! Set first record in the climatological case:    
    883           IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN 
    884             nbdy2_bt = 1 
    885           ELSEIF ((ln_bdy_clim).AND.(kbdy_bt==iman)) THEN 
    886             nbdy1_bt = 0 
    887             nbdy2_bt = imois 
    888           ELSE 
    889             nbdy2_bt = nbdy1_bt 
    890           END IF 
    891   
    892          ! Open Netcdf files: 
    893  
    894           CALL iom_open ( filbdy_data_bt_T, bdynumt_bt ) 
    895           CALL iom_open ( filbdy_data_bt_U, bdynumu_bt ) 
    896           CALL iom_open ( filbdy_data_bt_V, bdynumv_bt ) 
    897  
    898          ! Read first record: 
    899           ipj=1 
    900           jgrd=1 
    901           ipi=nblendta(jgrd) 
    902  
    903           ! ssh 
    904           jgrd=1 
    905           IF ( nblendta(jgrd) .le. 0 ) THEN  
    906             idvar = iom_varid( bdynumt_bt,'sossheig' ) 
    907             nblendta(jgrd) = iom_file(bdynumt_bt)%dimsz(1,idvar) 
    908           ENDIF 
    909           WRITE(numout,*) 'Dim size for sossheig is ',nblendta(jgrd) 
    910           ipi=nblendta(jgrd) 
    911  
    912           CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt ) 
    913  
    914           DO jb=1, nblen(jgrd) 
    915             sshbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1) 
    916           END DO 
    917   
    918           ! u-velocity 
    919           jgrd=2 
    920           IF ( nblendta(jgrd) .le. 0 ) THEN  
    921             idvar = iom_varid( bdynumu_bt,'vobtcrtx' ) 
    922             nblendta(jgrd) = iom_file(bdynumu_bt)%dimsz(1,idvar) 
    923           ENDIF 
    924           WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(jgrd) 
    925           ipi=nblendta(jgrd) 
    926  
    927           CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt ) 
    928  
    929           DO jb=1, nblen(jgrd) 
    930             ubtbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1) 
    931           END DO 
    932  
    933           ! v-velocity 
    934           jgrd=3 
    935           IF ( nblendta(jgrd) .le. 0 ) THEN  
    936             idvar = iom_varid( bdynumv_bt,'vobtcrty' ) 
    937             nblendta(jgrd) = iom_file(bdynumv_bt)%dimsz(1,idvar) 
    938           ENDIF 
    939           WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(jgrd) 
    940           ipi=nblendta(jgrd) 
    941  
    942           CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt ) 
    943  
    944           DO jb=1, nblen(jgrd) 
    945             vbtbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1) 
    946           END DO 
    947  
    948         END IF 
    949   
    950         ! In the case of constant boundary forcing fill bdy arrays once for all 
    951         IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN 
    952  
    953           ubtbdy  (:) = ubtbdydta  (:,2) 
    954           vbtbdy  (:) = vbtbdydta  (:,2) 
    955           sshbdy  (:) = sshbdydta  (:,2) 
    956  
    957           CALL iom_close( bdynumt_bt ) 
    958           CALL iom_close( bdynumu_bt ) 
    959           CALL iom_close( bdynumv_bt ) 
    960  
    961         END IF 
    962  
    963       ENDIF ! End if nit000 
    964  
    965       ! -------------------- ! 
    966       ! 2. At each time step ! 
    967       ! -------------------- ! 
    968  
    969       IF ((nbdy_dta==1).AND.(kbdy_bt>1)) THEN  
    970  
    971       ! 2.1 Read one more record if necessary 
    972       !************************************** 
    973  
    974         IF ( (ln_bdy_clim).AND.(imois/=nbdy1_bt) ) THEN ! remember that nbdy1_bt=0 for kt=nit000 
    975          nbdy1_bt = imois 
    976          nbdy2_bt = imois+1 
    977          nbdy1_bt = MOD( nbdy1_bt, iman ) 
    978          IF( nbdy1_bt == 0 ) nbdy1_bt = iman 
    979          nbdy2_bt = MOD( nbdy2_bt, iman ) 
    980          IF( nbdy2_bt == 0 ) nbdy2_bt = iman 
    981          lect=.true. 
    982  
    983         ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy2_bt))) THEN 
    984           nbdy1_bt=nbdy2_bt 
    985           nbdy2_bt=nbdy2_bt+1 
     996        ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN 
     997          nbdy_b_bt=nbdy_a_bt 
     998          nbdy_a_bt=nbdy_a_bt+1 
    986999          lect=.true. 
    9871000        END IF 
     
    9981011          ipj=1 
    9991012          ipk=jpk 
    1000           jgrd=1 
    1001           ipi=nblendta(jgrd) 
     1013          igrd=1 
     1014          ipi=nblendta(igrd) 
    10021015 
    10031016           
    10041017          ! ssh 
    1005           jgrd=1 
    1006           ipi=nblendta(jgrd) 
    1007  
    1008           CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt ) 
    1009  
    1010           DO jb=1, nblen(jgrd) 
    1011             sshbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1) 
     1018          igrd=1 
     1019          ipi=nblendta(igrd) 
     1020 
     1021          CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
     1022 
     1023          DO ib=1, nblen(igrd) 
     1024            sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    10121025          END DO 
    10131026 
    10141027          ! u-velocity 
    1015           jgrd=2 
    1016           ipi=nblendta(jgrd) 
    1017  
    1018           CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt ) 
    1019  
    1020           DO jb=1, nblen(jgrd) 
    1021             ubtbdydta(jb,2) =  pdta3(nbmap(jb,jgrd),1) 
     1028          igrd=2 
     1029          ipi=nblendta(igrd) 
     1030 
     1031          CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
     1032 
     1033          DO ib=1, nblen(igrd) 
     1034            ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
    10221035          END DO 
    10231036 
    10241037          ! v-velocity 
    1025           jgrd=3 
    1026           ipi=nblendta(jgrd) 
    1027  
    1028           CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt ) 
    1029  
    1030           DO jb=1, nblen(jgrd) 
    1031             vbtbdydta(jb,2) =  pdta3(nbmap(jb,jgrd),1) 
    1032           END DO 
    1033  
    1034  
    1035          IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy1_bt ',nbdy1_bt 
    1036          IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy2_bt ',nbdy2_bt 
     1038          igrd=3 
     1039          ipi=nblendta(igrd) 
     1040 
     1041          CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 
     1042 
     1043          DO ib=1, nblen(igrd) 
     1044            vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1) 
     1045          END DO 
     1046 
     1047 
     1048         IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b_bt ',nbdy_b_bt 
     1049         IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a_bt ',nbdy_a_bt 
    10371050         IF (.NOT.ln_bdy_clim) THEN 
    1038            IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy1_bt) 
     1051           IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy_b_bt) 
    10391052           IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer 
    1040            IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy2_bt) 
     1053           IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt) 
    10411054         ENDIF 
    10421055        END IF ! end lect=.true. 
     
    10471060     
    10481061        IF (ln_bdy_clim) THEN 
    1049           zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1_bt) ) + 0.5 - i15 
     1062          zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b_bt) ) + 0.5 - i15 
    10501063        ELSE           
    1051           zxy = FLOAT(istep_bt(nbdy1_bt)-itimer) / FLOAT(istep_bt(nbdy1_bt)-istep_bt(nbdy2_bt)) 
     1064          zxy = FLOAT(istep_bt(nbdy_b_bt)-itimer) / FLOAT(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt)) 
    10521065        END IF 
    10531066 
    1054           jgrd=1 
    1055           DO jb=1, nblen(jgrd) 
    1056             sshbdy(jb) = zxy      * sshbdydta(jb,2) + & 
    1057                        (1.-zxy) * sshbdydta(jb,1)    
    1058           END DO 
    1059  
    1060           jgrd=2 
    1061           DO jb=1, nblen(jgrd) 
    1062             ubtbdy(jb) = zxy      * ubtbdydta(jb,2) + & 
    1063                          (1.-zxy) * ubtbdydta(jb,1)    
    1064           END DO 
    1065  
    1066           jgrd=3 
    1067           DO jb=1, nblen(jgrd) 
    1068             vbtbdy(jb) = zxy      * vbtbdydta(jb,2) + & 
    1069                          (1.-zxy) * vbtbdydta(jb,1)    
    1070           END DO 
    1071  
    1072  
    1073       END IF !end if ((nbdy_dta==1).AND.(kbdy_bt>1)) 
     1067          igrd=1 
     1068          DO ib=1, nblen(igrd) 
     1069            sshbdy(ib) = zxy      * sshbdydta(ib,2) + & 
     1070                       (1.-zxy) * sshbdydta(ib,1)    
     1071          END DO 
     1072 
     1073          igrd=2 
     1074          DO ib=1, nblen(igrd) 
     1075            ubtbdy(ib) = zxy      * ubtbdydta(ib,2) + & 
     1076                         (1.-zxy) * ubtbdydta(ib,1)    
     1077          END DO 
     1078 
     1079          igrd=3 
     1080          DO ib=1, nblen(igrd) 
     1081            vbtbdy(ib) = zxy      * vbtbdydta(ib,2) + & 
     1082                         (1.-zxy) * vbtbdydta(ib,1)    
     1083          END DO 
     1084 
     1085 
     1086      END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) 
    10741087     
    10751088      ! ------------------- ! 
     
    10791092      ! Closing of the 3 files 
    10801093      IF( kt == nitend ) THEN 
    1081           CALL iom_close( bdynumt_bt ) 
    1082           CALL iom_close( bdynumu_bt ) 
    1083           CALL iom_close( bdynumv_bt ) 
     1094          CALL iom_close( numbdyt_bt ) 
     1095          CALL iom_close( numbdyu_bt ) 
     1096          CALL iom_close( numbdyv_bt ) 
    10841097      ENDIF 
    10851098 
    1086       ENDIF ! lk_bdy 
     1099      ENDIF ! ln_bdy_dyn_frs 
    10871100 
    10881101      END SUBROUTINE bdy_dta_bt 
     
    10901103 
    10911104#else 
    1092    !!------------------------------------------------------------------------------ 
    1093    !!   default option:           Dummy module NO Unstruct Open Boundary Conditions 
    1094    !!------------------------------------------------------------------------------ 
     1105   !!---------------------------------------------------------------------- 
     1106   !!   Dummy module                  NO Unstruct Open Boundary Conditions 
     1107   !!---------------------------------------------------------------------- 
    10951108CONTAINS 
    1096    SUBROUTINE bdy_dta( kt )             ! Dummy routine 
    1097       INTEGER, INTENT (in) :: kt 
     1109   SUBROUTINE bdy_dta( kt )              ! Empty routine 
    10981110      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt 
    10991111   END SUBROUTINE bdy_dta 
    1100  
    1101    SUBROUTINE bdy_dta_bt( kt, jit )             ! Dummy routine 
    1102       INTEGER, INTENT (in) :: kt, jit 
    1103       WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, jit 
     1112   SUBROUTINE bdy_dta_bt( kt, kit )      ! Empty routine 
     1113      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, kit 
    11041114   END SUBROUTINE bdy_dta_bt 
    11051115#endif 
  • trunk/NEMO/OPA_SRC/BDY/bdydyn.F90

    r911 r1125  
    11MODULE bdydyn 
    2    !!================================================================================= 
     2   !!====================================================================== 
    33   !!                       ***  MODULE  bdydyn  *** 
    4    !! Ocean dynamics:   Flow relaxation scheme of velocities on unstruc. open boundary 
    5    !!================================================================================= 
    6  
    7    !!--------------------------------------------------------------------------------- 
     4   !! Unstructured Open Boundary Cond. :   Flow relaxation scheme on velocities 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2005-02  (J. Chanut, A. Sellar)  Original code 
     7   !!             -   !  2007-07  (D. Storkey) Move Flather implementation to separate routine. 
     8   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     9   !!---------------------------------------------------------------------- 
     10#if defined key_bdy  
     11   !!---------------------------------------------------------------------- 
     12   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
     13   !!---------------------------------------------------------------------- 
    814   !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary 
    915   !!   bdy_dyn_fla    : Flather condition for barotropic solution 
    10    !!--------------------------------------------------------------------------------- 
    11 #if defined key_bdy || defined key_bdy_tides 
    12    !!--------------------------------------------------------------------------------- 
    13    !! * Modules used 
     16   !!---------------------------------------------------------------------- 
    1417   USE oce             ! ocean dynamics and tracers  
    1518   USE dom_oce         ! ocean space and time domain 
     
    1922   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2023   USE bdytides        ! for tidal harmonic forcing at boundary 
    21    USE in_out_manager 
     24   USE in_out_manager  ! 
    2225 
    2326   IMPLICIT NONE 
    2427   PRIVATE 
    2528 
    26    !! * Accessibility 
    27    PUBLIC bdy_dyn_frs  ! routine called in dynspg_flt (free surface case ONLY) 
    28 #if defined key_dynspg_exp || defined key_dynspg_ts 
    29    PUBLIC bdy_dyn_fla  ! routine called in dynspg_flt (free surface case ONLY) 
    30 #endif 
     29   PUBLIC   bdy_dyn_frs   ! routine called in dynspg_flt (free surface case ONLY) 
     30# if defined key_dynspg_exp || defined key_dynspg_ts 
     31   PUBLIC   bdy_dyn_fla   ! routine called in dynspg_flt (free surface case ONLY) 
     32# endif 
    3133 
    32    !!--------------------------------------------------------------------------------- 
     34   !!---------------------------------------------------------------------- 
     35   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     36   !! $Id: $  
     37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     38   !!---------------------------------------------------------------------- 
    3339 
    3440CONTAINS 
    3541 
    36    SUBROUTINE bdy_dyn_frs ( kt ) 
    37       !!------------------------------------------------------------------------------ 
    38       !!                      SUBROUTINE bdy_dyn_frs 
    39       !!                     *********************** 
     42   SUBROUTINE bdy_dyn_frs( kt ) 
     43      !!---------------------------------------------------------------------- 
     44      !!                  ***  SUBROUTINE bdy_dyn_frs  *** 
     45      !! 
    4046      !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the   
    4147      !!                case of unstructured open boundaries. 
     
    4450      !!               a three-dimensional baroclinic ocean model with realistic 
    4551      !!               topography. Tellus, 365-382. 
     52      !!---------------------------------------------------------------------- 
     53      INTEGER, INTENT( in ) ::   kt   ! Main time step counter 
    4654      !! 
    47       !! History : 
    48       !!   9.0  !  05-02 (J. Chanut, A. Sellar) Original 
    49       !!        !  07-07 (D. Storkey) Move Flather implementation to separate routine. 
    50       !!------------------------------------------------------------------------------ 
    51       !! * Arguments 
    52       INTEGER, INTENT( in ) ::   kt    ! Main time step counter 
     55      INTEGER  ::   ib, ik, igrd      ! dummy loop indices 
     56      INTEGER  ::   ii, ij            ! 2D addresses 
     57      REAL(wp) ::   zwgt              ! boundary weight 
     58      !!---------------------------------------------------------------------- 
     59      ! 
     60      IF(ln_bdy_dyn_frs) THEN ! If this is false, then this routine does nothing.  
    5361 
    54       !! * Local declarations 
    55       INTEGER ::   jb, jk, jgrd        ! dummy loop indices 
    56       INTEGER ::   ii, ij              ! 2D addresses 
    57       REAL(wp) :: zwgt                 ! boundary weight 
    58  
    59       !!------------------------------------------------------------------------------ 
    60       !!  OPA 9.0, LODYC-IPSL (2003) 
    61       !!------------------------------------------------------------------------------ 
    62  
    63       ! ---------------------------! 
    64       ! FRS on the total velocity :! 
    65       ! ---------------------------!   
    66    
    67       jgrd=2 !: Relaxation of zonal velocity 
    68       DO jb = 1, nblen(jgrd) 
    69         DO jk = 1, jpkm1 
    70           ii = nbi(jb,jgrd) 
    71           ij = nbj(jb,jgrd) 
    72           zwgt = nbw(jb,jgrd) 
    73  
    74           ua(ii,ij,jk) = ( ua(ii,ij,jk)*(1.-zwgt) + ubdy(jb,jk)*zwgt ) & 
    75                                                         * umask(ii,ij,jk) 
    76         END DO 
    77       END DO 
    78  
    79       jgrd=3 !: Relaxation of meridional velocity 
    80       DO jb = 1, nblen(jgrd) 
    81         DO jk = 1, jpkm1 
    82           ii = nbi(jb,jgrd) 
    83           ij = nbj(jb,jgrd) 
    84           zwgt = nbw(jb,jgrd) 
    85  
    86           va(ii,ij,jk) = ( va(ii,ij,jk)*(1.-zwgt) + vbdy(jb,jk)*zwgt ) & 
    87                                                         * vmask(ii,ij,jk) 
    88         END DO 
    89       END DO  
    90  
    91       CALL lbc_lnk( ua, 'U', 1. ) ! Boundary points should be updated 
    92       CALL lbc_lnk( va, 'V', 1. ) ! 
     62         IF( kt == nit000 ) THEN 
     63            IF(lwp) WRITE(numout,*) 
     64            IF(lwp) WRITE(numout,*) 'bdy_dyn : Flow Relaxation Scheme on momentum' 
     65            IF(lwp) WRITE(numout,*) '~~~~~~~' 
     66         ENDIF 
     67         ! 
     68         igrd = 2                      ! Relaxation of zonal velocity 
     69         DO ib = 1, nblen(igrd) 
     70            DO ik = 1, jpkm1 
     71               ii = nbi(ib,igrd) 
     72               ij = nbj(ib,igrd) 
     73               zwgt = nbw(ib,igrd) 
     74               ua(ii,ij,ik) = ( ua(ii,ij,ik) * ( 1.- zwgt ) + ubdy(ib,ik) * zwgt ) * umask(ii,ij,ik) 
     75            END DO 
     76         END DO 
     77         ! 
     78         igrd = 3                      ! Relaxation of meridional velocity 
     79         DO ib = 1, nblen(igrd) 
     80            DO ik = 1, jpkm1 
     81               ii = nbi(ib,igrd) 
     82               ij = nbj(ib,igrd) 
     83               zwgt = nbw(ib,igrd) 
     84               va(ii,ij,ik) = ( va(ii,ij,ik) * ( 1.- zwgt ) + vbdy(ib,ik) * zwgt ) * vmask(ii,ij,ik) 
     85            END DO 
     86         END DO  
     87         ! 
     88         CALL lbc_lnk( ua, 'U', 1. )   ! Boundary points should be updated 
     89         CALL lbc_lnk( va, 'V', 1. )   ! 
     90         ! 
     91      ENDIF ! ln_bdy_dyn_frs 
    9392 
    9493   END SUBROUTINE bdy_dyn_frs 
     94 
    9595 
    9696#if defined key_dynspg_exp || defined key_dynspg_ts 
    9797!! Option to use Flather with dynspg_flt not coded yet... 
    9898   SUBROUTINE bdy_dyn_fla 
    99       !!------------------------------------------------------------------------------ 
    100       !!                      SUBROUTINE bdy_dyn_fla 
    101       !!                     *********************** 
     99      !!---------------------------------------------------------------------- 
     100      !!                 ***  SUBROUTINE bdy_dyn_fla  *** 
     101      !!              
    102102      !!              - Apply Flather boundary conditions on normal barotropic velocities  
    103       !!                (ln_bdy_fla=.true.) 
     103      !!                (ln_bdy_dyn_fla=.true. or ln_bdy_tides=.true.) 
    104104      !! 
    105105      !! ** WARNINGS about FLATHER implementation: 
     
    113113      !!   ssh in the code is not updated). 
    114114      !! 
    115       !!             - Flather, R. A., 1976: A tidal model of the northwest European 
    116       !!               continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.      
    117       !! History : 
    118       !!   9.0  !  05-02 (J. Chanut, A. Sellar) Original 
    119       !!        !  07-07 (D. Storkey) Flather algorithm in separate routine. 
    120       !!------------------------------------------------------------------------------ 
    121       !! * Local declarations 
    122       INTEGER ::   jb, jk, jgrd, ji, jj, jim1, jip1, jjm1, jjp1   ! dummy loop indices 
    123       INTEGER ::   ii, ij                    ! 2D addresses 
    124       REAL(wp) ::  corr ! Flather correction 
    125       REAL(wp) :: zwgt                       ! boundary weight 
    126       REAL(wp) :: elapsed 
    127  
    128       !!------------------------------------------------------------------------------ 
    129       !!  OPA 9.0, LODYC-IPSL (2003) 
    130       !!------------------------------------------------------------------------------ 
     115      !! References:  Flather, R. A., 1976: A tidal model of the northwest European 
     116      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.      
     117      !!---------------------------------------------------------------------- 
     118      INTEGER  ::   ib, igrd                         ! dummy loop indices 
     119      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
     120      REAL(wp) ::   zcorr                            ! Flather correction 
     121      !!---------------------------------------------------------------------- 
    131122 
    132123      ! ---------------------------------! 
     
    134125      ! ---------------------------------!  
    135126      
     127      IF(ln_bdy_dyn_fla .OR. ln_bdy_tides) THEN ! If these are both false, then this routine does nothing.  
     128 
    136129      ! Fill temporary array with ssh data (here spgu): 
    137       jgrd = 1 
    138       DO jb = 1, nblenrim(jgrd) 
    139         ji = nbi(jb,jgrd) 
    140         jj = nbj(jb,jgrd) 
    141         spgu(ji, jj) = sshbdy(jb) 
    142 #if defined key_bdy_tides 
    143         spgu(ji, jj) = spgu(ji, jj) + sshtide(jb) 
    144 #endif 
     130      igrd = 1 
     131      spgu(:,:) = 0.0 
     132      DO ib = 1, nblenrim(igrd) 
     133         ii = nbi(ib,igrd) 
     134         ij = nbj(ib,igrd) 
     135         IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib) 
     136         IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(ib) 
    145137      END DO 
    146  
    147       jgrd = 2  !: Flather bc on u-velocity;  
    148                 ! remember that flagu=-1 if normal velocity direction is outward 
    149                 ! I think we should rather use after ssh ? 
    150  
    151       DO jb = 1, nblenrim(jgrd) 
    152         ji  = nbi(jb,jgrd) 
    153         jj  = nbj(jb,jgrd)  
    154         jim1 = ji+MAX(0, INT(flagu(jb))) ! T pts i-indice inside the boundary 
    155         jip1 = ji-MIN(0, INT(flagu(jb))) ! T pts i-indice outside the boundary  
    156       
    157         corr = - flagu(jb) * sqrt (grav / (hu_e(ji, jj) + 1.e-20) )  & 
    158                                       * ( sshn_e(jim1, jj) - spgu(jip1,jj) ) 
    159         ua_e(ji, jj) = ( ubtbdy(jb) + utide(jb) ) * hu_e(ji,jj) 
    160         if ( ln_bdy_fla ) then 
    161           ua_e(ji,jj) = ua_e(ji,jj) + corr * umask(ji,jj,1) * hu_e(ji,jj) 
    162         endif 
    163  
     138      ! 
     139      igrd = 2      ! Flather bc on u-velocity;  
     140      !             ! remember that flagu=-1 if normal velocity direction is outward 
     141      !             ! I think we should rather use after ssh ? 
     142      DO ib = 1, nblenrim(igrd) 
     143         ii  = nbi(ib,igrd) 
     144         ij  = nbj(ib,igrd)  
     145         iim1 = ii + MAX( 0, INT( flagu(ib) ) )   ! T pts i-indice inside the boundary 
     146         iip1 = ii - MIN( 0, INT( flagu(ib) ) )   ! T pts i-indice outside the boundary  
     147         ! 
     148         zcorr = - flagu(ib) * SQRT( grav / (hu_e(ii, ij) + 1.e-20) ) * ( sshn_e(iim1, ij) - spgu(iip1,ij) ) 
     149         ua_e(ii, ij) = ( ubtbdy(ib) + utide(ib) ) * hu_e(ii,ij) 
     150         ua_e(ii,ij) = ua_e(ii,ij) + zcorr * umask(ii,ij,1) * hu_e(ii,ij) 
    164151      END DO 
    165         
    166       jgrd = 3  !: Flather bc on v-velocity 
    167                 ! remember that flagv=-1 if normal velocity direction is outward 
    168                   
    169       DO jb = 1, nblenrim(jgrd) 
    170         ji  = nbi(jb,jgrd) 
    171         jj  = nbj(jb,jgrd)  
    172         jjm1 = jj+MAX(0, INT(flagv(jb))) ! T pts j-indice inside the boundary 
    173         jjp1 = jj-MIN(0, INT(flagv(jb))) ! T pts j-indice outside the boundary  
    174       
    175         corr = - flagv(jb) * sqrt (grav / (hv_e(ji, jj) + 1.e-20) )  & 
    176                                       * ( sshn_e(ji, jjm1) - spgu(ji,jjp1) ) 
    177         va_e(ji, jj) = ( vbtbdy(jb) + vtide(jb) ) * hv_e(ji,jj) 
    178         if ( ln_bdy_fla ) then 
    179           va_e(ji,jj) = va_e(ji,jj) + corr * vmask(ji,jj,1) * hv_e(ji,jj) 
    180         endif 
    181  
     152      ! 
     153      igrd = 3      ! Flather bc on v-velocity 
     154      !             ! remember that flagv=-1 if normal velocity direction is outward 
     155      DO ib = 1, nblenrim(igrd) 
     156         ii  = nbi(ib,igrd) 
     157         ij  = nbj(ib,igrd)  
     158         ijm1 = ij + MAX( 0, INT( flagv(ib) ) )   ! T pts j-indice inside the boundary 
     159         ijp1 = ij - MIN( 0, INT( flagv(ib) ) )   ! T pts j-indice outside the boundary  
     160         ! 
     161         zcorr = - flagv(ib) * SQRT( grav / (hv_e(ii, ij) + 1.e-20) ) * ( sshn_e(ii, ijm1) - spgu(ii,ijp1) ) 
     162         va_e(ii, ij) = ( vbtbdy(ib) + vtide(ib) ) * hv_e(ii,ij) 
     163         va_e(ii,ij) = va_e(ii,ij) + zcorr * vmask(ii,ij,1) * hv_e(ii,ij) 
    182164      END DO 
    183  
     165      ! 
    184166      CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated 
    185167      CALL lbc_lnk( va_e, 'V', 1. ) ! 
     168      ! 
     169      ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides 
    186170 
    187171   END SUBROUTINE bdy_dyn_fla 
     
    189173 
    190174#else 
    191    !!================================================================================= 
    192    !!                       ***  MODULE  bdydyn  *** 
    193    !! Ocean dynamics:   Flow relaxation scheme of velocities on unstruc. open boundary 
    194    !!================================================================================= 
     175   !!---------------------------------------------------------------------- 
     176   !!   Dummy module                   NO Unstruct Open Boundary Conditions 
     177   !!---------------------------------------------------------------------- 
    195178CONTAINS 
    196  
    197    SUBROUTINE bdy_dyn_frs 
    198                               ! No Unstructured open boundaries ==> empty routine 
     179   SUBROUTINE bdy_dyn_frs( kt )      ! Empty routine 
     180      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
    199181   END SUBROUTINE bdy_dyn_frs 
    200  
    201    SUBROUTINE bdy_dyn_fla 
    202                               ! No Unstructured open boundaries ==> empty routine 
     182   SUBROUTINE bdy_dyn_fla            ! Empty routine 
     183      WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?' 
    203184   END SUBROUTINE bdy_dyn_fla 
    204  
    205185#endif 
    206186 
     187   !!====================================================================== 
    207188END MODULE bdydyn 
  • trunk/NEMO/OPA_SRC/BDY/bdyini.F90

    r911 r1125  
    1  
    2  MODULE bdyini 
    3    !!================================================================================= 
     1MODULE bdyini 
     2   !!====================================================================== 
    43   !!                       ***  MODULE  bdyini  *** 
    5    !! Initialization of unstructured open boundaries 
    6    !!================================================================================= 
    7 #if defined key_bdy || defined key_bdy_tides 
    8    !!--------------------------------------------------------------------------------- 
    9    !!   'key_bdy'                                Unstructured Open Boundary Conditions 
    10    !!--------------------------------------------------------------------------------- 
     4   !! Unstructured open boundaries : initialisation 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
     7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module 
     8   !!             -   !  2007-01  (D. Storkey) Tidal forcing 
     9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     10   !!---------------------------------------------------------------------- 
     11#if defined key_bdy 
     12   !!---------------------------------------------------------------------- 
     13   !!   'key_bdy'                     Unstructured Open Boundary Conditions 
     14   !!---------------------------------------------------------------------- 
    1115   !!   bdy_init       : Initialization of unstructured open boundaries 
    12    !!--------------------------------------------------------------------------------- 
    13    !! * Modules used 
     16   !!---------------------------------------------------------------------- 
    1417   USE oce             ! ocean dynamics and tracers variables 
    1518   USE dom_oce         ! ocean space and time domain 
     
    2427   PRIVATE 
    2528 
    26    !! * Routine accessibility 
    27    PUBLIC bdy_init        ! routine called by opa.F90 
    28  
    29    !! * Substitutions 
    30  
    31    !!--------------------------------------------------------------------------------- 
    32    !!   OPA 9.0 , LODYC-IPSL  (2003) 
     29   PUBLIC   bdy_init   ! routine called by opa.F90 
     30 
     31   !!---------------------------------------------------------------------- 
     32   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     33   !! $Id: $  
     34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3335   !!--------------------------------------------------------------------------------- 
    3436 
     
    4749      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries 
    4850      !! 
    49       !! History : 
    50       !!  OPA 9.0 !  05-01  (J. Chanut, A. Sellar)  Original code 
    51       !!          !  07-01  (D. Storkey)  Update to use IOM module. 
    52       !!          !  07-01  (D. Storkey)  Tidal forcing. 
    5351      !!----------------------------------------------------------------------       
    54       !! * Local declarations 
    55       INTEGER ::   ji, jj, jk, jgrd,  & ! dummy loop indices 
    56                    jb, jr, icount,          & 
    57                    icountr, nb_rim, nb_len, nbr_max 
     52      INTEGER ::   ii, ij, ik, igrd, ib, ir   ! dummy loop indices 
     53      INTEGER ::   icount, icountr 
     54      INTEGER ::   ib_len, ibr_max 
    5855      INTEGER ::   iw, ie, is, in  
    5956      INTEGER ::   inum                 ! temporary logical unit 
    60       INTEGER ::   & ! temporary integers 
    61          dummy_id 
    62       INTEGER, DIMENSION (2)       ::   kdimsz 
    63       INTEGER, DIMENSION(jpbdta, jpbgrd) ::  &  !: Index arrays 
    64          nbidta, nbjdta, &              !: i and j indices of bdy dta 
    65          nbrdta                      !: Discrete distance from rim points 
    66       REAL(wp) :: & 
    67           efl, wfl, nfl, sfl            ! temporary scalars 
    68       REAL(wp) , DIMENSION(jpidta,jpjdta) ::   & 
    69          tmpmsk                         ! global domain mask 
    70       REAL(wp) , DIMENSION(jpbdta,1)      ::   & 
    71          ndta                           ! temporary array  
    72       CHARACTER(LEN=80),DIMENSION(3)     :: bdyfile 
    73  
    74       NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, & 
    75                       ln_bdy_clim, ln_bdy_vol, ln_bdy_fla, ln_bdy_mask,         & 
    76                       nbdy_dta, nb_rimwidth, volbdy 
    77  
     57      INTEGER ::   id_dummy             ! temporary integers 
     58      INTEGER ::   igrd_start, igrd_end ! start and end of loops on igrd 
     59      INTEGER, DIMENSION (2)             ::   kdimsz 
     60      INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
     61      INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbrdta           ! Discrete distance from rim points 
     62      REAL(wp) :: zefl, zwfl, znfl, zsfl                       ! temporary scalars 
     63      REAL(wp) , DIMENSION(jpidta,jpjdta) ::   zmask           ! global domain mask 
     64      REAL(wp) , DIMENSION(jpbdta,1)      ::   zdta            ! temporary array  
     65      CHARACTER(LEN=80),DIMENSION(3)      ::   clfile 
     66      !! 
     67      NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V,          & 
     68         &            ln_bdy_tides, ln_bdy_clim, ln_bdy_vol, ln_bdy_mask,                & 
     69         &            ln_bdy_dyn_fla, ln_bdy_dyn_frs, ln_bdy_tra_frs,                    & 
     70         &            nbdy_dta   , nb_rimwidth  , volbdy 
    7871      !!---------------------------------------------------------------------- 
    7972 
     
    8174      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries' 
    8275      IF(lwp) WRITE(numout,*) '~~~~~~~~' 
    83  
    84       IF( jperio /= 0 ) THEN 
    85          IF(lwp) WRITE(numout,*) 
    86          IF(lwp) WRITE(numout,*) ' E R R O R : Cyclic or symmetric,',   & 
    87             ' and unstructured open boundary condition are not compatible' 
    88          IF(lwp) WRITE(numout,*) ' ========== ' 
    89          IF(lwp) WRITE(numout,*) 
    90          nstop = nstop + 1 
    91       END IF 
     76      ! 
     77      IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,',   & 
     78           ' and unstructured open boundary condition are not compatible' ) 
    9279 
    9380#if defined key_obc 
    94       IF(lwp) WRITE(numout,*) 
    95       IF(lwp) WRITE(numout,*) ' E R R O R : Straight open boundaries,',   & 
    96             ' and unstructured open boundaries are not compatible' 
    97       IF(lwp) WRITE(numout,*) ' ========== ' 
    98       IF(lwp) WRITE(numout,*) 
    99       nstop = nstop + 1 
     81      CALL ctl_stop( 'Straight open boundaries,',   & 
     82           ' and unstructured open boundaries are not compatible' ) 
    10083#endif 
    10184 
    10285# if defined key_dynspg_rl 
    103       IF(lwp) WRITE(numout,*) 
    104       IF(lwp) WRITE(numout,*) ' E R R O R : Rigid lid,',   & 
    105             ' and unstructured open boundaries are not compatible' 
    106       IF(lwp) WRITE(numout,*) ' ========== ' 
    107       IF(lwp) WRITE(numout,*) 
    108       nstop = nstop + 1 
     86      CALL ctl_stop( 'Rigid lid,',   & 
     87           ' and unstructured open boundaries are not compatible' ) 
    10988#endif 
    11089 
    111       ! 0. Read namelist parameters 
     90      ! Read namelist parameters 
    11291      ! --------------------------- 
    113  
    11492      REWIND( numnam ) 
    11593      READ  ( numnam, nambdy ) 
     
    11896      IF(lwp) WRITE(numout,*) '         nambdy' 
    11997 
    120       IF ((nbdy_dta/=0).AND.(nbdy_dta/=1)) THEN ! Check nbdy_dta value 
    121         IF(lwp) WRITE(numout,*) 
    122         IF(lwp) WRITE(numout,*) ' E R R O R : nbdy_dta =',nbdy_dta,   & 
    123             'but it should have been 0 or 1' 
    124         IF(lwp) WRITE(numout,*) ' ========== ' 
    125         IF(lwp) WRITE(numout,*) 
    126         nstop = nstop + 1 
    127       ELSE 
    128         IF(lwp) WRITE(numout,*) ' ' 
    129         IF(lwp) WRITE(numout,*) '         data in file (=1) or     nbdy_dta = ', nbdy_dta 
    130         IF(lwp) WRITE(numout,*) '         initial state used (=0)' 
    131       END IF  
     98      ! Check nbdy_dta value 
     99      IF(lwp) WRITE(numout,*) 'nbdy_dta =', nbdy_dta       
     100      IF(lwp) WRITE(numout,*) ' ' 
     101      SELECT CASE( nbdy_dta ) 
     102      CASE( 0 ) 
     103        IF(lwp) WRITE(numout,*) '         initial state used for bdy data'         
     104      CASE( 1 ) 
     105        IF(lwp) WRITE(numout,*) '         boundary data taken from file' 
     106      CASE DEFAULT 
     107        CALL ctl_stop( 'nbdy_dta must be 0 or 1' ) 
     108      END SELECT 
    132109 
    133110      IF(lwp) WRITE(numout,*) ' ' 
    134111      IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth 
    135112 
     113      IF(lwp) WRITE(numout,*) ' ' 
     114      IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy 
     115 
    136116      IF (ln_bdy_vol) THEN 
    137         IF    (volbdy==1) THEN ! Check volbdy value 
    138           IF(lwp) WRITE(numout,*) ' ' 
    139           IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy 
     117        SELECT CASE ( volbdy ) ! Check volbdy value 
     118        CASE( 1 ) 
    140119          IF(lwp) WRITE(numout,*) '         The total volume will be constant' 
    141  
    142         ELSEIF (volbdy==0) THEN 
    143           IF(lwp) WRITE(numout,*) ' ' 
    144           IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy 
    145           IF(lwp) WRITE(numout,*) '         The total volume will vary according to & 
    146                                           &the surface E-P flux' 
    147         ELSE 
    148           IF(lwp) WRITE(numout,*) ' ' 
    149           IF(lwp) WRITE(numout,*) ' E R R O R : volbdy =',volbdy,   & 
    150             'but it should have been 0 or 1' 
    151           IF(lwp) WRITE(numout,*) ' ========== ' 
    152           IF(lwp) WRITE(numout,*) 
    153           nstop = nstop + 1 
    154         END IF 
     120        CASE( 0 ) 
     121          IF(lwp) WRITE(numout,*) '         The total volume will vary according' 
     122          IF(lwp) WRITE(numout,*) '         to the surface E-P flux' 
     123        CASE DEFAULT 
     124          CALL ctl_stop( 'volbdy must be 0 or 1' ) 
     125        END SELECT 
    155126      ELSE 
    156         IF(lwp) WRITE(numout,*) ' ' 
    157127        IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries' 
    158128        IF(lwp) WRITE(numout,*) ' ' 
    159129      ENDIF 
    160130 
    161       IF (ln_bdy_fla) THEN 
    162         IF(lwp) WRITE(numout,*) ' ' 
    163         IF(lwp) WRITE(numout,*) 'Flather bc with unstructured open boundaries' 
    164         IF(lwp) WRITE(numout,*) ' ' 
    165       ELSE 
    166         IF(lwp) WRITE(numout,*) ' ' 
    167         IF(lwp) WRITE(numout,*) 'NO Flather bc with unstructured open boundaries' 
    168         IF(lwp) WRITE(numout,*) ' ' 
    169       ENDIF 
    170  
    171       ! 0.5 Read tides namelist  
     131      IF (ln_bdy_tides) THEN 
     132        IF(lwp) WRITE(numout,*) ' ' 
     133        IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 
     134        IF(lwp) WRITE(numout,*) ' ' 
     135      ENDIF 
     136 
     137      IF (ln_bdy_dyn_fla) THEN 
     138        IF(lwp) WRITE(numout,*) ' ' 
     139        IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries' 
     140        IF(lwp) WRITE(numout,*) ' ' 
     141      ENDIF 
     142 
     143      IF (ln_bdy_dyn_frs) THEN 
     144        IF(lwp) WRITE(numout,*) ' ' 
     145        IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries' 
     146        IF(lwp) WRITE(numout,*) ' ' 
     147      ENDIF 
     148 
     149      IF (ln_bdy_tra_frs) THEN 
     150        IF(lwp) WRITE(numout,*) ' ' 
     151        IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries' 
     152        IF(lwp) WRITE(numout,*) ' ' 
     153      ENDIF 
     154 
     155      ! Read tides namelist  
    172156      ! ------------------------ 
    173  
    174       IF ( lk_bdy_tides )  CALL tide_init 
    175  
    176       ! 1. Read arrays defining unstructured open boundaries 
    177       ! ---------------------------------------------------- 
    178  
    179       ! 1.1 Read global 2D mask at T-points: bdytmask 
    180       ! ********************************************* 
    181       ! bdytmask=1 on the computational domain AND on open boundaries 
    182       !         =0 elsewhere    
     157      IF( ln_bdy_tides )   CALL tide_init 
     158 
     159      ! Read arrays defining unstructured open boundaries 
     160      ! ------------------------------------------------- 
     161 
     162      ! Read global 2D mask at T-points: bdytmask 
     163      ! ***************************************** 
     164      ! bdytmask = 1  on the computational domain AND on open boundaries 
     165      !          = 0  elsewhere    
    183166  
    184167      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 
    185         tmpmsk(: , : ) = 0.e0 
    186         tmpmsk(jpizoom+1:jpizoom+jpiglo-2,:  ) = 1.e0           
     168         zmask(         :                ,:) = 0.e0 
     169         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0           
    187170      ELSE IF ( ln_bdy_mask ) THEN 
    188         CALL iom_open( filbdy_mask, inum ) 
    189         CALL iom_get ( inum, jpdom_data, 'bdy_msk', tmpmsk(:,:) ) 
    190         CALL iom_close( inum ) 
     171         CALL iom_open( filbdy_mask, inum ) 
     172         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 
     173         CALL iom_close( inum ) 
    191174      ELSE 
    192         tmpmsk(:,:) = 1.0 
     175         zmask(:,:) = 1.e0 
    193176      ENDIF 
    194177 
    195178      ! Save mask over local domain       
    196       DO jj = 1, nlcj 
    197         DO ji = 1, nlci 
    198           bdytmask(ji,jj) = tmpmsk( mig(ji), mjg(jj)) 
    199         END DO 
     179      DO ij = 1, nlcj 
     180         DO ii = 1, nlci 
     181            bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 
     182         END DO 
    200183      END DO 
    201184 
    202185      ! Derive mask on U and V grid from mask on T grid 
    203       bdyumask(:,:)=0.e0 
    204       bdyvmask(:,:)=0.e0 
    205       DO jj=1, jpjm1 
    206         DO ji=1, jpim1 
    207           bdyumask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji+1, jj ) 
    208           bdyvmask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji  ,jj+1)   
    209         END DO 
     186      bdyumask(:,:) = 0.e0 
     187      bdyvmask(:,:) = 0.e0 
     188      DO ij=1, jpjm1 
     189         DO ii=1, jpim1 
     190            bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
     191            bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
     192         END DO 
    210193      END DO 
    211194 
     
    214197      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 
    215198 
    216       ! 1.2 Read discrete distance and mapping indices 
    217       ! ********************************************** 
    218         nbidta(:,:)=0. 
    219         nbjdta(:,:)=0. 
    220         nbrdta(:,:)=0. 
     199      ! Read discrete distance and mapping indices 
     200      ! ****************************************** 
     201      nbidta(:,:) = 0.e0 
     202      nbjdta(:,:) = 0.e0 
     203      nbrdta(:,:) = 0.e0 
    221204 
    222205      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 
    223  
    224         icount = 0 
    225       ! Define west boundary (from ji=2 to ji=1+nb_rimwidth): 
    226         DO jr=1,nb_rimwidth          
    227           DO jj=3,jpjglo-2 
    228             icount=icount+1 
    229             nbidta(icount,:) = jr + 1 + (jpizoom-1) 
    230             nbjdta(icount,:) = jj + (jpjzoom-1)  
    231             nbrdta(icount,:) = jr 
    232           END DO 
    233         END DO 
    234  
    235       ! Define east boundary (from ji=jpiglo-1 to ji=jpiglo-nb_rimwidth): 
    236         DO jr=1,nb_rimwidth          
    237           DO jj=3,jpjglo-2 
    238             icount=icount+1 
    239             nbidta(icount,:) = jpiglo-jr + (jpizoom-1) 
    240             nbidta(icount,2) = jpiglo-jr-1 + (jpizoom-1) ! special case for u points 
    241             nbjdta(icount,:) = jj + (jpjzoom-1) 
    242             nbrdta(icount,:) = jr 
    243           END DO 
    244         END DO 
     206         icount = 0 
     207         ! Define west boundary (from ii=2 to ii=1+nb_rimwidth): 
     208         DO ir = 1, nb_rimwidth          
     209            DO ij = 3, jpjglo-2 
     210               icount=icount+1 
     211               nbidta(icount,:) = ir + 1 + (jpizoom-1) 
     212               nbjdta(icount,:) = ij + (jpjzoom-1)  
     213               nbrdta(icount,:) = ir 
     214            END DO 
     215         END DO 
     216 
     217         ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth): 
     218         DO ir=1,nb_rimwidth          
     219            DO ij=3,jpjglo-2 
     220               icount=icount+1 
     221               nbidta(icount,:) = jpiglo-ir + (jpizoom-1) 
     222               nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points 
     223               nbjdta(icount,:) = ij + (jpjzoom-1) 
     224               nbrdta(icount,:) = ir 
     225            END DO 
     226         END DO 
    245227             
    246       ELSE 
    247  
    248       ! Read indices and distances in unstructured boundary data files  
    249  
    250         IF ( lk_bdy ) THEN  
    251           bdyfile(1) = filbdy_data_T 
    252           bdyfile(2) = filbdy_data_U 
    253           bdyfile(3) = filbdy_data_V 
    254         ELSE  
    255           ! In this case we have tides only at the boundaries 
    256           ! so read index arrays from tides files for first tidal component 
    257           bdyfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 
    258           bdyfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 
    259           bdyfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 
    260         ENDIF           
    261  
    262         DO jgrd = 1,3 
    263            CALL iom_open( bdyfile(jgrd), inum ) 
    264            dummy_id = iom_varid( inum, 'nbidta', kdimsz=kdimsz )   
    265            WRITE(numout,*) 'kdimsz : ',kdimsz 
    266            nb_len = kdimsz(1) 
    267            IF (nb_len > jpbdta) THEN 
    268              IF(lwp) WRITE(numout,*) 
    269              IF(lwp) WRITE(numout,*) ' E R R O R : jpbdta is too small:' 
    270              IF(lwp) WRITE(numout,*) ' ==========  Boundary array length in file is ', nb_len 
    271              IF(lwp) WRITE(numout,*) '             But jpbdta is ', jpbdta 
    272              IF(lwp) WRITE(numout,*) '             File :        ', bdyfile(jgrd) 
    273              IF(lwp) WRITE(numout,*) 
    274              nstop = nstop + 1 
    275            ENDIF 
    276            CALL iom_get ( inum, jpdom_unknown, 'nbidta', ndta(1:nb_len,:) ) 
    277            DO ji=1,nb_len 
    278              nbidta(ji,jgrd) = INT( ndta(ji,1) ) 
    279            ENDDO 
    280            CALL iom_get ( inum, jpdom_unknown, 'nbjdta', ndta(1:nb_len,:) ) 
    281            DO ji=1,nb_len 
    282              nbjdta(ji,jgrd) = INT( ndta(ji,1) ) 
    283            ENDDO 
    284            CALL iom_get ( inum, jpdom_unknown, 'nbrdta', ndta(1:nb_len,:) ) 
    285            DO ji=1,nb_len 
    286              nbrdta(ji,jgrd) = INT( ndta(ji,1) ) 
    287            ENDDO 
    288            CALL iom_close( inum ) 
    289  
    290            ! Check that rimwidth in file is big enough: 
    291            nbr_max = MAXVAL(nbrdta(:,jgrd)) 
    292            IF (nbr_max < nb_rimwidth) THEN 
    293            IF(lwp) WRITE(numout,*) 
    294            IF(lwp) WRITE(numout,*) ' E R R O R : Maximum rimwidth in file is ', nbr_max 
    295            IF(lwp) WRITE(numout,*) ' ==========  but nb_rimwidth is ', nb_rimwidth 
    296            IF(lwp) WRITE(numout,*) 
    297            nstop = nstop + 1 
    298            ELSE 
    299            IF(lwp) WRITE(numout,*) 
    300            IF(lwp) WRITE(numout,*) '             Maximum rimwidth in file is ', nbr_max 
    301            IF(lwp) WRITE(numout,*) 
    302            END IF            
    303  
    304         ENDDO 
    305  
    306       END IF  
    307  
    308       ! 1.3 Dispatch mapping indices and discrete distances on each processor 
    309       ! ********************************************************************* 
     228      ELSE            ! Read indices and distances in unstructured boundary data files  
     229 
     230         IF( ln_bdy_tides ) THEN  
     231            ! Read tides input files for preference in case there are 
     232            ! no bdydata files.  
     233            clfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 
     234            clfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 
     235            clfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 
     236         ELSE 
     237            clfile(1) = filbdy_data_T 
     238            clfile(2) = filbdy_data_U 
     239            clfile(3) = filbdy_data_V 
     240         ENDIF           
     241 
     242         ! how many files are we to read in? 
     243         igrd_start = 1 
     244         igrd_end   = 3 
     245         IF(.NOT. ln_bdy_tides ) THEN 
     246            IF(.NOT. (ln_bdy_dyn_fla) .AND..NOT. (ln_bdy_tra_frs)) THEN 
     247               ! No T-grid file. 
     248               igrd_start = 2 
     249            ELSEIF ( .NOT. ln_bdy_dyn_frs .AND..NOT. ln_bdy_dyn_fla ) THEN 
     250               ! No U-grid or V-grid file. 
     251               igrd_end   = 1          
     252            ENDIF 
     253         ENDIF 
     254 
     255         DO igrd = igrd_start, igrd_end 
     256            CALL iom_open( clfile(igrd), inum ) 
     257            id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz )   
     258            WRITE(numout,*) 'kdimsz : ',kdimsz 
     259            ib_len = kdimsz(1) 
     260            IF( ib_len > jpbdta) CALL ctl_stop(          & 
     261                'Boundary data array in file too long.', & 
     262                'File :', TRIM(clfile(igrd)),            & 
     263                'increase parameter jpbdta.' ) 
     264 
     265            CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) ) 
     266            DO ii = 1,ib_len 
     267               nbidta(ii,igrd) = INT( zdta(ii,1) ) 
     268            END DO 
     269            CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) ) 
     270            DO ii = 1,ib_len 
     271              nbjdta(ii,igrd) = INT( zdta(ii,1) ) 
     272            END DO 
     273            CALL iom_get ( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) ) 
     274            DO ii = 1,ib_len 
     275              nbrdta(ii,igrd) = INT( zdta(ii,1) ) 
     276            END DO 
     277            CALL iom_close( inum ) 
     278 
     279            ! Check that rimwidth in file is big enough: 
     280            ibr_max = MAXVAL( nbrdta(:,igrd) ) 
     281            IF(lwp) WRITE(numout,*) 
     282            IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 
     283            IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth 
     284            IF (ibr_max < nb_rimwidth) CALL ctl_stop( & 
     285                'nb_rimwidth is larger than maximum rimwidth in file' ) 
     286            ! 
     287         END DO 
     288         ! 
     289      ENDIF  
     290 
     291      ! Dispatch mapping indices and discrete distances on each processor 
     292      ! ***************************************************************** 
    310293      
    311       iw = mig(1)+1            ! if monotasking and no zoom, iw=2 
    312       ie = mig(1) + nlci-1-1   ! if monotasking and no zoom, ie=jpim1 
    313       is = mjg(1)+1            ! if monotasking and no zoom, is=2 
    314       in = mjg(1) + nlcj-1-1   ! if monotasking and no zoom, in=jpjm1 
    315  
    316       DO jgrd = 1, jpbgrd        
     294      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
     295      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
     296      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
     297      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
     298 
     299      DO igrd = igrd_start, igrd_end 
    317300        icount  = 0 
    318301        icountr = 0 
    319         DO nb_rim=1, nb_rimwidth 
    320           DO jb = 1, jpbdta 
    321           ! check if point is in local domain and equals nb_rim 
    322             IF ( (nbidta(jb,jgrd) >= iw ).AND. & 
    323                  (nbidta(jb,jgrd) <= ie ).AND. & 
    324                  (nbjdta(jb,jgrd) >= is ).AND. & 
    325                  (nbjdta(jb,jgrd) <= in ).AND. & 
    326                  (nbrdta(jb,jgrd) == nb_rim ) ) THEN 
    327  
    328               icount = icount  + 1 
    329              
    330               IF (nb_rim==1) icountr = icountr+1 
    331  
    332               IF (icount > jpbdim) THEN 
    333                 IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small' 
    334                 nstop = nstop + 1 
    335               ELSE 
    336                 nbi(icount, jgrd)  = nbidta(jb,jgrd)- mig(1)+1 
    337                 nbj(icount, jgrd)  = nbjdta(jb,jgrd)- mjg(1)+1 
    338                 nbr(icount, jgrd)  = nbrdta(jb,jgrd) 
    339                 nbmap(icount,jgrd) = jb 
    340               END IF             
    341             END IF 
    342           END DO 
    343         END DO 
    344         nblenrim(jgrd) = icountr !: length of rim boundary data on each proc 
    345         nblen   (jgrd) = icount  !: length of boundary data on each proc         
     302        nblen(igrd) = 0 
     303        nblenrim(igrd) = 0 
     304        nblendta(igrd) = 0 
     305        DO ir=1, nb_rimwidth 
     306          DO ib = 1, jpbdta 
     307          ! check if point is in local domain and equals ir 
     308            IF(  nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND.  & 
     309               & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND.   & 
     310               & nbrdta(ib,igrd) == ir  ) THEN 
     311               ! 
     312               icount = icount  + 1 
     313               ! 
     314               IF( ir == 1 )   icountr = icountr+1 
     315                  IF (icount > jpbdim) THEN 
     316                     IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small' 
     317                     nstop = nstop + 1 
     318                  ELSE 
     319                     nbi(icount, igrd)  = nbidta(ib,igrd)- mig(1)+1 
     320                     nbj(icount, igrd)  = nbjdta(ib,igrd)- mjg(1)+1 
     321                     nbr(icount, igrd)  = nbrdta(ib,igrd) 
     322                     nbmap(icount,igrd) = ib 
     323                  ENDIF             
     324               ENDIF 
     325            END DO 
     326         END DO 
     327         nblenrim(igrd) = icountr !: length of rim boundary data on each proc 
     328         nblen   (igrd) = icount  !: length of boundary data on each proc         
    346329      END DO  
    347330 
    348       ! 2. Compute rim weights 
    349       ! ---------------------- 
    350        
    351       DO  jgrd = 1, jpbgrd 
    352         DO jb = 1, nblen(jgrd) 
    353           ! tanh formulation 
    354           nbw(jb,jgrd) = 1.-TANH((FLOAT(nbr(jb,jgrd)-1))/2.) 
    355           ! quadratic 
    356 !          nbw(jb,jgrd) = (FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth))**2 
    357           ! linear 
    358 !          nbw(jb,jgrd) =  FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth) 
    359         END DO 
     331      ! Compute rim weights 
     332      ! ------------------- 
     333      DO igrd = igrd_start, igrd_end 
     334         DO ib = 1, nblen(igrd) 
     335            ! tanh formulation 
     336            nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 ) 
     337            ! quadratic 
     338!           nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2 
     339            ! linear 
     340!           nbw(ib,igrd) =  FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth) 
     341         END DO 
    360342      END DO  
    361343    
    362       ! 3. Mask corrections 
    363       ! ------------------- 
    364  
    365       DO jk=1, jpkm1 
    366         DO jj=1, jpj 
    367           DO ji=1, jpi 
    368             tmask(ji,jj,jk)=tmask(ji,jj,jk)*bdytmask(ji,jj) 
    369             umask(ji,jj,jk)=umask(ji,jj,jk)*bdyumask(ji,jj) 
    370             vmask(ji,jj,jk)=vmask(ji,jj,jk)*bdyvmask(ji,jj) 
    371             bmask(ji,jj)=bmask(ji,jj)*bdytmask(ji,jj) 
    372           END DO       
    373         END DO 
    374       END DO 
    375  
    376       ! I am not sure that it is useful:   
    377       DO jk=1, jpkm1 
    378         DO jj=2, jpjm1 
    379           DO ji=2, jpim1 
    380             fmask(ji,jj,jk) = fmask(ji,jj,jk) & 
    381                   &           * bdytmask(ji, jj ) * bdytmask(ji+1, jj )   & 
    382                   &           * bdytmask(ji,jj+1) * bdytmask(ji+1,jj+1) 
    383           END DO       
    384         END DO 
    385       END DO 
    386  
    387       tmask_i(:,:) = tmask(:,:,1)*tmask_i(:,:) 
    388  
    389       bdytmask(:,:)=tmask(:,:,1) 
     344      ! Mask corrections 
     345      ! ---------------- 
     346      DO ik = 1, jpkm1 
     347         DO ij = 1, jpj 
     348            DO ii = 1, jpi 
     349               tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
     350               umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
     351               vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
     352               bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij) 
     353            END DO       
     354         END DO 
     355      END DO 
     356 
     357      DO ik = 1, jpkm1 
     358         DO ij = 2, jpjm1 
     359            DO ii = 2, jpim1 
     360               fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
     361                  &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
     362            END DO       
     363         END DO 
     364      END DO 
     365 
     366      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)              
     367      bdytmask(:,:) = tmask(:,:,1) 
    390368 
    391369      ! bdy masks and bmask are now set to zero on boundary points: 
    392  
    393       jgrd=1 ! In the free surface case, bmask is at T-points 
    394       DO jb=1, nblenrim(jgrd)      
    395         bmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 
    396       END DO 
    397     
    398       jgrd=1 
    399       DO jb=1, nblenrim(jgrd)       
    400         bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 
    401       END DO 
    402  
    403       jgrd=2 
    404       DO jb=1, nblenrim(jgrd) 
    405         bdyumask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 
    406       END DO 
    407  
    408       jgrd=3 
    409       DO jb=1, nblenrim(jgrd) 
    410         bdyvmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 
    411       END DO 
    412  
    413        ! Lateral boundary conditions 
    414  
    415       CALL lbc_lnk( fmask, 'F', 1. ) 
     370      igrd = 1       ! In the free surface case, bmask is at T-points 
     371      DO ib = 1, nblenrim(igrd)      
     372        bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
     373      END DO 
     374      ! 
     375      igrd = 1 
     376      DO ib = 1, nblenrim(igrd)       
     377        bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
     378      END DO 
     379      ! 
     380      igrd = 2 
     381      DO ib = 1, nblenrim(igrd) 
     382        bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
     383      END DO 
     384      ! 
     385      igrd = 3 
     386      DO ib = 1, nblenrim(igrd) 
     387        bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 
     388      END DO 
     389 
     390      ! Lateral boundary conditions 
     391      CALL lbc_lnk( fmask        , 'F', 1. ) 
    416392      CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) 
    417393      CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) 
    418394      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 
    419395 
    420       IF ((ln_bdy_vol).OR.(ln_bdy_fla)) THEN  
    421  
    422       ! 4 Indices and directions of rim velocity components 
    423       ! --------------------------------------------------- 
    424         !flagu = -1 : u component is normal to the dynamical boundary  
    425         !             but its direction is outward 
    426         ! 
    427         !flagu =  0 : u is tangential 
    428         ! 
    429         !flagu =  1 : u is normal to the boundary  
    430         !             and is direction is inward 
     396      IF( ln_bdy_vol .OR. ln_bdy_dyn_fla ) THEN      ! Indices and directions of rim velocity components 
     397         ! 
     398         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 
     399         !flagu =  0 : u is tangential 
     400         !flagu =  1 : u is normal to the boundary and is direction is inward 
     401         icount = 0  
     402         flagu(:) = 0.e0 
    431403  
    432         icount = 0  
    433  
    434         flagu(:)=0.e0 
     404         igrd = 2      ! u-component  
     405         DO ib = 1, nblenrim(igrd)   
     406            zefl=bdytmask(nbi(ib,igrd)  , nbj(ib,igrd)) 
     407            zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd)) 
     408            IF( zefl + zwfl ==2 ) THEN 
     409               icount = icount +1 
     410            ELSE 
     411               flagu(ib)=-zefl+zwfl 
     412            ENDIF 
     413         END DO 
     414 
     415         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward 
     416         !flagv =  0 : u is tangential 
     417         !flagv =  1 : u is normal to the boundary and is direction is inward 
     418         flagv(:) = 0.e0 
     419 
     420         igrd = 3      ! v-component 
     421         DO ib = 1, nblenrim(igrd)   
     422            znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)) 
     423            zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1) 
     424            IF( znfl + zsfl ==2 ) THEN 
     425               icount = icount + 1 
     426            ELSE 
     427               flagv(ib) = -znfl + zsfl 
     428            END IF 
     429         END DO 
    435430  
    436         jgrd=2 ! u-component  
    437         DO jb=1, nblenrim(jgrd)   
    438           efl=bdytmask(nbi(jb,jgrd)  , nbj(jb,jgrd)) 
    439           wfl=bdytmask(nbi(jb,jgrd)+1, nbj(jb,jgrd)) 
    440           IF ((efl+wfl)==2) THEN 
    441             icount = icount +1 
    442           ELSE 
    443             flagu(jb)=-efl+wfl 
    444           END IF 
    445         END DO 
    446  
    447         !flagv = -1 : u component is normal to the dynamical boundary  
    448         !             but its direction is outward 
    449         ! 
    450         !flagv =  0 : u is tangential 
    451         ! 
    452         !flagv =  1 : u is normal to the boundary  
    453         !             and is direction is inward 
    454  
    455         flagv(:)=0.e0 
    456  
    457         jgrd=3 ! v-component 
    458         DO jb=1, nblenrim(jgrd)   
    459           nfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) 
    460           sfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)+1) 
    461           IF ((nfl+sfl)==2) THEN 
    462             icount = icount +1 
    463           ELSE 
    464             flagv(jb)=-nfl+sfl 
    465           END IF 
    466         END DO 
    467   
    468         IF( icount /= 0 ) THEN 
    469            IF(lwp) WRITE(numout,*) 
    470            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   & 
    471               ' are not boundary points. Check nbi, nbj, indices.' 
    472            IF(lwp) WRITE(numout,*) ' ========== ' 
    473            IF(lwp) WRITE(numout,*) 
    474            nstop = nstop + 1 
    475         END IF  
     431         IF( icount /= 0 ) THEN 
     432            IF(lwp) WRITE(numout,*) 
     433            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   & 
     434               ' are not boundary points. Check nbi, nbj, indices.' 
     435            IF(lwp) WRITE(numout,*) ' ========== ' 
     436            IF(lwp) WRITE(numout,*) 
     437            nstop = nstop + 1 
     438         ENDIF  
    476439     
    477       END IF 
    478  
    479       ! 5 Compute total lateral surface for volume correction: 
    480       ! ------------------------------------------------------ 
     440      ENDIF 
     441 
     442      ! Compute total lateral surface for volume correction: 
     443      ! ---------------------------------------------------- 
    481444  
    482445      bdysurftot = 0.e0  
    483      
    484       IF (ln_bdy_vol) THEN   
    485         jgrd=2 ! Lateral surface at U-points 
    486         DO jb=1, nblenrim(jgrd) 
    487           bdysurftot = bdysurftot + &  
    488                     hu(nbi(jb,jgrd), nbj(jb,jgrd)) * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    489                    * ABS(flagu(jb))*tmask_i(nbi(jb,jgrd)  , nbj(jb,jgrd)) & 
    490                                    *tmask_i(nbi(jb,jgrd)+1, nbj(jb,jgrd))                    
    491         END DO 
    492  
    493         jgrd=3 ! Add lateral surface at V-points 
    494         DO jb=1, nblenrim(jgrd) 
    495           bdysurftot = bdysurftot + &  
    496                     hv(nbi(jb,jgrd), nbj(jb,jgrd)) * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    497                    * ABS(flagv(jb))*tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    498                                    *tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)+1) 
    499         END DO 
    500  
    501         IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain 
     446      IF( ln_bdy_vol ) THEN   
     447         igrd = 2      ! Lateral surface at U-points 
     448         DO ib = 1, nblenrim(igrd) 
     449            bdysurftot = bdysurftot + hu     (nbi(ib,igrd)  ,nbj(ib,igrd))                      & 
     450               &                    * e2u    (nbi(ib,igrd)  ,nbj(ib,igrd)) * ABS( flagu(ib) )   & 
     451               &                    * tmask_i(nbi(ib,igrd)  ,nbj(ib,igrd))                      & 
     452               &                    * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd))                    
     453         END DO 
     454 
     455         igrd=3 ! Add lateral surface at V-points 
     456         DO ib = 1, nblenrim(igrd) 
     457            bdysurftot = bdysurftot + hv     (nbi(ib,igrd),nbj(ib,igrd)  )                      & 
     458               &                    * e1v    (nbi(ib,igrd),nbj(ib,igrd)  ) * ABS( flagv(ib) )   & 
     459               &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)  )                      & 
     460               &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1) 
     461         END DO 
     462 
     463         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain 
    502464      END IF    
    503465 
    504  
    505       ! 6. Initialise bdy data arrays 
    506       ! ----------------------------- 
    507  
     466      ! Initialise bdy data arrays 
     467      ! -------------------------- 
    508468      tbdy(:,:) = 0.e0 
    509469      sbdy(:,:) = 0.e0 
     
    514474      vbtbdy(:) = 0.e0 
    515475 
    516       ! 7. Read in tidal constituents and adjust for model start time 
    517       ! ------------------------------------------------------------- 
    518  
    519       IF ( lk_bdy_tides )  CALL tide_data 
    520  
     476      ! Read in tidal constituents and adjust for model start time 
     477      ! ---------------------------------------------------------- 
     478      IF( ln_bdy_tides )   CALL tide_data 
     479      ! 
    521480   END SUBROUTINE bdy_init 
    522481 
  • trunk/NEMO/OPA_SRC/BDY/bdytides.F90

    r911 r1125  
    11MODULE bdytides 
    2    !!================================================================================= 
     2   !!====================================================================== 
    33   !!                       ***  MODULE  bdytides  *** 
    44   !! Ocean dynamics:   Tidal forcing at open boundaries 
    5    !!================================================================================= 
    6 #if defined key_bdy_tides 
    7    !!--------------------------------------------------------------------------------- 
     5   !!====================================================================== 
     6   !! History :  2.0  !  2007-01  (D.Storkey)  Original code 
     7   !!            2.3  !  2008-01  (J.Holt)  Add date correction. Origins POLCOMS v6.3 2007 
     8   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     9   !!---------------------------------------------------------------------- 
     10#if defined key_bdy 
     11   !!---------------------------------------------------------------------- 
     12   !!   'key_bdy'     Unstructured Open Boundary Condition 
     13   !!---------------------------------------------------------------------- 
    814   !!   PUBLIC 
    9    !!   tide_init       : read of namelist 
    10    !!   tide_data       : read in and initialisation of tidal constituents at boundary 
    11    !!   tide_update     : calculation of tidal forcing at each timestep 
     15   !!      tide_init     : read of namelist 
     16   !!      tide_data     : read in and initialisation of tidal constituents at boundary 
     17   !!      tide_update   : calculation of tidal forcing at each timestep 
    1218   !!   PRIVATE 
    13    !!   uvset           :\ 
    14    !!   vday            : | Routines to correct tidal harmonics forcing for 
    15    !!   shpen           : | start time of integration 
    16    !!   ufset           : | 
    17    !!   vset            :/ 
    18    !!--------------------------------------------------------------------------------- 
    19  
    20    !!--------------------------------------------------------------------------------- 
    21    !! * Modules used 
     19   !!      uvset         :\ 
     20   !!      vday          : | Routines to correct tidal harmonics forcing for 
     21   !!      shpen         : | start time of integration 
     22   !!      ufset         : | 
     23   !!      vset          :/ 
     24   !!---------------------------------------------------------------------- 
    2225   USE oce             ! ocean dynamics and tracers  
    2326   USE dom_oce         ! ocean space and time domain 
     
    2932   USE bdy_oce         ! ocean open boundary conditions 
    3033   USE daymod          ! calendar 
     34 
    3135   IMPLICIT NONE 
    3236   PRIVATE 
    3337 
    34    !! * Accessibility 
    35    PUBLIC tide_init     ! routine called in bdyini 
    36    PUBLIC tide_data     ! routine called in bdyini 
    37    PUBLIC tide_update   ! routine called in bdydyn 
    38  
    39    !! * Module variables 
    40    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy_tides = .TRUE. !: tidal forcing at boundaries.  
    41  
    42    CHARACTER(len=80), PUBLIC :: & 
    43       filtide           !: Filename root for tidal input files 
    44  
    45    INTEGER, PUBLIC, PARAMETER ::  ntide_max = 15  ! Max number of tidal contituents 
    46    INTEGER, PUBLIC            ::  ntide           ! Actual number of tidal constituents 
    47  
    48    CHARACTER(len=4), PUBLIC, DIMENSION(ntide_max) :: & 
    49       tide_cpt         !: Names of tidal components used. 
    50  
    51    logical, PUBLIC      :: ln_tide_date ! if true correct tide phases and amplitude for model start date 
    52  
    53    REAL(wp), PUBLIC, DIMENSION(ntide_max) :: tide_speed ! Phase speed of tidal constituent (deg/hr) 
    54    INTEGER,          DIMENSION(ntide_max) :: indx 
    55    REAL(wp), DIMENSION(jpbdim,ntide_max) ::  & 
    56       ssh1, ssh2, & !: Tidal constituents : SSH 
    57       u1  , u2 ,       & !: Tidal constituents : U 
    58       v1  , v2           !: Tidal constituents : V 
     38   PUBLIC   tide_init     ! routine called in bdyini 
     39   PUBLIC   tide_data     ! routine called in bdyini 
     40   PUBLIC   tide_update   ! routine called in bdydyn 
     41 
     42   LOGICAL, PUBLIC            ::   ln_tide_date            !: =T correct tide phases and amplitude for model start date 
     43 
     44   INTEGER, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents 
     45   INTEGER            ::   ntide                 !: Actual number of tidal constituents 
     46 
     47   CHARACTER(len=80), PUBLIC                         ::   filtide    !: Filename root for tidal input files 
     48   CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used. 
     49 
     50   INTEGER , DIMENSION(jptides_max) ::   nindx        !: ??? 
     51   REAL(wp), DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
    5952    
     53   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   !: Tidal constituents : SSH 
     54   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     !: Tidal constituents : U 
     55   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     !: Tidal constituents : V 
    6056    
    61    !!--------------------------------------------------------------------------------- 
     57   !!---------------------------------------------------------------------- 
     58   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     59   !! $Id: $  
     60   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     61   !!---------------------------------------------------------------------- 
    6262 
    6363CONTAINS 
    6464 
    6565   SUBROUTINE tide_init 
    66       !!------------------------------------------------------------------------------ 
    67       !!                      SUBROUTINE tide_init  
    68       !!                     *********************** 
     66      !!---------------------------------------------------------------------- 
     67      !!                    ***  SUBROUTINE tide_init  *** 
     68      !!                      
    6969      !! ** Purpose : - Read in namelist for tides 
    7070      !! 
    71       !! History : 
    72       !!   NEMO v2.0  !  07-01 (D.Storkey) Original 
    73       !!------------------------------------------------------------------------------ 
    74       !! * Local declarations 
    75       INTEGER ::   jtide                  ! dummy loop index 
    76                                           ! different from nblendta!) 
    77   
    78       !!------------------------------------------------------------------------------ 
    79       !!  OPA 9.0, LODYC-IPSL (2007) 
    80       !!------------------------------------------------------------------------------ 
    81  
    82       NAMELIST/namtide/filtide, tide_cpt, tide_speed, ln_tide_date 
    83  
     71      !!---------------------------------------------------------------------- 
     72      INTEGER ::   itide                  ! dummy loop index  
     73      !! 
     74      NAMELIST/namtide/ln_tide_date, filtide, tide_cpt, tide_speed 
    8475      !!---------------------------------------------------------------------- 
    8576 
     
    8879      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
    8980 
    90       ! 0. Read namelist parameters 
    91       ! --------------------------- 
    92  
    93       do jtide = 1, ntide_max 
    94         tide_cpt(jtide) = '' 
    95       enddo 
    96  
    97       REWIND( numnam ) 
     81      ! Namelist namtide : tidal harmonic forcing at open boundaries 
     82      ln_tide_date = .false. 
     83      filtide(:) = '' 
     84      tide_speed(:) = 0.0 
     85      tide_cpt(:) = '' 
     86      REWIND( numnam )                                ! Read namelist parameters 
    9887      READ  ( numnam, namtide ) 
    99  
    100       ! control prints 
    101       IF(lwp) WRITE(numout,*) '         namtide' 
    102  
    103       ! Count number of components specified: 
    104       ntide=ntide_max 
    105       do jtide = 1, ntide_max 
    106         if ( tide_cpt(jtide) == '' ) then 
    107            ntide = jtide-1 
    108            exit 
    109         endif 
    110  
    111       enddo 
    112       ! find constients in standard list 
    113       do jtide = 1, ntide 
    114        indx(jtide) = 0 
    115        if ( TRIM(tide_cpt(jtide)) == 'Q1' ) indx(jtide) = 1 
    116        if ( TRIM(tide_cpt(jtide)) == 'O1' ) indx(jtide) = 2 
    117        if ( TRIM(tide_cpt(jtide)) == 'P1' ) indx(jtide) = 3 
    118        if ( TRIM(tide_cpt(jtide)) == 'S1' ) indx(jtide) = 4 
    119        if ( TRIM(tide_cpt(jtide)) == 'K1' ) indx(jtide) = 5 
    120        if ( TRIM(tide_cpt(jtide)) == '2N2' ) indx(jtide) = 6 
    121        if ( TRIM(tide_cpt(jtide)) == 'MU2' ) indx(jtide) = 7 
    122        if ( TRIM(tide_cpt(jtide)) == 'N2' ) indx(jtide) = 8 
    123        if ( TRIM(tide_cpt(jtide)) == 'NU2' ) indx(jtide) = 9 
    124        if ( TRIM(tide_cpt(jtide)) == 'M2' ) indx(jtide) = 10 
    125        if ( TRIM(tide_cpt(jtide)) == 'L2' ) indx(jtide) = 11 
    126        if ( TRIM(tide_cpt(jtide)) == 'T2' ) indx(jtide) = 12 
    127        if ( TRIM(tide_cpt(jtide)) == 'S2' ) indx(jtide) = 13 
    128        if ( TRIM(tide_cpt(jtide)) == 'K2' ) indx(jtide) = 14 
    129        if ( TRIM(tide_cpt(jtide)) == 'M4' ) indx(jtide) = 15 
    130       if (indx(jtide) == 0 ) then 
    131        if(lwp) write(numout,*) 'WARNING: constitunent', jtide,':', tide_cpt(jtide) & 
    132                                                          , 'not in standard list' 
    133       endif 
    134       enddo 
    135 ! 
    136       if ( ntide < 1 ) then  
    137          if(lwp) write(numout,*) 
    138          if(lwp) write(numout,*) ' E R R O R : Did not find any tidal components in namelist.' 
    139          if(lwp) write(numout,*) ' ========== ' 
    140          if(lwp) write(numout,*) 
    141          nstop = nstop + 1 
    142       else 
    143          if(lwp) write(numout,*) 
    144          if(lwp) write(numout,*) ntide,' tidal components specified : ' 
    145          if(lwp) write(numout,*) tide_cpt(1:ntide) 
    146          if(lwp) write(numout,*) ntide,' phase speeds (deg/hr) : ' 
    147          if(lwp) write(numout,*) tide_speed(1:ntide) 
    148       endif 
     88      !                                               ! Count number of components specified 
     89      ntide = jptides_max 
     90      itide = 1 
     91      DO WHILE( tide_cpt(itide) /= '' ) 
     92         ntide = itide 
     93         itide = itide + 1 
     94      END DO 
     95      !                                               ! find constituents in standard list 
     96      DO itide = 1, ntide 
     97         nindx(itide) = 0 
     98         IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
     99         IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
     100         IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
     101         IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
     102         IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
     103         IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
     104         IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
     105         IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
     106         IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
     107         IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
     108         IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
     109         IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
     110         IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
     111         IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
     112         IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
     113         IF( nindx(itide) == 0  .AND. lwp ) THEN 
     114            WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
     115            CALL ctl_warn( ctmp1 ) 
     116         ENDIF 
     117      END DO 
     118      !                                               ! Parameter control and print 
     119      IF( ntide < 1 ) THEN  
     120         CALL ctl_stop( '          Did not find any tidal components in namelist namtide' ) 
     121      ELSE 
     122         IF(lwp) WRITE(numout,*) '          Namelist namtide : tidal harmonic forcing at open boundaries' 
     123         IF(lwp) WRITE(numout,*) '             tidal components specified ', ntide 
     124         IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:ntide) 
     125         IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
     126         IF(lwp) WRITE(numout,*) '                ', tide_speed(1:ntide) 
     127      ENDIF 
    149128 
    150129      ! Initialisation of tidal harmonics arrays 
    151  
    152130      sshtide(:) = 0.e0 
    153       utide(:) = 0.e0 
    154       vtide(:) = 0.e0 
    155  
     131      utide  (:) = 0.e0 
     132      vtide  (:) = 0.e0 
     133      ! 
    156134   END SUBROUTINE tide_init 
    157135 
     136 
    158137   SUBROUTINE tide_data 
     138      !!---------------------------------------------------------------------- 
     139      !!                    ***  SUBROUTINE tide_data  *** 
     140      !!                     
     141      !! ** Purpose : - Read in tidal harmonics data and adjust for the start  
     142      !!                time of the model run.  
     143      !! 
     144      !!---------------------------------------------------------------------- 
     145      INTEGER ::   itide, igrd, ib        ! dummy loop indices 
     146      CHARACTER(len=80) :: clfile         ! full file name for tidal input file  
     147      INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read) 
     148      INTEGER, DIMENSION(3) :: lendta     ! length of data in the file (note may be different from nblendta!) 
     149      REAL(wp) ::  z_arg, z_atde, z_btde, z1t, z2t            
     150      REAL(wp), DIMENSION(jpbdta,1) ::   zdta   ! temporary array for data fields 
     151      REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc 
    159152      !!------------------------------------------------------------------------------ 
    160       !!                      SUBROUTINE tide_data 
    161       !!                     *********************** 
    162       !! ** Purpose : - Read in tidal harmonics data and adjust for the start time of 
    163       !!                the model run.  
    164       !! 
    165       !! History : 
    166       !!   NEMO v2.0  !  07-01 (D.Storkey) Original 
    167       !!              !  08-01 (J.Holt) Add date correction. 
    168       !!------------------------------------------------------------------------------ 
    169       !! * Local declarations 
    170       CHARACTER(len=80) :: tidefile       ! full file name for tidal input file  
    171       INTEGER ::   jtide, jgrd, jb        ! dummy loop indices 
    172       INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read) 
    173       INTEGER, DIMENSION(3) :: lendta     ! length of data in the file (note may be  
    174                                           ! different from nblendta!) 
    175       INTEGER :: iday,imonth,iyear 
    176       REAL(wp) ::  arg, atde, btde,z1t,z2t            
    177       REAL(wp), DIMENSION(jpbdta,1) ::   & 
    178          pdta                       ! temporary array for data fields 
    179       REAL(WP), DIMENSION(ntide_max) :: vplu, ftc 
    180   
    181       !!------------------------------------------------------------------------------ 
    182       !!  OPA 9.0, LODYC-IPSL (2007) 
    183       !!------------------------------------------------------------------------------ 
    184  
    185       ! 1. Open files and read in tidal forcing data 
    186       ! -------------------------------------------- 
     153 
     154      ! Open files and read in tidal forcing data 
     155      ! ----------------------------------------- 
    187156 
    188157      ipj = 1 
    189158 
    190       do jtide = 1 ,ntide 
    191  
    192          ! SSH fields 
    193          ! ---------- 
    194  
    195          tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_T.nc' 
    196          if(lwp) write(numout,*) 'Reading data from file ',tidefile 
    197          CALL iom_open( tidefile, inum ) 
    198  
    199          jgrd = 1 
    200          IF ( nblendta(jgrd) .le. 0 ) THEN  
    201            idvar = iom_varid( inum,'z1' ) 
    202            if(lwp) write(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 
    203            nblendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 
    204            WRITE(numout,*) 'Dim size for z1 is ',nblendta(jgrd) 
     159      DO itide = 1, ntide 
     160         !                                                              ! SSH fields 
     161         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
     162         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     163         CALL iom_open( clfile, inum ) 
     164         igrd = 1 
     165         IF( nblendta(igrd) <= 0 ) THEN  
     166            idvar = iom_varid( inum,'z1' ) 
     167            IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 
     168            nblendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
     169            WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd) 
    205170         ENDIF 
    206          ipi=nblendta(jgrd) 
    207  
    208          CALL iom_get ( inum, jpdom_unknown, 'z1', pdta(1:ipi,1:ipj) ) 
    209          DO jb=1, nblenrim(jgrd) 
    210             ssh1(jb,jtide) =  pdta(nbmap(jb,jgrd),1) 
    211          END DO 
    212  
    213          CALL iom_get ( inum, jpdom_unknown, 'z2', pdta(1:ipi,1:ipj) ) 
    214          DO jb=1, nblenrim(jgrd) 
    215             ssh2(jb,jtide) =  pdta(nbmap(jb,jgrd),1) 
    216          END DO 
    217  
     171         ipi = nblendta(igrd) 
     172         CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) ) 
     173         DO ib = 1, nblenrim(igrd) 
     174            ssh1(ib,itide) = zdta(nbmap(ib,igrd),1) 
     175         END DO 
     176         CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) ) 
     177         DO ib = 1, nblenrim(igrd) 
     178            ssh2(ib,itide) = zdta(nbmap(ib,igrd),1) 
     179         END DO 
    218180         CALL iom_close( inum ) 
    219  
    220          ! U fields 
    221          ! -------- 
    222  
    223          tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_U.nc' 
    224          if(lwp) write(numout,*) 'Reading data from file ',tidefile 
    225          CALL iom_open( tidefile, inum ) 
    226  
    227          jgrd = 2 
    228          IF ( lendta(jgrd) .le. 0 ) THEN  
    229            idvar = iom_varid( inum,'u1' ) 
    230            lendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 
    231            WRITE(numout,*) 'Dim size for u1 is ',lendta(jgrd) 
     181         ! 
     182         !                                                              ! U fields 
     183         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
     184         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     185         CALL iom_open( clfile, inum ) 
     186         igrd = 2 
     187         IF( lendta(igrd) <= 0 ) THEN  
     188            idvar = iom_varid( inum,'u1' ) 
     189            lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
     190            WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd) 
    232191         ENDIF 
    233          ipi=lendta(jgrd) 
    234  
    235          CALL iom_get ( inum, jpdom_unknown, 'u1', pdta(1:ipi,1:ipj) ) 
    236          DO jb=1, nblenrim(jgrd) 
    237             u1(jb,jtide) =  pdta(nbmap(jb,jgrd),1) 
    238          END DO 
    239  
    240          CALL iom_get ( inum, jpdom_unknown, 'u2', pdta(1:ipi,1:ipj) ) 
    241          DO jb=1, nblenrim(jgrd) 
    242             u2(jb,jtide) =  pdta(nbmap(jb,jgrd),1) 
    243          END DO 
    244  
     192         ipi = lendta(igrd) 
     193         CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) ) 
     194         DO ib = 1, nblenrim(igrd) 
     195            u1(ib,itide) = zdta(nbmap(ib,igrd),1) 
     196         END DO 
     197         CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) ) 
     198         DO ib = 1, nblenrim(igrd) 
     199            u2(ib,itide) = zdta(nbmap(ib,igrd),1) 
     200         END DO 
    245201         CALL iom_close( inum ) 
    246  
    247          ! V fields 
    248          ! -------- 
    249  
    250          tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_V.nc' 
    251          if(lwp) write(numout,*) 'Reading data from file ',tidefile 
    252          CALL iom_open( tidefile, inum ) 
    253  
    254          jgrd = 3 
    255          IF ( lendta(jgrd) .le. 0 ) THEN  
    256            idvar = iom_varid( inum,'v1' ) 
    257            lendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 
    258            WRITE(numout,*) 'Dim size for v1 is ',lendta(jgrd) 
     202         ! 
     203         !                                                              ! V fields 
     204         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
     205         if(lwp) write(numout,*) 'Reading data from file ', clfile 
     206         CALL iom_open( clfile, inum ) 
     207         igrd = 3 
     208         IF( lendta(igrd) <= 0 ) THEN  
     209            idvar = iom_varid( inum,'v1' ) 
     210            lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
     211            WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd) 
    259212         ENDIF 
    260          ipi=lendta(jgrd) 
    261  
    262          CALL iom_get ( inum, jpdom_unknown, 'v1', pdta(1:ipi,1:ipj) ) 
    263          DO jb=1, nblenrim(jgrd) 
    264             v1(jb,jtide) =  pdta(nbmap(jb,jgrd),1) 
    265          END DO 
    266  
    267          CALL iom_get ( inum, jpdom_unknown, 'v2', pdta(1:ipi,1:ipj) ) 
    268          DO jb=1, nblenrim(jgrd) 
    269             v2(jb,jtide) =  pdta(nbmap(jb,jgrd),1) 
    270          END DO 
    271  
     213         ipi = lendta(igrd) 
     214         CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) ) 
     215         DO ib = 1, nblenrim(igrd) 
     216            v1(ib,itide) = zdta(nbmap(ib,igrd),1) 
     217         END DO 
     218         CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) ) 
     219         DO ib=1, nblenrim(igrd) 
     220            v2(ib,itide) = zdta(nbmap(ib,igrd),1) 
     221         END DO 
    272222         CALL iom_close( inum ) 
    273          enddo ! end loop on tidal components 
    274 ! 
    275 !        correct for date factors 
    276 ! 
    277          if( ln_tide_date ) then 
    278  
    279  
    280 ! Calculate date corrects for 15 standard consituents 
    281           iyear  = int(ndate0 / 10000  )                           ! initial year 
    282           imonth = int((ndate0 - iyear * 10000 ) / 100 )          ! initial month 
    283           iday   = int(ndate0 - iyear * 10000 - imonth * 100)     ! initial day betwen 1 and 30 
    284  
    285             call uvset(0,iday,imonth,iyear,ftc,vplu) 
    286 ! 
    287             if(lwp) write(numout,*) 'Correcting tide for date:',iday,imonth,iyear 
    288  
    289             do jtide = 1, ntide 
    290 ! 
    291              if(indx(jtide) .ne. 0) then 
    292               arg=3.14159265d0*vplu(indx(jtide))/180.0d0 
    293               atde=ftc(indx(jtide))*cos(arg) 
    294               btde=ftc(indx(jtide))*sin(arg) 
    295               if(lwp) then 
    296                 write(numout,'(2i5,8f10.6)') jtide,indx(jtide),tide_speed(jtide), & 
    297                                    ftc(indx(jtide)),vplu(indx(jtide)) 
    298               endif             
    299              else 
    300              atde = 1.0_wp 
    301              btde = 0.0_wp 
    302              endif 
    303 !  elevation          
    304            jgrd = 1 
    305            do jb=1, nblenrim(jgrd)                 
    306             z1t=atde*ssh1(jb,jtide)+btde*ssh2(jb,jtide) 
    307             z2t=atde*ssh2(jb,jtide)-btde*ssh1(jb,jtide) 
    308             ssh1(jb,jtide) = z1t 
    309             ssh2(jb,jtide) = z2t 
    310            end do 
    311 !  u        
    312            jgrd = 2 
    313            do jb=1, nblenrim(jgrd)                 
    314             z1t=atde*u1(jb,jtide)+btde*u2(jb,jtide) 
    315             z2t=atde*u2(jb,jtide)-btde*u1(jb,jtide) 
    316             u1(jb,jtide) = z1t 
    317             u2(jb,jtide) = z2t 
    318            end do 
    319 !  v        
    320            jgrd = 3 
    321            do jb=1, nblenrim(jgrd)                 
    322             z1t=atde*v1(jb,jtide)+btde*v2(jb,jtide) 
    323             z2t=atde*v2(jb,jtide)-btde*v1(jb,jtide) 
    324             v1(jb,jtide) = z1t 
    325             v2(jb,jtide) = z2t 
    326            end do 
    327  
    328       enddo ! end loop on tidal components 
    329       endif ! date correction 
    330  
     223         ! 
     224      END DO ! end loop on tidal components 
     225 
     226      IF( ln_tide_date ) THEN      ! correct for date factors 
     227 
     228!! used nmonth, nyear and nday from daymod.... 
     229         ! Calculate date corrects for 15 standard consituents 
     230         ! This is the initialisation step, so nday, nmonth, nyear are the  
     231         ! initial date/time of the integration. 
     232 
     233         CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
     234 
     235         IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
     236 
     237         DO itide = 1, ntide       ! loop on tidal components 
     238            ! 
     239            IF( nindx(itide) /= 0 ) THEN 
     240!!gm use rpi  and rad global variable   
     241               z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
     242               z_atde=z_ftc(nindx(itide))*cos(z_arg) 
     243               z_btde=z_ftc(nindx(itide))*sin(z_arg) 
     244               IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide),   & 
     245                  &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
     246            ELSE 
     247               z_atde = 1.0_wp 
     248               z_btde = 0.0_wp 
     249            ENDIF 
     250            !                                         !  elevation          
     251            igrd = 1 
     252            DO ib = 1, nblenrim(igrd)                 
     253               z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide) 
     254               z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide) 
     255               ssh1(ib,itide) = z1t 
     256               ssh2(ib,itide) = z2t 
     257            END DO 
     258            !                                         !  u        
     259            igrd = 2 
     260            DO ib = 1, nblenrim(igrd)                 
     261               z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide) 
     262               z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide) 
     263               u1(ib,itide) = z1t 
     264               u2(ib,itide) = z2t 
     265            END DO 
     266            !                                         !  v        
     267            igrd = 3 
     268            DO ib = 1, nblenrim(igrd)                 
     269               z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide) 
     270               z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide) 
     271               v1(ib,itide) = z1t 
     272               v2(ib,itide) = z2t 
     273            END DO 
     274            ! 
     275         END DO   ! end loop on tidal components 
     276         ! 
     277      ENDIF ! date correction 
     278      ! 
    331279   END SUBROUTINE tide_data 
    332280 
     281 
    333282   SUBROUTINE tide_update ( kt, jit ) 
    334       !!------------------------------------------------------------------------------ 
    335       !!                      SUBROUTINE tide_update 
    336       !!                     ************************ 
     283      !!---------------------------------------------------------------------- 
     284      !!                 ***  SUBROUTINE tide_update  *** 
     285      !!                 
    337286      !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.  
    338287      !!                 
    339       !! 
    340       !! History : 
    341       !!   NEMO v2.0  !  06-12 (D.Storkey) Original 
    342       !!------------------------------------------------------------------------------ 
    343       !! * Arguments 
     288      !!---------------------------------------------------------------------- 
    344289      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter 
     290!!gm doctor jit ==> kit 
    345291      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option) 
    346  
    347       !! * Local declarations 
    348       INTEGER ::   jtide, jgrd, jb       ! dummy loop indices 
    349       REAL(wp) ::  arg, sarg                       
    350       REAL(wp), DIMENSION(ntide_max) :: sist, cost 
    351   
    352       !!------------------------------------------------------------------------------ 
    353       !!  OPA 9.0, LODYC-IPSL (2003) 
    354       !!------------------------------------------------------------------------------ 
     292      !! 
     293      INTEGER  ::   itide, igrd, ib      ! dummy loop indices 
     294      REAL(wp) ::   z_arg, z_sarg            !             
     295      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
     296      !!---------------------------------------------------------------------- 
    355297 
    356298      ! Note tide phase speeds are in deg/hour, so we need to convert the 
    357299      ! elapsed time in seconds to hours by dividing by 3600.0 
    358       if ( jit .eq. 0 ) then   
    359         arg = kt*rdt*rad/3600.0 
    360       else ! we are in a barotropic subcycle (for timesplitting option) 
    361         arg = ( (kt-1)*rdt + jit*rdtbt ) * rad/3600.0 
    362       endif 
    363  
    364       do jtide = 1,ntide 
    365         sarg = arg*tide_speed(jtide) 
    366         cost(jtide) = cos(sarg) 
    367         sist(jtide) = sin(sarg) 
    368       enddo 
    369  
    370       !! summing of tidal constituents into BDY arrays 
    371  
     300      IF( jit == 0 ) THEN   
     301         z_arg = kt * rdt * rad / 3600.0 
     302      ELSE                              ! we are in a barotropic subcycle (for timesplitting option) 
     303         z_arg = ( (kt-1) * rdt + jit * rdtbt ) * rad / 3600.0 
     304      ENDIF 
     305 
     306      DO itide = 1, ntide 
     307         z_sarg = z_arg * tide_speed(itide) 
     308         z_cost(itide) = COS( z_sarg ) 
     309         z_sist(itide) = SIN( z_sarg ) 
     310      END DO 
     311 
     312      ! summing of tidal constituents into BDY arrays 
    372313      sshtide(:) = 0.0 
    373       utide(:) = 0.0 
    374       vtide(:) = 0.0 
    375  
    376       do jtide = 1 ,ntide 
    377         jgrd=1 !: SSH on tracer grid. 
    378         do jb = 1, nblenrim(jgrd) 
    379           sshtide(jb) =sshtide(jb)+ ssh1(jb,jtide)*cost(jtide) + ssh2(jb,jtide)*sist(jtide) 
    380       !    if(lwp) write(numout,*) 'z',jb,jtide,sshtide(jb), ssh1(jb,jtide),ssh2(jb,jtide) 
    381         enddo 
    382         jgrd=2 !: U grid 
    383         do jb=1, nblenrim(jgrd) 
    384           utide(jb) = utide(jb)+ u1(jb,jtide)*cost(jtide) + u2(jb,jtide)*sist(jtide) 
    385       !    if(lwp) write(numout,*) 'u',jb,jtide,utide(jb), u1(jb,jtide),u2(jb,jtide) 
    386  
    387         enddo 
    388         jgrd=3 !: V grid 
    389         do jb=1, nblenrim(jgrd) 
    390           vtide(jb) = vtide(jb)+ v1(jb,jtide)*cost(jtide) + v2(jb,jtide)*sist(jtide) 
    391       !    if(lwp) write(numout,*) 'v',jb,jtide,vtide(jb), v1(jb,jtide),v2(jb,jtide) 
    392  
    393         enddo 
    394       enddo 
    395  
     314      utide (:) = 0.0 
     315      vtide (:) = 0.0 
     316      ! 
     317      DO itide = 1, ntide 
     318         igrd=1                              ! SSH on tracer grid. 
     319         DO ib = 1, nblenrim(igrd) 
     320            sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide) 
     321            !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide) 
     322         END DO 
     323         igrd=2                              ! U grid 
     324         DO ib=1, nblenrim(igrd) 
     325            utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide) 
     326            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide) 
     327         END DO 
     328         igrd=3                              ! V grid 
     329         DO ib=1, nblenrim(igrd) 
     330            vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide) 
     331            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide) 
     332         END DO 
     333      END DO 
     334      ! 
    396335   END SUBROUTINE tide_update 
    397336 
    398 ! 
    399 ! 
    400    SUBROUTINE uvset (ihs,iday,imnth,iyr,f,vplu) 
    401       !!------------------------------------------------------------------------------ 
    402       !!                      SUBROUTINE uvset 
    403       !!                     ************************ 
     337!!gm  doctor naming of dummy argument variables!!!   and all local variables 
     338 
     339   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 
     340      !!---------------------------------------------------------------------- 
     341      !!                   ***  SUBROUTINE uvset  *** 
     342      !!                      
    404343      !! ** Purpose : - adjust tidal forcing for date factors 
    405344      !!                 
    406       !! 
    407       !! History : 
    408       !! 
    409       !! Origins POLCOMS v6.3 2007 
    410       !!   NEMO v2.3  !  Jason Holt 
    411       !!------------------------------------------------------------------------------ 
     345      !!---------------------------------------------------------------------- 
    412346      implicit none 
    413       !! * Arguments 
    414347      INTEGER, INTENT( in ) ::   ihs     ! Start time hours 
    415348      INTEGER, INTENT( in ) ::   iday    ! start time days 
    416       INTEGER, INTENT( in ) ::   imnth  ! start time month 
    417       INTEGER, INTENT( in ) ::   iyr    ! start time year 
    418       INTEGER, PARAMETER ::  nc =15    ! maximum number of constituents 
    419       
    420       REAL(WP)              ::   f(nc)    ! nodal correction 
    421       REAL(WP)              ::   vplu(nc) ! phase correction 
    422 !     
    423       INTEGER         :: year,vd,ivdy,ndc,i,k 
    424       REAL(WP)        :: u(nc),v(nc),zig(nc),rtd 
    425       REAL(WP)        :: ss,h,p,en,p1 
     349      INTEGER, INTENT( in ) ::   imnth   ! start time month 
     350      INTEGER, INTENT( in ) ::   iyr     ! start time year 
     351      !! 
     352!!gm  nc is jptides_max ???? 
     353      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents 
    426354      CHARACTER(len=8), DIMENSION(nc) :: cname 
    427 ! 
    428       !!------------------------------------------------------------------------------ 
    429       !!  OPA 9.0, LODYC-IPSL (2007) 
    430       !!------------------------------------------------------------------------------ 
    431  
    432       data cname/ 'q1', 'o1', 'p1', 's1', 'k1', & 
    433                 '2n2','mu2', 'n2','nu2', 'm2', 'l2', 't2', 's2', 'k2', & 
    434                 'm4'/ 
    435       data zig/.2338507481,.2433518789,.2610826055,.2617993878 & 
    436      ,        .2625161701                                     & 
    437      ,        .4868657873,.4881373225,.4963669182,.4976384533 & 
    438      ,        .5058680490,.5153691799,.5228820265,.5235987756 & 
    439      ,        .5250323419                                     & 
    440      ,       1.011736098/ 
     355      INTEGER                 ::   year, vd, ivdy, ndc, i, k 
     356      REAL(wp)                ::   ss, h, p, en, p1, rtd 
     357      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction 
     358      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction 
     359      REAL(wp), DIMENSION(nc) ::   u, v, zig 
     360      !! 
     361      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   & 
     362         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   & 
     363         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      / 
     364      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   & 
     365         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   & 
     366         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   / 
     367      !!---------------------------------------------------------------------- 
    441368! 
    442369!     ihs             -  start time gmt on ... 
    443370!     iday/imnth/iyr  -  date   e.g.   12/10/87 
    444371! 
    445       call vday(iday,imnth,iyr,ivdy) 
    446       vd=ivdy 
     372      CALL vday(iday,imnth,iyr,ivdy) 
     373      vd = ivdy 
    447374! 
    448375!rp   note change of year number for d. blackman shpen 
     
    453380!.....year  =  year of required data 
    454381!.....vd    =  day of required data..set up for 0000gmt day year 
    455 ! 
    456382      ndc = nc 
    457383!.....ndc   =  number of constituents allowed 
    458 ! 
    459       rtd = 360.0/6.2831852 
    460       do i = 1,ndc 
    461       zig(i) = zig(i)*rtd 
    462 !      sigo(i)= zig(i) 
    463       enddo 
    464 ! 
    465       if(year == 0) then 
    466       return 
    467       endif 
     384!!gm use rpi ? 
     385      rtd = 360.0 / 6.2831852 
     386      DO i = 1, ndc 
     387         zig(i) = zig(i)*rtd 
     388         ! sigo(i)= zig(i) 
     389      END DO 
     390 
     391!!gm try to avoid RETURN  in F90 
     392      IF( year == 0 )   RETURN 
    468393       
    469       call shpen (year,vd,ss,h,p,en,p1) 
    470       call ufset (p,en,u,f) 
    471       call vset (ss,h,p,en,p1,v) 
    472 ! 
    473       do  k=1,nc 
    474        
    475       vplu(k)=v(k)+u(k) 
    476       vplu(k)=vplu(k)+dble(ihs)*zig(k) 
    477 ! 
    478       do while ( vplu(k) < 0 ) 
    479         vplu(k) = vplu(k) + 360.0 
    480       enddo 
    481 ! 
    482       do while (vplu(k) > 360. ) 
    483         vplu(k) = vplu(k) - 360.0 
    484       enddo 
    485 ! 
    486       enddo 
    487 ! 
     394         CALL shpen( year, vd, ss, h , p , en, p1 ) 
     395         CALL ufset( p   , en, u , f ) 
     396         CALL vset ( ss  , h , p , en, p1, v ) 
     397         ! 
     398         DO k = 1, nc 
     399            z_vplu(k) = v(k) + u(k) 
     400            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 
     401            DO WHILE( z_vplu(k) < 0    ) 
     402               z_vplu(k) = z_vplu(k) + 360.0 
     403            END DO 
     404            DO WHILE( z_vplu(k) > 360. ) 
     405               z_vplu(k) = z_vplu(k) - 360.0 
     406            END DO 
     407         END DO 
     408         ! 
    488409      END SUBROUTINE uvset 
    489410 
    490411 
    491       SUBROUTINE vday(iday,imnth,iy,ivdy) 
    492       !!------------------------------------------------------------------------------ 
    493       !!                      SUBROUTINE vday 
    494       !!                     **************** 
     412      SUBROUTINE vday( iday, imnth, iy, ivdy ) 
     413      !!---------------------------------------------------------------------- 
     414      !!                   *** SUBROUTINE vday  *** 
     415      !!                   
    495416      !! ** Purpose : - adjust tidal forcing for date factors 
    496417      !!                 
    497       !! 
    498       !! History : 
    499       !! 
    500       !! Origins POLCOMS v6.3 2007 
    501       !!   NEMO v2.3  !  Jason Holt 
     418      !!---------------------------------------------------------------------- 
     419      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ???? 
     420      INTEGER, INTENT(  out) ::   ivdy              ! ??? 
     421      !!  
     422      INTEGER ::   iyr 
    502423      !!------------------------------------------------------------------------------ 
    503       implicit none 
    504 ! 
    505 !     subroutine arguments 
    506       integer :: iday,imnth,iy,ivdy 
    507 ! 
    508 ! local variables 
    509       integer iyr 
    510       !!------------------------------------------------------------------------------ 
    511       !!  NEMO 2.3, LODYC-IPSL (2008) 
    512       !!------------------------------------------------------------------------------ 
    513  
    514 ! ================================================================= 
    515 ! 
    516 !     calculate day number in year from day/month/year 
    517 ! 
    518 ! ================================================================= 
     424 
     425!!gm   nday_year in day mode is the variable compiuted here, no? 
     426!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year 
     427 
     428      !calculate day number in year from day/month/year 
    519429       if(imnth.eq.1) ivdy=iday 
    520430       if(imnth.eq.2) ivdy=iday+31 
     
    529439       if(imnth.eq.11) ivdy=iday+304 
    530440       if(imnth.eq.12) ivdy=iday+334 
    531       iyr=iy 
     441       iyr=iy 
    532442       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
    533443       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 
    534444       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
    535  
    536       RETURN 
    537       END SUBROUTINE vday 
    538  
    539       SUBROUTINE shpen (year,vd,s,h,p,en,p1) 
    540       !!------------------------------------------------------------------------------ 
    541       !!                      SUBROUTINE shpen 
    542       !!                     ***************** 
     445      ! 
     446   END SUBROUTINE vday 
     447 
     448!!doctor norme for dummy arguments 
     449 
     450   SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 
     451      !!---------------------------------------------------------------------- 
     452      !!                    ***  SUBROUTINE shpen  *** 
     453      !!                    
    543454      !! ** Purpose : - calculate astronomical arguments for tides 
    544455      !!                this version from d. blackman 30 nove 1990 
    545       !!                 
    546       !! History : 
    547       !! 
    548       !! Origins POLCOMS v6.3 2007 
    549       !!   NEMO v2.3  !  Jason Holt 
    550       !!------------------------------------------------------------------------------ 
    551       implicit none 
    552  
    553 !     subroutine arguments 
    554       integer  ::year,vd 
    555       real(wp) :: s,h,p,en,p1 
    556  
    557 !     local variables 
    558  
    559       integer yr,ilc,icent,it,iday,ild,ipos,nn,iyd 
     456      !! 
     457      !!---------------------------------------------------------------------- 
     458!!gm add INTENT in, out or inout....    DOCTOR name.... 
     459!!gm please do not use variable name with a single letter:  impossible to search in a code 
     460      INTEGER  ::   year,vd 
     461      REAL(wp) ::   s,h,p,en,p1 
     462      !! 
     463      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd 
    560464      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat 
    561  
    562       !!------------------------------------------------------------------------------ 
    563       !!  NEMO 2.3, LODYC-IPSL (2008) 
    564       !!------------------------------------------------------------------------------ 
    565  
    566       data delt /-5.04,-3.90,-2.87,-0.58,0.71,1.80, & 
    567       3.08, 4.63, 5.86, 7.21, 8.58,10.50,12.10,    & 
    568      12.49,14.41,15.59,15.81,17.52,19.01,18.39,    & 
    569      19.55,20.36,21.01,21.81,21.76,22.35,22.68,    & 
    570      22.94,22.93,22.69,22.94,23.20,23.31,23.63,    & 
    571      23.47,23.68,23.62,23.53,23.59,23.99,23.80,    & 
    572      24.20,24.99,24.97,25.72,26.21,26.37,26.89,    & 
    573      27.68,28.13,28.94,29.42,29.66,30.29,30.96,    & 
    574      31.09,31.59,31.52,31.92,32.45,32.91,33.39,    & 
    575      33.80,34.23,34.73,35.40,36.14,36.99,37.87,    & 
    576      38.75,39.70,40.70,41.68,42.82,43.96,45.00,    & 
    577      45.98,47.00,48.03,49.10,50.10,50.97,51.81,    & 
    578      52.57/ 
    579 ! 
     465      !! 
     466      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           & 
     467         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    & 
     468         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    & 
     469         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    & 
     470         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    & 
     471         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    & 
     472         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    & 
     473         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    & 
     474         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    & 
     475         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    & 
     476         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    & 
     477         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    & 
     478         &       52.57 / 
     479      !!---------------------------------------------------------------------- 
     480 
    580481      cycle = 360.0 
    581482      ilc = 0 
    582       icent = year/100 
    583       yr = year - icent*100 
     483      icent = year / 100 
     484      yr = year - icent * 100 
    584485      t = icent - 20 
    585486! 
     
    588489!     see notes by cartwright 
    589490! 
     491!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90) 
     492!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    590493      it = icent - 20 
    591494      if (it) 1,2,2 
     
    602505      td = 0.0 
    603506! 
     507!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    604508      if (yr) 4,5,4 
    605509! 
     
    623527    7 deltat = delta * 1.0e-6 
    624528! 
     529!!gm   precision of the computation   :  example for s it should be replace by: 
     530!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results 
    625531      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 
    626       h   = 280.4661 + 36000.7698*t + 0.0003*t*t + 11.0*deltat 
    627       p   =  83.3535 + 4069.0139*t - 0.0103*t*t + deltat 
    628       en  = 234.9555 + 1934.1363*t - 0.0021*t*t + deltat 
    629       p1  = 282.9384 + 1.7195*t + 0.0005*t*t 
    630 ! 
    631       nn = s/cycle 
    632       s = s - nn*cycle 
    633       if ( s .lt. 0.0) s = s+cycle 
    634 ! 
    635       nn = h/cycle 
    636       h = h-cycle*nn 
    637       if (h .lt. 0.0) h = h+cycle 
    638 ! 
    639       nn = p/cycle 
    640       p = p- cycle*nn 
    641       if (p .lt. 0.0) p = p+cycle 
    642 ! 
    643       nn = en/cycle 
    644       en = en-cycle*nn 
    645       if(en .lt. 0.0) en = en + cycle 
     532      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t + 11.0*deltat 
     533      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +      deltat 
     534      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +      deltat 
     535      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t 
     536      ! 
     537      nn = s / cycle 
     538      s = s - nn * cycle 
     539      IF( s < 0.e0 )   s = s + cycle 
     540      ! 
     541      nn = h / cycle 
     542      h = h - cycle * nn 
     543      IF( h < 0.e0 )   h = h + cycle 
     544      ! 
     545      nn = p / cycle 
     546      p = p - cycle * nn 
     547      IF( p < 0.e0)   p = p + cycle 
     548      ! 
     549      nn = en / cycle 
     550      en = en - cycle * nn 
     551      IF( en < 0.e0 )  en = en + cycle 
    646552      en = cycle - en 
    647 ! 
    648       nn = p1/cycle 
    649       p1 = p1 - nn*cycle 
    650 ! 
    651       RETURN 
    652 ! 
    653       END SUBROUTINE shpen 
    654  
    655  
    656       SUBROUTINE ufset (p,cn,b,a) 
    657       !!------------------------------------------------------------------------------ 
    658       !!                      SUBROUTINE ufset 
    659       !!                     ***************** 
     553      ! 
     554      nn = p1 / cycle 
     555      p1 = p1 - nn * cycle 
     556      ! 
     557   END SUBROUTINE shpen 
     558 
     559 
     560   SUBROUTINE ufset( p, cn, b, a ) 
     561      !!---------------------------------------------------------------------- 
     562      !!                    ***  SUBROUTINE ufset  *** 
     563      !!                     
    660564      !! ** Purpose : - calculate nodal parameters for the tides 
    661565      !!                 
    662       !! 
    663       !! History : 
    664       !! 
    665       !! Origins POLCOMS v6.3 2007 
    666       !!   NEMO v2.3  !  Jason Holt 
    667       !!------------------------------------------------------------------------------ 
    668  
    669       implicit none 
     566      !!---------------------------------------------------------------------- 
     567!!gm  doctor naming of dummy argument variables!!!   and all local variables 
     568!!gm  nc is jptides_max ???? 
    670569      integer nc 
    671570      parameter (nc=15) 
    672 !     subroutine arguments 
    673       real(wp) p,cn 
    674 ! 
    675 ! 
    676 !     local variables 
    677       real(wp) ::   w1,w2,w3,w4,w5,w6,w7,w8,nw,pw,rad 
    678       real(wp) ::   a(nc),b(nc) 
    679       integer ndc,k 
    680 ! 
    681       !!------------------------------------------------------------------------------ 
    682       !!  NEMO 2.3, LODYC-IPSL (2008) 
    683       !!------------------------------------------------------------------------------ 
    684  
    685       ndc=nc 
    686 ! 
     571      REAL(wp) p,cn 
     572      !! 
     573       
     574!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z" 
     575      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 
     576      REAL(wp) ::   a(nc), b(nc) 
     577      INTEGER  ::   ndc, k 
     578      !!---------------------------------------------------------------------- 
     579 
     580      ndc = nc 
     581 
    687582!    a=f       ,  b =u 
    688583!    t is  zero as compared to tifa. 
     584!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module 
    689585      rad = 6.2831852d0/360.0 
    690       pw = p*rad 
    691       nw = cn*rad 
    692       w1 = cos(nw) 
    693       w2 = cos(2*nw) 
    694       w3 = cos(3*nw) 
    695       w4 = sin(nw) 
    696       w5 = sin(2*nw) 
    697       w6 = sin(3*nw) 
    698       w7 = 1 -0.2505*cos(2*pw) -0.1102*cos(2*pw-nw) & 
    699             -0.156*cos(2*pw-2*nw) -0.037*cos(nw) 
    700       w8 = -0.2505*sin(2*pw) -0.1102*sin(2*pw-nw) & 
    701           -0.0156*sin(2*pw-2*nw) -0.037*sin(nw) 
    702 ! 
    703       a(1) = 1.0089+0.1871*w1-0.0147*w2+0.0014*w3 
    704       b(1) = 0.1885*w4 - 0.0234*w5+.0033*w6 
    705 !   q1 
     586      pw = p  * rad 
     587      nw = cn * rad 
     588      w1 = cos(   nw ) 
     589      w2 = cos( 2*nw ) 
     590      w3 = cos( 3*nw ) 
     591      w4 = sin(   nw ) 
     592      w5 = sin( 2*nw ) 
     593      w6 = sin( 3*nw ) 
     594      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )  & 
     595         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw ) 
     596      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )  & 
     597         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw ) 
     598      ! 
     599      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 
     600      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 
     601      !   q1 
    706602      a(2) = a(1) 
    707603      b(2) = b(1) 
    708 !   o1 
     604      !   o1 
    709605      a(3) = 1.0 
    710606      b(3) = 0.0 
    711 !   p1 
     607      !   p1 
    712608      a(4) = 1.0 
    713609      b(4) = 0.0 
    714 !   s1 
     610      !   s1 
    715611      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 
    716612      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 
    717 !   k1 
     613      !   k1 
    718614      a(6) =1.0004 -0.0373*w1+ 0.0002*w2 
    719615      b(6) = -0.0374*w4 
    720 !  2n2 
     616      !  2n2 
    721617      a(7) = a(6) 
    722618      b(7) = b(6) 
    723 !  mu2 
     619      !  mu2 
    724620      a(8) = a(6) 
    725621      b(8) = b(6) 
    726 !   n2 
     622      !   n2 
    727623      a(9) = a(6) 
    728624      b(9) = b(6) 
    729 !  nu2 
     625      !  nu2 
    730626      a(10) = a(6) 
    731627      b(10) = b(6) 
    732 !   m2 
    733       a(11) = sqrt(w7*w7+w8*w8) 
    734       b(11) = atan(w8/w7) 
    735       if(w7.lt.0) b(11) = b(11) + 3.141992 
    736 !   l2 
     628      !   m2 
     629      a(11) = SQRT( w7 * w7 + w8 * w8 ) 
     630      b(11) = ATAN( w8 / w7 ) 
     631!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ???? 
     632      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992 
     633      !   l2 
    737634      a(12) = 1.0 
    738635      b(12) = 0.0 
    739 !   t2 
     636      !   t2 
    740637      a(13)= a(12) 
    741638      b(13)= b(12) 
    742 !   s2 
     639      !   s2 
    743640      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 
    744641      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 
    745 !   k2 
     642      !   k2 
    746643      a(15) = a(6)*a(6) 
    747644      b(15) = 2*b(6) 
    748 !   m4 
    749 ! 
    750       do 40 k = 1,ndc 
    751       b(k) = b(k)/rad 
    752 32    if (b(k)) 34,35,35 
    753 34    b(k) = b(k) + 360.0 
    754       go to 32 
    755 35    if (b(k)-360.0) 40,37,37 
    756 37    b(k) = b(k)-360.0 
    757       go to 35 
     645      !   m4 
     646!!gm  old coding,  remove GOTO and label of lines 
     647!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
     648      DO 40 k = 1,ndc 
     649         b(k) = b(k)/rad 
     65032       if (b(k)) 34,35,35 
     65134       b(k) = b(k) + 360.0 
     652         go to 32 
     65335       if (b(k)-360.0) 40,37,37 
     65437       b(k) = b(k)-360.0 
     655         go to 35 
    75865640    continue 
    759       RETURN 
    760       END SUBROUTINE ufset 
    761  
    762       SUBROUTINE vset( s,h,p,en,p1,v) 
    763       !!------------------------------------------------------------------------------ 
    764       !!                      SUBROUTINE vset 
    765       !!                     **************** 
     657      ! 
     658   END SUBROUTINE ufset 
     659    
     660 
     661   SUBROUTINE vset( s,h,p,en,p1,v) 
     662      !!---------------------------------------------------------------------- 
     663      !!                    ***  SUBROUTINE vset  *** 
     664      !!                    
    766665      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 
    767666      !!                 
    768       !! 
    769       !! History : 
    770       !! 
    771       !! Origins POLCOMS v6.3 2007 
    772       !!   NEMO v2.3  !  Jason Holt 
    773       !!------------------------------------------------------------------------------ 
    774       implicit none 
     667      !!---------------------------------------------------------------------- 
     668!!gm  doctor naming of dummy argument variables!!!   and all local variables 
     669!!gm  nc is jptides_max ???? 
     670!!gm  en argument is not used: suppress it ? 
    775671      integer nc 
    776672      parameter (nc=15) 
    777 !     subroutine arguments 
    778673      real(wp) s,h,p,en,p1 
    779674      real(wp)   v(nc) 
    780 ! 
    781       integer ndc,k 
    782       !!------------------------------------------------------------------------------ 
    783       !!  NEMO 2.3, LODYC-IPSL (2008) 
    784       !!------------------------------------------------------------------------------ 
     675      !! 
     676      integer ndc, k 
     677      !!---------------------------------------------------------------------- 
    785678 
    786679      ndc = nc 
    787 !   v s  are computed here. 
    788       v(1) =-3*s +h +p +270 
    789 !   q1 
    790       v(2) =-2*s +h +270.0 
    791 !   o1 
    792       v(3) =-h +270 
    793 !   p1 
    794       v(4) =180 
    795 !   s1 
    796       v(5) =h +90.0 
    797 !   k1 
    798       v(6) =-4*s +2*h +2*p 
    799 !   2n2 
    800       v(7) =-4*(s-h) 
    801 !   mu2 
    802       v(8) =-3*s +2*h +p 
    803 !   n2 
    804       v(9) =-3*s +4*h -p 
    805 !   mu2 
    806       v(10) =-2*s +2*h 
    807 !   m2 
    808       v(11) =-s +2*h -p +180 
    809 !   l2 
    810       v(12) =-h +p1 
    811 !   t2 
    812       v(13) =0 
    813 !   s2 
    814       v(14) =h+h 
    815 !   k2 
    816       v(15) =2*v(10) 
    817 !   m4 
    818 ! 
     680      !   v s  are computed here. 
     681      v(1) =-3*s +h +p +270      ! Q1 
     682      v(2) =-2*s +h +270.0       ! O1 
     683      v(3) =-h +270              ! P1 
     684      v(4) =180                  ! S1 
     685      v(5) =h +90.0              ! K1 
     686      v(6) =-4*s +2*h +2*p       ! 2N2 
     687      v(7) =-4*(s-h)             ! MU2 
     688      v(8) =-3*s +2*h +p         ! N2 
     689      v(9) =-3*s +4*h -p         ! MU2 
     690      v(10) =-2*s +2*h           ! M2 
     691      v(11) =-s +2*h -p +180     ! L2  
     692      v(12) =-h +p1              ! T2 
     693      v(13) =0                   ! S2 
     694      v(14) =h+h                 ! K2 
     695      v(15) =2*v(10)             ! M4 
     696! 
     697!!gm  old coding,  remove GOTO and label of lines 
     698!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    819699      do 72 k = 1, ndc 
    820 69    if (v(k) ) 70,71,71 
    821 70    v(k)=v(k)+360.0 
    822       go to 69 
    823 71    if ( v(k) - 360.0) 72,73,73 
    824 73    v(k)=v(k)-360.0 
    825       go to 71 
     70069       if( v(k) )  70,71,71 
     70170       v(k) = v(k)+360.0 
     702         go to 69 
     70371       if( v(k) - 360.0 )  72,73,73 
     70473       v(k) = v(k)-360.0 
     705         go to 71 
    82670672    continue 
    827       RETURN 
    828       END SUBROUTINE vset 
     707      ! 
     708   END SUBROUTINE vset 
    829709 
    830710#else 
    831    !!================================================================================= 
    832    !!                       ***  MODULE  bdytides *** 
    833    !!================================================================================= 
    834  
    835    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy_tides = .FALSE. !: tidal forcing at boundaries.  
    836  
    837    CHARACTER(len=80), PUBLIC :: & 
    838       filtide           !: Filename root for tidal input files 
    839  
    840    CHARACTER(len=4), PUBLIC, DIMENSION(1) :: & 
    841       tide_cpt         !: Names of tidal components used. 
     711   !!---------------------------------------------------------------------- 
     712   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides 
     713   !!---------------------------------------------------------------------- 
     714!!gm  are you sure we need to define filtide and tide_cpt ? 
     715   CHARACTER(len=80), PUBLIC               ::   filtide                !: Filename root for tidal input files 
     716   CHARACTER(len=4 ), PUBLIC, DIMENSION(1) ::   tide_cpt               !: Names of tidal components used. 
    842717 
    843718CONTAINS 
    844  
    845    SUBROUTINE tide_init 
    846                               ! No tidal forcing at boundaries ==> empty routine 
     719   SUBROUTINE tide_init                ! Empty routine 
    847720   END SUBROUTINE tide_init 
    848  
    849    SUBROUTINE tide_data 
    850                               ! No tidal forcing at boundaries ==> empty routine 
     721   SUBROUTINE tide_data                ! Empty routine 
    851722   END SUBROUTINE tide_data 
    852  
    853    SUBROUTINE tide_update( kt, jit ) 
    854       INTEGER :: kt, jit 
    855       WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit 
    856  
    857                               ! No tidal forcing at boundaries ==> empty routine 
     723   SUBROUTINE tide_update( kt, kit )   ! Empty routine 
     724      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit 
    858725   END SUBROUTINE tide_update 
    859  
    860726#endif 
    861727 
     728   !!====================================================================== 
    862729END MODULE bdytides 
  • trunk/NEMO/OPA_SRC/BDY/bdytra.F90

    r911 r1125  
    11MODULE bdytra 
    2    !!================================================================================= 
     2   !!====================================================================== 
    33   !!                       ***  MODULE  bdytra  *** 
    44   !! Ocean tracers:   Flow Relaxation Scheme of tracers on each open boundary 
    5    !!================================================================================= 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
     7   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     8   !!---------------------------------------------------------------------- 
    69#if defined key_bdy 
    7    !!--------------------------------------------------------------------------------- 
    8    !!   'key_bdy'      :                         Unstructured Open Boundary Conditions 
    9    !!--------------------------------------------------------------------------------- 
     10   !!---------------------------------------------------------------------- 
     11   !!   'key_bdy'                     Unstructured Open Boundary Conditions 
     12   !!---------------------------------------------------------------------- 
    1013   !!   bdy_tra        : Relaxation of tracers on unstructured open boundaries 
    11    !!--------------------------------------------------------------------------------- 
    12    !! * Modules used 
     14   !!---------------------------------------------------------------------- 
    1315   USE oce             ! ocean dynamics and tracers variables 
    1416   USE dom_oce         ! ocean space and time domain variables  
    1517   USE bdy_oce         ! ocean open boundary conditions 
    1618   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     19   USE in_out_manager  ! I/O manager 
    1720 
    1821   IMPLICIT NONE 
    1922   PRIVATE 
    2023 
    21    !! * Accessibility 
    2224   PUBLIC bdy_tra     ! routine called in tranxt.F90  
    2325 
    24    !! * Substitutions 
    25  
    26    !!--------------------------------------------------------------------------------- 
    27    !!   OPA 9.0 , LODYC-IPSL  (2003) 
    28    !!--------------------------------------------------------------------------------- 
     26   !!---------------------------------------------------------------------- 
     27   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     28   !! $Id: $  
     29   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     30   !!---------------------------------------------------------------------- 
    2931 
    3032CONTAINS 
    3133 
    3234   SUBROUTINE bdy_tra( kt ) 
    33       !!------------------------------------------------------------------------------ 
     35      !!---------------------------------------------------------------------- 
    3436      !!                 ***  SUBROUTINE bdy_tra  *** 
    3537      !!                     
     
    3739      !!              case of unstructured open boundaries. 
    3840      !!  
    39       !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in  
    40       !!             a three-dimensional baroclinic ocean model with realistic 
    41       !!             topography. Tellus, 365-382. 
    42       !!  History : 
    43       !!   9.0  !  05-01 (J. Chanut, A. Sellar) Original 
    44       !!------------------------------------------------------------------------------ 
    45       !! * Arguments 
     41      !! Reference : Engedahl H., 1995, Tellus, 365-382. 
     42      !!---------------------------------------------------------------------- 
    4643      INTEGER, INTENT( in ) ::   kt 
     44      !!  
     45      REAL(wp) ::   zwgt           ! boundary weight 
     46      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     47      INTEGER  ::   ii, ij         ! 2D addresses 
     48      !!---------------------------------------------------------------------- 
     49      ! 
     50      IF(ln_bdy_tra_frs) THEN ! If this is false, then this routine does nothing.  
    4751 
    48       !! * Local declarations 
    49       REAL(wp) :: zwgt                       ! boundary weight 
    50       INTEGER ::   jb, jk, jgrd              ! dummy loop indices 
    51       INTEGER ::   ii, ij                    ! 2D addresses 
    52       !!------------------------------------------------------------------------------ 
    53  
    54       jgrd=1 !: Everything is at T-points here 
    55   
    56       DO jb = 1, nblen(jgrd) 
    57         DO jk = 1, jpkm1 
    58           ii = nbi(jb,jgrd) 
    59           ij = nbj(jb,jgrd) 
    60           zwgt = nbw(jb,jgrd) 
    61  
    62           ! Temperature relaxation at the boundary    
    63           ta(ii,ij,jk) = ( ta(ii,ij,jk)*(1.-zwgt)  +  tbdy(jb,jk)*zwgt ) & 
    64                                                           * tmask(ii,ij,jk)          
    65  
    66           ! Salinity relaxation at the boundary    
    67           sa(ii,ij,jk) = ( sa(ii,ij,jk)*(1.-zwgt)  +  sbdy(jb,jk)*zwgt ) & 
    68                                                           * tmask(ii,ij,jk) 
    69     
     52      IF( kt == nit000 ) THEN 
     53         IF(lwp) WRITE(numout,*) 
     54         IF(lwp) WRITE(numout,*) 'bdy_tra : Flow Relaxation Scheme for tracers' 
     55         IF(lwp) WRITE(numout,*) '~~~~~~~' 
     56      ENDIF 
     57      ! 
     58      igrd = 1                       ! Everything is at T-points here 
     59      DO ib = 1, nblen(igrd) 
     60         DO ik = 1, jpkm1 
     61            ii = nbi(ib,igrd) 
     62            ij = nbj(ib,igrd) 
     63            zwgt = nbw(ib,igrd) 
     64            ta(ii,ij,ik) = ( ta(ii,ij,ik) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)          
     65            sa(ii,ij,ik) = ( sa(ii,ij,ik) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik) 
    7066        END DO 
    7167      END DO  
    72  
    73       CALL lbc_lnk( ta, 'T', 1. ) ! Boundary points should be updated 
    74       CALL lbc_lnk( sa, 'T', 1. ) ! 
     68      ! 
     69      CALL lbc_lnk( ta, 'T', 1. )   ! Boundary points should be updated 
     70      CALL lbc_lnk( sa, 'T', 1. )   ! 
     71      ! 
     72      ENDIF ! ln_bdy_tra_frs 
    7573 
    7674   END SUBROUTINE bdy_tra 
     75    
    7776#else 
    78    !!--------------------------------------------------------------------------------- 
    79    !!   Default option                                                    Empty module 
    80    !!--------------------------------------------------------------------------------- 
     77   !!---------------------------------------------------------------------- 
     78   !!   Dummy module                   NO Unstruct Open Boundary Conditions 
     79   !!---------------------------------------------------------------------- 
    8180CONTAINS 
    82    SUBROUTINE bdy_tra      ! Empty routine 
     81   SUBROUTINE bdy_tra(kt)      ! Empty routine 
     82      WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt 
    8383   END SUBROUTINE bdy_tra 
    8484#endif 
    8585 
    86    !!================================================================================= 
     86   !!====================================================================== 
    8787END MODULE bdytra 
  • trunk/NEMO/OPA_SRC/BDY/bdyvol.F90

    r911 r1125  
    11MODULE bdyvol 
    2    !!================================================================================= 
     2   !!====================================================================== 
    33   !!                       ***  MODULE  bdyvol  *** 
    44   !! Ocean dynamic :  Volume constraint when unstructured boundary  
    55   !!                  and Free surface are used 
    6    !!================================================================================= 
     6   !!====================================================================== 
     7   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code 
     8   !!             -   !  2006-01  (J. Chanut) Bug correction 
     9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     10   !!---------------------------------------------------------------------- 
    711#if   defined key_bdy   &&   defined key_dynspg_flt 
    8    !!--------------------------------------------------------------------------------- 
    9    !!   'key_bdy'               and              unstructured open boundary conditions 
    10    !!   'key_dynspg_flt'                                  constant volume free surface 
    11    !!--------------------------------------------------------------------------------- 
    12    !! * Modules used 
     12   !!---------------------------------------------------------------------- 
     13   !!   'key_bdy'            and      unstructured open boundary conditions 
     14   !!   'key_dynspg_flt'                              filtered free surface 
     15   !!---------------------------------------------------------------------- 
    1316   USE oce             ! ocean dynamics and tracers  
    1417   USE dom_oce         ! ocean space and time domain  
     
    2225   PRIVATE 
    2326 
    24    !! * Accessibility 
    2527   PUBLIC bdy_vol        ! routine called by dynspg_flt.h90 
    2628 
    2729   !! * Substitutions 
    2830#  include "domzgr_substitute.h90" 
    29    !!--------------------------------------------------------------------------------- 
    30    !!   OPA 9.0 , LODYC-IPSL  (2003) 
    31    !!--------------------------------------------------------------------------------- 
     31   !!---------------------------------------------------------------------- 
     32   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     33   !! $Id: $  
     34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     35   !!---------------------------------------------------------------------- 
    3236 
    3337CONTAINS 
    3438 
    35    SUBROUTINE bdy_vol ( kt ) 
    36       !!------------------------------------------------------------------------------ 
     39   SUBROUTINE bdy_vol( kt ) 
     40      !!---------------------------------------------------------------------- 
    3741      !!                      ***  ROUTINE bdyvol  *** 
    3842      !! 
    39       !! ** Purpose :  
    40       !!      This routine is called in dynspg_flt to control  
     43      !! ** Purpose :   This routine is called in dynspg_flt to control  
    4144      !!      the volume of the system. A correction velocity is calculated 
    4245      !!      to correct the total transport through the unstructured OBC.  
     
    4447      !!      linear free surface coded in OPA 8.2 
    4548      !! 
    46       !! ** Method :   
    47       !!      The correction velocity (zubtpecor here) is defined calculating 
     49      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating 
    4850      !!      the total transport through all open boundaries (trans_bdy) minus 
    49       !!      the cumulate E-P flux (zCflxemp) divided by the total lateral  
     51      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral  
    5052      !!      surface (bdysurftot) of the unstructured boundary.  
    51       !! 
    52       !!      zubtpecor = [trans_bdy - zCflxemp ]*(1./bdysurftot) 
    53       !! 
    54       !!      with zCflxemp => sum of (Evaporation minus Precipitation) 
     53      !!         zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot) 
     54      !!      with z_cflxemp => sum of (Evaporation minus Precipitation) 
    5555      !!                       over all the domain in m3/s at each time step. 
    56       !! 
    57       !!      zCflxemp < 0 when precipitation dominate 
    58       !!      zCflxemp > 0 when evaporation dominate 
     56      !!      z_cflxemp < 0 when precipitation dominate 
     57      !!      z_cflxemp > 0 when evaporation dominate 
    5958      !! 
    6059      !!      There are 2 options (user's desiderata):  
    61       !! 
    6260      !!         1/ The volume changes according to E-P, this is the default 
    6361      !!            option. In this case the cumulate E-P flux are setting to 
    64       !!            zero (zCflxemp=0) to calculate the correction velocity. So 
     62      !!            zero (z_cflxemp=0) to calculate the correction velocity. So 
    6563      !!            it will only balance the flux through open boundaries. 
    6664      !!            (set volbdy to 0 in tne namelist for this option) 
    67       !! 
    6865      !!         2/ The volume is constant even with E-P flux. In this case 
    6966      !!            the correction velocity must balance both the flux  
     
    7168      !!            surface.  
    7269      !!            (set volbdy to 1 in tne namelist for this option) 
     70      !!---------------------------------------------------------------------- 
     71      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    7372      !! 
    74       !! History : 
    75       !!   8.5  !  05-01 (J. Chanut, A. Sellar) Original code 
    76       !!        !  06-01 (J. Chanut) Bug correction 
    77       !!---------------------------------------------------------------------------- 
    78       !! * Arguments 
    79       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    80  
    81       !! * Local declarations 
    82       INTEGER ::   ji,jj,jb, jgrd, jk 
    83       REAL(wp) ::   zubtpecor 
    84       REAL(wp) ::   zCflxemp 
    85       REAL(wp) ::   ztranst 
     73      INTEGER  ::   ji, jj, jk, jb, jgrd 
     74      INTEGER  ::   ii, ij 
     75      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst 
    8676      !!----------------------------------------------------------------------------- 
    8777 
    8878      IF( kt == nit000 ) THEN  
    89          IF(lwp) WRITE(numout,*)'        ' 
     79         IF(lwp) WRITE(numout,*) 
    9080         IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC' 
    9181         IF(lwp) WRITE(numout,*)'~~~~~~~' 
    92          IF(lwp) WRITE(numout,*)'        ' 
    9382      END IF  
    9483 
    95       ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain. 
    96       ! --------------------------------------------------------------------------- 
    97   
    98       zCflxemp = 0.e0 
    99  
     84      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
     85      ! ----------------------------------------------------------------------- 
     86      z_cflxemp = 0.e0 
    10087      DO jj = 1, jpj 
    10188         DO ji = 1, jpi 
    102             zCflxemp = zCflxemp + ( (emp(ji,jj)*bdytmask(ji,jj)*tmask_i(ji,jj) )/rauw) & 
    103                                   *e1v(ji,jj)*e2u(ji,jj) 
     89            z_cflxemp = z_cflxemp + emp(ji,jj) * bdytmask(ji,jj) * tmask_i(ji,jj) / rauw * e1v(ji,jj) * e2u(ji,jj) 
    10490         END DO 
    10591      END DO 
    106       IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain 
     92      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )   ! sum over the global domain 
    10793 
    108       ! 2. Barotropic velocity through the unstructured open boundary 
    109       ! ------------------------------------------------------------- 
    110  
     94      ! Barotropic velocity through the unstructured open boundary 
     95      ! ---------------------------------------------------------- 
    11196      zubtpecor = 0.e0 
    112  
    113       jgrd = 2 ! cumulate u component contribution first  
     97      jgrd = 2                               ! cumulate u component contribution first  
    11498      DO jb = 1, nblenrim(jgrd) 
    115         DO jk = 1, jpkm1 
    116           zubtpecor = zubtpecor + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 
    117                                             * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    118                                             * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    119         END DO 
     99         DO jk = 1, jpkm1 
     100            ii = nbi(jb,jgrd) 
     101            ij = nbj(jb,jgrd) 
     102            zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     103         END DO 
    120104      END DO 
    121  
    122       jgrd = 3 ! then add v component contribution 
    123        DO jb = 1, nblenrim(jgrd) 
    124         DO jk = 1, jpkm1 
    125           zubtpecor = zubtpecor + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 
    126                                             * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    127                                             * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk)  
    128         END DO 
     105      jgrd = 3                               ! then add v component contribution 
     106      DO jb = 1, nblenrim(jgrd) 
     107         DO jk = 1, jpkm1 
     108            ii = nbi(jb,jgrd) 
     109            ij = nbj(jb,jgrd) 
     110            zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
     111         END DO 
    129112      END DO 
    130  
    131113      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
    132114 
    133  
    134       ! 3. The normal velocity correction 
    135       ! --------------------------------- 
    136  
    137       IF (volbdy==1) THEN 
    138         zubtpecor = (zubtpecor - zCflxemp)*(1./bdysurftot)  
    139       ELSE 
    140         zubtpecor =  zubtpecor*(1./bdysurftot) 
     115      ! The normal velocity correction 
     116      ! ------------------------------ 
     117      IF (volbdy==1) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot  
     118      ELSE                  ;   zubtpecor =   zubtpecor             / bdysurftot 
    141119      END IF 
    142120 
     121      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 
     122      ! ------------------------------------------------------------- 
     123      ztranst = 0.e0 
     124      jgrd = 2                               ! correct u component 
     125      DO jb = 1, nblenrim(jgrd) 
     126         DO jk = 1, jpkm1 
     127            ii = nbi(jb,jgrd) 
     128            ij = nbj(jb,jgrd) 
     129            ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 
     130            ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     131         END DO 
     132      END DO 
     133      jgrd = 3                              ! correct v component 
     134      DO jb = 1, nblenrim(jgrd) 
     135         DO jk = 1, jpkm1 
     136            ii = nbi(jb,jgrd) 
     137            ij = nbj(jb,jgrd) 
     138            va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
     139            ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
     140         END DO 
     141      END DO 
     142      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
     143  
     144      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 
     145      ! ------------------------------------------------------ 
     146 
    143147      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
    144          IF(lwp) WRITE(numout,*)'        ' 
     148         IF(lwp) WRITE(numout,*) 
    145149         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt 
    146150         IF(lwp) WRITE(numout,*)'~~~~~~~ ' 
    147          IF(lwp) WRITE(numout,*)'        ' 
    148          IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)' 
    149          IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',bdysurftot,'(m2)' 
    150          IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)' 
    151          IF(lwp) WRITE(numout,*)'        ' 
     151         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)' 
     152         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', bdysurftot, '(m2)' 
     153         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)' 
     154         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)' 
    152155      END IF  
    153  
    154       ! 4. Correction of the total velocity on the unstructured  
    155       !    boundary to respect the mass flux conservation 
    156       ! ------------------------------------------------------- 
    157  
    158       ztranst = 0.e0 
    159  
    160       jgrd = 2 ! correct u component 
    161       DO jb = 1, nblenrim(jgrd) 
    162         DO jk = 1, jpkm1 
    163           ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) = ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) & 
    164                             -flagu(jb) * zubtpecor * umask(nbi(jb, jgrd), nbj(jb, jgrd), jk) 
    165           ztranst = ztranst + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 
    166                                         * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    167                                         * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    168         END DO 
    169       END DO 
    170  
    171       jgrd = 3 ! correct v component 
    172        DO jb = 1, nblenrim(jgrd) 
    173         DO jk = 1, jpkm1 
    174           va(nbi(jb, jgrd), nbj(jb, jgrd), jk) = va(nbi(jb, jgrd), nbj(jb, jgrd), jk) & 
    175                             -flagv(jb) * zubtpecor * vmask(nbi(jb, jgrd), nbj(jb, jgrd), jk) 
    176           ztranst = ztranst + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 
    177                                         * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 
    178                                         * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
    179         END DO 
    180       END DO 
    181        
    182       IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
    183   
    184       ! 5. Check the cumulate transport through unstructured OBC 
    185       !    once barotropic velocities corrected 
    186       ! -------------------------------------------------------- 
    187  
    188       IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
    189          IF(lwp) WRITE(numout,*)'        ' 
    190          IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)' 
    191          IF(lwp) WRITE(numout,*)'        ' 
    192       END IF  
    193  
     156      ! 
    194157   END SUBROUTINE bdy_vol 
    195158 
    196159#else 
    197    !!--------------------------------------------------------------------------------- 
    198    !!  Default option :                                                   Empty module 
    199    !!--------------------------------------------------------------------------------- 
     160   !!---------------------------------------------------------------------- 
     161   !!   Dummy module                   NO Unstruct Open Boundary Conditions 
     162   !!---------------------------------------------------------------------- 
    200163CONTAINS 
    201    SUBROUTINE bdy_vol        ! Empty routine 
     164   SUBROUTINE bdy_vol( kt )        ! Empty routine 
     165      WRITE(*,*) 'bdy_vol: You should not have seen this print! error?', kt 
    202166   END SUBROUTINE bdy_vol 
    203167#endif 
    204168 
    205    !!================================================================================= 
     169   !!====================================================================== 
    206170END MODULE bdyvol 
  • trunk/NEMO/OPA_SRC/DOM/domvvl.F90

    r1058 r1125  
    302302#endif 
    303303 
    304 #if defined key_bdy || defined key_bdy_tides 
    305       DO jj = 1, jpj 
    306          DO ji = 1, jpi 
    307             zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 
    308          END DO 
    309       END DO 
     304#if defined key_bdy 
     305      zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) 
    310306#endif 
    311307 
  • trunk/NEMO/OPA_SRC/DYN/divcur.F90

    r1058 r1125  
    136136#endif 
    137137#endif          
    138 #if defined key_bdy || defined key_bdy_tides 
     138#if defined key_bdy 
    139139         ! unstructured open boundaries (div must be zero behind the open boundary) 
    140140         DO jj = 1, jpj 
     
    354354#endif 
    355355#endif          
    356 #if defined key_bdy || defined key_bdy_tides 
     356#if defined key_bdy 
    357357         ! unstructured open boundaries (div must be zero behind the open boundary) 
    358358         DO jj = 1, jpj 
  • trunk/NEMO/OPA_SRC/DYN/dynnxt.F90

    r911 r1125  
    238238      DO jk = 1, jpkm1                                 ! Horizontal slab 
    239239         !                                             ! =============== 
    240 # elif defined key_bdy || key_bdy_tides 
     240# elif defined key_bdy  
    241241         !                                             ! =============== 
    242242      END DO                                           !   End of slab 
  • trunk/NEMO/OPA_SRC/DYN/wzvmod.F90

    r1040 r1125  
    6161      REAL(wp) ::   z2dt, zraur          ! temporary scalar 
    6262      REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv 
    63 #if defined key_bdy || defined key_bdy_tides 
     63#if defined key_bdy 
    6464      INTEGER  ::     jgrd, jb           ! temporary scalars 
    6565#endif 
     
    113113#endif 
    114114 
    115 #if defined key_bdy || defined key_bdy_tides 
     115#if defined key_bdy 
    116116         jgrd=1 !: tracer grid. 
    117117         DO jb = 1, nblenrim(jgrd) 
  • trunk/NEMO/OPA_SRC/opa.F90

    r1037 r1125  
    5353   USE bdy_par         ! unstructured open boundary cond. parameters 
    5454   USE bdyini          ! unstructured open boundary cond. initialization (bdy_init routine) 
    55    USE bdytides        ! tides at open boundaries initialization (lk_bdy_tides) 
    5655   USE istate          ! initial state setting          (istate_init routine) 
    5756   USE eosbn2          ! equation of state            (eos bn2 routine) 
     
    272271      IF( lk_obc    )   CALL obc_init       ! Open boundaries  
    273272 
    274       IF( lk_bdy .OR. lk_bdy_tides )   & 
    275          &              CALL bdy_init       ! Unstructured open boundaries 
     273      IF( lk_bdy    )   CALL bdy_init       ! Unstructured open boundaries 
    276274 
    277275      CALL istate_init                      ! ocean initial state (Dynamics and tracers) 
Note: See TracChangeset for help on using the changeset viewer.