Changeset 367

Timestamp:

2005-12-28T10:25:10+01:00 (18 years ago)

Author:

opalod

Message:

nemo_v1_update_035 : CT : OBCs adapted to the new surface pressure gradient algorithms

Location:

trunk/NEMO/OPA_SRC

Files:

• DOM/dommsk.F90 (modified) (1 diff)
• DYN/dynnxt.F90 (modified) (2 diffs)
• DYN/dynspg.F90 (modified) (3 diffs)
• DYN/dynspg_exp.F90 (modified) (3 diffs)
• DYN/dynspg_oce.F90 (added)
• DYN/dynspg_ts.F90 (modified) (20 diffs)
• OBC/obc_oce.F90 (modified) (9 diffs)
• OBC/obccli.F90 (modified) (1 diff)
• OBC/obcdta.F90 (modified) (21 diffs)
• OBC/obcdyn.F90 (modified) (20 diffs)
• OBC/obcdyn_bt.F90 (added)
• OBC/obcfla.F90 (added)
• OBC/obcini.F90 (modified) (8 diffs)
• OBC/obcrad.F90 (modified) (8 diffs)
• OBC/obcrst.F90 (modified) (20 diffs)
• OBC/obcvol.F90 (modified) (3 diffs)
• SOL/solver.F90 (modified) (1 diff)
• TRA/traadv_cen2.F90 (modified) (1 diff)
• TRA/traadv_muscl.F90 (modified) (1 diff)
• TRA/traadv_muscl2.F90 (modified) (1 diff)
• TRA/traadv_tvd.F90 (modified) (1 diff)
• opa.F90 (modified) (1 diff)
• restart.F90 (modified) (1 diff)
• step.F90 (modified) (1 diff)
• stpctl.F90 (modified) (1 diff)

trunk/NEMO/OPA_SRC/DOM/dommsk.F90

@@ -18 +18 @@
    USE lib_mpp
    USE solisl          ! ???
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
 
    IMPLICIT NONE

trunk/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -12 +12 @@
    USE dom_oce         ! ocean space and time domain
    USE in_out_manager  ! I/O manager
+   USE obc_oce         ! ocean open boundary conditions
    USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
+   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
+   USE obcvol          ! ocean open boundary condition (obc_vol routines)
+   USE dynspg_oce      ! type of surface pressure gradient
    USE lbclnk          ! lateral boundary condition (or mpp link)
    USE prtctl          ! Print control
@@ -105 +109 @@
       ! Update (ua,va) along open boundaries (only in the rigid-lid case)
       CALL obc_dyn( kt )
+
+      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
+         !Flather boundary condition :
+         !        - Update sea surface height on each open boundary
+         !                 sshn (= after ssh) for explicit case
+         !                 sshn_b (= after ssha_b) for time-splitting case
+         !        - Correct the barotropic velocities
+         CALL obc_dyn_bt( kt )
+
+         !Boundary conditions on sshn ( after ssh)
+         CALL lbc_lnk( sshn, 'T', 1. )
+
+         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+            CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
+         ENDIF
+
+         IF ( ln_vol_cst ) CALL obc_vol( kt )
+
+      ENDIF
+
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab

trunk/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -12 +12 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
+   USE obc_oce        ! ocean open boundary conditions
+   USE dynspg_oce     ! surface pressure gradient variables
    USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
@@ -27 +29 @@
    !! * Accessibility
    PUBLIC dyn_spg         ! routine called by step module
-
-   !! * Public variables
-#if   defined key_dynspg_exp   ||  defined key_esopa
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .TRUE.  !: Explicit free surface flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .FALSE. !: Explicit free surface flag
-#endif
-#if   defined key_dynspg_ts   ||  defined key_esopa
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .TRUE.  !: Free surface with time splitting flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .FALSE. !: Free surface with time splitting flag
-#endif
-#if   defined key_dynspg_flt  ||  defined key_esopa
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .TRUE.  !: Filtered free surface cst volume flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .FALSE. !: Filtered free surface cst volume flag
-#endif
-#if   defined key_dynspg_rl
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl  = .TRUE.  !: Rigid-lid flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl  = .FALSE. !: Rigid-lid flag
-#endif
 
    !! * module variables
@@ -208 +188 @@
       ENDIF
 
+#if key_obc
+      ! Conservation of ocean volume (key_dynspg_flt)
+      ! ---------------------------------------------
+      IF( lk_dynspg_flt ) ln_vol_cst = .true.
+
+      ! Application of Flather's algorithm at open boundaries
+      ! -----------------------------------------------------
+      IF( lk_dynspg_flt ) ln_obc_fla = .false.
+      IF( lk_dynspg_exp ) ln_obc_fla = .true.
+      IF( lk_dynspg_ts  ) ln_obc_fla = .true.
+#endif
+
    END SUBROUTINE dyn_spg_ctl
 

trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90

@@ -20 +20 @@
    USE ocesbc          ! ocean surface boundary condition
    USE obc_oce         ! Lateral open boundary condition
+   USE obcdta          ! open boundary condition data     (obc_dta_bt routine)
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -101 +102 @@
       ENDIF
 
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! 0. Initialization
+      ! -----------------
+      ! read or estimate sea surface height and vertically integrated velocities
+      IF( lk_obc )   CALL obc_dta_bt( kt, 0 )
       z2dt = 2. * rdt                                       ! time step: leap-frog
       IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
@@ -179 +182 @@
       ENDIF
 
-      !Boundary conditions on sshn
-      CALL lbc_lnk( sshn, 'T', 1. )
+      ! Boundary conditions on sshn
+      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
  
       IF(ln_ctl) THEN         ! print sum trends (used for debugging)

trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -17 +17 @@
    USE phycst          ! physical constants
    USE ocesbc          ! ocean surface boundary condition
+   USE obcdta          ! open boundary condition data     
+   USE obcfla          ! Flather open boundary condition  
    USE dynvor          ! vorticity term
    USE obc_oce         ! Lateral open boundary condition
@@ -22 +24 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE dynspg_oce      ! surface pressure gradient variables
    USE in_out_manager  ! I/O manager
 
@@ -29 +32 @@
    !! * Accessibility
    PUBLIC dyn_spg_ts  ! routine called by step.F90
-
-   !! * Module variables
-      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop
-         sshn_b, sshb_b,               &  ! sea surface heigth (now, before)
-         un_b  , vn_b                     ! vertically integrated horizontal velocities (now)
 
    !! * Substitutions
@@ -61 +59 @@
       !!          surface gradient and the Coriolis force are updated within
       !!          the barotropic integration.
-      !!      -2- Barotropic loop : updates of sea surface height (zssha_e) and 
-      !!          barotropic transports (zua_e and zva_e) through barotropic 
+      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and 
+      !!          barotropic transports (ua_e and va_e) through barotropic 
       !!          momentum and continuity integration. Barotropic former 
       !!          variables are time averaging over the full barotropic cycle
@@ -97 +95 @@
          zssha_b, zua_b, zva_b,             &  !     "        "
          zsshb_e, zub_e, zvb_e,             &  !     "        "
-         zsshn_e, zun_e, zvn_e,             &  !     "        "
-         zssha_e, zua_e, zva_e                 !     "        "
+         zun_e, zvn_e                          !     "        "
       REAL(wp), DIMENSION(jpi,jpj),SAVE ::  &
          ztnw, ztne, ztsw, ztse
@@ -105 +102 @@
       ! Arrays initialization
       ! ---------------------
-      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0   ;   zua_e(:,:) = 0.e0
-      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0   ;   zva_e(:,:) = 0.e0
+      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0
+      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0
       zhdiv(:,:) = 0.e0
 
@@ -138 +135 @@
             ENDIF
          ENDIF
-         zssha_e(:,:) = sshn(:,:)
-         zua_e  (:,:) = un_b(:,:)
-         zva_e  (:,:) = vn_b(:,:)
+         ssha_e(:,:) = sshn(:,:)
+         ua_e(:,:)   = un_b(:,:)
+         va_e(:,:)   = vn_b(:,:)
 
          IF( ln_dynvor_een ) THEN
@@ -278 +275 @@
       ! variables for the barotropic equations
       zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
-      zsshn_e(:,:) = sshn_b(:,:)
-      zub_e(:,:) = un_b(:,:)       ! barotropic transports issued from the barotropic equations (before and now)
-      zvb_e(:,:) = vn_b(:,:)
-      zun_e(:,:) = un_b(:,:)
-      zvn_e(:,:) = vn_b(:,:)
-      zssha_b(:,:) = sshn(:,:)        ! time averaged variables over all sub-timesteps
-      zua_b(:,:) = un_b(:,:)   
-      zva_b(:,:) = vn_b(:,:)
+      sshn_e (:,:) = sshn_b(:,:)
+      zub_e  (:,:) = un_b  (:,:)       ! barotropic transports issued from the barotropic equations (before and now)
+      zvb_e  (:,:) = vn_b  (:,:)
+      zun_e  (:,:) = un_b  (:,:)
+      zvn_e  (:,:) = vn_b  (:,:)
+      zssha_b(:,:) = sshn  (:,:)       ! time averaged variables over all sub-timesteps
+      zua_b  (:,:) = un_b  (:,:)   
+      zva_b  (:,:) = vn_b  (:,:)
+
+      ! set ssh corrections to 0
+      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
+#if defined key_obc
+      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
+      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
+      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
+      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
+#endif
 
       ! Barotropic integration over 2 baroclinic time steps
@@ -296 +302 @@
          z2dt_e = 2. * rdtbt
          IF ( jit == 1 )   z2dt_e = rdtbt
+
+         ! Time interpolation of open boundary condition data
+         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
 
          ! Horizontal divergence of barotropic transports
@@ -312 +321 @@
          ! open boundaries (div must be zero behind the open boundary)
          !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
-         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
-         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
-         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
-         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
+         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0      ! east
+         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0      ! west
+         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
+         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0      ! south
 #endif
 
@@ -322 +331 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
+               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
             &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
             END DO
@@ -336 +345 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
+                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
                   ! energy conserving formulation for planetary vorticity term
                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
@@ -346 +355 @@
                   zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
                   ! after transports
-                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
@@ -355 +364 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
+                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
                   ! enstrophy conserving formulation for planetary vorticity term
                   zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
@@ -365 +374 @@
                   zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
                   ! after transports
