source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7900

Last change on this file since 7900 was 7900, checked in by glong, 5 years ago

Comment changes to dynhpg.F90, moved halo exchange to 2d method from 3d method in divcur.F90 and picked off spg and wind stress in dynspg_ts.F90

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