Changeset 3367 for branches

Timestamp:

2012-04-25T12:21:21+02:00 (12 years ago)

Author:

greffray

Message:

Merge tidal packages
Merge OBC-BDY packages
Add atmospheric pressure at the open boundaries
Modification of the namelist for AMM12 configuration

Location:

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM

Files:

• CONFIG/AMM12/EXP00/namelist (modified) (2 diffs)
• NEMO/OPA_SRC/BDY/bdy_oce.F90 (modified) (1 diff)
• NEMO/OPA_SRC/BDY/bdydta.F90 (modified) (24 diffs)
• NEMO/OPA_SRC/BDY/bdydyn3d.F90 (modified) (5 diffs)
• NEMO/OPA_SRC/BDY/bdyini.F90 (modified) (13 diffs)
• NEMO/OPA_SRC/BDY/bdytides.F90 (modified) (11 diffs)
• NEMO/OPA_SRC/BDY/bdytra.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/SBC/sbcapr.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/SBC/sbctide.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SBC/tideini.F90 (added)
• NEMO/OPA_SRC/nemogcm.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/step.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/step_oce.F90 (modified) (1 diff)

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/CONFIG/AMM12/EXP00/namelist

@@ -371 +371 @@
 &nam_tide      !   tide parameters (#ifdef key_tide)
 !-----------------------------------------------------------------------
-   ln_tide_pot   = .true.   !  use tidal potential forcing
-   nb_harmo      =    11    !  number of constituents used
-   clname(1)     =   'M2'   !  name of constituent
-   clname(2)     =   'S2'
-   clname(3)     =   'N2'
-   clname(4)     =   'K1'
-   clname(5)     =   'O1'
-   clname(6)     =   'Q1'
-   clname(7)     =   'M4'
-   clname(8)     =   'K2'
-   clname(9)     =   'P1'
-   clname(10)    =   'Mf'
-   clname(11)    =   'Mm'
+   ln_tide_pot   = .false.   !  use tidal potential forcing
+   clname(1)     =   'Q1'   !  name of constituent
+   clname(2)     =   'O1'
+   clname(3)     =   'P1'
+   clname(4)     =   'S1'
+   clname(5)     =   'K1'
+   clname(6)     =   '2N2'
+   clname(7)     =   'MU2'
+   clname(8)     =   'N2'
+   clname(9)     =   'NU2'
+   clname(10)    =   'M2'
+   clname(11)    =   'L2'
+   clname(12)    =   'T2'
+   clname(13)    =   'S2'
+   clname(14)    =   'K2'
+   clname(15)    =   'M4'
 /
 !-----------------------------------------------------------------------
@@ -427 +430 @@
 !-----------------------------------------------------------------------
    filtide      = 'bdydta/amm12_bdytide_'         !  file name root of tidal forcing files
-    tide_cpt(1)   ='Q1'  !  names of tidal components used
-    tide_cpt(2)   ='O1'  !  names of tidal components used
-    tide_cpt(3)   ='P1'  !  names of tidal components used
-    tide_cpt(4)   ='S1'  !  names of tidal components used
-    tide_cpt(5)   ='K1'  !  names of tidal components used
-    tide_cpt(6)   ='2N2' !  names of tidal components used
-    tide_cpt(7)   ='MU2' !  names of tidal components used
-    tide_cpt(8)   ='N2'  !  names of tidal components used
-    tide_cpt(9)   ='NU2' !  names of tidal components used
-    tide_cpt(10)   ='M2'  !  names of tidal components used
-    tide_cpt(11)   ='L2'  !  names of tidal components used
-    tide_cpt(12)   ='T2'  !  names of tidal components used
-    tide_cpt(13)   ='S2'  !  names of tidal components used
-    tide_cpt(14)   ='K2'  !  names of tidal components used
-    tide_cpt(15)   ='M4'  !  names of tidal components used
-    tide_speed(1)   = 13.398661 !  phase speeds of tidal components (deg/hour)
-    tide_speed(2)   = 13.943036 !  phase speeds of tidal components (deg/hour)
-    tide_speed(3)   = 14.958932 !  phase speeds of tidal components (deg/hour)
-    tide_speed(4)   = 15.000001 !  phase speeds of tidal components (deg/hour)
-    tide_speed(5)   = 15.041069 !  phase speeds of tidal components (deg/hour)
-    tide_speed(6)   = 27.895355 !  phase speeds of tidal components (deg/hour)
-    tide_speed(7)   = 27.968210 !  phase speeds of tidal components (deg/hour)
-    tide_speed(8)   = 28.439730 !  phase speeds of tidal components (deg/hour)
-    tide_speed(9)   = 28.512585 !  phase speeds of tidal components (deg/hour)
-    tide_speed(10)   = 28.984106 !  phase speeds of tidal components (deg/hour)
-    tide_speed(11)   = 29.528479 !  phase speeds of tidal components (deg/hour)
-    tide_speed(12)   = 29.958935 !  phase speeds of tidal components (deg/hour)
-    tide_speed(13)   = 30.000002 !  phase speeds of tidal components (deg/hour)
-    tide_speed(14)   = 30.082138 !  phase speeds of tidal components (deg/hour)
-    tide_speed(15)   = 57.968212 !  phase speeds of tidal components (deg/hour)
-    ln_tide_date = .true.               !  adjust tidal harmonics for start date of run
 /
 !!======================================================================

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

@@ -73 +73 @@
    INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta             !: = 0 use the initial state as bdy dta ; 
                                                             !: = 1 read it in a NetCDF file
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp
+   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp
+   LOGICAL                    ::   ln_tra_dmp_tot
+   LOGICAL                    ::   ln_dyn3d_dmp_tot
+
 #if defined key_lim2
    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables 

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

@@ -32 +32 @@
    USE ice_2
 #endif
+   USE sbcapr
 
    IMPLICIT NONE
@@ -108 +109 @@
 
             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 1
                DO ib = 1, ilen1(igrd)
@@ -134 +131 @@
 
             IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 2 
                DO ib = 1, ilen1(igrd)
@@ -158 +151 @@
 
             IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 1                       ! Everything is at T-points here
                DO ib = 1, ilen1(igrd)
@@ -176 +165 @@
 #if defined key_lim2
             IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 1                       ! Everything is at T-points here
                DO ib = 1, ilen1(igrd)
@@ -214 +199 @@
                      dta_bdy(ib_bdy)%v2d(:) = 0.0
                   ENDIF
-                  IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data
-                     jend = jstart + 2
-                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),   &
-                     &              jit=jit, time_offset=time_offset )
-                  ENDIF
-                  IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                     CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   & 
-                     &                 jit=jit, time_offset=time_offset )
+                  IF (nn_tra(ib_bdy).ne.4) THEN
+                   IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN
+                     ! For the runoff case, no need to update the forcing (already done in the baroclinic part)
+                     jend = nb_bdy_fld(ib_bdy)
+                     IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2
+                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset )
+                     IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2
+                     ! If full velocities in boundary data then split into barotropic and baroclinic data
+                     IF( ln_full_vel_array(ib_bdy) .and.                                             &
+                    &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN
+                        igrd = 2                      ! zonal velocity
+                        dta_bdy(ib_bdy)%u2d(:) = 0.0
+                        DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                           ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                           ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                           DO ik = 1, jpkm1
+                              dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
+                    &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                           END DO
+                           dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
+                           DO ik = 1, jpkm1
+                              dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)
+                           END DO
+                        END DO
+                        igrd = 3                      ! meridional velocity
+                        dta_bdy(ib_bdy)%v2d(:) = 0.0
+                        DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                           ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                           ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                           DO ik = 1, jpkm1
+                              dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
+                    &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
+                           END DO
+                           dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
+                           DO ik = 1, jpkm1
+                              dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)
+                           END DO
+                        END DO
+                     ENDIF                    
+                   ENDIF
+                   IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                      CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   & 
+                      &                 jit=jit, time_offset=time_offset )
+                   ENDIF
                   ENDIF
                ENDIF
             ELSE