-                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
@@ -375 +384 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
+                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
                   ! energy/enstrophy conserving formulation for planetary vorticity term
                   zcubt = + zfac25 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
@@ -383 +392 @@
                      &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
                   ! after transports
-                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
          ENDIF
 
-         ! ... Boundary conditions on zua_e, zva_e, zssha_e
-         CALL lbc_lnk( zua_e, 'U', -1. )
-         CALL lbc_lnk( zva_e, 'V', -1. )
-         CALL lbc_lnk( zssha_e, 'T', 1. )
+         ! Flather's boundary condition for the barotropic loop :
+         !         - Update sea surface height on each open boundary
+         !         - Correct the barotropic transports
+         IF( lk_obc )   CALL obc_fla_ts( kt )
+
+
+         ! ... Boundary conditions on ua_e, va_e, ssha_e
+         CALL lbc_lnk( ua_e  , 'U', -1. )
+         CALL lbc_lnk( va_e  , 'V', -1. )
+         CALL lbc_lnk( ssha_e, 'T',  1. )
 
          ! temporal sum
          !-------------
-         zssha_b(:,:) = zssha_b(:,:) + zssha_e(:,:)
-         zua_b  (:,:) = zua_b  (:,:) + zua_e  (:,:)
-         zva_b  (:,:) = zva_b  (:,:) + zva_e  (:,:) 
+         zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:)
+         zua_b  (:,:) = zua_b  (:,:) + ua_e  (:,:)
+         zva_b  (:,:) = zva_b  (:,:) + va_e  (:,:) 
 
          ! Time filter and swap of dynamics arrays
          ! ---------------------------------------
          IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
-            zsshb_e(:,:) = zsshn_e(:,:)
-            zub_e  (:,:) = zun_e  (:,:)
-            zvb_e  (:,:) = zvn_e  (:,:)
-            zsshn_e(:,:) = zssha_e(:,:)
-            zun_e  (:,:) = zua_e  (:,:)
-            zvn_e  (:,:) = zva_e  (:,:)
+            zsshb_e(:,:) = sshn_e(:,:)
+            zub_e  (:,:) = zun_e (:,:)
+            zvb_e  (:,:) = zvn_e (:,:)
+            sshn_e (:,:) = ssha_e(:,:)
+            zun_e  (:,:) = ua_e  (:,:)
+            zvn_e  (:,:) = va_e  (:,:)
          ELSE                                        ! Asselin filtering
-            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * zsshn_e(:,:)
-            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * zun_e  (:,:)
-            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * zvn_e  (:,:)
-            zsshn_e(:,:) = zssha_e(:,:)
-            zun_e  (:,:) = zua_e  (:,:)
-            zvn_e  (:,:) = zva_e  (:,:)
+            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
+            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:)
+            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:)
+            sshn_e (:,:) = ssha_e(:,:)
+            zun_e  (:,:) = ua_e  (:,:)
+            zvn_e  (:,:) = va_e  (:,:)
          ENDIF
 
@@ -426 +441 @@
       zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
       zssha_b(:,:) = zcoef * zssha_b(:,:) 
-      zua_b  (:,:) = zcoef * zua_b  (:,:) 
-      zva_b  (:,:) = zcoef * zva_b  (:,:) 
+      zua_b  (:,:) = zcoef *  zua_b (:,:) 
+      zva_b  (:,:) = zcoef *  zva_b (:,:) 
+#if defined key_obc
+         IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
+         IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
+         IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
+         IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
+#endif
      
 
@@ -448 +469 @@
       ! open boundaries (div must be zero behind the open boundary)
       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
-      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
-      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
+      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
+      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
-      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
+      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
 #endif
 
@@ -460 +481 @@
 
       ! ... Boundary conditions on sshn
-      CALL lbc_lnk( sshn, 'T', 1. )
+      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
 
 

trunk/NEMO/OPA_SRC/OBC/obc_oce.F90

@@ -40 +40 @@
    !!--------------------------------------
    INTEGER ::              & !: * namelist ??? *
-      nbobc    = 1  ,      & !: number of open boundaries ( 1=< nbobc =< 4 ) 
-      nobc_dta = 0  ,      & !:  = 0 use the initial state as obc data
+      nbobc    = 2  ,      & !: number of open boundaries ( 1=< nbobc =< 4 ) 
+      nobc_dta = 0           !:  = 0 use the initial state as obc data
       !                      !   = 1 read obc data in obcxxx.dta files
 
    LOGICAL ::  ln_obc_clim = .true.  !:  obc data files are climatological
+   LOGICAL ::  ln_obc_fla  = .false. !:  Flather open boundary condition not used
+   LOGICAL ::  ln_vol_cst  = .true.  !:  Conservation of the whole volume
 
    REAL(wp) ::             & !!: open boundary namelist (namobc)
@@ -117 +119 @@
 
    REAL(wp), DIMENSION(jpjed:jpjef) ::   &  !:
-      bfoe                !: now climatology of the east boundary barotropic stream function 
+      bfoe,             & !: now climatology of the east boundary barotropic stream function 
+      sshfoe,           & !: now climatology of the east boundary sea surface height
+      ubtfoe,vbtfoe       !: now climatology of the east boundary barotropic transport
      
    REAL(wp), DIMENSION(jpj,jpk) ::   &  !:
@@ -124 +128 @@
       uclie               !: baroclinic componant of the zonal velocity after radiation 
       !                   ! in the obcdyn.F90 routine
