source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4254

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

Branch 2013/dev_r3858_NOC_ZTC, #863. Merge in final changes from the dev_r3867_MERCATOR1_DYN branch; mainly AGRIF and BDY compatibility

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