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 branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/dynspg_ts.F90 @ 8403

Last change on this file since 8403 was 8403, checked in by deazer, 7 years ago

Add in ROMS WAD option ln_rwd+changes for implicit Bed Friction for ln_wd option
Note no ramp placed on ROMS bed friction yet
CS15mini case added as a Test CASE
at this revision AMM15 with Pure sigma coords barotorpic runs for 4 days without failure
in with ROMS option with 20cm min deoth and 50 vertical levels
Both run for CS15mini
In real domains nothing done on reference level yet so real domains
must have not negative depth points yet.
But a basic test has been done in WAD channel test cases (WAD7)

No changes in Main line source yet. See the MY_SRC sub dir of CS15 and TEST_CASES/WAD
for actual code changes.

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