+
+   REAL(wp), DIMENSION(jpjed:jpjef,jpj) ::   &  !:
+      sshfoe_b            !: east boundary ssh correction averaged over the barotropic loop
+                          !: (if Flather's algoritm applied at open boundary)
+
+   REAL(wp), DIMENSION(jpjed:jpjef,0:jptobc+1) ::   &  !:
+      sshedta, ubtedta    !: array used for interpolating monthly data on the east boundary
 
    REAL(wp), DIMENSION(jpjed:jpjef,jpk,jptobc) ::   &  !:
@@ -168 +179 @@
 
    REAL(wp), DIMENSION(jpjwd:jpjwf) ::   &  !:
-      bfow                !: now climatology of the west boundary barotropic stream function
+      bfow,             & !: now climatology of the west boundary barotropic stream function
+      sshfow,           & !: now climatology of the west boundary sea surface height
+      ubtfow,vbtfow       !: now climatology of the west boundary barotropic transport
 
    REAL(wp), DIMENSION(jpj,jpk) ::   &  !:
@@ -175 +188 @@
       ucliw               !: baroclinic componant of the zonal velocity after the radiation 
       !                   !  in the obcdyn.F90 routine
+
+   REAL(wp), DIMENSION(jpjwd:jpjwf,jpj) ::   &  !:
+      sshfow_b            !: west boundary ssh correction averaged over the barotropic loop
+                          !: (if Flather's algoritm applied at open boundary)
+
+   REAL(wp), DIMENSION(jpjwd:jpjwf,0:jptobc+1) ::   &  !:
+      sshwdta, ubtwdta    !: array used for interpolating monthly data on the west boundary
 
    REAL(wp), DIMENSION(jpjwd:jpjwf,jpk,jptobc) ::   &  !:
@@ -220 +240 @@
 
    REAL(wp), DIMENSION(jpind:jpinf) ::   &  !:
-      bfon                !: now climatology of the north boundary barotropic stream function
+      bfon,             & !: now climatology of the north boundary barotropic stream function
+      sshfon,           & !: now climatology of the north boundary sea surface height
+      ubtfon,vbtfon       !: now climatology of the north boundary barotropic transport
 
    REAL(wp), DIMENSION(jpi,jpk) ::   &    !:
@@ -227 +249 @@
       vclin               !: baroclinic componant of the meridian velocity after the radiation
       !                   !  in yhe obcdyn.F90 routine
+
+   REAL(wp), DIMENSION(jpind:jpinf,jpj) ::   &  !:
+      sshfon_b            !: north boundary ssh correction averaged over the barotropic loop
+                          !: (if Flather's algoritm applied at open boundary)
+
+   REAL(wp), DIMENSION(jpind:jpinf,0:jptobc+1) ::   &  !:
+      sshndta, vbtndta   !: array used for interpolating monthly data on the north boundary
 
    REAL(wp), DIMENSION(jpind:jpinf,jpk,jptobc) ::   &  !:
@@ -271 +300 @@
 
    REAL(wp), DIMENSION(jpisd:jpisf) ::    &   !:
-      bfos                !: now climatology of the south boundary barotropic stream function
+      bfos,             & !: now climatology of the south boundary barotropic stream function
+      sshfos,           & !: now climatology of the south boundary sea surface height
+      ubtfos,vbtfos       !: now climatology of the south boundary barotropic transport
 
    REAL(wp), DIMENSION(jpi,jpk) ::    &   !:
@@ -278 +309 @@
       vclis               !: baroclinic componant of the meridian velocity after the radiation 
       !                   !  in the obcdyn.F90 routine
+
+   REAL(wp), DIMENSION(jpisd:jpisf,jpj) ::   &  !:
+      sshfos_b            !: south boundary ssh correction averaged over the barotropic loop
+                          !: (if Flather's algoritm applied at open boundary)
+
+   REAL(wp), DIMENSION(jpisd:jpisf,0:jptobc+1) ::    &    !:
+      sshsdta, vbtsdta    !: array used for interpolating monthly data on the south boundary
 
    REAL(wp), DIMENSION(jpisd:jpisf,jpk,jptobc) ::    &    !:

trunk/NEMO/OPA_SRC/OBC/obccli.F90

@@ -4 +4 @@
    !! Ocean dynamics:   Baroclinic componant of velocities on each open boundary
    !!===================================================================================
-#if defined key_obc && ! defined key_dynspg_fsc
+#if defined key_obc && defined key_dynspg_rl
    !!-----------------------------------------------------------------------------------
    !!   'key_obc'               and 
-   !!   'key_dynspg_fsc'
+   !!   'key_dynspg_rl'
    !!-----------------------------------------------------------------------------------
    !!   obc_cli_dyn : Compute the baroclinic componant after the radiation phase

trunk/NEMO/OPA_SRC/OBC/obcdta.F90

@@ -20 +20 @@
    USE in_out_manager  ! I/O logical units
    USE lib_mpp         ! distributed memory computing
-   USE dynspg_rl       ! 
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE ioipsl
-
-
-#  if ! defined key_dynspg_fsc
+#  if defined key_dynspg_rl
    USE obccli
 #  endif
@@ -31 +29 @@
    PRIVATE
 
-!! * Accessibility
+   !! * Accessibility
    PUBLIC obc_dta        ! routines called by step.F90
+   PUBLIC obc_dta_bt     ! routines called by dynspg_ts.F90
 
    !! * Shared module variables
    INTEGER ::   &
-      nlecto =  0,    &  ! switch for the first read
-      ntobc1      ,   &  ! first record used
-      ntobc2             ! second record used
-   
+      nlecto,   &  ! switch for the first read
+      ntobc1,   &  ! first record used
+      ntobc2,   &  ! second record used
+      itobc        ! number of time steps in OBC files 
+
+   REAL(wp), DIMENSION(:), ALLOCATABLE :: ztcobc      ! time_counter variable of BCs
+
    !! * Substitutions
+#  include "domzgr_substitute.h90"                                            
 #  include "obc_vectopt_loop_substitute.h90"
    !!---------------------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
+   !! $Header$
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
    !!---------------------------------------------------------------------------------
 
@@ -82 +87 @@
       INTEGER ::  isrel              ! number of seconds since 1/1/1992
       INTEGER, SAVE ::  itobce, itobcw,  & ! number of time steps in OBC files
-                        itobcs, itobcn,  & !    "       "       "       "
-                        itobc
+                        itobcs, itobcn     !    "       "       "       "
       INTEGER ::  ikprint        ! frequency for printouts.
       INTEGER :: fid_e, fid_w, fid_n, fid_s, fid  ! file identifiers
@@ -91 +95 @@
             start,             & ! starting index read
             count                ! number of indices to be read
-      ! time_counter variable of BCs
-      REAL(wp),DIMENSION(:),ALLOCATABLE :: ztcobc
      
       CHARACTER(LEN=25) :: f_name,v_name
@@ -112 +114 @@
       IF( kt == nit000 )   THEN
       
+         nlecto =  0
+
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*)     'obc_dta : find boundary data'
@@ -219 +223 @@
        ! 1.1  Tangential velocities set to zero
        ! --------------------------------------
-         IF( lp_obc_east  ) vfoe = 0.0
-         IF( lp_obc_west  ) vfow = 0.0
-         IF( lp_obc_south ) ufos = 0.0
-         IF( lp_obc_north ) ufon = 0.0
+         IF( lp_obc_east  ) vfoe = 0.e0
+         IF( lp_obc_west  ) vfow = 0.e0
+         IF( lp_obc_south ) ufos = 0.e0
+         IF( lp_obc_north ) ufon = 0.e0
       
        ! 1.2  Data temperature, salinity, normal velocities set to zero
@@ -344 +348 @@
             itimo = imois   
          ELSE
-            IF(lwp) WRITE(numout,*) 'data other than constant or monthly not written yet'
-            STOP
+            IF(lwp) WRITE(numout,*) 'data other than constant or monthly',kt
+            iman  = itobc
+            itimo = FLOOR( kt*rdt / (ztcobc(2)-ztcobc(1)) )
+            isrel = kt*rdt
          ENDIF
       ENDIF
@@ -370 +376 @@
             ENDIF
          ELSE
-            !!!!!!!!!!!!!ATTENTION el: A verifier en fction de la convention choisie pour
-            !!!!!!!!!!!!! le codage de nyear, pour les runs interannuels !!!!!!!!!!!!!!
-            !!!  attention if ln_obc_clim is true, go back to jan 1st after december 31st
-            iyrel=nyear-1991
-            IF( ( iyrel < 1 ) .OR. ( iyrel > 13 ) )   THEN
-               IF( lwp ) WRITE(numout,*) 'Pb OBCDTA : iyrel'
-               STOP
-            ENDIF
-            ! Compute nb of seconds from 1/1/1992 00:00 :
-            isrel=(365*(iyrel-1)+nday_year)*86400
-            IF( lwp )   THEN
-               WRITE(numout,*)'Nbre de secondes ecoulees depuis le 1/1/1992:'
-               WRITE(numout,*) isrel
-            ENDIF
-      
-            !  need to calculate here ntobc1 and ntobc2, the two time steps to be read
-      
+            isrel=kt*rdt
+            ntobc1 = itimo         ! first file record used
+            ntobc2 = ntobc1 + 1    ! last  file record used
+            ntobc1 = MOD( ntobc1, iman )
+            IF( ntobc1 == 0 )   ntobc1 = iman
+            ntobc2 = MOD( ntobc2, iman )
+            IF( ntobc2 == 0 )   ntobc2 = iman
+            IF(lwp) WRITE(numout,*) ' read obc first record file used ntobc1 ', ntobc1
+            IF(lwp) WRITE(numout,*) ' read obc last  record file used ntobc2 ', ntobc2
          ENDIF
                                ! ======================= !
@@ -395 +393 @@
             ! ... initialise the sedta, tedta, uedta arrays
             CALL flioopfd ('obceast_TS.nc',fid_e)
-            CALL obc_dta_gv (fid_e,'y','vosaline',sedta(:,:,1),jpjef-jpjed+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_e,'y','vosaline',sedta(:,:,2),jpjef-jpjed+1,jpk,ntobc2)
-            CALL obc_dta_gv (fid_e,'y','votemper',tedta(:,:,1),jpjef-jpjed+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_e,'y','votemper',tedta(:,:,2),jpjef-jpjed+1,jpk,ntobc2)
-            CALL flioclo (fid_e)
-
-            CALL flioopfd ('obceast_U.nc',fid_e)
-            CALL obc_dta_gv (fid_e,'y','vozocrtx',uedta(:,:,1),jpjef-jpjed+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_e,'y','vozocrtx',uedta(:,:,2),jpjef-jpjed+1,jpk,ntobc2)
+            CALL obc_dta_gv (fid_e,'y','vosaline',jpjef-jpjed+1,ntobc1,pdta_3D=sedta(:,:,1))
+            CALL obc_dta_gv (fid_e,'y','vosaline',jpjef-jpjed+1,ntobc2,pdta_3D=sedta(:,:,2))
+            CALL obc_dta_gv (fid_e,'y','votemper',jpjef-jpjed+1,ntobc1,pdta_3D=tedta(:,:,1))
+            CALL obc_dta_gv (fid_e,'y','votemper',jpjef-jpjed+1,ntobc2,pdta_3D=tedta(:,:,2))
+            CALL flioclo (fid_e)                                                           
+                                                                                           
+            CALL flioopfd ('obceast_U.nc',fid_e)                                           
+            CALL obc_dta_gv (fid_e,'y','vozocrtx',jpjef-jpjed+1,ntobc1,pdta_3D=uedta(:,:,1))
+            CALL obc_dta_gv (fid_e,'y','vozocrtx',jpjef-jpjed+1,ntobc2,pdta_3D=uedta(:,:,2))
             CALL flioclo (fid_e)
             !  Usually printout is done only once at kt = nit000,
@@ -429 +427 @@
             ! ... initialise the swdta, twdta, uwdta arrays
             CALL flioopfd ('obcwest_TS.nc',fid_w)
-            CALL obc_dta_gv (fid_w,'y','vosaline',swdta(:,:,1),jpjwf-jpjwd+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_w,'y','vosaline',swdta(:,:,2),jpjwf-jpjwd+1,jpk,ntobc2)
-            CALL obc_dta_gv (fid_w,'y','votemper',twdta(:,:,1),jpjwf-jpjwd+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_w,'y','votemper',twdta(:,:,2),jpjwf-jpjwd+1,jpk,ntobc2)
-            CALL flioclo (fid_w)
-
-            CALL flioopfd ('obcwest_U.nc',fid_w)
-            CALL obc_dta_gv (fid_w,'y','vozocrtx',uwdta(:,:,1),jpjwf-jpjwd+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_w,'y','vozocrtx',uwdta(:,:,2),jpjwf-jpjwd+1,jpk,ntobc2)
+            CALL obc_dta_gv (fid_w,'y','vosaline',jpjwf-jpjwd+1,ntobc1,pdta_3D=swdta(:,:,1))
+            CALL obc_dta_gv (fid_w,'y','vosaline',jpjwf-jpjwd+1,ntobc2,pdta_3D=swdta(:,:,2))
+            CALL obc_dta_gv (fid_w,'y','votemper',jpjwf-jpjwd+1,ntobc1,pdta_3D=twdta(:,:,1))
+            CALL obc_dta_gv (fid_w,'y','votemper',jpjwf-jpjwd+1,ntobc2,pdta_3D=twdta(:,:,2))
+            CALL flioclo (fid_w)                                                           
+                                                                                           
+            CALL flioopfd ('obcwest_U.nc',fid_w)                                           
+            CALL obc_dta_gv (fid_w,'y','vozocrtx',jpjwf-jpjwd+1,ntobc1,pdta_3D=uwdta(:,:,1))
+            CALL obc_dta_gv (fid_w,'y','vozocrtx',jpjwf-jpjwd+1,ntobc2,pdta_3D=uwdta(:,:,2))
             CALL flioclo (fid_w)
 
@@ -460 +458 @@
          IF( lp_obc_north )   THEN      
             CALL flioopfd ('obcnorth_TS.nc',fid_n)
-            CALL obc_dta_gv (fid_n,'x','vosaline',sndta(:,:,1),jpinf-jpind+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_n,'x','vosaline',sndta(:,:,2),jpinf-jpind+1,jpk,ntobc2)
-            CALL obc_dta_gv (fid_n,'x','votemper',tndta(:,:,1),jpinf-jpind+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_n,'x','votemper',tndta(:,:,2),jpinf-jpind+1,jpk,ntobc2)
-            CALL flioclo (fid_n)
-
-            CALL flioopfd ('obcnorth_V.nc',fid_n)
-            CALL obc_dta_gv (fid_n,'x','vomecrty',vndta(:,:,1),jpinf-jpind+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_n,'x','vomecrty',vndta(:,:,2),jpinf-jpind+1,jpk,ntobc2)
+            CALL obc_dta_gv (fid_n,'x','vosaline',jpinf-jpind+1,ntobc1,pdta_3D=sndta(:,:,1))
+            CALL obc_dta_gv (fid_n,'x','vosaline',jpinf-jpind+1,ntobc2,pdta_3D=sndta(:,:,2))
+            CALL obc_dta_gv (fid_n,'x','votemper',jpinf-jpind+1,ntobc1,pdta_3D=tndta(:,:,1))
+            CALL obc_dta_gv (fid_n,'x','votemper',jpinf-jpind+1,ntobc2,pdta_3D=tndta(:,:,2))
+            CALL flioclo (fid_n)                                                           
+                                                                                           
+            CALL flioopfd ('obcnorth_V.nc',fid_n)                                          
+            CALL obc_dta_gv (fid_n,'x','vomecrty',jpinf-jpind+1,ntobc1,pdta_3D=vndta(:,:,1))
+            CALL obc_dta_gv (fid_n,'x','vomecrty',jpinf-jpind+1,ntobc2,pdta_3D=vndta(:,:,2))
             CALL flioclo (fid_n)
 
@@ -491 +489 @@
          IF( lp_obc_south )   THEN      
             CALL flioopfd ('obcsouth_TS.nc',fid_s)
-            CALL obc_dta_gv (fid_s,'x','vosaline',ssdta(:,:,1),jpisf-jpisd+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_s,'x','vosaline',ssdta(:,:,2),jpisf-jpisd+1,jpk,ntobc2)
-            CALL obc_dta_gv (fid_s,'x','votemper',tsdta(:,:,1),jpisf-jpisd+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_s,'x','votemper',tsdta(:,:,2),jpisf-jpisd+1,jpk,ntobc2)
-            CALL flioclo (fid_s)
-
-            CALL flioopfd ('obcsouth_V.nc',fid_s)
-            CALL obc_dta_gv (fid_s,'x','vomecrty',vsdta(:,:,1),jpisf-jpisd+1,jpk,ntobc1)
-            CALL obc_dta_gv (fid_s,'x','vomecrty',vsdta(:,:,2),jpisf-jpisd+1,jpk,ntobc2)
+            CALL obc_dta_gv (fid_s,'x','vosaline',jpisf-jpisd+1,ntobc1,pdta_3D=ssdta(:,:,1))
+            CALL obc_dta_gv (fid_s,'x','vosaline',jpisf-jpisd+1,ntobc2,pdta_3D=ssdta(:,:,2))
+            CALL obc_dta_gv (fid_s,'x','votemper',jpisf-jpisd+1,ntobc1,pdta_3D=tsdta(:,:,1))
+            CALL obc_dta_gv (fid_s,'x','votemper',jpisf-jpisd+1,ntobc2,pdta_3D=tsdta(:,:,2))
+            CALL flioclo (fid_s)                                                           
+                                                                                           
+            CALL flioopfd ('obcsouth_V.nc',fid_s)                                          
+            CALL obc_dta_gv (fid_s,'x','vomecrty',jpisf-jpisd+1,ntobc1,pdta_3D=vsdta(:,:,1))
+            CALL obc_dta_gv (fid_s,'x','vomecrty',jpisf-jpisd+1,ntobc2,pdta_3D=vsdta(:,:,2))
             CALL flioclo (fid_s)
 
@@ -519 +517 @@
             ENDIF
          ENDIF
-      
-      ENDIF        !      end of the test on the condition to read or not the files 
+
+      ELSE
+         
+         nlecto = 0        !      no reading of OBC barotropic data                         
+
+      ENDIF                !      end of the test on the condition to read or not the files 
       
       ! 3.  Call at every time step :
@@ -596 +598 @@
    END SUBROUTINE obc_dta
       
-# if defined key_dynspg_fsc
+# if defined key_dynspg_rl
    !!-----------------------------------------------------------------------------
-   !!   'key_dynspg_fsc'                    free surface with constant volume
-   !!-----------------------------------------------------------------------------
-   SUBROUTINE obc_dta_psi ( kt )       ! Empty routine
-      !! * Arguments
-      INTEGER,INTENT(in) :: kt 
-      WRITE(*,*) 'obc_dta_psi: You should not have seen this print! error?', kt
-   END SUBROUTINE obc_dta_psi
-#else
-   !!-----------------------------------------------------------------------------
-   !!   Default option                                                   Rigid-lid
+   !!   Rigid-lid
    !!-----------------------------------------------------------------------------
 
@@ -636 +629 @@
       !!        ! 97-08 (G. Madec, J.M. Molines)
       !!   8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
+      !!   9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
       !!----------------------------------------------------------------------------
       !! * Arguments
@@ -684 +678 @@
       IF( nbobc > 1 )   THEN
          DO jnic = 1,nbobc - 1
-            gcbic(jnic) = 0.
+            gcbic(jnic) = 0.e0
             ip=mnic(0,jnic)
             DO jip = 1,ip
@@ -742 +736 @@
       IF( lpsouthobc)   THEN
 
-         IF( kt == nit000.OR.kt <= kbsfstart )   THEN
+         IF( kt == nit000 .OR. kt <= kbsfstart )   THEN
             OPEN(inum,file='obcsouthbsf.dta')
             READ(inum,*)
@@ -759 +753 @@
 
       IF( lpnorthobc)   THEN
-         IF( kt == nit000.OR.kt <= kbsfstart )   THEN
+         IF( kt == nit000 .OR. kt <= kbsfstart )   THEN
             OPEN(inum,file='obcnorthbsf.dta')
             READ(inum,*)
@@ -776 +770 @@
 
    END SUBROUTINE obc_dta_psi
-
+#else
+   !!-----------------------------------------------------------------------------
+   !!   Default option                    
+   !!-----------------------------------------------------------------------------
+   SUBROUTINE obc_dta_psi ( kt )       ! Empty routine
+      !! * Arguments
+      INTEGER,INTENT(in) :: kt 
+      WRITE(*,*) 'obc_dta_psi: You should not have seen this print! error?', kt
+   END SUBROUTINE obc_dta_psi
 # endif
 
-   SUBROUTINE obc_dta_gv (ifid,cldim,clobc,pdta,kobcij,kobck,ktobc)
+
+#if defined key_dynspg_ts || defined key_dynspg_exp
+   SUBROUTINE obc_dta_bt( kt, kbt )
+      !!---------------------------------------------------------------------------
+      !!                      ***  SUBROUTINE obc_dta  ***
+      !!
+      !! ** Purpose :   time interpolation of barotropic data for time-splitting scheme
+      !!                Data at the boundary must be in m2/s 
+      !!
+      !! History :
+      !!   9.0  !  05-11 (V. garnier) Original code
+      !!---------------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk, ii, ij   ! dummy loop indices
+      INTEGER ::   fid_e, fid_w, fid_n, fid_s, fid  ! file identifiers
+      INTEGER ::   itimo, iman, imois, i15
+      INTEGER ::   ntobcm, ntobcp, itimom, itimop
+      REAL(wp) ::  zxy
+      INTEGER ::   isrel, ikt           ! number of seconds since 1/1/1992
+      INTEGER ::   ikprint              ! frequency for printouts.
+
+      !!---------------------------------------------------------------------------
+
+      ! 1.   First call: check time frames available in files.
+      ! -------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+
+         ! 1.1  Barotropic tangential velocities set to zero
+         ! -------------------------------------------------
+         IF( lp_obc_east  ) vbtfoe(:) = 0.e0
+         IF( lp_obc_west  ) vbtfow(:) = 0.e0
+         IF( lp_obc_south ) ubtfos(:) = 0.e0
+         IF( lp_obc_north ) ubtfon(:) = 0.e0
+
+         ! 1.2  Sea surface height and normal barotropic velocities set to zero
+         !                               or initial conditions if nobc_dta == 0
+         ! --------------------------------------------------------------------
+
+          IF( lp_obc_east ) THEN
+             ! initialisation to zero
+             sshedta(:,:) = 0.e0
+             ubtedta(:,:) = 0.e0
+             !                                        ! ================== !
+             IF( nobc_dta == 0 )   THEN               ! initial state used !
+                !                                     ! ================== !
+                !  Fills sedta, tedta, uedta (global arrays)
+                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
+                DO ji = nie0, nie1
+                   DO jj = nje0p1, nje1m1
+                      ij = jj -1 + njmpp
+                      sshedta(ij,1) = sshn(ji+1,jj) * tmask(ji+1,jj,1)
+                   END DO
+                END DO
+             ENDIF
+          ENDIF
+
+          IF( lp_obc_west) THEN
+             ! initialisation to zero
+             sshwdta(:,:) = 0.e0
+             ubtwdta(:,:) = 0.e0
+             !                                        ! ================== !
+             IF( nobc_dta == 0 )   THEN               ! initial state used !
+                !                                     ! ================== !
+                !  Fills swdta, twdta, uwdta (global arrays)
+                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
+                DO ji = niw0, niw1
+                   DO jj = njw0p1, njw1m1
+                      ij = jj -1 + njmpp
+                      sshwdta(ij,1) = sshn(ji,jj) * tmask(ji,jj,1)
+                   END DO
+                END DO
+             ENDIF
+          ENDIF
+
+          IF( lp_obc_north) THEN
+             ! initialisation to zero
+             sshndta(:,:) = 0.e0
+             vbtndta(:,:) = 0.e0
+             !                                        ! ================== !
+             IF( nobc_dta == 0 )   THEN               ! initial state used !
+                !                                     ! ================== !
+                !  Fills sndta, tndta, vndta (global arrays)
+                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
+                DO jj = njn0, njn1
+                   DO ji = nin0p1, nin1m1
+                      DO jk = 1, jpkm1
+                         ii = ji -1 + nimpp
+                         vbtndta(ii,1) = vbtndta(ii,1) + vndta(ii,jk,1)*fse3v(ji,jj,jk)
+                      END DO
+                      sshndta(ii,1) = sshn(ii,jj+1) * tmask(ji,jj+1,1)
+                   END DO
+                END DO
+             ENDIF
+          ENDIF
+
+          IF( lp_obc_south) THEN
+             ! initialisation to zero
+             ssdta(:,:,:) = 0.e0
+             tsdta(:,:,:) = 0.e0
+             vsdta(:,:,:) = 0.e0
+             sshsdta(:,:) = 0.e0
+             vbtsdta(:,:) = 0.e0
+             !                                        ! ================== !
+             IF( nobc_dta == 0 )   THEN               ! initial state used !
+                !                                     ! ================== !
+                !  Fills ssdta, tsdta, vsdta (global arrays)
+                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
+                DO jj = njs0, njs1
+                   DO ji = nis0p1, nis1m1
+                      DO jk = 1, jpkm1
+                         ii = ji -1 + nimpp
+                         vbtsdta(ii,1) = vbtsdta(ii,1) + vsdta(ii,jk,1)*fse3v(ji,jj,jk)
+                      END DO
+                      sshsdta(ii,1) = sshn(ji,jj) * tmask(ii,jj,1)
+                   END DO
+                END DO
+             ENDIF
+          ENDIF
+
+       ENDIF        !       END IF kt == nit000
+
+      !!------------------------------------------------------------------------------------
+      ! 2.      Initialize the time we are at. Does this every time the routine is called,
+      !         excepted when nobc_dta = 0
+      !
+      IF( nobc_dta == 0) THEN
+         itimo = 1
+         zxy   = 0
+      ELSE
+         IF(itobc == 1) THEN
+            itimo = 1
+         ELSE IF (itobc == 12) THEN      !   BC are monthly
+            ! we assume we have climatology in that case
+            iman  = 12
+            i15   = nday / 16
+            imois = nmonth + i15 - 1
+            IF( imois == 0 )   imois = iman
+            itimo = imois
+         ELSE
+            IF(lwp) WRITE(numout,*) 'data other than constant or monthly',kt
+            iman  = itobc
+            itimo = FLOOR( kt*rdt / ztcobc(1))
+            isrel=kt*rdt
+         ENDIF
+      ENDIF
+
+      ! 2. Read two records in the file if necessary
+      ! ---------------------------------------------
+
+      IF( nobc_dta == 1 .AND. nlecto == 1 ) THEN
+
+         IF( lp_obc_east ) THEN
+            ! ... Read datafile and set sea surface height and barotropic velocity
+            ! ... initialise the sshedta, ubtedta arrays
+            sshedta(:,0) = sshedta(:,1)
+            ubtedta(:,0) = ubtedta(:,1)
+            CALL flioopfd ('obceast_TS.nc',fid_e)
+            CALL obc_dta_gv (fid_e,'y','vossurfh',jpjef-jpjed+1,ntobc1,pdta_2D=sshedta(:,1))
+            CALL obc_dta_gv (fid_e,'y','vossurfh',jpjef-jpjed+1,ntobc2,pdta_2D=sshedta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_e,'y','vossurfh',jpjef-jpjed+1,ntobc2+1,pdta_2D=sshedta(:,3))
+            ENDIF
+            CALL flioclo (fid_e)
+
+            CALL flioopfd ('obceast_U.nc',fid_e)
+            CALL obc_dta_gv (fid_e,'y','vozoubt',jpjef-jpjed+1,ntobc1,pdta_2D=ubtedta(:,1))
+            CALL obc_dta_gv (fid_e,'y','vozoubt',jpjef-jpjed+1,ntobc2,pdta_2D=ubtedta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_e,'y','vozoubt',jpjef-jpjed+1,ntobc2+1,pdta_2D=ubtedta(:,3))
+            ENDIF
+            CALL flioclo (fid_e)
+
+            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
+            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) ' Read East OBC barotropic data records ', ntobc1, ntobc2
+               ikprint = (jpjef-jpjed+1)/20 +1
+               WRITE(numout,*)
+               WRITE(numout,*) ' Sea surface height record 1'
+               CALL prihre( sshedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, ikprint, 1, 1, -3, 1., numout )
+               WRITE(numout,*)
+               WRITE(numout,*) ' Normal transport (m2/s) record 1',ikprint
+               CALL prihre( ubtedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, ikprint, 1, 1, -3, 1., numout )
+            ENDIF
+         ENDIF
+
+         IF( lp_obc_west ) THEN
+            ! ... Read datafile and set temperature, salinity and normal velocity
+            ! ... initialise the swdta, twdta, uwdta arrays
+            sshwdta(:,0) = sshwdta(:,1)
+            ubtwdta(:,0) = ubtwdta(:,1)
+            CALL flioopfd ('obcwest_TS.nc',fid_w)
+            CALL obc_dta_gv (fid_w,'y','vossurfh',jpjwf-jpjwd+1,ntobc1,pdta_2D=sshwdta(:,1))
+            CALL obc_dta_gv (fid_w,'y','vossurfh',jpjwf-jpjwd+1,ntobc2,pdta_2D=sshwdta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_w,'y','vossurfh',jpjwf-jpjwd+1,ntobc2+1,pdta_2D=sshwdta(:,3))
+            ENDIF
+            CALL flioclo (fid_w)
+
+            CALL flioopfd ('obcwest_U.nc',fid_w)
+            CALL obc_dta_gv (fid_w,'y','vozoubt',jpjwf-jpjwd+1,ntobc1,pdta_2D=ubtwdta(:,1))
+            CALL obc_dta_gv (fid_w,'y','vozoubt',jpjwf-jpjwd+1,ntobc2,pdta_2D=ubtwdta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_w,'y','vozoubt',jpjwf-jpjwd+1,ntobc2+1,pdta_2D=ubtwdta(:,3))
+            ENDIF
+            CALL flioclo (fid_w)
+
+            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
+            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) ' Read West OBC barotropic data records ', ntobc1, ntobc2
+               ikprint = (jpjwf-jpjwd+1)/20 +1
+               WRITE(numout,*)
+               WRITE(numout,*) ' Sea surface height record 1  - printout surface level'
+               CALL prihre( sshwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, ikprint, 1, 1, -3, 1., numout )
+               WRITE(numout,*)
+               WRITE(numout,*) ' Normal transport (m2/s) record 1'
+               CALL prihre( ubtwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, ikprint, 1, 1, -3, 1., numout )
+            ENDIF
+         ENDIF
+
+         IF( lp_obc_north) THEN
+            ! ... Read datafile and set sea surface height and barotropic velocity
+            ! ... initialise the sshndta, ubtndta arrays
+            sshndta(:,0) = sshndta(:,1)
+            vbtndta(:,0) = vbtndta(:,1)
+            CALL flioopfd ('obcnorth_TS.nc',fid_n)
+            CALL obc_dta_gv (fid_n,'x','vossurfh',jpinf-jpind+1,ntobc1,pdta_2D=sshndta(:,1))
+            CALL obc_dta_gv (fid_n,'x','vossurfh',jpinf-jpind+1,ntobc2,pdta_2D=sshndta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_n,'x','vossurfh',jpinf-jpind+1,ntobc2+1,pdta_2D=sshndta(:,3))
+            ENDIF
+            CALL flioclo (fid_n)
+
+            CALL flioopfd ('obcnorth_V.nc',fid_n)
+            CALL obc_dta_gv (fid_n,'x','vomevbt',jpinf-jpind+1,ntobc1,pdta_2D=vbtndta(:,1))
+            CALL obc_dta_gv (fid_n,'x','vomevbt',jpinf-jpind+1,ntobc2,pdta_2D=vbtndta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_n,'x','vomevbt',jpinf-jpind+1,ntobc2+1,pdta_2D=vbtndta(:,3))
+            ENDIF
+            CALL flioclo (fid_n)
+
+            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
+            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) ' Read North OBC barotropic data records ', ntobc1, ntobc2
+               ikprint = (jpinf-jpind+1)/20 +1
+               WRITE(numout,*)
+               WRITE(numout,*) ' Sea surface height record 1  - printout surface level'
+               CALL prihre( sshndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, ikprint, 1, 1, -3, 1., numout )
+               WRITE(numout,*)
+               WRITE(numout,*) ' Normal transport (m2/s) record 1'
+               CALL prihre( vbtndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, ikprint, 1, 1, -3, 1., numout )
+            ENDIF
+         ENDIF
+
+         IF( lp_obc_south) THEN
+            ! ... Read datafile and set sea surface height and barotropic velocity
+            ! ... initialise the sshsdta, ubtsdta arrays
+            sshsdta(:,0) = sshsdta(:,1)
+            vbtsdta(:,0) = vbtsdta(:,1)
+            CALL flioopfd ('obcsouth_TS.nc',fid_s)
+            CALL obc_dta_gv (fid_s,'x','vossurfh',jpisf-jpisd+1,ntobc1,pdta_2D=sshsdta(:,1))
+            CALL obc_dta_gv (fid_s,'x','vossurfh',jpisf-jpisd+1,ntobc2,pdta_2D=sshsdta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_s,'x','vossurfh',jpisf-jpisd+1,ntobc2+1,pdta_2D=sshsdta(:,3))
+            ENDIF
+            CALL flioclo (fid_s)
+
+            CALL flioopfd ('obcsouth_V.nc',fid_s)
+            CALL obc_dta_gv (fid_s,'x','vomevbt',jpisf-jpisd+1,ntobc1,pdta_2D=vbtsdta(:,1))
+            CALL obc_dta_gv (fid_s,'x','vomevbt',jpisf-jpisd+1,ntobc2,pdta_2D=vbtsdta(:,2))
+            IF( lk_dynspg_ts ) THEN
+               CALL obc_dta_gv (fid_s,'x','vomevbt',jpisf-jpisd+1,ntobc2+1,pdta_2D=vbtsdta(:,3))
+            ENDIF
+            CALL flioclo (fid_s)
+                
+            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
+            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) ' Read South OBC barotropic data records ', ntobc1, ntobc2
+               ikprint = (jpisf-jpisd+1)/20 +1
+               WRITE(numout,*)
+               WRITE(numout,*) ' Sea surface height record 1  - printout surface level'
+               CALL prihre( sshsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, ikprint, 1, 1, -3, 1., numout )
+               WRITE(numout,*)
+               WRITE(numout,*) ' Normal transport (m2/s) record 1'
+               CALL prihre( vbtsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, ikprint, 1, 1, -3, 1., numout )
+            ENDIF
+         ENDIF
+
+       ENDIF        !      end of the test on the condition to read or not the files
+
+      ! 3.  Call at every time step : Linear interpolation of BCs to current time step
+      ! ----------------------------------------------------------------------
+
+       IF( lk_dynspg_ts ) THEN
+          isrel = (kt-1)*rdt + kbt*rdtbt
+
+          IF( nobc_dta == 1 ) THEN
+             isrel = (kt-1)*rdt + kbt*rdtbt
+             itimo  = FLOOR(  kt*rdt    / (ztcobc(2)-ztcobc(1)) )
+             itimom = FLOOR( (kt-1)*rdt / (ztcobc(2)-ztcobc(1)) )
+             itimop = FLOOR( (kt+1)*rdt / (ztcobc(2)-ztcobc(1)) )
+             IF( itimom == itimo .AND. itimop == itimo ) THEN
+                ntobcm = ntobc1
+                ntobcp = ntobc2
+
+             ELSEIF ( itimom <= itimo .AND. itimop == itimo ) THEN
+                IF(  FLOOR( isrel / (ztcobc(2)-ztcobc(1)) ) < itimo ) THEN
+                   ntobcm = ntobc1-1
+                   ntobcp = ntobc2-1
+                ELSE
+                   ntobcm = ntobc1
+                   ntobcp = ntobc2
+                ENDIF
+
+             ELSEIF ( itimom == itimo .AND. itimop >= itimo ) THEN
+                IF(  FLOOR( isrel / (ztcobc(2)-ztcobc(1)) ) < itimop ) THEN
+                   ntobcm = ntobc1
+                   ntobcp = ntobc2
+                ELSE
+                   ntobcm = ntobc1+1
+                   ntobcp = ntobc2+1
+                ENDIF
+
+             ELSEIF ( itimom == itimo-1 .AND. itimop == itimo+1 ) THEN
+                IF(  FLOOR( isrel / (ztcobc(2)-ztcobc(1)) ) < itimo ) THEN
+                   ntobcm = ntobc1-1
+                   ntobcp = ntobc2-1
+                ELSEIF (  FLOOR( isrel / (ztcobc(2)-ztcobc(1)) ) < itimop ) THEN
+                   ntobcm = ntobc1
+                   ntobcp = ntobc2
+                ELSEIF (  FLOOR( isrel / (ztcobc(2)-ztcobc(1)) ) == itimop ) THEN
+                   ntobcm = ntobc1+1
+                   ntobcp = ntobc2+2
+                ELSE
+                   IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 1?'
+                ENDIF
+             ELSE
+                IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 2?'
+             ENDIF
+
+          ENDIF
+
+       ELSE IF( lk_dynspg_exp ) THEN
+          isrel=kt*rdt
+          ntobcm = ntobc1
+          ntobcp = ntobc2
+       ENDIF
+
+       IF( itobc == 1 .OR. nobc_dta == 0 ) THEN
+          zxy = 0.e0
+       ELSE IF( itobc == 12 ) THEN
+          zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
+       ELSE
+          zxy = (ztcobc(ntobcm)-FLOAT(isrel)) / (ztcobc(ntobcm)-ztcobc(ntobcp))
+       ENDIF
+
+       IF( lp_obc_east ) THEN           !  fills sshfoe, ubtfoe (local to each processor)
+          DO jj = nje0p1, nje1m1
+             ij = jj -1 + njmpp
+             sshfoe(jj) = ( zxy * sshedta(ij,2) + (1.-zxy) * sshedta(ij,1) ) * temsk(jj,1)
+             ubtfoe(jj) = ( zxy * ubtedta(ij,2) + (1.-zxy) * ubtedta(ij,1) ) * uemsk(jj,1)
+          END DO
+       ENDIF
+
+       IF( lp_obc_west) THEN            !  fills sshfow, ubtfow (local to each processor)
+          DO jj = njw0p1, njw1m1
+             ij = jj -1 + njmpp
+             sshfow(jj) = ( zxy * sshwdta(ij,2) + (1.-zxy) * sshwdta(ij,1) ) * twmsk(jj,1)
+             ubtfow(jj) = ( zxy * ubtwdta(ij,2) + (1.-zxy) * ubtwdta(ij,1) ) * uwmsk(jj,1)
+          END DO
+       ENDIF
+
+       IF( lp_obc_north) THEN           !  fills sshfon, vbtfon (local to each processor)
+          DO ji = nin0p1, nin1m1
+             ii = ji -1 + nimpp
+             sshfon(ji) = ( zxy * sshndta(ii,2) + (1.-zxy) * sshndta(ii,1) ) * tnmsk(ji,1)
+             vbtfon(ji) = ( zxy * vbtndta(ii,2) + (1.-zxy) * vbtndta(ii,1) ) * vnmsk(ji,1)
+          END DO
+       ENDIF
+
+       IF( lp_obc_south) THEN           !  fills sshfos, vbtfos (local to each processor)
+          DO ji = nis0p1, nis1m1
+             ii = ji -1 + nimpp
+             sshfos(ji) = ( zxy * sshsdta(ii,2) + (1.-zxy) * sshsdta(ii,1) ) * tsmsk(ji,1)
+             vbtfos(ji) = ( zxy * vbtsdta(ii,2) + (1.-zxy) * vbtsdta(ii,1) ) * vsmsk(ji,1)
+          END DO
+       ENDIF
+
+   END SUBROUTINE obc_dta_bt
+
+#else
+   !!-----------------------------------------------------------------------------
+   !!   Default option
+   !!-----------------------------------------------------------------------------
+   SUBROUTINE obc_dta_bt ( kt, kbt )       ! Empty routine
+      !! * Arguments
+      INTEGER,INTENT(in) :: kt
+      INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index
+      WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt
+   END SUBROUTINE obc_dta_bt
+#endif
+
+
+   SUBROUTINE obc_dta_gv (ifid,cldim,clobc,kobcij,ktobc,pdta_2D,pdta_3D)
       !!-----------------------------------------------------------------------------
       !!                       ***  SUBROUTINE obc_dta_gv  ***
       !!
