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

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2236

Last change on this file since 2236 was 2236, checked in by cetlod, 14 years ago

First guess of NEMO_v3.3

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 33.6 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
4   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
5   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
6   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
7   !!              -   ! 2008-01  (R. Benshila)  change averaging method
8   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
10  !!---------------------------------------------------------------------
11#if defined key_dynspg_ts   ||   defined key_esopa
12   !!----------------------------------------------------------------------
13   !!   'key_dynspg_ts'         free surface cst volume with time splitting
14   !!----------------------------------------------------------------------
15   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
16   !!                 splitting scheme and add to the general trend
17   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! surface boundary condition: ocean
22   USE dynspg_oce      ! surface pressure gradient variables
23   USE phycst          ! physical constants
24   USE domvvl          ! variable volume
25   USE zdfbfr          ! bottom friction
26   USE obcdta          ! open boundary condition data     
27   USE obcfla          ! Flather open boundary condition 
28   USE dynvor          ! vorticity term
29   USE obc_oce         ! Lateral open boundary condition
30   USE obc_par         ! open boundary condition parameters
31   USE bdy_oce         ! unstructured open boundaries
32   USE bdy_par         ! unstructured open boundaries
33   USE bdydta          ! unstructured open boundaries
34   USE bdydyn          ! unstructured open boundaries
35   USE bdytides        ! tidal forcing at unstructured open boundaries.
36   USE lib_mpp         ! distributed memory computing library
37   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
38   USE prtctl          ! Print control
39   USE in_out_manager  ! I/O manager
40   USE iom
41   USE restart         ! only for lrst_oce
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC dyn_spg_ts  ! routine called by step.F90
47   PUBLIC ts_rst      ! routine called by istate.F90
48
49
50   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter
51   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
52
53   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! averaged velocity
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!-------------------------------------------------------------------------
59   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
60   !! $Id$
61   !! Software is governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt)
62   !!-------------------------------------------------------------------------
63
64CONTAINS
65
66   SUBROUTINE dyn_spg_ts( kt )
67      !!----------------------------------------------------------------------
68      !!                  ***  routine dyn_spg_ts  ***
69      !!
70      !! ** Purpose :   Compute the now trend due to the surface pressure
71      !!      gradient in case of free surface formulation with time-splitting.
72      !!      Add it to the general trend of momentum equation.
73      !!
74      !! ** Method  :   Free surface formulation with time-splitting
75      !!      -1- Save the vertically integrated trend. This general trend is
76      !!          held constant over the barotropic integration.
77      !!          The Coriolis force is removed from the general trend as the
78      !!          surface gradient and the Coriolis force are updated within
79      !!          the barotropic integration.
80      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
81      !!          barotropic velocity (ua_e and va_e) through barotropic
82      !!          momentum and continuity integration. Barotropic former
83      !!          variables are time averaging over the full barotropic cycle
84      !!          (= 2 * baroclinic time step) and saved in zuX_b
85      !!          and zvX_b (X specifying after, now or before).
86      !!      -3- The new general trend becomes :
87      !!          ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) )
88      !!
89      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
90      !!
91      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
92      !!---------------------------------------------------------------------
93      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
94      !!
95      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
96      INTEGER  ::   icycle           ! temporary scalar
97      INTEGER  ::   ikbu, ikbv        ! temporary scalar
98
99      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! temporary scalars
100      REAL(wp) ::   z1_8, zx1, zy1                   !    -         -
101      REAL(wp) ::   z1_4, zx2, zy2                   !     -         -
102      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !     -         -
103      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !     -         -
104      !!
105      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e
106      !!
107      REAL(wp), DIMENSION(jpi,jpj) ::   zbfru  , zbfrv     ! 2D workspace
108      !!
109      REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace
110      !!
111      REAL(wp), DIMENSION(jpi,jpj) ::   zcu, zwx, zua, zun, zub   ! 2D workspace
112      REAL(wp), DIMENSION(jpi,jpj) ::   zcv, zwy, zva, zvn, zvb   !  -      -
113      REAL(wp), DIMENSION(jpi,jpj) ::   zun_e, zub_e, zu_sum      ! 2D workspace
114      REAL(wp), DIMENSION(jpi,jpj) ::   zvn_e, zvb_e, zv_sum      !  -      -
115      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_sum                  !  -      -
116      !!----------------------------------------------------------------------
117
118      IF( kt == nit000 ) THEN             !* initialisation
119         !
120         IF(lwp) WRITE(numout,*)
121         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
122         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
123         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
124         !
125         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b 
126         !
127         ua_e  (:,:) = un_b (:,:)
128         va_e  (:,:) = vn_b (:,:)
129         hu_e  (:,:) = hu   (:,:)
130         hv_e  (:,:) = hv   (:,:)
131         hur_e (:,:) = hur  (:,:)
132         hvr_e (:,:) = hvr  (:,:)
133         IF( ln_dynvor_een ) THEN
134            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
135            DO jj = 2, jpj
136               DO ji = fs_2, jpi   ! vector opt.
137                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
138                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
139                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
140                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
141               END DO
142            END DO
143         ENDIF
144         !
145      ENDIF
146
147      !                                   !* Local constant initialization
148      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
149      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
150      z1_4 = 0.5 * 0.5
151      zraur  = 1. / rau0                                    ! 1 / volumic mass
152      !
153      zhdiv(:,:) = 0.e0                                     ! barotropic divergence
154      zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F)
155      zv_sld = 0.e0   ;   zv_asp = 0.e0
156
157      ! -----------------------------------------------------------------------------
158      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
159      ! -----------------------------------------------------------------------------
160      !     
161      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
162      !                                   ! --------------------------
163      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   zub(:,:) = 0.e0
164      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   zvb(:,:) = 0.e0
165      !
166      DO jk = 1, jpkm1
167#if defined key_vectopt_loop
168         DO jj = 1, 1         !Vector opt. => forced unrolling
169            DO ji = 1, jpij
170#else
171         DO jj = 1, jpj
172            DO ji = 1, jpi
173#endif
174               !                                                                              ! now trend
175               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
176               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
177               !                                                                              ! now velocity
178               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
179               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
180               !                                                                              ! before velocity
181               zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 
182               zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
183            END DO
184         END DO
185      END DO
186
187      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
188      DO jk = 1, jpkm1                    ! --------------------------
189         DO jj = 2, jpjm1
190            DO ji = fs_2, fs_jpim1   ! vector opt.
191               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
192               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
193            END DO
194         END DO
195      END DO
196
197      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
198      !                                   ! ---------------------------====
199      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport
200      zwy(:,:) = zvn(:,:) * e1v(:,:)
201      !
202      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
203         DO jj = 2, jpjm1
204            DO ji = fs_2, fs_jpim1   ! vector opt.
205               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
206               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
207               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
208               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
209               ! energy conserving formulation for planetary vorticity term
210               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
211               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
212            END DO
213         END DO
214         !
215      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
219               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
220               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
221               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
222            END DO
223         END DO
224         !
225      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
226         DO jj = 2, jpjm1
227            DO ji = fs_2, fs_jpim1   ! vector opt.
228               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
229                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
230               zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
231                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
232            END DO
233         END DO
234         !
235      ENDIF
236
237      !                                   !* Right-Hand-Side of the barotropic momentum equation
238      !                                   ! ----------------------------------------------------
239      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
240         DO jj = 2, jpjm1 
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
243                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
244               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
245                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
246            END DO
247         END DO
248      ENDIF
249
250      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
251         DO ji = fs_2, fs_jpim1
252            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
253            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
254         END DO
255      END DO
256
257                   
258      !                                             ! Remove barotropic contribution of bottom friction
259      !                                             ! from the barotropic transport trend
260      zcoef = -1. / z2dt_b
261# if defined key_vectopt_loop
262      DO jj = 1, 1
263         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
264# else
265      DO jj = 2, jpjm1
266         DO ji = 2, jpim1
267# endif
268            ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
269            ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
270            !
271            ! Apply stability criteria for bottom friction
272            !RBbug for vvl and external mode we may need to
273            ! use varying fse3
274            zbfru  (ji,jj) = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zcoef )
275            zbfrv  (ji,jj) = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zcoef )
276         END DO
277      END DO
278
279      IF( lk_vvl ) THEN
280         DO jj = 2, jpjm1
281            DO ji = fs_2, fs_jpim1   ! vector opt.
282               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj)   &
283                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
284               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj)   &
285                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
286            END DO
287         END DO
288      ELSE
289         DO jj = 2, jpjm1
290            DO ji = fs_2, fs_jpim1   ! vector opt.
291               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj) * hur(ji,jj)
292               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj) * hvr(ji,jj)
293            END DO
294         END DO
295      ENDIF
296
297      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
298      !                                   ! --------------------------
299      zua(:,:) = zua(:,:) * hur(:,:)
300      zva(:,:) = zva(:,:) * hvr(:,:)
301      !
302      IF( lk_vvl ) THEN
303         zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
304         zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
305      ELSE
306         zub(:,:) = zub(:,:) * hur(:,:)
307         zvb(:,:) = zvb(:,:) * hvr(:,:)
308      ENDIF
309
310      ! -----------------------------------------------------------------------
311      !  Phase 2 : Integration of the barotropic equations with time splitting
312      ! -----------------------------------------------------------------------
313      !
314      !                                             ! ==================== !
315      !                                             !    Initialisations   !
316      !                                             ! ==================== !
317      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
318     
319      !                                ! Start from NOW field
320      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
321      hv_e   (:,:) = hv   (:,:) 
322      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
323      hvr_e  (:,:) = hvr  (:,:)
324!RBbug     zsshb_e(:,:) = sshn (:,:) 
325      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
326      sshn_e (:,:) = sshn (:,:)
327     
328      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
329      zvn_e  (:,:) = vn_b (:,:)
330      zub_e  (:,:) = un_b (:,:)
331      zvb_e  (:,:) = vn_b (:,:)
332
333      zu_sum  (:,:) = un_b (:,:)           ! summation
334      zv_sum  (:,:) = vn_b (:,:)
335      zssh_sum(:,:) = sshn (:,:)
336
337#if defined key_obc
338      ! set ssh corrections to 0
339      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
340      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
341      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
342      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
343      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
344#endif
345
346      !                                             ! ==================== !
347      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
348         !                                          ! ==================== !
349         z2dt_e = 2. * ( rdt / nn_baro )
350         IF( jn == 1 )   z2dt_e = rdt / nn_baro
351
352         !                                                !* Update the forcing (OBC, BDY and tides)
353         !                                                !  ------------------
354         IF( lk_obc                     )   CALL obc_dta_bt( kt, jn   )
355         IF( lk_bdy  .OR.  ln_bdy_tides )   CALL bdy_dta_bt( kt, jn+1, icycle )
356
357         !                                                !* after ssh_e
358         !                                                !  -----------
359         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
360            DO ji = fs_2, fs_jpim1   ! vector opt.
361               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
362                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
363                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
364                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
365            END DO
366         END DO
367         !
368#if defined key_obc
369         !                                                     ! OBC : zhdiv must be zero behind the open boundary
370!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
371         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
372         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
373         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
374         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
375#endif
376#if defined key_bdy
377         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
378#endif
379         !
380         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
381            DO ji = fs_2, fs_jpim1   ! vector opt.
382               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
383            END DO
384         END DO
385
386         !                                                !* after barotropic velocities (vorticity scheme dependent)
387         !                                                !  --------------------------- 
388         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
389         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
390         !
391         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
392            DO jj = 2, jpjm1
393               DO ji = fs_2, fs_jpim1   ! vector opt.
394                  ! surface pressure gradient
395                  IF( lk_vvl) THEN
396                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
397                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
398                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
399                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
400                  ELSE
401                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
402                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
403                  ENDIF
404                  ! energy conserving formulation for planetary vorticity term
405                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
406                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
407                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
408                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
409                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
410                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
411                  ! after velocities with implicit bottom friction
412                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   &
413                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
414                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   &
415                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
416               END DO
417            END DO
418            !
419         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
420            DO jj = 2, jpjm1
421               DO ji = fs_2, fs_jpim1   ! vector opt.
422                   ! surface pressure gradient
423                  IF( lk_vvl) THEN
424                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
425                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
426                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
427                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
428                  ELSE
429                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
430                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
431                  ENDIF
432                  ! enstrophy conserving formulation for planetary vorticity term
433                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
434                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
435                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
436                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
437                  ! after velocities with implicit bottom friction
438                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   &
439                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
440                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   &
441                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
442               END DO
443            END DO
444            !
445         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
446            DO jj = 2, jpjm1
447               DO ji = fs_2, fs_jpim1   ! vector opt.
448                  ! surface pressure gradient
449                  IF( lk_vvl) THEN
450                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
451                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
452                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
453                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
454                  ELSE
455                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
456                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
457                  ENDIF
458                  ! energy/enstrophy conserving formulation for planetary vorticity term
459                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
460                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
461                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
462                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
463                  ! after velocities with implicit bottom friction
464                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   &
465                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
466                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   &
467                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
468               END DO
469            END DO
470            !
471         ENDIF
472         !                                                !* domain lateral boundary
473         !                                                !  -----------------------
474         !                                                      ! Flather's boundary condition for the barotropic loop :
475         !                                                      !         - Update sea surface height on each open boundary
476         !                                                      !         - Correct the velocity
477
478         IF( lk_obc                   )   CALL obc_fla_ts
479         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla( sshn_e ) 
480         !
481         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
482         CALL lbc_lnk( va_e  , 'V', -1. )
483         CALL lbc_lnk( ssha_e, 'T',  1. )
484
485         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
486         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
487         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
488
489         !                                                !* Time filter and swap
490         !                                                !  --------------------
491         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
492            zsshb_e(:,:) = sshn_e(:,:)
493            zub_e  (:,:) = zun_e (:,:)
494            zvb_e  (:,:) = zvn_e (:,:)
495            sshn_e (:,:) = ssha_e(:,:)
496            zun_e  (:,:) = ua_e  (:,:)
497            zvn_e  (:,:) = va_e  (:,:)
498         ELSE                                                   ! Swap + Filter
499            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
500            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
501            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
502            sshn_e (:,:) = ssha_e(:,:)
503            zun_e  (:,:) = ua_e  (:,:)
504            zvn_e  (:,:) = va_e  (:,:)
505         ENDIF
506
507         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
508            !                                             !  ------------------
509            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
510               DO ji = 1, fs_jpim1   ! Vector opt.
511                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
512                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
513                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
514                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
515                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
516                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
517               END DO
518            END DO
519            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
520            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
521            !
522            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
523            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
524            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
525            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
526            !
527         ENDIF
528         !                                                 ! ==================== !
529      END DO                                               !        end loop      !
530      !                                                    ! ==================== !
531
532#if defined key_obc
533      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
534      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
535      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
536      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
537#endif
538
539      ! -----------------------------------------------------------------------------
540      ! Phase 3. update the general trend with the barotropic trend
541      ! -----------------------------------------------------------------------------
542      !
543      !                                   !* Time average ==> after barotropic u, v, ssh
544      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
545      un_b  (:,:) = zcoef * zu_sum  (:,:) 
546      vn_b  (:,:) = zcoef * zv_sum  (:,:) 
547      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
548      !
549      !                                   !* update the general momentum trend
550      DO jk=1,jpkm1
551         ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b
552         va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b
553      END DO
554      !
555      !                                   !* write time-spliting arrays in the restart
556      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
557      !
558      !
559   END SUBROUTINE dyn_spg_ts
560
561
562   SUBROUTINE ts_rst( kt, cdrw )
563      !!---------------------------------------------------------------------
564      !!                   ***  ROUTINE ts_rst  ***
565      !!
566      !! ** Purpose : Read or write time-splitting arrays in restart file
567      !!----------------------------------------------------------------------
568      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
569      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
570      !
571      INTEGER ::  ji, jk        ! dummy loop indices
572      !!----------------------------------------------------------------------
573      !
574      IF( TRIM(cdrw) == 'READ' ) THEN
575         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
576            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
577            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
578         ELSE
579            un_b (:,:) = 0.e0
580            vn_b (:,:) = 0.e0
581            ! vertical sum
582            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
583               DO jk = 1, jpkm1
584                  DO ji = 1, jpij
585                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
586                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
587                  END DO
588               END DO
589            ELSE                             ! No  vector opt.
590               DO jk = 1, jpkm1
591                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
592                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
593               END DO
594            ENDIF
595            un_b (:,:) = un_b(:,:) * hur(:,:)
596            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
597         ENDIF
598         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
599            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered extrenal ssh
600         ELSE
601            sshn_b(:,:)=sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
602         ENDIF
603      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
604         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
605         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
606         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !
607      ENDIF
608      !
609   END SUBROUTINE ts_rst
610
611#else
612   !!----------------------------------------------------------------------
613   !!   Default case :   Empty module   No standart free surface cst volume
614   !!----------------------------------------------------------------------
615CONTAINS
616   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
617      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
618   END SUBROUTINE dyn_spg_ts
619   SUBROUTINE ts_rst( kt, cdrw )     ! Empty routine
620      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
621      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
622      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
623   END SUBROUTINE ts_rst   
624#endif
625   
626   !!======================================================================
627END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.