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 @ 294

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

nemo_v1_update_007 : CT : comments correction

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