-      !! ** Purpose :   Read a OBC forcing field from netcdf file 
+      !! ** Purpose :   Read an OBC forcing field from netcdf file 
       !!                Input file are supposed to be 3D e.g.
       !!                - for a South or North OB : longitude x depth x time
@@ -794 +1207 @@
       !! * Arguments
       INTEGER, INTENT(IN) ::   &
-         ifid,                 & ! netcdf file name identifier
+         ifid  ,               & ! netcdf file name identifier
          kobcij,               & ! Horizontal (i or j) dimension of the array
-         kobck,                & ! vertical dimension
          ktobc                   ! starting time index read
       CHARACTER(LEN=*), INTENT(IN)    ::   &
          cldim,                & ! dimension along which is the open boundary ('x' or 'y')
          clobc                   ! name of the netcdf variable read
-      REAL, DIMENSION(kobcij,kobck,1), INTENT(OUT) ::   &
-         pdta                    ! 3D array of OBC forcing field
+      REAL, DIMENSION(kobcij,jpk,1), INTENT(OUT), OPTIONAL ::   &
+         pdta_3D                 ! 3D array of OBC forcing field
+      REAL, DIMENSION(kobcij,1), INTENT(OUT), OPTIONAL ::   &
+         pdta_2D                 ! 3D array of OBC forcing field
       
       !! * Local declarations
