Changeset 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY

Timestamp:

2020-05-14T21:46:00+02:00 (4 years ago)

Author:

smueller

Message:

Synchronizing with ^{/NEMO/trunk@12925 (ticket #2170)}

Location:

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

Files:

• . (modified) (1 prop)
• src/OCE/BDY/bdy_oce.F90 (modified) (1 diff)
• src/OCE/BDY/bdydta.F90 (modified) (26 diffs)
• src/OCE/BDY/bdydyn.F90 (modified) (6 diffs)
• src/OCE/BDY/bdydyn3d.F90 (modified) (27 diffs)
• src/OCE/BDY/bdyice.F90 (modified) (1 diff)
• src/OCE/BDY/bdyini.F90 (modified) (7 diffs)
• src/OCE/BDY/bdylib.F90 (modified) (12 diffs)
• src/OCE/BDY/bdytides.F90 (modified) (19 diffs)
• src/OCE/BDY/bdytra.F90 (modified) (10 diffs)
• src/OCE/BDY/bdyvol.F90 (modified) (3 diffs)

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

@@ -6 +6 @@
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
+
+# SETTE
+^/utils/CI/sette@HEAD         sette

@@ -6 +6 @@
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
+
+# SETTE
+^/utils/CI/sette@HEAD         sette

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdy_oce.F90

@@ -141 +141 @@
    LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lsend_bdyext   !: mark needed communication for given boundary, grid and neighbour
    LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lrecv_bdyext   !:  when searching towards the exterior of the computational domain
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdydta.F90

@@ -23 +23 @@
    USE phycst         ! physical constants
    USE sbcapr         ! atmospheric pressure forcing
-   USE sbctide        ! Tidal forcing or not
+   USE tide_mod, ONLY: ln_tide ! tidal forcing
    USE bdy_oce        ! ocean open boundary conditions  
    USE bdytides       ! tidal forcing at boundaries
@@ -68 +68 @@
 !$AGRIF_END_DO_NOT_TREAT
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -75 +77 @@
 CONTAINS
 
-   SUBROUTINE bdy_dta( kt, kit, kt_offset )
+   SUBROUTINE bdy_dta( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta  ***
@@ -85 +87 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
-      INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option)
-      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit
-      !                                               ! is present then units = subcycle timesteps.
-      !                                               ! kt_offset = 0 => get data at "now" time level
-      !                                               ! kt_offset = -1 => get data at "before" time level
-      !                                               ! kt_offset = +1 => get data at "after" time level
-      !                                               ! etc.
+      INTEGER, INTENT(in)           ::   Kmm          ! ocean time level index
       !
       INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices
       INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers
-      INTEGER,   DIMENSION(jpbgrd)     ::   ilen1 
-      INTEGER,   DIMENSION(:), POINTER ::   nblen, nblenrim  ! short cuts
       TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut
       TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias
@@ -105 +99 @@
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN
+      IF( kt == nit000 ) THEN
 
          ! Calculate depth-mean currents
@@ -112 +106 @@
          DO jbdy = 1, nb_bdy
             !
-            nblen    => idx_bdy(jbdy)%nblen
-            nblenrim => idx_bdy(jbdy)%nblenrim
-            !
             IF( nn_dyn2d_dta(jbdy) == 0 ) THEN 
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_ssh ) THEN 
                   igrd = 1
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is allocated and used only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                     dta_bdy(jbdy)%ssh(ib) = ssh(ii,ij,Kmm) * tmask(ii,ij,1)         
                   END DO
                ENDIF
-               IF( dta_bdy(jbdy)%lneed_dyn2d) THEN 
+               IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
                   igrd = 2
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)      ! u2d is used either over the whole bdy or only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
+                     dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1)         
                   END DO
+               ENDIF
+               IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
                   igrd = 3
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)      ! v2d is used either over the whole bdy or only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
+                     dta_bdy(jbdy)%v2d(ib) = vv_b(ii,ij,Kmm) * vmask(ii,ij,1)         
                   END DO
                ENDIF
@@ -142 +134 @@
             !
             IF( nn_dyn3d_dta(jbdy) == 0 ) THEN 
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN 
                   igrd = 2 
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      DO ik = 1, jpkm1
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
+                        dta_bdy(jbdy)%u3d(ib,ik) =  ( uu(ii,ij,ik,Kmm) - uu_b(ii,ij,Kmm) ) * umask(ii,ij,ik)         
                      END DO
                   END DO
                   igrd = 3 
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      DO ik = 1, jpkm1
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vv(ii,ij,ik,Kmm) - vv_b(ii,ij,Kmm) ) * vmask(ii,ij,ik)         
                      END DO
                   END DO
@@ -164 +155 @@
 
             IF( nn_tra_dta(jbdy) == 0 ) THEN 
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_tra ) THEN
                   igrd = 1 
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      DO ik = 1, jpkm1
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_bdytem) * tmask(ii,ij,ik)         
-                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_bdysal) * tmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%tem(ib,ik) = ts(ii,ij,ik,jp_tem,Kmm) * tmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%sal(ib,ik) = ts(ii,ij,ik,jp_sal,Kmm) * tmask(ii,ij,ik)         
                      END DO
                   END DO
@@ -180 +170 @@
 #if defined key_si3
             IF( nn_ice_dta(jbdy) == 0 ) THEN    ! set ice to initial values
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_ice ) THEN
                   igrd = 1   
                   DO jl = 1, jpl
-                     DO ib = 1, ilen1(igrd)
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
@@ -216 +205 @@
          ! read/update all bdy data
          ! ------------------------
-         CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset )
-
+         ! BDY: use pt_offset=0.5 as applied at the end of the step and fldread is referenced at the middle of the step
+         CALL fld_read( kt, 1, bf_alias, pt_offset = 0.5_wp, Kmm = Kmm )
          ! apply some corrections in some specific cases...
          ! --------------------------------------------------
          !
          ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s)
-         IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN   ! runoff and we read u/v2d
+         IF( cn_tra(jbdy) == 'runoff' ) THEN   ! runoff
             !
-            igrd = 2                      ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
-            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-               dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
-            END DO
-            igrd = 3                      ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
-            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-               dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
-            END DO
+            IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
+               igrd = 2                         ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
+               DO ib = 1, SIZE(dta_alias%u2d)   ! u2d is used either over the whole bdy or only on the rim
+                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+               END DO
+            ENDIF
+            IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
+               igrd = 3                         ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
+               DO ib = 1, SIZE(dta_alias%v2d)   ! v2d is used either over the whole bdy or only on the rim
+                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+               END DO
+            ENDIF
          ENDIF
 
          ! tidal harmonic forcing ONLY: initialise arrays
          IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d 
