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/UKMO/NEMO_4.0_mirror/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/OCE/DYN/dynspg_ts.F90 @ 10888

Last change on this file since 10888 was 10888, checked in by davestorkey, 5 years ago

branches/UKMO/NEMO_4.0_mirror : clear SVN keywords

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