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

Last change on this file since 474 was 474, checked in by opalod, 15 years ago

nemo_v1_update_061: SM: end of ctl_stop + mpi optimization in _bilap

  • 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 sol_oce         ! ocean elliptic solver
21   USE solpcg          ! preconditionned conjugate gradient solver
22   USE solsor          ! Successive Over-relaxation solver
23   USE solfet          ! FETI solver
24   USE solsor_e        ! Successive Over-relaxation solver with MPP optimization
25   USE solisl          ! ???
26   USE obc_oce         ! Lateral open boundary condition
27   USE lib_mpp         ! distributed memory computing library
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29   USE in_out_manager  ! I/O manager
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Accessibility
35   PUBLIC dyn_spg_rl   ! called by step.F90
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40#  include "obc_vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !!   OPA 9.0 , LOCEAN-IPSL (2005)
43   !! $Header$
44   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE dyn_spg_rl( kt, kindic )
50      !!----------------------------------------------------------------------
51      !!                  ***  routine dyn_spg_rl  ***
52      !!
53      !! ** Purpose :   Compute the now trend due to the surface pressure
54      !!      gradient for the rigid-lid case,  add it to the general trend of
55      !!      momentum equation.
56      !!
57      !! ** Method  :   Rigid-lid appromimation: the surface pressure gradient
58      !!      is given by:
59      !!         spgu = 1/rau0 d/dx(ps) = Mu + 1/(hu e2u) dj-1(bsfd)
60      !!         spgv = 1/rau0 d/dy(ps) = Mv - 1/(hv e1v) di-1(bsfd)
61      !!      where (Mu,Mv) is the vertically averaged momentum trend (i.e.
62      !!      the vertical ponderated sum of the general momentum trend),
63      !!      and bsfd is the barotropic streamfunction trend.
64      !!      The trend is computed as follows:
65      !!         -1- compute the vertically averaged momentum trend (Mu,Mv)
66      !!         -2- compute the barotropic streamfunction trend by solving an
67      !!      ellipic equation using a diagonal preconditioned conjugate
68      !!      gradient or a successive-over-relaxation method (depending
69      !!      on nsolv, a namelist parameter).
70      !!         -3- add to bsfd the island trends if lk_isl=T
71      !!         -4- compute the after streamfunction is for further diagnos-
72      !!      tics using a leap-frog scheme.
73      !!         -5- add the momentum trend associated with the surface pres-
74      !!      sure gradient to the general trend.
75      !!
76      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
77      !!
78      !! References :
79      !!      Madec et al. 1988, ocean modelling, issue 78, 1-6.
80      !!
81      !! History :
82      !!        !  96-05  (G. Madec, M. Imbard, M. Guyon)  rewitting in 1
83      !!                  routine, without pointers, and with the same matrix
84      !!                  for sor and pcg, mpp exchange, and symmetric conditions
85      !!        !  96-07  (A. Weaver)  Euler forward step
86      !!        !  96-11  (A. Weaver)  correction to preconditioning:
87      !!        !  98-02  (M. Guyon)  FETI method
88      !!        !  98-05  (G. Roullet)  free surface
89      !!        !  97-09  (J.-M. Molines)  Open boundaries
90      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
91      !!        !  02-11  (C. Talandier, A-M Treguier) Open boundaries
92      !!   9.0  !  04-08  (C. Talandier)  New trends organization
93      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
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) ::  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      ! applied the lateral boundary conditions
199      IF( nsolv == 4)   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
200
201      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
202
203      ! Relative precision (computation on one processor)
204      rnorme = 0.e0
205      rnorme = SUM( gcb(1:nlci,1:nlcj) * gcdmat(1:nlci,1:nlcj) * gcb(1:nlci,1:nlcj) * bmask(:,:) )
206      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
207
208      epsr = eps*eps*rnorme
209      ncut = 0
210      ! if rnorme is 0, the solution is 0, the solver isn't called
211      IF( rnorme == 0.e0 ) THEN
212         bsfd (:,:) = 0.e0
213         res   = 0.e0
214         niter = 0
215         ncut  = 999
216      ENDIF
217
218      kindic = 0
219      ! solve the bsf system  ===> solution in gcx array
220      IF( ncut == 0 ) THEN
221         SELECT CASE ( nsolv )
222         CASE ( 1 )                    ! diagonal preconditioned conjuguate gradient
223            CALL sol_pcg( kindic )
224         CASE( 2 )                     ! successive-over-relaxation
225            CALL sol_sor( kindic )
226         CASE( 3 )                     ! FETI solver
227            CALL sol_fet( kindic )
228         CASE( 4 )                     ! successive-over-relaxation with extra outer halo
229            CALL sol_sor_e( kindic )
230         CASE DEFAULT                  ! e r r o r in nsolv namelist parameter
231            WRITE(ctmp1,*) ' ~~~~~~~~~~                not = ', nsolv
232            CALL ctl_stop( ' dyn_spg_rl : e r r o r, nsolv = 1, 2 ,3 or 4', ctmp1 )
233         END SELECT
234      ENDIF
235
236
237      ! bsf trend update
238      ! ----------------
239
240      bsfd(1:nlci,1:nlcj) = gcx(1:nlci,1:nlcj)
241
242     
243      ! update bsf trend with islands trend
244      ! -----------------------------------
245
246      IF( lk_isl )   CALL isl_dyn_spg         ! update bsfd
247
248
249# if defined key_obc
250      ! Compute bsf trend for OBC points (not masked)
251
252      IF( lp_obc_east ) THEN
253      ! compute bsf trend at the boundary from bsfeob, computed in obc_spg
254      IF( neuler == 0 .AND. kt == nit000 ) THEN
255         z2dtr = 1. / rdt
256      ELSE
257         z2dtr = 1. / (2. * rdt )
258      ENDIF
259      ! (jped,jpefm1),nieob
260      DO ji = fs_nie0, fs_nie1   ! vector opt.
261         DO jj = nje0m1, nje1m1
262            bsfd(ji,jj) = ( bsfeob(jj) - bsfb(ji,jj) ) * z2dtr
263         END DO
264      END DO
265      ENDIF
266
267      IF( lp_obc_west ) THEN
268      ! compute bsf trend at the boundary from bsfwob, computed in obc_spg
269      IF( neuler == 0 .AND. kt == nit000 ) THEN
270         z2dtr = 1. / rdt
271      ELSE
272         z2dtr = 1. / ( 2. * rdt )
273      ENDIF
274      ! (jpwd,jpwfm1),niwob
275      DO ji = fs_niw0, fs_niw1   ! vector opt.
276         DO jj = njw0m1, njw1m1
277            bsfd(ji,jj) = ( bsfwob(jj) - bsfb(ji,jj) ) * z2dtr
278         END DO
279      END DO
280      ENDIF
281
282      IF( lp_obc_north ) THEN
283      ! compute bsf trend at the boundary from bsfnob, computed in obc_spg
284      IF( neuler == 0 .AND. kt == nit000 ) THEN
285         z2dtr = 1. / rdt
286      ELSE
287         z2dtr = 1. / ( 2. * rdt )
288      ENDIF
289      ! njnob,(jpnd,jpnfm1)
290      DO jj = fs_njn0, fs_njn1   ! vector opt.
291         DO ji = nin0m1, nin1m1
292            bsfd(ji,jj) = ( bsfnob(ji) - bsfb(ji,jj) ) * z2dtr
293         END DO
294      END DO
295      ENDIF
296
297      IF( lp_obc_south ) THEN
298      ! compute bsf trend at the boundary from bsfsob, computed in obc_spg
299      IF( neuler == 0 .AND. kt == nit000 ) THEN
300         z2dtr = 1. / rdt
301      ELSE
302         z2dtr = 1. / ( 2. * rdt )
303      ENDIF
304      ! njsob,(jpsd,jpsfm1)
305      DO jj = fs_njs0, fs_njs1   ! vector opt.
306         DO ji = nis0m1, nis1m1
307            bsfd(ji,jj) = ( bsfsob(ji) - bsfb(ji,jj) ) * z2dtr
308         END DO
309      END DO
310      ENDIF
311
312      ! compute bsf trend for isolated coastline points
313
314      IF( neuler == 0 .AND. kt == nit000 ) THEN
315         z2dtr = 1. / rdt
316      ELSE
317         z2dtr = 1. /( 2. * rdt )
318      ENDIF
319
320      IF( nbobc > 1 ) THEN
321         DO jnic = 1,nbobc - 1
322            ip = mnic(0,jnic)
323            DO jip = 1,ip
324               ii = miic(jip,0,jnic)
325               ij = mjic(jip,0,jnic)
326               IF( ii >= 1 + nimpp - 1 .AND. ii <= jpi + nimpp -1 .AND. &
327                   ij >= 1 + njmpp - 1 .AND. ij <= jpj + njmpp -1 ) THEN
328                  iii = ii - nimpp + 1
329                  ijj = ij - njmpp + 1
330                  bsfd(iii,ijj) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
331               ENDIF
332            END DO
333         END DO
334      ENDIF
335# endif
336
337      ! 4. Barotropic stream function and array swap
338      ! --------------------------------------------
339
340      ! Leap-frog time scheme, time filter and array swap
341      IF( neuler == 0 .AND. kt == nit000 ) THEN
342         ! Euler time stepping (first time step, starting from rest)
343         z2dt = rdt
344         DO jj = 1, jpj
345            DO ji = 1, jpi
346               zbsfa       = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
347               bsfb(ji,jj) = bsfn(ji,jj)
348               bsfn(ji,jj) = zbsfa 
349            END DO
350         END DO
351      ELSE
352         ! Leap-frog time stepping - Asselin filter
353         z2dt = 2.*rdt
354         DO jj = 1, jpj
355            DO ji = 1, jpi
356               zbsfa       = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
357               bsfb(ji,jj) = atfp * ( bsfb(ji,jj) + zbsfa ) + atfp1 * bsfn(ji,jj)
358               bsfn(ji,jj) = zbsfa
359            END DO
360         END DO
361      ENDIF
362
363# if defined key_obc
364      ! Swap of boundary arrays
365      IF( lp_obc_east ) THEN
366      ! (jped,jpef),nieob
367      IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
368         DO jj = nje0m1, nje1
369            ! fields itm2 <== itm
370            bebnd(jj,ib  ,itm2) = bebnd(jj,ib  ,itm)
371            bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
372            bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
373            bebnd(jj,ib  ,itm ) = bebnd(jj,ib  ,it )
374         END DO
375      ELSE
376         ! fields itm <== it  plus time filter at the boundary
377         DO ji = fs_nie0, fs_nie1   ! vector opt.
378            DO jj = nje0m1, nje1
379               bebnd(jj,ib  ,itm2) = bebnd(jj,ib  ,itm)
380               bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
381               bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
382               bebnd(jj,ib  ,itm ) = atfp * ( bebnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bebnd(jj,ib,it)
383               bebnd(jj,ibm ,itm ) = bebnd(jj,ibm ,it )
384               bebnd(jj,ibm2,itm ) = bebnd(jj,ibm2,it )
385            END DO
386         END DO
387      ENDIF
388      ! fields it <== now (kt+1)
389      DO ji = fs_nie0, fs_nie1   ! vector opt.
390         DO jj = nje0m1, nje1
391            bebnd(jj,ib  ,it  ) = bsfn (ji  ,jj)
392            bebnd(jj,ibm ,it  ) = bsfn (ji-1,jj)
393            bebnd(jj,ibm2,it  ) = bsfn (ji-2,jj)
394         END DO
395      END DO
396      IF( lk_mpp )   CALL mppobc( bebnd, jpjed, jpjef, jpieob, 3*3, 2, jpj )
397      ENDIF
398
399      IF( lp_obc_west ) THEN
400      ! (jpwd,jpwf),niwob
401      IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
402         DO jj = njw0m1, njw1
403            ! fields itm2 <== itm
404            bwbnd(jj,ib  ,itm2) = bwbnd(jj,ib  ,itm)
405            bwbnd(jj,ibm ,itm2) = bwbnd(jj,ibm ,itm)
406            bwbnd(jj,ibm2,itm2) = bwbnd(jj,ibm2,itm)
407            bwbnd(jj,ib  ,itm ) = bwbnd(jj,ib  ,it )
408         END DO
409      ELSE
410         DO ji = fs_niw0, fs_niw1   ! Vector opt.
411            DO jj = njw0m1, njw1
412               bwbnd(jj,ib  ,itm2) = bwbnd(jj,ib  ,itm)
413               bwbnd(jj,ibm ,itm2) = bwbnd(jj,ibm ,itm)
414               bwbnd(jj,ibm2,itm2) = bwbnd(jj,ibm2,itm)
415               ! fields itm <== it  plus time filter at the boundary
416               bwbnd(jj,ib  ,itm ) = atfp * ( bwbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bwbnd(jj,ib,it)
417               bwbnd(jj,ibm ,itm ) = bwbnd(jj,ibm ,it )
418               bwbnd(jj,ibm2,itm ) = bwbnd(jj,ibm2,it )
419            END DO
420         END DO
421      ENDIF
422      ! fields it <== now (kt+1)
423      DO ji = fs_niw0, fs_niw1   ! Vector opt.
424         DO jj = njw0m1, njw1
425            bwbnd(jj,ib  ,it  ) = bsfn (ji  ,jj)
426            bwbnd(jj,ibm ,it  ) = bsfn (ji+1,jj)
427            bwbnd(jj,ibm2,it  ) = bsfn (ji+2,jj)
428         END DO
429      END DO
430      IF( lk_mpp )   CALL mppobc( bwbnd, jpjwd, jpjwf, jpiwob, 3*3, 2, jpj )
431      ENDIF
432
433      IF( lp_obc_north ) THEN
434   ! njnob,(jpnd,jpnf)
435      IF( kt < nit000 + 3 .AND. .NOT.ln_rstart ) THEN
436         DO ji = nin0m1, nin1
437            ! fields itm2 <== itm
438            bnbnd(ji,ib  ,itm2) = bnbnd(ji,ib  ,itm)
439            bnbnd(ji,ibm ,itm2) = bnbnd(ji,ibm ,itm)
440            bnbnd(ji,ibm2,itm2) = bnbnd(ji,ibm2,itm)
441            bnbnd(ji,ib  ,itm ) = bnbnd(ji,ib  ,it )
442         END DO
443      ELSE
444         DO jj = fs_njn0, fs_njn1   ! Vector opt.
445            DO ji = nin0m1, nin1
446               bnbnd(ji,ib  ,itm2) = bnbnd(ji,ib  ,itm)
447               bnbnd(ji,ibm ,itm2) = bnbnd(ji,ibm ,itm)
448               bnbnd(ji,ibm2,itm2) = bnbnd(ji,ibm2,itm)
449               ! fields itm <== it  plus time filter at the boundary
450               bnbnd(jj,ib  ,itm ) = atfp * ( bnbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bnbnd(jj,ib,it)
451               bnbnd(ji,ibm ,itm ) = bnbnd(ji,ibm ,it )
452               bnbnd(ji,ibm2,itm ) = bnbnd(ji,ibm2,it )
453            END DO
454          END DO
455      ENDIF
456      ! fields it <== now (kt+1)
457      DO jj = fs_njn0, fs_njn1   ! Vector opt.
458         DO ji = nin0m1, nin1
459            bnbnd(ji,ib  ,it  ) = bsfn (ji,jj  )
460            bnbnd(ji,ibm ,it  ) = bsfn (ji,jj-1)
461            bnbnd(ji,ibm2,it  ) = bsfn (ji,jj-2)
462         END DO
463      END DO
464      IF( lk_mpp )   CALL mppobc( bnbnd, jpind, jpinf, jpjnob, 3*3, 1, jpi )
465      ENDIF
466
467      IF( lp_obc_south ) THEN
468         ! njsob,(jpsd,jpsf)
469         IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
470            DO ji = nis0m1, nis1
471               ! fields itm2 <== itm
472               bsbnd(ji,ib  ,itm2) = bsbnd(ji,ib  ,itm)
473               bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
474               bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
475               bsbnd(ji,ib  ,itm ) = bsbnd(ji,ib  ,it )
476            END DO
477         ELSE
478            DO jj = fs_njs0, fs_njs1   ! vector opt.
479               DO ji = nis0m1, nis1
480                  bsbnd(ji,ib  ,itm2) = bsbnd(ji,ib  ,itm)
481                  bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
482                  bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
483                  ! fields itm <== it  plus time filter at the boundary
484                  bsbnd(jj,ib  ,itm ) = atfp * ( bsbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bsbnd(jj,ib,it)
485                  bsbnd(ji,ibm ,itm ) = bsbnd(ji,ibm ,it )
486                  bsbnd(ji,ibm2,itm ) = bsbnd(ji,ibm2,it )
487               END DO
488            END DO
489         ENDIF
490         DO jj = fs_njs0, fs_njs1   ! vector opt.
491            DO ji = nis0m1, nis1 
492               ! fields it <== now (kt+1)
493               bsbnd(ji,ib  ,it  ) = bsfn (ji,jj  )
494               bsbnd(ji,ibm ,it  ) = bsfn (ji,jj+1)
495               bsbnd(ji,ibm2,it  ) = bsfn (ji,jj+2)
496            END DO
497         END DO
498         IF( lk_mpp )   CALL mppobc( bsbnd, jpisd, jpisf, jpjsob, 3*3, 1, jpi )
499      ENDIF
500# endif
501      !
502      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
503     
504      !  add the surface pressure trend to the general trend
505      ! -----------------------------------------------------
506     
507      DO jj = 2, jpjm1
508
509         ! update the surface pressure gradient with the barotropic trend
510         DO ji = 2, jpim1
511            spgu(ji,jj) = spgu(ji,jj) + hur(ji,jj) / e2u(ji,jj) * ( bsfd(ji,jj) - bsfd(ji  ,jj-1) )
512            spgv(ji,jj) = spgv(ji,jj) - hvr(ji,jj) / e1v(ji,jj) * ( bsfd(ji,jj) - bsfd(ji-1,jj  ) )
513         END DO
514         ! add the surface pressure gradient trend to the general trend
515         DO jk = 1, jpkm1
516            DO ji = 2, jpim1
517               ua(ji,jj,jk) = ua(ji,jj,jk) - spgu(ji,jj)
518               va(ji,jj,jk) = va(ji,jj,jk) - spgv(ji,jj)
519            END DO
520         END DO
521
522      END DO
523
524   END SUBROUTINE dyn_spg_rl
525
526#else
527   !!----------------------------------------------------------------------
528   !!   'key_dynspg_rl'                                        NO rigid lid
529   !!----------------------------------------------------------------------
530CONTAINS
531   SUBROUTINE dyn_spg_rl( kt, kindic )       ! Empty routine
532      WRITE(*,*) 'dyn_spg_rl: You should not have seen this print! error?', kt, kindic
533   END SUBROUTINE dyn_spg_rl
534#endif
535   
536   !!======================================================================
537END MODULE dynspg_rl
Note: See TracBrowser for help on using the repository browser.