-            IF( dta_alias%lneed_ssh   ) dta_alias%ssh(:) = 0._wp
-            IF( dta_alias%lneed_dyn2d ) dta_alias%u2d(:) = 0._wp
-            IF( dta_alias%lneed_dyn2d ) dta_alias%v2d(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp
          ENDIF
 
@@ -249 +242 @@
             !
             igrd = 2                       ! zonal velocity
-            dta_alias%u2d(:) = 0._wp       ! compute barotrope zonal velocity and put it in u2d
             DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                ii   = idx_bdy(jbdy)%nbi(ib,igrd)
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               dta_alias%u2d(ib) = 0._wp   ! compute barotrope zonal velocity and put it in u2d
                DO ik = 1, jpkm1
-                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
                END DO
-               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu_n(ii,ij)
+               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm)
                DO ik = 1, jpkm1            ! compute barocline zonal velocity and put it in u3d
                   dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib)
@@ -262 +255 @@
             END DO
             igrd = 3                       ! meridional velocity
-            dta_alias%v2d(:) = 0._wp       ! compute barotrope meridional velocity and put it in v2d
             DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                ii   = idx_bdy(jbdy)%nbi(ib,igrd)
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               dta_alias%v2d(ib) = 0._wp   ! compute barotrope meridional velocity and put it in v2d
                DO ik = 1, jpkm1
-                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
                END DO
-               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv_n(ii,ij)
+               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm)
                DO ik = 1, jpkm1            ! compute barocline meridional velocity and put it in v3d
                   dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib)
@@ -275 +268 @@
             END DO
          ENDIF   ! ltotvel
-
-         ! update tidal harmonic forcing
-         IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
-            CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   & 
-               &                 kit = kit, kt_offset = kt_offset )
-         ENDIF
 
          !  atm surface pressure : add inverted barometer effect to ssh if it was read
@@ -293 +280 @@
 
 #if defined key_si3
-         IF( dta_alias%lneed_ice ) THEN
+         IF( dta_alias%lneed_ice .AND. idx_bdy(jbdy)%nblen(1) > 0 ) THEN
             ! fill temperature and salinity arrays
             IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem (jbdy)
@@ -302 +289 @@
                &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i )
             IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy)
+            
+            ! if T_i is read and not T_su, set T_su = T_i
+            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_i)%fnow(:,1,:)
+            ! if T_s is read and not T_su, set T_su = T_s
+            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_s)%fnow(:,1,:)
+            ! if T_i is read and not T_s, set T_s = T_i
+            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdyt_i)%fnow(:,1,:)
+            ! if T_su is read and not T_s, set T_s = T_su
+            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdytsu)%fnow(:,1,:)
             ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
             IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
                &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdytsu)%fnow(:,1,:) + 271.15 )
-            ! if T_su is read and not T_s, set T_s = T_su
-            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) &
-               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdytsu)%fnow(:,1,:)
-            ! if T_s is read and not T_su, set T_su = T_s
-            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
-               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_s)%fnow(:,1,:)
             ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
             IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
@@ -341 +335 @@
             DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop
                IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
-                  nblen => idx_bdy(jbdy)%nblen
-                  nblenrim => idx_bdy(jbdy)%nblenrim
-                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
-                     IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
-                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
-                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
-                  ENDIF
-               END DO
-            ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
-               !
-               CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset )
-            ENDIF
+                  IF( ASSOCIATED(dta_bdy(jbdy)%ssh) ) dta_bdy_s(jbdy)%ssh(:) = dta_bdy(jbdy)%ssh(:)
+                  IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) dta_bdy_s(jbdy)%u2d(:) = dta_bdy(jbdy)%u2d(:)
+                  IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) dta_bdy_s(jbdy)%v2d(:) = dta_bdy(jbdy)%v2d(:)
+               ENDIF
+            END DO
+         ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
+            !
+            CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp )
          ENDIF
-         !
-         IF( ln_timing )   CALL timing_stop('bdy_dta')
-         !
-      END SUBROUTINE bdy_dta
-
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('bdy_dta')
+      !
+   END SUBROUTINE bdy_dta
+   
 
    SUBROUTINE bdy_dta_init
@@ -373 +364 @@
       INTEGER ::   ierror, ios     ! 
       !
+      INTEGER ::   nbdy_rdstart, nbdy_loc
+      CHARACTER(LEN=50)                      ::   cerrmsg       ! error string
       CHARACTER(len=3)                       ::   cl3           ! 
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
@@ -388 +381 @@
       LOGICAL                                ::   llneed        !
       LOGICAL                                ::   llread        !
+      LOGICAL                                ::   llfullbdy     !
       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
@@ -415 +409 @@
       ! Read namelists
       ! --------------
-      REWIND(numnam_cfg)
+      nbdy_rdstart = 1
       DO jbdy = 1, nb_bdy
 
@@ -421 +415 @@
          WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy
 
-         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind 
-         REWIND(numnam_ref)
+         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we read from the beginning
          READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist' )
@@ -431 +424 @@
             & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   &
             & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN
-            ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another
-            READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
+            !
+            ! Need to support possibility of reading more than one
+            ! nambdy_dta from the namelist_cfg internal file.
+            ! Do this by finding the jbdy'th occurence of nambdy_dta in the
+            ! character buffer as the starting point.
+            !
+            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_dta' )
+            IF( nbdy_loc .GT. 0 ) THEN
+               nbdy_rdstart = nbdy_rdstart + nbdy_loc
+            ELSE
+               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',jbdy,' of nambdy_dta not found'
+               ios = -1
+               CALL ctl_nam ( ios , cerrmsg )
+            ENDIF
+            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_dta, IOSTAT = ios, ERR = 902)
 902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist' )
             IF(lwm) WRITE( numond, nambdy_dta )           
@@ -442 +448 @@
             IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file
                CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info
-               CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after
+               CALL fld_def( bf(jp_bdya_i,jbdy) )
+               CALL iom_open( bf(jp_bdya_i,jbdy)%clname, bf(jp_bdya_i,jbdy)%num )
                idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld )
                IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl
                ELSE                                                            ;   ipl = 1            ! xy or xyt
                ENDIF
+               CALL iom_close( bf(jp_bdya_i,jbdy)%num )
+               bf(jp_bdya_i,jbdy)%clrootname = 'NOT USED'   ! reset to default value as this subdomain may not need to read this bdy
             ENDIF
          ENDIF
@@ -487 +496 @@
                igrd = 2                                                    ! U point
                ipk = 1                                                     ! surface data
-               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%u2d will be needed
                llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file
                bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy
                bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta
-               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy
-               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
+               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need u2d over the whole bdy or only over the rim?
+               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)
+               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)
                ENDIF
             ENDIF
@@ -499 +509 @@
                igrd = 3                                                    ! V point
                ipk = 1                                                     ! surface data
-               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%v2d will be needed
                llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file
                bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy
                bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta 
