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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2621

Last change on this file since 2621 was 2618, checked in by gm, 13 years ago

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