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

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

nemo_v1_update_071:RB: add iom for restart and reorganization of restart

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