source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9261

Last change on this file since 9261 was 9261, checked in by andmirek, 3 years ago

#1978 only 1 and 2 for nn_timing

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