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

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

nemo_v1_update_064 : CT : general trends update including the addition of mean windows analysis possibility in the mixed layer

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