-               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy
-               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
+               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need v2d over the whole bdy or only over the rim?
+               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)
+               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)
                ENDIF
             ENDIF
@@ -615 +626 @@
             ENDIF
 
-            IF( llneed ) THEN                                              ! dta_bdy(jbdy)%xxx will be needed
+            IF( llneed .AND. iszdim > 0 ) THEN                             ! dta_bdy(jbdy)%xxx will be needed
                !                                                           !   -> must be associated with an allocated target
                ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) )              ! allocate the target
@@ -624 +635 @@
                   bf_alias(1)%imap    => idx_bdy(jbdy)%nbmap(1:iszdim,igrd)   ! associate the mapping used for this bdy
                   bf_alias(1)%igrd    = igrd                                  ! used only for vertical integration of 3D arrays
+                  bf_alias(1)%ibdy    = jbdy                                  !  "    "    "     "          "      "  "    "    
                   bf_alias(1)%ltotvel = ln_full_vel                           ! T if u3d is full velocity
                   bf_alias(1)%lzint   = ln_zinterp                            ! T if it requires a vertical interpolation

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdydyn.F90

@@ -37 +37 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn( kt, dyn3d_only )
+   SUBROUTINE bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn  ***
@@ -44 +44 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)           ::   kt           ! Main time step counter
-      LOGICAL, INTENT(in), OPTIONAL ::   dyn3d_only   ! T => only update baroclinic velocities
+      INTEGER                             , INTENT(in)    ::   kt           ! Main time step counter
+      INTEGER                             , INTENT(in)    ::   Kbb, Kaa     ! Ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv     ! Ocean velocities (to be updated at open boundaries)
+      LOGICAL, OPTIONAL                   , INTENT(in)    ::   dyn3d_only   ! T => only update baroclinic velocities
       !
       INTEGER ::   jk, ii, ij, ib_bdy, ib, igrd     ! Loop counter
       LOGICAL ::   ll_dyn2d, ll_dyn3d, ll_orlanski
-      REAL(wp), DIMENSION(jpi,jpj) :: pua2d, pva2d     ! after barotropic velocities
+      REAL(wp), DIMENSION(jpi,jpj) :: zua2d, zva2d     ! after barotropic velocities
       !!----------------------------------------------------------------------
       !
@@ -70 +72 @@
 
       !                          ! "After" velocities: 
-      pua2d(:,:) = 0._wp
-      pva2d(:,:) = 0._wp     
+      zua2d(:,:) = 0._wp
+      zva2d(:,:) = 0._wp     
       DO jk = 1, jpkm1
-         pua2d(:,:) = pua2d(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         pva2d(:,:) = pva2d(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+         zua2d(:,:) = zua2d(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+         zva2d(:,:) = zva2d(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
       END DO
-      pua2d(:,:) = pua2d(:,:) * r1_hu_a(:,:)
-      pva2d(:,:) = pva2d(:,:) * r1_hv_a(:,:)
+      zua2d(:,:) = zua2d(:,:) * r1_hu(:,:,Kaa)
+      zva2d(:,:) = zva2d(:,:) * r1_hv(:,:,Kaa)
 
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) - pua2d(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) - pva2d(:,:) ) * vmask(:,:,jk)
+         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zua2d(:,:) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zva2d(:,:) ) * vmask(:,:,jk)
       END DO
 
@@ -87 +89 @@
       IF( ll_orlanski ) THEN     ! "Before" velocities (Orlanski condition only) 
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = ( ub(:,:,jk) - ub_b(:,:) ) * umask(:,:,jk)
-            vb(:,:,jk) = ( vb(:,:,jk) - vb_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) - uu_b(:,:,Kbb) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) - vv_b(:,:,Kbb) ) * vmask(:,:,jk)
          END DO
       ENDIF
@@ -97 +99 @@
       !-------------------------------------------------------
 
-      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, r1_hu_a(:,:), r1_hv_a(:,:), ssha )
+      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, zua2d, zva2d, uu_b(:,:,Kbb), vv_b(:,:,Kbb), r1_hu(:,:,Kaa), r1_hv(:,:,Kaa), ssh(:,:,Kaa) )
 
-      IF( ll_dyn3d )   CALL bdy_dyn3d( kt )
+      IF( ll_dyn3d )   CALL bdy_dyn3d( kt, Kbb, puu, pvv, Kaa )
 
       !-------------------------------------------------------
@@ -106 +108 @@
       !
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk)
+         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) + zua2d(:,:) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) + zva2d(:,:) ) * vmask(:,:,jk)
       END DO
       !
       IF ( ll_orlanski ) THEN
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = ( ub(:,:,jk) + ub_b(:,:) ) * umask(:,:,jk)
-            vb(:,:,jk) = ( vb(:,:,jk) + vb_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) + uu_b(:,:,Kbb) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) + vv_b(:,:,Kbb) ) * vmask(:,:,jk)
          END DO
       END IF

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdydyn3d.F90

@@ -33 +33 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn3d( kt )
+   SUBROUTINE bdy_dyn3d( kt, Kbb, puu, pvv, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d  ***
@@ -40 +40 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! Main time step counter
+      INTEGER                             , INTENT( in    ) ::   kt        ! Main time step counter
+      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
       !
       INTEGER  ::   ib_bdy, ir     ! BDY set index, rim index
@@ -58 +60 @@
             CASE('none')        ;   CYCLE
             CASE('frs' )        ! treat the whole boundary at once
-               IF( ir == 0) CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+                       IF( ir == 0) CALL bdy_dyn3d_frs( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
             CASE('specified')   ! treat the whole rim      at once
-               IF( ir == 0) CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+                       IF( ir == 0) CALL bdy_dyn3d_spe( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
             CASE('zero')        ! treat the whole rim      at once
-               IF( ir == 0) CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-            CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.false. )
-            CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.true.  )
-            CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy, llrim0 )
-            CASE('neumann')     ;   CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy, llrim0 )
+                       IF( ir == 0) CALL bdy_dyn3d_zro( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+            CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.false. )
+            CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.true.  )
+            CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy, llrim0 )
+            CASE('neumann')     ;   CALL bdy_dyn3d_nmn( puu, pvv, Kaa, idx_bdy(ib_bdy), ib_bdy, llrim0 )
             CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
             END SELECT
@@ -97 +99 @@
          !
          IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', ua, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
