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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6049

Last change on this file since 6049 was 6049, checked in by acc, 8 years ago

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