-               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                  dta_bdy(ib_bdy)%ssh(:) = 0.0
-                  dta_bdy(ib_bdy)%u2d(:) = 0.0
-                  dta_bdy(ib_bdy)%v2d(:) = 0.0
+               IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition
+                  jend = nb_bdy_fld(ib_bdy)
+                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
+                  !
+                  igrd = 2                      ! zonal velocity
+                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+                  END DO
+                  !
+                  igrd = 3                      ! meridional velocity
+                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+                  END DO
+               ELSE
+                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     dta_bdy(ib_bdy)%ssh(:) = 0.0
+                     dta_bdy(ib_bdy)%u2d(:) = 0.0
+                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                  ENDIF
+                  IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data
+                     jend = nb_bdy_fld(ib_bdy)
+                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
+                  ENDIF
+                  ! If full velocities in boundary data then split into barotropic and baroclinic data
+                  IF( ln_full_vel_array(ib_bdy) .and.                                             &
+                 &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN
+                     igrd = 2                      ! zonal velocity
+                     dta_bdy(ib_bdy)%u2d(:) = 0.0
+                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        DO ik = 1, jpkm1
+                           dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
+                 &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                        END DO
+                        dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
+                        DO ik = 1, jpkm1
+                           dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)
+                        END DO
+                     END DO
+                     igrd = 3                      ! meridional velocity
+                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        DO ik = 1, jpkm1
+                           dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
+                 &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
+                        END DO
+                        dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
+                        DO ik = 1, jpkm1
+                           dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)
+                        END DO
+                     END DO
+                  ENDIF
+                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                     CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )
+                  ENDIF
                ENDIF
-               IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data
-                  jend = jstart + nb_bdy_fld(ib_bdy) - 1
-                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
-               ENDIF
-               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                  CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )
-               ENDIF
             ENDIF
             jstart = jend+1
-
-            ! If full velocities in boundary data then split into barotropic and baroclinic data
-            ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same
-            ! time as the dynspg_ts option). 
-
-            IF( ln_full_vel_array(ib_bdy) .and.                                             & 
-           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 
-
-               igrd = 2                      ! zonal velocity
-               dta_bdy(ib_bdy)%u2d(:) = 0.0
-               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+         END IF ! nn_dta(ib_bdy) = 1
+      END DO  ! ib_bdy
+
+      IF ( ln_apr_obc ) THEN
+         DO ib_bdy = 1, nb_bdy
+            IF (nn_tra(ib_bdy).NE.4)THEN
+               igrd = 1                      ! meridional velocity
+               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
                   ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                   ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  DO ik = 1, jpkm1
-                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
-              &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
-                  END DO
-                  dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
-                  DO ik = 1, jpkm1
-                     dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
-                  END DO
-               END DO
-
-               igrd = 3                      ! meridional velocity
-               dta_bdy(ib_bdy)%v2d(:) = 0.0
-               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  DO ik = 1, jpkm1
-                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
-              &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
-                  END DO
-                  dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
-                  DO ik = 1, jpkm1
-                     dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
-                  END DO
-               END DO
-    
-            ENDIF
-
-         END IF ! nn_dta(ib_bdy) = 1
-      END DO  ! ib_bdy
+                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij)
+               ENDDO
+            ENDIF
+         ENDDO
+      ENDIF
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dta')
@@ -326 +371 @@
       IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init')
 
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) ''
+
       ! Set nn_dta
       DO ib_bdy = 1, nb_bdy
-         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       &
-                               ,nn_dyn3d_dta(ib_bdy)       &
-                               ,nn_tra_dta(ib_bdy)         &
-#if defined key_lim2
-                               ,nn_ice_lim2_dta(ib_bdy)    &
-#endif
-                              )
-         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1
+         IF (nn_dyn2d_dta(ib_bdy).eq.1.OR.nn_dyn2d_dta(ib_bdy).eq.3) nn_dta(ib_bdy) = 1
+         IF (nn_dyn3d_dta(ib_bdy).ge.1) nn_dta(ib_bdy) = 1
+         IF (nn_tra_dta  (ib_bdy).ge.1) nn_dta(ib_bdy) = 1
+#if defined key_lim2
+         IF (nn_ice_lim2_dta(ib_bdy).eq.1) nn_dta(ib_bdy) = 1
+#endif
+         IF (lwp) THEN
+            WRITE(numout,*) 'Bdy package number ',ib_bdy
+            IF (nn_dta(ib_bdy).eq.1) THEN
+               WRITE(numout,*) 'We use external data'
+            ELSE
+               WRITE(numout,*) 'We do not use external data'
+            ENDIF
+         ENDIF
       END DO
 
@@ -357 +413 @@
          ENDIF
 #endif               
+         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy)
       ENDDO            
 
@@ -408 +465 @@
             ln_full_vel_array(ib_bdy) = ln_full_vel
 
-            IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN
-               CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',&
-            &                  'with dynspg_ts option' )   ;   RETURN  
-            ENDIF             
-
             nblen => idx_bdy(ib_bdy)%nblen
             nblenrim => idx_bdy(ib_bdy)%nblenrim
@@ -420 +472 @@
             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
 
-               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
+               IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
                   ibdy(jfld) = ib_bdy
                   igrid(jfld) = 1
-                  ilen1(jfld) = nblenrim(igrid(jfld))
+                  ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
                ENDIF
 
                IF( .not. ln_full_vel_array(ib_bdy) ) THEN
-
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u2d
                   ibdy(jfld) = ib_bdy
                   igrid(jfld) = 2
-                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                     ilen1(jfld) = nblen(igrid(jfld))
-                  ELSE
-                     ilen1(jfld) = nblenrim(igrid(jfld))
-                  ENDIF
+                  ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
 
@@ -446 +493 @@
                   ibdy(jfld) = ib_bdy
                   igrid(jfld) = 3
-                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                     ilen1(jfld) = nblen(igrid(jfld))
-                  ELSE
-                     ilen1(jfld) = nblenrim(igrid(jfld))
-                  ENDIF
+                  ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
-
                ENDIF
 
@@ -466 +508 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 2
-               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -477 +515 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 3
-               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -493 +527 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -504 +534 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -521 +547 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = 1
 
@@ -532 +554 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = 1
 
@@ -543 +561 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = 1
 
@@ -566 +580 @@
       ENDDO ! ib_bdy
 
-
       DO jfld = 1, nb_bdy_fld_sum
          ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) )
@@ -577 +590 @@
       jstart = 1
       DO ib_bdy = 1, nb_bdy
