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_jki.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynspg_ts_jki.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.1 KB
Line 
1MODULE dynspg_ts_jki
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_ts_jki  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if defined key_dynspg_ts   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_ts'                    free surface with time splitting
9   !!----------------------------------------------------------------------
10   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
11   !!                 splitting scheme and add to the general trend
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE phycst          ! physical constants
17   USE ocesbc          ! ocean surface boundary condition
18   USE obcdta          ! open boundary condition data     
19   USE obcfla          ! Flather open boundary condition 
20   USE dynvor          ! vorticity term
21   USE obc_oce         ! Lateral open boundary condition
22   USE obc_par         ! open boundary condition parameters
23   USE lib_mpp         ! distributed memory computing library
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE prtctl          ! Print control
26   USE dynspg_oce      ! surface pressure gradient variables
27   USE dynspg_ts       ! surface pressure gradient
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O library
30   USE restart         ! only for lrst_oce
31
32   IMPLICIT NONE
33   PRIVATE
34
35   !! * Accessibility
36   PUBLIC dyn_spg_ts_jki  ! routine called by step.F90
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !!   OPA 9.0 , LODYC-IPSL  (2005)
43   !! $Header$
44   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE dyn_spg_ts_jki( kt )
50      !!----------------------------------------------------------------------
51      !!                  ***  routine dyn_spg_ts_jki  ***
52      !!
53      !! ** Purpose :   Compute the now trend due to the surface pressure
54      !!      gradient in case of free surface formulation with time-splitting.
55      !!      Add it to the general trend of momentum equation.
56      !!      Compute the free surface.
57      !!
58      !! ** Method  :   Free surface formulation with time-splitting
59      !!      -1- Save the vertically integrated trend. This general trend is
60      !!          held constant over the barotropic integration.
61      !!          The Coriolis force is removed from the general trend as the
62      !!          surface gradient and the Coriolis force are updated within
63      !!          the barotropic integration.
64      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
65      !!          barotropic transports (ua_e and va_e) through barotropic
66      !!          momentum and continuity integration. Barotropic former
67      !!          variables are time averaging over the full barotropic cycle
68      !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b
69      !!          and zvX_b (X specifying after, now or before).
70      !!      -3- Update of sea surface height from time averaged barotropic
71      !!          variables.
72      !!        - apply lateral boundary conditions on sshn.
73      !!      -4- The new general trend becomes :
74      !!          ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H
75      !!
76      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
77      !!
78      !! References :
79      !!   Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
80      !!
81      !! History :
82      !!   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code
83      !!        !  05-11  (V. Garnier, G. Madec)  optimization
84      !!---------------------------------------------------------------------
85      !! * Arguments
86      INTEGER, INTENT( in )  ::   kt           ! ocean time-step index
87
88      !! * Local declarations
89      INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices
90      INTEGER  ::  icycle                      ! temporary scalar
91      REAL(wp) ::                           &
92         zraur, zcoef, z2dt_e, z2dt_b, zfac25,   &  ! temporary scalars
93         zfact1, zspgu, zcubt, zx1, zy1,    &  !     "        "
94         zfact2, zspgv, zcvbt, zx2, zy2        !     "        "
95      REAL(wp), DIMENSION(jpi,jpj) ::       &
96         zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays
97         zua, zva, zub, zvb,                &  !     "        "
98         zssha_b, zua_b, zva_b,             &  !     "        "
99         zsshb_e, zub_e, zvb_e,             &  !     "        "
100         zun_e, zvn_e                          !     "        "
101      REAL(wp), DIMENSION(jpi,jpj),SAVE ::  &
102         ztnw, ztne, ztsw, ztse
103      !!----------------------------------------------------------------------
104
105      ! Arrays initialization
106      ! ---------------------
107      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0
108      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0
109      zhdiv(:,:) = 0.e0
110
111      IF( kt == nit000 ) THEN
112
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_spg_ts_jki : surface pressure gradient trend'
115         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   free surface with time splitting with j-k-i loop'
116         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2*rdt/rdtbt )
117
118         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields:
119         !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
120
121         ssha_e(:,:) = sshn(:,:)
122         ua_e  (:,:) = un_b(:,:)
123         va_e  (:,:) = vn_b(:,:)
124
125         IF( ln_dynvor_een ) THEN
126            ztne(1,:) = 0.e0   ;   ztnw(1,:) = 0.e0   ;   ztse(1,:) = 0.e0   ;   ztsw(1,:) = 0.e0
127            DO jj = 2, jpj
128               DO ji = 2, jpi
129                  ztne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
130                  ztnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
131                  ztse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
132                  ztsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
133               END DO
134            END DO
135         ENDIF
136
137      ENDIF
138   
139      ! Local constant initialization
140      ! --------------------------------
141      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
142      IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
143      zfact1 = 0.5 * 0.25                                   ! coefficient for vorticity estimates
144      zfact2 = 0.5 * 0.5
145      zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water
146     
147      ! -----------------------------------------------------------------------------
148      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
149      ! -----------------------------------------------------------------------------
150
151      DO jj = 1, jpj
152
153         ! variables for the barotropic equations
154         zsshb_e(:,jj) = sshn_b(:,jj)       ! (barotropic) sea surface height (before and now)
155         sshn_e (:,jj) = sshn_b(:,jj)
156         zub_e  (:,jj) = un_b  (:,jj)       ! barotropic transports issued from the barotropic equations (before and now)
157         zvb_e  (:,jj) = vn_b  (:,jj)
158         zun_e  (:,jj) = un_b  (:,jj)
159         zvn_e  (:,jj) = vn_b  (:,jj)
160         zssha_b(:,jj) = sshn  (:,jj)        ! time averaged variables over all sub-timesteps
161         zua_b  (:,jj) = un_b  (:,jj)   
162         zva_b  (:,jj) = vn_b  (:,jj)
163
164         ! Vertically integrated quantities
165         ! --------------------------------
166         zua(:,jj) = 0.e0
167         zva(:,jj) = 0.e0
168         zub(:,jj) = 0.e0
169         zvb(:,jj) = 0.e0
170         zwx(:,jj) = 0.e0
171         zwy(:,jj) = 0.e0
172
173         ! vertical sum
174         DO jk = 1, jpkm1
175            !                                                           ! Vertically integrated momentum trends
176            zua(:,jj) = zua(:,jj) + fse3u(:,jj,jk) * umask(:,jj,jk) * ua(:,jj,jk)
177            zva(:,jj) = zva(:,jj) + fse3v(:,jj,jk) * vmask(:,jj,jk) * va(:,jj,jk)
178            !                                                           ! Vertically integrated transports (before)
179            zub(:,jj) = zub(:,jj) + fse3u(:,jj,jk) * ub(:,jj,jk)
180            zvb(:,jj) = zvb(:,jj) + fse3v(:,jj,jk) * vb(:,jj,jk)
181            !                                                           ! Planetary vorticity (now)
182            zwx(:,jj) = zwx(:,jj) + e2u(:,jj) * fse3u(:,jj,jk) * un(:,jj,jk)
183            zwy(:,jj) = zwy(:,jj) + e1v(:,jj) * fse3v(:,jj,jk) * vn(:,jj,jk)
184         END DO
185
186      END DO
187
188      DO jj = 2, jpjm1
189
190         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
191            DO ji = 2, jpim1
192               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
193               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
194               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
195               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
196               ! energy conserving formulation for planetary vorticity term
197               zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
198               zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
199            END DO
200
201         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
202            DO ji = 2, jpim1
203               zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
204                              + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
205               zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
206                              + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
207               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
208               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
209            END DO
210
211         ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
212         zfac25 = 0.25
213            DO ji = 2, jpim1
214               zcu(ji,jj) = + zfac25 / e1u(ji,jj)   &
215                  &       * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
216                  &          + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
217               zcv(ji,jj) = - zfac25 / e2v(ji,jj)   &
218                  &       * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
219                  &          + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
220            END DO
221
222         ENDIF
223
224
225         ! Remove barotropic trend from general momentum trend
226         DO jk = 1 , jpkm1
227            DO ji = 2, jpim1
228               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
229               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
230            END DO
231         END DO
232
233         ! Remove coriolis term from barotropic trend
234         ! ------------------------------------------
235         DO ji = 2, jpim1
236            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
237            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
238         END DO
239
240      END DO
241
242      ! -----------------------------------------------------------------------
243      !  Phase 2 : Integration of the barotropic equations with time splitting
244      ! -----------------------------------------------------------------------
245
246      ! Initialisations
247      !----------------
248      ! Number of iteration of the barotropic loop
249      icycle = FLOOR( z2dt_b / rdtbt )
250
251      ! set ssh corrections to 0
252      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
253#if defined key_obc
254      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
255      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
256      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
257      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
258#endif
259
260      ! Barotropic integration over 2 baroclinic time steps
261      ! ---------------------------------------------------
262
263      !                                                    ! ==================== !
264      DO jit = 1, icycle                                   !  sub-time-step loop  !
265         !                                                 ! ==================== !
266
267         z2dt_e = 2. * rdtbt
268         IF ( jit == 1 )   z2dt_e = rdtbt
269
270         ! Time interpolation of open boundary condition data
271         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
272
273         DO jj = 2, jpjm1
274
275         ! Horizontal divergence of barotropic transports
276         !--------------------------------------------------
277            DO ji = 2, jpim1
278               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              &
279                  &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              &
280                  &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              &
281                  &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          &
282                  &           / (e1t(ji,jj)*e2t(ji,jj))
283            END DO
284
285#if defined key_obc
286         ! open boundaries (div must be zero behind the open boundary)
287         !  mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column
288         IF( lp_obc_east  ) THEN
289            IF( nje0   <= jj .AND. jj <= nje1   )   zhdiv(nie0p1:nie1p1,jj) = 0.e0      ! east
290         ENDIF
291         IF( lp_obc_west  ) THEN
292            IF( njw0   <= jj .AND. jj <= njw1   )   zhdiv(niw0  :niw1  ,jj) = 0.e0      ! west
293         ENDIF
294         IF( lp_obc_north ) THEN
295            IF( njn0p1 <= jj .AND. jj <= njn1p1 )   zhdiv(nin0  :nin1  ,jj) = 0.e0      ! north
296         ENDIF
297         IF( lp_obc_south ) THEN
298            IF( njs0   <= jj .AND. jj <= njs1   )   zhdiv(nis0  :nis1  ,jj) = 0.e0      ! south
299         ENDIF
300#endif
301
302         ! Sea surface height from the barotropic system
303         !----------------------------------------------
304            DO ji = 2, jpim1
305               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
306                  &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
307            END DO
308
309         END DO
310
311         ! evolution of the barotropic transport ( following the vorticity scheme used)
312         ! ----------------------------------------------------------------------------
313         zwx(:,:) = e2u(:,:) * zun_e(:,:)
314         zwy(:,:) = e1v(:,:) * zvn_e(:,:)
315
316         DO jj = 2, jpjm1
317
318            IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
319               DO ji = 2, jpim1
320                  ! surface pressure gradient
321                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
322                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
323                  ! energy conserving formulation for planetary vorticity term
324                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
325                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
326                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
327                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
328                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
329                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
330                  ! after transports
331                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
332                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
333               END DO
334
335            ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
336               DO ji = 2, jpim1
337                  ! surface pressure gradient
338                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
339                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
340                  ! enstrophy conserving formulation for planetary vorticity term
341                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
342                                 + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
343                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
344                                 + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
345                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
346                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
347                  ! after transports
348                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
349                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
350               END DO
351
352            ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
353               zfac25 = 0.25
354               DO ji = 2, jpim1
355                  ! surface pressure gradient
356                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
357                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
358                  ! energy/enstrophy conserving formulation for planetary vorticity term
359                  zcubt = + zfac25 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
360                     &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
361                  zcvbt = - zfac25 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
362                     &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
363                  ! after transports
364                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
365                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
366               END DO
367            ENDIF
368
369         END DO
370
371
372         ! Flather's boundary condition for the barotropic loop :
373         !         - Update sea surface height on each open boundary
374         !         - Correct the barotropic transports
375         IF( lk_obc )   CALL obc_fla_ts
376
377         ! ... Boundary conditions on ua_e, va_e, ssha_e
378         CALL lbc_lnk( ua_e  , 'U', -1. )
379         CALL lbc_lnk( va_e  , 'V', -1. )
380         CALL lbc_lnk( ssha_e, 'T',  1. )
381
382         DO jj = 1, jpj
383
384            ! temporal sum
385            !-------------
386            zssha_b(:,jj) = zssha_b(:,jj) + ssha_e(:,jj)
387            zua_b  (:,jj) = zua_b  (:,jj) + ua_e  (:,jj)
388            zva_b  (:,jj) = zva_b  (:,jj) + va_e  (:,jj) 
389
390            ! Time filter and swap of dynamics arrays
391            ! ---------------------------------------
392            IF( neuler == 0 .AND. jit == 1 ) THEN       ! Euler (forward) time stepping
393               zsshb_e(:,jj) = sshn_e(:,jj)
394               zub_e  (:,jj) = zun_e (:,jj)
395               zvb_e  (:,jj) = zvn_e (:,jj)
396               sshn_e (:,jj) = ssha_e(:,jj)
397               zun_e  (:,jj) = ua_e  (:,jj)
398               zvn_e  (:,jj) = va_e  (:,jj)
399            ELSE                                        ! Asselin filtering
400               zsshb_e(:,jj) = atfp * ( zsshb_e(:,jj) + ssha_e(:,jj) ) + atfp1 * sshn_e (:,jj)
401               zub_e  (:,jj) = atfp * ( zub_e  (:,jj) + ua_e  (:,jj) ) + atfp1 * zun_e  (:,jj)
402               zvb_e  (:,jj) = atfp * ( zvb_e  (:,jj) + va_e  (:,jj) ) + atfp1 * zvn_e  (:,jj)
403               sshn_e (:,jj) = ssha_e(:,jj)
404               zun_e  (:,jj) = ua_e  (:,jj)
405               zvn_e  (:,jj) = va_e  (:,jj)
406            ENDIF
407
408         END DO
409
410         !                                                 ! ==================== !
411      END DO                                               !        end loop      !
412      !                                                    ! ==================== !
413
414
415      ! Time average of after barotropic variables
416      zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
417      zssha_b(:,:) = zcoef * zssha_b(:,:) 
418      zua_b  (:,:) = zcoef * zua_b  (:,:) 
419      zva_b  (:,:) = zcoef * zva_b  (:,:) 
420#if defined key_obc
421         IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
422         IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
423         IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
424         IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
425#endif
426     
427
428      ! ---------------------------------------------------------------------------
429      ! Phase 3 : Update sea surface height from time averaged barotropic variables
430      ! ---------------------------------------------------------------------------
431
432      sshb(:,:) = sshn(:,:)
433 
434      ! Horizontal divergence of time averaged barotropic transports
435      !-------------------------------------------------------------
436      DO jj = 2, jpjm1
437         DO ji = 2, jpim1
438            zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
439           &                +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
440           &             / ( e1t(ji,jj) * e2t(ji,jj) )
441         END DO
442
443#if defined key_obc
444         ! open boundaries (div must be zero behind the open boundary)
445         !  mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column
446         IF( lp_obc_east  ) THEN
447            IF( nje0   <= jj .AND. jj <= nje1   )   zhdiv(nie0p1:nie1p1,jj) = 0.e0      ! east
448         ENDIF
449         IF( lp_obc_west  ) THEN
450            IF( njw0   <= jj .AND. jj <= njw1   )   zhdiv(niw0  :niw1  ,jj) = 0.e0      ! west
451         ENDIF
452         IF( lp_obc_north ) THEN
453            IF( njn0p1 <= jj .AND. jj <= njn1p1 )   zhdiv(nin0  :nin1  ,jj) = 0.e0      ! north
454         ENDIF
455         IF( lp_obc_south ) THEN
456            IF( njs0   <= jj .AND. jj <= njs1   )   zhdiv(nis0  :nis1  ,jj) = 0.e0      ! south
457         ENDIF
458#endif
459
460         ! sea surface height
461         !-------------------
462         DO ji = 2, jpim1
463            sshn(ji,jj) = (  sshb_b(ji,jj) - z2dt_b * ( zraur * emp(ji,jj) + zhdiv(ji,jj) )  ) * tmask(ji,jj,1)
464         END DO
465
466      END DO
467
468      ! ... Boundary conditions on sshn
469      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
470
471
472      ! -----------------------------------------------------------------------------
473      ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
474      ! -----------------------------------------------------------------------------
475
476      DO jj = 1, jpj
477
478         ! Swap on time averaged barotropic variables
479         ! ------------------------------------------
480         sshb_b(:,jj) = sshn_b (:,jj)
481         sshn_b(:,jj) = zssha_b(:,jj)
482         un_b  (:,jj) = zua_b  (:,jj) 
483         vn_b  (:,jj) = zva_b  (:,jj) 
484   
485         ! add time averaged barotropic coriolis and surface pressure gradient
486         ! terms to the general momentum trend
487         ! --------------------------------------------------------------------
488         DO jk = 1, jpkm1
489            ua(:,jj,jk) = ua(:,jj,jk) + hur(:,jj) * ( zua_b(:,jj) - zub(:,jj) ) / z2dt_b
490            va(:,jj,jk) = va(:,jj,jk) + hvr(:,jj) * ( zva_b(:,jj) - zvb(:,jj) ) / z2dt_b
491         END DO
492      END DO
493
494      ! write filtered free surface arrays in restart file
495      ! --------------------------------------------------
496      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
497
498      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
499         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
500      ENDIF
501     
502   END SUBROUTINE dyn_spg_ts_jki
503#else
504   !!----------------------------------------------------------------------
505   !!   Default case :   Empty module   No standart free surface cst volume
506   !!----------------------------------------------------------------------
507CONTAINS
508   SUBROUTINE dyn_spg_ts_jki( kt )       ! Empty routine
509      WRITE(*,*) 'dyn_spg_ts_jki: You should not have seen this print! error?', kt
510   END SUBROUTINE dyn_spg_ts_jki
511#endif
512   
513   !!======================================================================
514END MODULE dynspg_ts_jki
Note: See TracBrowser for help on using the repository browser.