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 trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90 @ 247

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

CL : Add CVS Header and CeCILL licence information

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.2 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 obc_oce         ! Lateral open boundary condition
30   USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
31   USE obcvol          ! ocean open boundary condition (obc_vol routines)
32   USE lib_mpp         ! distributed memory computing library
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE cla_dynspg      ! cross land advection
35
36   IMPLICIT NONE
37   PRIVATE
38
39   !! * Accessibility
40   PUBLIC dyn_spg_fsc  ! routine called by step.F90
41
42   !! * Shared module variables
43   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_fsc = .TRUE.    !: free surface constant volume flag
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_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_fsc( kt, kindic )
57      !!----------------------------------------------------------------------
58      !!                  ***  routine dyn_spg_fsc  ***
59      !!
60      !! ** Purpose :   Compute the now trend due to the surface pressure
61      !!      gradient in case of free surface formulation with a constant
62      !!      ocean volume 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      !!
97      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
98      !!             - Save the trends in (ztdua,ztdva) ('key_trddyn')
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      !!---------------------------------------------------------------------
110      !! * Modules used     
111      USE oce, ONLY :    ztdua => ta,      & ! use ta as 3D workspace   
112                         ztdva => sa         ! use sa as 3D workspace   
113
114      !! * Arguments
115      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
116      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag
117                                             ! if the solver doesn t converge
118                                             ! the flag is < 0
119      !! * Local declarations
120      INTEGER  ::   ji, jj, jk               ! dummy loop indices
121      REAL(wp) ::                         & 
122         z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars
123         znurau, zssha, zspgu, zspgv,     &  !   "          "
124         zgcb, zbtd,                      &  !   "          "
125         ztdgu, ztdgv                        !   "          "
126      !!----------------------------------------------------------------------
127
128      IF( kt == nit000 ) THEN
129         IF(lwp) WRITE(numout,*)
130         IF(lwp) WRITE(numout,*) 'dyn_spg_fsc : surface pressure gradient trend'
131         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
132       
133         ! set to zero free surface specific arrays
134         spgu(:,:) = 0.e0                     ! surface pressur gradient (i-direction)
135         spgv(:,:) = 0.e0                     ! surface pressur gradient (j-direction)
136      ENDIF
137
138      ! 0. Local constant initialization
139      ! --------------------------------
140      ! time step: leap-frog
141      z2dt = 2. * rdt
142      ! time step: Euler if restart from rest
143      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
144      ! coefficients
145      z2dtg  = grav * z2dt
146      zraur  = 1. / rauw
147      znugdt =  rnu * grav * z2dt
148      znurau =  znugdt * zraur
149
150      ! 1. Surface pressure gradient (now)
151      ! ----------------------------
152      DO jj = 2, jpjm1
153         DO ji = fs_2, fs_jpim1   ! vector opt.
154            zspgu =      - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
155            zspgv =      - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
156            spgu(ji,jj) = zspgu
157            spgv(ji,jj) = zspgv
158         END DO
159      END DO 
160
161      ! 2. Add the surface pressure trend to the general trend
162      ! ------------------------------------------------------
163      DO jk = 1, jpkm1
164         DO jj = 2, jpjm1
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
167               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
168            END DO
169         END DO
170      END DO
171
172      ! Save the surface pressure gradient trend for diagnostics
173      IF( l_trddyn )   THEN
174         DO jk = 1, jpkm1
175            ztdua(:,:,jk) = spgu(:,:) 
176            ztdva(:,:,jk) = spgv(:,:) 
177         END DO
178      ENDIF
179     
180      ! 1. Evaluate the masked next velocity
181      ! ------------------------------------
182      !     (effect of the additional force not included)
183      DO jk = 1, jpkm1
184         DO jj = 2, jpjm1
185            DO ji = fs_2, fs_jpim1   ! vector opt.
186               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
187               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
188            END DO
189         END DO
190      END DO
191#if defined key_obc
192      ! Update velocities on each open boundary with the radiation algorithm
193      CALL obc_dyn( kt )
194      ! Correction of the barotropic componant velocity to control the volume of the system
195      CALL obc_vol( kt )
196#endif
197#if defined key_orca_r2
198      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (update (ua,va))
199#endif
200
201      ! 2. compute the next vertically averaged velocity
202      ! ------------------------------------------------
203      !     (effect of the additional force not included)
204      ! initialize to zero
205      DO jj = 2, jpjm1
206         DO ji = fs_2, fs_jpim1   ! vector opt.
207            spgu(ji,jj) = 0.e0
208            spgv(ji,jj) = 0.e0
209         END DO
210      END DO
211
212      ! vertical sum
213!CDIR NOLOOPCHG
214      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
215         DO jk = 1, jpkm1
216            DO ji = 1, jpij
217               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
218               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
219            END DO
220         END DO
221      ELSE                        ! No  vector opt.
222         DO jk = 1, jpkm1
223            DO jj = 2, jpjm1
224               DO ji = 2, jpim1
225                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
226                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
227               END DO
228            END DO
229         END DO
230      ENDIF
231
232      ! transport: multiplied by the horizontal scale factor
233      DO jj = 2, jpjm1
234         DO ji = fs_2, fs_jpim1   ! vector opt.
235            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
236            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
237         END DO
238      END DO
239
240      ! Boundary conditions on (spgu,spgv)
241      CALL lbc_lnk( spgu, 'U', -1. )
242      CALL lbc_lnk( spgv, 'V', -1. )
243
244      ! 3. Right hand side of the elliptic equation and first guess
245      ! -----------------------------------------------------------
246      DO jj = 2, jpjm1
247         DO ji = fs_2, fs_jpim1   ! vector opt.
248            ! Divergence of the after vertically averaged velocity
249            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
250                  + spgv(ji,jj) - spgv(ji,jj-1)
251            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
252            ! First guess of the after barotropic transport divergence
253            zbtd = gcx(ji,jj)
254            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
255            gcxb(ji,jj) =      zbtd
256         END DO
257      END DO
258
259      ! 4. Relative precision (computation on one processor)
260      ! ---------------------
261      rnorme =0.
262      rnorme = SUM( gcb(:,:) * gcdmat(:,:) * gcb(:,:) )
263      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
264
265      epsr = eps * eps * rnorme
266      ncut = 0
267      ! if rnorme is 0, the solution is 0, the solver isn't called
268      IF( rnorme == 0.e0 ) THEN
269         gcx(:,:) = 0.e0
270         res   = 0.e0
271         niter = 0
272         ncut  = 999
273      ENDIF
274
275      ! 5. Evaluate the next transport divergence
276      ! -----------------------------------------
277      !    Iterarive solver for the elliptic equation (except IF sol.=0)
278      !    (output in gcx with boundary conditions applied)
279      kindic = 0
280      IF( ncut == 0 ) THEN
281         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
282            CALL sol_pcg( kindic )
283         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
284            CALL sol_sor( kindic )
285         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
286            CALL sol_fet( kindic )
287         ELSE                          ! e r r o r in nsolv namelist parameter
288            IF(lwp) WRITE(numout,cform_err)
289            IF(lwp) WRITE(numout,*) ' dyn_spg_fsc : e r r o r, nsolv = 1, 2 or 3'
290            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~                not = ', nsolv
291            nstop = nstop + 1
292         ENDIF
293      ENDIF
294
295      ! 6. Transport divergence gradient multiplied by z2dt
296      ! -----------------------------------------------====
297      DO jj = 2, jpjm1
298         DO ji = fs_2, fs_jpim1   ! vector opt.
299            ! trend of Transport divergence gradient
300            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
301            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
302            ! multiplied by z2dt
303#if defined key_obc
304            ! caution : grad D = 0 along open boundaries
305            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
306            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
307#else
308            spgu(ji,jj) = z2dt * ztdgu
309            spgv(ji,jj) = z2dt * ztdgv
310#endif
311         END DO
312      END DO
313
314      ! 7.  Add the trends multiplied by z2dt to the after velocity
315      ! -----------------------------------------------------------
316      !     ( c a u t i o n : (ua,va) here are the after velocity not the
317      !                       trend, the leap-frog time stepping will not
318      !                       be done in dynnxt.F routine)
319      DO jk = 1, jpkm1
320         DO jj = 2, jpjm1
321            DO ji = fs_2, fs_jpim1   ! vector opt.
322               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
323               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
324            END DO
325         END DO
326      END DO
327
328      ! save the surface pressure gradient trends for diagnostic
329      ! momentum trends
330      IF( l_trddyn )   THEN
331         DO jk = 1, jpkm1
332            ztdua(:,:,jk) = ztdua(:,:,jk) + spgu(:,:)/z2dt 
333            ztdva(:,:,jk) = ztdva(:,:,jk) + spgv(:,:)/z2dt 
334         END DO
335
336         CALL trd_mod(ztdua, ztdva, jpdtdspg, 'DYN', kt)
337      ENDIF
338
339      IF(l_ctl) THEN         ! print sum trends (used for debugging)
340         WRITE(numout,*) ' spg  - Ua: ', SUM( ua(2:nictl,2:njctl,1:jpkm1)*umask(2:nictl,2:njctl,1:jpkm1) ),   &
341            &                   ' Va: ', SUM( va(2:nictl,2:njctl,1:jpkm1)*vmask(2:nictl,2:njctl,1:jpkm1) )
342      ENDIF
343
344
345      ! 8. Sea surface elevation time stepping
346      ! --------------------------------------
347      ! Euler (forward) time stepping, no time filter
348      IF( neuler == 0 .AND. kt == nit000 ) THEN
349         DO jj = 1, jpj
350            DO ji = 1, jpi
351               ! after free surface elevation
352               zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
353               ! swap of arrays
354               sshb(ji,jj) = sshn(ji,jj)
355               sshn(ji,jj) = zssha
356            END DO
357         END DO
358      ELSE
359         ! Leap-frog time stepping and time filter
360         DO jj = 1, jpj
361            DO ji = 1, jpi
362               ! after free surface elevation
363               zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
364               ! time filter and array swap
365               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
366               sshn(ji,jj) = zssha
367            END DO
368         END DO
369      ENDIF
370
371
372      IF(l_ctl) THEN         ! print sum trends (used for debugging)
373         WRITE(numout,*) ' spg - ssh:', SUM( sshn(2:nictl,2:njctl)*tmask(2:nictl,2:njctl,1) )
374      ENDIF
375
376   END SUBROUTINE dyn_spg_fsc
377
378#else
379   !!----------------------------------------------------------------------
380   !!   Default case :   Empty module   No standart free surface cst volume
381   !!----------------------------------------------------------------------
382   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_fsc = .FALSE.   !: free surface constant volume flag
383CONTAINS
384   SUBROUTINE dyn_spg_fsc( kt, kindic )       ! Empty routine
385      WRITE(*,*) 'dyn_spg_fsc: You should not have seen this print! error?', kt, kindic
386   END SUBROUTINE dyn_spg_fsc
387#endif
388   
389   !!======================================================================
390END MODULE dynspg_fsc
Note: See TracBrowser for help on using the repository browser.