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_atsk.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.7 KB
Line 
1MODULE dynspg_fsc_atsk
2   !!======================================================================
3   !!                  ***  MODULE  dynspg_fsc_atsk  ***
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' & 'key_autotasking'   free surface cst volume
9   !!                                                   j-k-i loop (j-slab)
10   !!----------------------------------------------------------------------
11   !!   dyn_spg_fsc_atsk : 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 , LODYC-IPSL   (2003)
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE trdtra_oce     ! ocean active tracer trend
21   USE trddyn_oce     ! ocean dynamics trend
22   USE in_out_manager  ! I/O manager
23   USE phycst          ! physical constant
24   USE ocesbc          ! Ocen Surface Boundary condition
25   USE flxrnf          ! ???
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 obc_oce         ! Lateral open boundary condition
31   USE obcdyn          ! open boudary condition
32   USE obcdyn          !    "              "
33   USE obcvol          !    "              "
34   USE lib_mpp         ! ???
35   USE lbclnk          ! ???
36   USE cla_dynspg      ! ???
37
38   IMPLICIT NONE
39   PRIVATE
40
41   !! * Accessibility
42   PUBLIC dyn_spg_fsc_atsk        ! routine called by step.F90
43
44   !! * Shares module variables
45   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_fsc_tsk = .TRUE.    ! free surf. cst vol. flag
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49   !!----------------------------------------------------------------------
50   !!   OPA 9.0 , LODYC-IPSL  (2003)
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE dyn_spg_fsc_atsk( kt, kindic )
56      !!----------------------------------------------------------------------
57      !!                  ***  routine dyn_spg_fsc_atsk  ***
58      !!
59      !! ** Purpose :   Compute the now trend due to the surface pressure
60      !!      gradient for free surface formulation with a constant ocean
61      !!      volume case, add it to the general trend of momentum equation.
62      !!
63      !! ** Method  :   Free surface formulation. The surface pressure gradient
64      !!      is given by:
65      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( etn + rnu btda )
66      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( etn + rnu btda )
67      !!      where etn is the free surface elevation and btda is the after
68      !!      of the free surface elevation
69      !!       -1- compute the after sea surface elevation from the cinematic
70      !!      surface boundary condition:
71      !!              zssha = sshb + 2 rdt ( wn - emp )
72      !!           Time filter applied on now sea surface elevation to avoid
73      !!      the divergence of two consecutive time-steps and swap of free
74      !!      surface arrays to start the next time step:
75      !!              sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
76      !!              sshn = zssha
77      !!       -2- evaluate the surface presure trend (including the addi-
78      !!      tional force) in three steps:
79      !!        a- compute the right hand side of the elliptic equation:
80      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
81      !!         where (spgu,spgv) are given by:
82      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
83      !!                 - g 2 rdt hu /e1u di[sshn + emp]
84      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
85      !!                 - g 2 rdt hv /e2v dj[sshn + emp]
86      !!         and define the first guess from previous computation :
87      !!            zbtd = btda
88      !!            btda = 2 zbtd - btdb
89      !!            btdb = zbtd
90      !!        b- compute the relative accuracy to be reached by the
91      !!         iterative solver
92      !!        c- apply the solver by a call to sol... routine
93      !!       -3- compute and add the free surface pressure gradient inclu-
94      !!      ding the additional force used to stabilize the equation.
95      !!      several slabs used ('key-autotasking')
96      !!
97      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
98      !!             - Save the trends in (utrd,vtrd) ('key_diatrends')
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      !!---------------------------------------------------------------------
109      !! * Arguments
110      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
111      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag
112                                             ! if the solver doesn t converge
113                                             ! the flag is < 0
114      !! * Local declarations
115      INTEGER  ::   ji, jj, jk               ! dummy loop indices
116      REAL(wp) ::   &              ! temporary scalars
117         z2dt, z2dtg, zraur, znugdt, znurau,   &
118         zssha, zspgu, zspgv, zegu, zegv,   &
119         zgcb, zbtd, ztdgu, ztdgv, zgwgt
120      !!----------------------------------------------------------------------
121
122      IF( kt == nit000 ) THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'dyn_spg_fsc_atsk : surface pressure gradient trend'
125         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~  (free surface constant volume, autotasking case)'
126
127         ! set to zero free surface specific arrays
128         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction)
129         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
130         IF( .NOT.ln_rstart ) THEN
131            sshb(:,:) = 0.e0      ! before sea-surface height
132            sshn(:,:) = 0.e0      ! now    sea-surface height
133         ENDIF
134      ENDIF
135
136      ! 0. Local constant initialization
137      ! --------------------------------
138      ! time step: leap-frog
139      z2dt = 2. * rdt
140      ! time step: Euler if restart from rest
141      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
142      ! coefficients
143      z2dtg  = g * z2dt
144      zraur  = 1. / rauw
145      znugdt =  rnu * g * z2dt
146      znurau =  znugdt * zraur
147#if defined key_mpp
148      ! Mpp : export boundary values of to neighboring processors
149      !!bug ???  why only in mpp?  is it really needed???
150      CALL lbc_lnk( ua, 'U' , -1. )
151      CALL lbc_lnk( va, 'V' , -1. )
152#endif
153      !                                                ! ===============
154      DO jj = 2, jpjm1                                 !  Vertical slab
155         !                                             ! ===============
156         ! 1. Surface pressure gradient (now)
157         ! ----------------------------
158         DO ji = 2, jpim1
159            zspgu =      - g * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
160            zspgv =      - g * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
161            zegu  = + znurau * ( emp (ji+1,jj) - emp (ji,jj) ) / e1u(ji,jj)
162            zegv  = + znurau * ( emp (ji,jj+1) - emp (ji,jj) ) / e2v(ji,jj)
163            spgu(ji,jj) = zspgu + zegu
164            spgv(ji,jj) = zspgv + zegv
165#if defined key_trddyn
166            ! save the surface pressure gradient trends
167            utrd2(ji,jj,1) = zspgu
168            vtrd2(ji,jj,1) = zspgv
169            ! save the first part of the additional force trend
170            utrd2(ji,jj,2) = zegu
171            vtrd2(ji,jj,2) = zegv
172#endif
173         END DO 
174
175         ! 2. Add the surface pressure trend to the general trend
176         ! ------------------------------------------------------
177         DO jk = 1, jpkm1
178            DO ji = 2, jpim1
179               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
180               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
181#if defined key_trddyn
182               ! save the surface pressure gradient trend for diagnostics
183               utrd(ji,jj,jk,8) = spgu(ji,jj)
184               vtrd(ji,jj,jk,8) = spgv(ji,jj) 
185#endif
186            END DO
187         END DO
188         !                                             ! ===============
189      END DO                                           !   End of slab
190      !                                                ! ===============
191
192      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
193
194      !                                                ! ===============
195      DO jj = 2, jpjm1                                 !  Vertical slab
196         !                                             ! ===============
197
198         ! 1. Evaluate the masked next velocity
199         ! ------------------------------------
200         !     (effect of the additional force not included)
201
202         DO jk = 1, jpkm1
203            DO ji = 2, jpim1
204               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
205               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
206            END DO
207         END DO
208
209         !                                             ! ===============
210      END DO                                           !   End of slab
211      !                                                ! ===============
212
213#if defined key_obc
214         ! Update velocities on each open boundary with the radiation algorithm
215         CALL obc_dyn( kt )
216         ! Correction of the barotropic componant velocity to control the volume of the system
217         CALL obc_vol( kt )
218#endif
219#if defined key_orca_r2
220      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (Update (ua,va))
221#endif
222
223      !                                                ! ===============
224      DO jj = 2, jpjm1                                 !  Vertical slab
225         !                                             ! ===============
226
227         ! 2. compute the next vertically averaged velocity
228         ! ------------------------------------------------
229         !     (effect of the additional force not included)
230         ! initialize to zero
231         DO ji = 2, jpim1
232            spgu(ji,jj) = 0.e0
233            spgv(ji,jj) = 0.e0
234         END DO
235
236         ! vertical sum
237         DO jk = 1, jpk
238            DO ji = 2, jpim1
239               spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
240               spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
241            END DO
242         END DO
243
244         ! transport: multiplied by the horizontal scale factor
245         DO ji = 2, jpim1
246            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
247            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
248         END DO
249
250         !                                             ! ===============
251      END DO                                           !   End of slab
252      !                                                ! ===============
253
254      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
255     
256      ! Boundary conditions on (spgu,spgv)
257      CALL lbc_lnk( spgu, 'U', -1. )
258      CALL lbc_lnk( spgv, 'V', -1. )
259
260      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
261
262      ! 3. Right hand side of the elliptic equation and first guess
263      ! -----------------------------------------------------------
264      DO jj = 2, jpjm1
265         DO ji = 2, jpim1
266            ! Divergence of the after vertically averaged velocity
267            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
268                  + spgv(ji,jj) - spgv(ji,jj-1)
269            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
270            ! First guess of the after barotropic transport divergence
271            zbtd = gcx(ji,jj)
272            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
273            gcxb(ji,jj) =      zbtd
274         END DO
275      END DO
276
277      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
278     
279      ! 4. Relative precision (computation on one processor)
280      ! ---------------------
281      rnorme =0.
282      DO jj = 1, jpj
283         DO ji = 1, jpi
284            zgwgt  = gcdmat(ji,jj) * gcb(ji,jj)
285            rnorme = rnorme + gcb(ji,jj) * zgwgt
286         END DO
287      END DO
288#if defined key_mpp
289      CALL mpp_sum( rnorme )
290#endif
291      epsr = eps * eps * rnorme
292      ncut = 0
293      ! if rnorme is 0, the solution is 0, the solver isn't called
294      IF( rnorme == 0.e0 ) THEN
295         gcx(:,:) = 0.e0
296         res   = 0.e0
297         niter = 0
298         ncut  = 999
299      ENDIF
300     
301      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
302     
303      ! 5. Evaluate the next transport divergence
304      ! -----------------------------------------
305      !    Iterarive solver for the elliptic equation (except IF sol.=0)
306      !    (output in gcx with boundary conditions applied)
307      kindic = 0
308      IF( ncut == 0 ) THEN
309         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
310            CALL sol_pcg( kindic )
311         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
312            CALL sol_sor( kt, kindic )
313         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
314            CALL sol_fet( kindic )
315         ELSE                          ! e r r o r in nsolv namelist parameter
316            IF(lwp) WRITE(numout,cform_err)
317            IF(lwp) WRITE(numout,*) ' dyn_spg_fsc_atsk : e r r o r, nsolv = 1, 2 or 3'
318            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~                not = ', nsolv
319            nstop = nstop + 1
320         ENDIF
321      ENDIF
322
323      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
324
325      !                                                ! ===============
326      DO jj = 2, jpjm1                                 !  Vertical slab
327         !                                             ! ===============
328         
329         ! 6. Transport divergence gradient multiplied by z2dt
330         ! -----------------------------------------------====
331         DO ji = 2, jpim1
332            ! trend of Transport divergence gradient
333            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
334            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
335            ! multiplied by z2dt
336#if defined key_obc
337            ! caution : grad D = 0 along open boundaries
338            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
339            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
340#else
341            spgu(ji,jj) = z2dt * ztdgu
342            spgv(ji,jj) = z2dt * ztdgv
343#endif
344#if defined key_trddyn
345            ! save the transport divergence gradient trends
346            utrd2(ji,jj,2) = utrd2(ji,jj,2) + ztdgu
347            vtrd2(ji,jj,2) = vtrd2(ji,jj,2) + ztdgv
348#endif
349         END DO
350         
351         ! 7.  Add the trends multiplied by z2dt to the after velocity
352         ! -----------------------------------------------------------
353         !     ( c a u t i o n : (ua,va) here are the after velocity not the
354         !                       trend, the leap-frog time stepping will not
355         !                       be done in dynnxt.F routine)
356         DO jk = 1, jpkm1
357            DO ji = 2, jpim1
358               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
359               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
360#if defined key_trddyn
361               ! save the surface pressure gradient trend for diagnostics
362               utrd(ji,jj,jk,8) = utrd(ji,jj,jk,8) + spgu(ji,jj)/z2dt
363               vtrd(ji,jj,jk,8) = vtrd(ji,jj,jk,8) + spgv(ji,jj)/z2dt
364#endif
365            END DO
366         END DO
367
368         ! 8. Sea surface elevation time stepping
369         ! --------------------------------------
370         ! Euler (forward) time stepping, no time filter
371         IF( neuler == 0 .AND. kt == nit000 ) THEN
372            DO ji = 1, jpi
373               ! after free surface elevation
374               zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
375               ! swap of arrays
376               sshb(ji,jj) = sshn(ji,jj)
377               sshn(ji,jj) = zssha
378            END DO
379         ELSE
380            ! Leap-frog time stepping and time filter
381            DO ji = 1, jpi
382               ! after free surface elevation
383               zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
384               ! time filter and array swap
385               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
386               sshn(ji,jj) = zssha
387            END DO
388         ENDIF
389         !                                             ! ===============
390      END DO                                           !   End of slab
391      !                                                ! ===============
392
393      !Boundary conditions on sshn
394      CALL lbc_lnk( sshn, 'T', 1. )
395
396   END SUBROUTINE dyn_spg_fsc_atsk
397
398#else
399   !!----------------------------------------------------------------------
400   !!   Default case :                                         Empty module
401   !!----------------------------------------------------------------------
402   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_fsc_tsk = .FALSE.   ! free surf. cst vol. flag
403CONTAINS
404   SUBROUTINE dyn_spg_fsc_atsk( kt, kindic )      ! Empty module
405      WRITE(*,*) kt, kindic
406   END SUBROUTINE dyn_spg_fsc_atsk
407#endif
408   
409   !!======================================================================
410END MODULE dynspg_fsc_atsk
Note: See TracBrowser for help on using the repository browser.