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

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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

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