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

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

Changeset 3651


Ignore:
Timestamp:
2012-11-26T11:46:39+01:00 (12 years ago)
Author:
cbricaud
Message:

merge dev_MERCATOR_2012_rev3555 into dev_NOC_MERCATOR_2012 ; see ticket 1020

Location:
branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM
Files:
32 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/ARCH/arch-ifort_MERCATOR_CLUSTER.fcm

    r3294 r3651  
    1818%NCDF_LIB            -L $(NETCDF_LIB) -lnetcdf 
    1919%FC                  mpif90 
    20 %FCFLAGS             -assume byterecl  -convert big_endian -i4 -r8 -automatic -align all -O3 
     20%FCFLAGS             -assume byterecl  -convert big_endian -i4 -r8 -automatic -align all -O0 
    2121%FFLAGS              %FCFLAGS 
    2222%LD                  mpif90 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/CONFIG/AMM12/EXP00/namelist

    r3633 r3651  
    404404&nam_tide      !   tide parameters (#ifdef key_tide) 
    405405!----------------------------------------------------------------------- 
    406    ln_tide_pot   = .true.   !  use tidal potential forcing 
    407    nb_harmo      =    11    !  number of constituents used 
    408    clname(1)     =   'M2'   !  name of constituent 
    409    clname(2)     =   'S2' 
    410    clname(3)     =   'N2' 
    411    clname(4)     =   'K1' 
    412    clname(5)     =   'O1' 
    413    clname(6)     =   'Q1' 
    414    clname(7)     =   'M4' 
    415    clname(8)     =   'K2' 
    416    clname(9)     =   'P1' 
    417    clname(10)    =   'Mf' 
    418    clname(11)    =   'Mm' 
     406   ln_tide_pot   = .false.   !  use tidal potential forcing 
     407   clname(1)     =   'Q1'   !  name of constituent 
     408   clname(2)     =   'O1' 
     409   clname(3)     =   'P1' 
     410   clname(4)     =   'S1' 
     411   clname(5)     =   'K1' 
     412   clname(6)     =   '2N2' 
     413   clname(7)     =   'MU2' 
     414   clname(8)     =   'N2' 
     415   clname(9)     =   'NU2' 
     416   clname(10)    =   'M2' 
     417   clname(11)    =   'L2' 
     418   clname(12)    =   'T2' 
     419   clname(13)    =   'S2' 
     420   clname(14)    =   'K2' 
     421   clname(15)    =   'M4' 
    419422/ 
    420423!----------------------------------------------------------------------- 
    … …  
    460463!----------------------------------------------------------------------- 
    461464   filtide      = 'bdydta/amm12_bdytide_'         !  file name root of tidal forcing files 
    462     tide_cpt(1)   ='Q1'  !  names of tidal components used 
    463     tide_cpt(2)   ='O1'  !  names of tidal components used 
    464     tide_cpt(3)   ='P1'  !  names of tidal components used 
    465     tide_cpt(4)   ='S1'  !  names of tidal components used 
    466     tide_cpt(5)   ='K1'  !  names of tidal components used 
    467     tide_cpt(6)   ='2N2' !  names of tidal components used 
    468     tide_cpt(7)   ='MU2' !  names of tidal components used 
    469     tide_cpt(8)   ='N2'  !  names of tidal components used 
    470     tide_cpt(9)   ='NU2' !  names of tidal components used 
    471     tide_cpt(10)   ='M2'  !  names of tidal components used 
    472     tide_cpt(11)   ='L2'  !  names of tidal components used 
    473     tide_cpt(12)   ='T2'  !  names of tidal components used 
    474     tide_cpt(13)   ='S2'  !  names of tidal components used 
    475     tide_cpt(14)   ='K2'  !  names of tidal components used 
    476     tide_cpt(15)   ='M4'  !  names of tidal components used 
    477     tide_speed(1)   = 13.398661 !  phase speeds of tidal components (deg/hour) 
    478     tide_speed(2)   = 13.943036 !  phase speeds of tidal components (deg/hour) 
    479     tide_speed(3)   = 14.958932 !  phase speeds of tidal components (deg/hour) 
    480     tide_speed(4)   = 15.000001 !  phase speeds of tidal components (deg/hour) 
    481     tide_speed(5)   = 15.041069 !  phase speeds of tidal components (deg/hour) 
    482     tide_speed(6)   = 27.895355 !  phase speeds of tidal components (deg/hour) 
    483     tide_speed(7)   = 27.968210 !  phase speeds of tidal components (deg/hour) 
    484     tide_speed(8)   = 28.439730 !  phase speeds of tidal components (deg/hour) 
    485     tide_speed(9)   = 28.512585 !  phase speeds of tidal components (deg/hour) 
    486     tide_speed(10)   = 28.984106 !  phase speeds of tidal components (deg/hour) 
    487     tide_speed(11)   = 29.528479 !  phase speeds of tidal components (deg/hour) 
    488     tide_speed(12)   = 29.958935 !  phase speeds of tidal components (deg/hour) 
    489     tide_speed(13)   = 30.000002 !  phase speeds of tidal components (deg/hour) 
    490     tide_speed(14)   = 30.082138 !  phase speeds of tidal components (deg/hour) 
    491     tide_speed(15)   = 57.968212 !  phase speeds of tidal components (deg/hour) 
    492     ln_tide_date = .true.               !  adjust tidal harmonics for start date of run 
    493465/ 
    494466!!====================================================================== 
    … …  
    932904                           !     ln_ssh                  Logical switch for SSH observations               
    933905 
    934    ln_sst     = .false.    ! Logical switch for SST observations               
    935                            !     ln_reysst               Logical switch for Reynolds observations        
    936                            !     ln_ghrsst               Logical switch for GHRSST observations           
     906   ln_sst     = .true.     ! Logical switch for SST observations 
     907   ln_reysst  = .true.     !     ln_reysst               Logical switch for Reynolds observations 
     908   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations       
    937909 
    938910   ln_sstfb   = .false.    ! Logical switch for feedback SST data           

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist

    r3625 r3651  
    931931                           !     ln_ssh                  Logical switch for SSH observations               
    932932 
    933    ln_sst     = .false.    ! Logical switch for SST observations               
    934                            !     ln_reysst               Logical switch for Reynolds observations        
    935                            !     ln_ghrsst               Logical switch for GHRSST observations           
     933   ln_sst     = .true.     ! Logical switch for SST observations 
     934   ln_reysst  = .true.     !     ln_reysst               Logical switch for Reynolds observations 
     935   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations       
    936936 
    937937   ln_sstfb   = .false.    ! Logical switch for feedback SST data           

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/CONFIG/GYRE/EXP00/namelist

    r3625 r3651  
    926926                           !     ln_ssh                  Logical switch for SSH observations               
    927927 
    928    ln_sst     = .false.    ! Logical switch for SST observations               
    929                            !     ln_reysst               Logical switch for Reynolds observations        
    930                            !     ln_ghrsst               Logical switch for GHRSST observations           
     928   ln_sst     = .true.     ! Logical switch for SST observations               
     929   ln_reysst  = .true.     !     ln_reysst               Logical switch for Reynolds observations        
     930   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations           
    931931 
    932932   ln_sstfb   = .false.    ! Logical switch for feedback SST data           

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3625 r3651  
    930930                           !     ln_ssh                  Logical switch for SSH observations               
    931931 
    932    ln_sst     = .false.    ! Logical switch for SST observations               
    933                            !     ln_reysst               Logical switch for Reynolds observations        
    934                            !     ln_ghrsst               Logical switch for GHRSST observations           
     932   ln_sst     = .true.     ! Logical switch for SST observations 
     933   ln_reysst  = .true.     !     ln_reysst               Logical switch for Reynolds observations 
     934   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations       
    935935 
    936936   ln_sstfb   = .false.    ! Logical switch for feedback SST data           

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist

    r3625 r3651  
    929929                           !     ln_ssh                  Logical switch for SSH observations               
    930930 
    931    ln_sst     = .false.    ! Logical switch for SST observations               
    932                            !     ln_reysst               Logical switch for Reynolds observations        
    933                            !     ln_ghrsst               Logical switch for GHRSST observations           
     931   ln_sst     = .true.     ! Logical switch for SST observations 
     932   ln_reysst  = .true.     !     ln_reysst               Logical switch for Reynolds observations 
     933   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations       
    934934 
    935935   ln_sstfb   = .false.    ! Logical switch for feedback SST data           

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r3294 r3651  
    2828      INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap 
    2929      REAL   , POINTER, DIMENSION(:,:)   ::  nbw 
     30      REAL   , POINTER, DIMENSION(:,:)   ::  nbd 
    3031      REAL   , POINTER, DIMENSION(:)     ::  flagu 
    3132      REAL   , POINTER, DIMENSION(:)     ::  flagv 
    … …  
    7374   INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta             !: = 0 use the initial state as bdy dta ;  
    7475                                                            !: = 1 read it in a NetCDF file 
     76   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
     77   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
     78   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
     79 
    7580#if defined key_lim2 
    7681   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables  
    … …  
    101106   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions 
    102107                                                                          !: =1 => some data to be read in from data files 
    103    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays 
     108   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
     109   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    104110   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
    105111   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r3294 r3651  
    3232   USE ice_2 
    3333#endif 
     34   USE sbcapr 
    3435 
    3536   IMPLICIT NONE 
    … …  
    108109 
    109110            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN  
    110                IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
    111                   ilen1(:) = nblen(:) 
    112                ELSE 
    113                   ilen1(:) = nblenrim(:) 
    114                ENDIF 
     111               ilen1(:) = nblen(:) 
    115112               igrd = 1 
    116113               DO ib = 1, ilen1(igrd) 
    … …  
    134131 
    135132            IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN  
    136                IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
    137                   ilen1(:) = nblen(:) 
    138                ELSE 
    139                   ilen1(:) = nblenrim(:) 
    140                ENDIF 
     133               ilen1(:) = nblen(:) 
    141134               igrd = 2  
    142135               DO ib = 1, ilen1(igrd) 
    … …  
    158151 
    159152            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN  
    160                IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
    161                   ilen1(:) = nblen(:) 
    162                ELSE 
    163                   ilen1(:) = nblenrim(:) 
    164                ENDIF 
     153               ilen1(:) = nblen(:) 
    165154               igrd = 1                       ! Everything is at T-points here 
    166155               DO ib = 1, ilen1(igrd) 
    … …  
    176165#if defined key_lim2 
    177166            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
    178                IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
    179                   ilen1(:) = nblen(:) 
    180                ELSE 
    181                   ilen1(:) = nblenrim(:) 
    182                ENDIF 
     167               ilen1(:) = nblen(:) 
    183168               igrd = 1                       ! Everything is at T-points here 
    184169               DO ib = 1, ilen1(igrd) 
    … …  
    207192            IF( PRESENT(jit) ) THEN 
    208193               ! Update barotropic boundary conditions only 
    209                ! jit is optional argument for fld_read and tide_update 
     194               ! jit is optional argument for fld_read and bdytide_update 
    210195               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
    211196                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    … …  
    214199                     dta_bdy(ib_bdy)%v2d(:) = 0.0 
    215200                  ENDIF 
    216                   IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data 
    217                      jend = jstart + 2 
    218                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),   & 
    219                      &              jit=jit, time_offset=time_offset ) 
    220                   ENDIF 
    221                   IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    222                      CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
    223                      &                 jit=jit, time_offset=time_offset ) 
     201                  IF (nn_tra(ib_bdy).ne.4) THEN 
     202                     IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN 
     203                        ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 
     204                        jend = nb_bdy_fld(ib_bdy) 
     205                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 
     206                        CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 
     207                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 
     208                        ! If full velocities in boundary data then split into barotropic and baroclinic data 
     209                        IF( ln_full_vel_array(ib_bdy) .and.                                             & 
     210                       &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 
     211                           igrd = 2                      ! zonal velocity 
     212                           dta_bdy(ib_bdy)%u2d(:) = 0.0 
     213                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     214                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     215                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     216                              DO ik = 1, jpkm1 
     217                                 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
     218                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     219                              END DO 
     220                              dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     221                              DO ik = 1, jpkm1 
     222                                 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
     223                              END DO 
     224                           END DO 
     225                           igrd = 3                      ! meridional velocity 
     226                           dta_bdy(ib_bdy)%v2d(:) = 0.0 
     227                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     228                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     229                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     230                              DO ik = 1, jpkm1 
     231                                 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
     232                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     233                              END DO 
     234                              dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     235                              DO ik = 1, jpkm1 
     236                                 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
     237                              END DO 
     238                           END DO 
     239                        ENDIF                     
     240                     ENDIF 
     241                     IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     242                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
     243                          &                 jit=jit, time_offset=time_offset ) 
     244                     ENDIF 
    224245                  ENDIF 
    225246               ENDIF 
    226247            ELSE 
    227                IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    228                   dta_bdy(ib_bdy)%ssh(:) = 0.0 
    229                   dta_bdy(ib_bdy)%u2d(:) = 0.0 
    230                   dta_bdy(ib_bdy)%v2d(:) = 0.0 
     248               IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition 
     249                  jend = nb_bdy_fld(ib_bdy) 
     250                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 
     251                  ! 
     252                  igrd = 2                      ! zonal velocity 
     253                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     254                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     255                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     256                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
     257                  END DO 
     258                  ! 
     259                  igrd = 3                      ! meridional velocity 
     260                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     261                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     262                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     263                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
     264                  END DO 
     265               ELSE 
     266                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     267                     dta_bdy(ib_bdy)%ssh(:) = 0.0 
     268                     dta_bdy(ib_bdy)%u2d(:) = 0.0 
     269                     dta_bdy(ib_bdy)%v2d(:) = 0.0 
     270                  ENDIF 
     271                  IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
     272                     jend = nb_bdy_fld(ib_bdy) 
     273                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 
     274                  ENDIF 
     275                  ! If full velocities in boundary data then split into barotropic and baroclinic data 
     276                  IF( ln_full_vel_array(ib_bdy) .and.                                             & 
     277                 &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 
     278                     igrd = 2                      ! zonal velocity 
     279                     dta_bdy(ib_bdy)%u2d(:) = 0.0 
     280                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     281                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     282                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     283                        DO ik = 1, jpkm1 
     284                           dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
     285                 &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     286                        END DO 
     287                        dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     288                        DO ik = 1, jpkm1 
     289                           dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
     290                        END DO 
     291                     END DO 
     292                     igrd = 3                      ! meridional velocity 
     293                     dta_bdy(ib_bdy)%v2d(:) = 0.0 
     294                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     295                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     296                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     297                        DO ik = 1, jpkm1 
     298                           dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
     299                 &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     300                        END DO 
     301                        dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     302                        DO ik = 1, jpkm1 
     303                           dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
     304                        END DO 
     305                     END DO 
     306                  ENDIF 
     307                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     308                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
     309                  ENDIF 
    231310               ENDIF 
    232                IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
    233                   jend = jstart + nb_bdy_fld(ib_bdy) - 1 
    234                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 
    235                ENDIF 
    236                IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    237                   CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
    238                ENDIF 
    239311            ENDIF 
    240312            jstart = jend+1 
    241  
    242             ! If full velocities in boundary data then split into barotropic and baroclinic data 
    243             ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same 
    244             ! time as the dynspg_ts option).  
    245  
    246             IF( ln_full_vel_array(ib_bdy) .and.                                             &  
    247            &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN  
    248  
    249                igrd = 2                      ! zonal velocity 
    250                dta_bdy(ib_bdy)%u2d(:) = 0.0 
    251                DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     313         END IF ! nn_dta(ib_bdy) = 1 
     314      END DO  ! ib_bdy 
     315 
     316      IF ( ln_apr_obc ) THEN 
     317         DO ib_bdy = 1, nb_bdy 
     318            IF (nn_tra(ib_bdy).NE.4)THEN 
     319               igrd = 1                      ! meridional velocity 
     320               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    252321                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    253322                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    254                   DO ik = 1, jpkm1 
    255                      dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
    256               &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
    257                   END DO 
    258                   dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
    259                   DO ik = 1, jpkm1 
    260                      dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)  
    261                   END DO 
    262                END DO 
    263  
    264                igrd = 3                      ! meridional velocity 
    265                dta_bdy(ib_bdy)%v2d(:) = 0.0 
    266                DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    267                   ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    268                   ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    269                   DO ik = 1, jpkm1 
    270                      dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
    271               &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
    272                   END DO 
    273                   dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
    274                   DO ik = 1, jpkm1 
    275                      dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)  
    276                   END DO 
    277                END DO 
    278      
    279             ENDIF 
    280  
    281          END IF ! nn_dta(ib_bdy) = 1 
    282       END DO  ! ib_bdy 
     323                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij) 
     324               ENDDO 
     325            ENDIF 
     326         ENDDO 
     327      ENDIF 
    283328 
    284329      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta') 
    … …  
    326371      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init') 
    327372 
     373      IF(lwp) WRITE(numout,*) 
     374      IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries' 
     375      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     376      IF(lwp) WRITE(numout,*) '' 
     377 
    328378      ! Set nn_dta 
    329379      DO ib_bdy = 1, nb_bdy 
    … …  
    357407         ENDIF 
    358408#endif                
     409         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy) 
    359410      ENDDO             
    360411 
    … …  
    408459            ln_full_vel_array(ib_bdy) = ln_full_vel 
    409460 
    410             IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN 
    411                CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',& 
    412             &                  'with dynspg_ts option' )   ;   RETURN   
    413             ENDIF              
    414  
    415461            nblen => idx_bdy(ib_bdy)%nblen 
    416462            nblenrim => idx_bdy(ib_bdy)%nblenrim 
    … …  
    420466            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
    421467 
    422                IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
     468               IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
    423469                  jfld = jfld + 1 
    424470                  blf_i(jfld) = bn_ssh 
    425471                  ibdy(jfld) = ib_bdy 
    426472                  igrid(jfld) = 1 
    427                   ilen1(jfld) = nblenrim(igrid(jfld)) 
     473                  ilen1(jfld) = nblen(igrid(jfld)) 
    428474                  ilen3(jfld) = 1 
    429475               ENDIF 
    430476 
    431477               IF( .not. ln_full_vel_array(ib_bdy) ) THEN 
    432  
    433478                  jfld = jfld + 1 
    434479                  blf_i(jfld) = bn_u2d 
    435480                  ibdy(jfld) = ib_bdy 
    436481                  igrid(jfld) = 2 
    437                   IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
    438                      ilen1(jfld) = nblen(igrid(jfld)) 
    439                   ELSE 
    440                      ilen1(jfld) = nblenrim(igrid(jfld)) 
    441                   ENDIF 
     482                  ilen1(jfld) = nblen(igrid(jfld)) 
    442483                  ilen3(jfld) = 1 
    443484 
    … …  
    446487                  ibdy(jfld) = ib_bdy 
    447488                  igrid(jfld) = 3 
    448                   IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
    449                      ilen1(jfld) = nblen(igrid(jfld)) 
    450                   ELSE 
    451                      ilen1(jfld) = nblenrim(igrid(jfld)) 
    452                   ENDIF 
     489                  ilen1(jfld) = nblen(igrid(jfld)) 
    453490                  ilen3(jfld) = 1 
    454  
    455491               ENDIF 
    456492 
    … …  
    466502               ibdy(jfld) = ib_bdy 
    467503               igrid(jfld) = 2 
    468                IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
    469                   ilen1(jfld) = nblen(igrid(jfld)) 
    470                ELSE 
    471                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    472                ENDIF 
     504               ilen1(jfld) = nblen(igrid(jfld)) 
    473505               ilen3(jfld) = jpk 
    474506 
    … …  
    477509               ibdy(jfld) = ib_bdy 
    478510               igrid(jfld) = 3 
    479                IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
    480                   ilen1(jfld) = nblen(igrid(jfld)) 
    481                ELSE 
    482                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    483                ENDIF 
     511               ilen1(jfld) = nblen(igrid(jfld)) 
    484512               ilen3(jfld) = jpk 
    485513 
    … …  
    493521               ibdy(jfld) = ib_bdy 
    494522               igrid(jfld) = 1 
    495                IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
    496                   ilen1(jfld) = nblen(igrid(jfld)) 
    497                ELSE 
    498                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    499                ENDIF 
     523               ilen1(jfld) = nblen(igrid(jfld)) 
    500524               ilen3(jfld) = jpk 
    501525 
    … …  
    504528               ibdy(jfld) = ib_bdy 
    505529               igrid(jfld) = 1 
    506                IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
    507                   ilen1(jfld) = nblen(igrid(jfld)) 
    508                ELSE 
    509                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    510                ENDIF 
     530               ilen1(jfld) = nblen(igrid(jfld)) 
    511531               ilen3(jfld) = jpk 
    512532 
    … …  
    521541               ibdy(jfld) = ib_bdy 
    522542               igrid(jfld) = 1 
    523                IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
    524                   ilen1(jfld) = nblen(igrid(jfld)) 
    525                ELSE 
    526                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    527                ENDIF 
     543               ilen1(jfld) = nblen(igrid(jfld)) 
    528544               ilen3(jfld) = 1 
    529545 
    … …  
    532548               ibdy(jfld) = ib_bdy 
    533549               igrid(jfld) = 1 
    534                IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
    535                   ilen1(jfld) = nblen(igrid(jfld)) 
    536                ELSE 
    537                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    538                ENDIF 
     550               ilen1(jfld) = nblen(igrid(jfld)) 
    539551               ilen3(jfld) = 1 
    540552 
    … …  
    543555               ibdy(jfld) = ib_bdy 
    544556               igrid(jfld) = 1 
    545                IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
    546                   ilen1(jfld) = nblen(igrid(jfld)) 
    547                ELSE 
    548                   ilen1(jfld) = nblenrim(igrid(jfld)) 
    549                ENDIF 
     557               ilen1(jfld) = nblen(igrid(jfld)) 
    550558               ilen3(jfld) = 1 
    551559 
    … …  
    566574      ENDDO ! ib_bdy 
    567575 
    568  
    569576      DO jfld = 1, nb_bdy_fld_sum 
    570577         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) 
    … …  
    577584      jstart = 1 
    578585      DO ib_bdy = 1, nb_bdy 
    579          jend = jstart + nb_bdy_fld(ib_bdy) - 1 
     586         jend = nb_bdy_fld(ib_bdy)  
    580587         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   & 
    581588         &              'open boundary conditions', 'nambdy_dta' ) 
    … …  
    596603         IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 
    597604            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 
    598                IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
    599                   ilen0(1:3) = nblen(1:3) 
    600                ELSE 
    601                   ilen0(1:3) = nblenrim(1:3) 
    602                ENDIF 
    603                ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) ) 
     605               ilen0(1:3) = nblen(1:3) 
    604606               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
    605607               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
     608               IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN 
     609                  jfld = jfld + 1 
     610                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
     611               ELSE 
     612                  ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) ) 
     613               ENDIF 
    606614            ELSE 
    607615               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
    … …  
    617625 
    618626         IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
    619             IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
    620                ilen0(1:3) = nblen(1:3) 
    621             ELSE 
    622                ilen0(1:3) = nblenrim(1:3) 
    623             ENDIF 
     627            ilen0(1:3) = nblen(1:3) 
    624628            ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 
    625629            ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) 
    … …  
    636640         IF (nn_tra(ib_bdy) .gt. 0) THEN 
    637641            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
    638                IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
    639                   ilen0(1:3) = nblen(1:3) 
    640                ELSE 
    641                   ilen0(1:3) = nblenrim(1:3) 
    642                ENDIF 
     642               ilen0(1:3) = nblen(1:3) 
    643643               ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 
    644644               ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) 
    … …  
    654654         IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
    655655            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
    656                IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
    657                   ilen0(1:3) = nblen(1:3) 
    658                ELSE 
    659                   ilen0(1:3) = nblenrim(1:3) 
    660                ENDIF 
     656               ilen0(1:3) = nblen(1:3) 
    661657               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 
    662658               ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r3294 r3651  
    1414   !!---------------------------------------------------------------------- 
    1515   USE timing          ! Timing 
     16   USE wrk_nemo        ! Memory Allocation 
    1617   USE oce             ! ocean dynamics and tracers  
    1718   USE dom_oce         ! ocean space and time domain 
    … …  
    1920   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2021   USE in_out_manager  ! 
     22   Use phycst 
    2123 
    2224   IMPLICIT NONE 
    … …  
    2426 
    2527   PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn 
    26  
     28   PUBLIC   bdy_dyn3d_dmp ! routine called by step 
     29 
     30   !! * Substitutions 
     31#  include "domzgr_substitute.h90" 
    2732   !!---------------------------------------------------------------------- 
    2833   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    5560         CASE(jp_frs) 
    5661            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     62         CASE(2) 
     63            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     64         CASE(3) 
     65            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    5766         CASE DEFAULT 
    5867            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 
    … …  
    6170 
    6271   END SUBROUTINE bdy_dyn3d 
     72 
     73   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt ) 
     74      !!---------------------------------------------------------------------- 
     75      !!                  ***  SUBROUTINE bdy_dyn3d_spe  *** 
     76      !! 
     77      !! ** Purpose : - Apply a specified value for baroclinic velocities 
     78      !!                at open boundaries. 
     79      !! 
     80      !!---------------------------------------------------------------------- 
     81      INTEGER                     ::   kt 
     82      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     83      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     84      !! 
     85      INTEGER  ::   jb, jk         ! dummy loop indices 
     86      INTEGER  ::   ii, ij, igrd   ! local integers 
     87      REAL(wp) ::   zwgt           ! boundary weight 
     88      !!---------------------------------------------------------------------- 
     89      ! 
     90      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe') 
     91      ! 
     92      igrd = 2                      ! Relaxation of zonal velocity 
     93      DO jb = 1, idx%nblenrim(igrd) 
     94         DO jk = 1, jpkm1 
     95            ii   = idx%nbi(jb,igrd) 
     96            ij   = idx%nbj(jb,igrd) 
     97            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk) 
     98         END DO 
     99      END DO 
     100      ! 
     101      igrd = 3                      ! Relaxation of meridional velocity 
     102      DO jb = 1, idx%nblenrim(igrd) 
     103         DO jk = 1, jpkm1 
     104            ii   = idx%nbi(jb,igrd) 
     105            ij   = idx%nbj(jb,igrd) 
     106            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk) 
     107         END DO 
     108      END DO 
     109      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
     110      ! 
     111      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     112 
     113      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe') 
     114 
     115   END SUBROUTINE bdy_dyn3d_spe 
     116 
     117   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt ) 
     118      !!---------------------------------------------------------------------- 
     119      !!                  ***  SUBROUTINE bdy_dyn3d_zro  *** 
     120      !! 
     121      !! ** Purpose : - baroclinic velocities = 0. at open boundaries. 
     122      !! 
     123      !!---------------------------------------------------------------------- 
     124      INTEGER                     ::   kt 
     125      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     126      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     127      !! 
     128      INTEGER  ::   ib, ik         ! dummy loop indices 
     129      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers 
     130      REAL(wp) ::   zwgt           ! boundary weight 
     131      !!---------------------------------------------------------------------- 
     132      ! 
     133      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro') 
     134      ! 
     135      igrd = 2                       ! Everything is at T-points here 
     136      DO ib = 1, idx%nblenrim(igrd) 
     137         ii = idx%nbi(ib,igrd) 
     138         ij = idx%nbj(ib,igrd) 
     139         DO ik = 1, jpkm1 
     140            ua(ii,ij,ik) = 0._wp 
     141         END DO 
     142      END DO 
     143 
     144      igrd = 3                       ! Everything is at T-points here 
     145      DO ib = 1, idx%nblenrim(igrd) 
     146         ii = idx%nbi(ib,igrd) 
     147         ij = idx%nbj(ib,igrd) 
     148         DO ik = 1, jpkm1 
     149            va(ii,ij,ik) = 0._wp 
     150         END DO 
     151      END DO 
     152      ! 
     153      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
     154      ! 
     155      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     156 
     157      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro') 
     158 
     159   END SUBROUTINE bdy_dyn3d_zro 
    63160 
    64161   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt ) 
    … …  
    111208   END SUBROUTINE bdy_dyn3d_frs 
    112209 
     210   SUBROUTINE bdy_dyn3d_dmp( kt ) 
     211      !!---------------------------------------------------------------------- 
     212      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  *** 
     213      !! 
     214      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries. 
     215      !! 
     216      !!---------------------------------------------------------------------- 
     217      INTEGER                     ::   kt 
     218      !! 
     219      INTEGER  ::   jb, jk         ! dummy loop indices 
     220      INTEGER  ::   ii, ij, igrd   ! local integers 
     221      REAL(wp) ::   zwgt           ! boundary weight 
     222      INTEGER  ::  ib_bdy          ! loop index 
     223      !!---------------------------------------------------------------------- 
     224      ! 
     225      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp') 
     226      ! 
     227      !------------------------------------------------------- 
     228      ! Remove barotropic part from before velocity 
     229      !------------------------------------------------------- 
     230      CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
     231 
     232      pu2d(:,:) = 0.e0 
     233      pv2d(:,:) = 0.e0 
     234 
     235      DO jk = 1, jpkm1 
     236#if defined key_vvl 
     237         pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk)  
     238         pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk) 
     239#else 
     240         pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk) 
     241         pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk) 
     242#endif 
     243      END DO 
     244 
     245      IF( lk_vvl ) THEN 
     246         pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     247         pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
     248      ELSE 
     249         pu2d(:,:) = pv2d(:,:) * hur(:,:) 
     250         pv2d(:,:) = pu2d(:,:) * hvr(:,:) 
     251      ENDIF 
     252 
     253      DO ib_bdy=1, nb_bdy 
     254         IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN 
     255            igrd = 2                      ! Relaxation of zonal velocity 
     256            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     257               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
     258               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
     259               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
     260               DO jk = 1, jpkm1 
     261                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk) 
     262               END DO 
     263            END DO 
     264            ! 
     265            igrd = 3                      ! Relaxation of meridional velocity 
     266            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     267               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
     268               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
     269               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
     270               DO jk = 1, jpkm1 
     271                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk) 
     272               END DO 
     273            END DO 
     274         ENDIF 
     275      ENDDO 
     276      ! 
     277      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     278      ! 
     279      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
     280      ! 
     281      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp') 
     282 
     283   END SUBROUTINE bdy_dyn3d_dmp 
    113284 
    114285#else 
    … …  
    118289CONTAINS 
    119290   SUBROUTINE bdy_dyn3d( kt )      ! Empty routine 
    120       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
     291      WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt 
    121292   END SUBROUTINE bdy_dyn3d 
     293 
     294   SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine 
     295      WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt 
     296   END SUBROUTINE bdy_dyn3d_dmp 
     297 
    122298#endif 
    123299 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r3632 r3651  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
    1212   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     13   !!            3.4  !  2012     (J. Chanut) straight open boundary case update 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_bdy 
    … …  
    2627   USE lib_mpp         ! for mpp_sum   
    2728   USE iom             ! I/O 
     29   USE sbctide, ONLY: lk_tide ! Tidal forcing or not 
     30   USE phycst, ONLY: rday 
    2831 
    2932   IMPLICIT NONE 
    … …  
    3235   PUBLIC   bdy_init   ! routine called in nemo_init 
    3336 
     37   INTEGER, PARAMETER          :: jp_nseg = 100 
     38   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured 
     39                                               ! open boundary data files 
     40   ! Straight open boundary segment parameters: 
     41   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
     42   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge 
     43   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw 
     44   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn 
     45   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs 
    3446   !!---------------------------------------------------------------------- 
    3547   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    … …  
    5365      ! namelist variables 
    5466      !------------------- 
    55       INTEGER, PARAMETER          :: jp_nseg = 100 
    56       INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
    57       INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft 
    58       INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft 
    59       INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft 
    60       INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft 
     67      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
     68      CHARACTER(LEN=1)   ::   ctypebdy 
     69      INTEGER :: nbdyind, nbdybeg, nbdyend 
    6170 
    6271      ! local variables 
    … …  
    6675      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       - 
    6776      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       - 
     77      INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
     78      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
    6879      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts 
    6980      REAL   , POINTER  ::  flagu, flagv                   !    -   - 
    7081      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    71       INTEGER, DIMENSION (2)                ::   kdimsz 
     82      INTEGER, DIMENSION (2)                  ::   kdimsz 
    7283      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
    7384      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    7485      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    75       REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask 
    76       CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
    77       CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid 
     86      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    7887      !! 
    7988      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             & 
    8089         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, & 
    8190         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &   
     91         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              & 
    8292#if defined key_lim2 
    8393         &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
    … …  
    8595         &             ln_vol, nn_volctl, nn_rimwidth 
    8696      !! 
    87       NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
    88                              nbdysegw, jpiwob, jpjwdt, jpjwft,             & 
    89                              nbdysegn, jpjnob, jpindt, jpinft,             & 
    90                              nbdysegs, jpjsob, jpisdt, jpisft 
     97      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
    9198 
    9299      !!---------------------------------------------------------------------- 
    … …  
    105112 
    106113      cgrid= (/'t','u','v'/) 
    107  
     114       
    108115      ! ----------------------------------------- 
    109116      ! Initialise and read namelist parameters 
    … …  
    121128      nn_tra(:)         = 0 
    122129      nn_tra_dta(:)     = -1  ! uninitialised flag 
     130      ln_tra_dmp(:)     = .false. 
     131      ln_dyn3d_dmp(:)   = .false. 
     132      rn_time_dmp(:)    = 1. 
    123133#if defined key_lim2 
    124134      nn_ice_lim2(:)    = 0 
    … …  
    135145      ! Check and write out namelist parameters 
    136146      ! ----------------------------------------- 
    137  
    138147      !                                   ! control prints 
    139148      IF(lwp) WRITE(numout,*) '         nambdy' 
    … …  
    158167        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    159168        SELECT CASE( nn_dyn2d(ib_bdy) )                   
    160           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    161           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    162           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
     169          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     170          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     171          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    163172          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
    164173        END SELECT 
    … …  
    171180              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    172181           END SELECT 
     182           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN 
     183             CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' ) 
     184           ENDIF 
    173185        ENDIF 
    174186        IF(lwp) WRITE(numout,*) 
    … …  
    176188        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    177189        SELECT CASE( nn_dyn3d(ib_bdy) )                   
    178           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    179           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     190          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     191          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     192          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
     193          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    180194          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 
    181195        END SELECT 
    … …  
    187201           END SELECT 
    188202        ENDIF 
     203 
     204        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
     205           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 
     206              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
     207              ln_dyn3d_dmp(ib_bdy)=.false. 
     208           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 
     209              CALL ctl_stop( 'Use FRS OR relaxation' ) 
     210           ELSE 
     211              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
     212              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     213              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     214           ENDIF 
     215        ELSE 
     216           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
     217        ENDIF 
    189218        IF(lwp) WRITE(numout,*) 
    190219 
    191220        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    192221        SELECT CASE( nn_tra(ib_bdy) )                   
    193           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    194           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     222          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     223          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     224          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
     225          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     226          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
    195227          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
    196228        END SELECT 
    … …  
    201233              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
    202234           END SELECT 
     235        ENDIF 
     236 
     237        IF ( ln_tra_dmp(ib_bdy) ) THEN 
     238           IF ( nn_tra(ib_bdy).EQ.0 ) THEN 
     239              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
     240              ln_tra_dmp(ib_bdy)=.false. 
     241           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 
     242              CALL ctl_stop( 'Use FRS OR relaxation' ) 
     243           ELSE 
     244              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
     245              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     246              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     247           ENDIF 
     248        ELSE 
     249           IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
    203250        ENDIF 
    204251        IF(lwp) WRITE(numout,*) 
    … …  
    221268#endif 
    222269 
    223         IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy) 
     270        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    224271        IF(lwp) WRITE(numout,*) 
    225272 
    226273      ENDDO 
    227274 
    228      IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    229        IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    230        IF(lwp) WRITE(numout,*) 
    231        SELECT CASE ( nn_volctl ) 
    232          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    233          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    234          CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    235        END SELECT 
    236        IF(lwp) WRITE(numout,*) 
    237      ELSE 
    238        IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    239        IF(lwp) WRITE(numout,*) 
     275     IF (nb_bdy .gt. 0) THEN 
     276        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     277          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
     278          IF(lwp) WRITE(numout,*) 
     279          SELECT CASE ( nn_volctl ) 
     280            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
     281            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
     282            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
     283          END SELECT 
     284          IF(lwp) WRITE(numout,*) 
     285        ELSE 
     286          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
     287          IF(lwp) WRITE(numout,*) 
     288        ENDIF 
    240289     ENDIF 
    241290 
    … …  
    247296      ! --------------------------------------------- 
    248297      REWIND( numnam )                     
     298                
     299      nblendta(:,:) = 0 
     300      nbdysege = 0 
     301      nbdysegw = 0 
     302      nbdysegn = 0 
     303      nbdysegs = 0 
     304      icount   = 0 ! count user defined segments 
     305      ! Dimensions below are used to allocate arrays to read external data 
     306      jpbdtas = 1 ! Maximum size of boundary data (structured case) 
     307      jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 
     308 
    249309      DO ib_bdy = 1, nb_bdy 
    250310 
    251          jpbdta = 1 
    252311         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    253312  
     313            icount = icount + 1 
    254314            ! No REWIND here because may need to read more than one nambdy_index namelist. 
    255315            READ  ( numnam, nambdy_index ) 
    256316 
    257             ! Automatic boundary definition: if nbdysegX = -1 
    258             ! set boundary to whole side of model domain. 
    259             IF( nbdysege == -1 ) THEN 
    260                nbdysege = 1 
    261                jpieob(1) = jpiglo - 1 
    262                jpjedt(1) = 2 
    263                jpjeft(1) = jpjglo - 1 
    264             ENDIF 
    265             IF( nbdysegw == -1 ) THEN 
    266                nbdysegw = 1 
    267                jpiwob(1) = 2 
    268                jpjwdt(1) = 2 
    269                jpjwft(1) = jpjglo - 1 
    270             ENDIF 
    271             IF( nbdysegn == -1 ) THEN 
    272                nbdysegn = 1 
    273                jpjnob(1) = jpjglo - 1 
    274                jpindt(1) = 2 
    275                jpinft(1) = jpiglo - 1 
    276             ENDIF 
    277             IF( nbdysegs == -1 ) THEN 
    278                nbdysegs = 1 
    279                jpjsob(1) = 2 
    280                jpisdt(1) = 2 
    281                jpisft(1) = jpiglo - 1 
    282             ENDIF 
    283  
    284             nblendta(:,ib_bdy) = 0 
    285             DO iseg = 1, nbdysege 
    286                igrd = 1 
    287                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
    288                igrd = 2 
    289                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
    290                igrd = 3 
    291                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)                
    292             ENDDO 
    293             DO iseg = 1, nbdysegw 
    294                igrd = 1 
    295                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
    296                igrd = 2 
    297                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
    298                igrd = 3 
    299                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)                
    300             ENDDO 
    301             DO iseg = 1, nbdysegn 
    302                igrd = 1 
    303                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1                
    304                igrd = 2 
    305                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)                
    306                igrd = 3 
    307                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 
    308             ENDDO 
    309             DO iseg = 1, nbdysegs 
    310                igrd = 1 
    311                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
    312                igrd = 2 
    313                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 
    314                igrd = 3 
    315                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
    316             ENDDO 
    317  
    318             nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
    319             jpbdta = MAXVAL(nblendta(:,ib_bdy))                
    320  
     317            SELECT CASE ( TRIM(ctypebdy) ) 
     318              CASE( 'N' ) 
     319                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     320                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
     321                    nbdybeg  = 2 
     322                    nbdyend  = jpiglo - 1 
     323                 ENDIF 
     324                 nbdysegn = nbdysegn + 1 
     325                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number 
     326                 jpjnob(nbdysegn) = nbdyind 
     327                 jpindt(nbdysegn) = nbdybeg 
     328                 jpinft(nbdysegn) = nbdyend 
     329                 ! 
     330              CASE( 'S' ) 
     331                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     332                    nbdyind  = 2           ! set boundary to whole side of model domain. 
     333                    nbdybeg  = 2 
     334                    nbdyend  = jpiglo - 1 
     335                 ENDIF 
     336                 nbdysegs = nbdysegs + 1 
     337                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number 
     338                 jpjsob(nbdysegs) = nbdyind 
     339                 jpisdt(nbdysegs) = nbdybeg 
     340                 jpisft(nbdysegs) = nbdyend 
     341                 ! 
     342              CASE( 'E' ) 
     343                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     344                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
     345                    nbdybeg  = 2 
     346                    nbdyend  = jpjglo - 1 
     347                 ENDIF 
     348                 nbdysege = nbdysege + 1  
     349                 npckge(nbdysege) = ib_bdy ! Save bdy package number 
     350                 jpieob(nbdysege) = nbdyind 
     351                 jpjedt(nbdysege) = nbdybeg 
     352                 jpjeft(nbdysege) = nbdyend 
     353                 ! 
     354              CASE( 'W' ) 
     355                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     356                    nbdyind  = 2           ! set boundary to whole side of model domain. 
     357                    nbdybeg  = 2 
     358                    nbdyend  = jpjglo - 1 
     359                 ENDIF 
     360                 nbdysegw = nbdysegw + 1 
     361                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number 
     362                 jpiwob(nbdysegw) = nbdyind 
     363                 jpjwdt(nbdysegw) = nbdybeg 
     364                 jpjwft(nbdysegw) = nbdyend 
     365                 ! 
     366              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
     367            END SELECT 
     368 
     369            ! For simplicity we assume that in case of straight bdy, arrays have the same length 
     370            ! (even if it is true that last tangential velocity points 
     371            ! are useless). This simplifies a little bit boundary data format (and agrees with format 
     372            ! used so far in obc package) 
     373 
     374            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 
     375            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 
     376            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 
     377            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 
    321378 
    322379         ELSE            ! Read size of arrays in boundary coordinates file. 
    323  
    324  
    325380            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    326             jpbdta = 1 
    327381            DO igrd = 1, jpbgrd 
    328382               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    329383               nblendta(igrd,ib_bdy) = kdimsz(1) 
    330                jpbdta = MAX(jpbdta, kdimsz(1)) 
    331             ENDDO 
     384               jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     385            ENDDO 
     386            CALL iom_close( inum ) 
    332387 
    333388         ENDIF  
    … …  
    335390      ENDDO ! ib_bdy 
    336391 
    337       ! Allocate arrays 
    338       !--------------- 
    339       ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    340          &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    341  
    342       ALLOCATE( dta_global(jpbdta, 1, jpk) ) 
     392      IF (nb_bdy>0) THEN 
     393         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
     394 
     395         ! Allocate arrays 
     396         !--------------- 
     397         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
     398            &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     399 
     400         ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
     401         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
     402         !  
     403      ENDIF 
     404 
     405      ! Now look for crossings in user (namelist) defined open boundary segments: 
     406      !-------------------------------------------------------------------------- 
     407      IF ( icount>0 ) CALL bdy_ctl_seg 
    343408 
    344409      ! Calculate global boundary index arrays or read in from file 
    345       !------------------------------------------------------------ 
    346       REWIND( numnam )                     
     410      !------------------------------------------------------------                
     411      ! 1. Read global index arrays from boundary coordinates file. 
    347412      DO ib_bdy = 1, nb_bdy 
    348413 
    349          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 
    350  
    351             ! No REWIND here because may need to read more than one nambdy_index namelist. 
    352             READ  ( numnam, nambdy_index ) 
    353  
    354             ! Automatic boundary definition: if nbdysegX = -1 
    355             ! set boundary to whole side of model domain. 
    356             IF( nbdysege == -1 ) THEN 
    357                nbdysege = 1 
    358                jpieob(1) = jpiglo - 1 
    359                jpjedt(1) = 2 
    360                jpjeft(1) = jpjglo - 1 
    361             ENDIF 
    362             IF( nbdysegw == -1 ) THEN 
    363                nbdysegw = 1 
    364                jpiwob(1) = 2 
    365                jpjwdt(1) = 2 
    366                jpjwft(1) = jpjglo - 1 
    367             ENDIF 
    368             IF( nbdysegn == -1 ) THEN 
    369                nbdysegn = 1 
    370                jpjnob(1) = jpjglo - 1 
    371                jpindt(1) = 2 
    372                jpinft(1) = jpiglo - 1 
    373             ENDIF 
    374             IF( nbdysegs == -1 ) THEN 
    375                nbdysegs = 1 
    376                jpjsob(1) = 2 
    377                jpisdt(1) = 2 
    378                jpisft(1) = jpiglo - 1 
    379             ENDIF 
    380  
    381             ! ------------ T points ------------- 
    382             igrd = 1   
    383             icount = 0 
    384             DO ir = 1, nn_rimwidth(ib_bdy) 
    385                ! east 
    386                DO iseg = 1, nbdysege 
    387                   DO ij = jpjedt(iseg), jpjeft(iseg) 
    388                      icount = icount + 1 
    389                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
    390                      nbjdta(icount, igrd, ib_bdy) = ij 
    391                      nbrdta(icount, igrd, ib_bdy) = ir 
    392                   ENDDO 
    393                ENDDO 
    394                ! west 
    395                DO iseg = 1, nbdysegw 
    396                   DO ij = jpjwdt(iseg), jpjwft(iseg) 
    397                      icount = icount + 1 
    398                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    399                      nbjdta(icount, igrd, ib_bdy) = ij 
    400                      nbrdta(icount, igrd, ib_bdy) = ir 
    401                   ENDDO 
    402                ENDDO 
    403                ! north 
    404                DO iseg = 1, nbdysegn 
    405                   DO ii = jpindt(iseg), jpinft(iseg) 
    406                      icount = icount + 1 
    407                      nbidta(icount, igrd, ib_bdy) = ii 
    408                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
    409                      nbrdta(icount, igrd, ib_bdy) = ir 
    410                   ENDDO 
    411                ENDDO 
    412                ! south 
    413                DO iseg = 1, nbdysegs 
    414                   DO ii = jpisdt(iseg), jpisft(iseg) 
    415                      icount = icount + 1 
    416                      nbidta(icount, igrd, ib_bdy) = ii 
    417                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    418                      nbrdta(icount, igrd, ib_bdy) = ir 
    419                   ENDDO 
    420                ENDDO 
    421             ENDDO 
    422  
    423             ! ------------ U points ------------- 
    424             igrd = 2   
    425             icount = 0 
    426             DO ir = 1, nn_rimwidth(ib_bdy) 
    427                ! east 
    428                DO iseg = 1, nbdysege 
    429                   DO ij = jpjedt(iseg), jpjeft(iseg) 
    430                      icount = icount + 1 
    431                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 
    432                      nbjdta(icount, igrd, ib_bdy) = ij 
    433                      nbrdta(icount, igrd, ib_bdy) = ir 
    434                   ENDDO 
    435                ENDDO 
    436                ! west 
    437                DO iseg = 1, nbdysegw 
    438                   DO ij = jpjwdt(iseg), jpjwft(iseg) 
    439                      icount = icount + 1 
    440                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    441                      nbjdta(icount, igrd, ib_bdy) = ij 
    442                      nbrdta(icount, igrd, ib_bdy) = ir 
    443                   ENDDO 
    444                ENDDO 
    445                ! north 
    446                DO iseg = 1, nbdysegn 
    447                   DO ii = jpindt(iseg), jpinft(iseg) - 1 
    448                      icount = icount + 1 
    449                      nbidta(icount, igrd, ib_bdy) = ii 
    450                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
    451                      nbrdta(icount, igrd, ib_bdy) = ir 
    452                   ENDDO 
    453                ENDDO 
    454                ! south 
    455                DO iseg = 1, nbdysegs 
    456                   DO ii = jpisdt(iseg), jpisft(iseg) - 1 
    457                      icount = icount + 1 
    458                      nbidta(icount, igrd, ib_bdy) = ii 
    459                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    460                      nbrdta(icount, igrd, ib_bdy) = ir 
    461                   ENDDO 
    462                ENDDO 
    463             ENDDO 
    464  
    465             ! ------------ V points ------------- 
    466             igrd = 3   
    467             icount = 0 
    468             DO ir = 1, nn_rimwidth(ib_bdy) 
    469                ! east 
    470                DO iseg = 1, nbdysege 
    471                   DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
    472                      icount = icount + 1 
    473                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
    474                      nbjdta(icount, igrd, ib_bdy) = ij 
    475                      nbrdta(icount, igrd, ib_bdy) = ir 
    476                   ENDDO 
    477                ENDDO 
    478                ! west 
    479                DO iseg = 1, nbdysegw 
    480                   DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
    481                      icount = icount + 1 
    482                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    483                      nbjdta(icount, igrd, ib_bdy) = ij 
    484                      nbrdta(icount, igrd, ib_bdy) = ir 
    485                   ENDDO 
    486                ENDDO 
    487                ! north 
    488                DO iseg = 1, nbdysegn 
    489                   DO ii = jpindt(iseg), jpinft(iseg) 
    490                      icount = icount + 1 
    491                      nbidta(icount, igrd, ib_bdy) = ii 
    492                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 
    493                      nbrdta(icount, igrd, ib_bdy) = ir 
    494                   ENDDO 
    495                ENDDO 
    496                ! south 
    497                DO iseg = 1, nbdysegs 
    498                   DO ii = jpisdt(iseg), jpisft(iseg) 
    499                      icount = icount + 1 
    500                      nbidta(icount, igrd, ib_bdy) = ii 
    501                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    502                      nbrdta(icount, igrd, ib_bdy) = ir 
    503                   ENDDO 
    504                ENDDO 
    505             ENDDO 
    506  
    507          ELSE            ! Read global index arrays from boundary coordinates file. 
    508  
     414         IF( ln_coords_file(ib_bdy) ) THEN 
     415 
     416            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    509417            DO igrd = 1, jpbgrd 
    510418               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
    … …  
    527435               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
    528436                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
    529  
    530437            END DO 
    531438            CALL iom_close( inum ) 
    … …  
    533440         ENDIF  
    534441 
    535       ENDDO  
     442      ENDDO       
     443     
     444      ! 2. Now fill indices corresponding to straight open boundary arrays: 
     445      ! East 
     446      !----- 
     447      DO iseg = 1, nbdysege 
     448         ib_bdy = npckge(iseg) 
     449         ! 
     450         ! ------------ T points ------------- 
     451         igrd=1 
     452         icount=0 
     453         DO ir = 1, nn_rimwidth(ib_bdy) 
     454            DO ij = jpjedt(iseg), jpjeft(iseg) 
     455               icount = icount + 1 
     456               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     457               nbjdta(icount, igrd, ib_bdy) = ij 
     458               nbrdta(icount, igrd, ib_bdy) = ir 
     459            ENDDO 
     460         ENDDO 
     461         ! 
     462         ! ------------ U points ------------- 
     463         igrd=2 
     464         icount=0 
     465         DO ir = 1, nn_rimwidth(ib_bdy) 
     466            DO ij = jpjedt(iseg), jpjeft(iseg) 
     467               icount = icount + 1 
     468               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
     469               nbjdta(icount, igrd, ib_bdy) = ij 
     470               nbrdta(icount, igrd, ib_bdy) = ir 
     471            ENDDO 
     472         ENDDO 
     473         ! 
     474         ! ------------ V points ------------- 
     475         igrd=3 
     476         icount=0 
     477         DO ir = 1, nn_rimwidth(ib_bdy) 
     478!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     479            DO ij = jpjedt(iseg), jpjeft(iseg) 
     480               icount = icount + 1 
     481               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     482               nbjdta(icount, igrd, ib_bdy) = ij 
     483               nbrdta(icount, igrd, ib_bdy) = ir 
     484            ENDDO 
     485            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     486            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     487         ENDDO 
     488      ENDDO 
     489      ! 
     490      ! West 
     491      !----- 
     492      DO iseg = 1, nbdysegw 
     493         ib_bdy = npckgw(iseg) 
     494         ! 
     495         ! ------------ T points ------------- 
     496         igrd=1 
     497         icount=0 
     498         DO ir = 1, nn_rimwidth(ib_bdy) 
     499            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     500               icount = icount + 1 
     501               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     502               nbjdta(icount, igrd, ib_bdy) = ij 
     503               nbrdta(icount, igrd, ib_bdy) = ir 
     504            ENDDO 
     505         ENDDO 
     506         ! 
     507         ! ------------ U points ------------- 
     508         igrd=2 
     509         icount=0 
     510         DO ir = 1, nn_rimwidth(ib_bdy) 
     511            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     512               icount = icount + 1 
     513               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     514               nbjdta(icount, igrd, ib_bdy) = ij 
     515               nbrdta(icount, igrd, ib_bdy) = ir 
     516            ENDDO 
     517         ENDDO 
     518         ! 
     519         ! ------------ V points ------------- 
     520         igrd=3 
     521         icount=0 
     522         DO ir = 1, nn_rimwidth(ib_bdy) 
     523!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     524            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     525               icount = icount + 1 
     526               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     527               nbjdta(icount, igrd, ib_bdy) = ij 
     528               nbrdta(icount, igrd, ib_bdy) = ir 
     529            ENDDO 
     530            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     531            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     532         ENDDO 
     533      ENDDO 
     534      ! 
     535      ! North 
     536      !----- 
     537      DO iseg = 1, nbdysegn 
     538         ib_bdy = npckgn(iseg) 
     539         ! 
     540         ! ------------ T points ------------- 
     541         igrd=1 
     542         icount=0 
     543         DO ir = 1, nn_rimwidth(ib_bdy) 
     544            DO ii = jpindt(iseg), jpinft(iseg) 
     545               icount = icount + 1 
     546               nbidta(icount, igrd, ib_bdy) = ii 
     547               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
     548               nbrdta(icount, igrd, ib_bdy) = ir 
     549            ENDDO 
     550         ENDDO 
     551         ! 
     552         ! ------------ U points ------------- 
     553         igrd=2 
     554         icount=0 
     555         DO ir = 1, nn_rimwidth(ib_bdy) 
     556!            DO ii = jpindt(iseg), jpinft(iseg) - 1 
     557            DO ii = jpindt(iseg), jpinft(iseg) 
     558               icount = icount + 1 
     559               nbidta(icount, igrd, ib_bdy) = ii 
     560               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
     561               nbrdta(icount, igrd, ib_bdy) = ir 
     562            ENDDO 
     563            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     564            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     565         ENDDO 
     566         ! 
     567         ! ------------ V points ------------- 
     568         igrd=3 
     569         icount=0 
     570         DO ir = 1, nn_rimwidth(ib_bdy) 
     571            DO ii = jpindt(iseg), jpinft(iseg) 
     572               icount = icount + 1 
     573               nbidta(icount, igrd, ib_bdy) = ii 
     574               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
     575               nbrdta(icount, igrd, ib_bdy) = ir 
     576            ENDDO 
     577         ENDDO 
     578      ENDDO 
     579      ! 
     580      ! South 
     581      !----- 
     582      DO iseg = 1, nbdysegs 
     583         ib_bdy = npckgs(iseg) 
     584         ! 
     585         ! ------------ T points ------------- 
     586         igrd=1 
     587         icount=0 
     588         DO ir = 1, nn_rimwidth(ib_bdy) 
     589            DO ii = jpisdt(iseg), jpisft(iseg) 
     590               icount = icount + 1 
     591               nbidta(icount, igrd, ib_bdy) = ii 
     592               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     593               nbrdta(icount, igrd, ib_bdy) = ir 
     594            ENDDO 
     595         ENDDO 
     596         ! 
     597         ! ------------ U points ------------- 
     598         igrd=2 
     599         icount=0 
     600         DO ir = 1, nn_rimwidth(ib_bdy) 
     601!            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     602            DO ii = jpisdt(iseg), jpisft(iseg) 
     603               icount = icount + 1 
     604               nbidta(icount, igrd, ib_bdy) = ii 
     605               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     606               nbrdta(icount, igrd, ib_bdy) = ir 
     607            ENDDO 
     608            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     609            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     610         ENDDO 
     611         ! 
     612         ! ------------ V points ------------- 
     613         igrd=3 
     614         icount=0 
     615         DO ir = 1, nn_rimwidth(ib_bdy) 
     616            DO ii = jpisdt(iseg), jpisft(iseg) 
     617               icount = icount + 1 
     618               nbidta(icount, igrd, ib_bdy) = ii 
     619               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     620               nbrdta(icount, igrd, ib_bdy) = ir 
     621            ENDDO 
     622         ENDDO 
     623      ENDDO 
     624 
     625      !  Deal with duplicated points 
     626      !----------------------------- 
     627      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated) 
     628      ! if their distance to the bdy is greater than the other 
     629      ! If their distance are the same, just keep only one to avoid updating a point twice 
     630      DO igrd = 1, jpbgrd 
     631         DO ib_bdy1 = 1, nb_bdy 
     632            DO ib_bdy2 = 1, nb_bdy 
     633               IF (ib_bdy1/=ib_bdy2) THEN 
     634                  DO ib1 = 1, nblendta(igrd,ib_bdy1) 
     635                     DO ib2 = 1, nblendta(igrd,ib_bdy2) 
     636                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 
     637                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
     638!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
     639!                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
     640!                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
     641                           ! keep only points with the lowest distance to boundary: 
     642                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 
     643                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     644                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     645                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 
     646                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     647                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     648                           ! Arbitrary choice if distances are the same: 
     649                           ELSE 
     650                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     651                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     652                           ENDIF 
     653                        END IF 
     654                     END DO 
     655                  END DO 
     656               ENDIF 
     657            END DO 
     658         END DO 
     659      END DO 
    536660 
    537661      ! Work out dimensions of boundary data on each processor 
    538662      ! ------------------------------------------------------ 
    539       
    540       iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
    541       ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
    542       is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
    543       in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
     663 
     664      ! Rather assume that boundary data indices are given on global domain 
     665      ! TO BE DISCUSSED ? 
     666!      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
     667!      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
     668!      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
     669!      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1       
     670      iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2 
     671      ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1 
     672      is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2 
     673      in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1 
    544674 
    545675      DO ib_bdy = 1, nb_bdy 
    … …  
    577707         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
    578708         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
     709         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 
    579710         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
    580711         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
    … …  
    596727                     ! 
    597728                     icount = icount  + 1 
    598                      idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
    599                      idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
     729 
     730                     ! Rather assume that boundary data indices are given on global domain 
     731                     ! TO BE DISCUSSED ? 
     732!                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
     733!                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
     734                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom 
     735                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom 
    600736                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
    601737                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
    … …  
    611747               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    612748               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
    613 !              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic 
    614 !              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear 
     749!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
     750!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear 
     751            END DO 
     752         END DO  
     753 
     754         ! Compute damping coefficients 
     755         ! ---------------------------- 
     756         DO igrd = 1, jpbgrd 
     757            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     758               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
     759               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &  
     760               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    615761            END DO 
    616762         END DO  
    … …  
    627773      !          = 0  elsewhere    
    628774  
    629       IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution 
    630          zmask(         :                ,:) = 0.e0 
    631          zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0           
    632       ELSE IF( ln_mask_file ) THEN 
     775      IF( ln_mask_file ) THEN 
    633776         CALL iom_open( cn_mask_file, inum ) 
    634          CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 
     777         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 
    635778         CALL iom_close( inum ) 
    636       ELSE 
    637          zmask(:,:) = 1.e0 
    638       ENDIF 
    639  
    640       DO ij = 1, nlcj      ! Save mask over local domain       
    641          DO ii = 1, nlci 
    642             bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 
     779 
     780         ! Derive mask on U and V grid from mask on T grid 
     781         bdyumask(:,:) = 0.e0 
     782         bdyvmask(:,:) = 0.e0 
     783         DO ij=1, jpjm1 
     784            DO ii=1, jpim1 
     785               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
     786               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
     787            END DO 
    643788         END DO 
    644       END DO 
    645  
    646       ! Derive mask on U and V grid from mask on T grid 
    647       bdyumask(:,:) = 0.e0 
    648       bdyvmask(:,:) = 0.e0 
    649       DO ij=1, jpjm1 
    650          DO ii=1, jpim1 
    651             bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
    652             bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
     789         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
     790 
     791 
     792         ! Mask corrections 
     793         ! ---------------- 
     794         DO ik = 1, jpkm1 
     795            DO ij = 1, jpj 
     796               DO ii = 1, jpi 
     797                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
     798                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
     799                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
     800                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij) 
     801               END DO       
     802            END DO 
    653803         END DO 
    654       END DO 
    655       CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
    656  
    657  
    658       ! Mask corrections 
    659       ! ---------------- 
    660       DO ik = 1, jpkm1 
    661          DO ij = 1, jpj 
    662             DO ii = 1, jpi 
    663                tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
    664                umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
    665                vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
    666                bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij) 
    667             END DO       
     804 
     805         DO ik = 1, jpkm1 
     806            DO ij = 2, jpjm1 
     807               DO ii = 2, jpim1 
     808                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
     809                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
     810               END DO       
     811            END DO 
    668812         END DO 
    669       END DO 
    670  
    671       DO ik = 1, jpkm1 
    672          DO ij = 2, jpjm1 
    673             DO ii = 2, jpim1 
    674                fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
    675                   &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
    676             END DO       
    677          END DO 
    678       END DO 
    679  
    680       tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)              
     813 
     814         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 
     815 
     816      ENDIF ! ln_mask_file=.TRUE. 
     817       
    681818      bdytmask(:,:) = tmask(:,:,1) 
    682819 
    … …  
    753890            END IF 
    754891         END DO 
    755   
     892 
    756893         IF( icount /= 0 ) THEN 
    757894            IF(lwp) WRITE(numout,*) 
    … …  
    767904      ! Compute total lateral surface for volume correction: 
    768905      ! ---------------------------------------------------- 
     906      ! JC: this must be done at each time step with key_vvl 
    769907      bdysurftot = 0.e0  
    770908      IF( ln_vol ) THEN   
    … …  
    800938      ! Tidy up 
    801939      !-------- 
    802       DEALLOCATE(nbidta, nbjdta, nbrdta) 
     940      IF (nb_bdy>0) THEN 
     941         DEALLOCATE(nbidta, nbjdta, nbrdta) 
     942      ENDIF 
    803943 
    804944      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
    805945 
    806946   END SUBROUTINE bdy_init 
     947 
     948   SUBROUTINE bdy_ctl_seg 
     949      !!---------------------------------------------------------------------- 
     950      !!                 ***  ROUTINE bdy_ctl_seg  *** 
     951      !! 
     952      !! ** Purpose :   Check straight open boundary segments location 
     953      !! 
     954      !! ** Method  :   - Look for open boundary corners 
     955      !!                - Check that segments start or end on land  
     956      !!---------------------------------------------------------------------- 
     957      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest   
     958      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns   
     959      REAL(wp), DIMENSION(2) ::   ztestmask 
     960      !!---------------------------------------------------------------------- 
     961      ! 
     962      IF (lwp) WRITE(numout,*) ' ' 
     963      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments' 
     964      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     965      ! 
     966      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege 
     967      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw 
     968      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn 
     969      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs 
     970      ! 1. Check bounds 
     971      !---------------- 
     972      DO ib = 1, nbdysegn 
     973         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib) 
     974         IF ((jpjnob(ib).ge.jpjglo-1).or.&  
     975            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     976         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     977         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     978         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
     979      END DO 
     980      ! 
     981      DO ib = 1, nbdysegs 
     982         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib) 
     983         IF ((jpjsob(ib).ge.jpjglo-1).or.&  
     984            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     985         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     986         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     987         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
     988      END DO 
     989      ! 
     990      DO ib = 1, nbdysege 
     991         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib) 
     992         IF ((jpieob(ib).ge.jpiglo-1).or.&  
     993            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     994         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     995         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     996         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
     997      END DO 
     998      ! 
     999      DO ib = 1, nbdysegw 
     1000         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib) 
     1001         IF ((jpiwob(ib).ge.jpiglo-1).or.&  
     1002            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     1003         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     1004         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     1005         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
     1006      ENDDO 
     1007      ! 
     1008      !       
     1009      ! 2. Look for segment crossings 
     1010      !------------------------------  
     1011      IF (lwp) WRITE(numout,*) '**Look for segments corners  :' 
     1012      ! 
     1013      itest = 0 ! corner number 
     1014      ! 
     1015      ! flag to detect if start or end of open boundary belongs to a corner 
     1016      ! if not (=0), it must be on land. 
     1017      ! if a corner is detected, save bdy package number for further tests 
     1018      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0. 
     1019      ! South/West crossings 
     1020      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN 
     1021         DO ib1 = 1, nbdysegw         
     1022            DO ib2 = 1, nbdysegs 
     1023               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. & 
     1024                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. & 
     1025                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. & 
     1026                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN 
     1027                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN  
     1028                     ! We have a possible South-West corner                       
     1029!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)  
     1030!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2) 
     1031                     icornw(ib1,1) = npckgs(ib2) 
     1032                     icorns(ib2,1) = npckgw(ib1) 
     1033                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 
     1034                     IF(lwp) WRITE(numout,*) 
     1035                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1036                     &                                     jpisft(ib2), jpjwft(ib1) 
     1037                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1038                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &  
     1039                     &                                                    ' and South segment: ',npckgs(ib2) 
     1040                     IF(lwp) WRITE(numout,*) 
     1041                     nstop = nstop + 1 
     1042                  ELSE 
     1043                     IF(lwp) WRITE(numout,*) 
     1044                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices' 
     1045                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , & 
     1046                     &                                                    ' and South segment: ',npckgs(ib2) 
     1047                     IF(lwp) WRITE(numout,*) 
     1048                     nstop = nstop+1 
     1049                  END IF 
     1050               END IF 
     1051            END DO 
     1052         END DO 
     1053      END IF 
     1054      ! 
     1055      ! South/East crossings 
     1056      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN 
     1057         DO ib1 = 1, nbdysege 
     1058            DO ib2 = 1, nbdysegs 
     1059               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. & 
     1060                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. & 
     1061                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. & 
     1062                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN 
     1063                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN 
     1064                     ! We have a possible South-East corner  
     1065!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)  
     1066!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2) 
     1067                     icorne(ib1,1) = npckgs(ib2) 
     1068                     icorns(ib2,2) = npckge(ib1) 
     1069                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 
     1070                     IF(lwp) WRITE(numout,*) 
     1071                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1072                     &                                     jpisdt(ib2), jpjeft(ib1) 
     1073                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1074                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), & 
     1075                     &                                                    ' and South segment: ',npckgs(ib2) 
     1076                     IF(lwp) WRITE(numout,*) 
     1077                     nstop = nstop + 1 
     1078                  ELSE 
     1079                     IF(lwp) WRITE(numout,*) 
     1080                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices' 
     1081                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
     1082                     &                                                    ' and South segment: ',npckgs(ib2) 
     1083                     IF(lwp) WRITE(numout,*) 
     1084                     nstop = nstop + 1 
     1085                  END IF 
     1086               END IF 
     1087            END DO 
     1088         END DO 
     1089      END IF 
     1090      ! 
     1091      ! North/West crossings 
     1092      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN 
     1093         DO ib1 = 1, nbdysegw         
     1094            DO ib2 = 1, nbdysegn 
     1095               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. & 
     1096                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. & 
     1097                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. & 
     1098                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN 
     1099                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN 
     1100                     ! We have a possible North-West corner  
     1101!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)  
     1102!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2) 
     1103                     icornw(ib1,2) = npckgn(ib2) 
     1104                     icornn(ib2,1) = npckgw(ib1) 
     1105                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 
     1106                     IF(lwp) WRITE(numout,*) 
     1107                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1108                     &                                     jpinft(ib2), jpjwdt(ib1) 
     1109                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1110                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
     1111                     &                                                    ' and North segment: ',npckgn(ib2) 
     1112                     IF(lwp) WRITE(numout,*) 
     1113                     nstop = nstop + 1 
     1114                  ELSE 
     1115                     IF(lwp) WRITE(numout,*) 
     1116                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices' 
     1117                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), & 
     1118                     &                                                    ' and North segment: ',npckgn(ib2) 
     1119                     IF(lwp) WRITE(numout,*) 
     1120                     nstop = nstop + 1 
     1121                  END IF 
     1122               END IF 
     1123            END DO 
     1124         END DO 
     1125      END IF 
     1126      ! 
     1127      ! North/East crossings 
     1128      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN 
     1129         DO ib1 = 1, nbdysege         
     1130            DO ib2 = 1, nbdysegn 
     1131               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. & 
     1132                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. & 
     1133                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. & 
     1134                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN 
     1135                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN 
     1136                     ! We have a possible North-East corner  
     1137!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1) 
     1138!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2) 
     1139                     icorne(ib1,2) = npckgn(ib2) 
     1140                     icornn(ib2,2) = npckge(ib1) 
     1141                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 
     1142                     IF(lwp) WRITE(numout,*) 
     1143                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1144                     &                                     jpindt(ib2), jpjedt(ib1) 
     1145                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1146                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), & 
     1147                     &                                                    ' and North segment: ',npckgn(ib2) 
     1148                     IF(lwp) WRITE(numout,*) 
     1149                     nstop = nstop + 1 
     1150                  ELSE 
     1151                     IF(lwp) WRITE(numout,*) 
     1152                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices' 
     1153                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
     1154                     &                                                    ' and North segment: ',npckgn(ib2) 
     1155                     IF(lwp) WRITE(numout,*) 
     1156                     nstop = nstop + 1 
     1157                  END IF 
     1158               END IF 
     1159            END DO 
     1160         END DO 
     1161      END IF 
     1162      ! 
     1163      ! 3. Check if segment extremities are on land 
     1164      !--------------------------------------------  
     1165      ! 
     1166      ! West segments 
     1167      DO ib = 1, nbdysegw 
     1168         ! get mask at boundary extremities: 
     1169         ztestmask(1:2)=0. 
     1170         DO ji = 1, jpi 
     1171            DO jj = 1, jpj              
     1172              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &  
     1173               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1174              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &  
     1175               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1176            END DO 
     1177         END DO 
     1178         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1179 
     1180         IF (ztestmask(1)==1) THEN  
     1181            IF (icornw(ib,1)==0) THEN 
     1182               IF(lwp) WRITE(numout,*) 
     1183               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
     1184               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
     1185               IF(lwp) WRITE(numout,*) 
     1186               nstop = nstop + 1 
     1187            ELSE 
     1188               ! This is a corner 
     1189               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib) 
     1190               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1)) 
     1191               itest=itest+1 
     1192            ENDIF 
     1193         ENDIF 
     1194         IF (ztestmask(2)==1) THEN 
     1195            IF (icornw(ib,2)==0) THEN 
     1196               IF(lwp) WRITE(numout,*) 
     1197               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
     1198               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
     1199               IF(lwp) WRITE(numout,*) 
     1200               nstop = nstop + 1 
     1201            ELSE 
     1202               ! This is a corner 
     1203               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib) 
     1204               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2)) 
     1205               itest=itest+1 
     1206            ENDIF 
     1207         ENDIF 
     1208      END DO 
     1209      ! 
     1210      ! East segments 
     1211      DO ib = 1, nbdysege 
     1212         ! get mask at boundary extremities: 
     1213         ztestmask(1:2)=0. 
     1214         DO ji = 1, jpi 
     1215            DO jj = 1, jpj              
     1216              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &  
     1217               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1218              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &  
     1219               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1220            END DO 
     1221         END DO 
     1222         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1223 
     1224         IF (ztestmask(1)==1) THEN 
     1225            IF (icorne(ib,1)==0) THEN 
     1226               IF(lwp) WRITE(numout,*) 
     1227               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 
     1228               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
     1229               IF(lwp) WRITE(numout,*) 
     1230               nstop = nstop + 1  
     1231            ELSE 
     1232               ! This is a corner 
     1233               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib) 
     1234               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1)) 
     1235               itest=itest+1 
     1236            ENDIF 
     1237         ENDIF 
     1238         IF (ztestmask(2)==1) THEN 
     1239            IF (icorne(ib,2)==0) THEN 
     1240               IF(lwp) WRITE(numout,*) 
     1241               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 
     1242               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
     1243               IF(lwp) WRITE(numout,*) 
     1244               nstop = nstop + 1 
     1245            ELSE 
     1246               ! This is a corner 
     1247               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib) 
     1248               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2)) 
     1249               itest=itest+1 
     1250            ENDIF 
     1251         ENDIF 
     1252      END DO 
     1253      ! 
     1254      ! South segments 
     1255      DO ib = 1, nbdysegs 
     1256         ! get mask at boundary extremities: 
     1257         ztestmask(1:2)=0. 
     1258         DO ji = 1, jpi 
     1259            DO jj = 1, jpj              
     1260              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &  
     1261               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1262              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &  
     1263               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1264            END DO 
     1265         END DO 
     1266         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1267 
     1268         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 
     1269            IF(lwp) WRITE(numout,*) 
     1270            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
     1271            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
     1272            IF(lwp) WRITE(numout,*) 
     1273            nstop = nstop + 1 
     1274         ENDIF 
     1275         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 
     1276            IF(lwp) WRITE(numout,*) 
     1277            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
     1278            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
     1279            IF(lwp) WRITE(numout,*) 
     1280            nstop = nstop + 1 
     1281         ENDIF 
     1282      END DO 
     1283      ! 
     1284      ! North segments 
     1285      DO ib = 1, nbdysegn 
     1286         ! get mask at boundary extremities: 
     1287         ztestmask(1:2)=0. 
     1288         DO ji = 1, jpi 
     1289            DO jj = 1, jpj              
     1290              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &  
     1291               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1292              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &  
     1293               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1294            END DO 
     1295         END DO 
     1296         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1297 
     1298         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 
     1299            IF(lwp) WRITE(numout,*) 
     1300            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
     1301            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                   
     1302            IF(lwp) WRITE(numout,*) 
     1303            nstop = nstop + 1 
     1304         ENDIF 
     1305         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 
     1306            IF(lwp) WRITE(numout,*) 
     1307            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
     1308            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                   
     1309            IF(lwp) WRITE(numout,*) 
     1310            nstop = nstop + 1 
     1311         ENDIF 
     1312      END DO 
     1313      ! 
     1314      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found' 
     1315      ! 
     1316      ! Other tests TBD:  
     1317      ! segments completly on land 
     1318      ! optimized open boundary array length according to landmask 
     1319      ! Nudging layers that overlap with interior domain 
     1320      ! 
     1321   END SUBROUTINE bdy_ctl_seg 
     1322 
     1323   SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 
     1324      !!---------------------------------------------------------------------- 
     1325      !!                 ***  ROUTINE bdy_ctl_corn  *** 
     1326      !! 
     1327      !! ** Purpose :   Check numerical schemes consistency between 
     1328      !!                segments having a common corner 
     1329      !! 
     1330      !! ** Method  :    
     1331      !!---------------------------------------------------------------------- 
     1332      INTEGER, INTENT(in)  ::   ib1, ib2 
     1333      INTEGER :: itest 
     1334      !!---------------------------------------------------------------------- 
     1335      itest = 0 
     1336 
     1337      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1 
     1338      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1 
     1339      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1 
     1340      ! 
     1341      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1 
     1342      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1 
     1343      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1 
     1344      ! 
     1345      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1    
     1346      ! 
     1347      IF ( itest>0 ) THEN 
     1348         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2 
     1349         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                   
     1350         IF(lwp) WRITE(numout,*) 
     1351         nstop = nstop + 1 
     1352      ENDIF 
     1353      ! 
     1354   END SUBROUTINE bdy_ctl_corn 
    8071355 
    8081356#else 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r3294 r3651  
    88   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
    10    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     10   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
    1111   !!---------------------------------------------------------------------- 
    1212#if defined key_bdy 
    … …  
    1515   !!---------------------------------------------------------------------- 
    1616   !!   PUBLIC 
    17    !!      tide_init     : read of namelist and initialisation of tidal harmonics data 
     17   !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data 
    1818   !!      tide_update   : calculation of tidal forcing at each timestep 
    19    !!   PRIVATE 
    20    !!      uvset         :\ 
    21    !!      vday          : | Routines to correct tidal harmonics forcing for 
    22    !!      shpen         : | start time of integration 
    23    !!      ufset         : | 
    24    !!      vset          :/ 
    2519   !!---------------------------------------------------------------------- 
    2620   USE timing          ! Timing 
    … …  
    3428   USE bdy_oce         ! ocean open boundary conditions 
    3529   USE daymod          ! calendar 
     30   USE wrk_nemo        ! Memory allocation 
     31   USE tideini 
     32!   USE tide_mod       ! Useless ?? 
    3633   USE fldread, ONLY: fld_map 
    3734 
    … …  
    3936   PRIVATE 
    4037 
    41    PUBLIC   tide_init     ! routine called in nemo_init 
    42    PUBLIC   tide_update   ! routine called in bdydyn 
     38   PUBLIC   bdytide_init     ! routine called in bdy_init 
     39   PUBLIC   bdytide_update   ! routine called in bdy_dta 
    4340 
    4441   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
    45       INTEGER                                ::   ncpt       !: Actual number of tidal components 
    46       REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr) 
    47       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH 
    48       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U 
    49       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V 
     42      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file) 
     43      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file) 
     44      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file) 
     45      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.) 
     46      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.) 
     47      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.) 
    5048   END TYPE TIDES_DATA 
    5149 
    52    INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents 
    53  
    54    TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data 
    55     
     50   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
     51 
    5652   !!---------------------------------------------------------------------- 
    5753   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    6157CONTAINS 
    6258 
    63    SUBROUTINE tide_init 
    64       !!---------------------------------------------------------------------- 
    65       !!                    ***  SUBROUTINE tide_init  *** 
     59   SUBROUTINE bdytide_init 
     60      !!---------------------------------------------------------------------- 
     61      !!                    ***  SUBROUTINE bdytide_init  *** 
    6662      !!                      
    6763      !! ** Purpose : - Read in namelist for tides and initialise external 
    … …  
    7167      !! namelist variables 
    7268      !!------------------- 
    73       LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date 
    74       CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files 
    75       CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used. 
    76       REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
     69      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files 
     70      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data 
     71      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data 
    7772      !! 
    78       INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components 
    79       INTEGER                                   ::   ib_bdy, itide, ib  !: dummy loop indices 
     73      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices 
     74      INTEGER                                   ::   ii, ij              !: dummy loop indices 
    8075      INTEGER                                   ::   inum, igrd 
    81       INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays) 
    82       CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file  
    83       REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t           
    84       REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc 
    85       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data 
     76      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays) 
     77      INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts 
     78      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file  
     79      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data 
     80      REAL(wp), POINTER, DIMENSION(:,:)         ::   ztr, zti            !:  "     "    "   "   "   "        "      "  
    8681      !! 
    87       TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut    
     82      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
    8883      !! 
    89       NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed 
    90       !!---------------------------------------------------------------------- 
    91  
    92       IF( nn_timing == 1 ) CALL timing_start('tide_init') 
    93  
    94       IF(lwp) WRITE(numout,*) 
    95       IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 
    96       IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
     84      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
     85      !!---------------------------------------------------------------------- 
     86 
     87      IF( nn_timing == 1 ) CALL timing_start('bdytide_init') 
     88 
     89      IF (nb_bdy>0) THEN 
     90         IF(lwp) WRITE(numout,*) 
     91         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 
     92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     93      ENDIF 
     94 
     95      ln_bdytide_2ddta = .FALSE. 
     96      ln_bdytide_conj  = .FALSE. 
    9797 
    9898      REWIND(numnam) 
    … …  
    101101 
    102102            td => tides(ib_bdy) 
     103            nblen => idx_bdy(ib_bdy)%nblen 
     104            nblenrim => idx_bdy(ib_bdy)%nblenrim 
    103105 
    104106            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
    105             ln_tide_date = .false. 
    106107            filtide(:) = '' 
    107             tide_speed(:) = 0.0 
    108             tide_cpt(:) = '' 
    109108 
    110109            ! Don't REWIND here - may need to read more than one of these namelists. 
    111110            READ  ( numnam, nambdy_tide ) 
    112             !                                               ! Count number of components specified 
    113             td%ncpt = 0 
    114             DO itide = 1, jptides_max 
    115               IF( tide_cpt(itide) /= '' ) THEN 
    116                  td%ncpt = td%ncpt + 1 
    117               ENDIF 
    118             END DO 
    119  
    120             ! Fill in phase speeds from namelist 
    121             ALLOCATE( td%speed(td%ncpt) ) 
    122             td%speed = tide_speed(1:td%ncpt) 
    123  
    124             ! Find constituents in standard list 
    125             DO itide = 1, td%ncpt 
    126                nindx(itide) = 0 
    127                IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
    128                IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
    129                IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
    130                IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
    131                IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
    132                IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
    133                IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
    134                IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
    135                IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
    136                IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
    137                IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
    138                IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
    139                IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
    140                IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
    141                IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
    142                IF( nindx(itide) == 0  .AND. lwp ) THEN 
    143                   WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
    144                   CALL ctl_warn( ctmp1 ) 
    145                ENDIF 
    146             END DO 
    147111            !                                               ! Parameter control and print 
    148             IF( td%ncpt < 1 ) THEN  
    149                CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' ) 
     112            IF(lwp) WRITE(numout,*) '  ' 
     113            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
     114            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta 
     115            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj 
     116            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
     117            IF(lwp) THEN  
     118                    WRITE(numout,*) '             Tidal cpt name    -     Phase speed (deg/hr)'             
     119               DO itide = 1, nb_harmo 
     120                  WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)    
     121               END DO 
     122            ENDIF  
     123            IF(lwp) WRITE(numout,*) ' ' 
     124 
     125            ! Allocate space for tidal harmonics data - get size from OBC data arrays 
     126            ! ----------------------------------------------------------------------- 
     127 
     128            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
     129            ! relaxation area       
     130            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     131               ilen0(:)=nblen(:) 
    150132            ELSE 
    151                IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
    152                IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt 
    153                IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt) 
    154                IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
    155                IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt) 
     133               ilen0(:)=nblenrim(:) 
    156134            ENDIF 
    157135 
    158  
    159             ! Allocate space for tidal harmonics data -  
    160             ! get size from OBC data arrays 
    161             ! --------------------------------------- 
    162  
    163             ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh )  
    164             ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 
    165  
    166             ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )  
    167             ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 
    168  
    169             ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d )  
    170             ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 
    171  
    172             ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 
    173  
    174  
    175             ! Open files and read in tidal forcing data 
    176             ! ----------------------------------------- 
    177  
    178             DO itide = 1, td%ncpt 
    179                !                                                              ! SSH fields 
    180                clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
    181                IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    182                CALL iom_open( clfile, inum ) 
    183                CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    184                td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) 
    185                CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    186                td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     136            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 
     137            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 
     138 
     139            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 
     140            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 
     141 
     142            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 
     143            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 
     144 
     145            td%ssh0(:,:,:) = 0.e0 
     146            td%ssh(:,:,:) = 0.e0 
     147            td%u0(:,:,:) = 0.e0 
     148            td%u(:,:,:) = 0.e0 
     149            td%v0(:,:,:) = 0.e0 
     150            td%v(:,:,:) = 0.e0 
     151 
     152            IF (ln_bdytide_2ddta) THEN 
     153               ! It is assumed that each data file contains all complex harmonic amplitudes 
     154               ! given on the data domain (ie global, jpidta x jpjdta) 
     155               ! 
     156               CALL wrk_alloc( jpi, jpj, zti, ztr ) 
     157               ! 
     158               ! SSH fields 
     159               clfile = TRIM(filtide)//'_grid_T.nc' 
     160               CALL iom_open (clfile , inum )  
     161               igrd = 1                       ! Everything is at T-points here 
     162               DO itide = 1, nb_harmo 
     163                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 
     164                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )  
     165                  DO ib = 1, ilen0(igrd) 
     166                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     167                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     168                     td%ssh0(ib,itide,1) = ztr(ii,ij) 
     169                     td%ssh0(ib,itide,2) = zti(ii,ij) 
     170                  END DO 
     171               END DO  
    187172               CALL iom_close( inum ) 
    188                !                                                              ! U fields 
    189                clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
    190                IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    191                CALL iom_open( clfile, inum ) 
    192                CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    193                td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) 
    194                CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    195                td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     173               ! 
     174               ! U fields 
     175               clfile = TRIM(filtide)//'_grid_U.nc' 
     176               CALL iom_open (clfile , inum )  
     177               igrd = 2                       ! Everything is at U-points here 
     178               DO itide = 1, nb_harmo 
     179                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) ) 
     180                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) ) 
     181                  DO ib = 1, ilen0(igrd) 
     182                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     183                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     184                     td%u0(ib,itide,1) = ztr(ii,ij) 
     185                     td%u0(ib,itide,2) = zti(ii,ij) 
     186                  END DO 
     187               END DO 
    196188               CALL iom_close( inum ) 
    197                !                                                              ! V fields 
    198                clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
    199                IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    200                CALL iom_open( clfile, inum ) 
    201                CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    202                td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) 
    203                CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    204                td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     189               ! 
     190               ! V fields 
     191               clfile = TRIM(filtide)//'_grid_V.nc' 
     192               CALL iom_open (clfile , inum )  
     193               igrd = 3                       ! Everything is at V-points here 
     194               DO itide = 1, nb_harmo 
     195                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) ) 
     196                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) ) 
     197                  DO ib = 1, ilen0(igrd) 
     198                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     199                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     200                     td%v0(ib,itide,1) = ztr(ii,ij) 
     201                     td%v0(ib,itide,2) = zti(ii,ij) 
     202                  END DO 
     203               END DO   
    205204               CALL iom_close( inum ) 
    206205               ! 
    207             END DO ! end loop on tidal components 
    208  
    209             IF( ln_tide_date ) THEN      ! correct for date factors 
    210  
    211 !! used nmonth, nyear and nday from daymod.... 
    212                ! Calculate date corrects for 15 standard consituents 
    213                ! This is the initialisation step, so nday, nmonth, nyear are the  
    214                ! initial date/time of the integration. 
    215                  print *, nday,nmonth,nyear 
    216                  nyear  = int(ndate0 / 10000  )                          ! initial year 
    217                  nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month 
    218                  nday   = int(ndate0 - nyear * 10000 - nmonth * 100) 
    219  
    220                CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
    221  
    222                IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
    223  
    224                DO itide = 1, td%ncpt       ! loop on tidal components 
     206               CALL wrk_dealloc( jpi, jpj, ztr, zti )  
     207               ! 
     208            ELSE             
     209               ! 
     210               ! Read tidal data only on bdy segments 
     211               !  
     212               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 
     213 
     214               ! Open files and read in tidal forcing data 
     215               ! ----------------------------------------- 
     216 
     217               DO itide = 1, nb_harmo 
     218                  !                                                              ! SSH fields 
     219                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
     220                  CALL iom_open( clfile, inum ) 
     221                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     222                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 
     223                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     224                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     225                  CALL iom_close( inum ) 
     226                  !                                                              ! U fields 
     227                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
     228                  CALL iom_open( clfile, inum ) 
     229                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     230                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 
     231                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     232                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     233                  CALL iom_close( inum ) 
     234                  !                                                              ! V fields 
     235                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
     236                  CALL iom_open( clfile, inum ) 
     237                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     238                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 
     239                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     240                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     241                  CALL iom_close( inum ) 
    225242                  ! 
    226                   IF( nindx(itide) /= 0 ) THEN 
    227 !!gm use rpi  and rad global variable   
    228                      z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
    229                      z_atde=z_ftc(nindx(itide))*cos(z_arg) 
    230                      z_btde=z_ftc(nindx(itide))*sin(z_arg) 
    231                      IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   & 
    232                         &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
    233                   ELSE 
    234                      z_atde = 1.0_wp 
    235                      z_btde = 0.0_wp 
    236                   ENDIF 
    237                   !                                         !  elevation          
    238                   igrd = 1 
    239                   DO ib = 1, ilen0(igrd)                 
    240                      z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2) 
    241                      z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1) 
    242                      td%ssh(ib,itide,1) = z1t 
    243                      td%ssh(ib,itide,2) = z2t 
    244                   END DO 
    245                   !                                         !  u        
    246                   igrd = 2 
    247                   DO ib = 1, ilen0(igrd)                 
    248                      z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2) 
    249                      z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1) 
    250                      td%u(ib,itide,1) = z1t 
    251                      td%u(ib,itide,2) = z2t 
    252                   END DO 
    253                   !                                         !  v        
    254                   igrd = 3 
    255                   DO ib = 1, ilen0(igrd)                 
    256                      z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2) 
    257                      z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1) 
    258                      td%v(ib,itide,1) = z1t 
    259                      td%v(ib,itide,2) = z2t 
    260                   END DO 
    261                   ! 
    262                END DO   ! end loop on tidal components 
    263                ! 
    264             ENDIF ! date correction 
     243               END DO ! end loop on tidal components 
     244               ! 
     245               DEALLOCATE( dta_read ) 
     246            ENDIF ! ln_bdytide_2ddta=.true. 
     247            ! 
     248            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files 
     249               td%ssh0(:,:,2) = - td%ssh0(:,:,2) 
     250               td%u0  (:,:,2) = - td%u0  (:,:,2) 
     251               td%v0  (:,:,2) = - td%v0  (:,:,2) 
     252            ENDIF 
    265253            ! 
    266254         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
    … …  
    268256      END DO ! loop on ib_bdy 
    269257 
    270       IF( nn_timing == 1 ) CALL timing_stop('tide_init') 
    271  
    272    END SUBROUTINE tide_init 
    273  
    274  
    275    SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 
    276       !!---------------------------------------------------------------------- 
    277       !!                 ***  SUBROUTINE tide_update  *** 
     258      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init') 
     259 
     260   END SUBROUTINE bdytide_init 
     261 
     262   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset ) 
     263      !!---------------------------------------------------------------------- 
     264      !!                 ***  SUBROUTINE bdytide_update  *** 
    278265      !!                 
    279266      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
    280267      !!                 
    281268      !!---------------------------------------------------------------------- 
    282       INTEGER, INTENT( in )          ::   kt      ! Main timestep counter 
    283 !!gm doctor jit ==> kit 
    284       TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices 
    285       TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data 
    286       TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data 
    287       INTEGER,INTENT(in),OPTIONAL    ::   jit     ! Barotropic timestep counter (for timesplitting option) 
    288       INTEGER,INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
    289                                                        ! is present then units = subcycle timesteps. 
    290                                                        ! time_offset = 0 => get data at "now" time level 
    291                                                        ! time_offset = -1 => get data at "before" time level 
    292                                                        ! time_offset = +1 => get data at "after" time level 
    293                                                        ! etc. 
     269      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter 
     270      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices 
     271      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data 
     272      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data 
     273      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option) 
     274      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit 
     275                                                        ! is present then units = subcycle timesteps. 
     276                                                        ! time_offset = 0  => get data at "now"    time level 
     277                                                        ! time_offset = -1 => get data at "before" time level 
     278                                                        ! time_offset = +1 => get data at "after"  time level 
     279                                                        ! etc. 
    294280      !! 
    295       INTEGER                          :: itide, igrd, ib      ! dummy loop indices 
    296       INTEGER                          :: time_add             ! time offset in units of timesteps 
    297       REAL(wp)                         :: z_arg, z_sarg       
    298       REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
    299       !!---------------------------------------------------------------------- 
    300  
    301       IF( nn_timing == 1 ) CALL timing_start('tide_update') 
     281      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays) 
     282      INTEGER                          :: itide, igrd, ib   ! dummy loop indices 
     283      INTEGER                          :: time_add          ! time offset in units of timesteps 
     284      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp       
     285      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 
     286      !!---------------------------------------------------------------------- 
     287 
     288      IF( nn_timing == 1 ) CALL timing_start('bdytide_update') 
     289 
     290      ilen0(1) =  SIZE(td%ssh(:,1,1)) 
     291      ilen0(2) =  SIZE(td%u(:,1,1)) 
     292      ilen0(3) =  SIZE(td%v(:,1,1)) 
     293 
     294      zflag=1 
     295      IF ( PRESENT(jit) ) THEN 
     296        IF ( jit /= 1 ) zflag=0 
     297      ENDIF 
     298 
     299      IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN 
     300        ! 
     301        kt_tide = kt 
     302        ! 
     303        IF(lwp) THEN 
     304           WRITE(numout,*) 
     305           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 
     306           WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
     307        ENDIF 
     308        ! 
     309        CALL tide_init_elevation ( idx, td ) 
     310        CALL tide_init_velocities( idx, td ) 
     311        ! 
     312      ENDIF  
    302313 
    303314      time_add = 0 
    … …  
    306317      ENDIF 
    307318          
    308       ! Note tide phase speeds are in deg/hour, so we need to convert the 
    309       ! elapsed time in seconds to hours by dividing by 3600.0 
    310319      IF( PRESENT(jit) ) THEN   
    311          z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
     320         z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) 
    312321      ELSE                               
    313          z_arg = (kt+time_add) * rdt * rad / 3600.0 
     322         z_arg = ((kt-kt_tide)+time_add) * rdt 
    314323      ENDIF 
    315324 
    316       DO itide = 1, td%ncpt 
    317          z_sarg = z_arg * td%speed(itide) 
     325      ! Linear ramp on tidal component at open boundaries  
     326      zramp = 1. 
     327      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
     328 
     329      DO itide = 1, nb_harmo 
     330         z_sarg = z_arg * omega_tide(itide) 
    318331         z_cost(itide) = COS( z_sarg ) 
    319332         z_sist(itide) = SIN( z_sarg ) 
    320333      END DO 
    321334 
    322       DO itide = 1, td%ncpt 
    323          igrd=1                              ! SSH on tracer grid. 
    324          DO ib = 1, idx%nblenrim(igrd) 
    325             dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 
    326             !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2) 
     335      DO itide = 1, nb_harmo 
     336         igrd=1                              ! SSH on tracer grid 
     337         DO ib = 1, ilen0(igrd) 
     338            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 
    327339         END DO 
    328340         igrd=2                              ! U grid 
    329          DO ib=1, idx%nblenrim(igrd) 
    330             dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 
    331             !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 
     341         DO ib = 1, ilen0(igrd) 
     342            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)) 
    332343         END DO 
    333344         igrd=3                              ! V grid 
    334          DO ib=1, idx%nblenrim(igrd) 
    335             dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 
    336             !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 
     345         DO ib = 1, ilen0(igrd)  
     346            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)) 
    337347         END DO 
    338348      END DO 
    339349      ! 
    340       IF( nn_timing == 1 ) CALL timing_stop('tide_update') 
     350      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update') 
    341351      ! 
    342    END SUBROUTINE tide_update 
    343  
    344 !!gm  doctor naming of dummy argument variables!!!   and all local variables 
    345  
    346    SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 
    347       !!---------------------------------------------------------------------- 
    348       !!                   ***  SUBROUTINE uvset  *** 
    349       !!                      
    350       !! ** Purpose : - adjust tidal forcing for date factors 
    351       !!                 
    352       !!---------------------------------------------------------------------- 
    353       implicit none 
    354       INTEGER, INTENT( in ) ::   ihs     ! Start time hours 
    355       INTEGER, INTENT( in ) ::   iday    ! start time days 
    356       INTEGER, INTENT( in ) ::   imnth   ! start time month 
    357       INTEGER, INTENT( in ) ::   iyr     ! start time year 
    358       !! 
    359 !!gm  nc is jptides_max ???? 
    360       INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents 
    361       CHARACTER(len=8), DIMENSION(nc) :: cname 
    362       INTEGER                 ::   year, vd, ivdy, ndc, i, k 
    363       REAL(wp)                ::   ss, h, p, en, p1, rtd 
    364       REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction 
    365       REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction 
    366       REAL(wp), DIMENSION(nc) ::   u, v, zig 
    367       !! 
    368       DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   & 
    369          &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   & 
    370          &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      / 
    371       DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   & 
    372          &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   & 
    373          &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   / 
    374       !!---------------------------------------------------------------------- 
    375 ! 
    376 !     ihs             -  start time gmt on ... 
    377 !     iday/imnth/iyr  -  date   e.g.   12/10/87 
    378 ! 
    379       CALL vday(iday,imnth,iyr,ivdy) 
    380       vd = ivdy 
    381 ! 
    382 !rp   note change of year number for d. blackman shpen 
    383 !rp   if(iyr.ge.1000) year=iyr-1900 
    384 !rp   if(iyr.lt.1000) year=iyr 
    385       year = iyr 
    386 ! 
    387 !.....year  =  year of required data 
    388 !.....vd    =  day of required data..set up for 0000gmt day year 
    389       ndc = nc 
    390 !.....ndc   =  number of constituents allowed 
    391 !!gm use rpi ? 
    392       rtd = 360.0 / 6.2831852 
    393       DO i = 1, ndc 
    394          zig(i) = zig(i)*rtd 
    395          ! sigo(i)= zig(i) 
     352   END SUBROUTINE bdytide_update 
     353 
     354   SUBROUTINE tide_init_elevation( idx, td ) 
     355      !!---------------------------------------------------------------------- 
     356      !!                 ***  ROUTINE tide_init_elevation  *** 
     357      !!---------------------------------------------------------------------- 
     358      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices 
     359      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data 
     360      !! * Local declarations 
     361      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays) 
     362      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide 
     363      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices 
     364 
     365      igrd=1    
     366                              ! SSH on tracer grid. 
     367    
     368      ilen0(1) =  SIZE(td%ssh0(:,1,1)) 
     369 
     370      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 
     371 
     372      DO itide = 1, nb_harmo 
     373         DO ib = 1, ilen0(igrd) 
     374            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.) 
     375            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 
     376         END DO 
     377         DO ib = 1 , ilen0(igrd) 
     378            mod_tide(ib)=mod_tide(ib)*ftide(itide) 
     379            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     380         ENDDO 
     381         DO ib = 1 , ilen0(igrd) 
     382            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
     383            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     384         ENDDO 
    396385      END DO 
    397386 
    398 !!gm try to avoid RETURN  in F90 
    399       IF( year == 0 )   RETURN 
    400        
    401          CALL shpen( year, vd, ss, h , p , en, p1 ) 
    402          CALL ufset( p   , en, u , f ) 
    403          CALL vset ( ss  , h , p , en, p1, v ) 
    404          ! 
    405          DO k = 1, nc 
    406             z_vplu(k) = v(k) + u(k) 
    407             z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 
    408             DO WHILE( z_vplu(k) < 0    ) 
    409                z_vplu(k) = z_vplu(k) + 360.0 
    410             END DO 
    411             DO WHILE( z_vplu(k) > 360. ) 
    412                z_vplu(k) = z_vplu(k) - 360.0 
    413             END DO 
    414          END DO 
    415          ! 
    416       END SUBROUTINE uvset 
    417  
    418  
    419       SUBROUTINE vday( iday, imnth, iy, ivdy ) 
    420       !!---------------------------------------------------------------------- 
    421       !!                   *** SUBROUTINE vday  *** 
    422       !!                   
    423       !! ** Purpose : - adjust tidal forcing for date factors 
    424       !!                 
    425       !!---------------------------------------------------------------------- 
    426       INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ???? 
    427       INTEGER, INTENT(  out) ::   ivdy              ! ??? 
    428       !!  
    429       INTEGER ::   iyr 
    430       !!------------------------------------------------------------------------------ 
    431  
    432 !!gm   nday_year in day mode is the variable compiuted here, no? 
    433 !!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year 
    434  
    435       !calculate day number in year from day/month/year 
    436        if(imnth.eq.1) ivdy=iday 
    437        if(imnth.eq.2) ivdy=iday+31 
    438        if(imnth.eq.3) ivdy=iday+59 
    439        if(imnth.eq.4) ivdy=iday+90 
    440        if(imnth.eq.5) ivdy=iday+120 
    441        if(imnth.eq.6) ivdy=iday+151 
    442        if(imnth.eq.7) ivdy=iday+181 
    443        if(imnth.eq.8) ivdy=iday+212 
    444        if(imnth.eq.9) ivdy=iday+243 
    445        if(imnth.eq.10) ivdy=iday+273 
    446        if(imnth.eq.11) ivdy=iday+304 
    447        if(imnth.eq.12) ivdy=iday+334 
    448        iyr=iy 
    449        if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
    450        if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 
    451        if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
    452       ! 
    453    END SUBROUTINE vday 
    454  
    455 !!doctor norme for dummy arguments 
    456  
    457    SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 
    458       !!---------------------------------------------------------------------- 
    459       !!                    ***  SUBROUTINE shpen  *** 
    460       !!                    
    461       !! ** Purpose : - calculate astronomical arguments for tides 
    462       !!                this version from d. blackman 30 nove 1990 
    463       !! 
    464       !!---------------------------------------------------------------------- 
    465 !!gm add INTENT in, out or inout....    DOCTOR name.... 
    466 !!gm please do not use variable name with a single letter:  impossible to search in a code 
    467       INTEGER  ::   year,vd 
    468       REAL(wp) ::   s,h,p,en,p1 
    469       !! 
    470       INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd 
    471       REAL(wp) ::   cycle,t,td,delt(84),delta,deltat 
    472       !! 
    473       DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           & 
    474          &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    & 
    475          &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    & 
    476          &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    & 
    477          &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    & 
    478          &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    & 
    479          &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    & 
    480          &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    & 
    481          &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    & 
    482          &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    & 
    483          &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    & 
    484          &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    & 
    485          &       52.57 / 
    486       !!---------------------------------------------------------------------- 
    487  
    488       cycle = 360.0 
    489       ilc = 0 
    490       icent = year / 100 
    491       yr = year - icent * 100 
    492       t = icent - 20 
    493 ! 
    494 !     for the following equations 
    495 !     time origin is fixed at 00 hr of jan 1st,2000. 
    496 !     see notes by cartwright 
    497 ! 
    498 !!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90) 
    499 !!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    500       it = icent - 20 
    501       if (it) 1,2,2 
    502     1 iday = it/4 -it 
    503       go to 3 
    504     2 iday = (it+3)/4 - it 
    505 ! 
    506 !     t is in julian century 
    507 !     correction in gegorian calander where only century year divisible 
    508 !     by 4 is leap year. 
    509 ! 
    510     3 continue 
    511 ! 
    512       td = 0.0 
    513 ! 
    514 !!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    515       if (yr) 4,5,4 
    516 ! 
    517     4 iyd = 365*yr 
    518       ild = (yr-1)/4 
    519       if((icent - (icent/4)*4) .eq. 0) ilc = 1 
    520       td = iyd + ild + ilc 
    521 ! 
    522     5 td = td + iday + vd -1.0 - 0.5 
    523       t = t + (td/36525.0) 
    524 ! 
    525       ipos=year-1899 
    526       if (ipos .lt. 0) go to 7 
    527       if (ipos .gt. 83) go to 6 
    528 ! 
    529       delta = (delt(ipos+1)+delt(ipos))/2.0 
    530       go to 7 
    531 ! 
    532     6 delta= (65.0-50.5)/20.0*(year-1980)+50.5 
    533 ! 
    534     7 deltat = delta * 1.0e-6 
    535 ! 
    536 !!gm   precision of the computation   :  example for s it should be replace by: 
    537 !!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results 
    538       s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 
    539       h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat 
    540       p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat 
    541       en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat 
    542       p1  = 282.9384 + 1.7195     *t + 0.0005*t*t 
    543       ! 
    544       nn = s / cycle 
    545       s = s - nn * cycle 
    546       IF( s < 0.e0 )   s = s + cycle 
    547       ! 
    548       nn = h / cycle 
    549       h = h - cycle * nn 
    550       IF( h < 0.e0 )   h = h + cycle 
    551       ! 
    552       nn = p / cycle 
    553       p = p - cycle * nn 
    554       IF( p < 0.e0)   p = p + cycle 
    555       ! 
    556       nn = en / cycle 
    557       en = en - cycle * nn 
    558       IF( en < 0.e0 )   en = en + cycle 
    559       en = cycle - en 
    560       ! 
    561       nn = p1 / cycle 
    562       p1 = p1 - nn * cycle 
    563       ! 
    564    END SUBROUTINE shpen 
    565  
    566  
    567    SUBROUTINE ufset( p, cn, b, a ) 
    568       !!---------------------------------------------------------------------- 
    569       !!                    ***  SUBROUTINE ufset  *** 
    570       !!                     
    571       !! ** Purpose : - calculate nodal parameters for the tides 
    572       !!                 
    573       !!---------------------------------------------------------------------- 
    574 !!gm  doctor naming of dummy argument variables!!!   and all local variables 
    575 !!gm  nc is jptides_max ???? 
    576       integer nc 
    577       parameter (nc=15) 
    578       REAL(wp) p,cn 
    579       !! 
    580        
    581 !!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z" 
    582       REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 
    583       REAL(wp) ::   a(nc), b(nc) 
    584       INTEGER  ::   ndc, k 
    585       !!---------------------------------------------------------------------- 
    586  
    587       ndc = nc 
    588  
    589 !    a=f       ,  b =u 
    590 !    t is  zero as compared to tifa. 
    591 !! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module 
    592       rad = 6.2831852d0/360.0 
    593       pw = p  * rad 
    594       nw = cn * rad 
    595       w1 = cos(   nw ) 
    596       w2 = cos( 2*nw ) 
    597       w3 = cos( 3*nw ) 
    598       w4 = sin(   nw ) 
    599       w5 = sin( 2*nw ) 
    600       w6 = sin( 3*nw ) 
    601       w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   & 
    602          &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw ) 
    603       w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   & 
    604          &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw ) 
    605       ! 
    606       a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 
    607       b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 
    608       !   q1 
    609       a(2) = a(1) 
    610       b(2) = b(1) 
    611       !   o1 
    612       a(3) = 1.0 
    613       b(3) = 0.0 
    614       !   p1 
    615       a(4) = 1.0 
    616       b(4) = 0.0 
    617       !   s1 
    618       a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 
    619       b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 
    620       !   k1 
    621       a(6) =1.0004 -0.0373*w1+ 0.0002*w2 
    622       b(6) = -0.0374*w4 
    623       !  2n2 
    624       a(7) = a(6) 
    625       b(7) = b(6) 
    626       !  mu2 
    627       a(8) = a(6) 
    628       b(8) = b(6) 
    629       !   n2 
    630       a(9) = a(6) 
    631       b(9) = b(6) 
    632       !  nu2 
    633       a(10) = a(6) 
    634       b(10) = b(6) 
    635       !   m2 
    636       a(11) = SQRT( w7 * w7 + w8 * w8 ) 
    637       b(11) = ATAN( w8 / w7 ) 
    638 !!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ???? 
    639       IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992 
    640       !   l2 
    641       a(12) = 1.0 
    642       b(12) = 0.0 
    643       !   t2 
    644       a(13)= a(12) 
    645       b(13)= b(12) 
    646       !   s2 
    647       a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 
    648       b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 
    649       !   k2 
    650       a(15) = a(6)*a(6) 
    651       b(15) = 2*b(6) 
    652       !   m4 
    653 !!gm  old coding,  remove GOTO and label of lines 
    654 !!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    655       DO 40 k = 1,ndc 
    656          b(k) = b(k)/rad 
    657 32       if (b(k)) 34,35,35 
    658 34       b(k) = b(k) + 360.0 
    659          go to 32 
    660 35       if (b(k)-360.0) 40,37,37 
    661 37       b(k) = b(k)-360.0 
    662          go to 35 
    663 40    continue 
    664       ! 
    665    END SUBROUTINE ufset 
    666     
    667  
    668    SUBROUTINE vset( s,h,p,en,p1,v) 
    669       !!---------------------------------------------------------------------- 
    670       !!                    ***  SUBROUTINE vset  *** 
    671       !!                    
    672       !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 
    673       !!                 
    674       !!---------------------------------------------------------------------- 
    675 !!gm  doctor naming of dummy argument variables!!!   and all local variables 
    676 !!gm  nc is jptides_max ???? 
    677 !!gm  en argument is not used: suppress it ? 
    678       integer nc 
    679       parameter (nc=15) 
    680       real(wp) s,h,p,en,p1 
    681       real(wp)   v(nc) 
    682       !! 
    683       integer ndc, k 
    684       !!---------------------------------------------------------------------- 
    685  
    686       ndc = nc 
    687       !   v s  are computed here. 
    688       v(1) =-3*s +h +p +270      ! Q1 
    689       v(2) =-2*s +h +270.0       ! O1 
    690       v(3) =-h +270              ! P1 
    691       v(4) =180                  ! S1 
    692       v(5) =h +90.0              ! K1 
    693       v(6) =-4*s +2*h +2*p       ! 2N2 
    694       v(7) =-4*(s-h)             ! MU2 
    695       v(8) =-3*s +2*h +p         ! N2 
    696       v(9) =-3*s +4*h -p         ! MU2 
    697       v(10) =-2*s +2*h           ! M2 
    698       v(11) =-s +2*h -p +180     ! L2  
    699       v(12) =-h +p1              ! T2 
    700       v(13) =0                   ! S2 
    701       v(14) =h+h                 ! K2 
    702       v(15) =2*v(10)             ! M4 
    703 ! 
    704 !!gm  old coding,  remove GOTO and label of lines 
    705 !!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    706       do 72 k = 1, ndc 
    707 69       if( v(k) )   70,71,71 
    708 70       v(k) = v(k)+360.0 
    709          go to 69 
    710 71       if( v(k) - 360.0 )   72,73,73 
    711 73       v(k) = v(k)-360.0 
    712          go to 71 
    713 72    continue 
    714       ! 
    715    END SUBROUTINE vset 
    716  
     387      DEALLOCATE(mod_tide,phi_tide) 
     388 
     389   END SUBROUTINE tide_init_elevation 
     390 
     391   SUBROUTINE tide_init_velocities( idx, td ) 
     392      !!---------------------------------------------------------------------- 
     393      !!                 ***  ROUTINE tide_init_elevation  *** 
     394      !!---------------------------------------------------------------------- 
     395      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices 
     396      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data 
     397      !! * Local declarations 
     398      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays) 
     399      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide 
     400      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices 
     401 
     402      ilen0(2) =  SIZE(td%u0(:,1,1)) 
     403      ilen0(3) =  SIZE(td%v0(:,1,1)) 
     404 
     405      igrd=2                                 ! U grid. 
     406 
     407      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 
     408 
     409      DO itide = 1, nb_harmo 
     410         DO ib = 1, ilen0(igrd) 
     411            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.) 
     412            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 
     413         END DO 
     414         DO ib = 1, ilen0(igrd) 
     415            mod_tide(ib)=mod_tide(ib)*ftide(itide) 
     416            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     417         ENDDO 
     418         DO ib = 1, ilen0(igrd) 
     419            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
     420            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     421         ENDDO 
     422      END DO 
     423 
     424      DEALLOCATE(mod_tide,phi_tide) 
     425 
     426      igrd=3                                 ! V grid. 
     427 
     428      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 
     429 
     430      DO itide = 1, nb_harmo 
     431         DO ib = 1, ilen0(igrd) 
     432            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.) 
     433            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 
     434         END DO 
     435         DO ib = 1, ilen0(igrd) 
     436            mod_tide(ib)=mod_tide(ib)*ftide(itide) 
     437            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     438         ENDDO 
     439         DO ib = 1, ilen0(igrd) 
     440            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
     441            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     442         ENDDO 
     443      END DO 
     444 
     445      DEALLOCATE(mod_tide,phi_tide) 
     446 
     447  END SUBROUTINE tide_init_velocities 
    717448#else 
    718449   !!---------------------------------------------------------------------- 
    719450   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides 
    720451   !!---------------------------------------------------------------------- 
    721 !!gm  are you sure we need to define filtide and tide_cpt ? 
    722    CHARACTER(len=80), PUBLIC               ::   filtide                !: Filename root for tidal input files 
    723    CHARACTER(len=4 ), PUBLIC, DIMENSION(1) ::   tide_cpt               !: Names of tidal components used. 
    724  
    725452CONTAINS 
    726    SUBROUTINE tide_init                ! Empty routine 
    727    END SUBROUTINE tide_init 
    728    SUBROUTINE tide_data                ! Empty routine 
    729    END SUBROUTINE tide_data 
    730    SUBROUTINE tide_update( kt, kit )   ! Empty routine 
    731       WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit 
    732    END SUBROUTINE tide_update 
     453   SUBROUTINE bdytide_init             ! Empty routine 
     454      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?' 
     455   END SUBROUTINE bdytide_init 
     456   SUBROUTINE bdytide_update( kt, jit )   ! Empty routine 
     457      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 
     458   END SUBROUTINE bdytide_update 
    733459#endif 
    734460 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r3294 r3651  
    2323   USE in_out_manager  ! I/O manager 
    2424 
     25 
    2526   IMPLICIT NONE 
    2627   PRIVATE 
    2728 
    2829   PUBLIC bdy_tra      ! routine called in tranxt.F90  
     30   PUBLIC bdy_tra_dmp  ! routine called in step.F90  
    2931 
    3032   !!---------------------------------------------------------------------- 
    … …  
    5355         CASE(jp_frs) 
    5456            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     57         CASE(2) 
     58            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     59         CASE(3) 
     60            CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
     61         CASE(4) 
     62            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    5563         CASE DEFAULT 
    5664            CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' ) 
    5765         END SELECT 
    5866      ENDDO 
     67      ! 
     68      ! Boundary points should be updated 
     69      IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) 
     70      IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    5971 
    6072   END SUBROUTINE bdy_tra 
    … …  
    90102      END DO  
    91103      ! 
    92       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
    93       ! 
    94104      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    95105      ! 
    … …  
    97107      ! 
    98108   END SUBROUTINE bdy_tra_frs 
    99     
     109   
     110   SUBROUTINE bdy_tra_spe( idx, dta, kt ) 
     111      !!---------------------------------------------------------------------- 
     112      !!                 ***  SUBROUTINE bdy_tra_frs  *** 
     113      !!                     
     114      !! ** Purpose : Apply a specified value for tracers at open boundaries. 
     115      !!  
     116      !!---------------------------------------------------------------------- 
     117      INTEGER,         INTENT(in) ::   kt 
     118      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     119      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     120      !!  
     121      REAL(wp) ::   zwgt           ! boundary weight 
     122      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     123      INTEGER  ::   ii, ij         ! 2D addresses 
     124      !!---------------------------------------------------------------------- 
     125      ! 
     126      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_spe') 
     127      ! 
     128      igrd = 1                       ! Everything is at T-points here 
     129      DO ib = 1, idx%nblenrim(igrd) 
     130         ii = idx%nbi(ib,igrd) 
     131         ij = idx%nbj(ib,igrd) 
     132         DO ik = 1, jpkm1 
     133            tsa(ii,ij,ik,jp_tem) = dta%tem(ib,ik) * tmask(ii,ij,ik) 
     134            tsa(ii,ij,ik,jp_sal) = dta%sal(ib,ik) * tmask(ii,ij,ik) 
     135         END DO 
     136      END DO 
     137      ! 
     138      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     139      ! 
     140      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_spe') 
     141      ! 
     142   END SUBROUTINE bdy_tra_spe 
     143 
     144   SUBROUTINE bdy_tra_nmn( idx, dta, kt ) 
     145      !!---------------------------------------------------------------------- 
     146      !!                 ***  SUBROUTINE bdy_tra_nmn  *** 
     147      !!                     
     148      !! ** Purpose : Duplicate the value for tracers at open boundaries. 
     149      !!  
     150      !!---------------------------------------------------------------------- 
     151      INTEGER,         INTENT(in) ::   kt 
     152      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     153      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     154      !!  
     155      REAL(wp) ::   zwgt           ! boundary weight 
     156      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     157      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses 
     158      !!---------------------------------------------------------------------- 
     159      ! 
     160      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_nmn') 
     161      ! 
     162      igrd = 1                       ! Everything is at T-points here 
     163      DO ib = 1, idx%nblenrim(igrd) 
     164         ii = idx%nbi(ib,igrd) 
     165         ij = idx%nbj(ib,igrd) 
     166         DO ik = 1, jpkm1 
     167            ! search the sense of the gradient 
     168            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
     169            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
     170            IF ( zcoef1+zcoef2 == 0) THEN 
     171               ! corner 
     172               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik) 
     173               tsa(ii,ij,ik,jp_tem) = tsa(ii-1,ij  ,ik,jp_tem) * tmask(ii-1,ij  ,ik) + & 
     174                 &                    tsa(ii+1,ij  ,ik,jp_tem) * tmask(ii+1,ij  ,ik) + & 
     175                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + & 
     176                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik) 
     177               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     178               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + & 
     179                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + & 
     180                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + & 
     181                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik) 
     182               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     183            ELSE 
     184               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
     185               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     186               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 
     187               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) 
     188            ENDIF 
     189         END DO 
     190      END DO 
     191      ! 
     192      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     193      ! 
     194      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_nmn') 
     195      ! 
     196   END SUBROUTINE bdy_tra_nmn 
     197 
     198   SUBROUTINE bdy_tra_rnf( idx, dta, kt ) 
     199      !!---------------------------------------------------------------------- 
     200      !!                 ***  SUBROUTINE bdy_tra_rnf  *** 
     201      !!                     
     202      !! ** Purpose : Apply the runoff values for tracers at open boundaries: 
     203      !!                  - specified to 0.1 PSU for the salinity 
     204      !!                  - duplicate the value for the temperature 
     205      !!  
     206      !!---------------------------------------------------------------------- 
     207      INTEGER,         INTENT(in) ::   kt 
     208      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     209      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     210      !!  
     211      REAL(wp) ::   zwgt           ! boundary weight 
     212      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     213      INTEGER  ::   ii, ij, ip, jp ! 2D addresses 
     214      !!---------------------------------------------------------------------- 
     215      ! 
     216      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_rnf') 
     217      ! 
     218      igrd = 1                       ! Everything is at T-points here 
     219      DO ib = 1, idx%nblenrim(igrd) 
     220         ii = idx%nbi(ib,igrd) 
     221         ij = idx%nbj(ib,igrd) 
     222         DO ik = 1, jpkm1 
     223            ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
     224            jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     225            tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii,ij,ik) 
     226            tsa(ii,ij,ik,jp_sal) =                        0.1 * tmask(ii,ij,ik) 
     227         END DO 
     228      END DO 
     229      ! 
     230      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     231      ! 
     232      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_rnf') 
     233      ! 
     234   END SUBROUTINE bdy_tra_rnf 
     235 
     236   SUBROUTINE bdy_tra_dmp( kt ) 
     237      !!---------------------------------------------------------------------- 
     238      !!                 ***  SUBROUTINE bdy_tra_dmp  *** 
     239      !!                     
     240      !! ** Purpose : Apply damping for tracers at open boundaries. 
     241      !!  
     242      !!---------------------------------------------------------------------- 
     243      INTEGER,         INTENT(in) ::   kt 
     244      !!  
     245      REAL(wp) ::   zwgt           ! boundary weight 
     246      REAL(wp) ::   zta, zsa, ztime 
     247      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     248      INTEGER  ::   ii, ij         ! 2D addresses 
     249      INTEGER  ::   ib_bdy         ! Loop index 
     250      !!---------------------------------------------------------------------- 
     251      ! 
     252      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_dmp') 
     253      ! 
     254      DO ib_bdy=1, nb_bdy 
     255         IF ( ln_tra_dmp(ib_bdy) ) THEN 
     256            igrd = 1                       ! Everything is at T-points here 
     257            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     258               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     259               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     260               zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd) 
     261               DO ik = 1, jpkm1 
     262                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik) 
     263                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik) 
     264                  tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta 
     265                  tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa 
     266               END DO 
     267            END DO 
     268         ENDIF 
     269      ENDDO 
     270      ! 
     271      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp') 
     272      ! 
     273   END SUBROUTINE bdy_tra_dmp 
     274  
    100275#else 
    101276   !!---------------------------------------------------------------------- 
    … …  
    106281      WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt 
    107282   END SUBROUTINE bdy_tra 
     283 
     284   SUBROUTINE bdy_tra_dmp(kt)      ! Empty routine 
     285      WRITE(*,*) 'bdy_tra_dmp: You should not have seen this print! error?', kt 
     286   END SUBROUTINE bdy_tra_dmp 
     287 
    108288#endif 
    109289 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r3294 r3651  
    403403         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    404404         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
    405          IF ( ln_tide_pot ) CALL upd_tide( kt, jn ) 
     405         IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 
    406406 
    407407         !                                                !* after ssh_e 
    … …  
    453453                  ENDIF 
    454454                  ! add tidal astronomical forcing 
    455                   IF ( ln_tide_pot ) THEN  
     455                  IF ( ln_tide_pot .AND. lk_tide ) THEN  
    456456                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    457457                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    … …  
    503503                  ENDIF 
    504504                  ! add tidal astronomical forcing 
    505                   IF ( ln_tide_pot ) THEN 
     505                  IF ( ln_tide_pot .AND. lk_tide ) THEN 
    506506                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    507507                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    … …  
    550550                  ENDIF 
    551551                  ! add tidal astronomical forcing 
    552                   IF ( ln_tide_pot ) THEN 
     552                  IF ( ln_tide_pot .AND. lk_tide ) THEN 
    553553                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    554554                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r3294 r3651  
    106106   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    107107      & ld_velav     !: Velocity data is daily averaged 
     108   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
     109      & ld_sstnight  !: SST observation corresponds to night mean 
    108110 
    109111   !!---------------------------------------------------------------------- 
    … …  
    737739         ALLOCATE(sstdata(nsstsets)) 
    738740         ALLOCATE(sstdatqc(nsstsets)) 
     741         ALLOCATE(ld_sstnight(nsstsets)) 
    739742         sstdata(:)%nsurf=0 
    740          sstdatqc(:)%nsurf=0          
     743         sstdatqc(:)%nsurf=0     
     744         ld_sstnight(:)=.false. 
    741745 
    742746         nsstsets = 0 
    … …  
    745749 
    746750            nsstsets = nsstsets + 1 
     751 
     752            ld_sstnight(nsstsets) = .TRUE. 
    747753 
    748754            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 
    … …  
    757763         
    758764            nsstsets = nsstsets + 1 
     765 
     766            ld_sstnight(nsstsets) = .TRUE. 
    759767           
    760768            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 
    … …  
    774782             
    775783               nsstsets = nsstsets + 1 
     784 
     785               ld_sstnight(nsstsets) = .TRUE. 
    776786             
    777787               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 
    … …  
    10921102      IF ( ln_sst ) THEN 
    10931103         DO jsstset = 1, nsstsets 
    1094             CALL obs_sst_opt( sstdatqc(jsstset),                 & 
    1095                &              kstp, jpi, jpj, nit000, tsn(:,:,1,jp_tem), & 
    1096                &              tmask(:,:,1), n2dint ) 
     1104            CALL obs_sst_opt( sstdatqc(jsstset),                & 
     1105               &              kstp, jpi, jpj, nit000, idaystp,  & 
     1106               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  & 
     1107               &              n2dint, ld_sstnight(jsstset) ) 
    10971108         END DO 
    10981109      ENDIF 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r2715 r3651  
    614614   END SUBROUTINE obs_sla_opt 
    615615 
    616    SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, & 
    617       &                    psstn, psstmask, k2dint ) 
    618  
     616   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 
     617      &                    psstn, psstmask, k2dint, ld_nightav ) 
    619618      !!----------------------------------------------------------------------- 
    620619      !! 
    … …  
    647646      !! * Modules used 
    648647      USE obs_surf_def  ! Definition of storage space for surface observations 
     648      USE sbcdcy 
    649649 
    650650      IMPLICIT NONE 
    … …  
    659659                                       !   (kit000-1 = restart time) 
    660660      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
     661      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day   
    661662      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    662663         & psstn,  &    ! Model SST field 
    663664         & psstmask     ! Land-sea mask 
    664           
     665 
    665666      !! * Local declarations 
    666667      INTEGER :: ji 
    … …  
    670671      INTEGER :: isst 
    671672      INTEGER :: iobs 
     673      INTEGER :: idayend 
    672674      REAL(KIND=wp) :: zlam 
    673675      REAL(KIND=wp) :: zphi 
    674676      REAL(KIND=wp) :: zext(1), zobsmask(1) 
     677      REAL(KIND=wp) :: zdaystp 
     678      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
     679         & icount_sstnight,      & 
     680         & imask_night 
     681      REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
     682         & zintmp, & 
     683         & zouttmp, &  
     684         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
    675685      REAL(kind=wp), DIMENSION(2,2,1) :: & 
    676686         & zweig 
    … …  
    678688         & zmask, & 
    679689         & zsstl, & 
     690         & zsstm, & 
    680691         & zglam, & 
    681692         & zgphi 
    … …  
    683694         & igrdi, & 
    684695         & igrdj 
     696      LOGICAL, INTENT(IN) :: ld_nightav 
    685697 
    686698      !----------------------------------------------------------------------- 
    … …  
    690702      inrc = kt - kit000 + 2 
    691703      isst = sstdatqc%nsstp(inrc) 
     704 
     705      IF ( ld_nightav ) THEN 
     706 
     707      ! Initialize array for night mean 
     708 
     709      IF ( kt .EQ. 0 ) THEN 
     710         ALLOCATE ( icount_sstnight(kpi,kpj) ) 
     711         ALLOCATE ( imask_night(kpi,kpj) ) 
     712         ALLOCATE ( zintmp(kpi,kpj) ) 
     713         ALLOCATE ( zouttmp(kpi,kpj) ) 
     714         ALLOCATE ( zmeanday(kpi,kpj) ) 
     715         nday_qsr = -1   ! initialisation flag for nbc_dcy 
     716      ENDIF 
     717 
     718      ! Initialize daily mean for first timestep 
     719      idayend = MOD( kt - kit000 + 1, kdaystp ) 
     720 
     721      ! Added kt == 0 test to catch restart case  
     722      IF ( idayend == 1 .OR. kt == 0) THEN 
     723         IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 
     724         DO jj = 1, jpj 
     725            DO ji = 1, jpi 
     726               sstdatqc%vdmean(ji,jj) = 0.0 
     727               zmeanday(ji,jj) = 0.0 
     728               icount_sstnight(ji,jj) = 0 
     729            END DO 
     730         END DO 
     731      ENDIF 
     732 
     733      zintmp(:,:) = 0.0 
     734      zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 
     735      imask_night(:,:) = INT( zouttmp(:,:) ) 
     736 
     737      DO jj = 1, jpj 
     738         DO ji = 1, jpi 
     739            ! Increment the temperature field for computing night mean and counter 
     740            sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  & 
     741                   &                        + psstn(ji,jj)*imask_night(ji,jj) 
     742            zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj) 
     743            icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 
     744         END DO 
     745      END DO 
     746    
     747      ! Compute the daily mean at the end of day 
     748 
     749      zdaystp = 1.0 / REAL( kdaystp ) 
     750 
     751      IF ( idayend == 0 ) THEN  
     752         DO jj = 1, jpj 
     753            DO ji = 1, jpi 
     754               ! Test if "no night" point 
     755               IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 
     756                  sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 
     757                    &                        / icount_sstnight(ji,jj)  
     758               ELSE 
     759                  sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
     760               ENDIF 
     761            END DO 
     762         END DO 
     763      ENDIF 
     764 
     765      ENDIF 
    692766 
    693767      ! Get the data for interpolation 
    … …  
    722796      CALL obs_int_comm_2d( 2, 2, isst, & 
    723797         &                  igrdi, igrdj, psstn, zsstl ) 
    724        
     798 
     799      ! At the end of the day get interpolated means 
     800      IF ( idayend == 0 .AND. ld_nightav ) THEN 
     801 
     802         ALLOCATE( & 
     803            & zsstm(2,2,isst)  & 
     804            & ) 
     805 
     806         CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 
     807            &               sstdatqc%vdmean(:,:), zsstm ) 
     808 
     809      ENDIF 
     810 
    725811      ! Loop over observations 
    726812 
    … …  
    756842             
    757843         ! Interpolate the model SST to the observation point  
    758          CALL obs_int_h2d( 1, 1,      & 
     844 
     845         IF ( ld_nightav ) THEN 
     846 
     847           IF ( idayend == 0 )  THEN 
     848               ! Daily averaged/diurnal cycle of SST  data 
     849               CALL obs_int_h2d( 1, 1,      &  
     850                     &              zweig, zsstm(:,:,iobs), zext ) 
     851            ELSE  
     852               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     & 
     853                     &           ' number of night SST data should' // & 
     854                     &           ' only occur at the end of a given day' ) 
     855            ENDIF 
     856 
     857         ELSE 
     858 
     859            CALL obs_int_h2d( 1, 1,      & 
    759860            &              zweig, zsstl(:,:,iobs),  zext ) 
     861 
     862         ENDIF 
    760863          
    761864         sstdatqc%rmod(jobs,1) = zext(1) 
    … …  
    772875         & zsstl  & 
    773876         & ) 
     877 
     878      ! At the end of the day also get interpolated means 
     879      IF ( idayend == 0 .AND. ld_nightav ) THEN 
     880         DEALLOCATE( & 
     881            & zsstm  & 
     882            & ) 
     883      ENDIF 
    774884       
    775885      sstdatqc%nsurfup = sstdatqc%nsurfup + isst 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_seaice.F90

    r2287 r3651  
    326326         &               iindx   ) 
    327327       
    328       CALL obs_surf_alloc( seaicedata, iobs, kvars, kextr, kstp ) 
     328      CALL obs_surf_alloc( seaicedata, iobs, &  
     329                           kvars, kextr, kstp, jpi, jpj ) 
    329330       
    330331      ! * Read obs/positions, QC, all variable and assign to seaicedata 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sla.F90

    r2287 r3651  
    391391         &               iindx   ) 
    392392       
    393       CALL obs_surf_alloc( sladata, iobs, kvars, kextr, kstp ) 
     393      CALL obs_surf_alloc( sladata, iobs, kvars, kextr, & 
     394         &                 jpi, jpj, kstp ) 
    394395       
    395396      ! * Read obs/positions, QC, all variable and assign to sladata 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sst.F90

    r2287 r3651  
    326326         &               iindx   ) 
    327327       
    328       CALL obs_surf_alloc( sstdata, iobs, kvars, kextr, kstp ) 
     328      CALL obs_surf_alloc( sstdata, iobs, kvars, kextr, kstp, jpi, jpj ) 
    329329       
    330330      ! * Read obs/positions, QC, all variable and assign to sstdata 
    … …  
    701701      ! Allocate obs_surf data structure for time sorted data 
    702702          
    703       CALL obs_surf_alloc( sstdata, inumobs, kvars, kextra, kstp ) 
     703      CALL obs_surf_alloc( sstdata, inumobs, kvars, kextra, kstp, jpi, jpj ) 
    704704 
    705705      pjul = pjulini + 1 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

    r2287 r3651  
    4747      INTEGER :: nextra     !: Number of extra fields at observation points 
    4848      INTEGER :: nstp       !: Number of time steps 
     49      INTEGER :: npi        !: Number of 3D grid points 
     50      INTEGER :: npj 
    4951      INTEGER :: nsurfup    !: Observation counter used in obs_oper 
    5052 
    … …  
    7981         & rext           !: Extra fields interpolated to observation points 
    8082 
     83      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
     84         & vdmean         !: Time averaged of model field 
     85 
    8186      ! Arrays with size equal to the number of time steps in the window 
    8287 
    … …  
    103108CONTAINS 
    104109    
    105    SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp ) 
     110   SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj ) 
    106111      !!---------------------------------------------------------------------- 
    107112      !!                     ***  ROUTINE obs_surf_alloc  *** 
    … …  
    120125      INTEGER, INTENT(IN) :: kextra  ! Number of extra fields at observation points 
    121126      INTEGER, INTENT(IN) :: kstp    ! Number of time steps 
     127      INTEGER, INTENT(IN) :: kpi     ! Number of 3D grid points 
     128      INTEGER, INTENT(IN) :: kpj 
    122129 
    123130      !!* Local variables 
    … …  
    131138      surf%nvar     = kvar 
    132139      surf%nstp     = kstp 
     140      surf%npi      = kpi 
     141      surf%npj      = kpj 
    133142       
    134143      ! Allocate arrays of number of surface data size 
    … …  
    174183         & ) 
    175184 
     185      ! Allocate arrays of size number of grid points 
     186 
     187      ALLOCATE( & 
     188         & surf%vdmean(kpi,kpj) & 
     189         & ) 
     190 
    176191      ! Set defaults for compression indices 
    177192       
    … …  
    242257         & ) 
    243258 
     259      ! Deallocate arrays of size number of grid points size times 
     260      ! number of variables 
     261 
     262      DEALLOCATE( & 
     263         & surf%vdmean & 
     264         & ) 
     265 
    244266      ! Deallocate arrays of number of time step size 
    245267 
    … …  
    300322      IF ( lallocate ) THEN 
    301323         CALL obs_surf_alloc( newsurf,  insurf, surf%nvar, & 
    302             & surf%nextra, surf%nstp ) 
     324            & surf%nextra, surf%nstp, surf%npi, surf%npj ) 
    303325      ENDIF 
    304326 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r3294 r3651  
    667667      !!---------------------------------------------------------------------- 
    668668#if defined key_bdy 
    669       USE bdy_oce, ONLY:  dta_global         ! workspace to read in global data arrays 
     669      USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays 
    670670#endif  
    671671 
    … …  
    681681      INTEGER                                 ::   ilendta  ! length of data in file 
    682682      INTEGER                                 ::   idvar    ! variable ID 
    683       INTEGER                                 ::   ib, ik   ! loop counters 
     683      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters 
    684684      INTEGER                                 ::   ierr 
    685       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read ! work space for global data 
     685      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data 
    686686      !!--------------------------------------------------------------------- 
    687687             
    688 #if defined key_bdy 
    689       dta_read => dta_global 
    690 #endif 
    691  
    692688      ipi = SIZE( dta, 1 ) 
    693689      ipj = 1 
    … …  
    696692      idvar   = iom_varid( num, clvar ) 
    697693      ilendta = iom_file(num)%dimsz(1,idvar) 
     694 
     695#if defined key_bdy 
     696      ipj = iom_file(num)%dimsz(2,idvar) 
     697      IF (ipj == 1) THEN ! we assume that this is a structured open boundary file 
     698         dta_read => dta_global 
     699      ELSE 
     700         dta_read => dta_global2 
     701      ENDIF 
     702#endif 
     703 
     704      IF ((lwp).AND.(ipi>ilendta*ipj)) CALL ctl_stop('fld_map problem with file dimension') 
     705 
    698706      IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta 
    699707      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 
     708 
    700709 
    701710      SELECT CASE( ipk ) 
    … …  
    706715      END SELECT 
    707716      ! 
    708       DO ib = 1, ipi 
    709          DO ik = 1, ipk 
    710             dta(ib,1,ik) =  dta_read(map(ib),1,ik) 
     717      IF (ipj==1) THEN 
     718         DO ib = 1, ipi 
     719            DO ik = 1, ipk 
     720               dta(ib,1,ik) =  dta_read(map(ib),1,ik) 
     721            END DO 
    711722         END DO 
    712       END DO 
     723      ELSE ! we assume that this is a structured open boundary file 
     724         DO ib = 1, ipi 
     725            jj=1+floor(REAL(map(ib)-1)/REAL(ilendta)) 
     726            ji=map(ib)-(jj-1)*ilendta 
     727            DO ik = 1, ipk 
     728               dta(ib,1,ik) =  dta_read(ji,jj,ik) 
     729            END DO 
     730         END DO 
     731      ENDIF 
    713732 
    714733   END SUBROUTINE fld_map 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90

    r3294 r3651  
    6666      !! 
    6767      INTEGER            ::   ierror  ! local integer  
    68       REAL(wp)           ::   zpref   ! local scalar 
    6968      !! 
    7069      CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files 
    7170      TYPE(FLD_N)        ::  sn_apr   ! informations about the fields to be read 
    7271      !! 
    73       NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr 
     72      NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr, rpref, ln_apr_obc 
    7473      !!---------------------------------------------------------------------- 
    7574      ! 
    … …  
    113112         ! 
    114113         !                                            !* control check 
    115          IF( ln_apr_obc  )   & 
    116             CALL ctl_stop( 'sbc_apr: inverse barometer added to OBC ssh data not yet implemented ' ) 
    117          IF( ln_apr_obc .AND. .NOT. lk_obc )   & 
    118             CALL ctl_stop( 'sbc_apr: add inverse barometer to OBC requires to use key_obc' ) 
     114         IF ( ln_apr_obc  ) THEN 
     115            IF(lwp) WRITE(numout,*) '         Inverse barometer added to OBC ssh data' 
     116         ENDIF 
    119117         IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts )   & 
    120118            CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' ) 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90

    r3294 r3651  
    4949 
    5050 
    51    FUNCTION sbc_dcy( pqsrin ) RESULT( zqsrout ) 
     51   FUNCTION sbc_dcy( pqsrin, l_mask ) RESULT( zqsrout ) 
    5252      !!---------------------------------------------------------------------- 
    5353      !!                  ***  ROUTINE sbc_dcy  *** 
    … …  
    6363      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590. 
    6464      !!---------------------------------------------------------------------- 
     65      LOGICAL, OPTIONAL, INTENT(in) :: l_mask ! use the routine for night mask computation 
    6566      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux  
    6667      !! 
    6768      INTEGER  ::   ji, jj                                       ! dummy loop indices 
     69      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 
    6870      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad  
    6971      REAL(wp) ::   zlo, zup, zlousd, zupusd 
    7072      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos 
    7173      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest 
     74      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2 
    7275      REAL(wp), DIMENSION(jpi,jpj) ::   zqsrout                  ! output QSR flux with diurnal cycle 
    7376      !---------------------------statement functions------------------------ 
    … …  
    9093      zup = zlo + ( REAL(nn_fsbc, wp)     * rdttra(1) ) / rday 
    9194      !                                           
    92       IF( nday_qsr == -1 ) THEN       ! first time step only                
     95      IF( nday_qsr == -1 ) THEN       ! first time step only   
    9396         IF(lwp) THEN 
    9497            WRITE(numout,*) 
    … …  
    120123         zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) ) 
    121124         ! Compute A and B needed to compute the time integral of the diurnal cycle 
    122          
     125 
    123126         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad ) 
    124127         DO jj = 1, jpj 
    … …  
    129132            END DO   
    130133         END DO   
    131  
    132134         ! Compute the time of dawn and dusk 
    133135 
    … …  
    156158         rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp ) 
    157159         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp ) 
    158  
    159160         !     2.2 Compute the scalling function: 
    160161         !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk 
    … …  
    185186         ! 
    186187      ENDIF  
    187  
    188188         !     3. update qsr with the diurnal cycle 
    189189         !     ------------------------------------ 
    190190 
     191      imask_night(:,:) = 0 
    191192      DO jj = 1, jpj 
    192193         DO ji = 1, jpi 
     194            ztmpm = 0.0 
    193195            IF( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h 
    194196               ! 
    … …  
    200202                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    201203                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     204                  ztmpm = zupusd - zlousd 
     205                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
    202206                  ! 
    203207               ELSE                                         ! day time in two parts 
    … …  
    205209                  zupusd = MIN(zup, rdusk(ji,jj)) 
    206210                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     211                  ztmpm1=zupusd-zlousd 
    207212                  zlousd = MAX(zlo, rdawn(ji,jj)) 
    208213                  zupusd = MAX(zup, rdawn(ji,jj)) 
    209214                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     215                  ztmpm2 =zupusd-zlousd 
    210216                  ztmp = ztmp1 + ztmp2 
     217                  ztmpm = ztmpm1 + ztmpm2 
    211218                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     219                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
    212220               ENDIF 
    213221            ELSE                                   ! 24h light or 24h night 
    … …  
    216224                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    217225                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     226                  imask_night(ji,jj) = 0 
    218227                  ! 
    219228               ELSE                                         ! No day 
    220229                  zqsrout(ji,jj) = 0.e0 
     230                  imask_night(ji,jj) = 1 
    221231               ENDIF 
    222232            ENDIF 
    223233         END DO   
    224234      END DO   
     235      ! 
     236      IF ( PRESENT(l_mask) .AND. l_mask ) THEN 
     237         zqsrout(:,:) = float(imask_night(:,:)) 
     238      ENDIF 
    225239      ! 
    226240      IF( nn_timing == 1 )  CALL timing_stop('sbc_dcy') 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90

    r3294 r3651  
    1414  USE daymod 
    1515  USE dynspg_oce 
    16   USE tide_mod 
     16  USE tideini 
    1717  USE iom 
    1818 
    … …  
    2121 
    2222  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro 
    23   LOGICAL, PUBLIC :: ln_tide_pot = .false. 
     23 
    2424#if defined key_tide 
    2525 
    2626  LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE. 
    27  
    28   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide  
    29  
    30   REAL(wp), ALLOCATABLE, DIMENSION(:) ::  & 
    31        v0tide,      & 
    32        utide,       & 
    33        ftide 
    34  
    3527  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot,phi_pot 
    36  
    37   INTEGER, PUBLIC :: nb_harmo 
    38   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide 
    39   INTEGER, PUBLIC :: nn_tide, kt_tide 
    40  
    4128  !!--------------------------------------------------------------------------------- 
    4229  !!   OPA 9.0 , LODYC-IPSL  (2003) 
    … …  
    5138    !! * Arguments 
    5239    INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
    53     !! * Local declarations 
    54     INTEGER  :: jk,ji 
    55     CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 
    5640    !!---------------------------------------------------------------------- 
    5741 
    58     NAMELIST/nam_tide/ln_tide_pot,nb_harmo,clname,nn_tide 
     42    IF ( kt == nit000 .AND. .NOT. lk_dynspg_ts )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' ) 
    5943 
    60     IF ( kt == nit000 ) THEN 
     44    IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN 
     45      ! 
     46      kt_tide = kt 
    6147 
    62        IF( .NOT. lk_dynspg_ts  )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' ) 
     48      IF(lwp) THEN 
     49         WRITE(numout,*) 
     50         WRITE(numout,*) 'sbc_tide : (re)Initialization of the tidal potential at kt=',kt 
     51         WRITE(numout,*) '~~~~~~~ ' 
     52      ENDIF 
    6353 
    64     ! Read Namelist nam_tide 
     54      IF(lwp) THEN 
     55         IF ( kt == nit000 ) WRITE(numout,*) 'Apply astronomical potential : ln_tide_pot =', ln_tide_pot 
     56         CALL flush(numout) 
     57      ENDIF 
    6558 
    66     nn_tide=INT(rday/rdt) 
     59      IF ( kt == nit000 ) ALLOCATE(amp_pot(jpi,jpj,nb_harmo)) 
     60      IF ( kt == nit000 ) ALLOCATE(phi_pot(jpi,jpj,nb_harmo)) 
     61      IF ( kt == nit000 ) ALLOCATE(pot_astro(jpi,jpj)) 
    6762 
    68     CALL tide_init_Wave 
     63      amp_pot(:,:,:) = 0.e0 
     64      phi_pot(:,:,:) = 0.e0 
     65      pot_astro(:,:) = 0.e0 
    6966 
    70     REWIND ( numnam ) 
    71     READ   ( numnam, nam_tide ) 
    72  
    73     IF(lwp) THEN 
    74        WRITE(numout,*) 
    75        WRITE(numout,*) 'sbc_tide : Initialization of the tidal components' 
    76        WRITE(numout,*) '~~~~~~~ ' 
     67      IF ( ln_tide_pot ) CALL tide_init_potential 
     68      ! 
    7769    ENDIF 
    78  
    79     IF(lwp) THEN 
    80        WRITE(numout,*) '        Namelist nam_tide' 
    81        WRITE(numout,*) '        Apply astronomical potential : ln_tide_pot =', ln_tide_pot 
    82        WRITE(numout,*) '        nb_harmo    = ', nb_harmo 
    83        CALL flush(numout) 
    84     ENDIF 
    85  
    86     ALLOCATE(ntide     (nb_harmo)) 
    87     DO jk=1,nb_harmo 
    88        DO ji=1,jpmax_harmo 
    89           IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN 
    90              ntide(jk) = ji 
    91              EXIT 
    92           END IF 
    93        END DO 
    94     END DO 
    95     ALLOCATE(omega_tide(nb_harmo)) 
    96     ALLOCATE(v0tide    (nb_harmo)) 
    97     ALLOCATE(utide     (nb_harmo)) 
    98     ALLOCATE(ftide     (nb_harmo)) 
    99     ALLOCATE(amp_pot(jpi,jpj,nb_harmo)) 
    100     ALLOCATE(phi_pot(jpi,jpj,nb_harmo)) 
    101     ALLOCATE(pot_astro(jpi,jpj)) 
    102     ENDIF 
    103  
    104     IF ( MOD( kt - 1, nn_tide ) == 0 ) THEN 
    105       kt_tide = kt 
    106       CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo) 
    107     ENDIF 
    108  
    109     amp_pot(:,:,:) = 0.e0 
    110     phi_pot(:,:,:) = 0.e0 
    111     pot_astro(:,:) = 0.e0 
    112  
    113     IF (ln_tide_pot          ) CALL tide_init_potential 
    11470 
    11571  END SUBROUTINE sbc_tide 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

    r3294 r3651  
    1313  USE sbctide 
    1414  USE dynspg_oce 
     15  USE tideini, ONLY: ln_tide_ramp, rdttideramp 
    1516 
    1617  IMPLICIT NONE 
    … …  
    3334    INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index 
    3435    INTEGER  :: ji,jj,jk 
     36    REAL (wp) :: zramp 
    3537    REAL (wp), DIMENSION(nb_harmo) :: zwt  
    36     !............................................................................... 
    37     ! Potentiel astronomique 
    3838    !............................................................................... 
    3939 
    4040    pot_astro(:,:)=0.e0 
     41    zramp = 1.e0 
    4142 
    4243    IF (lk_dynspg_ts) THEN 
    4344       zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp))) 
     45       IF (ln_tide_ramp) THEN 
     46          zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.) 
     47       ENDIF 
    4448    ELSE 
    4549       zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt 
     50       IF (ln_tide_ramp) THEN 
     51          zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.)  
     52       ENDIF   
    4653    ENDIF 
    4754 
    … …  
    4956       do ji=1,jpi 
    5057          do jj=1,jpj 
    51              pot_astro(ji,jj)=pot_astro(ji,jj) + (amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))       
     58             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))       
    5259          enddo 
    5360       enddo 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r3625 r3651  
    4949   USE ice_domain_size, only: nx_global, ny_global 
    5050#endif 
     51   USE tideini         ! tidal components initialization   (tide_ini routine) 
    5152   USE obcini          ! open boundary cond. initialization (obc_ini routine) 
    5253   USE bdyini          ! open boundary cond. initialization (bdy_init routine) 
    5354   USE bdydta          ! open boundary cond. initialization (bdy_dta_init routine) 
    54    USE bdytides        ! open boundary cond. initialization (tide_init routine) 
     55   USE bdytides        ! open boundary cond. initialization (bdytide_init routine) 
    5556   USE istate          ! initial state setting          (istate_init routine) 
    5657   USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine) 
    … …  
    8081   USE mod_ioclient 
    8182#endif 
     83   USE sbctide, ONLY: lk_tide 
     84 
    8285 
    8386   IMPLICIT NONE 
    … …  
    322325 
    323326      IF( lk_obc        )   CALL     obc_init   ! Open boundaries  
     327 
     328                            CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
     329 
     330      IF( lk_tide       )   CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics 
     331 
    324332      IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation 
    325333      IF( lk_bdy        )   CALL     bdy_dta_init   ! Open boundaries initialisation of external data arrays 
    326       IF( lk_bdy        )   CALL     tide_init      ! Open boundaries initialisation of tidal harmonic forcing 
     334      IF( lk_bdy        )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing 
    327335 
    328336                            CALL dyn_nept_init  ! simplified form of Neptune effect 
    329  
    330                             CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
    331337 
    332338      !                                     ! Ocean physics 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/step.F90

    r3634 r3651  
    9595      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9696                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     97      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp ) 
    9798      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    9899      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries 
    … …  
    194195      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 
    195196      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends 
     197      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends 
    196198                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection 
    197199      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes 
    … …  
    226228         & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
    227229      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified) 
     230      IF( lk_bdy           )   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends 
    228231                               CALL dyn_adv( kstp )         ! advection (vector or flux form) 
    229232                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r3634 r3651  
    5555 
    5656   USE bdy_par          ! for lk_bdy 
     57   USE bdy_oce          ! for dmp logical 
    5758   USE bdydta           ! open boundary condition data     (bdy_dta routine) 
     59   USE bdytra           ! bdy cond. for tracers            (bdy_tra routine) 
     60   USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine) 
    5861 
    5962   USE sshwzv           ! vertical velocity and ssh        (ssh_wzv routine) 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_MERCATOR_CLUSTER

    r3608 r3651  
    44#PBS -e sette.$PBS_JOBID.err 
    55#PBS -o sette.$PBS_JOBID.out 
    6 #PBS -l nodes=2:ppn=8 
    7 #PBS -q multi 
     6#PBS -l nodes=NODES:ppn=NPROCSNODE 
     7#PBS -q QUEUE 
    88#PBS -l walltime=03:00:00 
    9 #PBS -l mem=24gb 
     9#PBS -l mem=MEMgb 
    1010 
    1111# 
    … …  
    2222# Local settings for machine IBM Power6 (VARGAS at IDRIS France) 
    2323# 
    24 export MPIRUN="mpiexec -n $OCEANCORES" 
     24#cbr export MPIRUN="mpiexec -n $OCEANCORES" 
     25export MPIRUN="mpirun -np $OCEANCORES" 
    2526 
    2627# 
    … …  
    5051# 
    5152cd ${EXE_DIR} 
    52  
    5353  echo Running on host `hostname` 
    5454  echo Time is `date` 
    … …  
    5757#  Run the parallel MPI executable  
    5858# 
     59  echo "Running time ${MPIRUN} ./opa" 
     60# 
    5961  if [ MPI_FLAG == "yes" ]; then 
    60   echo "Running time ${MPIRUN} ./opa" 
    61      time ${MPIRUN} ./opa 
     62#cbr     time ${MPIRUN} ./opa 
     63     mpirun -np $OCEANCORES ./opa 
     64#cbr     mpirun -np $OCEANCORES ./opa 
    6265  else 
    63   echo "Running time ./opa" 
    6466     time ./opa 
    6567  fi 

  • branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/SETTE/iodef_sette.xml

    r3635 r3651  
    140140   <field id="utau"         description="Wind Stress along i-axis"                    unit="N/m2" axis_ref="none" /> 
    141141   <field id="uoce"         description="ocean current along i-axis"                  unit="m/s"                  /> 
    142    <field id="uocetr_eff"   description="Effective ocean transport along i-axis"      unit="m3/s"                  /> 
     142   <field id="uocetr_eff"     description="Effective ocean transport along i-axis"    unit="m3/s"                 /> 
    143143   <!-- uoce_eiv: available with key_traldf_eiv and key_diaeiv --> 
    144144   <field id="uoce_eiv"     description="EIV ocean current along i-axis"              unit="m/s"                  /> 
    … …  
    158158   <field id="vtau"         description="Wind Stress along j-axis"                    unit="N/m2" axis_ref="none" /> 
    159159   <field id="voce"         description="ocean current along j-axis"                  unit="m/s"                  /> 
    160    <field id="vocetr_eff"   description="Effective ocean transport along j-axis"      unit="m3/s"                  /> 
     160   <field id="vocetr_eff"     description="Effective ocean tra,sport along j-axis"    unit="m3/s"                 /> 
    161161   <!-- voce_eiv: available with key_traldf_eiv and key_diaeiv --> 
    162162   <field id="voce_eiv"     description="EIV ocean current along j-axis"              unit="m/s"                  /> 
    … …  
    175175      <group id="grid_W"  axis_ref="depthw" grid_ref="grid_W"> 
    176176   <field id="woce"         description="ocean vertical velocity"                     unit="m/s"                  /> 
    177    <field id="wocetr_eff"   description="effective ocean vertical transport"          unit="m3/s"                 /> 
     177   <field id="wocetr_eff"     description="effective ocean vertical transport"          unit="m3/s"                 /> 
    178178   <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv --> 
    179179   <field id="woce_eiv"     description="EIV ocean vertical velocity"                 unit="m/s"                  /> 
    … …  
    209209   <field id="fram_trans" description="Sea Ice Mass Transport Through Fram Strait" unit="kg/s" /> 
    210210      </group> 
    211           
     211 
     212      <!-- variables available with key_float --> 
     213      <group id="floatvar" axis_ref="nfloat" grid_ref="scalarpoint" zoom_ref="1point"> 
     214        <field id="traj_lon"   description="floats longitude"   unit="deg"  operation="inst(X)" /> 
     215        <field id="traj_lat"   description="floats latitude"    unit="deg"   /> 
     216        <field id="traj_dep"   description="floats depth"       unit="m"     /> 
     217        <field id="traj_temp"  description="floats temperature" unit="degC"  /> 
     218        <field id="traj_salt"  description="floats salinity"    unit="psu"   /> 
     219        <field id="traj_dens"  description="floats density"     unit="kg/m3" /> 
     220        <field id="traj_group" description="floats group"       unit="none"  /> 
     221      </group> 
     222      
    212223      <!-- ptrc on T grid --> 
    213224 
    214      <group id="ptrc_T" axis_ref="deptht" grid_ref="grid_T">   
     225     <group id="ptrc_T" axis_ref="deptht" grid_ref="grid_T"> 
    215226       <field id="DIC"      description="Dissolved inorganic Concentration"        unit="mmol/m3" /> 
    216227       <field id="Alkalini" description="Total Alkalinity Concentration"           unit="mmol/m3"    /> 
    … …  
    237248       <field id="NO3"      description="Nitrates Concentration"                   unit="mmol/m3" /> 
    238249       <field id="NH4"      description="Ammonium Concentration"                   unit="mmol/m3" /> 
     250       <field id="DET"      description="Detritus"                                      unit="mmole-N/m3" /> 
     251       <field id="ZOO"      description="Zooplankton Concentration"                     unit="mmole-N/m3" /> 
     252       <field id="PHY"      description="Phytoplankton Concentration"                   unit="mmole-N/m3" /> 
     253       <field id="NO3"      description="Nitrate Concentration"                         unit="mmole-N/m3" /> 
     254       <field id="NH4"      description="Ammonium Concentration"                        unit="mmole-N/m3" /> 
     255       <field id="DOM"      description="Dissolved Organic Matter"                      unit="mmole-N/m3" /> 
    239256     </group> 
    240257 
    … …  
    246263       <field id="CO3sat"      description="CO3 saturation"                           unit="mol/m3"     axis_ref="deptht" /> 
    247264       <field id="PAR"         description="Photosynthetically Available Radiation"   unit="W/m2"       axis_ref="deptht" /> 
    248        <field id="PPPHY"       description="Primary production of nanophyto"          unit="molC/m3/s"   axis_ref="deptht" /> 
    249        <field id="PPPHY2"      description="Primary production of diatoms"            unit="molC/m3/s"   axis_ref="deptht" /> 
    250        <field id="PPNEWN"      description="New Primary production of nanophyto"      unit="molC/m3/s"   axis_ref="deptht" /> 
    251        <field id="PPNEWD"      description="New Primary production of diatoms"        unit="molC/m3/s"   axis_ref="deptht" /> 
    252        <field id="PBSi"        description="Primary production of Si diatoms"         unit="molSi/m3/s"   axis_ref="deptht" /> 
    253        <field id="PFeN"        description="Primary production of nano iron"          unit="molFe/m3/s"   axis_ref="deptht" /> 
    254        <field id="PFeD"        description="Primary production of diatoms iron"       unit="molFe/m3/s"   axis_ref="deptht" /> 
    255        <field id="xfracal"     description="Calcifying fraction"                      unit="-"          axis_ref="deptht" /> 
    256        <field id="PCAL"        description="Calcite production"                       unit="molC/m3/s"   axis_ref="deptht" /> 
    257        <field id="DCAL"        description="Calcite dissolution"                      unit="molC/m3/s"   axis_ref="deptht" /> 
    258        <field id="GRAZ1"       description="Grazing by microzooplankton"              unit="molC/m3/s"   axis_ref="deptht" /> 
    259        <field id="GRAZ2"       description="Grazing by mesozooplankton"              unit="molC/m3/s"   axis_ref="deptht" /> 
    260        <field id="REMIN"       description="Oxic remineralization of OM"             unit="molC/m3/s"   axis_ref="deptht" /> 
    261        <field id="DENIT"       description="Anoxic remineralization of OM"           unit="molC/m3/s"   axis_ref="deptht" /> 
    262        <field id="Nfix"        description="Nitrogen fixation"                         unit="molN/m3/s"   axis_ref="deptht" /> 
     265       <field id="PPPHY"       description="Primary production of nanophyto"          unit="mol/m3/s"   axis_ref="deptht" /> 
     266       <field id="PPPHY2"      description="Primary production of diatoms"            unit="mol/m3/s"   axis_ref="deptht" /> 
     267       <field id="PPNEWN"      description="New Primary production of nanophyto"      unit="mol/m3/s"   axis_ref="deptht" /> 
     268       <field id="PPNEWD"      description="New Primary production of diatoms"        unit="mol/m3/s"   axis_ref="deptht" /> 
     269       <field id="PBSi"        description="Primary production of Si diatoms"         unit="mol/m3/s"   axis_ref="deptht" /> 
     270       <field id="PFeN"        description="Primary production of nano iron"          unit="mol/m3/s"   axis_ref="deptht" /> 
     271       <field id="PFeD"        description="Primary production of diatoms iron"       unit="mol/m3/s"   axis_ref="deptht" /> 
     272       <field id="PCAL"        description="Calcite production"                       unit="mol/m3/s"   axis_ref="deptht" /> 
     273       <field id="DCAL"        description="Calcite dissolution"                      unit="mol/m3/s"   axis_ref="deptht" /> 
     274       <field id="GRAZ"        description="Grazing by zooplankton"                   unit="mol/m3/s"   axis_ref="deptht" /> 
    263275       <field id="Mumax"       description="Maximum growth rate"                      unit="s-1"        axis_ref="deptht" /> 
    264276       <field id="MuN"         description="Realized growth rate for nanophyto"       unit="s-1"        axis_ref="deptht" /> 
    265277       <field id="MuD"         description="Realized growth rate for diatomes"        unit="s-1"        axis_ref="deptht" /> 
     278       <field id="MuNlight"    description="Light limited growth rate for nanophyto"  unit="s-1"        axis_ref="deptht" /> 
     279       <field id="MuDlight"    description="Light limited growth rate for diatomes"   unit="s-1"        axis_ref="deptht" /> 
    266280       <field id="LNnut"       description="Nutrient limitation term in Nanophyto"    unit="-"          axis_ref="deptht" /> 
    267281       <field id="LDnut"       description="Nutrient limitation term in Diatoms"      unit="-"          axis_ref="deptht" /> 
    … …  
    270284       <field id="LNlight"     description="Light limitation term in Nanophyto"       unit="-"          axis_ref="deptht" /> 
    271285       <field id="LDlight"     description="Light limitation term in Diatoms"         unit="-"          axis_ref="deptht" /> 
    272        <field id="Fe2"        description="Iron II concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    273        <field id="Fe3"        description="Iron III concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    274        <field id="FeL1"        description="Complexed Iron concentration with L1"       unit="nmol/L"   axis_ref="deptht" /> 
    275        <field id="FeL2"        description="Complexed Iron concentration with L2"       unit="nmol/L"   axis_ref="deptht" /> 
    276        <field id="FeP"        description="Precipitated Iron III"                      unit="nmol/L"   axis_ref="deptht" /> 
    277        <field id="TL1"        description="Total L1 concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    278        <field id="TL2"        description="Total L2 concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    279        <field id="pdust"        description="dust concentration"                      unit="g/L"    /> 
    280        <field id="Totlig"        description="Total ligand concentation"                unit="nmol/L"   axis_ref="deptht" /> 
    281        <field id="Biron"        description="Bioavailable iron"                       unit="nmol/L"   axis_ref="deptht" /> 
    282        <field id="Sdenit"        description="Nitrate reduction in the sediments"     unit="molN/m2/s"      /> 
    283        <field id="Ironice"     description="Iron input/uptake due to sea ice"             unit="molFe/m2/s"      /> 
    284        <field id="HYDR"       description="Iron input from hydrothemal vents"          unit="molFe/m2/s"      /> 
     286       <field id="Nfix"        description="Nitrogen fixation at surface"             unit="mol/m2/s"                     /> 
    285287       <field id="EPC100"      description="Export of carbon particles at 100 m"      unit="mol/m2/s"                     /> 
    286288       <field id="EPFE100"     description="Export of biogenic iron at 100 m"         unit="mol/m2/s"                     /> 
    … …  
    293295       <field id="Dpo2"        description="Delta O2"                                 unit="uatm"                         /> 
    294296       <field id="Heup"        description="Euphotic layer depth"                     unit="m"                            /> 
    295        <field id="Irondep"     description="Iron deposition from dust"                 unit="mol/m2/s"                     /> 
    296        <field id="Ironsed"     description="Iron deposition from sediment"             unit="mol/m2/s"  axis_ref="deptht"  /> 
     297       <field id="Irondep"     description="Iron deposition"                          unit="mol/m2/s"                     /> 
     298       <field id="FNO3PHY"     description="FNO3PHY"                                  unit="-"          axis_ref="deptht" /> 
     299       <field id="FNH4PHY"     description="FNH4PHY"                                  unit="-"          axis_ref="deptht" /> 
     300       <field id="FNH4NO3"     description="FNH4NO3"                                  unit="-"          axis_ref="deptht" /> 
     301       <field id="TNO3PHY"     description="TNO3PHY"                                  unit="-"  /> 
     302       <field id="TNH4PHY"     description="TNH4PHY"                                  unit="-"  /> 
     303       <field id="TPHYDOM"     description="TPHYDOM"                                  unit="-"  /> 
     304       <field id="TPHYNH4"     description="TPHYNH4"                                  unit="-"  /> 
     305       <field id="TPHYZOO"     description="TPHYZOO"                                  unit="-"  /> 
     306       <field id="TPHYDET"     description="TPHYDET"                                  unit="-"  /> 
     307       <field id="TDETZOO"     description="TDETZOO"                                  unit="-"  /> 
     308       <field id="TDETSED"     description="TDETSED"                                  unit="-"  /> 
     309       <field id="TZOODET"     description="TZOODET"                                  unit="-"  /> 
     310       <field id="TZOOBOD"     description="TZOOBOD"                                  unit="-"  /> 
     311       <field id="TZOONH4"     description="TZOONH4"                                  unit="-"  /> 
     312       <field id="TZOODOM"     description="TZOODOM"                                  unit="-"  /> 
     313       <field id="TNH4NO3"     description="TNH4NO3"                                  unit="-"  /> 
     314       <field id="TDOMNH4"     description="TDOMNH4"                                  unit="-"  /> 
     315       <field id="TDETNH4"     description="TDETNH4"                                  unit="-"  /> 
     316       <field id="TPHYTOT"     description="TPHYTOT"                                  unit="-"  /> 
     317       <field id="TZOOTOT"     description="TZOOTOT"                                  unit="-"  /> 
     318       <field id="TDETDOM"     description="TDETDOM"                                  unit="-"  /> 
     319       <field id="SEDPOC"      description="SEDPOC"                                   unit="-"  /> 
    297320     </group> 
    298321        
    … …  
    374397      <axis id="depthv" description="Vertical V levels" unit="m" positive=".false." /> 
    375398      <axis id="depthw" description="Vertical W levels" unit="m" positive=".false." /> 
     399      <axis id="nfloat" description="Number of float"   unit="no unit" positive=".false." /> 
    376400      <axis id="none" description="axe non defini" unit="none" size="1" /> 
    377401    </axis_definition>  
    … …  
    551575  </context> 
    552576   
    553  
    554577  <context id="1_nemo"> 
    555578     
    … …  
    690713   <field id="utau"         description="Wind Stress along i-axis"                    unit="N/m2" axis_ref="none" /> 
    691714   <field id="uoce"         description="ocean current along i-axis"                  unit="m/s"                  /> 
    692    <field id="uocetr_eff"   description="Effective ocean transport along i-axis"      unit="m3/s"                  /> 
     715   <field id="uocetr_eff"     description="Effective ocean transport along i-axis"    unit="m3/s"                 /> 
    693716   <!-- uoce_eiv: available with key_traldf_eiv and key_diaeiv --> 
    694717   <field id="uoce_eiv"     description="EIV ocean current along i-axis"              unit="m/s"                  /> 
    … …  
    708731   <field id="vtau"         description="Wind Stress along j-axis"                    unit="N/m2" axis_ref="none" /> 
    709732   <field id="voce"         description="ocean current along j-axis"                  unit="m/s"                  /> 
    710    <field id="vocetr_eff"   description="Effective ocean transport along j-axis"      unit="m3/s"                  /> 
     733   <field id="vocetr_eff"     description="Effective ocean tra,sport along j-axis"    unit="m3/s"                 /> 
    711734   <!-- voce_eiv: available with key_traldf_eiv and key_diaeiv --> 
    712735   <field id="voce_eiv"     description="EIV ocean current along j-axis"              unit="m/s"                  /> 
    … …  
    725748      <group id="grid_W"  axis_ref="depthw" grid_ref="grid_W"> 
    726749   <field id="woce"         description="ocean vertical velocity"                     unit="m/s"                  /> 
    727    <field id="wocetr_eff"   description="effective ocean vertical transport"          unit="m3/s"                 /> 
     750   <field id="wocetr_eff"     description="effective ocean vertical transport"          unit="m3/s"                 /> 
    728751   <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv --> 
    729752   <field id="woce_eiv"     description="EIV ocean vertical velocity"                 unit="m/s"                  /> 
    … …  
    759782   <field id="fram_trans" description="Sea Ice Mass Transport Through Fram Strait" unit="kg/s" /> 
    760783      </group> 
    761           
     784 
     785      <!-- variables available with key_float --> 
     786      <group id="floatvar" axis_ref="nfloat" grid_ref="scalarpoint" zoom_ref="1point"> 
     787        <field id="traj_lon"   description="floats longitude"   unit="deg"  operation="inst(X)" /> 
     788        <field id="traj_lat"   description="floats latitude"    unit="deg"   /> 
     789        <field id="traj_dep"   description="floats depth"       unit="m"     /> 
     790        <field id="traj_temp"  description="floats temperature" unit="degC"  /> 
     791        <field id="traj_salt"  description="floats salinity"    unit="psu"   /> 
     792        <field id="traj_dens"  description="floats density"     unit="kg/m3" /> 
     793        <field id="traj_group" description="floats group"       unit="none"  /> 
     794      </group> 
     795      
    762796      <!-- ptrc on T grid --> 
    763797 
    764      <group id="ptrc_T" axis_ref="deptht" grid_ref="grid_T">   
     798     <group id="ptrc_T" axis_ref="deptht" grid_ref="grid_T"> 
    765799       <field id="DIC"      description="Dissolved inorganic Concentration"        unit="mmol/m3" /> 
    766800       <field id="Alkalini" description="Total Alkalinity Concentration"           unit="mmol/m3"    /> 
    … …  
    787821       <field id="NO3"      description="Nitrates Concentration"                   unit="mmol/m3" /> 
    788822       <field id="NH4"      description="Ammonium Concentration"                   unit="mmol/m3" /> 
     823       <field id="DET"      description="Detritus"                                      unit="mmole-N/m3" /> 
     824       <field id="ZOO"      description="Zooplankton Concentration"                     unit="mmole-N/m3" /> 
     825       <field id="PHY"      description="Phytoplankton Concentration"                   unit="mmole-N/m3" /> 
     826       <field id="NO3"      description="Nitrate Concentration"                         unit="mmole-N/m3" /> 
     827       <field id="NH4"      description="Ammonium Concentration"                        unit="mmole-N/m3" /> 
     828       <field id="DOM"      description="Dissolved Organic Matter"                      unit="mmole-N/m3" /> 
    789829     </group> 
    790830 
    … …  
    796836       <field id="CO3sat"      description="CO3 saturation"                           unit="mol/m3"     axis_ref="deptht" /> 
    797837       <field id="PAR"         description="Photosynthetically Available Radiation"   unit="W/m2"       axis_ref="deptht" /> 
    798        <field id="PPPHY"       description="Primary production of nanophyto"          unit="molC/m3/s"   axis_ref="deptht" /> 
    799        <field id="PPPHY2"      description="Primary production of diatoms"            unit="molC/m3/s"   axis_ref="deptht" /> 
    800        <field id="PPNEWN"      description="New Primary production of nanophyto"      unit="molC/m3/s"   axis_ref="deptht" /> 
    801        <field id="PPNEWD"      description="New Primary production of diatoms"        unit="molC/m3/s"   axis_ref="deptht" /> 
    802        <field id="PBSi"        description="Primary production of Si diatoms"         unit="molSi/m3/s"   axis_ref="deptht" /> 
    803        <field id="PFeN"        description="Primary production of nano iron"          unit="molFe/m3/s"   axis_ref="deptht" /> 
    804        <field id="PFeD"        description="Primary production of diatoms iron"       unit="molFe/m3/s"   axis_ref="deptht" /> 
    805        <field id="xfracal"     description="Calcifying fraction"                      unit="-"          axis_ref="deptht" /> 
    806        <field id="PCAL"        description="Calcite production"                       unit="molC/m3/s"   axis_ref="deptht" /> 
    807        <field id="DCAL"        description="Calcite dissolution"                      unit="molC/m3/s"   axis_ref="deptht" /> 
    808        <field id="GRAZ1"       description="Grazing by microzooplankton"              unit="molC/m3/s"   axis_ref="deptht" /> 
    809        <field id="GRAZ2"       description="Grazing by mesozooplankton"              unit="molC/m3/s"   axis_ref="deptht" /> 
    810        <field id="REMIN"       description="Oxic remineralization of OM"             unit="molC/m3/s"   axis_ref="deptht" /> 
    811        <field id="DENIT"       description="Anoxic remineralization of OM"           unit="molC/m3/s"   axis_ref="deptht" /> 
    812        <field id="Nfix"        description="Nitrogen fixation"                         unit="molN/m3/s"   axis_ref="deptht" /> 
     838       <field id="PPPHY"       description="Primary production of nanophyto"          unit="mol/m3/s"   axis_ref="deptht" /> 
     839       <field id="PPPHY2"      description="Primary production of diatoms"            unit="mol/m3/s"   axis_ref="deptht" /> 
     840       <field id="PPNEWN"      description="New Primary production of nanophyto"      unit="mol/m3/s"   axis_ref="deptht" /> 
     841       <field id="PPNEWD"      description="New Primary production of diatoms"        unit="mol/m3/s"   axis_ref="deptht" /> 
     842       <field id="PBSi"        description="Primary production of Si diatoms"         unit="mol/m3/s"   axis_ref="deptht" /> 
     843       <field id="PFeN"        description="Primary production of nano iron"          unit="mol/m3/s"   axis_ref="deptht" /> 
     844       <field id="PFeD"        description="Primary production of diatoms iron"       unit="mol/m3/s"   axis_ref="deptht" /> 
     845       <field id="PCAL"        description="Calcite production"                       unit="mol/m3/s"   axis_ref="deptht" /> 
     846       <field id="DCAL"        description="Calcite dissolution"                      unit="mol/m3/s"   axis_ref="deptht" /> 
     847       <field id="GRAZ"        description="Grazing by zooplankton"                   unit="mol/m3/s"   axis_ref="deptht" /> 
    813848       <field id="Mumax"       description="Maximum growth rate"                      unit="s-1"        axis_ref="deptht" /> 
    814849       <field id="MuN"         description="Realized growth rate for nanophyto"       unit="s-1"        axis_ref="deptht" /> 
    815850       <field id="MuD"         description="Realized growth rate for diatomes"        unit="s-1"        axis_ref="deptht" /> 
     851       <field id="MuNlight"    description="Light limited growth rate for nanophyto"  unit="s-1"        axis_ref="deptht" /> 
     852       <field id="MuDlight"    description="Light limited growth rate for diatomes"   unit="s-1"        axis_ref="deptht" /> 
    816853       <field id="LNnut"       description="Nutrient limitation term in Nanophyto"    unit="-"          axis_ref="deptht" /> 
    817854       <field id="LDnut"       description="Nutrient limitation term in Diatoms"      unit="-"          axis_ref="deptht" /> 
    … …  
    820857       <field id="LNlight"     description="Light limitation term in Nanophyto"       unit="-"          axis_ref="deptht" /> 
    821858       <field id="LDlight"     description="Light limitation term in Diatoms"         unit="-"          axis_ref="deptht" /> 
    822        <field id="Fe2"        description="Iron II concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    823        <field id="Fe3"        description="Iron III concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    824        <field id="FeL1"        description="Complexed Iron concentration with L1"       unit="nmol/L"   axis_ref="deptht" /> 
    825        <field id="FeL2"        description="Complexed Iron concentration with L2"       unit="nmol/L"   axis_ref="deptht" /> 
    826        <field id="FeP"        description="Precipitated Iron III"                      unit="nmol/L"   axis_ref="deptht" /> 
    827        <field id="TL1"        description="Total L1 concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    828        <field id="TL2"        description="Total L2 concentration"                      unit="nmol/L"   axis_ref="deptht" /> 
    829        <field id="pdust"        description="dust concentration"                      unit="g/L"    /> 
    830        <field id="Totlig"        description="Total ligand concentation"                unit="nmol/L"   axis_ref="deptht" /> 
    831        <field id="Biron"        description="Bioavailable iron"                       unit="nmol/L"   axis_ref="deptht" /> 
    832        <field id="Sdenit"        description="Nitrate reduction in the sediments"     unit="molN/m2/s"      /> 
    833        <field id="Ironice"     description="Iron input/uptake due to sea ice"             unit="molFe/m2/s"      /> 
    834        <field id="HYDR"       description="Iron input from hydrothemal vents"          unit="molFe/m2/s"      /> 
     859       <field id="Nfix"        description="Nitrogen fixation at surface"             unit="mol/m2/s"                     /> 
    835860       <field id="EPC100"      description="Export of carbon particles at 100 m"      unit="mol/m2/s"                     /> 
    836861       <field id="EPFE100"     description="Export of biogenic iron at 100 m"         unit="mol/m2/s"                     /> 
    … …  
    843868       <field id="Dpo2"        description="Delta O2"                                 unit="uatm"                         /> 
    844869       <field id="Heup"        description="Euphotic layer depth"                     unit="m"                            /> 
    845        <field id="Irondep"     description="Iron deposition from dust"                 unit="mol/m2/s"                     /> 
    846        <field id="Ironsed"     description="Iron deposition from sediment"             unit="mol/m2/s"  axis_ref="deptht"  /> 
     870       <field id="Irondep"     description="Iron deposition"                          unit="mol/m2/s"                     /> 
     871       <field id="FNO3PHY"     description="FNO3PHY"                                  unit="-"          axis_ref="deptht" /> 
     872       <field id="FNH4PHY"     description="FNH4PHY"                                  unit="-"          axis_ref="deptht" /> 
     873       <field id="FNH4NO3"     description="FNH4NO3"                                  unit="-"          axis_ref="deptht" /> 
     874       <field id="TNO3PHY"     description="TNO3PHY"                                  unit="-"  /> 
     875       <field id="TNH4PHY"     description="TNH4PHY"                                  unit="-"  /> 
     876       <field id="TPHYDOM"     description="TPHYDOM"                                  unit="-"  /> 
     877       <field id="TPHYNH4"     description="TPHYNH4"                                  unit="-"  /> 
     878       <field id="TPHYZOO"     description="TPHYZOO"                                  unit="-"  /> 
     879       <field id="TPHYDET"     description="TPHYDET"                                  unit="-"  /> 
     880       <field id="TDETZOO"     description="TDETZOO"                                  unit="-"  /> 
     881       <field id="TDETSED"     description="TDETSED"                                  unit="-"  /> 
     882       <field id="TZOODET"     description="TZOODET"                                  unit="-"  /> 
     883       <field id="TZOOBOD"     description="TZOOBOD"                                  unit="-"  /> 
     884       <field id="TZOONH4"     description="TZOONH4"                                  unit="-"  /> 
     885       <field id="TZOODOM"     description="TZOODOM"                                  unit="-"  /> 
     886       <field id="TNH4NO3"     description="TNH4NO3"                                  unit="-"  /> 
     887       <field id="TDOMNH4"     description="TDOMNH4"                                  unit="-"  /> 
     888       <field id="TDETNH4"     description="TDETNH4"                                  unit="-"  /> 
     889       <field id="TPHYTOT"     description="TPHYTOT"                                  unit="-"  /> 
     890       <field id="TZOOTOT"     description="TZOOTOT"                                  unit="-"  /> 
     891       <field id="TDETDOM"     description="TDETDOM"                                  unit="-"  /> 
     892       <field id="SEDPOC"      description="SEDPOC"                                   unit="-"  /> 
    847893     </group>