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.
dynspg_ts.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN/dynspg_ts.F90 @ 10386

Last change on this file since 10386 was 10386, checked in by smasson, 5 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: safer use of l_full_nf_update (2 lines north fold exchange), see #2133

  • Property svn:keywords set to Id
File size: 77.2 KB
RevLine 
[358]1MODULE dynspg_ts
[9023]2
3   !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
4
[358]5   !!======================================================================
[6140]6   !!                   ***  MODULE  dynspg_ts  ***
7   !! Ocean dynamics:  surface pressure gradient trend, split-explicit scheme
8   !!======================================================================
[1502]9   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
10   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
11   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
12   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
13   !!              -   ! 2008-01  (R. Benshila)  change averaging method
14   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
[2528]15   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[2724]16   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
[4292]17   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
18   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
[5930]19   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
[7646]20   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
[9019]21   !!             4.0  ! 2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
[2724]22   !!---------------------------------------------------------------------
[6140]23
[358]24   !!----------------------------------------------------------------------
[6140]25   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
26   !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
27   !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
28   !!   ts_rst         : read/write time-splitting fields in restart file
[358]29   !!----------------------------------------------------------------------
30   USE oce             ! ocean dynamics and tracers
31   USE dom_oce         ! ocean space and time domain
[888]32   USE sbc_oce         ! surface boundary condition: ocean
[9019]33   USE zdf_oce         ! vertical physics: variables
34   USE zdfdrg          ! vertical physics: top/bottom drag coef.
[5120]35   USE sbcisf          ! ice shelf variable (fwfisf)
[6140]36   USE sbcapr          ! surface boundary condition: atmospheric pressure
37   USE dynadv    , ONLY: ln_dynadv_vec
[9528]38   USE dynvor          ! vortivity scheme indicators
[358]39   USE phycst          ! physical constants
40   USE dynvor          ! vorticity term
[6152]41   USE wet_dry         ! wetting/drying flux limter
[7646]42   USE bdy_oce         ! open boundary
[5930]43   USE bdytides        ! open boundary condition data
[3294]44   USE bdydyn2d        ! open boundary conditions on barotropic variables
[4292]45   USE sbctide         ! tides
46   USE updtide         ! tide potential
[7646]47   USE sbcwave         ! surface wave
[9019]48   USE diatmb          ! Top,middle,bottom output
49#if defined key_agrif
[9570]50   USE agrif_oce_interp ! agrif
[9124]51   USE agrif_oce
[9019]52#endif
53#if defined key_asminc   
54   USE asminc          ! Assimilation increment
55#endif
[6140]56   !
57   USE in_out_manager  ! I/O manager
[358]58   USE lib_mpp         ! distributed memory computing library
59   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl          ! Print control
[2715]61   USE iom             ! IOM library
[4292]62   USE restart         ! only for lrst_oce
[9023]63   USE diatmb          ! Top,middle,bottom output
[358]64
65   IMPLICIT NONE
66   PRIVATE
67
[9124]68   PUBLIC dyn_spg_ts        ! called by dyn_spg
69   PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
[358]70
[9019]71   !! Time filtered arrays at baroclinic time step:
72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step
[9124]73   !
[9023]74   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
75   REAL(wp),SAVE :: rdtbt       ! Barotropic time step
[9019]76   !
77   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields
[9124]78   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwz                ! ff_f/h at F points
79   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftnw, ftne         ! triad of coriolis parameter
80   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftsw, ftse         ! (only used with een vorticity scheme)
[4292]81
[9043]82   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios
83   REAL(wp) ::   r1_8  = 0.125_wp         !
84   REAL(wp) ::   r1_4  = 0.25_wp          !
85   REAL(wp) ::   r1_2  = 0.5_wp           !
[508]86
[358]87   !! * Substitutions
88#  include "vectopt_loop_substitute.h90"
[2715]89   !!----------------------------------------------------------------------
[9598]90   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]91   !! $Id$
[10068]92   !! Software governed by the CeCILL license (see ./LICENSE)
[2715]93   !!----------------------------------------------------------------------
[358]94CONTAINS
95
[2715]96   INTEGER FUNCTION dyn_spg_ts_alloc()
97      !!----------------------------------------------------------------------
98      !!                  ***  routine dyn_spg_ts_alloc  ***
99      !!----------------------------------------------------------------------
[6140]100      INTEGER :: ierr(3)
[4292]101      !!----------------------------------------------------------------------
102      ierr(:) = 0
[6140]103      !
104      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
105      !
[9528]106      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
107         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
108         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
[6140]109         !
110      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
111      !
112      dyn_spg_ts_alloc = MAXVAL( ierr(:) )
113      !
[10297]114      IF( lk_mpp                )   CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
[5930]115      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
[2715]116      !
117   END FUNCTION dyn_spg_ts_alloc
118
[5836]119
[358]120   SUBROUTINE dyn_spg_ts( kt )
121      !!----------------------------------------------------------------------
122      !!
[6140]123      !! ** Purpose : - Compute the now trend due to the explicit time stepping
124      !!              of the quasi-linear barotropic system, and add it to the
125      !!              general momentum trend.
[358]126      !!
[6140]127      !! ** Method  : - split-explicit schem (time splitting) :
[4374]128      !!      Barotropic variables are advanced from internal time steps
129      !!      "n"   to "n+1" if ln_bt_fw=T
130      !!      or from
131      !!      "n-1" to "n+1" if ln_bt_fw=F
132      !!      thanks to a generalized forward-backward time stepping (see ref. below).
[358]133      !!
[4374]134      !! ** Action :
135      !!      -Update the filtered free surface at step "n+1"      : ssha
136      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
[9023]137      !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
[4374]138      !!      These are used to advect tracers and are compliant with discrete
139      !!      continuity equation taken at the baroclinic time steps. This
140      !!      ensures tracers conservation.
[6140]141      !!      - (ua, va) momentum trend updated with barotropic component.
[358]142      !!
[6140]143      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
[358]144      !!---------------------------------------------------------------------
[1502]145      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
[2715]146      !
[9554]147      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
[9019]148      LOGICAL  ::   ll_fw_start           ! =T : forward integration
[9554]149      LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
[9019]150      LOGICAL  ::   ll_tmp1, ll_tmp2      ! local logical variables used in W/D
151      INTEGER  ::   ikbu, iktu, noffset   ! local integers
152      INTEGER  ::   ikbv, iktv            !   -      -
[9554]153      REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars
[9528]154      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      -
155      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      -
156      REAL(wp) ::   za0, za1, za2, za3              !   -      -
157      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      -
[9019]158      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
159      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
160      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
[9528]161      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e
[9019]162      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
163      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
[3294]164      !
[9023]165      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.
166
167      INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point
168
169      REAL(wp) ::   zepsilon, zgamma            !   -      -
[9019]170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy   ! Wetting/Dying gravity filter coef.
[9023]171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
[358]173      !!----------------------------------------------------------------------
[3294]174      !
[9023]175      IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
176      !                                         !* Allocate temporary arrays
177      IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
[3294]178      !
[6140]179      zmdi=1.e+20                               !  missing data indicator for masking
[9019]180      !
[9023]181      zwdramp = r_rn_wdmin1               ! simplest ramp
182!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
183      !                                         ! reciprocal of baroclinic time step
[6140]184      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
185      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt
[4292]186      ENDIF
[9043]187      r1_2dt_b = 1.0_wp / z2dt_bf 
[4292]188      !
[9023]189      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart
[4292]190      ll_fw_start = .FALSE.
[9023]191      !                                         ! time offset in steps for bdy data update
[6140]192      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro
193      ELSE                       ;   noffset =   0 
194      ENDIF
[4292]195      !
[9023]196      IF( kt == nit000 ) THEN                   !* initialisation
[508]197         !
[358]198         IF(lwp) WRITE(numout,*)
199         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
200         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
[4354]201         IF(lwp) WRITE(numout,*)
[1502]202         !
[6140]203         IF( neuler == 0 )   ll_init=.TRUE.
[1502]204         !
[6140]205         IF( ln_bt_fw .OR. neuler == 0 ) THEN
206            ll_fw_start =.TRUE.
207            noffset     = 0
[4292]208         ELSE
[6140]209            ll_fw_start =.FALSE.
[4292]210         ENDIF
211         !
212         ! Set averaging weights and cycle length:
[6140]213         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]214         !
215      ENDIF
216      !
[9019]217      IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities)
218         DO jj = 2, jpjm1
219            DO ji = fs_2, fs_jpim1   ! vector opt.
[9112]220               zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
221               zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
[9019]222            END DO
223         END DO
224      ELSE                    ! bottom friction only
225         DO jj = 2, jpjm1
226            DO ji = fs_2, fs_jpim1   ! vector opt.
[9112]227               zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
228               zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
[9019]229            END DO
230         END DO
231      ENDIF
232      !
[4292]233      ! Set arrays to remove/compute coriolis trend.
234      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
235      ! Note that these arrays are also used during barotropic loop. These are however frozen
[4374]236      ! although they should be updated in the variable volume case. Not a big approximation.
[4292]237      ! To remove this approximation, copy lines below inside barotropic loop
[4374]238      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
[4292]239      !
[6140]240      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
[9528]241         !
242         SELECT CASE( nvor_scheme )
243         CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
[7646]244            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
[5836]245            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
246               DO jj = 1, jpjm1
247                  DO ji = 1, jpim1
[6140]248                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
249                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp 
[7646]250                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
[5836]251                  END DO
[5032]252               END DO
[5836]253            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
254               DO jj = 1, jpjm1
255                  DO ji = 1, jpim1
[9528]256                     zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
257                        &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
258                        &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
259                        &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
[7646]260                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
[5836]261                  END DO
[4292]262               END DO
[5836]263            END SELECT
[10170]264            CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
[5836]265            !
[7753]266            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[358]267            DO jj = 2, jpj
[5836]268               DO ji = 2, jpi
[4292]269                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
270                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
271                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
272                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
[358]273               END DO
274            END DO
[5836]275            !
[9528]276         CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
277            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
278            DO jj = 2, jpj
279               DO ji = 2, jpi
280                  z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
281                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
282                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
283                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
284                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
285               END DO
286            END DO
287            !
288         CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
289            !
[7753]290            zwz(:,:) = 0._wp
291            zhf(:,:) = 0._wp
[7646]292           
[9124]293!!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
[7646]294!!gm    A priori a better value should be something like :
295!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
296!!gm                     divided by the sum of the corresponding mask
297!!gm
298!!           
[9528]299            IF( .NOT.ln_sco ) THEN
[7646]300 
301   !!gm  agree the JC comment  : this should be done in a much clear way
302 
303   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
304   !     Set it to zero for the time being
305   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
306   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
307   !              ENDIF
308   !              zhf(:,:) = gdepw_0(:,:,jk+1)
[9528]309               !
310            ELSE
311               !
312               !zhf(:,:) = hbatf(:,:)
313               DO jj = 1, jpjm1
314                  DO ji = 1, jpim1
315                     zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
316                        &              + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
317                        &       / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
318                        &              + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
319                  END DO
320               END DO
321            ENDIF
322            !
323            DO jj = 1, jpjm1
324               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
325            END DO
326            !
[4292]327            DO jk = 1, jpkm1
328               DO jj = 1, jpjm1
[7753]329                  zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
[4292]330               END DO
331            END DO
[10170]332            CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
[4292]333            ! JC: TBC. hf should be greater than 0
334            DO jj = 1, jpj
335               DO ji = 1, jpi
[4370]336                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
[4292]337               END DO
338            END DO
[7753]339            zwz(:,:) = ff_f(:,:) * zwz(:,:)
[9528]340         END SELECT
[508]341      ENDIF
[1502]342      !
[4292]343      ! If forward start at previous time step, and centered integration,
344      ! then update averaging weights:
[5836]345      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
[4292]346         ll_fw_start=.FALSE.
[9019]347         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]348      ENDIF
349                         
[358]350      ! -----------------------------------------------------------------------------
351      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
352      ! -----------------------------------------------------------------------------
[1502]353      !     
[4292]354      !
[4354]355      !                                   !* e3*d/dt(Ua) (Vertically integrated)
[4292]356      !                                   ! --------------------------------------------------
[7753]357      zu_frc(:,:) = 0._wp
358      zv_frc(:,:) = 0._wp
[1502]359      !
360      DO jk = 1, jpkm1
[7753]361         zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
362         zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
[1502]363      END DO
[4292]364      !
[7753]365      zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
366      zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
[4292]367      !
[7753]368      !
[1502]369      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
[4292]370      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
[1502]371         DO jj = 2, jpjm1
372            DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]373               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
374               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
[1502]375            END DO
[358]376         END DO
[1502]377      END DO
[7646]378     
379!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
380!!gm  Is it correct to do so ?   I think so...
381     
382     
[4292]383      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
384      !                                   ! --------------------------------------------------------
[9528]385      !
[7753]386      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes
387      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
[1502]388      !
[9528]389      SELECT CASE( nvor_scheme )
390      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
[358]391         DO jj = 2, jpjm1
[9528]392            DO ji = 2, jpim1   ! vector opt.
393               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    &
394                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
395                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
396                  !
397               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    &
398                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
399                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
400            END DO 
401         END DO 
402         !         
403      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
404         DO jj = 2, jpjm1
[358]405            DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]406               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
407               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
408               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
409               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
[358]410               ! energy conserving formulation for planetary vorticity term
[9043]411               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
412               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
[358]413            END DO
414         END DO
[508]415         !
[9528]416      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
[358]417         DO jj = 2, jpjm1
418            DO ji = fs_2, fs_jpim1   ! vector opt.
[9043]419               zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
[5836]420                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
[9043]421               zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
[5836]422                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
[4292]423               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
424               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
[358]425            END DO
426         END DO
[508]427         !
[9528]428      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
[358]429         DO jj = 2, jpjm1
430            DO ji = fs_2, fs_jpim1   ! vector opt.
[9043]431               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
[5836]432                &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
433                &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
434                &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
[9043]435               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
[5836]436                &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
437                &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
438                &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
[358]439            END DO
440         END DO
[508]441         !
[9528]442      END SELECT
[4292]443      !
[1502]444      !                                   !* Right-Hand-Side of the barotropic momentum equation
445      !                                   ! ----------------------------------------------------
[6140]446      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
[9528]447         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters
448            DO jj = 2, jpjm1
449               DO ji = 2, jpim1 
450                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
451                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
452                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
453                     &                                                         > rn_wdmin1 + rn_wdmin2
454                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
455                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
456                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
457                  IF(ll_tmp1) THEN
458                     zcpx(ji,jj) = 1.0_wp
459                  ELSEIF(ll_tmp2) THEN
460                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
461                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
462                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
463                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
464                  ELSE
465                     zcpx(ji,jj) = 0._wp
466                  ENDIF
467                  !
468                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
469                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
470                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
[9023]471                     &                                                       > rn_wdmin1 + rn_wdmin2
[9528]472                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
473                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
474                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
[9023]475 
[9528]476                  IF(ll_tmp1) THEN
477                     zcpy(ji,jj) = 1.0_wp
478                  ELSE IF(ll_tmp2) THEN
479                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
480                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
481                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
482                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
483                  ELSE
484                     zcpy(ji,jj) = 0._wp
485                  ENDIF
486               END DO
487            END DO
488            !
489            DO jj = 2, jpjm1
490               DO ji = 2, jpim1
491                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
492                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
493                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
494                     &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
495               END DO
496            END DO
497            !
[6152]498         ELSE
[9019]499            !
500            DO jj = 2, jpjm1
501               DO ji = fs_2, fs_jpim1   ! vector opt.
502                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
503                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
504               END DO
505            END DO
506         ENDIF
507         !
[1502]508      ENDIF
[9019]509      !
[4292]510      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
[358]511         DO ji = fs_2, fs_jpim1
[6140]512             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
513             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
[3294]514          END DO
[4292]515      END DO 
516      !
[9023]517      !                                         ! Add bottom stress contribution from baroclinic velocities:     
518      IF (ln_bt_fw) THEN
[4292]519         DO jj = 2, jpjm1                         
520            DO ji = fs_2, fs_jpim1   ! vector opt.
521               ikbu = mbku(ji,jj)       
522               ikbv = mbkv(ji,jj)   
523               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
524               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
525            END DO
526         END DO
[3294]527      ELSE
[4292]528         DO jj = 2, jpjm1
529            DO ji = fs_2, fs_jpim1   ! vector opt.
530               ikbu = mbku(ji,jj)       
531               ikbv = mbkv(ji,jj)   
532               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
533               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
534            END DO
535         END DO
536      ENDIF
[1502]537      !
[4292]538      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
[9023]539      IF( ln_wd_il ) THEN
[9075]540         zztmp = -1._wp / rdtbt
[9045]541         DO jj = 2, jpjm1
542            DO ji = fs_2, fs_jpim1   ! vector opt.
[9112]543               zu_frc(ji,jj) = zu_frc(ji,jj) + & 
544               & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) *  wdrampu(ji,jj)
545               zv_frc(ji,jj) = zv_frc(ji,jj) + & 
546               & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) *  wdrampv(ji,jj)
[9045]547            END DO
548         END DO
[9023]549      ELSE
[9045]550         DO jj = 2, jpjm1
551            DO ji = fs_2, fs_jpim1   ! vector opt.
[9112]552               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)
553               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)
[9045]554            END DO
555         END DO
[9023]556      END IF
557      !
[9019]558      IF( ln_isfcav ) THEN       ! Add TOP stress contribution from baroclinic velocities:     
559         IF( ln_bt_fw ) THEN
560            DO jj = 2, jpjm1
561               DO ji = fs_2, fs_jpim1   ! vector opt.
562                  iktu = miku(ji,jj)
563                  iktv = mikv(ji,jj)
564                  zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
565                  zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
566               END DO
567            END DO
568         ELSE
569            DO jj = 2, jpjm1
570               DO ji = fs_2, fs_jpim1   ! vector opt.
571                  iktu = miku(ji,jj)
572                  iktv = mikv(ji,jj)
573                  zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
574                  zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
575               END DO
576            END DO
577         ENDIF
578         !
579         ! Note that the "unclipped" top friction parameter is used even with explicit drag
580         DO jj = 2, jpjm1             
581            DO ji = fs_2, fs_jpim1   ! vector opt.
[9112]582               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)
583               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)
[9019]584            END DO
585         END DO
586      ENDIF
[6140]587      !       
[9019]588      IF( ln_bt_fw ) THEN                        ! Add wind forcing
589         DO jj = 2, jpjm1
590            DO ji = fs_2, fs_jpim1   ! vector opt.
591               zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
592               zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
593            END DO
594         END DO
[2724]595      ELSE
[9043]596         zztmp = r1_rau0 * r1_2
[9019]597         DO jj = 2, jpjm1
598            DO ji = fs_2, fs_jpim1   ! vector opt.
599               zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
600               zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
601            END DO
602         END DO
[4292]603      ENDIF 
604      !
[9019]605      IF( ln_apr_dyn ) THEN                     ! Add atm pressure forcing
606         IF( ln_bt_fw ) THEN
[4292]607            DO jj = 2, jpjm1             
608               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]609                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
610                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
[4292]611                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
612                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
613               END DO
614            END DO
615         ELSE
[9043]616            zztmp = grav * r1_2
[4292]617            DO jj = 2, jpjm1             
618               DO ji = fs_2, fs_jpim1   ! vector opt.
[9019]619                  zu_spg = zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
620                      &             + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
621                  zv_spg = zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
622                      &             + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
[4292]623                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
624                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
625               END DO
626            END DO
627         ENDIF
[2724]628      ENDIF
[4292]629      !                                   !* Right-Hand-Side of the barotropic ssh equation
630      !                                   ! -----------------------------------------------
631      !                                         ! Surface net water flux and rivers
632      IF (ln_bt_fw) THEN
[9019]633         zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
[4292]634      ELSE
[9043]635         zztmp = r1_rau0 * r1_2
[9019]636         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
637                &                 + fwfisf(:,:) + fwfisf_b(:,:)                     )
[4292]638      ENDIF
[7646]639      !
640      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary
[7753]641         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
[7646]642      ENDIF
643      !
[4292]644#if defined key_asminc
645      !                                         ! Include the IAU weighted SSH increment
646      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
[7753]647         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
[4292]648      ENDIF
649#endif
[5656]650      !                                   !* Fill boundary data arrays for AGRIF
651      !                                   ! ------------------------------------
[4486]652#if defined key_agrif
653         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
654#endif
[4292]655      !
[358]656      ! -----------------------------------------------------------------------
[4292]657      !  Phase 2 : Integration of the barotropic equations
[358]658      ! -----------------------------------------------------------------------
[1502]659      !
660      !                                             ! ==================== !
661      !                                             !    Initialisations   !
[4292]662      !                                             ! ==================== ! 
[4370]663      ! Initialize barotropic variables:     
[4770]664      IF( ll_init )THEN
[7753]665         sshbb_e(:,:) = 0._wp
666         ubb_e  (:,:) = 0._wp
667         vbb_e  (:,:) = 0._wp
668         sshb_e (:,:) = 0._wp
669         ub_e   (:,:) = 0._wp
670         vb_e   (:,:) = 0._wp
[4700]671      ENDIF
[6152]672
[4700]673      !
[4370]674      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
[7753]675         sshn_e(:,:) =    sshn(:,:)           
676         un_e  (:,:) =    un_b(:,:)           
677         vn_e  (:,:) =    vn_b(:,:)
678         !
679         hu_e  (:,:) =    hu_n(:,:)       
680         hv_e  (:,:) =    hv_n(:,:) 
681         hur_e (:,:) = r1_hu_n(:,:)   
682         hvr_e (:,:) = r1_hv_n(:,:)
[4370]683      ELSE                                ! CENTRED integration: start from BEFORE fields
[7753]684         sshn_e(:,:) =    sshb(:,:)
685         un_e  (:,:) =    ub_b(:,:)         
686         vn_e  (:,:) =    vb_b(:,:)
687         !
688         hu_e  (:,:) =    hu_b(:,:)       
689         hv_e  (:,:) =    hv_b(:,:) 
690         hur_e (:,:) = r1_hu_b(:,:)   
691         hvr_e (:,:) = r1_hv_b(:,:)
[4292]692      ENDIF
693      !
694      !
[4370]695      !
[4292]696      ! Initialize sums:
[7753]697      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
698      va_b  (:,:) = 0._wp
699      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
700      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
701      vn_adv(:,:) = 0._wp
[9528]702      !
703      IF( ln_wd_dl ) THEN
[9023]704         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)
705         zvwdmask(:,:) = 0._wp  !
[9528]706         zuwdav2 (:,:) = 0._wp 
707         zvwdav2 (:,:) = 0._wp   
[9023]708      END IF 
709
[9528]710      !                                             ! ==================== !
[4292]711      DO jn = 1, icycle                             !  sub-time-step loop  !
[1502]712         !                                          ! ==================== !
[10386]713         !
714         l_full_nf_update = jn == icycle   ! false: disable full North fold update (performances) for jn = 1 to icycle-1
715         !                                                !  ------------------
[3294]716         !                                                !* Update the forcing (BDY and tides)
[1502]717         !                                                !  ------------------
[4292]718         ! Update only tidal forcing at open boundaries
[7646]719         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
720         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
[4292]721         !
722         ! Set extrapolation coefficients for predictor step:
723         IF ((jn<3).AND.ll_init) THEN      ! Forward           
724           za1 = 1._wp                                         
725           za2 = 0._wp                       
726           za3 = 0._wp                       
727         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
728           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
729           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
730           za3 =  0.281105_wp              ! za3 = bet
731         ENDIF
[367]732
[4292]733         ! Extrapolate barotropic velocities at step jit+0.5:
[7753]734         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
735         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
[4292]736
[6140]737         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
[4292]738            !                                             !  ------------------
739            ! Extrapolate Sea Level at step jit+0.5:
[7753]740            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[9023]741           
742            ! set wetting & drying mask at tracer points for this barotropic sub-step
743            IF ( ln_wd_dl ) THEN 
[9528]744               !
[9023]745               IF ( ln_wd_dl_rmp ) THEN
746                  DO jj = 1, jpj                                 
747                     DO ji = 1, jpi   ! vector opt. 
748                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
749!                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN
750                           ztwdmask(ji,jj) = 1._wp
751                        ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
752                           ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) ) 
753                        ELSE
754                           ztwdmask(ji,jj) = 0._wp
755                        END IF
756                     END DO
757                  END DO
758               ELSE
759                  DO jj = 1, jpj                                 
760                     DO ji = 1, jpi   ! vector opt. 
761                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
762                           ztwdmask(ji,jj) = 1._wp
763                        ELSE
764                           ztwdmask(ji,jj) = 0._wp
[9528]765                        ENDIF
[9023]766                     END DO
767                  END DO
[9528]768               ENDIF 
769               !
770            ENDIF 
771            !
[4292]772            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
773               DO ji = 2, fs_jpim1   ! Vector opt.
[9043]774                  zwx(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
[5836]775                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
776                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
[9043]777                  zwy(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
[5836]778                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
779                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
[4292]780               END DO
781            END DO
[10170]782            CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp )
[4292]783            !
[9528]784            zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
785            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:)
786            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
[4370]787         ELSE
[9528]788            zhup2_e(:,:) = hu_n(:,:)
789            zhvp2_e(:,:) = hv_n(:,:)
790            zhtp2_e(:,:) = ht_n(:,:)
[4292]791         ENDIF
792         !                                                !* after ssh
[1502]793         !                                                !  -----------
[4292]794         ! One should enforce volume conservation at open boundaries here
795         ! considering fluxes below:
796         !
[7753]797         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
798         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
[4486]799         !
800#if defined key_agrif
[6140]801         ! Set fluxes during predictor step to ensure volume conservation
802         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
[4486]803            IF((nbondi == -1).OR.(nbondi == 2)) THEN
[9019]804               DO jj = 1, jpj
[9116]805                  zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
[9845]806                  zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
[4486]807               END DO
808            ENDIF
809            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
810               DO jj=1,jpj
[9116]811                  zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
[9845]812                  zwy(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj)
[4486]813               END DO
814            ENDIF
815            IF((nbondj == -1).OR.(nbondj == 2)) THEN
816               DO ji=1,jpi
[9116]817                  zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
[9845]818                  zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
[4486]819               END DO
820            ENDIF
821            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
822               DO ji=1,jpi
[9116]823                  zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
[9845]824                  zwx(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1)
[4486]825               END DO
826            ENDIF
827         ENDIF
828#endif
[9528]829         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
[9023]830
831         IF ( ln_wd_dl ) THEN 
[9528]832            !
833            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
834            !
[9023]835            DO jj = 1, jpjm1                                 
836               DO ji = 1, jpim1   
837                  IF ( zwx(ji,jj) > 0.0 ) THEN
838                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj) 
839                  ELSE
840                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 
841                  END IF
842                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
843                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
844
845                  IF ( zwy(ji,jj) > 0.0 ) THEN
846                     zvwdmask(ji, jj) = ztwdmask(ji, jj  )
847                  ELSE
848                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 
849                  END IF
850                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 
851                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
852               END DO
853            END DO
[9528]854            !
855         ENDIF   
[9023]856         
[4486]857         ! Sum over sub-time-steps to compute advective velocities
858         za2 = wgtbtp2(jn)
[7753]859         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
860         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
[9023]861         
862         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
863         IF ( ln_wd_dl_bc ) THEN
864            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
865            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
866         END IF 
867
[4486]868         ! Set next sea level:
[4292]869         DO jj = 2, jpjm1                                 
[358]870            DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]871               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
[5836]872                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
[358]873            END DO
874         END DO
[7753]875         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
876         
[10170]877         CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T',  1._wp )
[4292]878
[6140]879         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
[7646]880         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
[4292]881#if defined key_agrif
[6140]882         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
[4292]883#endif
884         
885         ! Sea Surface Height at u-,v-points (vvl case only)
[6140]886         IF( .NOT.ln_linssh ) THEN                               
[4292]887            DO jj = 2, jpjm1
888               DO ji = 2, jpim1      ! NO Vector Opt.
[9043]889                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
[6140]890                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
891                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
[9043]892                  zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
[6140]893                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
894                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
[4292]895               END DO
[358]896            END DO
[10170]897            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
[4292]898         ENDIF   
899         !                                 
900         ! Half-step back interpolation of SSH for surface pressure computation:
901         !----------------------------------------------------------------------
902         IF ((jn==1).AND.ll_init) THEN
903           za0=1._wp                        ! Forward-backward
904           za1=0._wp                           
905           za2=0._wp
906           za3=0._wp
907         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
908           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
909           za1=-0.1666666666666_wp          ! za1 = gam
910           za2= 0.0833333333333_wp          ! za2 = eps
911           za3= 0._wp             
912         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
[9023]913            IF (rn_bt_alpha==0._wp) THEN
914               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps
915               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps
916               za2=0.088_wp                 ! za2 = gam
917               za3=0.013_wp                 ! za3 = eps
918            ELSE
919               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
920               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
921               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
922               za1 = 1._wp - za0 - zgamma - zepsilon
923               za2 = zgamma
924               za3 = zepsilon
925            ENDIF
[4292]926         ENDIF
[6140]927         !
[9528]928         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
929            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
930         
[9023]931         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters
[6152]932           DO jj = 2, jpjm1
[7646]933              DO ji = 2, jpim1 
934                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
[9023]935                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            &
936                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
[7646]937                     &                                                             > rn_wdmin1 + rn_wdmin2
938                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
939                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
[9023]940                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
[7646]941   
942                IF(ll_tmp1) THEN
943                  zcpx(ji,jj) = 1.0_wp
944                ELSE IF(ll_tmp2) THEN
945                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here
[9023]946                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
[7646]947                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) )
948                ELSE
949                  zcpx(ji,jj) = 0._wp
[9528]950                ENDIF
951                !
[7646]952                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
[9023]953                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            &
954                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
[7646]955                     &                                                             > rn_wdmin1 + rn_wdmin2
956                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
957                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
[9023]958                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
[7646]959   
960                IF(ll_tmp1) THEN
961                  zcpy(ji,jj) = 1.0_wp
[9528]962                ELSEIF(ll_tmp2) THEN
[7646]963                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
[9023]964                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
[7646]965                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) )
966                ELSE
967                  zcpy(ji,jj) = 0._wp
[9528]968                ENDIF
[6152]969              END DO
[7646]970           END DO
[9528]971         ENDIF
[1502]972         !
[4292]973         ! Compute associated depths at U and V points:
[9023]974         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form
[4292]975            !                                       
976            DO jj = 2, jpjm1                           
977               DO ji = 2, jpim1
[9043]978                  zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
[5836]979                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
980                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
[9043]981                  zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
[5836]982                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
983                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
[4292]984                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
985                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
986               END DO
987            END DO
[9528]988            !
[4292]989         ENDIF
990         !
991         ! Add Coriolis trend:
[6140]992         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]993         ! at each time step. We however keep them constant here for optimization.
994         ! Recall that zwx and zwy arrays hold fluxes at this stage:
995         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
996         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
997         !
[9528]998         SELECT CASE( nvor_scheme )
999         CASE( np_ENT )             ! energy conserving scheme (t-point)
1000         DO jj = 2, jpjm1
1001            DO ji = 2, jpim1   ! vector opt.
1002
1003               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
1004               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
1005           
1006               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   &
1007                  &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   &
1008                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   )
1009                  !
1010               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1011                  &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   & 
1012                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   ) 
1013            END DO 
1014         END DO 
1015         !         
1016         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point)
[358]1017            DO jj = 2, jpjm1
1018               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]1019                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
1020                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
1021                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
1022                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
[9043]1023                  zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1024                  zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
[358]1025               END DO
1026            END DO
[508]1027            !
[9528]1028         CASE( np_ENS )             ! enstrophy conserving scheme (f-point)
[358]1029            DO jj = 2, jpjm1
1030               DO ji = fs_2, fs_jpim1   ! vector opt.
[9043]1031                  zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
[5836]1032                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
[9043]1033                  zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
[5836]1034                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
[4292]1035                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1036                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
[358]1037               END DO
1038            END DO
[508]1039            !
[9528]1040         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f)
[358]1041            DO jj = 2, jpjm1
1042               DO ji = fs_2, fs_jpim1   ! vector opt.
[9528]1043                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  &
1044                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  &
1045                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  & 
1046                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  )
1047                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  & 
1048                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  &
1049                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  & 
1050                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  )
[358]1051               END DO
1052            END DO
[508]1053            !
[9528]1054         END SELECT
[4292]1055         !
1056         ! Add tidal astronomical forcing if defined
[7646]1057         IF ( ln_tide .AND. ln_tide_pot ) THEN
[4292]1058            DO jj = 2, jpjm1
1059               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]1060                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1061                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
[4292]1062                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1063                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1064               END DO
1065            END DO
1066         ENDIF
1067         !
[9023]1068         ! Add bottom stresses:
1069!jth do implicitly instead
1070         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[9045]1071            DO jj = 2, jpjm1
1072               DO ji = fs_2, fs_jpim1   ! vector opt.
1073                  zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1074                  zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1075               END DO
1076            END DO
[9023]1077         ENDIF 
[4292]1078         !
1079         ! Surface pressure trend:
[9023]1080         IF( ln_wd_il ) THEN
[6152]1081           DO jj = 2, jpjm1
1082              DO ji = 2, jpim1 
1083                 ! Add surface pressure gradient
1084                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1085                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
[9023]1086                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 
1087                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
[6152]1088              END DO
1089           END DO
1090         ELSE
1091           DO jj = 2, jpjm1
1092              DO ji = fs_2, fs_jpim1   ! vector opt.
1093                 ! Add surface pressure gradient
1094                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1095                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
[9023]1096                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1097                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
[6152]1098              END DO
1099           END DO
1100         END IF
1101
[4292]1102         !
1103         ! Set next velocities:
[9023]1104         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
[4292]1105            DO jj = 2, jpjm1
1106               DO ji = fs_2, fs_jpim1   ! vector opt.
[5930]1107                  ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
[4292]1108                            &     + rdtbt * (                      zwx(ji,jj)   &
1109                            &                                 + zu_trd(ji,jj)   &
1110                            &                                 + zu_frc(ji,jj) ) & 
[6140]1111                            &   ) * ssumask(ji,jj)
[358]1112
[5930]1113                  va_e(ji,jj) = (                                 vn_e(ji,jj)   &
[4292]1114                            &     + rdtbt * (                      zwy(ji,jj)   &
1115                            &                                 + zv_trd(ji,jj)   &
1116                            &                                 + zv_frc(ji,jj) ) &
[6140]1117                            &   ) * ssvmask(ji,jj)
[9023]1118 
[4292]1119               END DO
1120            END DO
[6140]1121            !
[9023]1122         ELSE                           !* Flux form
[4292]1123            DO jj = 2, jpjm1
1124               DO ji = fs_2, fs_jpim1   ! vector opt.
[3294]1125
[9023]1126                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1127                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1128
[6152]1129                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1130                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1131
[5930]1132                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
[4292]1133                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
1134                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
[6140]1135                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
[4292]1136                            &   ) * zhura
[358]1137
[5930]1138                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
[4292]1139                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
1140                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
[6140]1141                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
[4292]1142                            &   ) * zhvra
[592]1143               END DO
1144            END DO
[4292]1145         ENDIF
[10345]1146!jth implicit bottom friction:
1147         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1148            DO jj = 2, jpjm1
1149               DO ji = fs_2, fs_jpim1   ! vector opt.
1150                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1151                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1152               END DO
1153            END DO
1154         ENDIF
[9023]1155
1156         
[6140]1157         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
[9023]1158            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1159            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
[7753]1160            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1161            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
[1502]1162            !
[1438]1163         ENDIF
[6140]1164         !                                             !* domain lateral boundary
[10170]1165         CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
[4292]1166         !
[6140]1167         !                                                 ! open boundaries
[7646]1168         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
[4486]1169#if defined key_agrif                                                           
1170         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
[4292]1171#endif
1172         !                                             !* Swap
1173         !                                             !  ----
[7753]1174         ubb_e  (:,:) = ub_e  (:,:)
1175         ub_e   (:,:) = un_e  (:,:)
1176         un_e   (:,:) = ua_e  (:,:)
1177         !
1178         vbb_e  (:,:) = vb_e  (:,:)
1179         vb_e   (:,:) = vn_e  (:,:)
1180         vn_e   (:,:) = va_e  (:,:)
1181         !
1182         sshbb_e(:,:) = sshb_e(:,:)
1183         sshb_e (:,:) = sshn_e(:,:)
1184         sshn_e (:,:) = ssha_e(:,:)
[4292]1185
1186         !                                             !* Sum over whole bt loop
1187         !                                             !  ----------------------
1188         za1 = wgtbtp1(jn)                                   
[6140]1189         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
[7753]1190            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
1191            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
[9023]1192         ELSE                                       ! Sum transports
1193            IF ( .NOT.ln_wd_dl ) THEN 
1194               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
1195               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
1196            ELSE
1197               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
1198               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
1199            END IF
[4292]1200         ENDIF
[9023]1201         !                                          ! Sum sea level
[7753]1202         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
[9023]1203
[358]1204         !                                                 ! ==================== !
1205      END DO                                               !        end loop      !
1206      !                                                    ! ==================== !
[1438]1207      ! -----------------------------------------------------------------------------
[1502]1208      ! Phase 3. update the general trend with the barotropic trend
[1438]1209      ! -----------------------------------------------------------------------------
[1502]1210      !
[4292]1211      ! Set advection velocity correction:
[9023]1212      IF (ln_bt_fw) THEN
1213         zwx(:,:) = un_adv(:,:)
1214         zwy(:,:) = vn_adv(:,:)
1215         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
[9043]1216            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1217            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
[9023]1218            !
1219            ! Update corrective fluxes for next time step:
1220            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1221            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1222         ELSE
1223            un_bf(:,:) = 0._wp
1224            vn_bf(:,:) = 0._wp 
1225         END IF         
1226         ! Save integrated transport for next computation
[7753]1227         ub2_b(:,:) = zwx(:,:)
1228         vb2_b(:,:) = zwy(:,:)
[4292]1229      ENDIF
[9023]1230
1231
[4292]1232      !
1233      ! Update barotropic trend:
[6140]1234      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]1235         DO jk=1,jpkm1
[9043]1236            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1237            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
[4292]1238         END DO
1239      ELSE
[5930]1240         ! At this stage, ssha has been corrected: compute new depths at velocity points
1241         DO jj = 1, jpjm1
1242            DO ji = 1, jpim1      ! NO Vector Opt.
[9554]1243               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
1244                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      &
[5930]1245                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
[9554]1246               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
1247                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
[5930]1248                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1249            END DO
1250         END DO
[10170]1251         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
[5930]1252         !
[4292]1253         DO jk=1,jpkm1
[9043]1254            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1255            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
[4292]1256         END DO
1257         ! Save barotropic velocities not transport:
[7753]1258         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1259         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
[4292]1260      ENDIF
[9023]1261
1262
1263      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
[4292]1264      DO jk = 1, jpkm1
[9023]1265         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)
1266         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
[358]1267      END DO
[9023]1268
1269      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1270         DO jk = 1, jpkm1
[9109]1271            un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1272                       & + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk) 
1273            vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) & 
1274                       & + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk) 
[9023]1275         END DO
1276      END IF
1277
1278     
1279      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current
1280      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current
[1502]1281      !
[4486]1282#if defined key_agrif
1283      ! Save time integrated fluxes during child grid integration
[5656]1284      ! (used to update coarse grid transports at next time step)
[4486]1285      !
[6140]1286      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1287         IF( Agrif_NbStepint() == 0 ) THEN
[7753]1288            ub2_i_b(:,:) = 0._wp
1289            vb2_i_b(:,:) = 0._wp
[4486]1290         END IF
1291         !
1292         za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]1293         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1294         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]1295      ENDIF
1296#endif     
[1502]1297      !                                   !* write time-spliting arrays in the restart
[6140]1298      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
[508]1299      !
[9023]1300      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
1301      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]1302      !
[9019]1303      IF( ln_diatmb ) THEN
[9554]1304         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
1305         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
[6140]1306      ENDIF
[2715]1307      !
[508]1308   END SUBROUTINE dyn_spg_ts
1309
[6140]1310
[4292]1311   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1312      !!---------------------------------------------------------------------
1313      !!                   ***  ROUTINE ts_wgt  ***
1314      !!
1315      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1316      !!----------------------------------------------------------------------
1317      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
1318      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
1319      INTEGER, INTENT(inout) :: jpit      ! cycle length   
1320      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
1321                                                         zwgt2    ! Secondary weights
1322     
1323      INTEGER ::  jic, jn, ji                      ! temporary integers
1324      REAL(wp) :: za1, za2
1325      !!----------------------------------------------------------------------
[508]1326
[4292]1327      zwgt1(:) = 0._wp
1328      zwgt2(:) = 0._wp
1329
1330      ! Set time index when averaged value is requested
1331      IF (ll_fw) THEN
1332         jic = nn_baro
1333      ELSE
1334         jic = 2 * nn_baro
1335      ENDIF
1336
1337      ! Set primary weights:
1338      IF (ll_av) THEN
1339           ! Define simple boxcar window for primary weights
1340           ! (width = nn_baro, centered around jic)     
1341         SELECT CASE ( nn_bt_flt )
1342              CASE( 0 )  ! No averaging
1343                 zwgt1(jic) = 1._wp
1344                 jpit = jic
1345
1346              CASE( 1 )  ! Boxcar, width = nn_baro
1347                 DO jn = 1, 3*nn_baro
1348                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1349                    IF (za1 < 0.5_wp) THEN
1350                      zwgt1(jn) = 1._wp
1351                      jpit = jn
1352                    ENDIF
1353                 ENDDO
1354
1355              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
1356                 DO jn = 1, 3*nn_baro
1357                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1358                    IF (za1 < 1._wp) THEN
1359                      zwgt1(jn) = 1._wp
1360                      jpit = jn
1361                    ENDIF
1362                 ENDDO
1363              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1364         END SELECT
1365
1366      ELSE ! No time averaging
1367         zwgt1(jic) = 1._wp
1368         jpit = jic
1369      ENDIF
1370   
1371      ! Set secondary weights
1372      DO jn = 1, jpit
1373        DO ji = jn, jpit
1374             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1375        END DO
1376      END DO
1377
1378      ! Normalize weigths:
1379      za1 = 1._wp / SUM(zwgt1(1:jpit))
1380      za2 = 1._wp / SUM(zwgt2(1:jpit))
1381      DO jn = 1, jpit
1382        zwgt1(jn) = zwgt1(jn) * za1
1383        zwgt2(jn) = zwgt2(jn) * za2
1384      END DO
1385      !
1386   END SUBROUTINE ts_wgt
1387
[6140]1388
[508]1389   SUBROUTINE ts_rst( kt, cdrw )
1390      !!---------------------------------------------------------------------
1391      !!                   ***  ROUTINE ts_rst  ***
1392      !!
1393      !! ** Purpose : Read or write time-splitting arrays in restart file
1394      !!----------------------------------------------------------------------
[9528]1395      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1396      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
[508]1397      !!----------------------------------------------------------------------
1398      !
[9506]1399      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1400         !                                   ! ---------------
[10345]1401         IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN    !* Read the restart file
[9506]1402            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
1403            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
1404            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
1405            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
1406            IF( .NOT.ln_bt_av ) THEN
1407               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
1408               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
1409               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
1410               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
1411               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
1412               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
1413            ENDIF
[4486]1414#if defined key_agrif
[9506]1415            ! Read time integrated fluxes
1416            IF ( .NOT.Agrif_Root() ) THEN
1417               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
1418               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
1419            ENDIF
1420#endif
1421         ELSE                                   !* Start from rest
1422            IF(lwp) WRITE(numout,*)
1423            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
1424            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
1425            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
1426            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
1427#if defined key_agrif
1428            IF ( .NOT.Agrif_Root() ) THEN
1429               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
1430            ENDIF
1431#endif
[4486]1432         ENDIF
[9506]1433         !
1434      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1435         !                                   ! -------------------
1436         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[9367]1437         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1438         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
1439         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
1440         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
1441         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
[4292]1442         !
1443         IF (.NOT.ln_bt_av) THEN
[9367]1444            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
1445            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
1446            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
1447            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
1448            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
1449            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
[4292]1450         ENDIF
[4486]1451#if defined key_agrif
1452         ! Save time integrated fluxes
1453         IF ( .NOT.Agrif_Root() ) THEN
[9367]1454            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
1455            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
[4486]1456         ENDIF
1457#endif
[9367]1458         IF( lwxios ) CALL iom_swap(      cxios_context          )
[4292]1459      ENDIF
1460      !
1461   END SUBROUTINE ts_rst
[2528]1462
[6140]1463
1464   SUBROUTINE dyn_spg_ts_init
[4292]1465      !!---------------------------------------------------------------------
1466      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1467      !!
1468      !! ** Purpose : Set time splitting options
1469      !!----------------------------------------------------------------------
[6140]1470      INTEGER  ::   ji ,jj              ! dummy loop indices
1471      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
[9019]1472      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
[4292]1473      !!----------------------------------------------------------------------
[4370]1474      !
[5930]1475      ! Max courant number for ext. grav. waves
[4370]1476      !
[5930]1477      DO jj = 1, jpj
1478         DO ji =1, jpi
1479            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1480            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
[7646]1481            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
[4370]1482         END DO
[5930]1483      END DO
1484      !
[5836]1485      zcmax = MAXVAL( zcu(:,:) )
[10297]1486      IF( lk_mpp )   CALL mpp_max( 'dynspg_ts', zcmax )
[2528]1487
[4370]1488      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[6140]1489      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
[4292]1490     
[5836]1491      rdtbt = rdt / REAL( nn_baro , wp )
[4292]1492      zcmax = zcmax * rdtbt
[9023]1493      ! Print results
[4292]1494      IF(lwp) WRITE(numout,*)
[9169]1495      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1496      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]1497      IF( ln_bt_auto ) THEN
[9169]1498         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
[4370]1499         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
[4292]1500      ELSE
[9169]1501         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
[358]1502      ENDIF
[4292]1503
1504      IF(ln_bt_av) THEN
[9169]1505         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
[4292]1506      ELSE
[9169]1507         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]1508      ENDIF
[508]1509      !
[4292]1510      !
1511      IF(ln_bt_fw) THEN
[4370]1512         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
[4292]1513      ELSE
[4370]1514         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]1515      ENDIF
1516      !
[4486]1517#if defined key_agrif
1518      ! Restrict the use of Agrif to the forward case only
[9023]1519!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]1520#endif
1521      !
[4370]1522      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]1523      SELECT CASE ( nn_bt_flt )
[6140]1524         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1525         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1526         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
[9169]1527         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]1528      END SELECT
1529      !
[4370]1530      IF(lwp) WRITE(numout,*) ' '
1531      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1532      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1533      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1534      !
[9023]1535      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1536      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1537         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1538      ENDIF
1539      !
[6140]1540      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]1541         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1542      ENDIF
[6140]1543      IF( zcmax>0.9_wp ) THEN
[4292]1544         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1545      ENDIF
1546      !
[9124]1547      !                             ! Allocate time-splitting arrays
1548      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1549      !
1550      !                             ! read restart when needed
[9506]1551      CALL ts_rst( nit000, 'READ' )
[9124]1552      !
[9367]1553      IF( lwxios ) THEN
1554! define variables in restart file when writing with XIOS
1555         CALL iom_set_rstw_var_active('ub2_b')
1556         CALL iom_set_rstw_var_active('vb2_b')
1557         CALL iom_set_rstw_var_active('un_bf')
1558         CALL iom_set_rstw_var_active('vn_bf')
1559         !
1560         IF (.NOT.ln_bt_av) THEN
1561            CALL iom_set_rstw_var_active('sshbb_e')
1562            CALL iom_set_rstw_var_active('ubb_e')
1563            CALL iom_set_rstw_var_active('vbb_e')
1564            CALL iom_set_rstw_var_active('sshb_e')
1565            CALL iom_set_rstw_var_active('ub_e')
1566            CALL iom_set_rstw_var_active('vb_e')
1567         ENDIF
1568#if defined key_agrif
1569         ! Save time integrated fluxes
1570         IF ( .NOT.Agrif_Root() ) THEN
1571            CALL iom_set_rstw_var_active('ub2_i_b')
1572            CALL iom_set_rstw_var_active('vb2_i_b')
1573         ENDIF
1574#endif
1575      ENDIF
1576      !
[4292]1577   END SUBROUTINE dyn_spg_ts_init
[508]1578
[358]1579   !!======================================================================
1580END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.