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_fsc.F90 in tags/nemo_v1_06/NEMO/OPA_SRC/DYN – NEMO

source: tags/nemo_v1_06/NEMO/OPA_SRC/DYN/dynspg_fsc.F90 @ 7389

Last change on this file since 7389 was 314, checked in by opalod, 19 years ago

nemo_v1_update_017:RB: added call to sol_sor_e for free surface and rigid lid.

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