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

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

CL : Add CVS Header and CeCILL licence information

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