@@ -814 +1228 @@
       IF( l_exv )   THEN
          ! checks the number of dimensions
-         IF( indim == 3 )   THEN
-            istart(1:3) = (/ 1, 1, ktobc /)
-            icount(1:3) = (/ kobcij, kobck, 1 /)
-            CALL fliogetv (ifid,TRIM(clobc),pdta,start=istart(1:3),count=icount(1:3))
+         IF( indim == 2 )   THEN
+            istart(1:2) = (/ 1     , ktobc /)
+            icount(1:2) = (/ kobcij, 1     /)
+            CALL fliogetv (ifid,TRIM(clobc),pdta_2D,start=istart(1:2),count=icount(1:2))
+         ELSE IF( indim == 3 )   THEN
+            istart(1:3) = (/ 1     , 1    , ktobc /)
+            icount(1:3) = (/ kobcij, jpk  , 1     /)
+            CALL fliogetv (ifid,TRIM(clobc),pdta_3D,start=istart(1:3),count=icount(1:3))
          ELSE IF( indim == 4 )   THEN
             istart(1:4) = (/ 1, 1, 1, ktobc /)
             IF( TRIM(cldim) == 'y' )   THEN
-               icount(1:4) = (/ 1, kobcij, kobck, 1 /)
+               icount(1:4) = (/ 1     , kobcij, jpk  , 1 /)
             ELSE
