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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago

#1613: vvl by default, step III: Merge with the trunk (free surface simplification) (see wiki)

  • Property svn:keywords set to Id
File size: 54.9 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_ts  ***
4   !! Ocean dynamics:  surface pressure gradient trend, split-explicit scheme
5   !!======================================================================
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
12   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
13   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
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
16   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
17   !!---------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
21   !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
22   !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
23   !!   ts_rst         : read/write time-splitting fields in restart file
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE sbc_oce         ! surface boundary condition: ocean
28   USE zdf_oce         ! Bottom friction coefts
29   USE sbcisf          ! ice shelf variable (fwfisf)
30   USE sbcapr          ! surface boundary condition: atmospheric pressure
31   USE dynadv    , ONLY: ln_dynadv_vec
32   USE phycst          ! physical constants
33   USE dynvor          ! vorticity term
34   USE bdy_par         ! for lk_bdy
35   USE bdytides        ! open boundary condition data
36   USE bdydyn2d        ! open boundary conditions on barotropic variables
37   USE sbctide         ! tides
38   USE updtide         ! tide potential
39   !
40   USE in_out_manager  ! I/O manager
41   USE lib_mpp         ! distributed memory computing library
42   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
43   USE prtctl          ! Print control
44   USE iom             ! IOM library
45   USE restart         ! only for lrst_oce
46   USE wrk_nemo        ! Memory Allocation
47   USE timing          ! Timing   
48#if defined key_agrif
49   USE agrif_opa_interp ! agrif
50#endif
51#if defined key_asminc   
52   USE asminc          ! Assimilation increment
53#endif
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90
59   PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
60   PUBLIC dyn_spg_ts_init   !    "      "     "    "
61   PUBLIC ts_rst            !    "      "     "    "
62
63   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
64   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
65
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   wgtbtp1, wgtbtp2   !: 1st & 2nd weights used in time filtering of barotropic fields
67
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff/h at F points
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   !: triad of coriolis parameter
70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   !: (only used with een vorticity scheme)
71
72   !! Time filtered arrays at baroclinic time step:
73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv     !: Advection vel. at "now" barocl. step
74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b  , vb2_b      !: Half step fluxes (ln_bt_fw=T)
75#if defined key_agrif
76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_i_b,  vb2_i_b   !: Half step time integrated fluxes
77#endif
78
79   !! Arrays at barotropic time step:                   ! bef before !   before   !    now     !   after    !
80   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ubb_e    ,    ub_e    ,    un_e    ,    ua_e    !: u-external velocity
81   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vbb_e    ,    vb_e    ,    vn_e    ,    va_e    !: v-external velocity
82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshbb_e  ,    sshb_e  ,    sshn_e  ,    ssha_e  !: external ssh
83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hu_e                 !: external u-depth
84   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hv_e                 !: external v-depth
85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hur_e                !: inverse of u-depth
86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hvr_e                !: inverse of v-depth
87
88   !! * Substitutions
89#  include "vectopt_loop_substitute.h90"
90   !!----------------------------------------------------------------------
91   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
92   !! $Id$
93   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
94   !!----------------------------------------------------------------------
95CONTAINS
96
97   INTEGER FUNCTION dyn_spg_ts_alloc()
98      !!----------------------------------------------------------------------
99      !!                  ***  routine dyn_spg_ts_alloc  ***
100      !!----------------------------------------------------------------------
101      INTEGER :: ierr(5)
102      !!----------------------------------------------------------------------
103      ierr(:) = 0
104      !
105      ALLOCATE( ssha_e(jpi,jpj),  sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
106         &        ua_e(jpi,jpj),    un_e(jpi,jpj),   ub_e(jpi,jpj),   ubb_e(jpi,jpj), &
107         &        va_e(jpi,jpj),    vn_e(jpi,jpj),   vb_e(jpi,jpj),   vbb_e(jpi,jpj), &
108         &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT=ierr(1) )
109         !
110      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj)                , STAT=ierr(2) )
111      !
112      IF( ln_dynvor_een )   ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj)                 , & 
113         &                            ftsw(jpi,jpj) , ftse(jpi,jpj)                 , STAT=ierr(3) )
114         !
115      ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj)    , STAT=ierr(4) )
116#if defined key_agrif
117      ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                                  , STAT= ierr(5))
118#endif
119      !
120      dyn_spg_ts_alloc = MAXVAL( ierr(:) )
121      !
122      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
123      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
124      !
125   END FUNCTION dyn_spg_ts_alloc
126
127
128   SUBROUTINE dyn_spg_ts( kt )
129      !!----------------------------------------------------------------------
130      !!
131      !! ** Purpose : - Compute the now trend due to the explicit time stepping
132      !!              of the quasi-linear barotropic system, and add it to the
133      !!              general momentum trend.
134      !!
135      !! ** Method  : - split-explicit schem (time splitting) :
136      !!      Barotropic variables are advanced from internal time steps
137      !!      "n"   to "n+1" if ln_bt_fw=T
138      !!      or from
139      !!      "n-1" to "n+1" if ln_bt_fw=F
140      !!      thanks to a generalized forward-backward time stepping (see ref. below).
141      !!
142      !! ** Action :
143      !!      -Update the filtered free surface at step "n+1"      : ssha
144      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
145      !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
146      !!      These are used to advect tracers and are compliant with discrete
147      !!      continuity equation taken at the baroclinic time steps. This
148      !!      ensures tracers conservation.
149      !!      - (ua, va) momentum trend updated with barotropic component.
150      !!
151      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
152      !!---------------------------------------------------------------------
153      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
154      !
155      LOGICAL  ::   ll_fw_start      ! if true, forward integration
156      LOGICAL  ::   ll_init          ! if true, special startup of 2d equations
157      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
158      INTEGER  ::   ikbu, ikbv, noffset        ! local integers
159      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf   ! local scalars
160      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      -
161      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      -
162      REAL(wp) ::   zu_spg, zv_spg             !   -      -
163      REAL(wp) ::   zhura, zhvra               !   -      -
164      REAL(wp) ::   za0, za1, za2, za3         !   -      -
165      !
166      REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e
167      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
168      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv
169      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
170      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
171      REAL(wp), POINTER, DIMENSION(:,:) :: zhf
172      !!----------------------------------------------------------------------
173      !
174      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
175      !
176      !                                         !* Allocate temporary arrays
177      CALL wrk_alloc( jpi,jpj,   zsshp2_e, zhdiv )
178      CALL wrk_alloc( jpi,jpj,   zu_trd, zv_trd)
179      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc)
180      CALL wrk_alloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e)
181      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                                   )
182      CALL wrk_alloc( jpi,jpj,   zhf )
183      !
184      z1_12 = 1._wp / 12._wp                 !* Local constant initialization
185      z1_8  = 0.125_wp                                   
186      z1_4  = 0.25_wp
187      z1_2  = 0.5_wp     
188      zraur = 1._wp / rau0
189      !                                            ! reciprocal of baroclinic time step
190      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
191      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt
192      ENDIF
193      z1_2dt_b = 1.0_wp / z2dt_bf 
194      !
195      ll_init     = ln_bt_av                       ! if no time averaging, then no specific restart
196      ll_fw_start = .FALSE.
197      !                                            ! time offset in steps for bdy data update
198      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro
199      ELSE                       ;   noffset =   0 
200      ENDIF
201      !
202      IF( kt == nit000 ) THEN                !* initialisation
203         !
204         IF(lwp) WRITE(numout,*)
205         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
206         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
207         IF(lwp) WRITE(numout,*)
208         !
209         IF( neuler == 0 )   ll_init=.TRUE.
210         !
211         IF( ln_bt_fw .OR. neuler == 0 ) THEN
212            ll_fw_start =.TRUE.
213            noffset     = 0
214         ELSE
215            ll_fw_start =.FALSE.
216         ENDIF
217         !
218         ! Set averaging weights and cycle length:
219         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
220         !
221      ENDIF
222      !
223      ! Set arrays to remove/compute coriolis trend.
224      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
225      ! Note that these arrays are also used during barotropic loop. These are however frozen
226      ! although they should be updated in the variable volume case. Not a big approximation.
227      ! To remove this approximation, copy lines below inside barotropic loop
228      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
229      !
230      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
231         IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==!
232            SELECT CASE( nn_een_e3f )              !* ff/e3 at F-point
233            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
234               DO jj = 1, jpjm1
235                  DO ji = 1, jpim1
236                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
237                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp 
238                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
239                  END DO
240               END DO
241            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
242               DO jj = 1, jpjm1
243                  DO ji = 1, jpim1
244                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     &
245                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   &
246                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
247                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
248                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
249                  END DO
250               END DO
251            END SELECT
252            CALL lbc_lnk( zwz, 'F', 1._wp )
253            !
254            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
255            DO jj = 2, jpj
256               DO ji = 2, jpi
257                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
258                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
259                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
260                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
261               END DO
262            END DO
263            !
264         ELSE                                !== all other schemes (ENE, ENS, MIX)
265            zwz(:,:) = 0._wp
266            zhf(:,:) = 0._wp
267            IF ( .not. ln_sco ) THEN
268
269!!gm  agree the JC comment  : this should be done in a much clear way
270
271! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
272!     Set it to zero for the time being
273!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
274!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
275!              ENDIF
276!              zhf(:,:) = gdepw_0(:,:,jk+1)
277            ELSE
278               zhf(:,:) = hbatf(:,:)
279            END IF
280
281            DO jj = 1, jpjm1
282               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
283            END DO
284
285            DO jk = 1, jpkm1
286               DO jj = 1, jpjm1
287                  zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
288               END DO
289            END DO
290            CALL lbc_lnk( zhf, 'F', 1._wp )
291            ! JC: TBC. hf should be greater than 0
292            DO jj = 1, jpj
293               DO ji = 1, jpi
294                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
295               END DO
296            END DO
297            zwz(:,:) = ff(:,:) * zwz(:,:)
298         ENDIF
299      ENDIF
300      !
301      ! If forward start at previous time step, and centered integration,
302      ! then update averaging weights:
303      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
304         ll_fw_start=.FALSE.
305         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
306      ENDIF
307                         
308      ! -----------------------------------------------------------------------------
309      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
310      ! -----------------------------------------------------------------------------
311      !     
312      !
313      !                                   !* e3*d/dt(Ua) (Vertically integrated)
314      !                                   ! --------------------------------------------------
315      zu_frc(:,:) = 0._wp
316      zv_frc(:,:) = 0._wp
317      !
318      DO jk = 1, jpkm1
319         zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
320         zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
321      END DO
322      !
323      zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
324      zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
325      !
326      !
327      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
328      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
329         DO jj = 2, jpjm1
330            DO ji = fs_2, fs_jpim1   ! vector opt.
331               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
332               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
333            END DO
334         END DO
335      END DO
336      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
337      !                                   ! --------------------------------------------------------
338      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes
339      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
340      !
341      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
342         DO jj = 2, jpjm1
343            DO ji = fs_2, fs_jpim1   ! vector opt.
344               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
345               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
346               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
347               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
348               ! energy conserving formulation for planetary vorticity term
349               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
350               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
351            END DO
352         END DO
353         !
354      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
355         DO jj = 2, jpjm1
356            DO ji = fs_2, fs_jpim1   ! vector opt.
357               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
358                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
359               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
360                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
361               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
362               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
363            END DO
364         END DO
365         !
366      ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme
367         DO jj = 2, jpjm1
368            DO ji = fs_2, fs_jpim1   ! vector opt.
369               zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
370                &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
371                &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
372                &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
373               zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
374                &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
375                &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
376                &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
377            END DO
378         END DO
379         !
380      ENDIF 
381      !
382      !                                   !* Right-Hand-Side of the barotropic momentum equation
383      !                                   ! ----------------------------------------------------
384      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
385         DO jj = 2, jpjm1 
386            DO ji = fs_2, fs_jpim1   ! vector opt.
387               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
388               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)
389            END DO
390         END DO
391      ENDIF
392
393      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
394         DO ji = fs_2, fs_jpim1
395             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
396             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
397          END DO
398      END DO 
399      !
400      !                 ! Add bottom stress contribution from baroclinic velocities:     
401      IF (ln_bt_fw) THEN
402         DO jj = 2, jpjm1                         
403            DO ji = fs_2, fs_jpim1   ! vector opt.
404               ikbu = mbku(ji,jj)       
405               ikbv = mbkv(ji,jj)   
406               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
407               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
408            END DO
409         END DO
410      ELSE
411         DO jj = 2, jpjm1
412            DO ji = fs_2, fs_jpim1   ! vector opt.
413               ikbu = mbku(ji,jj)       
414               ikbv = mbkv(ji,jj)   
415               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
416               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
417            END DO
418         END DO
419      ENDIF
420      !
421      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
422      zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:)
423      zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:)
424      !       
425      IF( ln_bt_fw ) THEN                        ! Add wind forcing
426         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:)
427         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:)
428      ELSE
429         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:)
430         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:)
431      ENDIF 
432      !
433      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
434         IF (ln_bt_fw) THEN
435            DO jj = 2, jpjm1             
436               DO ji = fs_2, fs_jpim1   ! vector opt.
437                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
438                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
439                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
440                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
441               END DO
442            END DO
443         ELSE
444            DO jj = 2, jpjm1             
445               DO ji = fs_2, fs_jpim1   ! vector opt.
446                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
447                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
448                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
449                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
450                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
451                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
452               END DO
453            END DO
454         ENDIF
455      ENDIF
456      !                                   !* Right-Hand-Side of the barotropic ssh equation
457      !                                   ! -----------------------------------------------
458      !                                         ! Surface net water flux and rivers
459      IF (ln_bt_fw) THEN
460         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
461      ELSE
462         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
463                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
464      ENDIF
465#if defined key_asminc
466      !                                         ! Include the IAU weighted SSH increment
467      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
468         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
469      ENDIF
470#endif
471      !                                   !* Fill boundary data arrays for AGRIF
472      !                                   ! ------------------------------------
473#if defined key_agrif
474         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
475#endif
476      !
477      ! -----------------------------------------------------------------------
478      !  Phase 2 : Integration of the barotropic equations
479      ! -----------------------------------------------------------------------
480      !
481      !                                             ! ==================== !
482      !                                             !    Initialisations   !
483      !                                             ! ==================== ! 
484      ! Initialize barotropic variables:     
485      IF( ll_init )THEN
486         sshbb_e(:,:) = 0._wp
487         ubb_e  (:,:) = 0._wp
488         vbb_e  (:,:) = 0._wp
489         sshb_e (:,:) = 0._wp
490         ub_e   (:,:) = 0._wp
491         vb_e   (:,:) = 0._wp
492      ENDIF
493      !
494      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
495         sshn_e(:,:) =    sshn(:,:)           
496         un_e  (:,:) =    un_b(:,:)           
497         vn_e  (:,:) =    vn_b(:,:)
498         !
499         hu_e  (:,:) =    hu_n(:,:)       
500         hv_e  (:,:) =    hv_n(:,:) 
501         hur_e (:,:) = r1_hu_n(:,:)   
502         hvr_e (:,:) = r1_hv_n(:,:)
503      ELSE                                ! CENTRED integration: start from BEFORE fields
504         sshn_e(:,:) =    sshb(:,:)
505         un_e  (:,:) =    ub_b(:,:)         
506         vn_e  (:,:) =    vb_b(:,:)
507         !
508         hu_e  (:,:) =    hu_b(:,:)       
509         hv_e  (:,:) =    hv_b(:,:) 
510         hur_e (:,:) = r1_hu_b(:,:)   
511         hvr_e (:,:) = r1_hv_b(:,:)
512      ENDIF
513      !
514      !
515      !
516      ! Initialize sums:
517      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
518      va_b  (:,:) = 0._wp
519      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
520      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
521      vn_adv(:,:) = 0._wp
522      !                                             ! ==================== !
523      DO jn = 1, icycle                             !  sub-time-step loop  !
524         !                                          ! ==================== !
525         !                                                !* Update the forcing (BDY and tides)
526         !                                                !  ------------------
527         ! Update only tidal forcing at open boundaries
528#if defined key_tide
529         IF( lk_bdy      .AND. lk_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
530         IF( ln_tide_pot .AND. lk_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
531#endif
532         !
533         ! Set extrapolation coefficients for predictor step:
534         IF ((jn<3).AND.ll_init) THEN      ! Forward           
535           za1 = 1._wp                                         
536           za2 = 0._wp                       
537           za3 = 0._wp                       
538         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
539           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
540           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
541           za3 =  0.281105_wp              ! za3 = bet
542         ENDIF
543
544         ! Extrapolate barotropic velocities at step jit+0.5:
545         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
546         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
547
548         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
549            !                                             !  ------------------
550            ! Extrapolate Sea Level at step jit+0.5:
551            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
552            !
553            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
554               DO ji = 2, fs_jpim1   ! Vector opt.
555                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)     &
556                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
557                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
558                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)     &
559                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
560                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
561               END DO
562            END DO
563            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
564            !
565            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
566            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
567         ELSE
568            zhup2_e (:,:) = hu_n(:,:)
569            zhvp2_e (:,:) = hv_n(:,:)
570         ENDIF
571         !                                                !* after ssh
572         !                                                !  -----------
573         ! One should enforce volume conservation at open boundaries here
574         ! considering fluxes below:
575         !
576         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
577         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
578         !
579#if defined key_agrif
580         ! Set fluxes during predictor step to ensure volume conservation
581         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
582            IF((nbondi == -1).OR.(nbondi == 2)) THEN
583               DO jj=1,jpj
584                  zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
585               END DO
586            ENDIF
587            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
588               DO jj=1,jpj
589                  zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
590               END DO
591            ENDIF
592            IF((nbondj == -1).OR.(nbondj == 2)) THEN
593               DO ji=1,jpi
594                  zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
595               END DO
596            ENDIF
597            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
598               DO ji=1,jpi
599                  zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
600               END DO
601            ENDIF
602         ENDIF
603#endif
604         !
605         ! Sum over sub-time-steps to compute advective velocities
606         za2 = wgtbtp2(jn)
607         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
608         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
609         !
610         ! Set next sea level:
611         DO jj = 2, jpjm1                                 
612            DO ji = fs_2, fs_jpim1   ! vector opt.
613               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
614                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
615            END DO
616         END DO
617         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
618         CALL lbc_lnk( ssha_e, 'T',  1._wp )
619
620#if defined key_bdy
621         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
622         IF( lk_bdy )   CALL bdy_ssh( ssha_e )
623#endif
624#if defined key_agrif
625         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
626#endif
627         
628         ! Sea Surface Height at u-,v-points (vvl case only)
629         IF( .NOT.ln_linssh ) THEN                               
630            DO jj = 2, jpjm1
631               DO ji = 2, jpim1      ! NO Vector Opt.
632                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)  &
633                     &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
634                     &              +   e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
635                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)  &
636                     &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
637                     &              +   e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
638               END DO
639            END DO
640            CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
641         ENDIF   
642         !                                 
643         ! Half-step back interpolation of SSH for surface pressure computation:
644         !----------------------------------------------------------------------
645         IF ((jn==1).AND.ll_init) THEN
646           za0=1._wp                        ! Forward-backward
647           za1=0._wp                           
648           za2=0._wp
649           za3=0._wp
650         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
651           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
652           za1=-0.1666666666666_wp          ! za1 = gam
653           za2= 0.0833333333333_wp          ! za2 = eps
654           za3= 0._wp             
655         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
656           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
657           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
658           za2=0.088_wp                     ! za2 = gam
659           za3=0.013_wp                     ! za3 = eps
660         ENDIF
661         !
662         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
663          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
664         !
665         ! Compute associated depths at U and V points:
666         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN   !* Vector form
667            !                                       
668            DO jj = 2, jpjm1                           
669               DO ji = 2, jpim1
670                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e1e2u(ji  ,jj)    &
671                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
672                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
673                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e1e2v(ji  ,jj  )  &
674                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
675                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
676                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
677                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
678               END DO
679            END DO
680         ENDIF
681         !
682         ! Add Coriolis trend:
683         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
684         ! at each time step. We however keep them constant here for optimization.
685         ! Recall that zwx and zwy arrays hold fluxes at this stage:
686         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
687         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
688         !
689         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN     !==  energy conserving or mixed scheme  ==!
690            DO jj = 2, jpjm1
691               DO ji = fs_2, fs_jpim1   ! vector opt.
692                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
693                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
694                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
695                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
696                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
697                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
698               END DO
699            END DO
700            !
701         ELSEIF ( ln_dynvor_ens ) THEN                   !==  enstrophy conserving scheme  ==!
702            DO jj = 2, jpjm1
703               DO ji = fs_2, fs_jpim1   ! vector opt.
704                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
705                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
706                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
707                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
708                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
709                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
710               END DO
711            END DO
712            !
713         ELSEIF ( ln_dynvor_een ) THEN                   !==  energy and enstrophy conserving scheme  ==!
714            DO jj = 2, jpjm1
715               DO ji = fs_2, fs_jpim1   ! vector opt.
716                  zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
717                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
718                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
719                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
720                  zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
721                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
722                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
723                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
724               END DO
725            END DO
726            !
727         ENDIF
728         !
729         ! Add tidal astronomical forcing if defined
730         IF ( lk_tide.AND.ln_tide_pot ) THEN
731            DO jj = 2, jpjm1
732               DO ji = fs_2, fs_jpim1   ! vector opt.
733                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
734                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
735                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
736                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
737               END DO
738            END DO
739         ENDIF
740         !
741         ! Add bottom stresses:
742         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:)
743         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
744         !
745         ! Surface pressure trend:
746         DO jj = 2, jpjm1
747            DO ji = fs_2, fs_jpim1   ! vector opt.
748               ! Add surface pressure gradient
749               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
750               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
751               zwx(ji,jj) = zu_spg
752               zwy(ji,jj) = zv_spg
753            END DO
754         END DO
755         !
756         ! Set next velocities:
757         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   !* Vector form
758            DO jj = 2, jpjm1
759               DO ji = fs_2, fs_jpim1   ! vector opt.
760                  ua_e(ji,jj) = (                                 un_e(ji,jj)    & 
761                     &            + rdtbt * (                      zwx(ji,jj)    &
762                     &                                        + zu_trd(ji,jj)    &
763                     &                                        + zu_frc(ji,jj) )  ) * umask(ji,jj,1)
764                     !
765                  va_e(ji,jj) = (                                 vn_e(ji,jj)    &
766                     &            + rdtbt * (                      zwy(ji,jj)    &
767                     &                                        + zv_trd(ji,jj)    &
768                     &                                        + zv_frc(ji,jj) )  ) * vmask(ji,jj,1)
769               END DO
770            END DO
771            !
772         ELSE                                      !* Flux form
773            DO jj = 2, jpjm1
774               DO ji = fs_2, fs_jpim1   ! vector opt.
775                  zhura = umask(ji,jj,1) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1) )
776                  zhvra = vmask(ji,jj,1) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1) )
777                  !
778                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)    & 
779                     &            + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)    & 
780                     &                      + zhup2_e(ji,jj)  * zu_trd(ji,jj)    &
781                     &                      +    hu_n(ji,jj)  * zu_frc(ji,jj) )  ) * zhura
782                     !
783                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)    &
784                     &            + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)    &
785                     &                      + zhvp2_e(ji,jj)  * zv_trd(ji,jj)    &
786                     &                      +    hv_n(ji,jj)  * zv_frc(ji,jj) )  ) * zhvra
787               END DO
788            END DO
789         ENDIF
790         !
791         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
792            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
793            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
794            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
795            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
796            !
797         ENDIF
798         !                                             !* domain lateral boundary
799         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
800         !
801#if defined key_bdy 
802         !                                                 ! open boundaries
803         IF( lk_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
804#endif
805#if defined key_agrif                                                           
806         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
807#endif
808         !                                             !* Swap
809         !                                             !  ----
810         ubb_e  (:,:) = ub_e  (:,:)
811         ub_e   (:,:) = un_e  (:,:)
812         un_e   (:,:) = ua_e  (:,:)
813         !
814         vbb_e  (:,:) = vb_e  (:,:)
815         vb_e   (:,:) = vn_e  (:,:)
816         vn_e   (:,:) = va_e  (:,:)
817         !
818         sshbb_e(:,:) = sshb_e(:,:)
819         sshb_e (:,:) = sshn_e(:,:)
820         sshn_e (:,:) = ssha_e(:,:)
821
822         !                                             !* Sum over whole bt loop
823         !                                             !  ----------------------
824         za1 = wgtbtp1(jn)                                   
825         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
826            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
827            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
828         ELSE                                ! Sum transports
829            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
830            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
831         ENDIF
832         !                                   ! Sum sea level
833         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
834         !                                                 ! ==================== !
835      END DO                                               !        end loop      !
836      !                                                    ! ==================== !
837      ! -----------------------------------------------------------------------------
838      ! Phase 3. update the general trend with the barotropic trend
839      ! -----------------------------------------------------------------------------
840      !
841      ! Set advection velocity correction:
842      zwx(:,:) = un_adv(:,:)
843      zwy(:,:) = vn_adv(:,:)
844      IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN     
845         un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
846         vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
847      ELSE
848         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
849         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
850      END IF
851
852      IF( ln_bt_fw ) THEN ! Save integrated transport for next computation
853         ub2_b(:,:) = zwx(:,:)
854         vb2_b(:,:) = zwy(:,:)
855      ENDIF
856      !
857      ! Update barotropic trend:
858      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
859         DO jk=1,jpkm1
860            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
861            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
862         END DO
863      ELSE
864         ! At this stage, ssha has been corrected: compute new depths at velocity points
865         DO jj = 1, jpjm1
866            DO ji = 1, jpim1      ! NO Vector Opt.
867               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
868                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
869                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
870               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
871                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
872                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
873            END DO
874         END DO
875         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
876         !
877         DO jk=1,jpkm1
878            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
879            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
880         END DO
881         ! Save barotropic velocities not transport:
882         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
883         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
884      ENDIF
885      !
886      DO jk = 1, jpkm1
887         ! Correct velocities:
888         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
889         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
890         !
891      END DO
892      !
893      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
894      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic i-current
895      !
896#if defined key_agrif
897      ! Save time integrated fluxes during child grid integration
898      ! (used to update coarse grid transports at next time step)
899      !
900      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
901         IF( Agrif_NbStepint() == 0 ) THEN
902            ub2_i_b(:,:) = 0._wp
903            vb2_i_b(:,:) = 0._wp
904         END IF
905         !
906         za1 = 1._wp / REAL(Agrif_rhot(), wp)
907         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
908         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
909      ENDIF
910#endif     
911      !                                   !* write time-spliting arrays in the restart
912      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
913      !
914      CALL wrk_dealloc( jpi,jpj,   zsshp2_e, zhdiv )
915      CALL wrk_dealloc( jpi,jpj,   zu_trd, zv_trd )
916      CALL wrk_dealloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc )
917      CALL wrk_dealloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e )
918      CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   )
919      CALL wrk_dealloc( jpi,jpj,   zhf )
920      !
921      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
922      !
923   END SUBROUTINE dyn_spg_ts
924
925
926   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
927      !!---------------------------------------------------------------------
928      !!                   ***  ROUTINE ts_wgt  ***
929      !!
930      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
931      !!----------------------------------------------------------------------
932      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
933      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
934      INTEGER, INTENT(inout) :: jpit      ! cycle length   
935      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
936                                                         zwgt2    ! Secondary weights
937     
938      INTEGER ::  jic, jn, ji                      ! temporary integers
939      REAL(wp) :: za1, za2
940      !!----------------------------------------------------------------------
941
942      zwgt1(:) = 0._wp
943      zwgt2(:) = 0._wp
944
945      ! Set time index when averaged value is requested
946      IF (ll_fw) THEN
947         jic = nn_baro
948      ELSE
949         jic = 2 * nn_baro
950      ENDIF
951
952      ! Set primary weights:
953      IF (ll_av) THEN
954           ! Define simple boxcar window for primary weights
955           ! (width = nn_baro, centered around jic)     
956         SELECT CASE ( nn_bt_flt )
957              CASE( 0 )  ! No averaging
958                 zwgt1(jic) = 1._wp
959                 jpit = jic
960
961              CASE( 1 )  ! Boxcar, width = nn_baro
962                 DO jn = 1, 3*nn_baro
963                    za1 = ABS(float(jn-jic))/float(nn_baro) 
964                    IF (za1 < 0.5_wp) THEN
965                      zwgt1(jn) = 1._wp
966                      jpit = jn
967                    ENDIF
968                 ENDDO
969
970              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
971                 DO jn = 1, 3*nn_baro
972                    za1 = ABS(float(jn-jic))/float(nn_baro) 
973                    IF (za1 < 1._wp) THEN
974                      zwgt1(jn) = 1._wp
975                      jpit = jn
976                    ENDIF
977                 ENDDO
978              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
979         END SELECT
980
981      ELSE ! No time averaging
982         zwgt1(jic) = 1._wp
983         jpit = jic
984      ENDIF
985   
986      ! Set secondary weights
987      DO jn = 1, jpit
988        DO ji = jn, jpit
989             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
990        END DO
991      END DO
992
993      ! Normalize weigths:
994      za1 = 1._wp / SUM(zwgt1(1:jpit))
995      za2 = 1._wp / SUM(zwgt2(1:jpit))
996      DO jn = 1, jpit
997        zwgt1(jn) = zwgt1(jn) * za1
998        zwgt2(jn) = zwgt2(jn) * za2
999      END DO
1000      !
1001   END SUBROUTINE ts_wgt
1002
1003
1004   SUBROUTINE ts_rst( kt, cdrw )
1005      !!---------------------------------------------------------------------
1006      !!                   ***  ROUTINE ts_rst  ***
1007      !!
1008      !! ** Purpose : Read or write time-splitting arrays in restart file
1009      !!----------------------------------------------------------------------
1010      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1011      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1012      !
1013      !!----------------------------------------------------------------------
1014      !
1015      IF( TRIM(cdrw) == 'READ' ) THEN
1016         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1017         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1018         IF( .NOT.ln_bt_av ) THEN
1019            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1020            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1021            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1022            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1023            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1024            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1025         ENDIF
1026#if defined key_agrif
1027         ! Read time integrated fluxes
1028         IF ( .NOT.Agrif_Root() ) THEN
1029            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
1030            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
1031         ENDIF
1032#endif
1033      !
1034      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1035         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1036         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1037         !
1038         IF (.NOT.ln_bt_av) THEN
1039            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1040            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1041            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1042            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1043            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1044            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1045         ENDIF
1046#if defined key_agrif
1047         ! Save time integrated fluxes
1048         IF ( .NOT.Agrif_Root() ) THEN
1049            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
1050            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
1051         ENDIF
1052#endif
1053      ENDIF
1054      !
1055   END SUBROUTINE ts_rst
1056
1057
1058   SUBROUTINE dyn_spg_ts_init
1059      !!---------------------------------------------------------------------
1060      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1061      !!
1062      !! ** Purpose : Set time splitting options
1063      !!----------------------------------------------------------------------
1064      INTEGER  ::   ji ,jj              ! dummy loop indices
1065      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1066      REAL(wp), POINTER, DIMENSION(:,:) ::   zcu
1067      !!----------------------------------------------------------------------
1068      !
1069      ! Max courant number for ext. grav. waves
1070      !
1071      CALL wrk_alloc( jpi,jpj,   zcu )
1072      !
1073      DO jj = 1, jpj
1074         DO ji =1, jpi
1075            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1076            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1077            zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) )
1078         END DO
1079      END DO
1080      !
1081      zcmax = MAXVAL( zcu(:,:) )
1082      IF( lk_mpp )   CALL mpp_max( zcmax )
1083
1084      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1085      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1086     
1087      rdtbt = rdt / REAL( nn_baro , wp )
1088      zcmax = zcmax * rdtbt
1089                     ! Print results
1090      IF(lwp) WRITE(numout,*)
1091      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1092      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1093      IF( ln_bt_auto ) THEN
1094         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.true. Automatically set nn_baro '
1095         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1096      ELSE
1097         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist '
1098      ENDIF
1099
1100      IF(ln_bt_av) THEN
1101         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1102      ELSE
1103         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
1104      ENDIF
1105      !
1106      !
1107      IF(ln_bt_fw) THEN
1108         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1109      ELSE
1110         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1111      ENDIF
1112      !
1113#if defined key_agrif
1114      ! Restrict the use of Agrif to the forward case only
1115      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1116#endif
1117      !
1118      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1119      SELECT CASE ( nn_bt_flt )
1120         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1121         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1122         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1123         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1124      END SELECT
1125      !
1126      IF(lwp) WRITE(numout,*) ' '
1127      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1128      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1129      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1130      !
1131      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1132         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1133      ENDIF
1134      IF( zcmax>0.9_wp ) THEN
1135         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1136      ENDIF
1137      !
1138      CALL wrk_dealloc( jpi,jpj,   zcu )
1139      !
1140   END SUBROUTINE dyn_spg_ts_init
1141
1142   !!======================================================================
1143END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.