+            CALL lbc_lnk( 'bdydyn2d', puu(:,:,:,Kaa), 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
          END IF
          IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', va, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
+            CALL lbc_lnk( 'bdydyn2d', pvv(:,:,:,Kaa), 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
          END IF
       END DO   ! ir
@@ -107 +109 @@
 
 
-   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
+   SUBROUTINE bdy_dyn3d_spe( puu, pvv, Kaa, idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
@@ -115 +117 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER        , INTENT(in) ::   kt      ! time step index
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
-      INTEGER        , INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   kt        ! Time step
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -129 +133 @@
             ii   = idx%nbi(jb,igrd)
             ij   = idx%nbj(jb,igrd)
-            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
+            puu(ii,ij,jk,Kaa) = dta%u3d(jb,jk) * umask(ii,ij,jk)
          END DO
       END DO
@@ -138 +142 @@
             ii   = idx%nbi(jb,igrd)
             ij   = idx%nbj(jb,igrd)
-            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
+            pvv(ii,ij,jk,Kaa) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
          END DO
       END DO
@@ -145 +149 @@
 
 
-   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt, ib_bdy, llrim0 )
+   SUBROUTINE bdy_dyn3d_zgrad( puu, pvv, Kaa, idx, dta, kt, ib_bdy, llrim0 )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
@@ -152 +156 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER                     ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta      ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy   ! BDY set index
-      LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   kt
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
+      LOGICAL                             , INTENT( in    ) ::   llrim0   ! indicate if rim 0 is treated
       !!
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -178 +184 @@
             !
             DO jk = 1, jpkm1
-               ua(ii,ij,jk) = ua(ii,ij+flagv,jk) * umask(ii,ij+flagv,jk)
+               puu(ii,ij,jk,Kaa) = puu(ii,ij+flagv,jk,Kaa) * umask(ii,ij+flagv,jk)
             END DO
             !
@@ -198 +204 @@
             !
             DO jk = 1, jpkm1
-               va(ii,ij,jk) = va(ii+flagu,ij,jk) * vmask(ii+flagu,ij,jk)
+               pvv(ii,ij,jk,Kaa) = pvv(ii+flagu,ij,jk,Kaa) * vmask(ii+flagu,ij,jk)
             END DO
             !
@@ -207 +213 @@
 
 
-   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
+   SUBROUTINE bdy_dyn3d_zro( puu, pvv, Kaa, idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
@@ -214 +220 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER        , INTENT(in) ::   kt      ! time step index
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER                             , INTENT( in    ) ::   kt        ! time step index
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   ib, ik         ! dummy loop indices
@@ -228 +236 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            ua(ii,ij,ik) = 0._wp
+            puu(ii,ij,ik,Kaa) = 0._wp
          END DO
       END DO
@@ -237 +245 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            va(ii,ij,ik) = 0._wp
+            pvv(ii,ij,ik,Kaa) = 0._wp
          END DO
       END DO
@@ -244 +252 @@
 
 
-   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
+   SUBROUTINE bdy_dyn3d_frs( puu, pvv, Kaa, idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
@@ -255 +263 @@
       !!               topography. Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER        , INTENT(in) ::   kt      ! time step index
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER                             , INTENT( in    ) ::   kt        ! time step index
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -271 +281 @@
             ij   = idx%nbj(jb,igrd)
             zwgt = idx%nbw(jb,igrd)
-            ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta%u3d(jb,jk) - ua(ii,ij,jk) ) ) * umask(ii,ij,jk)
+            puu(ii,ij,jk,Kaa) = ( puu(ii,ij,jk,Kaa) + zwgt * ( dta%u3d(jb,jk) - puu(ii,ij,jk,Kaa) ) ) * umask(ii,ij,jk)
          END DO
       END DO
@@ -281 +291 @@
             ij   = idx%nbj(jb,igrd)
             zwgt = idx%nbw(jb,igrd)
-            va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta%v3d(jb,jk) - va(ii,ij,jk) ) ) * vmask(ii,ij,jk)
+            pvv(ii,ij,jk,Kaa) = ( pvv(ii,ij,jk,Kaa) + zwgt * ( dta%v3d(jb,jk) - pvv(ii,ij,jk,Kaa) ) ) * vmask(ii,ij,jk)
          END DO
       END DO   
@@ -288 +298 @@
 
 
-   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, llrim0, ll_npo )
+   SUBROUTINE bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx, dta, ib_bdy, llrim0, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
@@ -298 +308 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
-      INTEGER,                      INTENT(in) ::   ib_bdy   ! BDY set index
-      LOGICAL,                      INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
-      LOGICAL,                      INTENT(in) ::   ll_npo   ! switch for NPO version
+      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
+      LOGICAL                             , INTENT( in    ) ::   llrim0    ! indicate if rim 0 is treated
+      LOGICAL                             , INTENT( in    ) ::   ll_npo    ! switch for NPO version
 
       INTEGER  ::   jb, igrd                               ! dummy loop indices
       !!----------------------------------------------------------------------
       !
-      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
+      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities. 
       !
       igrd = 2      ! Orlanski bc on u-velocity; 
       !            
-      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo, llrim0 )
+      CALL bdy_orlanski_3d( idx, igrd, puu(:,:,:,Kbb), puu(:,:,:,Kaa), dta%u3d, ll_npo, llrim0 )
 
       igrd = 3      ! Orlanski bc on v-velocity
       !  
-      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo, llrim0 )
+      CALL bdy_orlanski_3d( idx, igrd, pvv(:,:,:,Kbb), pvv(:,:,:,Kaa), dta%v3d, ll_npo, llrim0 )
       !
    END SUBROUTINE bdy_dyn3d_orlanski
 
 
-   SUBROUTINE bdy_dyn3d_dmp( kt )
+   SUBROUTINE bdy_dyn3d_dmp( kt, Kbb, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
@@ -327 +339 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step index
+      INTEGER                             , INTENT( in    ) ::   kt        ! time step
+      INTEGER                             , INTENT( in    ) ::   Kbb, Krhs ! Time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities and trends (to be updated at open boundaries)
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -345 +359 @@
                zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
                DO jk = 1, jpkm1
-                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
-                                   ub(ii,ij,jk) + ub_b(ii,ij)) ) * umask(ii,ij,jk)
+                  puu(ii,ij,jk,Krhs) = ( puu(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
+                                   puu(ii,ij,jk,Kbb) + uu_b(ii,ij,Kbb)) ) * umask(ii,ij,jk)
                END DO
             END DO
@@ -356 +370 @@
                zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
                DO jk = 1, jpkm1
-                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
-                                   vb(ii,ij,jk) + vb_b(ii,ij)) ) * vmask(ii,ij,jk)
+                  pvv(ii,ij,jk,Krhs) = ( pvv(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
+                                   pvv(ii,ij,jk,Kbb) + vv_b(ii,ij,Kbb)) ) * vmask(ii,ij,jk)
                END DO
             END DO
@@ -368 +382 @@
 
 
-   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy, llrim0 )
+   SUBROUTINE bdy_dyn3d_nmn( puu, pvv, Kaa, idx, ib_bdy, llrim0 )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
@@ -377 +391 @@
       !!
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices
-      INTEGER,         INTENT(in) ::   ib_bdy   ! BDY set index
-      LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
+      LOGICAL                             , INTENT( in    ) ::   llrim0    ! indicate if rim 0 is treated
       INTEGER  ::   igrd                        ! dummy indice
       !!----------------------------------------------------------------------
       !
