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

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

RB:nemo_v1_update_038: first integration of Agrif :

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