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 tags/nemo_v1_07/NEMO/OPA_SRC/DYN – NEMO

source: tags/nemo_v1_07/NEMO/OPA_SRC/DYN/dynspg_rl.F90 @ 8864

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

nemo_v1_update_017:RB: added call to sol_sor_e for free surface and rigid lid.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.2 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 solsor_e        ! Successive Over-relaxation solver with MPP optimization
28   USE solisl          ! ???
29   USE obc_oce         ! Lateral open boundary condition
30   USE lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   !! * Accessibility
37   PUBLIC dyn_spg_rl   ! called by step.F90
38
39   !! * Shared module variables
40   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl = .TRUE.    !: rigid-lid flag
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45#  include "obc_vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!   OPA 9.0 , LOCEAN-IPSL (2005)
48   !! $Header$
49   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
50   !!----------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE dyn_spg_rl( kt, kindic )
55      !!----------------------------------------------------------------------
56      !!                  ***  routine dyn_spg_rl  ***
57      !!
58      !! ** Purpose :   Compute the now trend due to the surface pressure
59      !!      gradient for the rigid-lid case,  add it to the general trend of
60      !!      momentum equation.
61      !!
62      !! ** Method  :   Rigid-lid appromimation: the surface pressure gradient
63      !!      is given by:
64      !!         spgu = 1/rau0 d/dx(ps) = Mu + 1/(hu e2u) dj-1(bsfd)
65      !!         spgv = 1/rau0 d/dy(ps) = Mv - 1/(hv e1v) di-1(bsfd)
66      !!      where (Mu,Mv) is the vertically averaged momentum trend (i.e.
67      !!      the vertical ponderated sum of the general momentum trend),
68      !!      and bsfd is the barotropic streamfunction trend.
69      !!      The trend is computed as follows:
70      !!         -1- compute the vertically averaged momentum trend (Mu,Mv)
71      !!         -2- compute the barotropic streamfunction trend by solving an
72      !!      ellipic equation using a diagonal preconditioned conjugate
73      !!      gradient or a successive-over-relaxation method (depending
74      !!      on nsolv, a namelist parameter).
75      !!         -3- add to bsfd the island trends if lk_isl=T
76      !!         -4- compute the after streamfunction is for further diagnos-
77      !!      tics using a leap-frog scheme.
78      !!         -5- add the momentum trend associated with the surface pres-
79      !!      sure gradient to the general trend.
80      !!
81      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
82      !!             - Save the trends in (ztdua,ztdva) ('key_trddyn')
83      !!
84      !! References :
85      !!      Madec et al. 1988, ocean modelling, issue 78, 1-6.
86      !!
87      !! History :
88      !!        !  96-05  (G. Madec, M. Imbard, M. Guyon)  rewitting in 1
89      !!                  routine, without pointers, and with the same matrix
90      !!                  for sor and pcg, mpp exchange, and symmetric conditions
91      !!        !  96-07  (A. Weaver)  Euler forward step
92      !!        !  96-11  (A. Weaver)  correction to preconditioning:
93      !!        !  98-02  (M. Guyon)  FETI method
94      !!        !  98-05  (G. Roullet)  free surface
95      !!        !  97-09  (J.-M. Molines)  Open boundaries
96      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
97      !!        !  02-11  (C. Talandier, A-M Treguier) Open boundaries
98      !!   9.0  !  04-08  (C. Talandier)  New trends organization
99      !!---------------------------------------------------------------------
100      !! * Modules used     
101      USE oce, ONLY :    ztdua => ta,    & ! use ta as 3D workspace   
102                         ztdva => sa       ! use sa as 3D workspace   
103      !! * Arguments
104      INTEGER, INTENT( in  ) ::   kt       ! ocean time-step index
105      INTEGER, INTENT( out ) ::   kindic   ! solver flag, take a negative value
106      !                                    ! when the solver doesnot converge
107      !! * Local declarations
108      INTEGER ::   ji, jj, jk    ! dummy loop indices
109      REAL(wp) ::  zbsfa, zgcx, z2dt
110# if defined key_obc
111      INTEGER ::   ip, ii, ij
112      INTEGER ::   iii, ijj, jip, jnic
113      INTEGER ::   it, itm, itm2, ib, ibm, ibm2
114      REAL(wp) ::   z2dtr
115# endif
116      !!----------------------------------------------------------------------
117     
118      IF( kt == nit000 ) THEN
119         IF(lwp) WRITE(numout,*)
120         IF(lwp) WRITE(numout,*) 'dyn_spg_rl : surface pressure gradient trend'
121         IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
122
123         ! set to zero rigid-lid specific arrays
124         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction)
125         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
126      ENDIF
127
128      ! 0. Initializations:
129      ! -------------------
130# if defined key_obc
131      ! space index on boundary arrays
132      ib = 1
133      ibm = 2
134      ibm2 = 3
135      ! time index on boundary arrays
136      it = 1
137      itm = 2
138      itm2 = 3
139# endif
140
141      !                                                ! ===============
142      DO jj = 2, jpjm1                                 !  Vertical slab
143         !                                             ! ===============
144
145         ! 1. Vertically averaged momentum trend
146         ! -------------------------------------
147         ! initialization to zero
148         spgu(:,jj) = 0.
149         spgv(:,jj) = 0.
150         ! vertical sum
151         DO jk = 1, jpk
152            DO ji = 2, jpim1
153               spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
154               spgv(ji,jj) = spgv(ji,jj) + va(ji,jj,jk) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)
155            END DO
156         END DO 
157         ! divide by the depth
158         spgu(:,jj) = spgu(:,jj) * hur(:,jj)
159         spgv(:,jj) = spgv(:,jj) * hvr(:,jj)
160
161         !                                             ! ===============
162      END DO                                           !   End of slab
163      !                                                ! ===============
164
165      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
166
167      ! Boundary conditions on (spgu,spgv)
168      CALL lbc_lnk( spgu, 'U', -1. )
169      CALL lbc_lnk( spgv, 'V', -1. )
170
171      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
172
173      ! 2. Barotropic streamfunction trend (bsfd)
174      ! ----------------------------------
175
176      ! Curl of the vertically averaged velocity
177      DO jj = 2, jpjm1
178         DO ji = 2, jpim1
179            gcb(ji,jj) = -gcdprc(ji,jj)   &
180                       *(  ( e2v(ji+1,jj  )*spgv(ji+1,jj  ) - e2v(ji,jj)*spgv(ji,jj) )   &
181                          -( e1u(ji  ,jj+1)*spgu(ji  ,jj+1) - e1u(ji,jj)*spgu(ji,jj) )   ) 
182         END DO
183      END DO
184
185# if defined key_obc
186      ! Open boundary contribution
187      DO jj = 2, jpjm1
188         DO ji = 2, jpim1
189           gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj)
190         END DO
191      END DO
192# else
193      ! No open boundary contribution
194# endif
195
196      ! First guess using previous solution of the elliptic system and
197      ! not bsfd since the system is solved with 0 as coastal boundary
198      ! condition. Also include a swap array (gcx,gxcb)
199      DO jj = 2, jpjm1
200         DO ji = 2, jpim1
201            zgcx        = gcx(ji,jj)
202            gcx (ji,jj) = 2.*zgcx - gcxb(ji,jj)
203            gcxb(ji,jj) =    zgcx
204         END DO
205      END DO
206      ! applied the lateral boundary conditions
207      IF( nsolv == 4)   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
208
209      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
210
211      ! Relative precision (computation on one processor)
212      rnorme = 0.e0
213      rnorme = SUM( gcb(1:nlci,1:nlcj) * gcdmat(1:nlci,1:nlcj) * gcb(1:nlci,1:nlcj) * bmask(:,:) )
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( 4 )                     ! successive-over-relaxation with extra outer halo
237            CALL sol_sor_e( 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 ,3 or 4'
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(1:nlci,1:nlcj) = gcx(1:nlci,1:nlcj)
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.