-      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
+      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities. 
       !
       igrd = 2      ! Neumann bc on u-velocity; 
       !            
-      CALL bdy_nmn( idx, igrd, ua, llrim0 )   ! ua is masked
+      CALL bdy_nmn( idx, igrd, puu(:,:,:,Kaa), llrim0 )
 
       igrd = 3      ! Neumann bc on v-velocity
       !  
-      CALL bdy_nmn( idx, igrd, va, llrim0 )   ! va is masked
+      CALL bdy_nmn( idx, igrd, pvv(:,:,:,Kaa), llrim0 )
       !
    END SUBROUTINE bdy_dyn3d_nmn

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdyice.F90

@@ -179 +179 @@
 
             ! Then, a) transfer the snow excess into the ice (different from icethd_dh)
-            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )
+            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )
             ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
             !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdyini.F90

@@ -19 +19 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain
+   USE sbc_oce , ONLY: nn_ice
    USE bdy_oce        ! unstructured open boundary conditions
    USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
    USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
-   USE sbctide        ! Tidal forcing or not
-   USE phycst   , ONLY: rday
+   USE tide_mod, ONLY: ln_tide ! tidal forcing
+   USE phycst  , ONLY: rday
    !
    USE in_out_manager ! I/O units
@@ -75 +76 @@
       ! Read namelist parameters
       ! ------------------------
-      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries
       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' )
@@ -93 +93 @@
       cn_ice         (2:jp_bdy) = cn_ice         (1)
       nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
-      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
       READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
@@ -317 +316 @@
 
          dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
+
+         IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN
+            WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice
+            CALL ctl_stop( ctmp1 )
+         ENDIF
 
          IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN 
@@ -364 +368 @@
       ! -------------------------------------------------
 
-      REWIND( numnam_cfg )     
       nblendta(:,:) = 0
       nbdysege = 0
@@ -398 +401 @@
       jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
       ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
-    
+      nbrdta(:,:,:) = 0   ! initialize nbrdta as it may not be completely defined for each bdy
+      
       ! Calculate global boundary index arrays or read in from file
       !------------------------------------------------------------               
@@ -1079 +1083 @@
       INTEGER          ::   ios                 ! Local integer output status for namelist read
       INTEGER          ::   nbdyind, nbdybeg, nbdyend
+      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc
       CHARACTER(LEN=1) ::   ctypebdy   !     -        - 
+      CHARACTER(LEN=50)::   cerrmsg    !     -        - 
       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
       !!----------------------------------------------------------------------
-
-      ! No REWIND here because may need to read more than one nambdy_index namelist.
-      ! Read only namelist_cfg to avoid unseccessfull overwrite 
-      ! keep full control of the configuration namelist
-      READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
+      ! Need to support possibility of reading more than one nambdy_index from
+      ! the namelist_cfg internal file.
+      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the
+      ! character buffer as the starting point.
+      nbdy_rdstart = 1
+      DO nbdy_count = 1, kb_bdy
+       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
+       IF( nbdy_loc .GT. 0 ) THEN
+          nbdy_rdstart = nbdy_rdstart + nbdy_loc
+       ELSE
+          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
+          ios = -1
+          CALL ctl_nam ( ios , cerrmsg )
+       ENDIF
+      END DO
+      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
+      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
 904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
       IF(lwm) WRITE ( numond, nambdy_index )

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdylib.F90

@@ -35 +35 @@
 CONTAINS
 
-   SUBROUTINE bdy_frs( idx, pta, dta )
+   SUBROUTINE bdy_frs( idx, phia, dta )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_frs  ***
@@ -45 +45 @@
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
       REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       !!
       REAL(wp) ::   zwgt           ! boundary weight
@@ -58 +58 @@
             ij = idx%nbj(ib,igrd)
             zwgt = idx%nbw(ib,igrd)
-            pta(ii,ij,ik) = ( pta(ii,ij,ik) + zwgt * (dta(ib,ik) - pta(ii,ij,ik) ) ) * tmask(ii,ij,ik)
+            phia(ii,ij,ik) = ( phia(ii,ij,ik) + zwgt * (dta(ib,ik) - phia(ii,ij,ik) ) ) * tmask(ii,ij,ik)
          END DO
       END DO
@@ -65 +65 @@
 
 
-   SUBROUTINE bdy_spe( idx, pta, dta )
+   SUBROUTINE bdy_spe( idx, phia, dta )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_spe  ***
@@ -74 +74 @@
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
       REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       !!
       INTEGER  ::   ib, ik, igrd   ! dummy loop indices
@@ -85 +85 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            pta(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik)
+            phia(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik)
          END DO
       END DO
@@ -92 +92 @@
 
 
-   SUBROUTINE bdy_orl( idx, ptb, pta, dta, lrim0, ll_npo )
+   SUBROUTINE bdy_orl( idx, phib, phia, dta, lrim0, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_orl  ***
@@ -102 +102 @@
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
       REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptb  ! before tracer field
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phib  ! before tracer field
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       LOGICAL                 , OPTIONAL,  INTENT(in) ::   lrim0   ! indicate if rim 0 is treated
       LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version
@@ -112 +112 @@
       igrd = 1                       ! Everything is at T-points here
       !
-      CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, lrim0, ll_npo )
+      CALL bdy_orlanski_3d( idx, igrd, phib(:,:,:), phia(:,:,:), dta, lrim0, ll_npo )
       !
    END SUBROUTINE bdy_orl
@@ -240 +240 @@
          ! Centred derivative is calculated as average of "left" and "right" derivatives for 
          ! this reason. 
-         ! Note no rdt factor in expression for zdt because it cancels in the expressions for 
+         ! Note no rn_Dt factor in expression for zdt because it cancels in the expressions for 
          ! zrx and zry.
          zdt   =     phia(iibm1   ,ijbm1   ) - phib(iibm1   ,ijbm1   )
@@ -259 +259 @@
          zout = sign( 1., zrx )
          zout = 0.5*( zout + abs(zout) )
-         zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
+         zwgt = 2.*rn_Dt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
          ! only apply radiation on outflow points 
          if( ll_npo ) then     !! NPO version !!
@@ -425 +425 @@
             zout = sign( 1., zrx )
             zout = 0.5*( zout + abs(zout) )
-            zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
+            zwgt = 2.*rn_Dt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
             ! only apply radiation on outflow points 
             if( ll_npo ) then     !! NPO version !!

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdytides.F90

@@ -18 +18 @@
    USE phycst         ! physical constants
    USE bdy_oce        ! ocean open boundary conditions
