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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization

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