New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4687

Last change on this file since 4687 was 4687, checked in by acc, 10 years ago

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