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

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5232

Last change on this file since 5232 was 5232, checked in by davestorkey, 9 years ago

Svn keywords deactivated using "svn propdel" in
branch 2015/dev_r5021_UKMO1_CICE_coupling.

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