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 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 36.5 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
4   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
5   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
6   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
7   !!              -   ! 2008-01  (R. Benshila)  change averaging method
8   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
10  !!---------------------------------------------------------------------
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 dynspg.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_in_use, wrk_not_released
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(wrk_in_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         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b 
147         !
148         ua_e  (:,:) = un_b (:,:)
149         va_e  (:,:) = vn_b (:,:)
150         hu_e  (:,:) = hu   (:,:)
151         hv_e  (:,:) = hv   (:,:)
152         hur_e (:,:) = hur  (:,:)
153         hvr_e (:,:) = hvr  (:,:)
154         IF( ln_dynvor_een ) THEN
155            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
156            DO jj = 2, jpj
157               DO ji = fs_2, jpi   ! vector opt.
158                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
159                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
160                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
161                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
162               END DO
163            END DO
164         ENDIF
165         !
166      ENDIF
167
168      !                                   !* Local constant initialization
169      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
170      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
171      z1_4 = 0.5 * 0.5
172      zraur  = 1. / rau0                                    ! 1 / volumic mass
173      !
174      zhdiv(:,:) = 0.e0                                     ! barotropic divergence
175      zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F)
176      zv_sld = 0.e0   ;   zv_asp = 0.e0
177
178      ! -----------------------------------------------------------------------------
179      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
180      ! -----------------------------------------------------------------------------
181      !     
182      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
183      !                                   ! --------------------------
184      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0 
185      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0 
186      !
187      DO jk = 1, jpkm1
188#if defined key_vectopt_loop
189         DO jj = 1, 1         !Vector opt. => forced unrolling
190            DO ji = 1, jpij
191#else
192         DO jj = 1, jpj
193            DO ji = 1, jpi
194#endif
195               !                                                                              ! now trend
196               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
197               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
198               !                                                                              ! now velocity
199               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
200               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
201            END DO
202         END DO
203      END DO
204
205      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
206      DO jk = 1, jpkm1                    ! --------------------------
207         DO jj = 2, jpjm1
208            DO ji = fs_2, fs_jpim1   ! vector opt.
209               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
210               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
211            END DO
212         END DO
213      END DO
214
215      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
216      !                                   ! ---------------------------====
217      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport
218      zwy(:,:) = zvn(:,:) * e1v(:,:)
219      !
220      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
221         DO jj = 2, jpjm1
222            DO ji = fs_2, fs_jpim1   ! vector opt.
223               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
224               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
225               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
226               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
227               ! energy conserving formulation for planetary vorticity term
228               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
229               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
230            END DO
231         END DO
232         !
233      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
234         DO jj = 2, jpjm1
235            DO ji = fs_2, fs_jpim1   ! vector opt.
236               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
237               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
238               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
239               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
240            END DO
241         END DO
242         !
243      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
244         DO jj = 2, jpjm1
245            DO ji = fs_2, fs_jpim1   ! vector opt.
246               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
247                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
248               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)   &
249                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
250            END DO
251         END DO
252         !
253      ENDIF
254
255      !                                   !* Right-Hand-Side of the barotropic momentum equation
256      !                                   ! ----------------------------------------------------
257      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
258         DO jj = 2, jpjm1 
259            DO ji = fs_2, fs_jpim1   ! vector opt.
260               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
261                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
262               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
263                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
264            END DO
265         END DO
266      ENDIF
267
268      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
269         DO ji = fs_2, fs_jpim1
270            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
271            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
272         END DO
273      END DO
274
275                   
276      !                                             ! Remove barotropic contribution of bottom friction
277      !                                             ! from the barotropic transport trend
278      zcoef = -1. / z2dt_b
279# if defined key_vectopt_loop
280      DO jj = 1, 1
281         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
282# else
283      DO jj = 2, jpjm1
284         DO ji = 2, jpim1
285# endif
286            ! Apply stability criteria for bottom friction
287            !RBbug for vvl and external mode we may need to use varying fse3
288            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
289            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
290            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
291         END DO
292      END DO
293
294      IF( lk_vvl ) THEN
295         DO jj = 2, jpjm1
296            DO ji = fs_2, fs_jpim1   ! vector opt.
297               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
298                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
299               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
300                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
301            END DO
302         END DO
303      ELSE
304         DO jj = 2, jpjm1
305            DO ji = fs_2, fs_jpim1   ! vector opt.
306               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
307               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
308            END DO
309         END DO
310      ENDIF
311
312      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
313      !                                   ! --------------------------
314      zua(:,:) = zua(:,:) * hur(:,:)
315      zva(:,:) = zva(:,:) * hvr(:,:)
316      !
317
318      ! -----------------------------------------------------------------------
319      !  Phase 2 : Integration of the barotropic equations with time splitting
320      ! -----------------------------------------------------------------------
321      !
322      !                                             ! ==================== !
323      !                                             !    Initialisations   !
324      !                                             ! ==================== !
325      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
326     
327      !                                ! Start from NOW field
328      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
329      hv_e   (:,:) = hv   (:,:) 
330      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
331      hvr_e  (:,:) = hvr  (:,:)
332!RBbug     zsshb_e(:,:) = sshn (:,:) 
333      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
334      sshn_e (:,:) = sshn (:,:)
335     
336      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
337      zvn_e  (:,:) = vn_b (:,:)
338      zub_e  (:,:) = un_b (:,:)
339      zvb_e  (:,:) = vn_b (:,:)
340
341      zu_sum  (:,:) = un_b (:,:)           ! summation
342      zv_sum  (:,:) = vn_b (:,:)
343      zssh_sum(:,:) = sshn (:,:)
344
345#if defined key_obc
346      ! set ssh corrections to 0
347      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
348      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
349      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
350      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
351      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
352#endif
353
354      !                                             ! ==================== !
355      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
356         !                                          ! ==================== !
357         z2dt_e = 2. * ( rdt / nn_baro )
358         IF( jn == 1 )   z2dt_e = rdt / nn_baro
359
360         !                                                !* Update the forcing (OBC, BDY and tides)
361         !                                                !  ------------------
362         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
363         IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle )
364
365         !                                                !* after ssh_e
366         !                                                !  -----------
367         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
368            DO ji = fs_2, fs_jpim1   ! vector opt.
369               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
370                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
371                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
372                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
373            END DO
374         END DO
375         !
376#if defined key_obc
377         !                                                     ! OBC : zhdiv must be zero behind the open boundary
378!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
379         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
380         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
381         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
382         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
383#endif
384#if defined key_bdy
385         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
386#endif
387         !
388         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
389            DO ji = fs_2, fs_jpim1   ! vector opt.
390               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
391            END DO
392         END DO
393
394         !                                                !* after barotropic velocities (vorticity scheme dependent)
395         !                                                !  --------------------------- 
396         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
397         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
398         !
399         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
400            DO jj = 2, jpjm1
401               DO ji = fs_2, fs_jpim1   ! vector opt.
402                  ! surface pressure gradient
403                  IF( lk_vvl) THEN
404                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
405                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
406                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
407                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
408                  ELSE
409                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
410                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
411                  ENDIF
412                  ! energy conserving formulation for planetary vorticity term
413                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
414                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
415                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
416                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
417                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
418                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
419                  ! after velocities with implicit bottom friction
420                  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)   &
421                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
422                  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)   &
423                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
424               END DO
425            END DO
426            !
427         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
428            DO jj = 2, jpjm1
429               DO ji = fs_2, fs_jpim1   ! vector opt.
430                   ! surface pressure gradient
431                  IF( lk_vvl) THEN
432                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
433                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
434                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
435                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
436                  ELSE
437                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
438                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
439                  ENDIF
440                  ! enstrophy conserving formulation for planetary vorticity term
441                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
442                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
443                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
444                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
445                  ! after velocities with implicit bottom friction
446                  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)   &
447                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
448                  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)   &
449                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
450               END DO
451            END DO
452            !
453         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
454            DO jj = 2, jpjm1
455               DO ji = fs_2, fs_jpim1   ! vector opt.
456                  ! surface pressure gradient
457                  IF( lk_vvl) THEN
458                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
459                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
460                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
461                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
462                  ELSE
463                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
464                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
465                  ENDIF
466                  ! energy/enstrophy conserving formulation for planetary vorticity term
467                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
468                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
469                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
470                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
471                  ! after velocities with implicit bottom friction
472                  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)   &
473                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
474                  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)   &
475                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
476               END DO
477            END DO
478            !
479         ENDIF
480         !                                                !* domain lateral boundary
481         !                                                !  -----------------------
482         !                                                      ! Flather's boundary condition for the barotropic loop :
483         !                                                      !         - Update sea surface height on each open boundary
484         !                                                      !         - Correct the velocity
485
486         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
487         IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e ) 
488         !
489         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
490         CALL lbc_lnk( va_e  , 'V', -1. )
491         CALL lbc_lnk( ssha_e, 'T',  1. )
492
493         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
494         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
495         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
496
497         !                                                !* Time filter and swap
498         !                                                !  --------------------
499         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
500            zsshb_e(:,:) = sshn_e(:,:)
501            zub_e  (:,:) = zun_e (:,:)
502            zvb_e  (:,:) = zvn_e (:,:)
503            sshn_e (:,:) = ssha_e(:,:)
504            zun_e  (:,:) = ua_e  (:,:)
505            zvn_e  (:,:) = va_e  (:,:)
506         ELSE                                                   ! Swap + Filter
507            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
508            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
509            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
510            sshn_e (:,:) = ssha_e(:,:)
511            zun_e  (:,:) = ua_e  (:,:)
512            zvn_e  (:,:) = va_e  (:,:)
513         ENDIF
514
515         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
516            !                                             !  ------------------
517            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
518               DO ji = 1, fs_jpim1   ! Vector opt.
519                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
520                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
521                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
522                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
523                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
524                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
525               END DO
526            END DO
527            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
528            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
529            !
530            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
531            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
532            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
533            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
534            !
535         ENDIF
536         !                                                 ! ==================== !
537      END DO                                               !        end loop      !
538      !                                                    ! ==================== !
539
540#if defined key_obc
541      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
542      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
543      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
544      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
545#endif
546
547      ! -----------------------------------------------------------------------------
548      ! Phase 3. update the general trend with the barotropic trend
549      ! -----------------------------------------------------------------------------
550      !
551      !                                   !* Time average ==> after barotropic u, v, ssh
552      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
553      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
554      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
555      !
556      !                                   !* update the general momentum trend
557      DO jk=1,jpkm1
558         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
559         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
560      END DO
561      ub_b (:,:) = un_b(:,:)
562      vb_b (:,:) = vn_b(:,:)
563      un_b  (:,:) =  zu_sum(:,:) 
564      vn_b  (:,:) =  zv_sum(:,:) 
565      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
566      !
567      !                                   !* write time-spliting arrays in the restart
568      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
569      !
570      !
571      IF(wrk_not_released(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,         &
572                              11,12,13,14,15,16,17,18,19,20,21))THEN
573         CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays.')
574      END IF
575      !
576   END SUBROUTINE dyn_spg_ts
577
578
579   SUBROUTINE ts_rst( kt, cdrw )
580      !!---------------------------------------------------------------------
581      !!                   ***  ROUTINE ts_rst  ***
582      !!
583      !! ** Purpose : Read or write time-splitting arrays in restart file
584      !!----------------------------------------------------------------------
585      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
586      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
587      !
588      INTEGER ::  ji, jk        ! dummy loop indices
589      !!----------------------------------------------------------------------
590      !
591      IF( TRIM(cdrw) == 'READ' ) THEN
592         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
593            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
594            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
595         ELSE
596            un_b (:,:) = 0.e0
597            vn_b (:,:) = 0.e0
598            ! vertical sum
599            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
600               DO jk = 1, jpkm1
601                  DO ji = 1, jpij
602                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
603                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
604                  END DO
605               END DO
606            ELSE                             ! No  vector opt.
607               DO jk = 1, jpkm1
608                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
609                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
610               END DO
611            ENDIF
612            un_b (:,:) = un_b(:,:) * hur(:,:)
613            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
614         ENDIF
615
616         ! Vertically integrated velocity (before)
617         IF (neuler/=0) THEN
618            ub_b (:,:) = 0.e0
619            vb_b (:,:) = 0.e0
620
621            ! vertical sum
622            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
623               DO jk = 1, jpkm1
624                  DO ji = 1, jpij
625                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
626                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
627                  END DO
628               END DO
629            ELSE                             ! No  vector opt.
630               DO jk = 1, jpkm1
631                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
632                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
633               END DO
634            ENDIF
635
636            IF( lk_vvl ) THEN
637               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
638               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
639            ELSE
640               ub_b(:,:) = ub_b(:,:) * hur(:,:)
641               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
642            ENDIF
643         ELSE                                 ! neuler==0
644            ub_b (:,:) = un_b (:,:)
645            vb_b (:,:) = vn_b (:,:)
646         ENDIF
647
648         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
649            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered extrenal ssh
650         ELSE
651            sshn_b(:,:)=sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
652         ENDIF
653      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
654         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
655         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
656         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !
657      ENDIF
658      !
659   END SUBROUTINE ts_rst
660
661#else
662   !!----------------------------------------------------------------------
663   !!   Default case :   Empty module   No standart free surface cst volume
664   !!----------------------------------------------------------------------
665CONTAINS
666   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
667      dyn_spg_ts_alloc = 0
668   END FUNCTION dyn_spg_ts_alloc
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.