-         jend = jstart + nb_bdy_fld(ib_bdy) - 1
+         jend = nb_bdy_fld(ib_bdy) 
          CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   &
          &              'open boundary conditions', 'nambdy_dta' )
@@ -596 +609 @@
          IF (nn_dyn2d(ib_bdy) .gt. 0) THEN
             IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN
-               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen0(1:3) = nblen(1:3)
-               ELSE
-                  ilen0(1:3) = nblenrim(1:3)
-               ENDIF
-               ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) )
+               ilen0(1:3) = nblen(1:3)
                ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) )
                ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) )
+               IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN
+                  jfld = jfld + 1
+                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)
+               ELSE
+                  ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) )
+               ENDIF
             ELSE
                IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
@@ -617 +631 @@
 
          IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN
-            IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-               ilen0(1:3) = nblen(1:3)
-            ELSE
-               ilen0(1:3) = nblenrim(1:3)
-            ENDIF
+            ilen0(1:3) = nblen(1:3)
             ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) )
             ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) )
@@ -636 +646 @@
          IF (nn_tra(ib_bdy) .gt. 0) THEN
             IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen0(1:3) = nblen(1:3)
-               ELSE
-                  ilen0(1:3) = nblenrim(1:3)
-               ENDIF
+               ilen0(1:3) = nblen(1:3)
                ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) )
                ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) )
@@ -654 +660 @@
          IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN
             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen0(1:3) = nblen(1:3)
-               ELSE
-                  ilen0(1:3) = nblenrim(1:3)
-               ENDIF
+               ilen0(1:3) = nblen(1:3)
                ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) )
                ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) )

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

@@ -19 +19 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
+   Use phycst
 
    IMPLICIT NONE
@@ -24 +25 @@
 
    PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn
+   PUBLIC   bdy_dyn3d_dmp ! routine called by step
 
    !!----------------------------------------------------------------------
@@ -55 +57 @@
          CASE(jp_frs)
             CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(2)
+            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(3)
+            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
          CASE DEFAULT
             CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
@@ -62 +68 @@
    END SUBROUTINE bdy_dyn3d
 
+   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
+      !!
+      !! ** Purpose : - Apply a specified value for baroclinic velocities
+      !!                at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !!
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe')
+      !
+      igrd = 2                      ! Relaxation of zonal velocity
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            ij   = idx%nbj(jb,igrd)
+            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
+         END DO
+      END DO
+      !
+      igrd = 3                      ! Relaxation of meridional velocity
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            ij   = idx%nbj(jb,igrd)
+            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
+         END DO
+      END DO
+      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe')
+
+   END SUBROUTINE bdy_dyn3d_spe
+
+   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
+      !!
+      !! ** Purpose : - baroclinic velocities = 0. at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !!
+      INTEGER  ::   ib, ik         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro')
+      !
+      igrd = 2                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ua(ii,ij,ik) = 0._wp
+         END DO
+      END DO
+
+      igrd = 3                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            va(ii,ij,ik) = 0._wp
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro')
+
+   END SUBROUTINE bdy_dyn3d_zro
+
    SUBROUTINE bdy_dyn3d_frs( idx, dta, kt )
       !!----------------------------------------------------------------------
@@ -111 +205 @@
    END SUBROUTINE bdy_dyn3d_frs
 
+   SUBROUTINE bdy_dyn3d_dmp( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
+      !!
+      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      !!
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::  ib_bdy          ! loop index
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
+      !
+      DO ib_bdy=1, nb_bdy
+         IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN
+            igrd = 2                      ! Relaxation of zonal velocity
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+               zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday )
+               DO jk = 1, jpkm1
+                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) ) ) * umask(ii,ij,jk)
+               END DO
+            END DO
+            !
+            igrd = 3                      ! Relaxation of meridional velocity
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+               zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday )
+               DO jk = 1, jpkm1
+                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) ) ) * vmask(ii,ij,jk)
+               END DO
+            END DO
+         ENDIF
+      ENDDO
+      !
+      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
+
+   END SUBROUTINE bdy_dyn3d_dmp
 
 #else

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

@@ -80 +80 @@
          &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
          &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &  
+         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              &
 #if defined key_lim2
          &             nn_ice_lim2, nn_ice_lim2_dta,                       &
@@ -121 +122 @@
       nn_tra(:)         = 0
       nn_tra_dta(:)     = -1  ! uninitialised flag
+      ln_tra_dmp(:)     = .false.
+      ln_dyn3d_dmp(:)   = .false.
+      rn_time_dmp(:)    = 1.
 #if defined key_lim2
       nn_ice_lim2(:)    = 0
@@ -135 +139 @@
       ! Check and write out namelist parameters
       ! -----------------------------------------
-
+ 
+      ln_tra_dmp_tot=.false.
+      ln_dyn3d_dmp_tot=.false.
       !                                   ! control prints
       IF(lwp) WRITE(numout,*) '         nambdy'
@@ -178 +184 @@
           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
+          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
         END SELECT
@@ -189 +197 @@
         IF(lwp) WRITE(numout,*)
 
+        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
+           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN
+              IF(lwp) WRITE(numout,*) 'No open boundary condition for the baroclinic velocitues : we force ln_dyn3d_dmp=.false.'
+              ln_dyn3d_dmp(ib_bdy)=.false.
+           ELSE
+              IF(lwp) WRITE(numout,*) '      Damping of the baroclinic velocities along the boundaries'
+              IF(lwp) WRITE(numout,*) '         Time damping :',rn_time_dmp(ib_bdy),' days'
+              ln_dyn3d_dmp_tot=.true.
+           ENDIF
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+
         IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
         SELECT CASE( nn_tra(ib_bdy) )                  
           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
+          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
+          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
         END SELECT
@@ -201 +224 @@
               CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
            END SELECT
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+
+        IF ( ln_tra_dmp(ib_bdy) ) THEN
+           IF ( nn_tra(ib_bdy).EQ.0 ) THEN
+              IF(lwp) WRITE(numout,*) 'No open boundary condition for the tracer : we force ln_tra_dmp=.false.'
+              ln_tra_dmp(ib_bdy)=.false.
+           ELSE
+              IF(lwp) WRITE(numout,*) '      Damping of the scalars along the boundaries'
+              IF(lwp) WRITE(numout,*) '         Time dumping :',rn_time_dmp(ib_bdy),' days'
+              ln_tra_dmp_tot=.true.
+           ENDIF
         ENDIF
         IF(lwp) WRITE(numout,*)
@@ -247 +282 @@
       ! ---------------------------------------------
       REWIND( numnam )                    
+      jpbdta = 1 
+      nblendta(:,:) = 0
       DO ib_bdy = 1, nb_bdy
 
-         jpbdta = 1
          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
  
@@ -282 +318 @@
             ENDIF
 
-            nblendta(:,ib_bdy) = 0
             DO iseg = 1, nbdysege
                igrd = 1
@@ -317 +352 @@
 
             nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
-            jpbdta = MAXVAL(nblendta(:,ib_bdy))               
-
 
          ELSE            ! Read size of arrays in boundary coordinates file.
@@ -324 +357 @@
 
             CALL iom_open( cn_coords_file(ib_bdy), inum )
-            jpbdta = 1
             DO igrd = 1, jpbgrd
                id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
                nblendta(igrd,ib_bdy) = kdimsz(1)
