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

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

CT : UPDATE151 : New trends organization

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