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

source: tags/nemo_v3_2_2/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7661

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

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