-   USE tideini        ! 
+   USE tide_mod       ! 
    USE daymod         ! calendar
    !
@@ -30 +30 @@
 
    PUBLIC   bdytide_init     ! routine called in bdy_init
-   PUBLIC   bdytide_update   ! routine called in bdy_dta
    PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
 
@@ -45 +44 @@
    TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
 
+   INTEGER ::   kt_tide
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -62 +65 @@
       !! namelist variables
       !!-------------------
-      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
-      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
-      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
+      CHARACTER(len=80)                         ::   filtide             ! Filename root for tidal input files
+      LOGICAL                                   ::   ln_bdytide_2ddta    ! If true, read 2d harmonic data
       !!
-      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
-      INTEGER                                   ::   ii, ij              !: dummy loop indices
+      INTEGER                                   ::   ib_bdy, itide, ib   ! dummy loop indices
+      INTEGER                                   ::   ii, ij              ! dummy loop indices
       INTEGER                                   ::   inum, igrd
-      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
+      INTEGER                                   ::   isz                 ! bdy data size
       INTEGER                                   ::   ios                 ! Local integer output status for namelist read
-      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file 
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      " 
+      INTEGER                                   ::   nbdy_rdstart, nbdy_loc
+      CHARACTER(LEN=50)                         ::   cerrmsg             ! error string
+      CHARACTER(len=80)                         ::   clfile              ! full file name for tidal input file 
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            ! work space to read in tidal harmonics data
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !  "     "    "   "   "   "        "      " 
       !!
-      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
+      TYPE(TIDES_DATA), POINTER                 ::   td                  ! local short cut   
+      TYPE(  OBC_DATA), POINTER                 ::   dta                 ! local short cut
       !!
-      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
+      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
       !!----------------------------------------------------------------------
       !
@@ -84 +89 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
 
-      REWIND(numnam_cfg)
-
+
+      nbdy_rdstart = 1
       DO ib_bdy = 1, nb_bdy
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
             !
-            td => tides(ib_bdy)
-
+            td  => tides(ib_bdy)
+            dta => dta_bdy(ib_bdy)
+         
             ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
             filtide(:) = ''
 
-            REWIND( numnam_ref )
             READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
-            ! Don't REWIND here - may need to read more than one of these namelists. 
-            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
+            !
+            ! Need to support possibility of reading more than one
+            ! nambdy_tide from the namelist_cfg internal file.
+            ! Do this by finding the ib_bdy'th occurence of nambdy_tide in the
+            ! character buffer as the starting point.
+            !
+            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' )
+            IF( nbdy_loc .GT. 0 ) THEN
+               nbdy_rdstart = nbdy_rdstart + nbdy_loc
+            ELSE
+               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found'
+               ios = -1
+               CALL ctl_nam ( ios , cerrmsg )
+            ENDIF
+            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_tide, IOSTAT = ios, ERR = 902)
 902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
             IF(lwm) WRITE ( numond, nambdy_tide )
@@ -105 +123 @@
             IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
             IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
-            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
             IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
             IF(lwp) THEN 
                     WRITE(numout,*) '             Tidal components: ' 
                DO itide = 1, nb_harmo
-                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
+                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide 
                END DO
             ENDIF 
             IF(lwp) WRITE(numout,*) ' '
 
-            ! Allocate space for tidal harmonics data - get size from OBC data arrays
+            ! Allocate space for tidal harmonics data - get size from BDY data arrays
+            ! Allocate also slow varying data in the case of time splitting:
+            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
             ! -----------------------------------------------------------------------
-
-            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
-            ! relaxation area      
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
-            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
-            ENDIF
-
-            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
-            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
-
-            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
-            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
-
-            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
-            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
-
-            td%ssh0(:,:,:) = 0._wp
-            td%ssh (:,:,:) = 0._wp
-            td%u0  (:,:,:) = 0._wp
-            td%u   (:,:,:) = 0._wp
-            td%v0  (:,:,:) = 0._wp
-            td%v   (:,:,:) = 0._wp
-
+            IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+               isz = SIZE(dta%ssh)
+               ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) )
+               dta_bdy_s(ib_bdy)%ssh(:) = 0._wp   ! needed?
+            ENDIF
+            IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+               isz = SIZE(dta%u2d)
+               ALLOCATE( td%u0  ( isz, nb_harmo, 2 ), td%u  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) )
+               dta_bdy_s(ib_bdy)%u2d(:) = 0._wp   ! needed?
+            ENDIF
+            IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+               isz = SIZE(dta%v2d)
+               ALLOCATE( td%v0  ( isz, nb_harmo, 2 ), td%v  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) )
+               dta_bdy_s(ib_bdy)%v2d(:) = 0._wp   ! needed?
+            ENDIF
+
+            ! fill td%ssh0, td%u0, td%v0
+            ! -----------------------------------------------------------------------
             IF( ln_bdytide_2ddta ) THEN
+               !
                ! It is assumed that each data file contains all complex harmonic amplitudes
                ! given on the global domain (ie global, jpiglo x jpjglo)
@@ -147 +162 @@
                !
                ! SSH fields
-               clfile = TRIM(filtide)//'_grid_T.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 1                       ! Everything is at T-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%ssh0(ib,itide,1) = ztr(ii,ij)
-                     td%ssh0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO 
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_T.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 1                       ! Everything is at T-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
+                     CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
+                     DO ib = 1, SIZE(dta%ssh)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%ssh0(ib,itide,1) = ztr(ii,ij)
+                        td%ssh0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                ! U fields
-               clfile = TRIM(filtide)//'_grid_U.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 2                       ! Everything is at U-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%u0(ib,itide,1) = ztr(ii,ij)
-                     td%u0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_U.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 2                       ! Everything is at U-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
+                     DO ib = 1, SIZE(dta%u2d)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%u0(ib,itide,1) = ztr(ii,ij)
+                        td%u0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                ! V fields
-               clfile = TRIM(filtide)//'_grid_V.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 3                       ! Everything is at V-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%v0(ib,itide,1) = ztr(ii,ij)
-                     td%v0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO  
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_V.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 3                       ! Everything is at V-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
+                     DO ib = 1, SIZE(dta%v2d)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%v0(ib,itide,1) = ztr(ii,ij)
+                        td%v0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                DEALLOCATE( ztr, zti ) 
@@ -203 +221 @@
                ! Read tidal data only on bdy segments
                ! 
-               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
+               ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) )
                !
                ! Open files and read in tidal forcing data
@@ -210 +228 @@
                DO itide = 1, nb_harmo
                   !                                                              ! SSH fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