-               jpbdta = MAX(jpbdta, kdimsz(1))
             ENDDO
 
@@ -334 +365 @@
 
       ENDDO ! ib_bdy
+
+      jpbdta = MAXVAL(nblendta(:,:))
 
       ! Allocate arrays
@@ -610 +643 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
-               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
-!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic
-!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
+               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.      ! quadratic
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
             END DO
          END DO 
@@ -753 +786 @@
             END IF
          END DO
- 
+
          IF( icount /= 0 ) THEN
             IF(lwp) WRITE(numout,*)

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

@@ -8 +8 @@
    !!            3.0  !  2008-04  (NEMO team)  add in the reference version
    !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
-   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -15 +14 @@
    !!----------------------------------------------------------------------
    !!   PUBLIC
-   !!      tide_init     : read of namelist and initialisation of tidal harmonics data
+   !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data
    !!      tide_update   : calculation of tidal forcing at each timestep
-   !!   PRIVATE
-   !!      uvset         :\
-   !!      vday          : | Routines to correct tidal harmonics forcing for
-   !!      shpen         : | start time of integration
-   !!      ufset         : |
-   !!      vset          :/
    !!----------------------------------------------------------------------
    USE timing          ! Timing
@@ -34 +27 @@
    USE bdy_oce         ! ocean open boundary conditions
    USE daymod          ! calendar
+   USE tideini
+   USE tide_mod
    USE fldread, ONLY: fld_map
 
@@ -39 +34 @@
    PRIVATE
 
-   PUBLIC   tide_init     ! routine called in nemo_init
+   PUBLIC   bdytide_init  ! routine called in nemo_init
    PUBLIC   tide_update   ! routine called in bdydyn
 
    TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
-      INTEGER                                ::   ncpt       !: Actual number of tidal components
-      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.)
    END TYPE TIDES_DATA
 
-   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents
-
-   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data
+   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
    
    !!----------------------------------------------------------------------
@@ -61 +55 @@
 CONTAINS
 
-   SUBROUTINE tide_init
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE tide_init  ***
+   SUBROUTINE bdytide_init
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE bdytide_init  ***
       !!                     
       !! ** Purpose : - Read in namelist for tides and initialise external
@@ -71 +65 @@
       !! namelist variables
       !!-------------------
-      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date
-      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files
-      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used.
-      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
-      !!
-      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components
-      INTEGER                                   ::   ib_bdy, itide, ib  !: dummy loop indices
+      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
+      LOGICAL                                   ::   ln_conjug = .FALSE. !: F/T : tidal data in complex/complex conjugate
+      !!
+      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
       INTEGER                                   ::   inum, igrd
-      INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays)
-      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file 
-      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t           
-      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data
-      !!
-      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut   
-      !!
-      NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed
-      !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('tide_init')
+      INTEGER, DIMENSION(3)                     ::   ilen0               !: length of boundary data (from OBC arrays)
+      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
+      !!
+      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
+      !!
+      NAMELIST/nambdy_tide/filtide, ln_conjug
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
+      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~'
 
@@ -103 +92 @@
 
             ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
-            ln_tide_date = .false.
             filtide(:) = ''
-            tide_speed(:) = 0.0
-            tide_cpt(:) = ''
 
             ! Don't REWIND here - may need to read more than one of these namelists.
             READ  ( numnam, nambdy_tide )
-            !                                               ! Count number of components specified
-            td%ncpt = 0
-            DO itide = 1, jptides_max
-              IF( tide_cpt(itide) /= '' ) THEN
-                 td%ncpt = td%ncpt + 1
-              ENDIF
-            END DO
-
-            ! Fill in phase speeds from namelist
-            ALLOCATE( td%speed(td%ncpt) )
-            td%speed = tide_speed(1:td%ncpt)
-
-            ! Find constituents in standard list
-            DO itide = 1, td%ncpt
-               nindx(itide) = 0
-               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
-               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
-               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
-               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
-               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
-               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
-               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
-               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
-               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
-               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
-               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
-               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
-               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
-               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
-               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
-               IF( nindx(itide) == 0  .AND. lwp ) THEN
-                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
-                  CALL ctl_warn( ctmp1 )
-               ENDIF
-            END DO
             !                                               ! Parameter control and print
-            IF( td%ncpt < 1 ) THEN 
-               CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' )
-            ELSE
-               IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
-               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt
-               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt)
-               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
-               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt)
-            ENDIF
-
-
-            ! Allocate space for tidal harmonics data - 
-            ! get size from OBC data arrays
-            ! ---------------------------------------
+            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
+            IF(lwp) WRITE(numout,*) '             tidal components specified ', nb_harmo
+            IF(lwp) WRITE(numout,*) '                ', Wave(ntide(1:nb_harmo))%cname_tide
+            IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
+            IF(lwp) WRITE(numout,*) '                ', omega_tide(1:nb_harmo)
+
+            ! Allocate space for tidal harmonics data - get size from OBC data arrays
+            ! -----------------------------------------------------------------------
 
             ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 
-            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) )
+            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
+            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
 
             ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 
-            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) )
+            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
+            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
 
             ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 
-            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) )
+            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
+            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
 
             ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) )
-
 
             ! Open files and read in tidal forcing data
             ! -----------------------------------------
 
-            DO itide = 1, td%ncpt
+            DO itide = 1, nb_harmo
                !                                                              ! SSH fields
-               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
-               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               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%ssh(:,itide,1) = dta_read(1:ilen0(1),1,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%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1)
+               td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
                CALL iom_close( inum )
                !                                                              ! U fields
-               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
-               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               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%u(:,itide,1) = dta_read(1:ilen0(2),1,1)
+               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%u(:,itide,2) = dta_read(1:ilen0(2),1,1)
+               td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
                CALL iom_close( inum )
                !                                                              ! V fields
-               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
-               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               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%v(:,itide,1) = dta_read(1:ilen0(3),1,1)
+               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%v(:,itide,2) = dta_read(1:ilen0(3),1,1)
+               td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
                CALL iom_close( inum )
+               IF ( ln_conjug ) THEN
+                  IF(lwp) WRITE(numout,*) '       The tidal input data are written in complex conjugate'
+                  td%ssh0(:,itide,2) = - td%ssh0(:,itide,2)
+                  td%u0  (:,itide,2) = - td%u0  (:,itide,2)
+                  td%v0  (:,itide,2) = - td%v0  (:,itide,2)
+               ENDIF
                !
             END DO ! end loop on tidal components
