source: trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90 @ 358

Last change on this file since 358 was 358, checked in by opalod, 15 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.4 KB
Line 
1MODULE dynspg_flt_jki
2   !!======================================================================
3   !!                  ***  MODULE  dynspg_flt_jki  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if ( defined key_dynspg_flt && defined key_autotasking )   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_flt' & 'key_autotasking'          filtered free surface
9   !!                                                   j-k-i loop (j-slab)
10   !!----------------------------------------------------------------------
11   !!   dyn_spg_flt_jki : Update the momentum trend with the surface pressure
12   !!                     gradient for the free surf. constant volume case
13   !!                     with auto-tasking (j-slab) (no vectior opt.)
14   !!----------------------------------------------------------------------
15   !!   OPA 9.0 , LOCEAN-IPSL (2005)
16   !! $Header$
17   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
18   !!----------------------------------------------------------------------
19   !! * Modules used
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE zdf_oce         ! ocean vertical physics
23   USE phycst          ! physical constant
24   USE ocesbc          ! Ocean Surface Boundary condition
25   USE flxrnf          ! ocean runoffs
26   USE sol_oce         ! ocean elliptic solver
27   USE solpcg          ! preconditionned conjugate gradient solver
28   USE solsor          ! Successive Over-relaxation solver
29   USE solfet          ! FETI solver
30   USE solsor_e        ! Successive Over-relaxation solver with MPP optimization
31   USE obc_oce         ! Lateral open boundary condition
32   USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
33   USE obcvol          ! ocean open boundary condition (obc_vol routines)
34   USE lib_mpp         ! distributed memory computing library
35   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
36   USE cla_dynspg      ! cross land advection
37   USE prtctl          ! Print control
38   USE in_out_manager  ! I/O manager
39
40   IMPLICIT NONE
41   PRIVATE
42
43   !! * Accessibility
44   PUBLIC dyn_spg_flt_jki        ! routine called by step.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48   !!----------------------------------------------------------------------
49   !!   OPA 9.0 , LOCEAN-IPSL (2005)
50   !! $Header$
51   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE dyn_spg_flt_jki( kt, kindic )
57      !!----------------------------------------------------------------------
58      !!                  ***  routine dyn_spg_flt_jki  ***
59      !!
60      !! ** Purpose :   Compute the now trend due to the surface pressure
61      !!      gradient for free surface formulation with a constant ocean
62      !!      volume case, add it to the general trend of momentum equation.
63      !!
64      !! ** Method  :   Free surface formulation. The surface pressure gradient
65      !!      is given by:
66      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( etn + rnu btda )
67      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( etn + rnu btda )
68      !!      where etn is the free surface elevation and btda is the after
69      !!      of the free surface elevation
70      !!       -1- compute the after sea surface elevation from the cinematic
71      !!      surface boundary condition:
72      !!              zssha = sshb + 2 rdt ( wn - emp )
73      !!           Time filter applied on now sea surface elevation to avoid
74      !!      the divergence of two consecutive time-steps and swap of free
75      !!      surface arrays to start the next time step:
76      !!              sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
77      !!              sshn = zssha
78      !!       -2- evaluate the surface presure trend (including the addi-
79      !!      tional force) in three steps:
80      !!        a- compute the right hand side of the elliptic equation:
81      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
82      !!         where (spgu,spgv) are given by:
83      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
84      !!                 - grav 2 rdt hu /e1u di[sshn + emp]
85      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
86      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
87      !!         and define the first guess from previous computation :
88      !!            zbtd = btda
89      !!            btda = 2 zbtd - btdb
90      !!            btdb = zbtd
91      !!        b- compute the relative accuracy to be reached by the
92      !!         iterative solver
93      !!        c- apply the solver by a call to sol... routine
94      !!       -3- compute and add the free surface pressure gradient inclu-
95      !!      ding the additional force used to stabilize the equation.
96      !!      several slabs used ('key-autotasking')
97      !!
98      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
99      !!
100      !! References :
101      !!      Roullet and Madec 1999, JGR.
102      !!
103      !! History :
104      !!        !  98-05 (G. Roullet)  Original code
105      !!        !  98-10 (G. Madec, M. Imbard)  release 8.2
106      !!   8.5  !  02-08 (G. Madec)  F90: Free form and module
107      !!        !  02-11 (C. Talandier, A-M Treguier) Open boundaries
108      !!   9.0  !  04-08 (C. Talandier) New trends organization
109      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
110      !!---------------------------------------------------------------------
111      !! * Arguments
112      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
113      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag
114                                             ! if the solver doesn t converge
115                                             ! the flag is < 0
116      !! * Local declarations
117      INTEGER  ::   ji, jj, jk               ! dummy loop indices
118      REAL(wp) ::   &              ! temporary scalars
119         z2dt, z2dtg, zraur, znugdt, znurau,   &
120         zssha, zgcb, zbtd, ztdgu, ztdgv, zgwgt
121      !!----------------------------------------------------------------------
122
123      IF( kt == nit000 ) THEN
124         IF(lwp) WRITE(numout,*)
125         IF(lwp) WRITE(numout,*) 'dyn_spg_flt_jki : surface pressure gradient trend'
126         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~  (free surface constant volume, autotasking case)'
127
128         ! set to zero free surface specific arrays
129         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction)
130         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
131      ENDIF
132
133      ! 0. Local constant initialization
134      ! --------------------------------
135      ! time step: leap-frog
136      z2dt = 2. * rdt
137      ! time step: Euler if restart from rest
138      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
139      ! coefficients
140      z2dtg  = grav * z2dt
141      zraur  = 1. / rauw
142      znugdt =  rnu * grav * z2dt
143      znurau =  znugdt * zraur
144
145      !                                                ! ===============
146      DO jj = 2, jpjm1                                 !  Vertical slab
147         !                                             ! ===============
148         ! Surface pressure gradient (now)
149         DO ji = 2, jpim1
150            spgu(ji,jj) = - grav * ( sshn(ji+1,jj  ) - sshn(ji,jj) ) / e1u(ji,jj)
151            spgv(ji,jj) = - grav * ( sshn(ji  ,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
152         END DO 
153
154         !  Add the surface pressure trend to the general trend
155         DO jk = 1, jpkm1
156            DO ji = 2, jpim1
157               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
158               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
159            END DO
160         END DO
161
162         ! Evaluate the masked next velocity (effect of the additional force not included)
163         DO jk = 1, jpkm1
164            DO ji = 2, jpim1
165               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
166               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
167            END DO
168         END DO
169
170         !                                             ! ===============
171      END DO                                           !   End of slab
172      !                                                ! ===============
173
174#if defined key_obc
175         ! Update velocities on each open boundary with the radiation algorithm
176         CALL obc_dyn( kt )
177         ! Correction of the barotropic componant velocity to control the volume of the system
178         CALL obc_vol( kt )
179#endif
180#if defined key_orca_r2
181      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (Update (ua,va))
182#endif
183
184      !                                                ! ===============
185      DO jj = 2, jpjm1                                 !  Vertical slab
186         !                                             ! ===============
187
188         ! 2. compute the next vertically averaged velocity
189         ! ------------------------------------------------
190         !     (effect of the additional force not included)
191         ! initialize to zero
192         DO ji = 2, jpim1
193            spgu(ji,jj) = 0.e0
194            spgv(ji,jj) = 0.e0
195         END DO
196
197         ! vertical sum
198         DO jk = 1, jpk
199            DO ji = 2, jpim1
200               spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
201               spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
202            END DO
203         END DO
204
205         ! transport: multiplied by the horizontal scale factor
206         DO ji = 2, jpim1
207            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
208            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
209         END DO
210
211         !                                             ! ===============
212      END DO                                           !   End of slab
213      !                                                ! ===============
214
215      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
216     
217      ! Boundary conditions on (spgu,spgv)
218      CALL lbc_lnk( spgu, 'U', -1. )
219      CALL lbc_lnk( spgv, 'V', -1. )
220
221      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
222
223      ! 3. Right hand side of the elliptic equation and first guess
224      ! -----------------------------------------------------------
225      DO jj = 2, jpjm1
226         DO ji = 2, jpim1
227            ! Divergence of the after vertically averaged velocity
228            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
229                  + spgv(ji,jj) - spgv(ji,jj-1)
230            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
231            ! First guess of the after barotropic transport divergence
232            zbtd = gcx(ji,jj)
233            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
234            gcxb(ji,jj) =      zbtd
235         END DO
236      END DO
237      ! applied the lateral boundary conditions
238      IF( nsolv == 4)   CALL lbc_lnk_e( gcb, c_solver_pt, 1. ) 
239
240      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
241     
242      ! 4. Relative precision (computation on one processor)
243      ! ---------------------
244      rnorme =0.
245      DO jj = 1, jpj
246         DO ji = 1, jpi
247            zgwgt  = gcdmat(ji,jj) * gcb(ji,jj)
248            rnorme = rnorme + gcb(ji,jj) * zgwgt * bmask(ji,jj)
249         END DO
250      END DO
251      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
252
253      epsr = eps * eps * rnorme
254      ncut = 0
255      ! if rnorme is 0, the solution is 0, the solver isn't called
256      IF( rnorme == 0.e0 ) THEN
257         gcx(:,:) = 0.e0
258         res   = 0.e0
259         niter = 0
260         ncut  = 999
261      ENDIF
262     
263      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
264     
265      ! 5. Evaluate the next transport divergence
266      ! -----------------------------------------
267      !    Iterarive solver for the elliptic equation (except IF sol.=0)
268      !    (output in gcx with boundary conditions applied)
269      kindic = 0
270      IF( ncut == 0 ) THEN
271         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
272            CALL sol_pcg( kindic )
273         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
274            CALL sol_sor( kindic )
275         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
276            CALL sol_fet( kindic )
277         ELSEIF( nsolv == 4 ) THEN     ! successive-over-relaxation with extra outer halo
278            CALL sol_sor_e( kindic )
279         ELSE                          ! e r r o r in nsolv namelist parameter
280            IF(lwp) WRITE(numout,cform_err)
281            IF(lwp) WRITE(numout,*) ' dyn_spg_flt_jki : e r r o r, nsolv = 1, 2, 3 or 4'
282            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~                not = ', nsolv
283            nstop = nstop + 1
284         ENDIF
285      ENDIF
286
287      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
288
289!CDIR PARALLEL DO
290!$OMP PARALLEL DO
291      !                                                ! ===============
292      DO jj = 2, jpjm1                                 !  Vertical slab
293         !                                             ! ===============
294         
295         ! 6. Transport divergence gradient multiplied by z2dt
296         ! -----------------------------------------------====
297         DO ji = 2, jpim1
298            ! trend of Transport divergence gradient
299            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
300            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
301            ! multiplied by z2dt
302#if defined key_obc
303            ! caution : grad D = 0 along open boundaries
304            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
305            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
306#else
307            spgu(ji,jj) = z2dt * ztdgu
308            spgv(ji,jj) = z2dt * ztdgv
309#endif
310         END DO
311         
312         ! 7.  Add the trends multiplied by z2dt to the after velocity
313         ! -----------------------------------------------------------
314         !     ( c a u t i o n : (ua,va) here are the after velocity not the
315         !                       trend, the leap-frog time stepping will not
316         !                       be done in dynnxt.F routine)
317         DO jk = 1, jpkm1
318            DO ji = 2, jpim1
319               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
320               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
321            END DO
322         END DO
323
324         ! 8. Sea surface elevation time stepping
325         ! --------------------------------------
326         ! Euler (forward) time stepping, no time filter
327         IF( neuler == 0 .AND. kt == nit000 ) THEN
328            DO ji = 1, jpi
329               ! after free surface elevation
330               zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
331               ! swap of arrays
332               sshb(ji,jj) = sshn(ji,jj)
333               sshn(ji,jj) = zssha
334            END DO
335         ELSE
336            ! Leap-frog time stepping and time filter
337            DO ji = 1, jpi
338               ! after free surface elevation
339               zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
340               ! time filter and array swap
341               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
342               sshn(ji,jj) = zssha
343            END DO
344         ENDIF
345         !                                             ! ===============
346      END DO                                           !   End of slab
347      !                                                ! ===============
348
349      !Boundary conditions on sshn
350      CALL lbc_lnk( sshn, 'T', 1. )
351
352      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
353         CALL prt_ctl( tab3d_1=ua  , clinfo1=' spg  - Ua : ', mask1=umask,   &
354            &          tab3d_2=va  , clinfo2='        Va : ', mask2=vmask )
355         CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask)
356      ENDIF
357
358
359   END SUBROUTINE dyn_spg_flt_jki
360
361#else
362   !!----------------------------------------------------------------------
363   !!   Default case :                                         Empty module
364   !!----------------------------------------------------------------------
365CONTAINS
366   SUBROUTINE dyn_spg_flt_jki( kt, kindic )      ! Empty module
367      WRITE(*,*) 'dyn_spg_flt_jki: You should not have seen this print! error?', kt, kindic
368   END SUBROUTINE dyn_spg_flt_jki
369#endif
370   
371   !!======================================================================
372END MODULE dynspg_flt_jki
Note: See TracBrowser for help on using the repository browser.