-                  CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
-                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
-                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
-                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+                     isz = SIZE(dta%ssh)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                     td%ssh0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                     td%ssh0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !                                                              ! U fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
-                  CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
-                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
-                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
-                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+                     isz = SIZE(dta%u2d)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                     td%u0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                     td%u0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !                                                              ! V fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
-                  CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
-                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
-                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
-                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+                     isz = SIZE(dta%v2d)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                     td%v0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                     td%v0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !
                END DO ! end loop on tidal components
@@ -240 +267 @@
             ENDIF ! ln_bdytide_2ddta=.true.
             !
-            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
-               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
-               td%u0  (:,:,2) = - td%u0  (:,:,2)
-               td%v0  (:,:,2) = - td%v0  (:,:,2)
-            ENDIF
-            !
-            ! Allocate slow varying data in the case of time splitting:
-            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
-            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
-            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
-            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
-            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
-            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
-            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
-            !
          ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
          !
@@ -262 +274 @@
 
 
-   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdytide_update  ***
-      !!                
-      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
-      !!                
-      !!----------------------------------------------------------------------
-      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
-      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
-      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
-      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
-      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
-      !                                                 ! is present then units = subcycle timesteps.
-      !                                                 ! kt_offset = 0  => get data at "now"    time level
-      !                                                 ! kt_offset = -1 => get data at "before" time level
-      !                                                 ! kt_offset = +1 => get data at "after"  time level
-      !                                                 ! etc.
-      !
-      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER  ::   time_add              ! time offset in units of timesteps
-      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
-      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars    
-      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
-      !!----------------------------------------------------------------------
-      !
-      ilen0(1) =  SIZE(td%ssh(:,1,1))
-      ilen0(2) =  SIZE(td%u(:,1,1))
-      ilen0(3) =  SIZE(td%v(:,1,1))
-
-      zflag=1
-      IF ( PRESENT(kit) ) THEN
-        IF ( kit /= 1 ) zflag=0
-      ENDIF
-
-      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
-        !
-        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
-        !
-        IF(lwp) THEN
-           WRITE(numout,*)
-           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
-           WRITE(numout,*) '~~~~~~~~~~~~~~ '
-        ENDIF
-        !
-        CALL tide_init_elevation ( idx, td )
-        CALL tide_init_velocities( idx, td )
-        !
-      ENDIF 
-
-      time_add = 0
-      IF( PRESENT(kt_offset) ) THEN
-         time_add = kt_offset
-      ENDIF
-         
-      IF( PRESENT(kit) ) THEN  
-         z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
-      ELSE                              
-         z_arg = ((kt-kt_tide)+time_add) * rdt
-      ENDIF
-
-      ! Linear ramp on tidal component at open boundaries 
-      zramp = 1._wp
-      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
-
-      DO itide = 1, nb_harmo
-         z_sarg = z_arg * omega_tide(itide)
-         z_cost(itide) = COS( z_sarg )
-         z_sist(itide) = SIN( z_sarg )
-      END DO
-
-      DO itide = 1, nb_harmo
-         igrd=1                              ! SSH on tracer grid
-         DO ib = 1, ilen0(igrd)
-            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
-         END DO
-         igrd=2                              ! U grid
-         DO ib = 1, ilen0(igrd)
-            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
-         END DO
-         igrd=3                              ! V grid
-         DO ib = 1, ilen0(igrd) 
-            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
-         END DO
-      END DO
-      !
-   END SUBROUTINE bdytide_update
-
-
-   SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
+   SUBROUTINE bdy_dta_tides( kt, kit, pt_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dta_tides  ***
@@ -360 +283 @@
       INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
       INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
-      !                                              ! is present then units = subcycle timesteps.
-      !                                              ! kt_offset = 0  => get data at "now"    time level
-      !                                              ! kt_offset = -1 => get data at "before" time level
-      !                                              ! kt_offset = +1 => get data at "after"  time level
-      !                                              ! etc.
+      REAL(wp),OPTIONAL, INTENT(in) ::   pt_offset   ! time offset in units of timesteps
       !
       LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
-      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
-      INTEGER  ::   time_add                  ! time offset in units of timesteps
-      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
-      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
-      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist      
+      INTEGER  ::   itide, ib_bdy, ib         ! loop indices
+      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset   
       !!----------------------------------------------------------------------
       !
@@ -378 +293 @@
       IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
 
-      time_add = 0
-      IF( PRESENT(kt_offset) ) THEN
-         time_add = kt_offset
-      ENDIF
+      zt_offset = 0._wp
+      IF( PRESENT(pt_offset) )   zt_offset = pt_offset
       
       ! Absolute time from model initialization:   
       IF( PRESENT(kit) ) THEN  
-         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
+         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_e, wp) ) * rn_Dt
       ELSE                              
-         z_arg = ( kt + time_add ) * rdt
+         z_arg = ( REAL(kt, wp) + zt_offset ) * rn_Dt
       ENDIF
 
       ! Linear ramp on tidal component at open boundaries 
       zramp = 1.
-      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - REAL(nit000,wp)*rn_Dt)/(rn_tide_ramp_dt*rday),0.),1.)
 
       DO ib_bdy = 1,nb_bdy
          !
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
-            !
-            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
-            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
-            !
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
-            ELSE                                   ;   ilen0(:) = nblenrim(:)
-            ENDIF     
             !
             ! We refresh nodal factors every day below
             ! This should be done somewhere else
-            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
-               !
-               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
+            IF ( ( nsec_day == NINT(0.5_wp * rn_Dt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
+               !
+               kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rn_Dt)/rn_Dt)
                !
                IF(lwp) THEN
@@ -421 +327 @@
                !
             ENDIF
-            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
+            zoff = REAL(-kt_tide,wp) * rn_Dt ! time offset relative to nodal factor computation time
             !
             ! If time splitting, initialize arrays from slow varying open boundary data:
             IF ( PRESENT(kit) ) THEN           
-               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:)
             ENDIF
             !
@@ -433 +339 @@
             DO itide = 1, nb_harmo
                !
-               z_sarg = (z_arg + zoff) * omega_tide(itide)
+               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega
                z_cost = zramp * COS( z_sarg )
                z_sist = zramp * SIN( z_sarg )
                !
-               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
-                  igrd=1                              ! SSH on tracer grid
-                  DO ib = 1, ilen0(igrd)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN   ! SSH on tracer grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh)
                      dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
                         &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
@@ -446 +351 @@
                ENDIF
                !
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
-                  igrd=2                              ! U grid
-                  DO ib = 1, ilen0(igrd)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN  ! U grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d)
                      dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
                         &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
                         &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
                   END DO
-                  igrd=3                              ! V grid
-                  DO ib = 1, ilen0(igrd) 
+               ENDIF
+               !
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN   ! V grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d)
                      dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
                         &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
@@ -460 +366 @@
                   END DO
                ENDIF
+               !
             END DO             
-         END IF
+         ENDIF
       END DO
       !
@@ -474 +381 @@
       TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
       !