-
-            IF( ln_tide_date ) THEN      ! correct for date factors
-
-!! used nmonth, nyear and nday from daymod....
-               ! Calculate date corrects for 15 standard consituents
-               ! This is the initialisation step, so nday, nmonth, nyear are the 
-               ! initial date/time of the integration.
-                 print *, nday,nmonth,nyear
-                 nyear  = int(ndate0 / 10000  )                          ! initial year
-                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month
-                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100)
-
-               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
-
-               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
-
-               DO itide = 1, td%ncpt       ! loop on tidal components
-                  !
-                  IF( nindx(itide) /= 0 ) THEN
-!!gm use rpi  and rad global variable  
-                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
-                     z_atde=z_ftc(nindx(itide))*cos(z_arg)
-                     z_btde=z_ftc(nindx(itide))*sin(z_arg)
-                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   &
-                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
-                  ELSE
-                     z_atde = 1.0_wp
-                     z_btde = 0.0_wp
-                  ENDIF
-                  !                                         !  elevation         
-                  igrd = 1
-                  DO ib = 1, ilen0(igrd)                
-                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2)
-                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1)
-                     td%ssh(ib,itide,1) = z1t
-                     td%ssh(ib,itide,2) = z2t
-                  END DO
-                  !                                         !  u       
-                  igrd = 2
-                  DO ib = 1, ilen0(igrd)                
-                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2)
-                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1)
-                     td%u(ib,itide,1) = z1t
-                     td%u(ib,itide,2) = z2t
-                  END DO
-                  !                                         !  v       
-                  igrd = 3
-                  DO ib = 1, ilen0(igrd)                
-                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2)
-                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1)
-                     td%v(ib,itide,1) = z1t
-                     td%v(ib,itide,2) = z2t
-                  END DO
-                  !
-               END DO   ! end loop on tidal components
-               !
-            ENDIF ! date correction
+            !
+            DEALLOCATE( dta_read )
             !
          ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
@@ -268 +163 @@
       END DO ! loop on ib_bdy
 
-      IF( nn_timing == 1 ) CALL timing_stop('tide_init')
-
-   END SUBROUTINE tide_init
-
+      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
+
+   END SUBROUTINE bdytide_init
 
    SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )
@@ -280 +174 @@
       !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter
-!!gm doctor jit ==> kit
-      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices
-      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data
-      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data
-      INTEGER,INTENT(in),OPTIONAL    ::   jit     ! Barotropic timestep counter (for timesplitting option)
-      INTEGER,INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
-                                                       ! is present then units = subcycle timesteps.
-                                                       ! time_offset = 0 => get data at "now" time level
-                                                       ! time_offset = -1 => get data at "before" time level
-                                                       ! time_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
-      REAL(wp)                         :: z_arg, z_sarg      
-      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
+      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,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
+                                                        ! is present then units = subcycle timesteps.
+                                                        ! time_offset = 0  => get data at "now"    time level
+                                                        ! time_offset = -1 => get data at "before" time level
+                                                        ! time_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
+      REAL(wp)                         :: z_arg, z_sarg, zflag      
+      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
       !!----------------------------------------------------------------------
 
       IF( nn_timing == 1 ) CALL timing_start('tide_update')
+
+      zflag=1
+      IF ( PRESENT(jit) ) THEN
+        IF ( jit /= 1 ) zflag=0
+      ENDIF
+
+      IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN
+        !
+        kt_tide = kt
+        !
+        IF(lwp) THEN
+           WRITE(numout,*)
+           WRITE(numout,*) 'bdy_tide : (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
@@ -306 +219 @@
       ENDIF
          
-      ! Note tide phase speeds are in deg/hour, so we need to convert the
-      ! elapsed time in seconds to hours by dividing by 3600.0
       IF( PRESENT(jit) ) THEN  
-         z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
+         z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) )
       ELSE                              
-         z_arg = (kt+time_add) * rdt * rad / 3600.0
+         z_arg = ((kt-kt_tide)+time_add) * rdt
       ENDIF
 
-      DO itide = 1, td%ncpt
-         z_sarg = z_arg * td%speed(itide)
+      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, td%ncpt
-         igrd=1                              ! SSH on tracer grid.
+      DO itide = 1, nb_harmo
+         igrd=1                              ! SSH on tracer grid
          DO ib = 1, idx%nblenrim(igrd)
             dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)
-            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2)
          END DO
          igrd=2                              ! U grid
          DO ib=1, idx%nblenrim(igrd)
-            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide)
-            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2)
+            dta%u2d(ib) = dta%u2d(ib) + 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, idx%nblenrim(igrd)
-            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide)
-            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2)
+            dta%v2d(ib) = dta%v2d(ib) + td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)
          END DO
       END DO
@@ -342 +250 @@
    END SUBROUTINE tide_update
 