-               icount(1:4) = (/ kobcij, 1, kobck, 1 /)
+               icount(1:4) = (/ kobcij, 1     , jpk  , 1 /)
             ENDIF
             ALLOCATE (v_tmp_4(icount(1),icount(2),icount(3),icount(4)))
             CALL fliogetv (ifid,TRIM(clobc),v_tmp_4,start=istart(1:4),count=icount(1:4))
             IF( TRIM(cldim) == 'y' )   THEN
-               pdta(1:kobcij,1:kobck,1:1) = v_tmp_4(1,1:kobcij,1:kobck,1:1)
+               pdta_3D(1:kobcij,1:jpk,1:1) = v_tmp_4(1,1:kobcij,1:jpk,1:1)
             ELSE
-               pdta(1:kobcij,1:kobck,1:1) = v_tmp_4(1:kobcij,1,1:kobck,1:1)
+               pdta_3D(1:kobcij,1:jpk,1:1) = v_tmp_4(1:kobcij,1,1:jpk,1:1)
             ENDIF
             DEALLOCATE (v_tmp_4)

trunk/NEMO/OPA_SRC/OBC/obcdyn.F90

@@ -22 +22 @@
    USE lbclnk          ! ???
    USE lib_mpp         ! ???
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE obccli          ! ocean open boundary conditions: climatology
    USE in_out_manager  ! I/O manager
@@ -29 +30 @@
 
    !! * Accessibility
-   PUBLIC obc_dyn     ! routine called in dynspg_fsc (free surface case)
+   PUBLIC obc_dyn     ! routine called in dynspg_flt (free surface case)
                       ! routine called in dynnxt.F90 (rigid lid case)
 
@@ -57 +58 @@
       !! ** Purpose :
       !!      Compute  dynamics (u,v) at the open boundaries.
-      !!      if defined key_dynspg_fsc: 
-      !!                 this routine is called by dynspg_fsc and updates
+      !!      if defined key_dynspg_flt: 
+      !!                 this routine is called by dynspg_flt and updates
       !!                 ua, va which are the actual velocities (not trends)
       !!      else  (rigid lid case) , 
@@ -74 +75 @@
       !!        !  97-07 (G. Madec, J.-M. Molines) addition
       !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90
+      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -132 +134 @@
       !!         ! 00-06 (J.-M. Molines) 
       !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
+      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
       !!------------------------------------------------------------------------------
       !! * Arguments
@@ -144 +147 @@
       ! --------------------------------------------------------
 
-      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbceast ) THEN
+      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbceast .OR. lk_dynspg_exp ) THEN
 
          ! 1.1 U zonal velocity    
