source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5901

Last change on this file since 5901 was 5901, checked in by jamesharle, 6 years ago

merging branch with head of the trunk

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