source: trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 389

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

RB:nemo_v1_update_038: first integration of Agrif :

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