-      INTEGER ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
+      INTEGER ::   itide, isz, ib       ! dummy loop indices
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       !!----------------------------------------------------------------------
       !
-      igrd=1   
-                              ! SSH on tracer grid.
-      ilen0(1) =  SIZE(td%ssh0(:,1,1))
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
-      !
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
+      IF( ASSOCIATED(td%ssh0) ) THEN   ! SSH on tracer grid.
+         !
+         isz = SIZE( td%ssh0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1 , ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
-         ENDDO
-         DO ib = 1 , ilen0(igrd)
-            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide, phi_tide )
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
       !
    END SUBROUTINE tide_init_elevation
@@ -512 +419 @@
       TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
       !
-      INTEGER ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
+      INTEGER ::   itide, isz, ib        ! dummy loop indices
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       !!----------------------------------------------------------------------
       !
-      ilen0(2) =  SIZE(td%u0(:,1,1))
-      ilen0(3) =  SIZE(td%v0(:,1,1))
-      !
-      igrd=2                                 ! U grid.
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
-      !
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
+      IF( ASSOCIATED(td%u0) ) THEN   ! U grid. we use bdy u2d on this mpi subdomain
+         !
+         isz = SIZE( td%u0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
-         ENDDO
-         DO ib = 1, ilen0(igrd)
-            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide , phi_tide )
-      !
-      igrd=3                                 ! V grid.
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
-
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
+      !
+      IF( ASSOCIATED(td%v0) ) THEN   ! V grid. we use bdy u2d on this mpi subdomain
+         !
+         isz = SIZE( td%v0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
-         ENDDO
-         DO ib = 1, ilen0(igrd)
-            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide, phi_tide )
-      !
-  END SUBROUTINE tide_init_velocities
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
+      !
+   END SUBROUTINE tide_init_velocities
 
    !!======================================================================

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdytra.F90

@@ -40 +40 @@
 CONTAINS
 
-   SUBROUTINE bdy_tra( kt )
+   SUBROUTINE bdy_tra( kt, Kbb, pts, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_tra  ***
@@ -47 +47 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! Main time step counter
+      INTEGER                                  , INTENT(in)    :: kt        ! Main time step counter
+      INTEGER                                  , INTENT(in)    :: Kbb, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! tracer fields
       !
       INTEGER                        :: ib_bdy, jn, igrd, ir   ! Loop indeces
@@ -70 +72 @@
                CASE('none'        )   ;   CYCLE
                CASE('frs'         )   ! treat the whole boundary at once
-                  IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra )
+                  IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
                CASE('specified'   )   ! treat the whole rim      at once
-                  IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra )
-               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , tsa(:,:,:,jn), llrim0 )   ! tsa masked
-               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), &
+                  IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
+               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , pts(:,:,:,jn,Kaa), llrim0 )   ! tsa masked
+               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), &
                     & zdta(jn)%tra, llrim0, ll_npo=.false. )
-               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), &
+               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), &
                     & zdta(jn)%tra, llrim0, ll_npo=.true.  )
-               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), jn, llrim0 )
+               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), jn, llrim0 )
                CASE DEFAULT           ;   CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
                END SELECT
@@ -98 +100 @@
          END DO
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdytra', tsa, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+            CALL lbc_lnk( 'bdytra', pts(:,:,:,jn,Kaa), 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
          END IF
          !
@@ -106 +108 @@
 
 
-   SUBROUTINE bdy_rnf( idx, pta, jpa, llrim0 )
+   SUBROUTINE bdy_rnf( idx, pt, jpa, llrim0 )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_rnf  ***
@@ -116 +118 @@
       !!----------------------------------------------------------------------
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx      ! OBC indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta      ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pt       ! tracer trend
       INTEGER,                             INTENT(in) ::   jpa      ! TRA index
       LOGICAL,                             INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
       !
       INTEGER  ::   ib, ii, ij, igrd   ! dummy loop indices
-      INTEGER  ::   ik, ip, jp ! 2D addresses
       !!----------------------------------------------------------------------
       !
       igrd = 1                       ! Everything is at T-points here
       IF(      jpa == jp_tem ) THEN
-         CALL bdy_nmn( idx, igrd, pta, llrim0 )
+         CALL bdy_nmn( idx, igrd, pt, llrim0 )
       ELSE IF( jpa == jp_sal ) THEN
          IF( .NOT. llrim0 )   RETURN
@@ -132 +133 @@
             ii = idx%nbi(ib,igrd)
             ij = idx%nbj(ib,igrd)
-            pta(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1)
+            pt(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1)
          END DO
       END IF
@@ -139 +140 @@
 
 
-   SUBROUTINE bdy_tra_dmp( kt )
+   SUBROUTINE bdy_tra_dmp( kt, Kbb, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_tra_dmp  ***
@@ -146 +147 @@
       !! 
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   !
+      INTEGER                                  , INTENT(in)    :: kt        ! time step
+      INTEGER                                  , INTENT(in)    :: Kbb, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation
       !
       REAL(wp) ::   zwgt           ! boundary weight
@@ -165 +168 @@
                zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
                DO ik = 1, jpkm1
-                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik)
-                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik)
-                  tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta
-                  tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa
+                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - pts(ii,ij,ik,jp_tem,Kbb) ) * tmask(ii,ij,ik)
+                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - pts(ii,ij,ik,jp_sal,Kbb) ) * tmask(ii,ij,ik)
+                  pts(ii,ij,ik,jp_tem,Krhs) = pts(ii,ij,ik,jp_tem,Krhs) + zta
+                  pts(ii,ij,ik,jp_sal,Krhs) = pts(ii,ij,ik,jp_sal,Krhs) + zsa
                END DO
             END DO

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/BDY/bdyvol.F90

@@ -14 +14 @@
    USE bdy_oce        ! ocean open boundary conditions
    USE sbc_oce        ! ocean surface boundary conditions
+   USE isf_oce, ONLY : fwfisf_cav, fwfisf_par  ! ice shelf
    USE dom_oce        ! ocean space and time domain 
    USE phycst         ! physical constants
-   USE sbcisf         ! ice shelf
    !
    USE in_out_manager ! I/O manager
@@ -77 +77 @@
       ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
       ! -----------------------------------------------------------------------
-      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
+      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rho0
 
       ! Compute bdy surface each cycle if non linear free surface
@@ -143 +143 @@
       ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
       ! ------------------------------------------------------
-      IF( MOD( kt, nn_write ) == 0 .AND. ( kc == 1 ) ) THEN
+      IF( MOD( kt, MAX(nn_write,1) ) == 0 .AND. ( kc == 1 ) ) THEN
          !
          ! compute residual transport across boundary

@@ -6 +6 @@
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
+
+# SETTE
+^/utils/CI/sette@HEAD         sette