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

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

Initial revision

  • 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   !!----------------------------------------------------------------------
44   !!   OPA 9.0 , LODYC-IPSL  (2003)
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 l_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      !!             - Save the trends in (utrd,vtrd) ('key_diatrends')
78      !!
79      !! References :
80      !!      Madec et al. 1988, ocean modelling, issue 78, 1-6.
81      !!
82      !! History :
83      !!        !  96-05  (G. Madec, M. Imbard, M. Guyon)  rewitting in 1
84      !!                  routine, without pointers, and with the same matrix
85      !!                  for sor and pcg, mpp exchange, and symmetric conditions
86      !!        !  96-07  (A. Weaver)  Euler forward step
87      !!        !  96-11  (A. Weaver)  correction to preconditioning:
88      !!        !  98-02  (M. Guyon)  FETI method
89      !!        !  98-05  (G. Roullet)  free surface
90      !!        !  97-09  (J.-M. Molines)  Open boundaries
91      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
92      !!        !  02-11  (C. Talandier, A-M Treguier) Open boundaries
93      !!---------------------------------------------------------------------
94      !! * Arguments
95      INTEGER, INTENT( in  ) ::   kt       ! ocean time-step index
96      INTEGER, INTENT( out ) ::   kindic   ! solver flag, take a negative value
97      !                                    ! when the solver doesnot converge
98      !! * Local declarations
99      INTEGER ::   ji, jj, jk    ! dummy loop indices
100      REAL(wp) ::   zgwgt, zbsfa, zgcx, z2dt
101# if defined key_obc
102      INTEGER ::   ip, ii, ij
103      INTEGER ::   iii, ijj, jip, jnic
104      INTEGER ::   it, itm, itm2, ib, ibm, ibm2
105      REAL(wp) ::   z2dtr
106# endif
107      !!----------------------------------------------------------------------
108     
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'dyn_spg_rl : surface pressure gradient trend'
112         IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
113
114         ! set to zero rigid-lid specific arrays
115         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction)
116         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
117         bsfb(:,:) = 0.e0      ! before barotropic stream-function
118         bsfn(:,:) = 0.e0      ! now    barotropic stream-function
119         bsfd(:,:) = 0.e0      ! barotropic stream-function trend
120      ENDIF
121
122      ! 0. Initializations:
123      ! -------------------
124# if defined key_obc
125      ! space index on boundary arrays
126      ib = 1
127      ibm = 2
128      ibm2 = 3
129      ! time index on boundary arrays
130      it = 1
131      itm = 2
132      itm2 = 3
133# endif
134
135      !                                                ! ===============
136      DO jj = 2, jpjm1                                 !  Vertical slab
137         !                                             ! ===============
138
139         ! 1. Vertically averaged momentum trend
140         ! -------------------------------------
141         ! initialization to zero
142         spgu(:,jj) = 0.
143         spgv(:,jj) = 0.
144         ! vertical sum
145         DO jk = 1, jpk
146            DO ji = 2, jpim1
147               spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
148               spgv(ji,jj) = spgv(ji,jj) + va(ji,jj,jk) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)
149            END DO
150         END DO 
151         ! divide by the depth
152         spgu(:,jj) = spgu(:,jj) * hur(:,jj)
153         spgv(:,jj) = spgv(:,jj) * hvr(:,jj)
154
155         !                                             ! ===============
156      END DO                                           !   End of slab
157      !                                                ! ===============
158
159      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
160
161      ! Boundary conditions on (spgu,spgv)
162      CALL lbc_lnk( spgu, 'U', -1. )
163      CALL lbc_lnk( spgv, 'V', -1. )
164
165      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
166
167      ! 2. Barotropic streamfunction trend (bsfd)
168      ! ----------------------------------
169
170      ! Curl of the vertically averaged velocity
171      DO jj = 2, jpjm1
172         DO ji = 2, jpim1
173            gcb(ji,jj) = -gcdprc(ji,jj)   &
174                       *(  ( e2v(ji+1,jj  )*spgv(ji+1,jj  ) - e2v(ji,jj)*spgv(ji,jj) )   &
175                          -( e1u(ji  ,jj+1)*spgu(ji  ,jj+1) - e1u(ji,jj)*spgu(ji,jj) )   ) 
176         END DO
177      END DO
178
179# if defined key_obc
180      ! Open boundary contribution
181      DO jj = 2, jpjm1
182         DO ji = 2, jpim1
183           gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj)
184         END DO
185      END DO
186# else
187      ! No open boundary contribution
188# endif
189
190      ! First guess using previous solution of the elliptic system and
191      ! not bsfd since the system is solved with 0 as coastal boundary
192      ! condition. Also include a swap array (gcx,gxcb)
193      DO jj = 2, jpjm1
194         DO ji = 2, jpim1
195            zgcx        = gcx(ji,jj)
196            gcx (ji,jj) = 2.*zgcx - gcxb(ji,jj)
197            gcxb(ji,jj) =    zgcx
198         END DO
199      END DO
200
201      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
202
203      ! Relative precision (computation on one processor)
204      rnorme = 0.e0
205      DO jj = 1, jpj
206         DO ji = 1, jpi
207            zgwgt =  gcdmat(ji,jj) * gcb(ji,jj)
208            rnorme = rnorme + gcb(ji,jj) * zgwgt
209         END DO
210      END DO
211# if defined key_mpp
212      CALL mpp_sum( rnorme )
213# endif
214      epsr = eps*eps*rnorme
215      ncut = 0
216      ! if rnorme is 0, the solution is 0, the solver isn't called
217      IF( rnorme == 0.e0 ) THEN
218         bsfd (:,:) = 0.e0
219         res   = 0.e0
220         niter = 0
221         ncut  = 999
222      ENDIF
223
224      kindic = 0
225      ! solve the bsf system  ===> solution in gcx array
226      IF( ncut == 0 ) THEN
227         SELECT CASE ( nsolv )
228         CASE ( 1 )                    ! diagonal preconditioned conjuguate gradient
229            CALL sol_pcg( kindic )
230         CASE( 2 )                     ! successive-over-relaxation
231            CALL sol_sor( kt, kindic )
232         CASE( 3 )                     ! FETI solver
233            CALL sol_fet( kindic )
234         CASE DEFAULT                  ! e r r o r in nsolv namelist parameter
235            IF(lwp) WRITE(numout,cform_err)
236            IF(lwp) WRITE(numout,*) ' dyn_spg_rl : e r r o r, nsolv = 1, 2 or 3'
237            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~                not = ', nsolv
238            nstop = nstop + 1
239         END SELECT
240      ENDIF
241
242
243      ! bsf trend update
244      ! ----------------
245
246      bsfd(:,:) = gcx(:,:)
247
248     
249      ! update bsf trend with islands trend
250      ! -----------------------------------
251
252      IF( l_isl )   CALL isl_dyn_spg         ! update bsfd
253
254
255# if defined key_obc
256      ! Compute bsf trend for OBC points (not masked)
257
258      IF( lpeastobc ) THEN
259      ! compute bsf trend at the boundary from bsfeob, computed in obc_spg
260      IF( neuler == 0 .AND. kt == nit000 ) THEN
261         z2dtr = 1. / rdt
262      ELSE
263         z2dtr = 1. / (2. * rdt )
264      ENDIF
265      ! (jped,jpefm1),nieob
266      DO ji = fs_nie0, fs_nie1   ! vector opt.
267         DO jj = nje0m1, nje1m1
268            bsfd(ji,jj) = ( bsfeob(jj) - bsfb(ji,jj) ) * z2dtr
269         END DO
270      END DO
271      ENDIF
272
273      IF( lpwestobc ) THEN
274      ! compute bsf trend at the boundary from bsfwob, computed in obc_spg
275      IF( neuler == 0 .AND. kt == nit000 ) THEN
276         z2dtr = 1. / rdt
277      ELSE
278         z2dtr = 1. / ( 2. * rdt )
279      ENDIF
280      ! (jpwd,jpwfm1),niwob
281      DO ji = fs_niw0, fs_niw1   ! vector opt.
282         DO jj = njw0m1, njw1m1
283            bsfd(ji,jj) = ( bsfwob(jj) - bsfb(ji,jj) ) * z2dtr
284         END DO
285      END DO
286      ENDIF
287
288      IF( lpnorthobc ) THEN
289      ! compute bsf trend at the boundary from bsfnob, computed in obc_spg
290      IF( neuler == 0 .AND. kt == nit000 ) THEN
291         z2dtr = 1. / rdt
292      ELSE
293         z2dtr = 1. / ( 2. * rdt )
294      ENDIF
295      ! njnob,(jpnd,jpnfm1)
296      DO jj = fs_njn0, fs_njn1   ! vector opt.
297         DO ji = nin0m1, nin1m1
298            bsfd(ji,jj) = ( bsfnob(ji) - bsfb(ji,jj) ) * z2dtr
299         END DO
300      END DO
301      ENDIF
302
303      IF( lpsouthobc ) THEN
304      ! compute bsf trend at the boundary from bsfsob, computed in obc_spg
305      IF( neuler == 0 .AND. kt == nit000 ) THEN
306         z2dtr = 1. / rdt
307      ELSE
308         z2dtr = 1. / ( 2. * rdt )
309      ENDIF
310      ! njsob,(jpsd,jpsfm1)
311      DO jj = fs_njs0, fs_njs1   ! vector opt.
312         DO ji = nis0m1, nis1m1
313            bsfd(ji,jj) = ( bsfsob(ji) - bsfb(ji,jj) ) * z2dtr
314         END DO
315      END DO
316      ENDIF
317
318      ! compute bsf trend for isolated coastline points
319
320      IF( neuler == 0 .AND. kt == nit000 ) THEN
321         z2dtr = 1. / rdt
322      ELSE
323         z2dtr = 1. /( 2. * rdt )
324      ENDIF
325
326      IF( nbobc > 1 ) THEN
327         DO jnic = 1,nbobc - 1
328            ip = mnic(0,jnic)
329            DO jip = 1,ip
330               ii = miic(jip,0,jnic)
331               ij = mjic(jip,0,jnic)
332               IF( ii >= 1 + nimpp - 1 .AND. ii <= jpi + nimpp -1 .AND. &
333                   ij >= 1 + njmpp - 1 .AND. ij <= jpj + njmpp -1 ) THEN
334                  iii = ii - nimpp + 1
335                  ijj = ij - njmpp + 1
336                  bsfd(iii,ijj) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
337               ENDIF
338            END DO
339         END DO
340      ENDIF
341# endif
342
343      ! 4. Barotropic stream function and array swap
344      ! --------------------------------------------
345
346      ! Leap-frog time scheme, time filter and array swap
347      IF( neuler == 0 .AND. kt == nit000 ) THEN
348         ! Euler time stepping (first time step, starting from rest)
349         z2dt = rdt
350         DO jj = 1, jpj
351            DO ji = 1, jpi
352               zbsfa       = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
353               bsfb(ji,jj) = bsfn(ji,jj)
354               bsfn(ji,jj) = zbsfa 
355            END DO
356         END DO
357      ELSE
358         ! Leap-frog time stepping - Asselin filter
359         z2dt = 2.*rdt
360         DO jj = 1, jpj
361            DO ji = 1, jpi
362               zbsfa       = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
363               bsfb(ji,jj) = atfp * ( bsfb(ji,jj) + zbsfa ) + atfp1 * bsfn(ji,jj)
364               bsfn(ji,jj) = zbsfa
365            END DO
366         END DO
367      ENDIF
368
369# if defined key_obc
370      ! Swap of boundary arrays
371      IF( lpeastobc ) THEN
372      ! (jped,jpef),nieob
373      IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
374         DO jj = nje0m1, nje1
375            ! fields itm2 <== itm
376            bebnd(jj,ib  ,itm2) = bebnd(jj,ib  ,itm)
377            bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
378            bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
379            bebnd(jj,ib  ,itm ) = bebnd(jj,ib  ,it )
380         END DO
381      ELSE
382         ! fields itm <== it  plus time filter at the boundary
383         DO ji = fs_nie0, fs_nie1   ! vector opt.
384            DO jj = nje0m1, nje1
385               bebnd(jj,ib  ,itm2) = bebnd(jj,ib  ,itm)
386               bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
387               bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
388               bebnd(jj,ib  ,itm ) = atfp * ( bebnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bebnd(jj,ib,it)
389               bebnd(jj,ibm ,itm ) = bebnd(jj,ibm ,it )
390               bebnd(jj,ibm2,itm ) = bebnd(jj,ibm2,it )
391            END DO
392         END DO
393      ENDIF
394      ! fields it <== now (kt+1)
395      DO ji = fs_nie0, fs_nie1   ! vector opt.
396         DO jj = nje0m1, nje1
397            bebnd(jj,ib  ,it  ) = bsfn (ji  ,jj)
398            bebnd(jj,ibm ,it  ) = bsfn (ji-1,jj)
399            bebnd(jj,ibm2,it  ) = bsfn (ji-2,jj)
400         END DO
401      END DO
402#   if defined key_mpp
403      CALL mppobc( bebnd, jpjed, jpjef, jpieob, 3*3, 2, jpj )
404#   endif
405      ENDIF
406
407      IF( lpwestobc ) 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 defined key_mpp
439      CALL mppobc( bwbnd, jpjwd, jpjwf, jpiwob, 3*3, 2, jpj )
440#   endif
441      ENDIF
442
443      IF( lpnorthobc ) 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 defined key_mpp
475      CALL mppobc( bnbnd, jpind, jpinf, jpjnob, 3*3, 1, jpi )
476#   endif
477      ENDIF
478
479      IF( lpsouthobc ) THEN
480      ! njsob,(jpsd,jpsf)
481      IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
482         DO ji = nis0m1, nis1
483            ! fields itm2 <== itm
484            bsbnd(ji,ib  ,itm2) = bsbnd(ji,ib  ,itm)
485            bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
486            bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
487            bsbnd(ji,ib  ,itm ) = bsbnd(ji,ib  ,it )
488         END DO
489      ELSE
490         DO jj = fs_njs0, fs_njs1   ! vector opt.
491            DO ji = nis0m1, nis1
492               bsbnd(ji,ib  ,itm2) = bsbnd(ji,ib  ,itm)
493               bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
494               bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
495               ! fields itm <== it  plus time filter at the boundary
496               bsbnd(jj,ib  ,itm ) = atfp * ( bsbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bsbnd(jj,ib,it)
497               bsbnd(ji,ibm ,itm ) = bsbnd(ji,ibm ,it )
498               bsbnd(ji,ibm2,itm ) = bsbnd(ji,ibm2,it )
499            END DO
500         END DO
501      ENDIF
502      DO jj = fs_njs0, fs_njs1   ! vector opt.
503         DO ji = nis0m1, nis1 
504            ! fields it <== now (kt+1)
505            bsbnd(ji,ib  ,it  ) = bsfn (ji,jj  )
506            bsbnd(ji,ibm ,it  ) = bsfn (ji,jj+1)
507            bsbnd(ji,ibm2,it  ) = bsfn (ji,jj+2)
508         END DO
509      END DO
510#   if defined key_mpp
511      CALL mppobc( bsbnd, jpisd, jpisf, jpjsob, 3*3, 1, jpi )
512#   endif
513      ENDIF
514# endif
515      !
516      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
517     
518      !  add the surface pressure trend to the general trend
519      ! -----------------------------------------------------
520     
521      DO jj = 2, jpjm1
522
523         ! update the surface pressure gradient with the barotropic trend
524         DO ji = 2, jpim1
525            spgu(ji,jj) = spgu(ji,jj) + hur(ji,jj) / e2u(ji,jj) * ( bsfd(ji,jj) - bsfd(ji  ,jj-1) )
526            spgv(ji,jj) = spgv(ji,jj) - hvr(ji,jj) / e1v(ji,jj) * ( bsfd(ji,jj) - bsfd(ji-1,jj  ) )
527         END DO
528         ! add the surface pressure gradient trend to the general trend
529         DO jk = 1, jpkm1
530            DO ji = 2, jpim1
531               ua(ji,jj,jk) = ua(ji,jj,jk) - spgu(ji,jj)
532               va(ji,jj,jk) = va(ji,jj,jk) - spgv(ji,jj)
533# if defined key_trddyn
534               ! save the surface pressure gradient trend for diagnostics
535               utrd(ji,jj,jk,8) = -spgu(ji,jj)
536               vtrd(ji,jj,jk,8) = -spgv(ji,jj)
537# endif
538            END DO
539         END DO
540
541      END DO
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(*,*) 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.