-!!gm  doctor naming of dummy argument variables!!!   and all local variables
-
-   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
-      !!----------------------------------------------------------------------
-      !!                   ***  SUBROUTINE uvset  ***
-      !!                     
-      !! ** Purpose : - adjust tidal forcing for date factors
-      !!                
-      !!----------------------------------------------------------------------
-      implicit none
-      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
-      INTEGER, INTENT( in ) ::   iday    ! start time days
-      INTEGER, INTENT( in ) ::   imnth   ! start time month
-      INTEGER, INTENT( in ) ::   iyr     ! start time year
-      !!
-!!gm  nc is jptides_max ????
-      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents
-      CHARACTER(len=8), DIMENSION(nc) :: cname
-      INTEGER                 ::   year, vd, ivdy, ndc, i, k
-      REAL(wp)                ::   ss, h, p, en, p1, rtd
-      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction
-      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction
-      REAL(wp), DIMENSION(nc) ::   u, v, zig
-      !!
-      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   &
-         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   &
-         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      /
-      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   &
-         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   &
-         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   /
-      !!----------------------------------------------------------------------
-!
-!     ihs             -  start time gmt on ...
-!     iday/imnth/iyr  -  date   e.g.   12/10/87
-!
-      CALL vday(iday,imnth,iyr,ivdy)
-      vd = ivdy
-!
-!rp   note change of year number for d. blackman shpen
-!rp   if(iyr.ge.1000) year=iyr-1900
-!rp   if(iyr.lt.1000) year=iyr
-      year = iyr
-!
-!.....year  =  year of required data
-!.....vd    =  day of required data..set up for 0000gmt day year
-      ndc = nc
-!.....ndc   =  number of constituents allowed
-!!gm use rpi ?
-      rtd = 360.0 / 6.2831852
-      DO i = 1, ndc
-         zig(i) = zig(i)*rtd
-         ! sigo(i)= zig(i)
-      END DO
-
-!!gm try to avoid RETURN  in F90
-      IF( year == 0 )   RETURN
-      
-         CALL shpen( year, vd, ss, h , p , en, p1 )
-         CALL ufset( p   , en, u , f )
-         CALL vset ( ss  , h , p , en, p1, v )
-         !
-         DO k = 1, nc
-            z_vplu(k) = v(k) + u(k)
-            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
-            DO WHILE( z_vplu(k) < 0    )
-               z_vplu(k) = z_vplu(k) + 360.0
-            END DO
-            DO WHILE( z_vplu(k) > 360. )
-               z_vplu(k) = z_vplu(k) - 360.0
-            END DO
-         END DO
-         !
-      END SUBROUTINE uvset
-
-
-      SUBROUTINE vday( iday, imnth, iy, ivdy )
-      !!----------------------------------------------------------------------
-      !!                   *** SUBROUTINE vday  ***
-      !!                  
-      !! ** Purpose : - adjust tidal forcing for date factors
-      !!                
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ????
-      INTEGER, INTENT(  out) ::   ivdy              ! ???
-      !! 
-      INTEGER ::   iyr
-      !!------------------------------------------------------------------------------
-
-!!gm   nday_year in day mode is the variable compiuted here, no?
-!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year
-
-      !calculate day number in year from day/month/year
-       if(imnth.eq.1) ivdy=iday
-       if(imnth.eq.2) ivdy=iday+31
-       if(imnth.eq.3) ivdy=iday+59
-       if(imnth.eq.4) ivdy=iday+90
-       if(imnth.eq.5) ivdy=iday+120
-       if(imnth.eq.6) ivdy=iday+151
-       if(imnth.eq.7) ivdy=iday+181
-       if(imnth.eq.8) ivdy=iday+212
-       if(imnth.eq.9) ivdy=iday+243
-       if(imnth.eq.10) ivdy=iday+273
-       if(imnth.eq.11) ivdy=iday+304
-       if(imnth.eq.12) ivdy=iday+334
-       iyr=iy
-       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
-       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
-       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
-      !
-   END SUBROUTINE vday
-
-!!doctor norme for dummy arguments
-
-   SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE shpen  ***
-      !!                   
-      !! ** Purpose : - calculate astronomical arguments for tides
-      !!                this version from d. blackman 30 nove 1990
-      !!
-      !!----------------------------------------------------------------------
-!!gm add INTENT in, out or inout....    DOCTOR name....
-!!gm please do not use variable name with a single letter:  impossible to search in a code
-      INTEGER  ::   year,vd
-      REAL(wp) ::   s,h,p,en,p1
-      !!
-      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd
-      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
-      !!
-      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           &
-         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    &
-         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    &
-         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    &
-         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    &
-         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    &
-         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    &
-         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    &
-         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    &
-         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    &
-         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    &
-         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    &
-         &       52.57 /
-      !!----------------------------------------------------------------------
-
-      cycle = 360.0
-      ilc = 0
-      icent = year / 100
-      yr = year - icent * 100
-      t = icent - 20
-!
-!     for the following equations
-!     time origin is fixed at 00 hr of jan 1st,2000.
-!     see notes by cartwright
-!
-!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90)
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      it = icent - 20
-      if (it) 1,2,2
-    1 iday = it/4 -it
-      go to 3
-    2 iday = (it+3)/4 - it
-!
-!     t is in julian century
-!     correction in gegorian calander where only century year divisible
-!     by 4 is leap year.
-!
-    3 continue
-!
-      td = 0.0
-!
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      if (yr) 4,5,4
-!
-    4 iyd = 365*yr
-      ild = (yr-1)/4
-      if((icent - (icent/4)*4) .eq. 0) ilc = 1
-      td = iyd + ild + ilc
-!
-    5 td = td + iday + vd -1.0 - 0.5
-      t = t + (td/36525.0)
-!
-      ipos=year-1899
-      if (ipos .lt. 0) go to 7
-      if (ipos .gt. 83) go to 6
-!
-      delta = (delt(ipos+1)+delt(ipos))/2.0
-      go to 7
-!
-    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
-!
-    7 deltat = delta * 1.0e-6
-!
-!!gm   precision of the computation   :  example for s it should be replace by:
-!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results
-      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
-      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat
-      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat
-      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat
-      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t
-      !
-      nn = s / cycle
-      s = s - nn * cycle
-      IF( s < 0.e0 )   s = s + cycle
-      !
-      nn = h / cycle
-      h = h - cycle * nn
-      IF( h < 0.e0 )   h = h + cycle
-      !
-      nn = p / cycle
-      p = p - cycle * nn
-      IF( p < 0.e0)   p = p + cycle
-      !
-      nn = en / cycle
-      en = en - cycle * nn
-      IF( en < 0.e0 )   en = en + cycle
-      en = cycle - en
-      !
-      nn = p1 / cycle
-      p1 = p1 - nn * cycle
-      !
-   END SUBROUTINE shpen
-
-
-   SUBROUTINE ufset( p, cn, b, a )
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE ufset  ***
-      !!                    
-      !! ** Purpose : - calculate nodal parameters for the tides
-      !!                
-      !!----------------------------------------------------------------------
-!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!gm  nc is jptides_max ????
-      integer nc
-      parameter (nc=15)
-      REAL(wp) p,cn
-      !!
-      
-!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z"
-      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
-      REAL(wp) ::   a(nc), b(nc)
-      INTEGER  ::   ndc, k
-      !!----------------------------------------------------------------------
-
-      ndc = nc
-
-!    a=f       ,  b =u
-!    t is  zero as compared to tifa.
-!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module
-      rad = 6.2831852d0/360.0
-      pw = p  * rad
-      nw = cn * rad
-      w1 = cos(   nw )
-      w2 = cos( 2*nw )
-      w3 = cos( 3*nw )
-      w4 = sin(   nw )
-      w5 = sin( 2*nw )
-      w6 = sin( 3*nw )
-      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   &
-         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw )
-      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   &
-         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw )
-      !
-      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
-      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
-      !   q1
-      a(2) = a(1)
-      b(2) = b(1)
-      !   o1
-      a(3) = 1.0
-      b(3) = 0.0
-      !   p1
-      a(4) = 1.0
-      b(4) = 0.0
-      !   s1
-      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
-      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
-      !   k1
-      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
-      b(6) = -0.0374*w4
-      !  2n2
-      a(7) = a(6)
-      b(7) = b(6)
-      !  mu2
-      a(8) = a(6)
-      b(8) = b(6)
-      !   n2
-      a(9) = a(6)
-      b(9) = b(6)
-      !  nu2
-      a(10) = a(6)
-      b(10) = b(6)
-      !   m2
-      a(11) = SQRT( w7 * w7 + w8 * w8 )
-      b(11) = ATAN( w8 / w7 )
-!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ????
-      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992
-      !   l2
-      a(12) = 1.0
-      b(12) = 0.0
-      !   t2
-      a(13)= a(12)
-      b(13)= b(12)
-      !   s2
-      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
-      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
-      !   k2
-      a(15) = a(6)*a(6)
-      b(15) = 2*b(6)
-      !   m4
-!!gm  old coding,  remove GOTO and label of lines
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      DO 40 k = 1,ndc
-         b(k) = b(k)/rad
-32       if (b(k)) 34,35,35
-34       b(k) = b(k) + 360.0
-         go to 32
-35       if (b(k)-360.0) 40,37,37
-37       b(k) = b(k)-360.0
-         go to 35
-40    continue
-      !
-   END SUBROUTINE ufset
-   
-
-   SUBROUTINE vset( s,h,p,en,p1,v)
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE vset  ***
-      !!                   
-      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
-      !!                
-      !!----------------------------------------------------------------------
-!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!gm  nc is jptides_max ????
-!!gm  en argument is not used: suppress it ?
-      integer nc
-      parameter (nc=15)
-      real(wp) s,h,p,en,p1
-      real(wp)   v(nc)
-      !!
-      integer ndc, k
-      !!----------------------------------------------------------------------
-
-      ndc = nc
-      !   v s  are computed here.
-      v(1) =-3*s +h +p +270      ! Q1
-      v(2) =-2*s +h +270.0       ! O1
-      v(3) =-h +270              ! P1
-      v(4) =180                  ! S1
-      v(5) =h +90.0              ! K1
-      v(6) =-4*s +2*h +2*p       ! 2N2
-      v(7) =-4*(s-h)             ! MU2
-      v(8) =-3*s +2*h +p         ! N2
-      v(9) =-3*s +4*h -p         ! MU2
-      v(10) =-2*s +2*h           ! M2
-      v(11) =-s +2*h -p +180     ! L2 
-      v(12) =-h +p1              ! T2
-      v(13) =0                   ! S2
-      v(14) =h+h                 ! K2
-      v(15) =2*v(10)             ! M4
-!
-!!gm  old coding,  remove GOTO and label of lines
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      do 72 k = 1, ndc
-69       if( v(k) )   70,71,71
-70       v(k) = v(k)+360.0
-         go to 69
-71       if( v(k) - 360.0 )   72,73,73
-73       v(k) = v(k)-360.0
-         go to 71
-72    continue
-      !
-   END SUBROUTINE vset
-
+   SUBROUTINE tide_init_elevation( idx, td )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init_elevation  ***
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
+      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data
+      !! * Local declarations
+      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
+      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
+
+      igrd=1                                 ! SSH on tracer grid.
+
+      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
+
+      DO itide = 1, nb_harmo
+         DO ib = 1, idx%nblenrim(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))
+         END DO
+         DO ib = 1, idx%nblenrim(igrd)
+            mod_tide(ib)=mod_tide(ib)*ftide(itide)
+            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+         ENDDO
+         DO ib = 1, idx%nblenrim(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)
+
+   END SUBROUTINE tide_init_elevation
+
+   SUBROUTINE tide_init_velocities( idx, td )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init_elevation  ***
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
+      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
+      !! * Local declarations
+      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
+      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
+
+      igrd=2                                 ! U grid.
+
+      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
+
+      DO itide = 1, nb_harmo
+         DO ib = 1, idx%nblenrim(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))
+         END DO
+         DO ib = 1, idx%nblenrim(igrd)
+            mod_tide(ib)=mod_tide(ib)*ftide(itide)
+            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+         ENDDO
+         DO ib = 1, idx%nblenrim(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                                 ! U grid.
+
+      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
+
+      DO itide = 1, nb_harmo
+         DO ib = 1, idx%nblenrim(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))
+         END DO
+         DO ib = 1, idx%nblenrim(igrd)
+            mod_tide(ib)=mod_tide(ib)*ftide(itide)
+            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+         ENDDO
+         DO ib = 1, idx%nblenrim(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
 #else
    !!----------------------------------------------------------------------
    !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
    !!----------------------------------------------------------------------
-!!gm  are you sure we need to define filtide and tide_cpt ?
-   CHARACTER(len=80), PUBLIC               ::   filtide                !: Filename root for tidal input files
-   CHARACTER(len=4 ), PUBLIC, DIMENSION(1) ::   tide_cpt               !: Names of tidal components used.
-
 CONTAINS
-   SUBROUTINE tide_init                ! Empty routine
-   END SUBROUTINE tide_init
+   SUBROUTINE bdytide_init                ! Empty routine
+   END SUBROUTINE bdytide_init
    SUBROUTINE tide_data                ! Empty routine
    END SUBROUTINE tide_data
-   SUBROUTINE tide_update( kt, kit )   ! Empty routine
-      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit
+   SUBROUTINE tide_update( kt, jit )   ! Empty routine
+      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit
    END SUBROUTINE tide_update
 #endif

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

@@ -22 +22 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
+   USE phycst
 
    IMPLICIT NONE
@@ -27 +28 @@
 
    PUBLIC bdy_tra      ! routine called in tranxt.F90 
+   PUBLIC bdy_tra_dmp  ! routine called in tranxt.F90 
 
    !!----------------------------------------------------------------------
@@ -53 +55 @@
          CASE(jp_frs)
             CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(2)
+            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(3)
+            CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(4)
+            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
          CASE DEFAULT
             CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
@@ -97 +105 @@
       !
    END SUBROUTINE bdy_tra_frs
-   
+  
+   SUBROUTINE bdy_tra_spe( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_frs  ***
+      !!                    
+      !! ** Purpose : Apply a specified value for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_spe')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            tsa(ii,ij,ik,jp_tem) = dta%tem(ib,ik) * tmask(ii,ij,ik)
+            tsa(ii,ij,ik,jp_sal) = dta%sal(ib,ik) * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_spe')
+      !
+   END SUBROUTINE bdy_tra_spe
+
+   SUBROUTINE bdy_tra_nmn( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_nmn  ***
+      !!                    
+      !! ** Purpose : Duplicate the value for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_nmn')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ! search the sense of the gradient
+            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
+            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
+            IF ( zcoef1+zcoef2 == 0) THEN
+               ! corner
+               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
+               tsa(ii,ij,ik,jp_tem) = tsa(ii-1,ij  ,ik,jp_tem) * tmask(ii-1,ij  ,ik) + &
+                 &                    tsa(ii+1,ij  ,ik,jp_tem) * tmask(ii+1,ij  ,ik) + &
+                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + &
+                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik)
+               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + &
+                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + &
+                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + &
+                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik)
+               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+            ELSE
+               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik)
+               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik)
+            ENDIF
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_nmn')
+      !
+   END SUBROUTINE bdy_tra_nmn
+
+   SUBROUTINE bdy_tra_rnf( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_rnf  ***
+      !!                    
+      !! ** Purpose : Apply the runoff values for tracers at open boundaries:
+      !!                  - specified to 0.1 PSU for the salinity
+      !!                  - duplicate the value for the temperature
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_rnf')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+            jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+            tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii,ij,ik)
+            tsa(ii,ij,ik,jp_sal) =                        0.1 * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_rnf')
+      !
+   END SUBROUTINE bdy_tra_rnf
+
+   SUBROUTINE bdy_tra_dmp( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_dmp  ***
+      !!                    
+      !! ** Purpose : Apply damping for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      REAL(wp) ::   zta, zsa, ztime
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      INTEGER  ::   ib_bdy         ! Loop index
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_dmp')
+      !
+      DO ib_bdy=1, nb_bdy
+         IF ( ln_tra_dmp(ib_bdy) ) THEN
+            igrd = 1                       ! Everything is at T-points here
+            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+               zwgt = idx_bdy(ib_bdy)%nbw(ib,igrd) / ( rn_time_dmp(ib_bdy) * rday )
+               DO ik = 1, jpkm1
+                  tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik)
+                  tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik)
+               END DO
+            END DO
+         ENDIF
+      ENDDO
+      !
+      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )  ;  CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp')
+      !
+   END SUBROUTINE bdy_tra_dmp
+ 
 #else
    !!----------------------------------------------------------------------

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -403 +403 @@
          IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
          IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
