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

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

Changeset 1502


Ignore:
Timestamp:
2009-07-20T17:20:23+02:00 (15 years ago)
Author:
rblod
Message:

Update dynnxt and dynspg_ts for variable volume, see ticket #474

Location:
trunk/NEMO/OPA_SRC
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/OPA_SRC/BDY/bdydyn.F90

    r1146 r1502  
    77   !!             -   !  2007-07  (D. Storkey) Move Flather implementation to separate routine. 
    88   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
     9   !!            3.2  !  2008-04  (R. Benshila) consider velocity instead of transport  
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_bdy  
    … …  
    8687         END DO  
    8788         ! 
    88          CALL lbc_lnk( ua, 'U', 1. )   ! Boundary points should be updated 
    89          CALL lbc_lnk( va, 'V', 1. )   ! 
     89         CALL lbc_lnk( ua, 'U', -1. )   ! Boundary points should be updated 
     90         CALL lbc_lnk( va, 'V', -1. )   ! 
    9091         ! 
    9192      ENDIF ! ln_bdy_dyn_frs 
    … …  
    9697#if defined key_dynspg_exp || defined key_dynspg_ts 
    9798!! Option to use Flather with dynspg_flt not coded yet... 
    98    SUBROUTINE bdy_dyn_fla 
     99   SUBROUTINE bdy_dyn_fla( pssh ) 
    99100      !!---------------------------------------------------------------------- 
    100101      !!                 ***  SUBROUTINE bdy_dyn_fla  *** 
    … …  
    116117      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.      
    117118      !!---------------------------------------------------------------------- 
     119      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh 
     120 
    118121      INTEGER  ::   ib, igrd                         ! dummy loop indices 
    119122      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
    120123      REAL(wp) ::   zcorr                            ! Flather correction 
     124      REAL(wp) ::   zforc                            ! temporary scalar 
    121125      !!---------------------------------------------------------------------- 
    122126 
    … …  
    127131      IF(ln_bdy_dyn_fla .OR. ln_bdy_tides) THEN ! If these are both false, then this routine does nothing.  
    128132 
    129       ! Fill temporary array with ssh data (here spgu): 
    130       igrd = 1 
    131       spgu(:,:) = 0.0 
    132       DO ib = 1, nblenrim(igrd) 
    133          ii = nbi(ib,igrd) 
    134          ij = nbj(ib,igrd) 
    135          IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib) 
    136          IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(ib) 
    137       END DO 
     133         ! Fill temporary array with ssh data (here spgu): 
     134         igrd = 1 
     135         spgu(:,:) = 0.0 
     136         DO ib = 1, nblenrim(igrd) 
     137            ii = nbi(ib,igrd) 
     138            ij = nbj(ib,igrd) 
     139            IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib) 
     140            IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(ib) 
     141         END DO 
     142         ! 
     143         igrd = 2      ! Flather bc on u-velocity;  
     144         !             ! remember that flagu=-1 if normal velocity direction is outward 
     145         !             ! I think we should rather use after ssh ? 
     146         DO ib = 1, nblenrim(igrd) 
     147            ii  = nbi(ib,igrd) 
     148            ij  = nbj(ib,igrd)  
     149            iim1 = ii + MAX( 0, INT( flagu(ib) ) )   ! T pts i-indice inside the boundary 
     150            iip1 = ii - MIN( 0, INT( flagu(ib) ) )   ! T pts i-indice outside the boundary  
     151            ! 
     152            zcorr = - flagu(ib) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
     153            zforc = ubtbdy(ib) + utide(ib) 
     154            ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
     155         END DO 
     156         ! 
     157         igrd = 3      ! Flather bc on v-velocity 
     158         !             ! remember that flagv=-1 if normal velocity direction is outward 
     159         DO ib = 1, nblenrim(igrd) 
     160            ii  = nbi(ib,igrd) 
     161            ij  = nbj(ib,igrd)  
     162            ijm1 = ij + MAX( 0, INT( flagv(ib) ) )   ! T pts j-indice inside the boundary 
     163            ijp1 = ij - MIN( 0, INT( flagv(ib) ) )   ! T pts j-indice outside the boundary  
     164            ! 
     165            zcorr = - flagv(ib) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
     166            zforc = vbtbdy(ib) + vtide(ib) 
     167            va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
     168         END DO 
     169         CALL lbc_lnk( ua_e, 'U', 1. )   ! Boundary points should be updated 
     170         CALL lbc_lnk( va_e, 'V', 1. )   ! 
     171         ! 
     172      ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides 
    138173      ! 
    139       igrd = 2      ! Flather bc on u-velocity;  
    140       !             ! remember that flagu=-1 if normal velocity direction is outward 
    141       !             ! I think we should rather use after ssh ? 
    142       DO ib = 1, nblenrim(igrd) 
    143          ii  = nbi(ib,igrd) 
    144          ij  = nbj(ib,igrd)  
    145          iim1 = ii + MAX( 0, INT( flagu(ib) ) )   ! T pts i-indice inside the boundary 
    146          iip1 = ii - MIN( 0, INT( flagu(ib) ) )   ! T pts i-indice outside the boundary  
    147          ! 
    148          zcorr = - flagu(ib) * SQRT( grav / (hu_e(ii, ij) + 1.e-20) ) * ( sshn_e(iim1, ij) - spgu(iip1,ij) ) 
    149          ua_e(ii, ij) = ( ubtbdy(ib) + utide(ib) ) * hu_e(ii,ij) 
    150          ua_e(ii,ij) = ua_e(ii,ij) + zcorr * umask(ii,ij,1) * hu_e(ii,ij) 
    151       END DO 
    152       ! 
    153       igrd = 3      ! Flather bc on v-velocity 
    154       !             ! remember that flagv=-1 if normal velocity direction is outward 
    155       DO ib = 1, nblenrim(igrd) 
    156          ii  = nbi(ib,igrd) 
    157          ij  = nbj(ib,igrd)  
    158          ijm1 = ij + MAX( 0, INT( flagv(ib) ) )   ! T pts j-indice inside the boundary 
    159          ijp1 = ij - MIN( 0, INT( flagv(ib) ) )   ! T pts j-indice outside the boundary  
    160          ! 
    161          zcorr = - flagv(ib) * SQRT( grav / (hv_e(ii, ij) + 1.e-20) ) * ( sshn_e(ii, ijm1) - spgu(ii,ijp1) ) 
    162          va_e(ii, ij) = ( vbtbdy(ib) + vtide(ib) ) * hv_e(ii,ij) 
    163          va_e(ii,ij) = va_e(ii,ij) + zcorr * vmask(ii,ij,1) * hv_e(ii,ij) 
    164       END DO 
    165       ! 
    166       CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated 
    167       CALL lbc_lnk( va_e, 'V', 1. ) ! 
    168       ! 
    169       ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides 
    170  
    171174   END SUBROUTINE bdy_dyn_fla 
    172175#endif 

  • trunk/NEMO/OPA_SRC/DYN/dynnxt.F90

    r1438 r1502  
    11MODULE dynnxt 
    2    !!====================================================================== 
     2   !!========================================================================= 
    33   !!                       ***  MODULE  dynnxt  *** 
    44   !! Ocean dynamics: time stepping 
    5    !!====================================================================== 
    6    !!====================================================================== 
     5   !!========================================================================= 
    76   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code 
    87   !!                 !  1990-10  (C. Levy, G. Madec) 
    … …  
    1514   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
    1615   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.  
    17    !!            3.2  !  2009-04  (G. Madec, R.Benshila))  re-introduce the vvl option 
    18    !!---------------------------------------------------------------------- 
     16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option 
     17   !!------------------------------------------------------------------------- 
    1918   
    20    !!---------------------------------------------------------------------- 
    21    !!   dyn_nxt      : update the horizontal velocity from the momentum trend 
    22    !!---------------------------------------------------------------------- 
     19   !!------------------------------------------------------------------------- 
     20   !!   dyn_nxt      : obtain the next (after) horizontal velocity 
     21   !!------------------------------------------------------------------------- 
    2322   USE oce             ! ocean dynamics and tracers 
    2423   USE dom_oce         ! ocean space and time domain 
    25    USE in_out_manager  ! I/O manager 
     24   USE dynspg_oce      ! type of surface pressure gradient 
     25   USE dynadv          ! dynamics: vector invariant versus flux form 
     26   USE domvvl          ! variable volume 
    2627   USE obc_oce         ! ocean open boundary conditions 
    2728   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine) 
    … …  
    3132   USE bdydta          ! unstructured open boundary conditions 
    3233   USE bdydyn          ! unstructured open boundary conditions 
    33    USE dynspg_oce      ! type of surface pressure gradient 
     34   USE agrif_opa_update 
     35   USE agrif_opa_interp 
     36   USE in_out_manager  ! I/O manager 
    3437   USE lbclnk          ! lateral boundary condition (or mpp link) 
    3538   USE prtctl          ! Print control 
    36    USE agrif_opa_update 
    37    USE agrif_opa_interp 
    38    USE domvvl          ! variable volume 
    3939 
    4040   IMPLICIT NONE 
    … …  
    5757      !!                  ***  ROUTINE dyn_nxt  *** 
    5858      !!                    
    59       !! ** Purpose :   Compute the after horizontal velocity from the  
    60       !!      momentum trend. 
    61       !! 
    62       !! ** Method  :   Apply lateral boundary conditions on the trends (ua,va)  
    63       !!      through calls to routine lbc_lnk. 
    64       !!      After velocity is compute using a leap-frog scheme environment: 
    65       !!         (ua,va) = (ub,vb) + 2 rdt (ua,va) 
    66       !!      Note that if lk_dynspg_flt=T, the time stepping has already been 
    67       !!      performed in dynspg module 
    68       !!      Time filter applied on now horizontal velocity to avoid the 
    69       !!      divergence of two consecutive time-steps and swap of dynamics 
    70       !!      arrays to start the next time step: 
    71       !!         (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 
    72       !!         (un,vn) = (ua,va)  
    73       !! 
    74       !! ** Action : - Update ub,vb arrays, the before horizontal velocity 
    75       !!             - Update un,vn arrays, the now horizontal velocity 
    76       !! 
     59      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary  
     60      !!             condition on the after velocity, achieved the time stepping  
     61      !!             by applying the Asselin filter on now fields and swapping  
     62      !!             the fields. 
     63      !! 
     64      !! ** Method  : * After velocity is compute using a leap-frog scheme: 
     65      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va) 
     66      !!             Note that with flux form advection and variable volume layer 
     67      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted 
     68      !!             velocity. 
     69      !!             Note also that in filtered free surface (lk_dynspg_flt=T), 
     70      !!             the time stepping has already been done in dynspg module 
     71      !! 
     72      !!              * Apply lateral boundary conditions on after velocity  
     73      !!             at the local domain boundaries through lbc_lnk call, 
     74      !!             at the radiative open boundaries (lk_obc=T), 
     75      !!             at the relaxed   open boundaries (lk_bdy=T), and 
     76      !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
     77      !! 
     78      !!              * Apply the time filter applied and swap of the dynamics 
     79      !!             arrays to start the next time step: 
     80      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 
     81      !!                (un,vn) = (ua,va). 
     82      !!             Note that with flux form advection and variable volume layer 
     83      !!             (lk_vvl=T), the time filter is applied on thickness weighted 
     84      !!             velocity. 
     85      !! 
     86      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step 
     87      !!               un,vn   now horizontal velocity of next time-step 
    7788      !!---------------------------------------------------------------------- 
    7889      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    … …  
    8596      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         -  
    8697      REAL(wp) ::   zuf   , zvf              !    -         -  
    87      !!---------------------------------------------------------------------- 
     98      !!---------------------------------------------------------------------- 
    8899 
    89100      IF( kt == nit000 ) THEN 
    … …  
    93104      ENDIF 
    94105 
    95       ! Local constant initialization 
    96       z2dt = 2. * rdt 
     106#if defined key_dynspg_flt 
     107      ! 
     108      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine 
     109      ! ------------- 
     110 
     111      ! Update after velocity on domain lateral boundaries      (only local domain required) 
     112      ! -------------------------------------------------- 
     113      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries 
     114      CALL lbc_lnk( va, 'V', -1. )  
     115      ! 
     116#else 
     117      ! Next velocity :   Leap-frog time stepping 
     118      ! ------------- 
     119      z2dt = 2. * rdt                                 ! Euler or leap-frog time step  
    97120      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt 
    98  
    99       !! Explicit physics with thickness weighted updates 
    100  
    101       ! Lateral boundary conditions on ( ua, va ) 
    102       CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )  
    103  
    104       ! Next velocity 
    105       ! ------------- 
    106 #if defined key_dynspg_flt 
    107       ! Leap-frog time stepping already done in dynspg_flt.F routine 
    108 #else 
    109       IF( lk_vvl ) THEN                                  ! Varying levels 
     121      ! 
     122      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
    110123         DO jk = 1, jpkm1 
    111             ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)     & 
    112                &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) )   & 
    113                &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
    114             va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)     & 
    115                &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) )   & 
    116                &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
    117          END DO 
    118       ELSE 
    119          DO jk = 1, jpkm1 
    120                   ! Leap-frog time stepping 
    121124            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    122125            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    123126         END DO 
    124       ENDIF 
    125  
     127      ELSE                                            ! applied on thickness weighted velocity 
     128         DO jk = 1, jpkm1 
     129            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
     130               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
     131               &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
     132            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
     133               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
     134               &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
     135         END DO 
     136      ENDIF 
     137 
     138 
     139      ! Update after velocity on domain lateral boundaries 
     140      ! --------------------------------------------------       
     141      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries 
     142      CALL lbc_lnk( va, 'V', -1. )  
     143      ! 
    126144# if defined key_obc 
    127       ! Update (ua,va) along open boundaries (only in the rigid-lid case) 
     145      !                                !* OBC open boundaries 
    128146      CALL obc_dyn( kt ) 
    129  
     147      ! 
    130148      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 
    131          ! Flather boundary condition : 
    132          !        - Update sea surface height on each open boundary 
    133          !                 sshn (= after ssh) for explicit case 
    134          !                 sshn_b (= after ssha_b) for time-splitting case 
    135          !        - Correct the barotropic velocities 
     149         ! Flather boundary condition : - Update sea surface height on each open boundary 
     150         !                                       sshn   (= after ssh   ) for explicit case 
     151         !                                       sshn_b (= after ssha_b) for time-splitting case 
     152         !                              - Correct the barotropic velocities 
    136153         CALL obc_dyn_bt( kt ) 
    137154         ! 
    138          ! Boundary conditions on sshn ( after ssh) 
    139          CALL lbc_lnk( sshn, 'T', 1. ) 
     155!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called 
     156         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn 
    140157         ! 
    141158         IF( ln_vol_cst )   CALL obc_vol( kt ) 
    … …  
    143160         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask ) 
    144161      ENDIF 
    145  
     162      ! 
    146163# elif defined key_bdy  
    147       ! Update (ua,va) along open boundaries (for exp or ts options). 
    148       IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 
     164      !                                !* BDY open boundaries 
     165      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option 
    149166         ! 
    150167         CALL bdy_dyn_frs( kt ) 
    151168         ! 
    152          IF( ln_bdy_fla ) THEN 
     169         IF( ln_bdy_dyn_fla ) THEN 
    153170            ua_e(:,:) = 0.e0 
    154171            va_e(:,:) = 0.e0 
    155172            ! Set these variables for use in bdy_dyn_fla 
    156             hu_e(:,:) = hu(:,:) 
    157             hv_e(:,:) = hv(:,:) 
     173            hur_e(:,:) = hur(:,:) 
     174            hvr_e(:,:) = hvr(:,:) 
    158175            DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    159176               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    160177               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    161178            END DO 
     179            ua_e(:,:) = ua_e(:,:) * hur(:,:) 
     180            va_e(:,:) = va_e(:,:) * hvr(:,:) 
    162181            DO jk = 1 , jpkm1 
    163                ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:) 
    164                va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:) 
     182               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) 
     183               va(:,:,jk) = va(:,:,jk) - va_e(:,:) 
    165184            END DO 
    166185            CALL bdy_dta_bt( kt+1, 0) 
    167             CALL bdy_dyn_fla 
     186            CALL bdy_dyn_fla( sshn_b ) 
     187            DO jk = 1 , jpkm1 
     188               ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk) 
     189               va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk) 
     190            END DO 
    168191         ENDIF 
    169192         ! 
    170          DO jk = 1 , jpkm1 
    171             ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:) 
    172             va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:) 
    173          END DO 
    174193      ENDIF 
    175194# endif 
    176  
     195      ! 
    177196# if defined key_agrif 
    178       CALL Agrif_dyn( kt ) 
     197      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries 
    179198# endif 
    180199#endif 
    … …  
    182201      ! Time filter and swap of dynamics arrays 
    183202      ! ------------------------------------------ 
    184       IF( neuler == 0 .AND. kt == nit000 ) THEN 
    185          DO jk = 1, jpkm1                                  
    186             un(:,:,jk) = ua(:,:,jk) 
     203      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap 
     204         DO jk = 1, jpkm1 
     205            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua 
    187206            vn(:,:,jk) = va(:,:,jk) 
    188207         END DO 
    189       ELSE 
    190          IF( lk_vvl ) THEN                                ! Varying levels 
    191             DO jk = 1, jpkm1                              ! filter applied on thickness weighted velocities 
     208      ELSE                                             !* Leap-Frog : Asselin filter and swap 
     209         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity 
     210            DO jk = 1, jpkm1                               
     211               DO jj = 1, jpj 
     212                  DO ji = 1, jpi     
     213                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 
     214                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 
     215                     ! 
     216                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     217                     vb(ji,jj,jk) = zvf 
     218                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     219                     vn(ji,jj,jk) = va(ji,jj,jk) 
     220                  END DO 
     221               END DO 
     222            END DO 
     223         ELSE                                                ! applied on thickness weighted velocity 
     224            DO jk = 1, jpkm1 
    192225               DO jj = 1, jpj 
    193226                  DO ji = 1, jpi 
    … …  
    206239                     zve3b = vb(ji,jj,jk) * ze3v_b 
    207240                     ! 
    208                      ub(ji,jj,jk) = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   & 
    209                         &         / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 
    210                      vb(ji,jj,jk) = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   & 
    211                         &         / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 
    212                      un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk) 
    213                      vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk) 
    214                   END DO 
    215                END DO 
    216             END DO 
    217          ELSE                                             ! Fixed levels 
    218             DO jk = 1, jpkm1                              ! filter applied on velocities 
    219                DO jj = 1, jpj 
    220                   DO ji = 1, jpi 
    221                      zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 
    222                      zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 
    223                      ub(ji,jj,jk) = zuf 
     241                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   & 
     242                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 
     243                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   & 
     244                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 
     245                     ! 
     246                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    224247                     vb(ji,jj,jk) = zvf 
    225                      un(ji,jj,jk) = ua(ji,jj,jk) 
     248                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
    226249                     vn(ji,jj,jk) = va(ji,jj,jk) 
    227250                  END DO 
    … …  
    232255 
    233256#if defined key_agrif 
     257      ! Update velocity at AGRIF zoom boundaries 
    234258      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt ) 
    235259#endif       
    … …  
    240264   END SUBROUTINE dyn_nxt 
    241265 
    242    !!====================================================================== 
     266   !!========================================================================= 
    243267END MODULE dynnxt 

  • trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90

    r1438 r1502  
    4040 
    4141#if   defined key_dynspg_ts   ||  defined key_vvl   ||  defined key_esopa 
    42    !! Time splitting variables 
    43    !! ------------------------ 
    44       REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop 
    45          un_e  , vn_e,                        & ! vertically integrated horizontal velocities (now) 
    46          ua_e  , va_e                           ! vertically integrated horizontal velocities (after) 
    47       REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop 
    48          sshn_e, ssha_e,                      & ! sea surface heigth (now, after) 
    49          hu_e  , hv_e                           ! depth arrays for the barotropic solution 
     42  !! Time splitting variables   ! variables of the explicit barotropic loop 
     43  !! ------------------------ 
     44     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ua_e  , va_e             ! barotropic velocities (after) 
     45     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshn_e, ssha_e, sshn_b   ! sea surface heigth (now, after, average) 
     46     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hu_e  , hv_e             ! depth arrays for the barotropic solution 
     47     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hur_e , hvr_e            ! inverse of depth arrays  
    5048#endif 
    5149 

  • trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r1438 r1502  
    11MODULE dynspg_ts 
    22   !!====================================================================== 
    3    !! History :   1.0  !  04-12  (L. Bessieres, G. Madec)  Original code 
    4    !!              -   !  05-11  (V. Garnier, G. Madec)  optimization 
    5    !!              -   !  06-08  (S. Masson)  distributed restart using iom 
    6    !!             2.0  !  07-07  (D. Storkey) calls to BDY routines 
    7    !!              -   !  08-01  (R. Benshila)  change averaging method 
     3   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code 
     4   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization 
     5   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom 
     6   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines 
     7   !!              -   ! 2008-01  (R. Benshila)  change averaging method 
     8   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation 
    89  !!--------------------------------------------------------------------- 
    910#if defined key_dynspg_ts   ||   defined key_esopa 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   'key_dynspg_ts'         free surface cst volume with time splitting 
    12    !!---------------------------------------------------------------------- 
    1313   !!---------------------------------------------------------------------- 
    1414   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
    … …  
    4545   PUBLIC ts_rst      ! routine called by istate.F90 
    4646 
     47 
    4748   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter 
    4849   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme) 
     50 
     51   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! averaged velocity 
    4952 
    5053   !! * Substitutions 
    … …  
    6669      !!      gradient in case of free surface formulation with time-splitting. 
    6770      !!      Add it to the general trend of momentum equation. 
    68       !!      Compute the free surface. 
    6971      !! 
    7072      !! ** Method  :   Free surface formulation with time-splitting 
    … …  
    7577      !!          the barotropic integration. 
    7678      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and  
    77       !!          barotropic transports (ua_e and va_e) through barotropic  
     79      !!          barotropic velocity (ua_e and va_e) through barotropic  
    7880      !!          momentum and continuity integration. Barotropic former  
    7981      !!          variables are time averaging over the full barotropic cycle 
    80       !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b  
     82      !!          (= 2 * baroclinic time step) and saved in zuX_b  
    8183      !!          and zvX_b (X specifying after, now or before). 
    8284      !!      -3- The new general trend becomes : 
    83       !!          ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H 
     85      !!          ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) ) 
    8486      !! 
    8587      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
    … …  
    8789      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 
    8890      !!--------------------------------------------------------------------- 
    89       INTEGER, INTENT( in )  ::   kt           ! ocean time-step index 
    90       !! 
    91       INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices 
    92       INTEGER  ::  icycle                      ! temporary scalar 
    93       REAL(wp) ::                           & 
    94          zraur, zcoef, z2dt_e, z2dt_b, zfac25,   &  ! temporary scalars 
    95          zfact1, zspgu, zcubt, zx1, zy1,    &  !     "        " 
    96          zfact2, zspgv, zcvbt, zx2, zy2        !     "        " 
    97       REAL(wp), DIMENSION(jpi,jpj) ::       & 
    98          zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays 
    99          zua, zva, zub, zvb,                &  !     "        " 
    100          zua_e, zva_e, zssha_e,             &  !     "        " 
    101          zub_e, zvb_e, zsshb_e,             &  !     "        " 
    102          zu_sum, zv_sum, zssh_sum 
    103       !! Variable volume 
    104       REAL(wp), DIMENSION(jpi,jpj)     ::   & 
    105          zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace 
     91      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
     92      !! 
     93      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices 
     94      INTEGER  ::   icycle                      ! temporary scalar 
     95      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b,   &  ! temporary scalars 
     96         z1_8, zspgu, zcubt, zx1, zy1,    &  !     "        " 
     97         z1_4, zspgv, zcvbt, zx2, zy2        !     "        " 
     98      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e 
     99      !!  
     100      REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace 
     101      !! 
     102      REAL(wp), DIMENSION(jpi,jpj) ::   zcu, zwx, zua, zun, zub   ! 2D workspace 
     103      REAL(wp), DIMENSION(jpi,jpj) ::   zcv, zwy, zva, zvn, zvb   !  -      - 
     104      REAL(wp), DIMENSION(jpi,jpj) ::   zun_e, zub_e, zu_sum      ! 2D workspace 
     105      REAL(wp), DIMENSION(jpi,jpj) ::   zvn_e, zvb_e, zv_sum      !  -      - 
     106      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_sum                  !  -      - 
    106107      !!---------------------------------------------------------------------- 
    107108 
    108       ! Arrays initialization 
    109       ! --------------------- 
    110       zhdiv(:,:) = 0.e0 
    111  
    112  
    113       IF( kt == nit000 ) THEN 
     109      IF( kt == nit000 ) THEN             !* initialisation 
    114110         ! 
    115111         IF(lwp) WRITE(numout,*) 
    … …  
    117113         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting' 
    118114         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro 
    119  
     115         ! 
    120116         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: 
    121          !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 
    122  
    123          zssha_e(:,:) = sshn(:,:) 
    124          zua_e  (:,:) = un_e(:,:) 
    125          zva_e  (:,:) = vn_e(:,:) 
    126          hu_e   (:,:) = hu  (:,:) 
    127          hv_e   (:,:) = hv  (:,:) 
     117         !                               ! un_b, vn_b   
     118         ! 
     119         ua_e  (:,:) = un_b (:,:) 
     120         va_e  (:,:) = vn_b (:,:) 
     121         hu_e  (:,:) = hu   (:,:) 
     122         hv_e  (:,:) = hv   (:,:) 
     123         hur_e (:,:) = hur  (:,:) 
     124         hvr_e (:,:) = hvr  (:,:) 
    128125         IF( ln_dynvor_een ) THEN 
    129126            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0 
    … …  
    140137      ENDIF 
    141138 
    142       ! Substract the surface pressure gradient from the momentum 
    143       ! --------------------------------------------------------- 
    144       IF( lk_vvl ) THEN    ! Variable volume 
    145  
    146          ! 0. Local initialization 
    147          ! ----------------------- 
    148          zspgu_1(:,:) = 0.e0 
    149          zspgv_1(:,:) = 0.e0 
    150  
    151          ! 1. Surface pressure gradient (now) 
    152          ! ---------------------------- 
    153          DO jj = 2, jpjm1 
    154             DO ji = fs_2, fs_jpim1   ! vector opt. 
    155                zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  ) & 
    156                   &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e1u(ji,jj) 
    157                zspgv_1(ji,jj) = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1) & 
    158                   &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e2v(ji,jj) 
    159             END DO 
    160          END DO 
    161  
    162          ! 2. Substract the surface pressure trend from the general trend 
    163          ! ------------------------------------------------------ 
    164          DO jk = 1, jpkm1 
    165             DO jj = 2, jpjm1 
    166                DO ji = fs_2, fs_jpim1   ! vector opt. 
    167                   ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj) 
    168                   va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj) 
    169                END DO 
    170             END DO 
    171          END DO 
    172  
    173       ENDIF 
    174  
    175       ! Local constant initialization 
    176       ! -------------------------------- 
     139      !                                   !* Local constant initialization 
    177140      z2dt_b = 2.0 * rdt                                    ! baroclinic time step 
    178       IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt 
    179       zfact1 = 0.5 * 0.25                                   ! coefficient for vorticity estimates 
    180       zfact2 = 0.5 * 0.5 
     141      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates 
     142      z1_4 = 0.5 * 0.5 
    181143      zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water 
     144      ! 
     145      zhdiv(:,:) = 0.e0                                     ! barotropic divergence 
    182146 
    183147      ! ----------------------------------------------------------------------------- 
    184148      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    185149      ! ----------------------------------------------------------------------------- 
    186  
    187       ! Vertically integrated quantities 
    188       ! -------------------------------- 
    189       zua(:,:) = 0.e0 
    190       zva(:,:) = 0.e0 
    191       zub(:,:) = 0.e0 
    192       zvb(:,:) = 0.e0 
    193       zwx(:,:) = 0.e0 
    194       zwy(:,:) = 0.e0 
    195  
    196       ! vertical sum 
    197       IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    198          DO jk = 1, jpkm1 
     150      !       
     151      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
     152      !                                   ! -------------------------- 
     153      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   zub(:,:) = 0.e0 
     154      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   zvb(:,:) = 0.e0 
     155      ! 
     156      DO jk = 1, jpkm1 
     157#if defined key_vectopt_loop 
     158         DO jj = 1, 1         !Vector opt. => forced unrolling 
    199159            DO ji = 1, jpij 
    200                !                                                           ! Vertically integrated momentum trends 
    201                zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk) 
    202                zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 
    203                !                                                           ! Vertically integrated transports (before) 
    204                zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 
    205                zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 
    206                !                                                           ! Planetary vorticity transport fluxes (now) 
    207                zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) 
    208                zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk) 
    209             END DO 
    210          END DO 
    211       ELSE                             ! No  vector opt. 
    212          DO jk = 1, jpkm1 
    213             !                                                           ! Vertically integrated momentum trends 
    214             zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    215             zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    216             !                                                           ! Vertically integrated transports (before) 
    217             zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 
    218             zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 
    219             !                                                           ! Planetary vorticity (now) 
    220             zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    221             zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    222          END DO 
    223       ENDIF 
    224  
     160#else  
     161         DO jj = 1, jpj 
     162            DO ji = 1, jpi 
     163#endif 
     164               !                                                                              ! now trend 
     165               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
     166               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
     167               !                                                                              ! now velocity  
     168               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk) 
     169               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)                
     170               !                                                                              ! before velocity 
     171               zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  
     172               zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) 
     173            END DO 
     174         END DO 
     175      END DO 
     176 
     177      !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
     178      DO jk = 1, jpkm1                    ! -------------------------- 
     179         DO jj = 2, jpjm1 
     180            DO ji = fs_2, fs_jpim1   ! vector opt. 
     181               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
     182               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
     183            END DO 
     184         END DO 
     185      END DO 
     186 
     187      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent) 
     188      !                                   ! ---------------------------==== 
     189      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport  
     190      zwy(:,:) = zvn(:,:) * e1v(:,:) 
     191      ! 
    225192      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
    226193         DO jj = 2, jpjm1 
    … …  
    231198               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    232199               ! energy conserving formulation for planetary vorticity term 
    233                zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
    234                zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
     200               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
     201               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    235202            END DO 
    236203         END DO 
    … …  
    239206         DO jj = 2, jpjm1 
    240207            DO ji = fs_2, fs_jpim1   ! vector opt. 
    241                zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
    242                               + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    243                zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
    244                               + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     208               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     209               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    245210               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    246211               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    … …  
    249214         ! 
    250215      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme 
    251          zfac25 = 0.25 
    252216         DO jj = 2, jpjm1 
    253217            DO ji = fs_2, fs_jpim1   ! vector opt. 
    254                zcu(ji,jj) = + zfac25 / e1u(ji,jj)   & 
    255                   &       * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    256                   &          + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    257                zcv(ji,jj) = - zfac25 / e2v(ji,jj)   & 
    258                   &       * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    259                   &          + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     218               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     219                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     220               zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     221                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    260222            END DO 
    261223         END DO 
    … …  
    263225      ENDIF 
    264226 
    265  
    266       ! Remove barotropic trend from general momentum trend 
    267       ! --------------------------------------------------- 
    268       DO jk = 1 , jpkm1 
    269          DO jj = 2, jpjm1 
    270             DO ji = fs_2, fs_jpim1   ! vector opt. 
    271                ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
    272                va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
    273             END DO 
    274          END DO 
    275       END DO 
    276  
    277       ! Remove coriolis term from barotropic trend 
    278       ! ------------------------------------------ 
    279       DO jj = 2, jpjm1 
     227      !                                   !* Right-Hand-Side of the barotropic momentum equation 
     228      !                                   ! ---------------------------------------------------- 
     229      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient 
     230         DO jj = 2, jpjm1  
     231            DO ji = fs_2, fs_jpim1   ! vector opt. 
     232               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    & 
     233                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
     234               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    & 
     235                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
     236            END DO 
     237         END DO 
     238      ENDIF 
     239 
     240      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
    280241         DO ji = fs_2, fs_jpim1 
    281242            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    … …  
    284245      END DO 
    285246 
     247      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
     248      !                                   ! --------------------------  
     249      zua(:,:) = zua(:,:) * hur(:,:) 
     250      zva(:,:) = zva(:,:) * hvr(:,:) 
     251      ! 
     252      IF( lk_vvl ) THEN 
     253         zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
     254         zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
     255      ELSE 
     256         zub(:,:) = zub(:,:) * hur(:,:) 
     257         zvb(:,:) = zvb(:,:) * hvr(:,:) 
     258      ENDIF 
     259 
    286260      ! ----------------------------------------------------------------------- 
    287261      !  Phase 2 : Integration of the barotropic equations with time splitting 
    288262      ! ----------------------------------------------------------------------- 
    289  
    290       ! Initialisations 
    291       !---------------- 
    292       ! Number of iteration of the barotropic loop 
    293       icycle = 2  * nn_baro + 1 
    294  
    295       ! variables for the barotropic equations 
    296       zu_sum  (:,:) = 0.e0 
    297       zv_sum  (:,:) = 0.e0 
    298       zssh_sum(:,:) = 0.e0 
    299       hu_e    (:,:) = hu    (:,:)      ! (barotropic) ocean depth at u-point 
    300       hv_e    (:,:) = hv    (:,:)      ! (barotropic) ocean depth at v-point 
    301       zsshb_e (:,:) = sshn_e(:,:)      ! (barotropic) sea surface height (before and now) 
    302       ! vertical sum 
    303       un_e  (:,:) = 0.e0 
    304       vn_e  (:,:) = 0.e0 
    305       IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    306          DO jk = 1, jpkm1 
    307             DO ji = 1, jpij 
    308                un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
    309                vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
    310             END DO 
    311          END DO 
    312       ELSE                             ! No  vector opt. 
    313          DO jk = 1, jpkm1 
    314             un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    315             vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    316          END DO 
    317       ENDIF 
    318       zub_e (:,:) = un_e(:,:) 
    319       zvb_e (:,:) = vn_e(:,:) 
    320  
    321  
     263      ! 
     264      !                                             ! ==================== ! 
     265      !                                             !    Initialisations   ! 
     266      !                                             ! ==================== ! 
     267      icycle = 2  * nn_baro            ! Number of barotropic sub time-step 
     268       
     269      !                                ! Start from NOW field 
     270      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points 
     271      hv_e   (:,:) = hv   (:,:)  
     272      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points 
     273      hvr_e  (:,:) = hvr  (:,:) 
     274      zsshb_e(:,:) = sshn (:,:)            ! sea surface height (before and now) 
     275      sshn_e (:,:) = sshn (:,:) 
     276       
     277      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external) 
     278      zvn_e  (:,:) = vn_b (:,:) 
     279      zub_e  (:,:) = un_b (:,:) 
     280      zvb_e  (:,:) = vn_b (:,:) 
     281 
     282      zu_sum  (:,:) = un_b (:,:)           ! summation 
     283      zv_sum  (:,:) = vn_b (:,:) 
     284      zssh_sum(:,:) = sshn (:,:) 
     285 
     286#if defined key_obc 
    322287      ! set ssh corrections to 0 
    323288      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 
    324 #if defined key_obc 
    325289      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0 
    326290      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0 
    … …  
    329293#endif 
    330294 
    331       ! Barotropic integration over 2 baroclinic time steps 
    332       ! --------------------------------------------------- 
    333  
    334       !                                                    ! ==================== ! 
    335       DO jit = 1, icycle                                   !  sub-time-step loop  ! 
    336          !                                                 ! ==================== ! 
     295      !                                             ! ==================== ! 
     296      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1) 
     297         !                                          ! ==================== ! 
    337298         z2dt_e = 2. * ( rdt / nn_baro ) 
    338          IF ( jit == 1 )   z2dt_e = rdt / nn_baro 
    339  
    340          ! Time interpolation of open boundary condition data 
    341          IF( lk_obc )   CALL obc_dta_bt( kt, jit ) 
    342          IF( lk_bdy .OR. ln_bdy_tides)   CALL bdy_dta_bt( kt, jit+1 ) 
    343  
    344          ! Horizontal divergence of barotropic transports 
    345          !-------------------------------------------------- 
    346          zhdiv(:,:) = 0.e0 
    347          DO jj = 2, jpjm1 
    348             DO ji = fs_2, fs_jpim1   ! vector opt. 
    349                zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * un_e(ji  ,jj)              & 
    350                   &            -e2u(ji-1,jj  ) * un_e(ji-1,jj)              & 
    351                   &            +e1v(ji  ,jj  ) * vn_e(ji  ,jj)              & 
    352                   &            -e1v(ji  ,jj-1) * vn_e(ji  ,jj-1) )          & 
    353                   &           / (e1t(ji,jj)*e2t(ji,jj)) 
    354             END DO 
    355          END DO 
    356  
     299         IF( jn == 1 )   z2dt_e = rdt / nn_baro 
     300 
     301         !                                                !* Update the forcing (OBC, BDY and tides) 
     302         !                                                !  ------------------ 
     303         IF( lk_obc                     )   CALL obc_dta_bt( kt, jn   ) 
     304         IF( lk_bdy  .OR.  ln_bdy_tides )   CALL bdy_dta_bt( kt, jn+1 ) 
     305 
     306         !                                                !* after ssh_e 
     307         !                                                !  ----------- 
     308         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports 
     309            DO ji = fs_2, fs_jpim1   ! vector opt. 
     310               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
     311                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     & 
     312                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     & 
     313                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     314            END DO 
     315         END DO 
     316         ! 
    357317#if defined key_obc 
    358          ! open boundaries (div must be zero behind the open boundary) 
    359          !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    360          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)  = 0.e0      ! east 
    361          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)  = 0.e0      ! west 
     318         !                                                     ! OBC : zhdiv must be zero behind the open boundary 
     319!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
     320         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east 
     321         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west 
    362322         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north 
    363          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)  = 0.e0      ! south 
     323         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south 
    364324#endif 
    365  
    366325#if defined key_bdy 
    367          DO jj = 1, jpj 
    368             DO ji = 1, jpi 
    369                zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 
    370             END DO 
    371          END DO 
     326         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask 
    372327#endif 
    373  
    374          ! Sea surface height from the barotropic system 
    375          !---------------------------------------------- 
    376          DO jj = 2, jpjm1 
    377             DO ji = fs_2, fs_jpim1   ! vector opt. 
    378                zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  & 
    379                   &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
    380             END DO 
    381          END DO 
    382  
    383          ! evolution of the barotropic transport ( following the vorticity scheme used) 
    384          ! ---------------------------------------------------------------------------- 
    385          zwx(:,:) = e2u(:,:) * un_e(:,:) 
    386          zwy(:,:) = e1v(:,:) * vn_e(:,:) 
    387  
    388          IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
     328         ! 
     329         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
     330            DO ji = fs_2, fs_jpim1   ! vector opt. 
     331               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
     332            END DO 
     333         END DO 
     334 
     335         !                                                !* after barotropic velocities (vorticity scheme dependent) 
     336         !                                                !  ---------------------------   
     337         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport 
     338         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
     339         ! 
     340         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==! 
    389341            DO jj = 2, jpjm1 
    390342               DO ji = fs_2, fs_jpim1   ! vector opt. 
    391343                  ! surface pressure gradient 
    392344                  IF( lk_vvl) THEN 
    393                      zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) & 
    394                         &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj) 
    395                      zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) & 
    396                         &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj) 
     345                     zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
     346                        &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
     347                     zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
     348                        &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    397349                  ELSE 
    398                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 
    399                      zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 
     350                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     351                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    400352                  ENDIF 
    401353                  ! energy conserving formulation for planetary vorticity term 
    … …  
    404356                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    405357                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    406                   zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
    407                   zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    408                   ! after transports 
    409                   zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    410                   zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     358                  zcubt = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 
     359                  zcvbt =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
     360                  ! after barotropic velocity 
     361                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
     362                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
    411363               END DO 
    412364            END DO 
    413365            ! 
    414          ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme 
     366         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==! 
     367            DO jj = 2, jpjm1 
     368               DO ji = fs_2, fs_jpim1   ! vector opt. 
     369                   ! surface pressure gradient 
     370                  IF( lk_vvl) THEN 
     371                     zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
     372                        &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
     373                     zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
     374                        &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
     375                  ELSE 
     376                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     377                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
     378                  ENDIF 
     379                  ! enstrophy conserving formulation for planetary vorticity term 
     380                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     381                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     382                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 
     383                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
     384                  ! after barotropic velocity 
     385                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
     386                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     387               END DO 
     388            END DO 
     389            ! 
     390         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==! 
    415391            DO jj = 2, jpjm1 
    416392               DO ji = fs_2, fs_jpim1   ! vector opt. 
    417393                  ! surface pressure gradient 
    418394                  IF( lk_vvl) THEN 
    419                      zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) & 
    420                         &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj) 
    421                      zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) & 
    422                         &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj) 
     395                     zspgu = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
     396                        &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
     397                     zspgv = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
     398                        &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    423399                  ELSE 
    424                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 
    425                      zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 
    426                   ENDIF 
    427                   ! enstrophy conserving formulation for planetary vorticity term 
    428                   zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
    429                      + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    430                   zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
    431                      + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    432                   zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    433                   zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    434                   ! after transports 
    435                   zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    436                   zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
    437                END DO 
    438             END DO 
    439             ! 
    440          ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme 
    441             zfac25 = 0.25 
    442             DO jj = 2, jpjm1 
    443                DO ji = fs_2, fs_jpim1   ! vector opt. 
    444                   ! surface pressure gradient 
    445                   IF( lk_vvl) THEN 
    446                      zspgu = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  ) & 
    447                         &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj) 
    448                      zspgv = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) & 
    449                         &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj) 
    450                   ELSE 
    451                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 
    452                      zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 
     400                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     401                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    453402                  ENDIF 
    454403                  ! energy/enstrophy conserving formulation for planetary vorticity term 
    455                   zcubt = + zfac25 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    456                      &                             + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    457                   zcvbt = - zfac25 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    458                      &                             + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    459                   ! after transports 
    460                   zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    461                   zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     404                  zcubt = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     405                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 
     406                  zcvbt = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     407                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
     408                  ! after barotropic velocity 
     409                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
     410                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
    462411               END DO 
    463412            END DO 
    464413            !  
    465414         ENDIF 
    466  
    467          ! Flather's boundary condition for the barotropic loop : 
    468          !         - Update sea surface height on each open boundary 
    469          !         - Correct the barotropic transports 
    470          IF( lk_obc )   CALL obc_fla_ts 
    471          IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla 
    472  
    473          ! ... Boundary conditions on ua_e, va_e, ssha_e 
    474          CALL lbc_lnk( zua_e  , 'U', -1. ) 
    475          CALL lbc_lnk( zva_e  , 'V', -1. ) 
    476          CALL lbc_lnk( zssha_e, 'T',  1. ) 
    477  
    478          ! temporal sum 
    479          !------------- 
    480          zu_sum  (:,:) = zu_sum  (:,:) + zua_e  (:,:) 
    481          zv_sum  (:,:) = zv_sum  (:,:) + zva_e  (:,:)  
    482          zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:)  
    483  
    484          ! Time filter and swap of dynamics arrays 
    485          ! --------------------------------------- 
    486          IF( jit == 1 ) THEN   ! Euler (forward) time stepping 
    487             zsshb_e(:,:) = sshn_e (:,:) 
    488             zub_e  (:,:) = un_e   (:,:) 
    489             zvb_e  (:,:) = vn_e   (:,:) 
    490             sshn_e (:,:) = zssha_e(:,:) 
    491             un_e   (:,:) = zua_e  (:,:) 
    492             vn_e   (:,:) = zva_e  (:,:) 
    493          ELSE                                        ! Asselin filtering 
    494             zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
    495             zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * un_e  (:,:) 
    496             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * vn_e  (:,:) 
    497             sshn_e (:,:) = zssha_e(:,:) 
    498             un_e   (:,:) = zua_e  (:,:) 
    499             vn_e   (:,:) = zva_e  (:,:) 
     415         !                                                !* domain lateral boundary 
     416         !                                                !  ----------------------- 
     417         !                                                      ! Flather's boundary condition for the barotropic loop : 
     418         !                                                      !         - Update sea surface height on each open boundary 
     419         !                                                      !         - Correct the velocity 
     420 
     421         IF( lk_obc                   )   CALL obc_fla_ts 
     422         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla( sshn_e )  
     423         ! 
     424         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
     425         CALL lbc_lnk( va_e  , 'V', -1. ) 
     426         CALL lbc_lnk( ssha_e, 'T',  1. ) 
     427 
     428         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps 
     429         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:)  
     430         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)  
     431 
     432         !                                                !* Time filter and swap 
     433         !                                                !  -------------------- 
     434         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step) 
     435            zsshb_e(:,:) = sshn_e(:,:) 
     436            zub_e  (:,:) = zun_e (:,:) 
     437            zvb_e  (:,:) = zvn_e (:,:) 
     438            sshn_e (:,:) = ssha_e(:,:) 
     439            zun_e  (:,:) = ua_e  (:,:) 
     440            zvn_e  (:,:) = va_e  (:,:) 
     441         ELSE                                                   ! Swap + Filter 
     442            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
     443            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:) 
     444            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:) 
     445            sshn_e (:,:) = ssha_e(:,:) 
     446            zun_e  (:,:) = ua_e  (:,:) 
     447            zvn_e  (:,:) = va_e  (:,:) 
    500448         ENDIF 
    501449 
    502          IF( lk_vvl ) THEN ! Variable volume   ! needed for BDY ??? 
    503  
    504             ! Sea surface elevation 
    505             ! --------------------- 
    506             DO jj = 1, jpjm1 
    507                DO ji = 1,jpim1 
    508  
    509                   ! Sea Surface Height at u-point before 
    510                   zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )   & 
    511                      &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
    512                      &                + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
    513  
    514                   ! Sea Surface Height at v-point before 
    515                   zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )   & 
    516                      &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
    517                      &                + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
    518  
    519                END DO 
    520             END DO 
    521  
    522             ! Boundaries conditions 
    523             CALL lbc_lnk( zsshun_e, 'U', 1. )    ;  CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
    524  
    525             ! Ocean depth at U- and V-points 
    526             ! ------------------------------------------- 
    527             hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 
    528             hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 
    529  
     450         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
     451            !                                             !  ------------------ 
     452            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points 
     453               DO ji = 1, fs_jpim1   ! Vector opt. 
     454                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     455                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
     456                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
     457                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     458                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
     459                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
     460               END DO 
     461            END DO 
     462            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions 
     463            CALL lbc_lnk( zsshvn_e, 'V', 1. )  
     464            ! 
     465            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points 
     466            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 
     467            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) ) 
     468            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) ) 
    530469            ! 
    531470         ENDIF 
    … …  
    534473      !                                                    ! ==================== ! 
    535474 
    536       ! Time average of after barotropic variables 
    537       zcoef =  1.e0 / ( 2 * nn_baro  + 1 )  
    538       un_e  (:,:) = zcoef *  zu_sum  (:,:)  
    539       vn_e  (:,:) = zcoef *  zv_sum  (:,:)  
    540       sshn_e(:,:) = zcoef *  zssh_sum(:,:)  
    541  
    542475#if defined key_obc 
    543       IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) 
     476      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ????? 
    544477      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:) 
    545478      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:) 
    … …  
    548481 
    549482      ! ----------------------------------------------------------------------------- 
    550       ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step) 
     483      ! Phase 3. update the general trend with the barotropic trend 
    551484      ! ----------------------------------------------------------------------------- 
    552  
    553  
    554  
    555       ! add time averaged barotropic coriolis and surface pressure gradient 
    556       ! terms to the general momentum trend 
    557       ! -------------------------------------------------------------------- 
     485      ! 
     486      !                                   !* Time average ==> after barotropic u, v, ssh 
     487      zcoef =  1.e0 / ( 2 * nn_baro  + 1 )  
     488      un_b  (:,:) = zcoef * zu_sum  (:,:)  
     489      vn_b  (:,:) = zcoef * zv_sum  (:,:)  
     490      sshn_b(:,:) = zcoef * zssh_sum(:,:)  
     491      !  
     492      !                                   !* update the general momentum trend 
    558493      DO jk=1,jpkm1 
    559          ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b 
    560          va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b 
     494         ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b 
     495         va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b 
    561496      END DO 
    562  
    563       ! write filtered free surface arrays in restart file 
    564       ! -------------------------------------------------- 
     497      ! 
     498      !                                   !* write time-spliting arrays in the restart 
    565499      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    566  
    567500      ! 
    568501   END SUBROUTINE dyn_spg_ts 
    … …  
    582515      ! 
    583516      IF( TRIM(cdrw) == 'READ' ) THEN 
    584          IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 
    585             CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) )   ! free surface and 
    586             CALL iom_get( numror, jpdom_autoglo, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued 
    587             CALL iom_get( numror, jpdom_autoglo, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop 
     517         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 
     518            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued 
     519            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    588520         ELSE 
    589             sshn_e(:,:) = sshn(:,:) 
    590             un_e  (:,:) = 0.e0 
    591             vn_e  (:,:) = 0.e0 
     521            un_b (:,:) = 0.e0 
     522            vn_b (:,:) = 0.e0 
    592523            ! vertical sum 
    593524            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    594525               DO jk = 1, jpkm1 
    595526                  DO ji = 1, jpij 
    596                      un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
    597                      vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
     527                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
     528                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
    598529                  END DO 
    599530               END DO 
    600531            ELSE                             ! No  vector opt. 
    601532               DO jk = 1, jpkm1 
    602                   un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    603                   vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
     533                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
     534                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    604535               END DO 
    605536            ENDIF 
     537            un_b (:,:) = un_b(:,:) * hur(:,:) 
     538            vn_b (:,:) = vn_b(:,:) * hvr(:,:) 
    606539         ENDIF 
    607540      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    608          CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) )   ! free surface and 
    609          CALL iom_rstput( kt, nitrst, numrow, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued 
    610          CALL iom_rstput( kt, nitrst, numrow, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop 
     541         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! external velocity issued 
     542         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    611543      ENDIF 
    612544      ! 
Note: See TracChangeset for help on using the changeset viewer.