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

source: trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90 @ 109

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

CT : UPDATE069 : Vorticity diagnostics has been added

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.4 KB
Line 
1MODULE dynspg_rl
2   !!======================================================================
3   !!                    ***  MODULE  dynspg_rl  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if   defined key_dynspg_rl   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_rl'                                           rigid lid
9   !!----------------------------------------------------------------------
10   !!   dyn_spg_rl   : update the momentum trend with the surface pressure
11   !!                  for the rigid-lid case.
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE phycst          ! physical constants
17   USE ldftra_oce      ! ocean active tracers: lateral physics
18   USE ldfdyn_oce      ! ocean dynamics: lateral physics
19   USE zdf_oce         ! ocean vertical physics
20   USE trddyn_oce      ! ocean dynamics trends
21   USE in_out_manager  ! I/O manager
22   USE sol_oce         ! ocean elliptic solver
23   USE solpcg          ! preconditionned conjugate gradient solver
24   USE solsor          ! Successive Over-relaxation solver
25   USE solfet          ! FETI solver
26   USE solisl          ! ???
27   USE obc_oce         ! Lateral open boundary condition
28   USE lib_mpp         ! ???
29   USE lbclnk          ! ???
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Accessibility
35   PUBLIC dyn_spg_rl   ! called by step.F90
36
37   !! * Shared module variables
38   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl = .TRUE.    !: rigid-lid flag
39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42#  include "vectopt_loop_substitute.h90"
43#  include "obc_vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !!   OPA 9.0 , LODYC-IPSL  (2003)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE dyn_spg_rl( kt, kindic )
51      !!----------------------------------------------------------------------
52      !!                  ***  routine dyn_spg_rl  ***
53      !!
54      !! ** Purpose :   Compute the now trend due to the surface pressure
55      !!      gradient for the rigid-lid case,  add it to the general trend of
56      !!      momentum equation.
57      !!
58      !! ** Method  :   Rigid-lid appromimation: the surface pressure gradient
59      !!      is given by:
60      !!         spgu = 1/rau0 d/dx(ps) = Mu + 1/(hu e2u) dj-1(bsfd)
61      !!         spgv = 1/rau0 d/dy(ps) = Mv - 1/(hv e1v) di-1(bsfd)
62      !!      where (Mu,Mv) is the vertically averaged momentum trend (i.e.
63      !!      the vertical ponderated sum of the general momentum trend),
64      !!      and bsfd is the barotropic streamfunction trend.
65      !!      The trend is computed as follows:
66      !!         -1- compute the vertically averaged momentum trend (Mu,Mv)
67      !!         -2- compute the barotropic streamfunction trend by solving an
68      !!      ellipic equation using a diagonal preconditioned conjugate
69      !!      gradient or a successive-over-relaxation method (depending
70      !!      on nsolv, a namelist parameter).
71      !!         -3- add to bsfd the island trends if lk_isl=T
72      !!         -4- compute the after streamfunction is for further diagnos-
73      !!      tics using a leap-frog scheme.
74      !!         -5- add the momentum trend associated with the surface pres-
75      !!      sure gradient to the general trend.
76      !!
77      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
78      !!             - Save the trends in (utrd,vtrd) ('key_diatrends')
79      !!
80      !! References :
81      !!      Madec et al. 1988, ocean modelling, issue 78, 1-6.
82      !!
83      !! History :
84      !!        !  96-05  (G. Madec, M. Imbard, M. Guyon)  rewitting in 1
85      !!                  routine, without pointers, and with the same matrix
86      !!                  for sor and pcg, mpp exchange, and symmetric conditions
87      !!        !  96-07  (A. Weaver)  Euler forward step
88      !!        !  96-11  (A. Weaver)  correction to preconditioning:
89      !!        !  98-02  (M. Guyon)  FETI method
90      !!        !  98-05  (G. Roullet)  free surface
91      !!        !  97-09  (J.-M. Molines)  Open boundaries
92      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
93      !!        !  02-11  (C. Talandier, A-M Treguier) Open boundaries
94      !!---------------------------------------------------------------------
95      !! * Arguments
96      INTEGER, INTENT( in  ) ::   kt       ! ocean time-step index
97      INTEGER, INTENT( out ) ::   kindic   ! solver flag, take a negative value
98      !                                    ! when the solver doesnot converge
99      !! * Local declarations
100      INTEGER ::   ji, jj, jk    ! dummy loop indices
101      REAL(wp) ::   zgwgt, zbsfa, zgcx, z2dt
102# if defined key_obc
103      INTEGER ::   ip, ii, ij
104      INTEGER ::   iii, ijj, jip, jnic
105      INTEGER ::   it, itm, itm2, ib, ibm, ibm2
106      REAL(wp) ::   z2dtr
107# endif
108      !!----------------------------------------------------------------------
109     
110      IF( kt == nit000 ) THEN
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) 'dyn_spg_rl : surface pressure gradient trend'
113         IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
114
115         ! set to zero rigid-lid specific arrays
116         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction)
117         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
118      ENDIF
119
120      ! 0. Initializations:
121      ! -------------------
122# if defined key_obc
123      ! space index on boundary arrays
124      ib = 1
125      ibm = 2
126      ibm2 = 3
127      ! time index on boundary arrays
128      it = 1
129      itm = 2
130      itm2 = 3
131# endif
132
133      !                                                ! ===============
134      DO jj = 2, jpjm1                                 !  Vertical slab
135         !                                             ! ===============
136
137         ! 1. Vertically averaged momentum trend
138         ! -------------------------------------
139         ! initialization to zero
140         spgu(:,jj) = 0.
141         spgv(:,jj) = 0.
142         ! vertical sum
143         DO jk = 1, jpk
144            DO ji = 2, jpim1
145               spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
146               spgv(ji,jj) = spgv(ji,jj) + va(ji,jj,jk) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)
147            END DO
148         END DO 
149         ! divide by the depth
150         spgu(:,jj) = spgu(:,jj) * hur(:,jj)
151         spgv(:,jj) = spgv(:,jj) * hvr(:,jj)
152
153         !                                             ! ===============
154      END DO                                           !   End of slab
155      !                                                ! ===============
156
157      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
158
159      ! Boundary conditions on (spgu,spgv)
160      CALL lbc_lnk( spgu, 'U', -1. )
161      CALL lbc_lnk( spgv, 'V', -1. )
162
163      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
164
165      ! 2. Barotropic streamfunction trend (bsfd)
166      ! ----------------------------------
167
168      ! Curl of the vertically averaged velocity
169      DO jj = 2, jpjm1
170         DO ji = 2, jpim1
171            gcb(ji,jj) = -gcdprc(ji,jj)   &
172                       *(  ( e2v(ji+1,jj  )*spgv(ji+1,jj  ) - e2v(ji,jj)*spgv(ji,jj) )   &
173                          -( e1u(ji  ,jj+1)*spgu(ji  ,jj+1) - e1u(ji,jj)*spgu(ji,jj) )   ) 
174         END DO
175      END DO
176
177# if defined key_obc
178      ! Open boundary contribution
179      DO jj = 2, jpjm1
180         DO ji = 2, jpim1
181           gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj)
182         END DO
183      END DO
184# else
185      ! No open boundary contribution
186# endif
187
188      ! First guess using previous solution of the elliptic system and
189      ! not bsfd since the system is solved with 0 as coastal boundary
190      ! condition. Also include a swap array (gcx,gxcb)
191      DO jj = 2, jpjm1
192         DO ji = 2, jpim1
193            zgcx        = gcx(ji,jj)
194            gcx (ji,jj) = 2.*zgcx - gcxb(ji,jj)
195            gcxb(ji,jj) =    zgcx
196         END DO
197      END DO
198
199      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
200
201      ! Relative precision (computation on one processor)
202      rnorme = 0.e0
203      DO jj = 1, jpj
204         DO ji = 1, jpi
205            zgwgt =  gcdmat(ji,jj) * gcb(ji,jj)
206            rnorme = rnorme + gcb(ji,jj) * zgwgt
207         END DO
208      END DO
209      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
210
211      epsr = eps*eps*rnorme
212      ncut = 0
213      ! if rnorme is 0, the solution is 0, the solver isn't called
214      IF( rnorme == 0.e0 ) THEN
215         bsfd (:,:) = 0.e0
216         res   = 0.e0
217         niter = 0
218         ncut  = 999
219      ENDIF
220
221      kindic = 0
222      ! solve the bsf system  ===> solution in gcx array
223      IF( ncut == 0 ) THEN
224         SELECT CASE ( nsolv )
225         CASE ( 1 )                    ! diagonal preconditioned conjuguate gradient
226            CALL sol_pcg( kindic )
227         CASE( 2 )                     ! successive-over-relaxation
228            CALL sol_sor( kindic )
229         CASE( 3 )                     ! FETI solver
230            CALL sol_fet( kindic )
231         CASE DEFAULT                  ! e r r o r in nsolv namelist parameter
232            IF(lwp) WRITE(numout,cform_err)
233            IF(lwp) WRITE(numout,*) ' dyn_spg_rl : e r r o r, nsolv = 1, 2 or 3'
234            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~                not = ', nsolv
235            nstop = nstop + 1
236         END SELECT
237      ENDIF
238
239
240      ! bsf trend update
241      ! ----------------
242
243      bsfd(:,:) = gcx(:,:)
244
245     
246      ! update bsf trend with islands trend
247      ! -----------------------------------
248
249      IF( lk_isl )   CALL isl_dyn_spg         ! update bsfd
250
251
252# if defined key_obc
253      ! Compute bsf trend for OBC points (not masked)
254
255      IF( lp_obc_east ) THEN
256      ! compute bsf trend at the boundary from bsfeob, computed in obc_spg
257      IF( neuler == 0 .AND. kt == nit000 ) THEN
258         z2dtr = 1. / rdt
259      ELSE
260         z2dtr = 1. / (2. * rdt )
261      ENDIF
262      ! (jped,jpefm1),nieob
263      DO ji = fs_nie0, fs_nie1   ! vector opt.
264         DO jj = nje0m1, nje1m1
265            bsfd(ji,jj) = ( bsfeob(jj) - bsfb(ji,jj) ) * z2dtr
266         END DO
267      END DO
268      ENDIF
269
270      IF( lp_obc_west ) THEN
271      ! compute bsf trend at the boundary from bsfwob, computed in obc_spg
272      IF( neuler == 0 .AND. kt == nit000 ) THEN
273         z2dtr = 1. / rdt
274      ELSE
275         z2dtr = 1. / ( 2. * rdt )
276      ENDIF
277      ! (jpwd,jpwfm1),niwob
278      DO ji = fs_niw0, fs_niw1   ! vector opt.
279         DO jj = njw0m1, njw1m1
280            bsfd(ji,jj) = ( bsfwob(jj) - bsfb(ji,jj) ) * z2dtr
281         END DO
282      END DO
283      ENDIF
284
285      IF( lp_obc_north ) THEN
286      ! compute bsf trend at the boundary from bsfnob, computed in obc_spg
287      IF( neuler == 0 .AND. kt == nit000 ) THEN
288         z2dtr = 1. / rdt
289      ELSE
290         z2dtr = 1. / ( 2. * rdt )
291      ENDIF
292      ! njnob,(jpnd,jpnfm1)
293      DO jj = fs_njn0, fs_njn1   ! vector opt.
294         DO ji = nin0m1, nin1m1
295            bsfd(ji,jj) = ( bsfnob(ji) - bsfb(ji,jj) ) * z2dtr
296         END DO
297      END DO
298      ENDIF
299
300      IF( lp_obc_south ) THEN
301      ! compute bsf trend at the boundary from bsfsob, computed in obc_spg
302      IF( neuler == 0 .AND. kt == nit000 ) THEN
303         z2dtr = 1. / rdt
304      ELSE
305         z2dtr = 1. / ( 2. * rdt )
306      ENDIF
307      ! njsob,(jpsd,jpsfm1)
308      DO jj = fs_njs0, fs_njs1   ! vector opt.
309         DO ji = nis0m1, nis1m1
310            bsfd(ji,jj) = ( bsfsob(ji) - bsfb(ji,jj) ) * z2dtr
311         END DO
312      END DO
313      ENDIF
314
315      ! compute bsf trend for isolated coastline points
316
317      IF( neuler == 0 .AND. kt == nit000 ) THEN
318         z2dtr = 1. / rdt
319      ELSE
320         z2dtr = 1. /( 2. * rdt )
321      ENDIF
322
323      IF( nbobc > 1 ) THEN
324         DO jnic = 1,nbobc - 1
325            ip = mnic(0,jnic)
326            DO jip = 1,ip
327               ii = miic(jip,0,jnic)
328               ij = mjic(jip,0,jnic)
329               IF( ii >= 1 + nimpp - 1 .AND. ii <= jpi + nimpp -1 .AND. &
330                   ij >= 1 + njmpp - 1 .AND. ij <= jpj + njmpp -1 ) THEN
331                  iii = ii - nimpp + 1
332                  ijj = ij - njmpp + 1
333                  bsfd(iii,ijj) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
334               ENDIF
335            END DO
336         END DO
337      ENDIF
338# endif
339
340      ! 4. Barotropic stream function and array swap
341      ! --------------------------------------------
342
343      ! Leap-frog time scheme, time filter and array swap
344      IF( neuler == 0 .AND. kt == nit000 ) THEN
345         ! Euler time stepping (first time step, starting from rest)
346         z2dt = rdt
347         DO jj = 1, jpj
348            DO ji = 1, jpi
349               zbsfa       = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
350               bsfb(ji,jj) = bsfn(ji,jj)
351               bsfn(ji,jj) = zbsfa 
352            END DO
353         END DO
354      ELSE
355         ! Leap-frog time stepping - Asselin filter
356         z2dt = 2.*rdt
357         DO jj = 1, jpj
358            DO ji = 1, jpi
359               zbsfa       = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
360               bsfb(ji,jj) = atfp * ( bsfb(ji,jj) + zbsfa ) + atfp1 * bsfn(ji,jj)
361               bsfn(ji,jj) = zbsfa
362            END DO
363         END DO
364      ENDIF
365
366# if defined key_obc
367      ! Swap of boundary arrays
368      IF( lp_obc_east ) THEN
369      ! (jped,jpef),nieob
370      IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
371         DO jj = nje0m1, nje1
372            ! fields itm2 <== itm
373            bebnd(jj,ib  ,itm2) = bebnd(jj,ib  ,itm)
374            bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
375            bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
376            bebnd(jj,ib  ,itm ) = bebnd(jj,ib  ,it )
377         END DO
378      ELSE
379         ! fields itm <== it  plus time filter at the boundary
380         DO ji = fs_nie0, fs_nie1   ! vector opt.
381            DO jj = nje0m1, nje1
382               bebnd(jj,ib  ,itm2) = bebnd(jj,ib  ,itm)
383               bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
384               bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
385               bebnd(jj,ib  ,itm ) = atfp * ( bebnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bebnd(jj,ib,it)
386               bebnd(jj,ibm ,itm ) = bebnd(jj,ibm ,it )
387               bebnd(jj,ibm2,itm ) = bebnd(jj,ibm2,it )
388            END DO
389         END DO
390      ENDIF
391      ! fields it <== now (kt+1)
392      DO ji = fs_nie0, fs_nie1   ! vector opt.
393         DO jj = nje0m1, nje1
394            bebnd(jj,ib  ,it  ) = bsfn (ji  ,jj)
395            bebnd(jj,ibm ,it  ) = bsfn (ji-1,jj)
396            bebnd(jj,ibm2,it  ) = bsfn (ji-2,jj)
397         END DO
398      END DO
399      IF( lk_mpp )   CALL mppobc( bebnd, jpjed, jpjef, jpieob, 3*3, 2, jpj )
400      ENDIF
401
402      IF( lp_obc_west ) THEN
403      ! (jpwd,jpwf),niwob
404      IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
405         DO jj = njw0m1, njw1
406            ! fields itm2 <== itm
407            bwbnd(jj,ib  ,itm2) = bwbnd(jj,ib  ,itm)
408            bwbnd(jj,ibm ,itm2) = bwbnd(jj,ibm ,itm)
409            bwbnd(jj,ibm2,itm2) = bwbnd(jj,ibm2,itm)
410            bwbnd(jj,ib  ,itm ) = bwbnd(jj,ib  ,it )
411         END DO
412      ELSE
413         DO ji = fs_niw0, fs_niw1   ! Vector opt.
414            DO jj = njw0m1, njw1
415               bwbnd(jj,ib  ,itm2) = bwbnd(jj,ib  ,itm)
416               bwbnd(jj,ibm ,itm2) = bwbnd(jj,ibm ,itm)
417               bwbnd(jj,ibm2,itm2) = bwbnd(jj,ibm2,itm)
418               ! fields itm <== it  plus time filter at the boundary
419               bwbnd(jj,ib  ,itm ) = atfp * ( bwbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bwbnd(jj,ib,it)
420               bwbnd(jj,ibm ,itm ) = bwbnd(jj,ibm ,it )
421               bwbnd(jj,ibm2,itm ) = bwbnd(jj,ibm2,it )
422            END DO
423         END DO
424      ENDIF
425      ! fields it <== now (kt+1)
426      DO ji = fs_niw0, fs_niw1   ! Vector opt.
427         DO jj = njw0m1, njw1
428            bwbnd(jj,ib  ,it  ) = bsfn (ji  ,jj)
429            bwbnd(jj,ibm ,it  ) = bsfn (ji+1,jj)
430            bwbnd(jj,ibm2,it  ) = bsfn (ji+2,jj)
431         END DO
432      END DO
433      IF( lk_mpp )   CALL mppobc( bwbnd, jpjwd, jpjwf, jpiwob, 3*3, 2, jpj )
434      ENDIF
435
436      IF( lp_obc_north ) THEN
437   ! njnob,(jpnd,jpnf)
438      IF( kt < nit000 + 3 .AND. .NOT.ln_rstart ) THEN
439         DO ji = nin0m1, nin1
440            ! fields itm2 <== itm
441            bnbnd(ji,ib  ,itm2) = bnbnd(ji,ib  ,itm)
442            bnbnd(ji,ibm ,itm2) = bnbnd(ji,ibm ,itm)
443            bnbnd(ji,ibm2,itm2) = bnbnd(ji,ibm2,itm)
444            bnbnd(ji,ib  ,itm ) = bnbnd(ji,ib  ,it )
445         END DO
446      ELSE
447         DO jj = fs_njn0, fs_njn1   ! Vector opt.
448            DO ji = nin0m1, nin1
449               bnbnd(ji,ib  ,itm2) = bnbnd(ji,ib  ,itm)
450               bnbnd(ji,ibm ,itm2) = bnbnd(ji,ibm ,itm)
451               bnbnd(ji,ibm2,itm2) = bnbnd(ji,ibm2,itm)
452               ! fields itm <== it  plus time filter at the boundary
453               bnbnd(jj,ib  ,itm ) = atfp * ( bnbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bnbnd(jj,ib,it)
454               bnbnd(ji,ibm ,itm ) = bnbnd(ji,ibm ,it )
455               bnbnd(ji,ibm2,itm ) = bnbnd(ji,ibm2,it )
456            END DO
457          END DO
458      ENDIF
459      ! fields it <== now (kt+1)
460      DO jj = fs_njn0, fs_njn1   ! Vector opt.
461         DO ji = nin0m1, nin1
462            bnbnd(ji,ib  ,it  ) = bsfn (ji,jj  )
463            bnbnd(ji,ibm ,it  ) = bsfn (ji,jj-1)
464            bnbnd(ji,ibm2,it  ) = bsfn (ji,jj-2)
465         END DO
466      END DO
467      IF( lk_mpp )   CALL mppobc( bnbnd, jpind, jpinf, jpjnob, 3*3, 1, jpi )
468      ENDIF
469
470      IF( lp_obc_south ) THEN
471         ! njsob,(jpsd,jpsf)
472         IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
473            DO ji = nis0m1, nis1
474               ! fields itm2 <== itm
475               bsbnd(ji,ib  ,itm2) = bsbnd(ji,ib  ,itm)
476               bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
477               bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
478               bsbnd(ji,ib  ,itm ) = bsbnd(ji,ib  ,it )
479            END DO
480         ELSE
481            DO jj = fs_njs0, fs_njs1   ! vector opt.
482               DO ji = nis0m1, nis1
483                  bsbnd(ji,ib  ,itm2) = bsbnd(ji,ib  ,itm)
484                  bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
485                  bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
486                  ! fields itm <== it  plus time filter at the boundary
487                  bsbnd(jj,ib  ,itm ) = atfp * ( bsbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bsbnd(jj,ib,it)
488                  bsbnd(ji,ibm ,itm ) = bsbnd(ji,ibm ,it )
489                  bsbnd(ji,ibm2,itm ) = bsbnd(ji,ibm2,it )
490               END DO
491            END DO
492         ENDIF
493         DO jj = fs_njs0, fs_njs1   ! vector opt.
494            DO ji = nis0m1, nis1 
495               ! fields it <== now (kt+1)
496               bsbnd(ji,ib  ,it  ) = bsfn (ji,jj  )
497               bsbnd(ji,ibm ,it  ) = bsfn (ji,jj+1)
498               bsbnd(ji,ibm2,it  ) = bsfn (ji,jj+2)
499            END DO
500         END DO
501         IF( lk_mpp )   CALL mppobc( bsbnd, jpisd, jpisf, jpjsob, 3*3, 1, jpi )
502      ENDIF
503# endif
504      !
505      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
506     
507      !  add the surface pressure trend to the general trend
508      ! -----------------------------------------------------
509     
510      DO jj = 2, jpjm1
511
512         ! update the surface pressure gradient with the barotropic trend
513         DO ji = 2, jpim1
514            spgu(ji,jj) = spgu(ji,jj) + hur(ji,jj) / e2u(ji,jj) * ( bsfd(ji,jj) - bsfd(ji  ,jj-1) )
515            spgv(ji,jj) = spgv(ji,jj) - hvr(ji,jj) / e1v(ji,jj) * ( bsfd(ji,jj) - bsfd(ji-1,jj  ) )
516         END DO
517         ! add the surface pressure gradient trend to the general trend
518         DO jk = 1, jpkm1
519            DO ji = 2, jpim1
520               ua(ji,jj,jk) = ua(ji,jj,jk) - spgu(ji,jj)
521               va(ji,jj,jk) = va(ji,jj,jk) - spgv(ji,jj)
522# if defined key_trddyn || defined key_trd_vor
523               ! save the surface pressure gradient trend for diagnostics
524               utrd(ji,jj,jk,8) = -spgu(ji,jj)
525               vtrd(ji,jj,jk,8) = -spgv(ji,jj)
526# endif
527            END DO
528         END DO
529
530      END DO
531
532   END SUBROUTINE dyn_spg_rl
533
534#else
535   !!----------------------------------------------------------------------
536   !!   'key_dynspg_rl'                                        NO rigid lid
537   !!----------------------------------------------------------------------
538   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl = .FALSE.   !: rigid-lid flag
539CONTAINS
540   SUBROUTINE dyn_spg_rl( kt, kindic )       ! Empty routine
541      WRITE(*,*) 'dyn_spg_rl: You should not have seen this print! error?', kt, kindic
542   END SUBROUTINE dyn_spg_rl
543#endif
544   
545   !!======================================================================
546END MODULE dynspg_rl
Note: See TracBrowser for help on using the repository browser.