@@ -151 +154 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-# if defined key_dynspg_fsc
-                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uemsk(jj,jk)) + &
-                                 uemsk(jj,jk)*ufoe(jj,jk)
-# else
+# if defined key_dynspg_rl
                   ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uemsk(jj,jk)) +                     &
                                  uemsk(jj,jk)*( ufoe(jj,jk) - hur (ji,jj) / e2u (ji,jj) &
                                  * ( bsfn(ji,jj) - bsfn(ji,jj-1) ) )
+# else
+                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uemsk(jj,jk)) + &
+                                 uemsk(jj,jk)*ufoe(jj,jk)
 # endif 
                END DO
@@ -220 +223 @@
             END DO
          END DO
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
          ! ... ua must be a baroclinic velocity uclie() 
          CALL obc_cli( ua, uclie, nie0, nie1, 0, jpj )  
@@ -294 +297 @@
       !!         ! 00-06 (J.-M. Molines) 
       !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
+      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
       !!------------------------------------------------------------------------------
       !! * Arguments
@@ -306 +310 @@
       ! --------------------------------------------------------
 
-      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcwest ) THEN
+      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcwest .OR. lk_dynspg_exp ) THEN
 
          ! 1.1 U zonal velocity
@@ -313 +317 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-# if defined key_dynspg_fsc
-                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uwmsk(jj,jk)) + &
-                                 uwmsk(jj,jk)*ufow(jj,jk)
-# else 
+# if defined key_dynspg_rl
                   ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uwmsk(jj,jk)) +                     &
                                  uwmsk(jj,jk)*( ufow(jj,jk) - hur (ji,jj) / e2u (ji,jj) &
                                  * ( bsfn(ji,jj) - bsfn(ji,jj-1) ) )
+# else 
+                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uwmsk(jj,jk)) + &
+                                 uwmsk(jj,jk)*ufow(jj,jk)
 # endif
                END DO
@@ -381 +385 @@
             END DO
          END DO
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
          ! ... ua must be a baroclinic velocity ucliw() 
          CALL obc_cli( ua, ucliw, niw0, niw1, 0, jpj )  
@@ -454 +458 @@
       !!         ! 00-06 (J.-M. Molines) 
       !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
+      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
       !!------------------------------------------------------------------------------
       !! * Arguments
@@ -466 +471 @@
       ! ---------------------------------------------------------
  
-      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcnorth ) THEN
+      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcnorth  .OR. lk_dynspg_exp ) THEN
 
          ! 1.1 U zonal velocity
@@ -484 +489 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-# if defined key_dynspg_fsc
-                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vnmsk(ji,jk)) + &
-                                vfon(ji,jk)*vnmsk(ji,jk)
-# else
+# if defined key_dynspg_rl
                   va(ji,jj,jk)= va(ji,jj,jk) * (1.-vnmsk(ji,jk)) +                       &
                                 vnmsk(ji,jk) * ( vfon(ji,jk) + hvr (ji,jj) / e1v (ji,jj) &
                                 * ( bsfn(ji,jj) - bsfn(ji-1,jj) ) )
+# else
+                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vnmsk(ji,jk)) + &
+                                vfon(ji,jk)*vnmsk(ji,jk)
 # endif
                END DO
@@ -590 +595 @@
             END DO
          END DO
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
          ! ... va must be a baroclinic velocity vclin() 
          CALL obc_cli( va, vclin, njn0, njn1, 1, jpi )  
@@ -625 +630 @@
       !!         ! 00-06 (J.-M. Molines) 
       !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
+      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
       !!------------------------------------------------------------------------------
       !! * Arguments
@@ -640 +646 @@
       ! ---------------------------------------------------------
 
-      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcsouth ) THEN
+      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcsouth  .OR. lk_dynspg_exp ) THEN
 
          ! 1.1 U zonal velocity
@@ -658 +664 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-# if defined key_dynspg_fsc
-                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vsmsk(ji,jk)) + &
-                                vsmsk(ji,jk) * vfos(ji,jk)
-# else 
+# if defined key_dynspg_rl
                   va(ji,jj,jk)= va(ji,jj,jk) * (1.-vsmsk(ji,jk)) +                      &
                                 vsmsk(ji,jk) * (vfos(ji,jk) + hvr (ji,jj) / e1v (ji,jj) &
                                 * ( bsfn(ji,jj) - bsfn(ji-1,jj) ) )
+# else 
+                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vsmsk(ji,jk)) + &
+                                vsmsk(ji,jk) * vfos(ji,jk)
 # endif
                END DO
@@ -758 +764 @@
             END DO
          END DO
-# if ! defined key_dynspg_fsc 
+# if defined key_dynspg_rl 
          ! ... va must be a baroclinic velocity vclis() 
          CALL obc_cli( va, vclis, njs0, njs1, 1, jpi )  

trunk/NEMO/OPA_SRC/OBC/obcini.F90

@@ -55 +55 @@
       !!        !  97-11  (J.M. Molines)
       !!   8.5  !  02-11  (C. Talandier, A-M. Treguier) Free surface, F90
+      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
       !!----------------------------------------------------------------------
       !! * Modules used
@@ -69 +70 @@
          &             rdpeob, rdpwob, rdpnob, rdpsob,   &
          &             zbsic1, zbsic2, zbsic3,           &
-         &             nbic, volemp, nobc_dta, ln_obc_clim
+         &             nbic, volemp, nobc_dta,           &
+         &             ln_obc_clim, ln_vol_cst, ln_obc_fla
       !!----------------------------------------------------------------------
 
@@ -135 +137 @@
       IF(lwp) WRITE(numout,*) '         initial state used (=0)             '
       IF(lwp) WRITE(numout,*) '         climatology (true) or not:', ln_obc_clim
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) '                                 WARNING                                                  '
+      IF(lwp) WRITE(numout,*) '         Flather"s algorithm is applied with explicit free surface scheme                 '
+      IF(lwp) WRITE(numout,*) '         or with free surface time-splitting scheme                                       '
+      IF(lwp) WRITE(numout,*) '         Nor radiation neither relaxation is allowed with explicit free surface scheme:   '
+      IF(lwp) WRITE(numout,*) '         Radiation and/or relaxation is allowed with free surface time-splitting scheme '
+      IF(lwp) WRITE(numout,*) '         depending of the choice of rdpXin = rdpXob  = 0. for open boundaries             '
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) '         For the rigid-lid case or the filtered free surface case,                        '
+      IF(lwp) WRITE(numout,*) '         radiation, relaxation or presciption of data can be applied                      '
       IF( lwp.AND.lp_obc_east ) THEN
          WRITE(numout,*) '         East open boundary :'
@@ -317 +329 @@
          !... (jpjed,jpjefm1),jpieob
          DO jj = nje0, nje1m1
-# if defined key_dynspg_fsc
+# if defined key_dynspg_rl
+            DO ji = nie0, nie1
+# else
             DO ji = nie0p1, nie1p1
-# else
-            DO ji = nie0, nie1
 # endif
                bmask(ji,jj) = 0.e0
@@ -368 +380 @@
       IF( lp_obc_north ) THEN
          ! ... jpjnob,(jpind,jpisfm1)
-# if defined key_dynspg_fsc
+# if defined key_dynspg_rl
+         DO jj = njn0, njn1
+# else
          DO jj = njn0p1, njn1p1
-# else
-         DO jj = njn0, njn1
 # endif
             DO ji = nin0, nin1m1
@@ -418 +430 @@
       END IF
 
-# if defined key_dynspg_fsc
+# if defined key_dynspg_flt
 
       ! ... Initialize obcumask and obcvmask for the Force filtering 
-      !     boundary condition in dynspg_fsc
+      !     boundary condition in dynspg_flt
       obcumask(:,:) = umask(:,:,1)
       obcvmask(:,:) = vmask(:,:,1)
@@ -502 +514 @@
       END IF
 
-      ! 3.1 Total lateral surface for each open boundary
-      ! ------------------------------------------------
-
-      obcsurftot = 0.e0
-      
-      IF( lp_obc_west ) THEN    ! ... West open boundary vertical surface
-         DO ji = niw0, niw1
-            DO jj = 1, jpj 
-               obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1)
-            END DO
-         END DO
-       END IF
-
-      IF( lp_obc_east ) THEN    ! ... East open boundary vertical surface
-         DO ji = nie0, nie1
-            DO jj = 1, jpj 
-               obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1)
-            END DO
-         END DO
-      END IF
-
-      IF( lp_obc_north ) THEN   ! ... North open boundary vertical surface
-         DO jj = njn0, njn1
-            DO ji = 1, jpi
-               obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vnmsk(ji,1)
-            END DO
-         END DO
-      END IF
-
-      IF( lp_obc_south ) THEN   ! ... South open boundary vertical surface
-         DO jj = njs0, njs1
-            DO ji = 1, jpi
-               obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vsmsk(ji,1)
-            END DO
-         END DO
-      END IF
-      IF( lk_mpp )   CALL mpp_sum( obcsurftot )   ! sum over the global domain
-
+# endif
+
+# if ! defined key_dynspg_rl
+
+      IF ( ln_vol_cst ) THEN
+
+         ! 3.1 Total lateral surface for each open boundary
+         ! ------------------------------------------------
+
+         ! ... West open boundary surface
+         IF( lp_obc_west ) THEN
+            DO ji = niw0, niw1
+               DO jj = 1, jpj 
+                  obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1)
+               END DO
+            END DO
+         END IF
+
+         ! ... East open boundary surface
+         IF( lp_obc_east ) THEN
+            DO ji = nie0, nie1
+               DO jj = 1, jpj 
+                  obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1)
+               END DO
+            END DO
+         END IF
+
+         ! ... North open boundary vertical surface
+         IF( lp_obc_north ) THEN
+            DO jj = njn0, njn1
+               DO ji = 1, jpi
+                  obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vnmsk(ji,1)
+               END DO
+            END DO
+         END IF
+
+         ! ... South open boundary vertical surface
+         IF( lp_obc_south ) THEN
+            DO jj = njs0, njs1
+               DO ji = 1, jpi
+                  obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vsmsk(ji,1)
+               END DO
+            END DO
+         END IF
+         IF( lk_mpp )   CALL mpp_sum( obcsurftot )   ! sum over the global domain
+      ENDIF
 # endif
 
@@ -712 +732 @@
       END IF
 
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
       ! 7. Isolated coastline arrays initialization (rigid lid case only)
       ! -----------------------------------------------------------------

trunk/NEMO/OPA_SRC/OBC/obcrad.F90

