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 4292 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Ignore:
Timestamp:
2013-11-20T17:28:04+01:00 (10 years ago)
Author:
cetlod
Message:

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

Location:
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4147 r4292  
    767767      DO jj = 2, jpjm1 
    768768        DO ji = 2, jpim1 
    769           zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 
    770           zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 
     769          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation 
     770          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation 
    771771        END DO 
    772772      END DO 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r3764 r4292  
    1717   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module 
    1818   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 
     19   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1920   !!------------------------------------------------------------------------- 
    2021   
    … …  
    4243   USE wrk_nemo        ! Memory Allocation 
    4344   USE prtctl          ! Print control 
     45   USE dynspg_ts       ! Barotropic velocities 
     46 
    4447#if defined key_agrif 
    4548   USE agrif_opa_interp 
    … …  
    103106      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars 
    104107      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      - 
     108      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva 
    105109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f  
    106110      !!---------------------------------------------------------------------- 
    … …  
    109113      ! 
    110114      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     115      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva ) 
    111116      ! 
    112117      IF( kt == nit000 ) THEN 
    … …  
    127132      ! 
    128133#else 
     134 
     135# if defined key_dynspg_exp 
    129136      ! Next velocity :   Leap-frog time stepping 
    130137      ! ------------- 
    … …  
    147154         END DO 
    148155      ENDIF 
    149  
     156# endif 
     157 
     158# if defined key_dynspg_ts 
     159      ! Ensure below that barotropic velocities match time splitting estimate 
     160      ! Compute actual transport and replace it with ts estimate at "after" time step 
     161      zua(:,:) = 0._wp 
     162      zva(:,:) = 0._wp 
     163      IF (lk_vvl) THEN 
     164         DO jk = 1, jpkm1 
     165            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     166            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)     
     167         END DO 
     168         DO jk = 1, jpkm1 
     169            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_e(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
     170            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_e(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
     171         END DO 
     172      ELSE 
     173         DO jk = 1, jpkm1 
     174            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     175            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)     
     176         END DO 
     177         DO jk = 1, jpkm1 
     178            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk) 
     179            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk) 
     180         END DO 
     181      ENDIF 
     182 
     183      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 
     184         ! Remove advective velocity from "now velocities"  
     185         ! prior to asselin filtering      
     186         ! In the forward case, this is done below after asselin filtering     
     187         DO jk = 1, jpkm1 
     188            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
     189            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     190         END DO   
     191      ENDIF 
     192# endif 
    150193 
    151194      ! Update after velocity on domain lateral boundaries 
    … …  
    194237            vn(:,:,jk) = va(:,:,jk) 
    195238         END DO 
     239         IF (lk_vvl) THEN 
     240            DO jk = 1, jpkm1 
     241               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     242               fse3u_b(:,:,jk) = fse3u_n(:,:,jk) 
     243               fse3v_b(:,:,jk) = fse3v_n(:,:,jk) 
     244            ENDDO 
     245         ENDIF 
    196246      ELSE                                             !* Leap-Frog : Asselin filter and swap 
    197247         !                                ! =============! 
    … …  
    201251               DO jj = 1, jpj 
    202252                  DO ji = 1, jpi     
    203                      zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    204                      zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 
     253                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
     254                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
    205255                     ! 
    206256                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    … …  
    214264         ELSE                             ! Variable volume ! 
    215265            !                             ! ================! 
     266            ! Before scale factor at t-points 
     267            ! (used as a now filtered scale factor until the swap) 
     268            ! ---------------------------------------------------- 
     269            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     270               ! Remove asselin filtering on thicknesses if forward time splitting 
     271                  fse3t_b(:,:,:) = fse3t_n(:,:,:) 
     272            ELSE 
     273               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 
     274               ! Add volume filter correction: compatibility with tracer advection scheme 
     275               ! => time filter + conservation correction (only at the first level) 
     276               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    216277            ! 
    217             DO jk = 1, jpkm1                 ! Before scale factor at t-points 
    218                fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
    219                   &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     & 
    220                   &                         - 2._wp * fse3t_n(:,:,jk)            ) 
    221             END DO 
    222             zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors 
    223             fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     278            ENDIF 
    224279            ! 
    225             IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation) 
    226                ! 
    227                !                                      ! before scale factors at u- & v-pts (computed from fse3t_b) 
    228                CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 
    229                ! 
    230                DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity 
    231                   DO jj = 1, jpj                      !                                                 -------- 
     280            IF( ln_dynadv_vec ) THEN 
     281               ! Before scale factor at (u/v)-points 
     282               ! ----------------------------------- 
     283               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 
     284               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 
     285               ! Leap-Frog - Asselin filter and swap: applied on velocity 
     286               ! ----------------------------------- 
     287               DO jk = 1, jpkm1 
     288                  DO jj = 1, jpj 
    232289                     DO ji = 1, jpi 
    233                         zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    234                         zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 
     290                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
     291                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
    235292                        ! 
    236293                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    … …  
    242299               END DO 
    243300               ! 
    244             ELSE                             ! flux form (thickness weighted calulation) 
    245                ! 
    246                CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b) 
    247                ! 
    248                DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:  
    249                   DO jj = 1, jpj                      !                   applied on thickness weighted velocity 
     301            ELSE 
     302               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 
     303               !------------------------------------------------ 
     304               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 
     305               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 
     306               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 
     307               ! -----------------------------------             =========================== 
     308               DO jk = 1, jpkm1 
     309                  DO jj = 1, jpj 
    250310                     DO ji = 1, jpi                   !                              --------------------------- 
    251311                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 
    … …  
    272332         ENDIF 
    273333         ! 
    274       ENDIF 
     334         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     335         ! Remove asselin filtering of barotropic velocities if forward time splitting 
     336         ! note that we replace barotropic velocities by advective velocities        
     337            zua(:,:) = 0._wp 
     338            zva(:,:) = 0._wp 
     339            IF (lk_vvl) THEN 
     340               DO jk = 1, jpkm1 
     341                  zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
     342                  zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     343               END DO 
     344            ELSE 
     345               DO jk = 1, jpkm1 
     346                  zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
     347                  zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     348               END DO 
     349            ENDIF 
     350            DO jk = 1, jpkm1 
     351               ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 
     352               vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
     353            END DO 
     354         ENDIF 
     355         ! 
     356      ENDIF ! neuler =/0 
    275357 
    276358      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    … …  
    278360      !  
    279361      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     362      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva ) 
    280363      ! 
    281364      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r4245 r4292  
    2323   USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine) 
    2424   USE dynadv         ! dynamics: vector invariant versus flux form 
     25   USE dynhpg, ONLY: ln_dynhpg_imp 
     26   USE sbctide 
     27   USE updtide 
    2528   USE trdmod         ! ocean dynamics trends 
    2629   USE trdmod_oce     ! ocean variables trends 
    … …  
    101104      ENDIF 
    102105 
    103       IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==! 
    104          zg_2 = grav * 0.5 
    105          DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
     106      IF(      ln_apr_dyn                                                &   ! atmos. pressure 
     107         .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting) 
     108         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice 
     109         ! 
     110         DO jj = 2, jpjm1 
    106111            DO ji = fs_2, fs_jpim1   ! vector opt. 
    107                spgu(ji,jj) =  zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    & 
    108                   &                   + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
    109                spgv(ji,jj) =  zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    & 
    110                   &                   + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
    111             END DO 
    112          END DO 
    113          DO jk = 1, jpkm1                          ! Add the apg to the general trend 
     112               spgu(ji,jj) = 0._wp 
     113               spgv(ji,jj) = 0._wp 
     114            END DO 
     115         END DO          
     116         ! 
     117         IF( ln_apr_dyn .AND. (.NOT. lk_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==! 
     118            zg_2 = grav * 0.5 
     119            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
     120               DO ji = fs_2, fs_jpim1   ! vector opt. 
     121                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    & 
     122                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     123                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    & 
     124                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     125               END DO 
     126            END DO 
     127         ENDIF 
     128         ! 
     129         !                                    !==  tide potential forcing term  ==! 
     130         IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case 
     131            ! 
     132            CALL upd_tide( kt )                      ! update tide potential 
     133            ! 
     134            DO jj = 2, jpjm1                         ! add tide potential forcing 
     135               DO ji = fs_2, fs_jpim1   ! vector opt. 
     136                  spgv(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     137                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     138               END DO  
     139            END DO 
     140         ENDIF 
     141         ! 
     142         IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==! 
     143            CALL wrk_alloc( jpi, jpj, zpice ) 
     144            !                                             
     145            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) 
     146            zgrau0r     = - grav * r1_rau0 
     147            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r 
     148            DO jj = 2, jpjm1 
     149               DO ji = fs_2, fs_jpim1   ! vector opt. 
     150                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
     151                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     152               END DO 
     153            END DO 
     154            ! 
     155            CALL wrk_dealloc( jpi, jpj, zpice )          
     156         ENDIF 
     157         ! 
     158         DO jk = 1, jpkm1                     !== Add all terms to the general trend 
    114159            DO jj = 2, jpjm1 
    115160               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    118163               END DO 
    119164            END DO 
    120          END DO 
    121       ENDIF 
    122  
    123       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: Pressure gradient due to snow-ice mass ==! 
    124          CALL wrk_alloc( jpi, jpj, zpice ) 
    125          !                                             
    126          zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) 
    127          zgrau0r     = - grav * r1_rau0 
    128          zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r 
    129          DO jj = 2, jpjm1 
    130             DO ji = fs_2, fs_jpim1   ! vector opt. 
    131                spgu(ji,jj) = ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    132                spgv(ji,jj) = ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
    133             END DO 
    134          END DO 
    135          DO jk = 1, jpkm1                             ! Add the surface pressure trend to the general trend 
    136             DO jj = 2, jpjm1 
    137                DO ji = fs_2, fs_jpim1   ! vector opt. 
    138                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    139                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
    140                END DO 
    141             END DO 
    142          END DO 
    143          ! 
    144          CALL wrk_dealloc( jpi, jpj, zpice ) 
    145       ENDIF 
    146  
     165         END DO          
     166      ENDIF 
    147167 
    148168      SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend 
    … …  
    209229      ENDIF 
    210230 
     231      IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 
     232      ! (do it now, to set nn_baro, used to allocate some arrays later on) 
    211233      !                        ! allocate dyn_spg arrays 
    212234      IF( lk_dynspg_ts ) THEN 
    … …  
    248270      ENDIF 
    249271 
    250       !                        ! Control of momentum formulation 
    251       IF( lk_dynspg_ts .AND. lk_vvl ) THEN 
    252          IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' ) 
     272      !               ! Control of hydrostatic pressure choice 
     273      IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 
     274         CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 
    253275      ENDIF 
    254276      ! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r3680 r4292  
    9191               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 
    9292            END DO  
    93          END DO  
     93         END DO 
     94         ! 
    9495         DO jk = 1, jpkm1                    ! Add it to the general trend 
    9596            DO jj = 2, jpjm1 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90

    r3294 r4292  
    3939   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur_e , hvr_e    ! inverse of hu_e and hv_e 
    4040   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_b           ! before field without time-filter 
     41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b     ! after  averaged velocities 
     42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b     ! now    averaged velocities 
     43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b     ! before averaged velocities 
     44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv ! Advection vel. at "now" barocl. step 
     45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b  ! Advection vel. at "now-0.5" barocl. step 
    4146 
    4247   !!---------------------------------------------------------------------- 
    … …  
    5358      ALLOCATE( sshn_e(jpi,jpj) , ua_e(jpi,jpj) , hu_e(jpi,jpj) , hur_e(jpi,jpj) ,      & 
    5459         &      ssha_e(jpi,jpj) , va_e(jpi,jpj) , hv_e(jpi,jpj) , hvr_e(jpi,jpj) ,      & 
     60         &      ub_b(jpi,jpj)   , vb_b(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj)  ,      & 
     61         &      ua_b(jpi,jpj)   , va_b(jpi,jpj)                                  ,      &  
     62         &      ub2_b(jpi,jpj)  , vb2_b(jpi,jpj)                                 ,      & 
     63         &      un_adv(jpi,jpj) , vn_adv(jpi,jpj)                                ,      & 
    5564         &      sshn_b(jpi,jpj)                                                  , STAT = dynspg_oce_alloc ) 
    5665         ! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r3680 r4292  
    99   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations 
    1010   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b 
     11   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
     12   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1113   !!--------------------------------------------------------------------- 
    1214#if defined key_dynspg_ts   ||   defined key_esopa 
    … …  
    1618   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
    1719   !!                 splitting scheme and add to the general trend  
    18    !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file 
    1920   !!---------------------------------------------------------------------- 
    2021   USE oce             ! ocean dynamics and tracers 
    … …  
    2425   USE phycst          ! physical constants 
    2526   USE domvvl          ! variable volume 
    26    USE zdfbfr          ! bottom friction 
    2727   USE dynvor          ! vorticity term 
    28    USE obc_oce         ! Lateral open boundary condition 
    29    USE obc_par         ! open boundary condition parameters 
    30    USE obcdta          ! open boundary condition data      
    31    USE obcfla          ! Flather open boundary condition   
    3228   USE bdy_par         ! for lk_bdy 
    3329   USE bdy_oce         ! Lateral open boundary condition 
    34    USE bdydta          ! open boundary condition data      
     30   USE bdytides        ! open boundary condition data      
    3531   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    36    USE sbctide 
    37    USE updtide 
     32   USE sbctide         ! tides 
     33   USE updtide         ! tide potential 
    3834   USE lib_mpp         ! distributed memory computing library 
    3935   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    4137   USE in_out_manager  ! I/O manager 
    4238   USE iom             ! IOM library 
     39   USE restart         ! only for lrst_oce 
    4340   USE zdf_oce         ! Vertical diffusion 
    4441   USE wrk_nemo        ! Memory Allocation 
    45    USE timing          ! Timing 
     42   USE timing          ! Timing     
     43   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     44   USE dynadv, ONLY: ln_dynadv_vec 
     45#if defined key_agrif 
     46   USE agrif_opa_interp ! agrif 
     47#endif 
    4648 
    4749 
    … …  
    4951   PRIVATE 
    5052 
    51    PUBLIC dyn_spg_ts        ! routine called by step.F90 
    52    PUBLIC ts_rst            ! routine called by istate.F90 
    53    PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90 
    54  
    55  
     53   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90  
     54   PUBLIC dyn_spg_ts_alloc  !    "      "     "    " 
     55   PUBLIC dyn_spg_ts_init   !    "      "     "    " 
     56 
     57   ! Potential namelist parameters below to be read in dyn_spg_ts_init 
     58   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping 
     59   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables 
     60   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically 
     61   INTEGER,  PRIVATE, PARAMETER :: nn_bt_flt=1            !: Filter choice 
     62   REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp      !: Max. courant number (used if ln_bt_nn_auto=T) 
     63   ! End namelist parameters 
     64 
     65   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 
     66   REAL(wp),SAVE :: rdtbt   ! Barotropic time step 
     67 
     68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 
     69                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables 
     70                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables 
     71 
     72   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points 
    5673   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter 
    5774   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme) 
    5875 
    59    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b   ! now    averaged velocity 
    60    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b   ! before averaged velocity 
     76   ! Arrays below are saved to allow testing of the "no time averaging" option 
     77   ! If this option is not retained, these could be replaced by temporary arrays 
     78   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 
     79                                                   ubb_e, ub_e,     & 
     80                                                   vbb_e, vb_e 
    6181 
    6282   !! * Substitutions 
    … …  
    6484#  include "vectopt_loop_substitute.h90" 
    6585   !!---------------------------------------------------------------------- 
    66    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    67    !! $Id$ 
     86   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
     87   !! $Id: dynspg_ts.F90 
    6888   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6989   !!---------------------------------------------------------------------- 
    … …  
    7494      !!                  ***  routine dyn_spg_ts_alloc  *** 
    7595      !!---------------------------------------------------------------------- 
    76       ALLOCATE( ftnw  (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) ,     & 
    77          &      ftsw  (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc ) 
    78          ! 
     96      INTEGER :: ierr(3) 
     97      !!---------------------------------------------------------------------- 
     98      ierr(:) = 0 
     99 
     100      ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 
     101         &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , & 
     102         &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , STAT= ierr(1) ) 
     103 
     104      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 
     105 
     106      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     107                             &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
     108 
     109      dyn_spg_ts_alloc = MAXVAL(ierr(:)) 
     110 
    79111      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc ) 
    80112      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') 
    … …  
    82114   END FUNCTION dyn_spg_ts_alloc 
    83115 
    84  
    85116   SUBROUTINE dyn_spg_ts( kt ) 
    86117      !!---------------------------------------------------------------------- 
    87       !!                  ***  routine dyn_spg_ts  *** 
    88118      !! 
    89       !! ** Purpose :   Compute the now trend due to the surface pressure 
    90       !!      gradient in case of free surface formulation with time-splitting. 
    91       !!      Add it to the general trend of momentum equation. 
     119      !! ** Purpose :    
     120      !!      -Compute the now trend due to the explicit time stepping 
     121      !!      of the quasi-linear barotropic system. Barotropic variables are 
     122      !!      advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F) 
     123      !!      or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a 
     124      !!      generalized forward-backward (see ref. below) time stepping. 
     125      !!      -Update the free surface at step "n+1" (ssha, zsshu_a, zsshv_a). 
     126      !!      -Compute barotropic advective velocities at step "n" to be used  
     127      !!      to advect tracers latter on. These are compliant with discrete 
     128      !!      continuity equation taken at the baroclinic time steps, thus  
     129      !!      ensuring tracers conservation. 
    92130      !! 
    93       !! ** Method  :   Free surface formulation with time-splitting 
    94       !!      -1- Save the vertically integrated trend. This general trend is 
    95       !!          held constant over the barotropic integration. 
    96       !!          The Coriolis force is removed from the general trend as the 
    97       !!          surface gradient and the Coriolis force are updated within 
    98       !!          the barotropic integration. 
    99       !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and  
    100       !!          barotropic velocity (ua_e and va_e) through barotropic  
    101       !!          momentum and continuity integration. Barotropic former  
    102       !!          variables are time averaging over the full barotropic cycle 
    103       !!          (= 2 * baroclinic time step) and saved in uX_b  
    104       !!          and vX_b (X specifying after, now or before). 
    105       !!      -3- The new general trend becomes : 
    106       !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b ) 
     131      !! ** Method  :   
    107132      !! 
    108       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
     133      !! ** Action : - Update barotropic velocities: ua_b, va_b 
     134      !!             - Update trend (ua,va) with barotropic component 
     135      !!             - Update ssha, zsshu_a, zsshv_a 
     136      !!             - Update barotropic advective velocity at kt=now 
    109137      !! 
    110       !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 
     138      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:  
     139      !!              The regional oceanic modeling system (ROMS):  
     140      !!              a split-explicit, free-surface, 
     141      !!              topography-following-coordinate oceanic model.  
     142      !!              Ocean Modelling, 9, 347-404.  
    111143      !!--------------------------------------------------------------------- 
    112144      ! 
    113145      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    114146      ! 
    115       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    116       INTEGER  ::   icycle           ! local scalar 
    117       INTEGER  ::   ikbu, ikbv       ! local scalar 
    118       REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars 
    119       REAL(wp) ::   z1_8, zx1, zy1                            !   -      - 
    120       REAL(wp) ::   z1_4, zx2, zy2                            !   -      - 
    121       REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
    122       REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
    123       REAL(wp) ::   ua_btm, va_btm                            !   -      - 
    124       ! 
    125       REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv  
    126       REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e  
    127       REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 
     147      LOGICAL  ::   ll_fw_start        ! if true, forward integration  
     148      LOGICAL  ::   ll_init         ! if true, special startup of 2d equations 
     149      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
     150      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     151      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
     152      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     153      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      - 
     154      REAL(wp) ::   zu_spg, zv_spg     !   -      - 
     155      REAL(wp) ::   zhura, zhvra          !   -      - 
     156      REAL(wp) ::   za0, za1, za2, za3    !   -      - 
     157      ! 
     158      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     159      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 
     160      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 
     161      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 
     162      REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b 
     163      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
     164      REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf 
    128165      !!---------------------------------------------------------------------- 
    129166      ! 
    130167      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
    131168      ! 
    132       CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
    133       CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
    134       CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    135       ! 
    136       IF( kt == nit000 ) THEN             !* initialisation 
     169      !                                         !* Allocate temporay arrays 
     170      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 
     171      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  ) 
     172      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 
     173      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
     174      CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b                                     ) 
     175      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
     176      CALL wrk_alloc( jpi, jpj, zht, zhf ) 
     177      ! 
     178      !                                         !* Local constant initialization 
     179      z1_12 = 1._wp / 12._wp  
     180      z1_8  = 0.125_wp                                    
     181      z1_4  = 0.25_wp 
     182      z1_2  = 0.5_wp      
     183      zraur = 1._wp / rau0 
     184      ! 
     185      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
     186        z2dt_bf = rdt 
     187      ELSE 
     188        z2dt_bf = 2.0_wp * rdt 
     189      ENDIF 
     190      z1_2dt_b = 1.0_wp / z2dt_bf  
     191      ! 
     192      ll_init = ln_bt_av                           ! if no time averaging, then no specific restart  
     193      ll_fw_start = .FALSE. 
     194      ! 
     195                                                       ! time offset in steps for bdy data update 
     196      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
     197      ! 
     198      IF( kt == nit000 ) THEN                !* initialisation 
    137199         ! 
    138200         IF(lwp) WRITE(numout,*) 
    … …  
    141203         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro 
    142204         ! 
    143          CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b   
    144          ! 
    145          ua_e  (:,:) = un_b (:,:) 
    146          va_e  (:,:) = vn_b (:,:) 
    147          hu_e  (:,:) = hu   (:,:) 
    148          hv_e  (:,:) = hv   (:,:) 
    149          hur_e (:,:) = hur  (:,:) 
    150          hvr_e (:,:) = hvr  (:,:) 
    151          IF( ln_dynvor_een ) THEN 
    152             ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
     205         IF (neuler==0) ll_init=.TRUE. 
     206         ! 
     207         IF (ln_bt_fw.OR.(neuler==0)) THEN 
     208           ll_fw_start=.TRUE. 
     209           noffset = 0 
     210         ELSE 
     211           ll_fw_start=.FALSE. 
     212         ENDIF 
     213         ! 
     214         ! Set averaging weights and cycle length: 
     215         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 
     216         ! 
     217         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' )  
     218         ! 
     219      ENDIF 
     220      ! 
     221      ! Set arrays to remove/compute coriolis trend. 
     222      ! Do it once at kt=nit000 if volume is fixed, else at each long time step. 
     223      ! Note that these arrays are also used during barotropic loop. These are however frozen 
     224      ! although they should be updated in variable volume case. Not a big approximation. 
     225      ! To remove this approximation, copy lines below inside barotropic loop 
     226      ! and update depths at T-F points (ht and hf resp.) at each barotropic time step 
     227      ! 
     228      IF ( kt == nit000 .OR. lk_vvl ) THEN 
     229         IF ( ln_dynvor_een ) THEN 
     230            ! JC: Simplification needed below: define ht_0 even when volume is fixed 
     231            IF (lk_vvl) THEN 
     232               zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1)  
     233            ELSE 
     234               zht(:,:) = 0. 
     235               DO jk = 1, jpkm1 
     236                  zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     237               END DO 
     238            ENDIF 
     239 
     240            DO jj = 1, jpjm1 
     241               DO ji = 1, jpim1 
     242                  zwz(ji,jj) =   ( zht(ji  ,jj+1) + zht(ji+1,jj+1) +                     & 
     243                        &          zht(ji  ,jj  ) + zht(ji+1,jj  )   )                   & 
     244                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
     245                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     246                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
     247               END DO 
     248            END DO 
     249            CALL lbc_lnk( zwz, 'F', 1._wp ) 
     250            zwz(:,:) = ff(:,:) * zwz(:,:) 
     251 
     252            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    153253            DO jj = 2, jpj 
    154254               DO ji = fs_2, jpi   ! vector opt. 
    155                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp 
    156                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp 
    157                   ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 
    158                   ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp 
    159                END DO 
    160             END DO 
    161          ENDIF 
    162          ! 
    163       ENDIF 
    164  
    165       !                                                     !* Local constant initialization 
    166       z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step 
    167       IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic  
    168                                                                         ! time step (euler timestep) 
    169       z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates 
    170       z1_4     = 0.25_wp         
    171       zraur    = 1._wp / rau0                               ! 1 / volumic mass 
    172       ! 
    173       zhdiv(:,:) = 0._wp                                    ! barotropic divergence 
    174       zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
    175       zv_sld = 0._wp   ;   zv_asp = 0._wp 
    176  
    177       IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction 
    178         z2dt_bf = rdt 
    179       ELSE 
    180         z2dt_bf = 2.0_wp * rdt 
    181       ENDIF 
    182  
     255                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     256                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     257                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
     258                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
     259               END DO 
     260            END DO 
     261         ELSE 
     262            zwz(:,:) = 0._wp 
     263            zht(:,:) = 0. 
     264            IF ( .not. ln_sco ) THEN 
     265!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     266!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     267!              ENDIF 
     268!              zht(:,:) = gdepw_0(:,:,jk+1) 
     269            ELSE 
     270               zht(:,:) = hbatf(:,:) 
     271            END IF 
     272 
     273            DO jj = 1, jpjm1 
     274               zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     275            END DO 
     276 
     277            DO jk = 1, jpkm1 
     278               DO jj = 1, jpjm1 
     279                  zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     280               END DO 
     281            END DO 
     282            CALL lbc_lnk( zht, 'F', 1._wp ) 
     283            ! JC: TBC. hf should be greater than 0  
     284            DO jj = 1, jpj 
     285               DO ji = 1, jpi 
     286                  IF( zht(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array 
     287               END DO 
     288            END DO 
     289            zwz(:,:) = ff(:,:) * zwz(:,:) 
     290         ENDIF 
     291      ENDIF 
     292      ! 
     293      ! If forward start at previous time step, and centered integration,  
     294      ! then update averaging weights: 
     295      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN 
     296         ll_fw_start=.FALSE. 
     297         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 
     298      ENDIF 
     299 
     300      ! before inverse water column height at u- and v- points 
     301      IF( lk_vvl ) THEN 
     302         zhur_b(:,:) = 0. 
     303         zhvr_b(:,:) = 0. 
     304         DO jk = 1, jpk 
     305            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 
     306            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
     307         END DO 
     308         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) ) 
     309         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) ) 
     310      ELSE 
     311         zhur_b(:,:) = hur(:,:) 
     312         zhvr_b(:,:) = hvr(:,:) 
     313      ENDIF 
     314                           
    183315      ! ----------------------------------------------------------------------------- 
    184316      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    185317      ! ----------------------------------------------------------------------------- 
    186318      !       
     319      ! Some vertical sums (at now and before time steps) below could be suppressed  
     320      ! if one swap barotropic arrays somewhere 
     321      ! 
    187322      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
    188       !                                   ! -------------------------- 
    189       zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp 
    190       zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
     323      !                                   ! -------------------------------------------------- 
     324      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp 
     325      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp 
    191326      ! 
    192327      DO jk = 1, jpkm1 
    … …  
    198333            DO ji = 1, jpi 
    199334#endif 
    200                !                                                                              ! now trend 
    201                zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    202                zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    203                !                                                                              ! now velocity  
    204                zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk) 
    205                zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)                
    206                ! 
    207 #if defined key_vvl 
    208                ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk)  
    209                vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk) 
    210 #else 
    211                ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    212                vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
    213 #endif 
     335               !        ! now trend:                                                                    
     336               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
     337               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)          
     338               !        ! now bt transp:                    
     339               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)         
     340               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     341               !  ! before bt transp: 
     342               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
     343               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
    214344            END DO 
    215345         END DO 
    216346      END DO 
    217  
     347      ! 
     348      zu_frc(:,:) = zu_frc(:,:) * hur(:,:) 
     349      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:) 
     350      ! 
     351      IF( lk_vvl ) THEN 
     352          ub_b(:,:) = ub_b(:,:) * zhur_b(:,:) 
     353          vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:) 
     354      ELSE 
     355          ub_b(:,:) = ub_b(:,:) * hur(:,:) 
     356          vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
     357      ENDIF 
     358      ! 
    218359      !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    219       DO jk = 1, jpkm1                    ! -------------------------- 
     360      DO jk = 1, jpkm1                    ! ----------------------------------------------------------- 
    220361         DO jj = 2, jpjm1 
    221362            DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
    223                va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
     363               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 
     364               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 
    224365            END DO 
    225366         END DO 
    226367      END DO 
    227  
    228       !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent) 
    229       !                                   ! ---------------------------==== 
    230       zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport  
    231       zwy(:,:) = zvn(:,:) * e1v(:,:) 
     368      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
     369      !                                   ! -------------------------------------------------------- 
     370      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport  
     371      zwy(:,:) = vn_b(:,:) * e1v(:,:) 
    232372      ! 
    233373      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
    … …  
    239379               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    240380               ! energy conserving formulation for planetary vorticity term 
    241                zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
    242                zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    243             END DO 
    244          END DO 
    245          ! 
    246       ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme 
     381               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     382               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
     383            END DO 
     384         END DO 
     385         ! 
     386      ELSEIF ( ln_dynvor_ens ) THEN             ! enstrophy conserving scheme 
    247387         DO jj = 2, jpjm1 
    248388            DO ji = fs_2, fs_jpim1   ! vector opt. 
    249                zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    250                zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    251                zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    252                zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    253             END DO 
    254          END DO 
    255          ! 
    256       ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme 
     389               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
     390                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     391               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
     392                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     393               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     394               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     395            END DO 
     396         END DO 
     397         ! 
     398      ELSEIF ( ln_dynvor_een ) THEN             ! enstrophy and energy conserving scheme 
    257399         DO jj = 2, jpjm1 
    258400            DO ji = fs_2, fs_jpim1   ! vector opt. 
    259                zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    260                   &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    261                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)   & 
    262                   &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    263             END DO 
    264          END DO 
    265          ! 
    266       ENDIF 
    267  
     401               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
     402                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
     403                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
     404                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     405               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
     406                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
     407                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
     408                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     409            END DO 
     410         END DO 
     411         ! 
     412      ENDIF  
     413      ! 
     414      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities 
     415      vn_b (:,:) = vn_b(:,:) * hvr(:,:)   
    268416      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    269417      !                                   ! ---------------------------------------------------- 
    270       IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient 
     418      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    271419         DO jj = 2, jpjm1  
    272420            DO ji = fs_2, fs_jpim1   ! vector opt. 
    273                zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    & 
    274                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
    275                zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    & 
    276                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
    277             END DO 
    278          END DO 
    279       ENDIF 
    280  
    281       DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
     421               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     422               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     423            END DO 
     424         END DO 
     425      ENDIF 
     426 
     427      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    282428         DO ji = fs_2, fs_jpim1 
    283              zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    284              zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
     429             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 
     430             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
    285431          END DO 
    286       END DO 
    287  
    288                      
    289       !                                             ! Remove barotropic contribution of bottom friction  
    290       !                                             ! from the barotropic transport trend 
    291       zcoef = -1._wp * z1_2dt_b 
    292  
    293       IF(ln_bfrimp) THEN 
    294       !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
    295       !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction  
    296         DO jj = 2, jpjm1          
    297            DO ji = fs_2, fs_jpim1 
    298               ikbu = mbku(ji,jj) 
    299               ikbv = mbkv(ji,jj) 
    300               ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 
    301               va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 
    302  
    303               zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
    304               zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
    305            END DO 
    306         END DO 
    307  
    308       ELSE 
    309  
    310 # if defined key_vectopt_loop 
    311         DO jj = 1, 1 
    312            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    313 # else 
    314         DO jj = 2, jpjm1 
    315            DO ji = 2, jpim1 
    316 # endif 
    317             ! Apply stability criteria for bottom friction 
    318             !RBbug for vvl and external mode we may need to use varying fse3 
    319             !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    320               zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    321               zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    322            END DO 
    323         END DO 
    324  
    325         IF( lk_vvl ) THEN 
    326            DO jj = 2, jpjm1 
    327               DO ji = fs_2, fs_jpim1   ! vector opt. 
    328                  zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
    329                     &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 
    330                  zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
    331                     &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 
    332               END DO 
    333            END DO 
    334         ELSE 
    335            DO jj = 2, jpjm1 
    336               DO ji = fs_2, fs_jpim1   ! vector opt. 
    337                  zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
    338                  zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
    339               END DO 
    340            END DO 
    341         ENDIF 
    342       END IF    ! end (ln_bfrimp) 
    343  
    344                      
    345       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    346       !                                   ! --------------------------  
    347       zua(:,:) = zua(:,:) * hur(:,:) 
    348       zva(:,:) = zva(:,:) * hvr(:,:) 
    349       ! 
    350       IF( lk_vvl ) THEN 
    351          ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
    352          vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    353       ELSE 
    354          ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    355          vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    356       ENDIF 
    357  
     432      END DO  
     433      ! 
     434      !                 ! Add bottom stress contribution from baroclinic velocities:       
     435      IF (ln_bt_fw) THEN 
     436         DO jj = 2, jpjm1                           
     437            DO ji = fs_2, fs_jpim1   ! vector opt. 
     438               ikbu = mbku(ji,jj)        
     439               ikbv = mbkv(ji,jj)     
     440               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 
     441               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
     442            END DO 
     443         END DO 
     444      ELSE 
     445         DO jj = 2, jpjm1 
     446            DO ji = fs_2, fs_jpim1   ! vector opt. 
     447               ikbu = mbku(ji,jj)        
     448               ikbv = mbkv(ji,jj)     
     449               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 
     450               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
     451            END DO 
     452         END DO 
     453      ENDIF 
     454      ! 
     455      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
     456      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
     457      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     458      !        
     459      IF (ln_bt_fw) THEN                        ! Add wind forcing 
     460         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 
     461         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:) 
     462      ELSE 
     463         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 
     464         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 
     465      ENDIF   
     466      ! 
     467      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing 
     468         IF (ln_bt_fw) THEN 
     469            DO jj = 2, jpjm1               
     470               DO ji = fs_2, fs_jpim1   ! vector opt. 
     471                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj) 
     472                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj) 
     473                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
     474                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
     475               END DO 
     476            END DO 
     477         ELSE 
     478            DO jj = 2, jpjm1               
     479               DO ji = fs_2, fs_jpim1   ! vector opt. 
     480                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    & 
     481                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     482                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    & 
     483                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     484                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
     485                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
     486               END DO 
     487            END DO 
     488         ENDIF  
     489      ENDIF 
     490      !                                   !* Right-Hand-Side of the barotropic ssh equation 
     491      !                                   ! ----------------------------------------------- 
     492      !                                         ! Surface net water flux and rivers 
     493      IF (ln_bt_fw) THEN 
     494         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 
     495      ELSE 
     496         zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 
     497      ENDIF 
     498#if defined key_asminc 
     499      !                                         ! Include the IAU weighted SSH increment 
     500      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     501         zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:) 
     502      ENDIF 
     503#endif 
     504      ! 
    358505      ! ----------------------------------------------------------------------- 
    359       !  Phase 2 : Integration of the barotropic equations with time splitting 
     506      !  Phase 2 : Integration of the barotropic equations  
    360507      ! ----------------------------------------------------------------------- 
    361508      ! 
    362509      !                                             ! ==================== ! 
    363510      !                                             !    Initialisations   ! 
     511      !                                             ! ==================== !   
     512      ! Initialize barotropic variables:     
     513      IF (ln_bt_fw) THEN                  ! FORWARD integration:  start from NOW fields                              
     514         sshn_e (:,:) = sshn (:,:)             
     515         zun_e  (:,:) = un_b (:,:)             
     516         zvn_e  (:,:) = vn_b (:,:) 
     517      ELSE                                ! CENTERED integration: start from BEFORE fields 
     518         sshn_e (:,:) = sshb (:,:) 
     519         zun_e  (:,:) = ub_b (:,:)          
     520         zvn_e  (:,:) = vb_b (:,:) 
     521      ENDIF 
     522      ! 
     523      ! Initialize depths: 
     524      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 
     525         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 
     526         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 
     527         hur_e (:,:) = zhur_b(:,:) 
     528         hvr_e (:,:) = zhvr_b(:,:) 
     529      ELSE 
     530         hu_e  (:,:) = hu   (:,:)        
     531         hv_e  (:,:) = hv   (:,:)  
     532         hur_e (:,:) = hur  (:,:)     
     533         hvr_e (:,:) = hvr  (:,:) 
     534      ENDIF 
     535      ! 
     536      IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5: 
     537         zhup2_e (:,:) = hu(:,:) 
     538         zhvp2_e (:,:) = hv(:,:) 
     539      ENDIF 
     540      ! 
     541      ! Initialize sums: 
     542      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     543      va_b  (:,:) = 0._wp 
     544      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level 
     545      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
     546      zv_sum(:,:) = 0._wp 
    364547      !                                             ! ==================== ! 
    365       icycle = 2  * nn_baro            ! Number of barotropic sub time-step 
    366        
    367       !                                ! Start from NOW field 
    368       hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points 
    369       hv_e   (:,:) = hv   (:,:)  
    370       hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points 
    371       hvr_e  (:,:) = hvr  (:,:) 
    372 !RBbug     zsshb_e(:,:) = sshn (:,:)   
    373       zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now) 
    374       sshn_e (:,:) = sshn (:,:) 
    375        
    376       zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external) 
    377       zvn_e  (:,:) = vn_b (:,:) 
    378       zub_e  (:,:) = un_b (:,:) 
    379       zvb_e  (:,:) = vn_b (:,:) 
    380  
    381       zu_sum  (:,:) = un_b (:,:)           ! summation 
    382       zv_sum  (:,:) = vn_b (:,:) 
    383       zssh_sum(:,:) = sshn (:,:) 
    384  
    385 #if defined key_obc 
    386       ! set ssh corrections to 0 
    387       ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 
    388       IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp 
    389       IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp 
    390       IF( lp_obc_south )   sshfos_b(:,:) = 0._wp 
    391       IF( lp_obc_north )   sshfon_b(:,:) = 0._wp 
    392 #endif 
    393  
    394       !                                             ! ==================== ! 
    395       DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1) 
     548      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    396549         !                                          ! ==================== ! 
    397          z2dt_e = 2. * ( rdt / nn_baro ) 
    398          IF( jn == 1 )   z2dt_e = rdt / nn_baro 
    399  
    400550         !                                                !* Update the forcing (BDY and tides) 
    401551         !                                                !  ------------------ 
    402          IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    403          IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
    404          IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 
    405  
    406          !                                                !* after ssh_e 
     552         ! Update only tidal forcing at open boundaries 
     553#if defined key_tide 
     554         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
     555         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset ) 
     556#endif 
     557         ! 
     558         ! Set extrapolation coefficients for predictor step: 
     559         IF ((jn<3).AND.ll_init) THEN      ! Forward            
     560           za1 = 1._wp                                           
     561           za2 = 0._wp                         
     562           za3 = 0._wp                         
     563         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105  
     564           za1 =  1.781105_wp              ! za1 =   3/2 +   bet 
     565           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet) 
     566           za3 =  0.281105_wp              ! za3 = bet 
     567         ENDIF 
     568 
     569         ! Extrapolate barotropic velocities at step jit+0.5: 
     570         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
     571         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
     572 
     573         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
     574            !                                             !  ------------------ 
     575            ! Extrapolate Sea Level at step jit+0.5: 
     576            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
     577            ! 
     578            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
     579               DO ji = 2, fs_jpim1   ! Vector opt. 
     580                  zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     581                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     582                     &              +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 
     583                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     584                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     585                     &              +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
     586               END DO 
     587            END DO 
     588            CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp ) 
     589            ! 
     590            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)               ! Ocean depth at U- and V-points 
     591            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     592         ENDIF 
     593         !                                                !* after ssh 
    407594         !                                                !  ----------- 
    408          DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
     595         ! One should enforce volume conservation at open boundaries here 
     596         ! considering fluxes below: 
     597         ! 
     598         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
     599         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
     600         DO jj = 2, jpjm1                                  
    409601            DO ji = fs_2, fs_jpim1   ! vector opt. 
    410                zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
    411                   &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     & 
    412                   &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     & 
    413                   &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
    414             END DO 
    415          END DO 
    416          ! 
    417 #if defined key_obc 
    418          !                                                     ! OBC : zhdiv must be zero behind the open boundary 
    419 !!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    420          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east 
    421          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west 
    422          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north 
    423          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south 
     602               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   & 
     603                  &             + zwy(ji,jj) - zwy(ji,jj-1)   & 
     604                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     605            END DO 
     606         END DO 
     607         ! 
     608         ! Sum over sub-time-steps to compute advective velocities 
     609         za2 = wgtbtp2(jn) 
     610         zu_sum  (:,:) = zu_sum  (:,:) + za2 * ua_e  (:,:) * zhup2_e  (:,:) 
     611         zv_sum  (:,:) = zv_sum  (:,:) + za2 * va_e  (:,:) * zhvp2_e  (:,:) 
     612         ! 
     613         ! Set next sea level: 
     614         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     615         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
     616 
     617#if defined key_bdy 
     618         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 
     619         IF (lk_bdy) CALL bdy_ssh( ssha_e ) 
    424620#endif 
    425 #if defined key_bdy 
    426          zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask 
     621#if defined key_agrif 
     622         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 
    427623#endif 
    428          ! 
    429          DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
    430             DO ji = fs_2, fs_jpim1   ! vector opt. 
    431                ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)  
    432             END DO 
    433          END DO 
    434  
    435          !                                                !* after barotropic velocities (vorticity scheme dependent) 
    436          !                                                !  ---------------------------   
    437          zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport 
    438          zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
     624         !   
     625         ! Sea Surface Height at u-,v-points (vvl case only) 
     626         IF ( lk_vvl ) THEN                                 
     627            DO jj = 2, jpjm1 
     628               DO ji = 2, jpim1      ! NO Vector Opt. 
     629                  zsshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 & 
     630                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) & 
     631                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) ) 
     632                  zsshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 & 
     633                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) & 
     634                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) ) 
     635               END DO 
     636            END DO 
     637            CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) 
     638         ENDIF    
     639         !                                  
     640         ! Half-step back interpolation of SSH for surface pressure computation: 
     641         !---------------------------------------------------------------------- 
     642         IF ((jn==1).AND.ll_init) THEN 
     643           za0=1._wp                        ! Forward-backward 
     644           za1=0._wp                            
     645           za2=0._wp 
     646           za3=0._wp 
     647         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 
     648           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps 
     649           za1=-0.1666666666666_wp          ! za1 = gam 
     650           za2= 0.0833333333333_wp          ! za2 = eps 
     651           za3= 0._wp               
     652         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
     653           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps     
     654           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps  
     655           za2=0.088_wp                     ! za2 = gam 
     656           za3=0.013_wp                     ! za3 = eps 
     657         ENDIF 
     658 
     659         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
     660          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     661 
     662         ! 
     663         ! Compute associated depths at U and V points: 
     664         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN       
     665            !                                         
     666            DO jj = 2, jpjm1                             
     667               DO ji = 2, jpim1 
     668                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     669                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     670                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                  
     671                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     672                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     673                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
     674                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1  
     675                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 
     676               END DO 
     677            END DO 
     678         ENDIF 
     679         ! 
     680         ! Add Coriolis trend: 
     681         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated 
     682         ! at each time step. We however keep them constant here for optimization. 
     683         ! Recall that zwx and zwy arrays hold fluxes at this stage: 
     684         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5 
     685         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    439686         ! 
    440687         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==! 
    441688            DO jj = 2, jpjm1 
    442689               DO ji = fs_2, fs_jpim1   ! vector opt. 
    443                   ! surface pressure gradient 
    444                   IF( lk_vvl) THEN 
    445                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    446                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    447                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    448                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    449                   ELSE 
    450                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    451                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    452                   ENDIF 
    453                   ! add tidal astronomical forcing 
    454                   IF ( ln_tide_pot .AND. lk_tide ) THEN  
    455                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    456                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    457                   ENDIF 
    458                   ! energy conserving formulation for planetary vorticity term 
    459690                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
    460691                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    461692                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    462693                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    463                   zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 
    464                   zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    465                   ! after velocities with implicit bottom friction 
    466  
    467                   IF( ln_bfrimp ) THEN      ! implicit bottom friction 
    468                      !   A new method to implement the implicit bottom friction.  
    469                      !   H. Liu 
    470                      !   Sept 2011 
    471                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    472                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    473                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    474                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    475                      ! 
    476                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    477                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    478                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    479                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    480                      ! 
    481                   ELSE 
    482                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    483                       &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    484                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    485                       &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    486                   ENDIF 
     694                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     695                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    487696               END DO 
    488697            END DO 
    … …  
    491700            DO jj = 2, jpjm1 
    492701               DO ji = fs_2, fs_jpim1   ! vector opt. 
    493                    ! surface pressure gradient 
    494                   IF( lk_vvl) THEN 
    495                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    496                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    497                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    498                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    499                   ELSE 
    500                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    501                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    502                   ENDIF 
    503                   ! add tidal astronomical forcing 
    504                   IF ( ln_tide_pot .AND. lk_tide ) THEN 
    505                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    506                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    507                   ENDIF 
    508                   ! enstrophy conserving formulation for planetary vorticity term 
    509                   zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    510                   zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    511                   zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 
    512                   zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    513                   ! after velocities with implicit bottom friction 
    514                   IF( ln_bfrimp ) THEN 
    515                      !   A new method to implement the implicit bottom friction.  
    516                      !   H. Liu 
    517                      !   Sept 2011 
    518                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    519                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    520                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    521                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    522                      ! 
    523                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    524                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    525                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    526                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    527                      ! 
    528                   ELSE 
    529                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    530                      &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    531                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    532                      &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    533                   ENDIF 
     702                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
     703                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     704                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
     705                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     706                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     707                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    534708               END DO 
    535709            END DO 
    … …  
    538712            DO jj = 2, jpjm1 
    539713               DO ji = fs_2, fs_jpim1   ! vector opt. 
    540                   ! surface pressure gradient 
    541                   IF( lk_vvl) THEN 
    542                      zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    543                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    544                      zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    545                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    546                   ELSE 
    547                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    548                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    549                   ENDIF 
    550                   ! add tidal astronomical forcing 
    551                   IF ( ln_tide_pot .AND. lk_tide ) THEN 
    552                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    553                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    554                   ENDIF 
    555                   ! energy/enstrophy conserving formulation for planetary vorticity term 
    556                   zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    557                      &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 
    558                   zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    559                      &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    560                   ! after velocities with implicit bottom friction 
    561                   IF( ln_bfrimp ) THEN 
    562                      !   A new method to implement the implicit bottom friction.  
    563                      !   H. Liu 
    564                      !   Sept 2011 
    565                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    566                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    567                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    568                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    569                      ! 
    570                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    571                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    572                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    573                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    574                      ! 
    575                   ELSE 
    576                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    577                      &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    578                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    579                      &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    580                   ENDIF 
     714                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
     715                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
     716                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) &  
     717                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     718                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &  
     719                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
     720                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &  
     721                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    581722               END DO 
    582723            END DO 
    583724            !  
    584725         ENDIF 
    585          !                                                     !* domain lateral boundary 
    586          !                                                     !  ----------------------- 
    587  
    588                                                                ! OBC open boundaries 
    589          IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 
    590  
    591                                                                ! BDY open boundaries 
    592 #if defined key_bdy 
    593          pssh => sshn_e 
    594          phur => hur_e 
    595          phvr => hvr_e 
    596          pu2d => ua_e 
    597          pv2d => va_e 
    598  
    599          IF( lk_bdy )   CALL bdy_dyn2d( kt )  
    600 #endif 
    601  
    602          ! 
    603          CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
    604          CALL lbc_lnk( va_e  , 'V', -1. ) 
    605          CALL lbc_lnk( ssha_e, 'T',  1. ) 
    606  
    607          zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps 
    608          zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:)  
    609          zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)  
    610  
    611          !                                                !* Time filter and swap 
    612          !                                                !  -------------------- 
    613          IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step) 
    614             zsshb_e(:,:) = sshn_e(:,:) 
    615             zub_e  (:,:) = zun_e (:,:) 
    616             zvb_e  (:,:) = zvn_e (:,:) 
    617             sshn_e (:,:) = ssha_e(:,:) 
    618             zun_e  (:,:) = ua_e  (:,:) 
    619             zvn_e  (:,:) = va_e  (:,:) 
    620          ELSE                                                   ! Swap + Filter 
    621             zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
    622             zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:) 
    623             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:) 
    624             sshn_e (:,:) = ssha_e(:,:) 
    625             zun_e  (:,:) = ua_e  (:,:) 
    626             zvn_e  (:,:) = va_e  (:,:) 
    627          ENDIF 
    628  
    629          IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
    630             !                                             !  ------------------ 
    631             DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points 
    632                DO ji = 1, fs_jpim1   ! Vector opt. 
    633                   zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
    634                      &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
    635                      &                     +   e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
    636                   zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
    637                      &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
    638                      &                     +   e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
    639                END DO 
    640             END DO 
    641             CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions 
    642             CALL lbc_lnk( zsshvn_e, 'V', 1. )  
    643             ! 
    644             hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points 
    645             hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 
     726         ! 
     727         ! Add tidal astronomical forcing if defined 
     728         IF ( lk_tide.AND.ln_tide_pot ) THEN 
     729            DO jj = 2, jpjm1 
     730               DO ji = fs_2, fs_jpim1   ! vector opt. 
     731                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     732                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     733                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 
     734                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 
     735               END DO 
     736            END DO 
     737         ENDIF 
     738         ! 
     739         ! Add bottom stresses: 
     740         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
     741         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
     742         ! 
     743         ! Surface pressure trend: 
     744         DO jj = 2, jpjm1 
     745            DO ji = fs_2, fs_jpim1   ! vector opt. 
     746               ! Add surface pressure gradient 
     747               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     748               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     749               zwx(ji,jj) = zu_spg 
     750               zwy(ji,jj) = zv_spg 
     751            END DO 
     752         END DO 
     753         ! 
     754         ! Set next velocities: 
     755         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form 
     756            DO jj = 2, jpjm1 
     757               DO ji = fs_2, fs_jpim1   ! vector opt. 
     758                  ua_e(ji,jj) = (                                zun_e(ji,jj)   &  
     759                            &     + rdtbt * (                      zwx(ji,jj)   & 
     760                            &                                 + zu_trd(ji,jj)   & 
     761                            &                                 + zu_frc(ji,jj) ) &  
     762                            &   ) * umask(ji,jj,1) 
     763 
     764                  va_e(ji,jj) = (                                zvn_e(ji,jj)   & 
     765                            &     + rdtbt * (                      zwy(ji,jj)   & 
     766                            &                                 + zv_trd(ji,jj)   & 
     767                            &                                 + zv_frc(ji,jj) ) & 
     768                            &   ) * vmask(ji,jj,1) 
     769               END DO 
     770            END DO 
     771 
     772         ELSE                 ! Flux form 
     773            DO jj = 2, jpjm1 
     774               DO ji = fs_2, fs_jpim1   ! vector opt. 
     775 
     776                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
     777                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
     778 
     779                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
     780                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     781                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
     782                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     783                            &   ) * zhura 
     784 
     785                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   & 
     786                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   & 
     787                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
     788                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) & 
     789                            &   ) * zhvra 
     790               END DO 
     791            END DO 
     792         ENDIF 
     793         ! 
     794         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
     795            !                                          !  ----------------------------------------------         
     796            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     797            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    646798            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    647799            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
    648800            ! 
    649801         ENDIF 
     802         !                                                 !* domain lateral boundary 
     803         !                                                 !  ----------------------- 
     804         ! 
     805         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries  
     806         CALL lbc_lnk( va_e  , 'V', -1._wp ) 
     807 
     808#if defined key_bdy   
     809 
     810         pssh => ssha_e 
     811         phur => hur_e 
     812         phvr => hvr_e 
     813         pua2d => ua_e 
     814         pva2d => va_e 
     815         pub2d => zun_e 
     816         pvb2d => zvn_e 
     817                                       
     818         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries 
     819#endif 
     820#if defined key_agrif 
     821         IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( kt, jn ) ! Agrif 
     822#endif 
     823         !                                             !* Swap 
     824         !                                             !  ---- 
     825         ubb_e  (:,:) = ub_e  (:,:) 
     826         ub_e   (:,:) = zun_e (:,:) 
     827         zun_e  (:,:) = ua_e  (:,:) 
     828         ! 
     829         vbb_e  (:,:) = vb_e  (:,:) 
     830         vb_e   (:,:) = zvn_e (:,:) 
     831         zvn_e  (:,:) = va_e  (:,:) 
     832         ! 
     833         sshbb_e(:,:) = sshb_e(:,:) 
     834         sshb_e (:,:) = sshn_e(:,:) 
     835         sshn_e (:,:) = ssha_e(:,:) 
     836 
     837         !                                             !* Sum over whole bt loop 
     838         !                                             !  ---------------------- 
     839         za1 = wgtbtp1(jn)                                     
     840         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities 
     841            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
     842            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
     843         ELSE                                ! Sum transports 
     844            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
     845            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     846         ENDIF 
     847         !                                   ! Sum sea level 
     848         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
    650849         !                                                 ! ==================== ! 
    651850      END DO                                               !        end loop      ! 
    652851      !                                                    ! ==================== ! 
    653  
    654 #if defined key_obc 
    655       IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ????? 
    656       IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:) 
    657       IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:) 
    658       IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:) 
    659 #endif 
    660  
    661852      ! ----------------------------------------------------------------------------- 
    662853      ! Phase 3. update the general trend with the barotropic trend 
    663854      ! ----------------------------------------------------------------------------- 
    664855      ! 
    665       !                                   !* Time average ==> after barotropic u, v, ssh 
    666       zcoef =  1._wp / ( 2 * nn_baro  + 1 )  
    667       zu_sum(:,:) = zcoef * zu_sum  (:,:)  
    668       zv_sum(:,:) = zcoef * zv_sum  (:,:)  
    669       !  
    670       !                                   !* update the general momentum trend 
    671       DO jk=1,jpkm1 
    672          ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 
    673          va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 
     856      ! At this stage ssha holds a time averaged value 
     857      !                                                ! Sea Surface Height at u-,v- and f-points 
     858      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     859         DO jj = 1, jpjm1 
     860            DO ji = 1, jpim1      ! NO Vector Opt. 
     861               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
     862                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
     863                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     864               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
     865                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
     866                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     867            END DO 
     868         END DO 
     869         CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions 
     870      ENDIF 
     871      ! 
     872      ! Set advection velocity correction: 
     873      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
     874         un_adv(:,:) = zu_sum(:,:)*hur(:,:) 
     875         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 
     876      ELSE 
     877         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 
     878         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 
     879      END IF 
     880 
     881      IF (ln_bt_fw) THEN ! Save integrated transport for next computation 
     882         ub2_b(:,:) = zu_sum(:,:) 
     883         vb2_b(:,:) = zv_sum(:,:) 
     884      ENDIF 
     885      ! 
     886      ! Update barotropic trend: 
     887      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 
     888         DO jk=1,jpkm1 
     889            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
     890            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
     891         END DO 
     892      ELSE 
     893         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 
     894         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 
     895         DO jk=1,jpkm1 
     896            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b 
     897            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b 
     898         END DO 
     899         ! Save barotropic velocities not transport: 
     900         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 
     901         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     902      ENDIF 
     903      ! 
     904      DO jk = 1, jpkm1 
     905         ! Correct velocities: 
     906         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 
     907         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 
     908         ! 
    674909      END DO 
    675       un_b  (:,:) =  zu_sum(:,:)  
    676       vn_b  (:,:) =  zv_sum(:,:)  
    677       sshn_b(:,:) = zcoef * zssh_sum(:,:)  
    678910      ! 
    679911      !                                   !* write time-spliting arrays in the restart 
    680       IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    681       ! 
    682       CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
    683       CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
    684       CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
     912      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' ) 
     913      ! 
     914      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 
     915      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 
     916      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 
     917      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
     918      CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b                                     ) 
     919      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
     920      CALL wrk_dealloc( jpi, jpj, zht, zhf ) 
    685921      ! 
    686922      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    … …  
    688924   END SUBROUTINE dyn_spg_ts 
    689925 
     926   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) 
     927      !!--------------------------------------------------------------------- 
     928      !!                   ***  ROUTINE ts_wgt  *** 
     929      !! 
     930      !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 
     931      !!---------------------------------------------------------------------- 
     932      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true. 
     933      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true. 
     934      INTEGER, INTENT(inout) :: jpit      ! cycle length     
     935      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights 
     936                                                         zwgt2    ! Secondary weights 
     937       
     938      INTEGER ::  jic, jn, ji                      ! temporary integers 
     939      REAL(wp) :: za1, za2 
     940      !!---------------------------------------------------------------------- 
     941 
     942      zwgt1(:) = 0._wp 
     943      zwgt2(:) = 0._wp 
     944 
     945      ! Set time index when averaged value is requested 
     946      IF (ll_fw) THEN  
     947         jic = nn_baro 
     948      ELSE 
     949         jic = 2 * nn_baro 
     950      ENDIF 
     951 
     952      ! Set primary weights: 
     953      IF (ll_av) THEN 
     954           ! Define simple boxcar window for primary weights  
     955           ! (width = nn_baro, centered around jic)      
     956         SELECT CASE ( nn_bt_flt ) 
     957              CASE( 0 )  ! No averaging 
     958                 zwgt1(jic) = 1._wp 
     959                 jpit = jic 
     960 
     961              CASE( 1 )  ! Boxcar, width = nn_baro 
     962                 DO jn = 1, 3*nn_baro 
     963                    za1 = ABS(float(jn-jic))/float(nn_baro)  
     964                    IF (za1 < 0.5_wp) THEN 
     965                      zwgt1(jn) = 1._wp 
     966                      jpit = jn 
     967                    ENDIF 
     968                 ENDDO 
     969 
     970              CASE( 2 )  ! Boxcar, width = 2 * nn_baro 
     971                 DO jn = 1, 3*nn_baro 
     972                    za1 = ABS(float(jn-jic))/float(nn_baro)  
     973                    IF (za1 < 1._wp) THEN 
     974                      zwgt1(jn) = 1._wp 
     975                      jpit = jn 
     976                    ENDIF 
     977                 ENDDO 
     978              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 
     979         END SELECT 
     980 
     981      ELSE ! No time averaging 
     982         zwgt1(jic) = 1._wp 
     983         jpit = jic 
     984      ENDIF 
     985     
     986      ! Set secondary weights 
     987      DO jn = 1, jpit 
     988        DO ji = jn, jpit 
     989             zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 
     990        END DO 
     991      END DO 
     992 
     993      ! Normalize weigths: 
     994      za1 = 1._wp / SUM(zwgt1(1:jpit)) 
     995      za2 = 1._wp / SUM(zwgt2(1:jpit)) 
     996      DO jn = 1, jpit 
     997        zwgt1(jn) = zwgt1(jn) * za1 
     998        zwgt2(jn) = zwgt2(jn) * za2 
     999      END DO 
     1000      ! 
     1001   END SUBROUTINE ts_wgt 
    6901002 
    6911003   SUBROUTINE ts_rst( kt, cdrw ) 
    … …  
    6981010      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    6991011      ! 
    700       INTEGER ::  ji, jk        ! dummy loop indices 
    7011012      !!---------------------------------------------------------------------- 
    7021013      ! 
    7031014      IF( TRIM(cdrw) == 'READ' ) THEN 
    704          IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 
    705             CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued 
    706             CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
     1015         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )    
     1016         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) )  
     1017         IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN 
     1018            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )    
     1019            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )    
     1020            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) ) 
     1021            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) )  
     1022            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )    
     1023            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) ) 
    7071024         ELSE 
    708             un_b (:,:) = 0._wp 
    709             vn_b (:,:) = 0._wp 
    710             ! vertical sum 
    711             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    712                DO jk = 1, jpkm1 
    713                   DO ji = 1, jpij 
    714                      un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
    715                      vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
    716                   END DO 
    717                END DO 
    718             ELSE                             ! No  vector opt. 
    719                DO jk = 1, jpkm1 
    720                   un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    721                   vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    722                END DO 
    723             ENDIF 
    724             un_b (:,:) = un_b(:,:) * hur(:,:) 
    725             vn_b (:,:) = vn_b(:,:) * hvr(:,:) 
    726          ENDIF 
    727  
    728          ! Vertically integrated velocity (before) 
    729          IF (neuler/=0) THEN 
    730             ub_b (:,:) = 0._wp 
    731             vb_b (:,:) = 0._wp 
    732  
    733             ! vertical sum 
    734             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    735                DO jk = 1, jpkm1 
    736                   DO ji = 1, jpij 
    737                      ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 
    738                      vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 
    739                   END DO 
    740                END DO 
    741             ELSE                             ! No  vector opt. 
    742                DO jk = 1, jpkm1 
    743                   ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 
    744                   vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 
    745                END DO 
    746             ENDIF 
    747  
    748             IF( lk_vvl ) THEN 
    749                ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
    750                vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    751             ELSE 
    752                ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    753                vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    754             ENDIF 
    755          ELSE                                 ! neuler==0 
    756             ub_b (:,:) = un_b (:,:) 
    757             vb_b (:,:) = vn_b (:,:) 
    758          ENDIF 
    759  
    760          IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 
    761             CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh 
    762          ELSE 
    763             sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value    
    764          ENDIF  
     1025            sshbb_e = sshn_b                                                ! ACC GUESS WORK 
     1026            ubb_e   = ub_b 
     1027            vbb_e   = vb_b 
     1028            sshb_e  = sshn_b 
     1029            ub_e    = ub_b 
     1030            vb_e    = vb_b 
     1031         ENDIF 
     1032      ! 
    7651033      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    766          CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh 
    767          CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop 
    768          CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !  
     1034         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) ) 
     1035         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) ) 
     1036         ! 
     1037         IF (.NOT.ln_bt_av) THEN 
     1038            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )  
     1039            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) ) 
     1040            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) ) 
     1041            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) ) 
     1042            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) ) 
     1043            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) ) 
     1044         ENDIF 
    7691045      ENDIF 
    7701046      ! 
    7711047   END SUBROUTINE ts_rst 
    7721048 
     1049   SUBROUTINE dyn_spg_ts_init( kt ) 
     1050      !!--------------------------------------------------------------------- 
     1051      !!                   ***  ROUTINE dyn_spg_ts_init  *** 
     1052      !! 
     1053      !! ** Purpose : Set time splitting options 
     1054      !!---------------------------------------------------------------------- 
     1055      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     1056      ! 
     1057      INTEGER  :: ji ,jj, jk 
     1058      REAL(wp) :: zxr2, zyr2, zcmax 
     1059      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht 
     1060      !! 
     1061!      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 
     1062!      &                  nn_baro, rn_bt_cmax, nn_bt_flt 
     1063      !!---------------------------------------------------------------------- 
     1064!      REWIND( numnam )              !* Namelist namsplit: split-explicit free surface 
     1065!      READ  ( numnam, namsplit ) 
     1066      !         ! Max courant number for ext. grav. waves 
     1067      ! 
     1068      CALL wrk_alloc( jpi, jpj, zcu, zht ) 
     1069      ! 
     1070      ! JC: Simplification needed below: define ht_0 even when volume is fixed 
     1071      IF (lk_vvl) THEN  
     1072         zht(:,:) = ht_0(:,:) * tmask(:,:,1) 
     1073      ELSE 
     1074         zht(:,:) = 0. 
     1075         DO jk = 1, jpkm1 
     1076            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     1077         END DO 
     1078      ENDIF 
     1079 
     1080      DO jj = 1, jpj 
     1081         DO ji =1, jpi 
     1082            zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 
     1083            zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 
     1084            zcu(ji,jj) = sqrt(grav*zht(ji,jj)*(zxr2 + zyr2) ) 
     1085         END DO 
     1086      END DO 
     1087 
     1088      zcmax = MAXVAL(zcu(:,:)) 
     1089      IF( lk_mpp )   CALL mpp_max( zcmax ) 
     1090 
     1091      ! Estimate number of iterations to satisfy a max courant number=0.8  
     1092      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
     1093       
     1094      rdtbt = rdt / FLOAT(nn_baro) 
     1095      zcmax = zcmax * rdtbt 
     1096                     ! Print results 
     1097      IF(lwp) WRITE(numout,*) 
     1098      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 
     1099      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     1100      IF( ln_bt_nn_auto ) THEN 
     1101         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro ' 
     1102         IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 
     1103      ELSE 
     1104         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist ' 
     1105      ENDIF 
     1106      IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro 
     1107      IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt 
     1108      IF(lwp) WRITE(numout,*) ' Maximum Courant number is   :', zcmax 
     1109 
     1110      IF(ln_bt_av) THEN 
     1111         IF(lwp) WRITE(numout,*) ' ln_bt_av=.true.  => Time averaging over nn_baro time steps is on ' 
     1112      ELSE 
     1113         IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables ' 
     1114      ENDIF 
     1115      ! 
     1116      ! 
     1117      IF(ln_bt_fw) THEN 
     1118         IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables ' 
     1119      ELSE 
     1120         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables ' 
     1121      ENDIF 
     1122      ! 
     1123      IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 
     1124      SELECT CASE ( nn_bt_flt ) 
     1125           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '      Dirac' 
     1126           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = nn_baro' 
     1127           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = 2*nn_baro'  
     1128           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 
     1129      END SELECT 
     1130      ! 
     1131      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN 
     1132         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 
     1133      ENDIF 
     1134      IF ( zcmax>0.9_wp ) THEN 
     1135         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )           
     1136      ENDIF 
     1137      ! 
     1138      CALL wrk_dealloc( jpi, jpj, zcu, zht ) 
     1139      ! 
     1140   END SUBROUTINE dyn_spg_ts_init 
     1141 
    7731142#else 
    774    !!---------------------------------------------------------------------- 
    775    !!   Default case :   Empty module   No standart free surface cst volume 
    776    !!---------------------------------------------------------------------- 
     1143   !!--------------------------------------------------------------------------- 
     1144   !!   Default case :   Empty module   No standard free surface constant volume 
     1145   !!--------------------------------------------------------------------------- 
     1146 
     1147   USE par_kind 
     1148   LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping 
    7771149CONTAINS 
    7781150   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function 
    … …  
    7871159      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    7881160      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw 
    789    END SUBROUTINE ts_rst     
     1161   END SUBROUTINE ts_rst   
     1162   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine 
     1163      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     1164      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt 
     1165   END SUBROUTINE dyn_spg_ts_init 
    7901166#endif 
    7911167    

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r4147 r4292  
    572572      INTEGER  ::   ierr                                          ! local integer 
    573573      REAL(wp) ::   zfac12, zua, zva                              ! local scalars 
     574      REAL(wp) ::   zmsk, ze3                                     ! local scalars 
    574575      !                                                           !  3D workspace  
    575576      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz 
    … …  
    577578#if defined key_vvl 
    578579      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T) 
    579 #endif 
    580 #if ! defined key_vvl 
     580#else 
    581581      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all 
    582582#endif 
    … …  
    604604      ENDIF 
    605605 
    606       IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t) 
     606      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 
    607607         DO jk = 1, jpk 
    608608            DO jj = 1, jpjm1 
    609609               DO ji = 1, jpim1 
    610                   ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    611                      &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25 
    612                   IF( ze3f(ji,jj,jk) /= 0._wp )   ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk) 
     610                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     611                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
     612                  zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
     613                     &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
     614                  IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
    613615               END DO 
    614616            END DO 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r3625 r4292  
    1616   USE oce             ! ocean dynamics and tracers 
    1717   USE dom_oce         ! ocean space and time domain 
     18   USE domvvl          ! variable volume 
    1819   USE sbc_oce         ! surface boundary condition: ocean 
    1920   USE zdf_oce         ! ocean vertical physics 
    … …  
    2425   USE wrk_nemo        ! Memory Allocation 
    2526   USE timing          ! Timing 
     27   USE dynadv          ! dynamics: vector invariant versus flux form 
     28   USE dynspg_oce, ONLY: lk_dynspg_ts, ua_b, va_b 
     29   USE dynspg_ts 
    2630 
    2731   IMPLICIT NONE 
    … …  
    2933 
    3034   PUBLIC   dyn_zdf_imp   ! called by step.F90 
     35 
     36   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise  
    3137 
    3238   !! * Substitutions 
    … …  
    6470      INTEGER  ::   ikbu, ikbv   ! local integers 
    6571      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
     72      REAL(wp) ::   ze3ua, ze3va 
    6673      !!---------------------------------------------------------------------- 
    6774 
    6875      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
    69       REAL(wp), POINTER, DIMENSION(:,:)   ::  zavmu, zavmv 
    7076      !!---------------------------------------------------------------------- 
    7177      ! 
    … …  
    7379      ! 
    7480      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )  
    75       CALL wrk_alloc( jpi,jpj, zavmu, zavmv )  
    7681      ! 
    7782      IF( kt == nit000 ) THEN 
    … …  
    7984         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
    8085         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     86         ! 
     87         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator 
     88         ELSE                ;    r_vvl = 0._wp        
     89         ENDIF 
    8190      ENDIF 
    8291 
    … …  
    94103      IF( ln_bfrimp ) THEN 
    95104# if defined key_vectopt_loop 
    96       DO jj = 1, 1 
    97          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     105         DO jj = 1, 1 
     106            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    98107# else 
    99       DO jj = 2, jpjm1 
    100          DO ji = 2, jpim1 
     108         DO jj = 2, jpjm1 
     109            DO ji = 2, jpim1 
    101110# endif 
    102             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    103             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    104             zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 
    105             zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 
    106             avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)  
    107             avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    108          END DO 
    109       END DO 
    110       ENDIF 
     111               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
     112               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     113               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
     114               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     115            END DO 
     116         END DO 
     117      ENDIF 
     118 
     119#if defined key_dynspg_ts 
     120      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
     121         DO jk = 1, jpkm1 
     122            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     123            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     124         END DO 
     125      ELSE                                            ! applied on thickness weighted velocity 
     126         DO jk = 1, jpkm1 
     127            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
     128               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
     129               &                               / fse3u_a(:,:,jk) * umask(:,:,jk) 
     130            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
     131               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
     132               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk) 
     133         END DO 
     134      ENDIF 
     135 
     136      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 
     137         ! remove barotropic velocities: 
     138         DO jk = 1, jpkm1 
     139            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 
     140            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 
     141         ENDDO 
     142         ! Add bottom stress due to barotropic component only: 
     143         DO jj = 2, jpjm1         
     144            DO ji = fs_2, fs_jpim1   ! vector opt. 
     145               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     146               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     147               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
     148               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     149               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
     150               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
     151            END DO 
     152         END DO 
     153      ENDIF 
     154#endif 
    111155 
    112156      ! 2. Vertical diffusion on u 
    … …  
    119163         DO jj = 2, jpjm1  
    120164            DO ji = fs_2, fs_jpim1   ! vector opt. 
    121                zcoef = - p2dt / fse3u(ji,jj,jk) 
     165               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point 
     166               zcoef = - p2dt / ze3ua       
    122167               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    123168               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    … …  
    128173         END DO 
    129174      END DO 
    130       DO jj = 2, jpjm1        ! Surface boudary conditions 
     175      DO jj = 2, jpjm1        ! Surface boundary conditions 
    131176         DO ji = fs_2, fs_jpim1   ! vector opt. 
    132177            zwi(ji,jj,1) = 0._wp 
    … …  
    160205      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    161206         DO ji = fs_2, fs_jpim1   ! vector opt. 
    162             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    163                &                                                       * r1_rau0 / fse3u(ji,jj,1)       ) 
     207            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
     208#if defined key_dynspg_ts 
     209            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     210               &                                      / ( ze3ua * rau0 )  
     211#else 
     212            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     213               &                                                     / ( fse3u(ji,jj,1) * rau0     ) )  
     214#endif 
    164215         END DO 
    165216      END DO 
    … …  
    167218         DO jj = 2, jpjm1    
    168219            DO ji = fs_2, fs_jpim1   ! vector opt. 
    169                zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side 
     220#if defined key_dynspg_ts 
     221               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
     222#else 
     223               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) 
     224#endif 
    170225               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    171226            END DO 
    … …  
    186241      END DO 
    187242 
     243#if ! defined key_dynspg_ts 
    188244      ! Normalization to obtain the general momentum trend ua 
    189245      DO jk = 1, jpkm1 
    … …  
    194250         END DO 
    195251      END DO 
    196  
     252#endif 
    197253 
    198254      ! 3. Vertical diffusion on v 
    … …  
    205261         DO jj = 2, jpjm1    
    206262            DO ji = fs_2, fs_jpim1   ! vector opt. 
    207                zcoef = -p2dt / fse3v(ji,jj,jk) 
     263               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point 
     264               zcoef = - p2dt / ze3va 
    208265               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    209266               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    … …  
    214271         END DO 
    215272      END DO 
    216       DO jj = 2, jpjm1        ! Surface boudary conditions 
     273      DO jj = 2, jpjm1        ! Surface boundary conditions 
    217274         DO ji = fs_2, fs_jpim1   ! vector opt. 
    218275            zwi(ji,jj,1) = 0._wp 
    … …  
    246303      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    247304         DO ji = fs_2, fs_jpim1   ! vector opt. 
    248             va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    249                &                                                       * r1_rau0 / fse3v(ji,jj,1)       ) 
     305            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
     306#if defined key_dynspg_ts             
     307            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     308               &                                      / ( ze3va * rau0 )  
     309#else 
     310            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     311               &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
     312#endif 
    250313         END DO 
    251314      END DO 
    … …  
    253316         DO jj = 2, jpjm1 
    254317            DO ji = fs_2, fs_jpim1   ! vector opt. 
    255                zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side 
     318#if defined key_dynspg_ts 
     319               zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
     320#else 
     321               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) 
     322#endif 
    256323               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    257324            END DO 
    … …  
    259326      END DO 
    260327      ! 
    261       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
     328      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
    262329         DO ji = fs_2, fs_jpim1   ! vector opt. 
    263330            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    … …  
    273340 
    274341      ! Normalization to obtain the general momentum trend va 
     342#if ! defined key_dynspg_ts 
    275343      DO jk = 1, jpkm1 
    276344         DO jj = 2, jpjm1    
    … …  
    280348         END DO 
    281349      END DO 
    282  
     350#endif 
     351 
     352      ! J. Chanut: Lines below are useless ? 
    283353      !! restore bottom layer avmu(v)  
    284354      IF( ln_bfrimp ) THEN 
    … …  
    292362            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    293363            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    294             avmu(ji,jj,ikbu+1) = zavmu(ji,jj) 
    295             avmv(ji,jj,ikbv+1) = zavmv(ji,jj) 
     364            avmu(ji,jj,ikbu+1) = 0.e0 
     365            avmv(ji,jj,ikbv+1) = 0.e0 
    296366         END DO 
    297367      END DO 
    … …  
    299369      ! 
    300370      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)  
    301       CALL wrk_dealloc( jpi,jpj, zavmu, zavmv)  
    302371      ! 
    303372      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp') 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r3764 r4292  
    88   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    10    !!---------------------------------------------------------------------- 
    11  
    12    !!---------------------------------------------------------------------- 
    13    !!   ssh_wzv        : after ssh & now vertical velocity 
    14    !!   ssh_nxt        : filter ans swap the ssh arrays 
     10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
     11   !!---------------------------------------------------------------------- 
     12 
     13   !!---------------------------------------------------------------------- 
     14   !!   ssh_nxt        : after ssh 
     15   !!   ssh_swp        : filter ans swap the ssh arrays 
     16   !!   wzv            : compute now vertical velocity 
    1517   !!---------------------------------------------------------------------- 
    1618   USE oce             ! ocean dynamics and tracers variables 
    … …  
    2022   USE divcur          ! hor. divergence and curl      (div & cur routines) 
    2123   USE iom             ! I/O library 
     24   USE restart         ! only for lrst_oce 
    2225   USE in_out_manager  ! I/O manager 
    2326   USE prtctl          ! Print control 
    … …  
    2831   USE obc_oce 
    2932   USE bdy_oce 
     33   USE bdy_par          
     34   USE bdydyn2d        ! bdy_ssh routine 
    3035   USE diaar5, ONLY:   lk_diaar5 
    3136   USE iom 
    32    USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff  
     37   USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div   ! River runoff  
     38   USE dynspg_ts,   ONLY: ln_bt_fw 
     39   USE dynspg_oce, ONLY: lk_dynspg_ts 
    3340#if defined key_agrif 
    3441   USE agrif_opa_update 
    … …  
    4451   PRIVATE 
    4552 
    46    PUBLIC   ssh_wzv    ! called by step.F90 
    4753   PUBLIC   ssh_nxt    ! called by step.F90 
     54   PUBLIC   wzv        ! called by step.F90 
     55   PUBLIC   ssh_swp    ! called by step.F90 
    4856 
    4957   !! * Substitutions 
    … …  
    5765CONTAINS 
    5866 
    59    SUBROUTINE ssh_wzv( kt )  
    60       !!---------------------------------------------------------------------- 
    61       !!                ***  ROUTINE ssh_wzv  *** 
     67   SUBROUTINE ssh_nxt( kt ) 
     68      !!---------------------------------------------------------------------- 
     69      !!                ***  ROUTINE ssh_nxt  *** 
    6270      !!                    
    63       !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity 
    64       !!              and update the now vertical coordinate (lk_vvl=T). 
    65       !! 
    66       !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    67       !!      velocity is computed by integrating the horizontal divergence   
    68       !!      from the bottom to the surface minus the scale factor evolution. 
    69       !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     71      !! ** Purpose :   compute the after ssh (ssha) 
     72      !! 
     73      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment 
     74      !!      is computed by integrating the horizontal divergence and multiply by 
     75      !!      by the time step. 
    7076      !! 
    7177      !! ** action  :   ssha    : after sea surface height 
    72       !!                wn      : now vertical velocity 
    73       !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T) 
    74       !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points 
    7578      !! 
    7679      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7780      !!---------------------------------------------------------------------- 
    78       INTEGER, INTENT(in) ::   kt   ! time step 
    79       ! 
    80       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    81       REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars 
    82       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv 
    83       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d 
    84       !!---------------------------------------------------------------------- 
    85       ! 
    86       IF( nn_timing == 1 )  CALL timing_start('ssh_wzv') 
    87       ! 
    88       CALL wrk_alloc( jpi, jpj, z2d, zhdiv )  
     81      ! 
     82      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv 
     83      INTEGER, INTENT(in) ::   kt                      ! time step 
     84      !  
     85      INTEGER             ::   jk                      ! dummy loop indice 
     86      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
     87      !!---------------------------------------------------------------------- 
     88      ! 
     89      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
     90      ! 
     91      CALL wrk_alloc( jpi, jpj, zhdiv )  
    8992      ! 
    9093      IF( kt == nit000 ) THEN 
    9194         ! 
    9295         IF(lwp) WRITE(numout,*) 
    93          IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 
     96         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 
    9497         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    9598         ! 
    96          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    97          ! 
    98          IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only) 
    99             DO jj = 1, jpjm1 
    100                DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
    101                   zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
    102                   zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
    103                   zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
    104                   sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    105                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    106                   sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    107                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    108                   sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    109                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    110                   sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    111                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    112                END DO 
    113             END DO 
    114             CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
    115             CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
    116             DO jj = 1, jpjm1 
    117                DO ji = 1, jpim1      ! NO Vector Opt. 
    118                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
    119                        &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    120                        &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    121                        &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    122                END DO 
    123             END DO 
    124             CALL lbc_lnk( sshf_n, 'F', 1. ) 
    125          ENDIF 
    126          ! 
    127       ENDIF 
    128  
    129       !                                           !------------------------------------------! 
    130       IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case) 
    131          !                                        !------------------------------------------! 
    132          DO jk = 1, jpkm1 
    133             fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays 
    134             fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 
    135             fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 
    136             ! 
    137             fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays 
    138             fse3u (:,:,jk) = fse3u_n (:,:,jk) 
    139             fse3v (:,:,jk) = fse3v_n (:,:,jk) 
    140             fse3f (:,:,jk) = fse3f_n (:,:,jk) 
    141             fse3w (:,:,jk) = fse3w_n (:,:,jk) 
    142             fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 
    143             fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
    144          END DO 
    145          ! 
    146          hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points) 
    147          hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
    148          !                                            ! now masked inverse of the ocean depth (at u- and v-points) 
    149          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) ) 
    150          hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) ) 
    151          !  
    15299      ENDIF 
    153100      ! 
    … …  
    162109      zhdiv(:,:) = 0._wp 
    163110      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    164         zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 
     111        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 
    165112      END DO 
    166113      !                                                ! Sea surface elevation time stepping 
    167114      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 
    168115      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 
    169       z1_rau0 = 0.5 / rau0 
     116      z1_rau0 = 0.5_wp * r1_rau0 
    170117      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    171118 
    … …  
    180127#endif 
    181128#if defined key_bdy 
    182       ssha(:,:) = ssha(:,:) * bdytmask(:,:) 
    183       CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm) 
    184 #endif 
     129      ! bg jchanut tschanges 
     130      ! These lines are not necessary with time splitting since 
     131      ! boundary condition on sea level is set during ts loop 
     132      IF (lk_bdy) THEN 
     133         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 
     134         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 
     135      ENDIF 
     136#endif 
     137      ! end jchanut tschanges 
    185138#if defined key_asminc 
    186139      !                                                ! Include the IAU weighted SSH increment 
    … …  
    190143      ENDIF 
    191144#endif 
    192       !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
    193       IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
    194          DO jj = 1, jpjm1 
    195             DO ji = 1, jpim1      ! NO Vector Opt. 
    196                sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
    197                   &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
    198                   &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
    199                sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
    200                   &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
    201                   &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     145 
     146      !                                           !------------------------------! 
     147      !                                           !           outputs            ! 
     148      !                                           !------------------------------! 
     149      CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
     150      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     151      ! 
     152      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
     153      ! 
     154      CALL wrk_dealloc( jpi, jpj, zhdiv )  
     155      ! 
     156      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt') 
     157      ! 
     158   END SUBROUTINE ssh_nxt 
     159 
     160    
     161   SUBROUTINE wzv( kt ) 
     162      !!---------------------------------------------------------------------- 
     163      !!                ***  ROUTINE wzv  *** 
     164      !!                    
     165      !! ** Purpose :   compute the now vertical velocity 
     166      !! 
     167      !! ** Method  : - Using the incompressibility hypothesis, the vertical  
     168      !!      velocity is computed by integrating the horizontal divergence   
     169      !!      from the bottom to the surface minus the scale factor evolution. 
     170      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     171      !! 
     172      !! ** action  :   wn      : now vertical velocity 
     173      !! 
     174      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     175      !!---------------------------------------------------------------------- 
     176      ! 
     177      INTEGER, INTENT(in) ::   kt           ! time step 
     178      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
     179      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
     180      ! 
     181      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     182      REAL(wp)            ::   z1_2dt       ! local scalars 
     183      !!---------------------------------------------------------------------- 
     184       
     185      IF( nn_timing == 1 )  CALL timing_start('wzv') 
     186      ! 
     187      IF( kt == nit000 ) THEN 
     188         ! 
     189         IF(lwp) WRITE(numout,*) 
     190         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
     191         IF(lwp) WRITE(numout,*) '~~~~~ ' 
     192         ! 
     193         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     194         ! 
     195      ENDIF 
     196      !                                           !------------------------------! 
     197      !                                           !     Now Vertical Velocity    ! 
     198      !                                           !------------------------------! 
     199      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog) 
     200      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt 
     201      ! 
     202      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
     203         CALL wrk_alloc( jpi, jpj, jpk, zhdiv )  
     204         ! 
     205         DO jk = 1, jpkm1 
     206            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
     207            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
     208            DO jj = 2, jpjm1 
     209               DO ji = fs_2, fs_jpim1   ! vector opt. 
     210                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
     211               END DO 
    202212            END DO 
    203213         END DO 
    204          CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions 
    205       ENDIF 
    206  
    207       !                                           !------------------------------! 
    208       !                                           !     Now Vertical Velocity    ! 
    209       !                                           !------------------------------! 
    210       z1_2dt = 1.e0 / z2dt 
    211       DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    212          ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    213          wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    214             &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    215             &                         * tmask(:,:,jk) * z1_2dt 
     214         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
     215         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
     216         !                             ! Same question holds for hdivn. Perhaps just for security 
     217         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     218            ! computation of w 
     219            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    & 
     220               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     221         END DO 
     222         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     223         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )  
     224      ELSE   ! z_star and linear free surface cases 
     225         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     226            ! computation of w 
     227            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   & 
     228               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     229         END DO 
     230      ENDIF 
     231 
    216232#if defined key_bdy 
    217233         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
    218234#endif 
    219       END DO 
    220  
     235      ! 
    221236      !                                           !------------------------------! 
    222237      !                                           !           outputs            ! 
    223238      !                                           !------------------------------! 
    224       CALL iom_put( "woce", wn                    )   ! vertical velocity 
    225       CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
    226       CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     239      CALL iom_put( "woce", wn )   ! vertical velocity 
    227240      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
     241         CALL wrk_alloc( jpi, jpj, z2d )  
     242         CALL wrk_alloc( jpi, jpj, jpk, z3d )  
    228243         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    229          CALL wrk_alloc( jpi,jpj,jpk, z3d ) 
    230          z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 
     244         z2d(:,:) = rau0 * e12t(:,:) 
    231245         DO jk = 1, jpk 
    232246            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
    … …  
    234248         CALL iom_put( "w_masstr" , z3d                     )   
    235249         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    236          CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 
    237       ENDIF 
    238       ! 
    239       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
    240       ! 
    241       CALL wrk_dealloc( jpi, jpj, z2d, zhdiv )  
    242       ! 
    243       IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv') 
    244       ! 
    245    END SUBROUTINE ssh_wzv 
    246  
    247  
    248    SUBROUTINE ssh_nxt( kt ) 
     250         CALL wrk_dealloc( jpi, jpj, z2d  )  
     251         CALL wrk_dealloc( jpi, jpj, jpk, z3d )  
     252      ENDIF 
     253      ! 
     254      IF( nn_timing == 1 )  CALL timing_stop('wzv') 
     255 
     256 
     257   END SUBROUTINE wzv 
     258 
     259   SUBROUTINE ssh_swp( kt ) 
    249260      !!---------------------------------------------------------------------- 
    250261      !!                    ***  ROUTINE ssh_nxt  *** 
    … …  
    252263      !! ** Purpose :   achieve the sea surface  height time stepping by  
    253264      !!              applying Asselin time filter and swapping the arrays 
    254       !!              ssha  already computed in ssh_wzv   
     265      !!              ssha  already computed in ssh_nxt   
    255266      !! 
    256267      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
    … …  
    266277      !!---------------------------------------------------------------------- 
    267278      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    268       !! 
    269       INTEGER  ::   ji, jj   ! dummy loop indices 
    270       REAL(wp) ::   zec      ! temporary scalar 
    271       !!---------------------------------------------------------------------- 
    272       ! 
    273       IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
     279      !!---------------------------------------------------------------------- 
     280      ! 
     281      IF( nn_timing == 1 )  CALL timing_start('ssh_swp') 
    274282      ! 
    275283      IF( kt == nit000 ) THEN 
    276284         IF(lwp) WRITE(numout,*) 
    277          IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 
     285         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
    278286         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    279287      ENDIF 
    280288 
    281       !                       !--------------------------! 
    282       IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points) 
    283          !                    !--------------------------! 
    284          ! 
    285          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
    286             sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now) 
    287             sshu_n(:,:) = sshu_a(:,:) 
    288             sshv_n(:,:) = sshv_a(:,:) 
    289             DO jj = 1, jpjm1                                ! ssh now at f-point 
    290                DO ji = 1, jpim1      ! NO Vector Opt. 
    291                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 & 
    292                      &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    293                      &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    294                      &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    295                END DO 
    296             END DO 
    297             CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions 
    298             ! 
    299          ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
    300             zec = atfp * rdt / rau0 
    301             DO jj = 1, jpj 
    302                DO ji = 1, jpi                               ! before <-- now filtered 
    303                   sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   & 
    304                      &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 
    305                   sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after 
    306                   sshu_n(ji,jj) = sshu_a(ji,jj) 
    307                   sshv_n(ji,jj) = sshv_a(ji,jj) 
    308                END DO 
    309             END DO 
    310             DO jj = 1, jpjm1                                ! ssh now at f-point 
    311                DO ji = 1, jpim1      ! NO Vector Opt. 
    312                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 & 
    313                      &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    314                      &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    315                      &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    316                END DO 
    317             END DO 
    318             CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions 
    319             ! 
    320             DO jj = 1, jpjm1                                ! ssh before at u- & v-points 
    321                DO ji = 1, jpim1      ! NO Vector Opt. 
    322                   sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
    323                      &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    324                      &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    325                   sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
    326                      &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    327                      &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    328                END DO 
    329             END DO 
    330             CALL lbc_lnk( sshu_b, 'U', 1. ) 
    331             CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions 
    332             ! 
    333          ENDIF 
    334          !                    !--------------------------! 
    335       ELSE                    !        fixed levels      !     (ssh at t-point only) 
    336          !                    !--------------------------! 
    337          ! 
    338          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
    339             sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now) 
    340             ! 
    341          ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap 
    342             DO jj = 1, jpj 
    343                DO ji = 1, jpi                               ! before <-- now filtered 
    344                   sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 
    345                   sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after 
    346                END DO 
    347             END DO 
    348          ENDIF 
    349          ! 
     289# if defined key_dynspg_ts 
     290      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter 
     291# else 
     292      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter 
     293#endif 
     294         sshb(:,:) = sshn(:,:)                           ! before <-- now 
     295         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now) 
     296      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
     297         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered 
     298         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     299         sshn(:,:) = ssha(:,:)                           ! now <-- after 
    350300      ENDIF 
    351301      ! 
    … …  
    357307      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    358308      ! 
    359       IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt') 
    360       ! 
    361    END SUBROUTINE ssh_nxt 
     309      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp') 
     310      ! 
     311   END SUBROUTINE ssh_swp 
    362312 
    363313   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.