-         IF ( ln_tide_pot ) CALL upd_tide( kt, jn )
+         IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn )
 
          !                                                !* after ssh_e
@@ -453 +453 @@
                   ENDIF
                   ! add tidal astronomical forcing
-                  IF ( ln_tide_pot ) THEN 
+                  IF ( ln_tide_pot .AND. lk_tide ) THEN 
                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
@@ -503 +503 @@
                   ENDIF
                   ! add tidal astronomical forcing
-                  IF ( ln_tide_pot ) THEN
+                  IF ( ln_tide_pot .AND. lk_tide ) THEN
                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
@@ -550 +550 @@
                   ENDIF
                   ! add tidal astronomical forcing
-                  IF ( ln_tide_pot ) THEN
+                  IF ( ln_tide_pot .AND. lk_tide ) THEN
                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90

@@ -66 +66 @@
       !!
       INTEGER            ::   ierror  ! local integer 
-      REAL(wp)           ::   zpref   ! local scalar
       !!
       CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files
       TYPE(FLD_N)        ::  sn_apr   ! informations about the fields to be read
       !!
-      NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr
+      NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr, rpref, ln_apr_obc
       !!----------------------------------------------------------------------
       !
@@ -113 +112 @@
          !
          !                                            !* control check
-         IF( ln_apr_obc  )   &
-            CALL ctl_stop( 'sbc_apr: inverse barometer added to OBC ssh data not yet implemented ' )
-         IF( ln_apr_obc .AND. .NOT. lk_obc )   &
-            CALL ctl_stop( 'sbc_apr: add inverse barometer to OBC requires to use key_obc' )
+         IF ( ln_apr_obc  ) THEN
+            IF(lwp) WRITE(numout,*) '         Inverse barometer added to OBC ssh data'
+         ENDIF
          IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts )   &
             CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' )

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90