@@ -134 +134 @@
          ! ... fields nit <== now (kt+1) 
          ! ... Total or baroclinic velocity at b, bm and bm2
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zucb   = uclie(jj,jk)
 # else
                   zucb   = un(ji,jj,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zucbm  = un(ji-1,jj,jk) + hur(ji-1,jj) / e2u(ji-1,jj) &
                            * ( bsfn(ji-1,jj) - bsfn(ji-1,jj-1) )
@@ -145 +145 @@
                   zucbm  = un(ji-1,jj,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zucbm2 = un(ji-2,jj,jk) + hur(ji-2,jj) / e2u(ji-2,jj) &
                            * ( bsfn(ji-2,jj) - bsfn(ji-2,jj-1) )
@@ -412 +412 @@
                   uwbnd(jj,jk,nibm2,nitm) = uwbnd(jj,jk,nibm2,nit)*uwmsk(jj,jk)
          ! ... total or baroclinic velocity at b, bm and bm2
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zucb   = ucliw(jj,jk) 
 # else
                   zucb   = un (ji,jj,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zucbm  = un (ji+1,jj,jk) + hur (ji+1,jj) / e2u (ji+1,jj) &
                            * ( bsfn(ji+1,jj) - bsfn(ji+1,jj-1) )
@@ -423 +423 @@
                   zucbm  = un (ji+1,jj,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zucbm2 = un (ji+2,jj,jk) + hur (ji+2,jj) / e2u (ji+2,jj) &
                            * ( bsfn(ji+2,jj) - bsfn(ji+2,jj-1) )
@@ -738 +738 @@
          ! ... fields nit <== now (kt+1)
          ! ... total or baroclinic velocity at b, bm and bm2
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zvcb   = vclin(ji,jk)
 # else
                   zvcb   = vn (ji,jj,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zvcbm  = vn (ji,jj-1,jk) - hvr (ji,jj-1) / e1v (ji,jj-1) &
                            * ( bsfn(ji,jj-1) - bsfn(ji-1,jj-1) )
@@ -749 +749 @@
                   zvcbm  = vn (ji,jj-1,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zvcbm2 = vn (ji,jj-2,jk) - hvr (ji,jj-2) / e1v (ji,jj-2) &
                            * ( bsfn(ji,jj-2) - bsfn(ji-1,jj-2) )
@@ -1026 +1026 @@
                   vsbnd(ji,jk,nibm2,nitm) = vsbnd(ji,jk,nibm2,nit)*vsmsk(ji,jk)
          ! ... total or baroclinic velocity at b, bm and bm2
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zvcb   = vclis(ji,jk)
 # else
                   zvcb   = vn (ji,jj,jk)
 # endif
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                   zvcbm  = vn (ji,jj+1,jk) - hvr (ji,jj+1) / e1v (ji,jj+1) & 
                            * ( bsfn(ji,jj+1) - bsfn(ji-1,jj+1) )
@@ -1037 +1037 @@
                   zvcbm  = vn (ji,jj+1,jk)
 # endif
-# if ! defined key_dynspg_fsc  
+# if defined key_dynspg_rl  
                   zvcbm2 = vn (ji,jj+2,jk) - hvr (ji,jj+2) / e1v (ji,jj+2) &
                            * ( bsfn(ji,jj+2) - bsfn(ji-1,jj+2) )

trunk/NEMO/OPA_SRC/OBC/obcrst.F90

@@ -130 +130 @@
                      PRINT *,'Narea =',narea,' write jrec =2 east'
                      WRITE(inum,REC=jrec)                                    &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                            ((  bebnd(jfoe,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -145 +145 @@
                         jrec = 2 + jj + njmpp -1 -jpjed
                         WRITE (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                               ((  bebnd(jfoe,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -175 +175 @@
                      PRINT *,'Narea =',narea,' write jrec =',jrec,' west'
                      WRITE (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                            ((  bwbnd(jfow,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -190 +190 @@
                         jrec = 3 + jpjef -jpjed + jj + njmpp -1 -jpjwd
                         WRITE (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                               ((  bwbnd(jfow,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -219 +219 @@
                      ifon = jpind -nimpp +1
                      WRITE (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                            ((  bnbnd(ifon,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -234 +234 @@
                         jrec = 4 + jpjef -jpjed + jpjwf -jpjwd +ji + nimpp -1  -jpind
                         WRITE (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                               ((  bnbnd(ifon,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -264 +264 @@
                      ifos = jpisd -nimpp + 1
                      WRITE (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                            ((  bsbnd(ifos,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -280 +280 @@
                               ji + nimpp -1 -jpisd
                         WRITE (inum,REC=jrec) &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                               ((  bsbnd(ifos,   jb,jt),          jb=1,3),jt=1,3), &
 # endif 
@@ -520 +520 @@
                jfoe = jpjed -njmpp + 1
                READ (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                     ((  bebnd(jfoe,   jb,jt),          jb=1,3),jt=1,3), &
 # endif 
@@ -535 +535 @@
                   jrec = 2 + jj + njmpp -1 -jpjed
                   READ (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                        ((  bebnd(jfoe,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -562 +562 @@
                jfow = jpjwd -njmpp + 1
                READ (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                     ((  bwbnd(jfow,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -577 +577 @@
                   jrec = 3 + jpjef -jpjed + jj + njmpp -1 -jpjwd
                   READ (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                        ((  bwbnd(jfow,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -604 +604 @@
                ifon = jpind -nimpp +1
                READ (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                     ((  bnbnd(ifon,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -619 +619 @@
                   jrec = 4 + jpjef -jpjed + jpjwf -jpjwd +ji + nimpp -1  -jpind
                   READ (inum,REC=jrec)                                   & 
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                        ((  bnbnd(ifon,   jb,jt),          jb=1,3),jt=1,3), &
 # endif  
@@ -646 +646 @@
                ifos = jpisd -nimpp + 1
                READ (inum,REC=jrec)                                   &
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                     ((  bsbnd(ifos,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -662 +662 @@
                         ji + nimpp -1 -jpisd
                   READ (inum,REC=jrec)                                   & 
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
                        ((  bsbnd(ifos,   jb,jt),          jb=1,3),jt=1,3), &
 # endif
@@ -680 +680 @@
       IF( lk_mpp ) THEN
          IF( lp_obc_east ) THEN
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
             CALL mppobc(bebnd,jpjed,jpjef,jpieob,3*3,2,jpj)
 # endif
@@ -689 +689 @@
          ENDIF
          IF( lp_obc_west ) THEN
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
             CALL mppobc(bwbnd,jpjwd,jpjwf,jpiwob,3*3,2,jpj)
 # endif
@@ -698 +698 @@
          ENDIF
          IF( lp_obc_north ) THEN 
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
             CALL mppobc(bnbnd,jpind,jpinf,jpjnob  ,3*3    ,1,jpi)
 # endif
@@ -707 +707 @@
          ENDIF
          IF( lp_obc_south ) THEN
-# if ! defined key_dynspg_fsc
+# if defined key_dynspg_rl
             CALL mppobc(bsbnd,jpisd,jpisf,jpjsob,    3*3,1,jpi)
 # endif

trunk/NEMO/OPA_SRC/OBC/obcvol.F90

@@ -4 +4 @@
    !! Ocean dynamic :  Volume constraint when OBC and Free surface are used
    !!=================================================================================
-#if   defined key_obc   &&   defined key_dynspg_fsc
+#if   defined key_obc   &&   ! defined key_dynspg_rl
    !!---------------------------------------------------------------------------------
    !!   'key_obc'               and                           open boundary conditions
-   !!   'key_dynspg_fsc'                                  constant volume free surface
+   !!   'key_dynspg_flt'                                  constant volume free surface
    !!---------------------------------------------------------------------------------
    !! * Modules used
@@ -22 +22 @@
 
    !! * Accessibility
-   PUBLIC obc_vol        ! routine called by dynspg_fsc.h90
+   PUBLIC obc_vol        ! routine called by dynspg_flt
 
    !! * Substitutions
@@ -40 +40 @@
       !!
       !! ** Purpose : 
-      !!      This routine is called in dynspg_fsc to control 
+      !!      This routine is called in dynspg_flt to control 
       !!      the volume of the system. A correction velocity is calculated
       !!      to correct the total transport through the OBC. 

trunk/NEMO/OPA_SRC/SOL/solver.F90

@@ -19 +19 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
 
    IMPLICIT NONE

trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90

@@ -22 +22 @@
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE prtctl          ! Print control
 

trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90

@@ -15 +15 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE trabbl          ! tracers: bottom boundary layer
    USE lib_mpp         ! distribued memory computing

trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90

@@ -15 +15 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE trabbl          ! tracers: bottom boundary layer
    USE lib_mpp

trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90

@@ -17 +17 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE trabbl          ! Advective term of BBL
    USE lib_mpp

trunk/NEMO/OPA_SRC/opa.F90

@@ -47 +47 @@
 
    USE step            ! OPA time-stepping                  (stp     routine)
-   USE dynspg          ! Control choice of surface pressure gradient schemes
+   USE dynspg_oce      ! Control choice of surface pressure gradient schemes
    USE prtctl          ! Print control                 (prt_ctl_init routine)
    USE ini1d           ! re-initialization of u-v mask for the 1D configuration

trunk/NEMO/OPA_SRC/restart.F90

@@ -21 +21 @@
    USE blk_oce         ! bulk variables
    USE flx_oce         ! sea-ice/ocean forcings variables
-   USE dynspg          ! choice/control of key cpp for surface pressure gradient
-   USE dynspg_ts       ! free surface time splitting scheme variables
+   USE dynspg_oce      ! free surface time splitting scheme variables
    USE cpl_oce,         ONLY : lk_cpl              !
 

trunk/NEMO/OPA_SRC/step.F90

@@ -35 +35 @@
    USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine)
    USE dynhpg_atsk     ! hydrostatic pressure grad.  (dyn_hpg_atsk routine)
+   USE dynspg_oce      ! surface pressure gradient        (dyn_spg routine)
    USE dynspg          ! surface pressure gradient        (dyn_spg routine)
    USE dynkeg          ! kinetic energy gradient          (dyn_keg routine)

trunk/NEMO/OPA_SRC/stpctl.F90

@@ -17 +17 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! distributed memory computing
-   USE dynspg          ! pressure gradient schemes 
+   USE dynspg_oce      ! pressure gradient schemes 
 
    IMPLICIT NONE