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

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_sit/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12381

Last change on this file since 12381 was 12381, checked in by dcarneir, 2 years ago

Updating GO6 branch with G06 version in the apps trunk

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