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

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

CT : UPDATE152 : optimization: reduce the MPI communications in doing only one CALL lbc_lnk() with a 3D array instead of jpk CALL lbc_lnk() with a 2D array; or suppressing useless CALL lbc_lnk()

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