@@ -14 +14 @@
   USE daymod
   USE dynspg_oce
-  USE tide_mod
+  USE tideini
   USE iom
 
@@ -21 +21 @@
 
   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro
-  LOGICAL, PUBLIC :: ln_tide_pot = .false.
+
 #if defined key_tide
 
   LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE.
-
-  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide 
-
-  REAL(wp), ALLOCATABLE, DIMENSION(:) ::  &
-       v0tide,      &
-       utide,       &
-       ftide
-
   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot,phi_pot
-
-  INTEGER, PUBLIC :: nb_harmo
-  INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide
-  INTEGER, PUBLIC :: nn_tide, kt_tide
-
   !!---------------------------------------------------------------------------------
   !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -51 +38 @@
     !! * Arguments
     INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-    !! * Local declarations
-    INTEGER  :: jk,ji
-    CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
     !!----------------------------------------------------------------------
 
-    NAMELIST/nam_tide/ln_tide_pot,nb_harmo,clname,nn_tide
+    IF ( kt == nit000 .AND. .NOT. lk_dynspg_ts )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' )
 
-    IF ( kt == nit000 ) THEN
+    IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN
+      !
+      kt_tide = kt
 
-       IF( .NOT. lk_dynspg_ts  )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' )
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_tide : (re)Initialization of the tidal potential at kt=',kt
+         WRITE(numout,*) '~~~~~~~ '
+      ENDIF
 
-    ! Read Namelist nam_tide
+      IF(lwp) THEN
+         IF ( kt == nit000 ) WRITE(numout,*) 'Apply astronomical potential : ln_tide_pot =', ln_tide_pot
+         CALL flush(numout)
+      ENDIF
 
-    nn_tide=INT(rday/rdt)
+      IF ( kt == nit000 ) ALLOCATE(amp_pot(jpi,jpj,nb_harmo))
+      IF ( kt == nit000 ) ALLOCATE(phi_pot(jpi,jpj,nb_harmo))
+      IF ( kt == nit000 ) ALLOCATE(pot_astro(jpi,jpj))
 
-    CALL tide_init_Wave
+      amp_pot(:,:,:) = 0.e0
+      phi_pot(:,:,:) = 0.e0
+      pot_astro(:,:) = 0.e0
 
-    REWIND ( numnam )
-    READ   ( numnam, nam_tide )
-
-    IF(lwp) THEN
-       WRITE(numout,*)
-       WRITE(numout,*) 'sbc_tide : Initialization of the tidal components'
-       WRITE(numout,*) '~~~~~~~ '
+      IF ( ln_tide_pot ) CALL tide_init_potential
+      !
     ENDIF
-
-    IF(lwp) THEN
-       WRITE(numout,*) '        Namelist nam_tide'
-       WRITE(numout,*) '        Apply astronomical potential : ln_tide_pot =', ln_tide_pot
-       WRITE(numout,*) '        nb_harmo    = ', nb_harmo
-       CALL flush(numout)
-    ENDIF
-
-    ALLOCATE(ntide     (nb_harmo))
-    DO jk=1,nb_harmo
-       DO ji=1,jpmax_harmo
-          IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN
-             ntide(jk) = ji
-             EXIT
-          END IF
-       END DO
-    END DO
-    ALLOCATE(omega_tide(nb_harmo))
-    ALLOCATE(v0tide    (nb_harmo))
-    ALLOCATE(utide     (nb_harmo))
-    ALLOCATE(ftide     (nb_harmo))
-    ALLOCATE(amp_pot(jpi,jpj,nb_harmo))
-    ALLOCATE(phi_pot(jpi,jpj,nb_harmo))
-    ALLOCATE(pot_astro(jpi,jpj))
-    ENDIF
-
-    IF ( MOD( kt - 1, nn_tide ) == 0 ) THEN
-      kt_tide = kt
-      CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo)
-    ENDIF
-
-    amp_pot(:,:,:) = 0.e0
-    phi_pot(:,:,:) = 0.e0
-    pot_astro(:,:) = 0.e0
-
-    IF (ln_tide_pot          ) CALL tide_init_potential
 
   END SUBROUTINE sbc_tide

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -46 +46 @@
    USE mppini          ! shared/distributed memory setting (mpp_init routine)
    USE domain          ! domain initialization             (dom_init routine)
+   USE tideini         ! tidal components initialization   (tide_ini routine)
    USE obcini          ! open boundary cond. initialization (obc_ini routine)
    USE bdyini          ! open boundary cond. initialization (bdy_init routine)
    USE bdydta          ! open boundary cond. initialization (bdy_dta_init routine)
-   USE bdytides        ! open boundary cond. initialization (tide_init routine)
+   USE bdytides        ! open boundary cond. initialization (bdytide_init routine)
    USE istate          ! initial state setting          (istate_init routine)
    USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
@@ -314 +315 @@
 
       IF( lk_obc        )   CALL     obc_init   ! Open boundaries 
+
+                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
+
+                            CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics
+
       IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation
       IF( lk_bdy        )   CALL     bdy_dta_init   ! Open boundaries initialisation of external data arrays
-      IF( lk_bdy        )   CALL     tide_init      ! Open boundaries initialisation of tidal harmonic forcing
+      IF( lk_bdy        )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing
 
                             CALL flush(numout)
                             CALL dyn_nept_init  ! simplified form of Neptune effect
                             CALL flush(numout)
-
-                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
       !                                     ! Ocean physics

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -95 +95 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+      IF( kstp /= nit000 )   CALL tide_init ( kstp )
       IF( lk_tide    )   CALL sbc_tide( kstp )
       IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
@@ -187 +188 @@
       IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
       IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
+      IF( lk_bdy .AND.     &
+         & ln_tra_dmp_tot)   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
       IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
@@ -219 +222 @@
          & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
       IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
+      IF( lk_bdy .AND.     &
+         & ln_dyn3d_dmp_tot)   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
                                CALL dyn_adv( kstp )         ! advection (vector or flux form)
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis

branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

@@ -55 +55 @@
 
    USE bdy_par          ! for lk_bdy
+   USE bdy_oce          ! for dmp logical
    USE bdydta           ! open boundary condition data     (bdy_dta routine)
+   USE bdytra           ! bdy cond. for tracers            (bdy_tra routine)
+   USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine)
 
    USE sshwzv           ! vertical velocity and ssh        (ssh_wzv routine)