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.
diaspr.F90 in trunk/NEMO/OPA_SRC/DIA – NEMO

source: trunk/NEMO/OPA_SRC/DIA/diaspr.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.4 KB
Line 
1MODULE diaspr
2   !!======================================================================
3   !!                       ***  MODULE  diaspr  ***
4   !! Ocean diagnostics:  surface pressure (rigid-lid case)
5   !!=====================================================================
6#if   defined key_diaspr   &&   defined key_dynspg_rl
7   !!----------------------------------------------------------------------
8   !!   'key_diaspr'        and                surface pressure diagnostics
9   !!   'key_dynspg_rl'                                      rigid-lid case
10   !!----------------------------------------------------------------------
11   !!   dia_spr      : update momentum and tracer Kz from a tke scheme
12   !!   sprmat       : initialization, namelist read, and parameters control
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! physical constants
18   USE in_out_manager  ! I/O manager
19   USE sol_oce         ! ocean elliptic solver
20   USE solpcg          ! preconditionned conjugate gradient solver
21   USE solsor          ! Successive Over-relaxation solver
22   USE solfet          ! FETI solver
23   USE lib_mpp         ! distributed memory computing library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Routine accessibility
29   PUBLIC dia_spr   ! routine called by step.F90
30
31   !! * Shared module variables
32   LOGICAL, PUBLIC, PARAMETER ::   lk_diaspr = .TRUE.    !: surface pressure diag. flag
33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gps         !: surface pressure
34
35   !! * Module variables
36   INTEGER ::                 &
37      nmoyps,                 &  ! time step for average
38      nindic,                 &  ! indicator of convergence of the solver
39      !                          ! namspr  surface pressure diagnostic
40      nmaxp ,                 &  ! maximum of iterations for the solver
41      niterp                     ! number of iteration done by the solver
42
43   REAL(wp) ::     &
44      ! namspr  surface pressure diagnostic
45      epsp                       ! absolute precision of the solver
46
47      !! * Namelist
48      NAMELIST/namspr/ nmaxp, epsp, niterp
49
50   REAL(wp) ::     &
51      e1e2t                      ! ???
52
53   REAL(wp), PUBLIC DIMENSION(jpi,jpj) ::   &
54      spgum, spgvm,           &  ! average value of the surface pressure gradients
55      gpsuu, gpsvv,           &  ! surface pressure gradients computed from comp. PS
56      gcdpsc,                 &  ! inverse diagonal preconditioning matrix
57      gcsmat,                 &  ! diagonal preconditioning matrix
58      spmsk                      ! surface pressure Mask
59
60   REAL(wp), DIMENSION(jpi,jpj,4) ::   &
61      gcps                       ! extra-diagonal elements of SPG matrix
62   !!----------------------------------------------------------------------
63   !!   OPA 9.0 , LOCEAN-IPSL (2005)
64   !! $Header$
65   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
66   !!----------------------------------------------------------------------
67
68CONTAINS
69
70   SUBROUTINE dia_spr( kt )
71      !!---------------------------------------------------------------------
72      !!                  ***  ROUTINE dia_spr  ***
73      !!
74      !! ** Purpose :   compute the surface pressure from its gradient
75      !!
76      !! ** Method  :   rigid-lid appromimation: the surface pressure
77      !!      gradient is given by:
78      !!           spgu = 1/rau0 d/dx(ps) = Mu + 1/(hu e2u) dj-1(bsfd)
79      !!           spgv = 1/rau0 d/dy(ps) = Mv - 1/(hv e1v) di-1(bsfd)
80      !!
81      !!      where (Mu,Mv) is the vertically averaged momentum trend, i.e.
82      !!      the vertical ponderated sum of the general momentum trend.
83      !!      where bsfd is the trend of the barotropic stream function.
84      !!
85      !!       taking the divergence of the surface pressure gradient provides
86      !!      an elliptic equation for ps which is solved using either a
87      !!      diagonal preconditioned conjugate gradient method (solpcg.f) or
88      !!      an successive-over-relaxation method (solsor.f) or FETI method
89      !!      (solfet.F).
90      !!
91      !!      n.b. this resolution is valid with topography, cyclic east-west
92      !!      boundary conditions and islands.
93      !!
94      !! History :
95      !!        !  98-01  (G. Madec & M. Ioualalen)  Original code
96      !!        !  98-02  (M. Guyon)  FETI method
97      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
98      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
99      !!----------------------------------------------------------------------
100      !! * Arguments
101      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
102
103      !! * Local declarations
104      INTEGER  ::  ji, jj
105      INTEGER  ::   imax, ijt, iju
106      REAL(wp) ::   zpsmea, zeps, zmoyr
107      REAL(wp) ::   ztab(jpi,jpj,8)
108      REAL(wp) ::   zemin1, zemax1, zemin2, zemax2, zgwgt
109      REAL(wp) ::   z1, z2, zcompt,z3,z4
110      REAL(wp) ::   zdif1, zdif2, zvar1, zvar2
111      !!----------------------------------------------------------------------
112
113
114      ! 0. initialisation (the first time step)
115      ! ---------------------------------------
116     
117      IF( kt == nit000 ) THEN
118
119         ! Namelist namspr : surface pressure
120
121         nmaxp  = 2000
122         epsp   = 1.e-6
123         niterp = 16
124
125         ! Read Namelist namspr : surface pressure diagnostics
126         REWIND ( numnam )
127         READ(numnam,namspr)
128
129         IF(lwp) THEN
130            WRITE(numout,*) 'dia_spr : surface pressure diagnostic (rigid-lid case)'
131            WRITE(numout,*) '~~~~~~~'
132            WRITE(numout,*)
133            WRITE(numout,*) '          Namelist namspr : set solver parameters'
134            WRITE(numout,*)
135            WRITE(numout,*) '             maximum iterations for solver  nmaxp  = ', nmaxp
136            WRITE(numout,*) '             absolute precision of solver   epsp   = ', epsp
137            WRITE(numout,*) '             number of solver iterations    niterp = ', niterp
138            WRITE(numout,*) '             frequeny of averaged output    nwrite = ', nwrite
139            WRITE(numout,*)
140         ENDIF
141
142         ! control
143# if ! defined key_dynspg_rl
144      CALL ctl_stop( '          surface pressure already explicitly computed !!' )
145# endif
146
147         ! compute the ocean surface
148         e1e2t = 0.e0
149         DO jj = 2, jpjm1
150            DO ji = 2, jpim1
151               e1e2t = e1e2t + e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,1)
152            END DO
153         END DO
154         IF( lk_mpp )   CALL  mpp_sum( e1e2t )   ! sum over the global domain
155         
156         ! build the matrix for the surface pressure
157         CALL sprmat
158         
159         ! set to zero the mean surface pressure gradient
160         nmoyps = 0
161         spgum(:,:) = 0.e0
162         spgvm(:,:) = 0.e0
163
164      ENDIF
165
166      ! 1. cumulate the surface pressure gradient (at each time step)
167      ! -----------------------------------------
168
169      nmoyps = nmoyps + 1
170      spgum(:,:) = spgum(:,:) + spgu(:,:)
171      spgvm(:,:) = spgvm(:,:) + spgv(:,:)
172     
173
174      ! 2. ps computation each nwrite time step
175      ! ---------------------------------------
176     
177      ! RETURN IF not the right time to compute ps
178      IF ( MOD(kt-nit000+1,nwrite) /= 0 ) RETURN
179     
180     
181      ! mean surface pressure gradient
182      !   averaging and mask
183      zmoyr = 1./float(nmoyps)
184      DO jj = 2, jpjm1
185         DO ji = 2, jpim1
186            spgum(ji,jj) = spgum(ji,jj) * zmoyr * umask(ji,jj,1)
187            spgvm(ji,jj) = spgvm(ji,jj) * zmoyr * vmask(ji,jj,1)
188         END DO
189      END DO
190
191      CALL  lbc_lnk(spgum, 'U', -1. )
192      CALL  lbc_lnk(spgvm, 'V', -1. )
193
194     
195      ! SAVE in local arrays and variables of solver informations
196      zeps   = eps
197      imax  = nmax 
198      ztab(:,:,1) = gcp   (:,:,1)
199      ztab(:,:,2) = gcp   (:,:,2)
200      ztab(:,:,3) = gcp   (:,:,3)
201      ztab(:,:,4) = gcp   (:,:,4)
202      ztab(:,:,5) = gcdprc(:,:  )
203      ztab(:,:,6) = gcdmat(:,:  )
204      ztab(:,:,7) = gcx   (:,:  )
205      ztab(:,:,8) = bmask (:,:  )
206
207      ! replace bsf solver informations by ps solver one
208      eps    = epsp
209      nmax   = nmaxp
210      gcp   (:,:,1) = gcps  (:,:,1)
211      gcp   (:,:,2) = gcps  (:,:,2)
212      gcp   (:,:,3) = gcps  (:,:,3)
213      gcp   (:,:,4) = gcps  (:,:,4)
214      gcdprc(:,:  ) = gcdpsc(:,:  )
215      gcdmat(:,:  ) = gcsmat(:,:  )
216      bmask (:,:  ) = spmsk (:,:  )
217      !    first guess: ps
218      gcx   (:,:  ) = gps   (:,:  )
219
220      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
221
222      ! right hand side: 2d div. of the surface pressure gradient
223      DO jj = 2, jpjm1
224         DO ji = 2, jpim1
225            gcb(ji,jj) = -gcdpsc(ji,jj)*   &
226               (  e2u(ji,jj)*spgum(ji,jj) - e2u(ji-1,jj)*spgum(ji-1,jj)   &
227               + e1v(ji,jj)*spgvm(ji,jj) - e1v(ji,jj-1)*spgvm(ji,jj-1) )
228         END DO
229      END DO
230     
231      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
232     
233      ! relative PRECISION
234      rnorme = 0.
235      DO jj = 1, jpj
236         DO ji = 1, jpi
237            rnorme = rnorme + gcb(ji,jj) * gcsmat(ji,jj) * gcb(ji,jj)
238         END DO
239      END DO
240      IF( lk_mpp )   CALL  mpp_sum( rnorme )   ! sum over the global domain
241
242      epsr=eps*eps*rnorme
243      ncut=0
244      !   IF the second member is 0 the solution is 0, solpcg isn't called
245      IF ( rnorme == 0.e0 ) THEN
246         gps(:,:) = 0.e0
247         res   = 0.e0
248         niter = 0
249         ncut  = 999
250      ENDIF
251     
252      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
253     
254      nindic = 0
255
256      ! iterarive solver of the spg system (except IF sol.=0)
257      !     (OUTPUT in gcx with boundary conditions applied)
258      IF ( ncut == 0 ) THEN
259         IF ( nsolv == 1 ) THEN
260            CALL sol_pcg( nindic )         !   diagonal preconditioned conjuguate gradient
261         ELSE IF ( nsolv == 2 ) THEN
262            CALL sol_sor( nindic )     !   successive-over-relaxation
263         ELSE IF(nsolv == 3) THEN
264            CALL sol_fet( nindic )         !   FETI solver
265         ELSE
266            !   e r r o r  in nsolv namelist PARAMETER
267            IF(lwp) THEN
268               WRITE(numout,*) ' dia_spr: e r r o r, nsolv = 1 or 2'
269               WRITE(numout,*) ' ******               not = ',nsolv
270            ENDIF
271            STOP 'dia_spr'
272         ENDIF
273      ENDIF
274     
275      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
276
277     
278      ! sp solver statistics  (i.e. problem for the solver)
279      IF ( epsr < 0.) THEN
280         IF(lwp) THEN
281            WRITE(numout,*)'rrrrrrrrrrrrrrr'
282            IF(lwp)WRITE(numout,*)'dia_spr-1:',epsr
283            IF(lwp)WRITE(numout,*)'rrrrrrrrrrrrrrr'
284         ENDIF
285      ENDIF
286      IF(lwp)WRITE(numout,9300) kt, niter, res, SQRT(epsr)/eps
287      IF (nindic < 0) THEN
288         IF(lwp) THEN
289            WRITE(numout,9100)
290            WRITE(numout,*) ' dia_spr : the surface pressure solver DO not converge'
291            WRITE(numout,*) ' ====== ' 
292            WRITE(numout,*) 
293         ENDIF
294      ENDIF
2959100  FORMAT( /,' ===>>>> : w a r n i n g',/,'          ===============',/ )
2969300  FORMAT(' it :', i8, ' niter :', i4, ' res :',e20.10,' b :', e20.10)
297
298      ! recover bsf solver informations and SAVE ps for next computation
299      eps    = zeps
300      nmax   = imax 
301      gps   (:,:  ) = gcx (:,:)
302      gcp   (:,:,1) = ztab(:,:,1)
303      gcp   (:,:,2) = ztab(:,:,2)
304      gcp   (:,:,3) = ztab(:,:,3)
305      gcp   (:,:,4) = ztab(:,:,4)
306      gcdprc(:,:  ) = ztab(:,:,5)
307      gcdmat(:,:  ) = ztab(:,:,6)
308      gcx   (:,:  ) = ztab(:,:,7)
309      bmask (:,:  ) = ztab(:,:,8)
310     
311      ! compute and substract the mean value
312     
313      zpsmea = 0.e0
314      DO jj=2,jpjm1
315         DO ji=2,jpim1
316            zpsmea = zpsmea + gps(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1)
317         END DO
318      END DO
319      IF( lk_mpp )   CALL  mpp_sum( zpsmea )   ! sum over the global domain
320
321      zpsmea = zpsmea / e1e2t
322      gps(:,:) = ( gps(:,:) - zpsmea ) * tmask(:,:,1)
323 
324      IF(lwp)WRITE(numout,*) ' mean value of ps = ',zpsmea,' is substracted'
325      ! ----------------------------------------
326      ! i. compute the surface pressure gradient
327      !    from the computed surface pressure
328      ! ----------------------------------------
329
330      DO jj=2,jpjm1
331         DO ji=2,jpim1
332            gpsuu(ji,jj)=(gps(ji+1,jj)-gps(ji,jj))/e1u(ji,jj) * umask(ji,jj,1)
333            gpsvv(ji,jj)=(gps(ji,jj+1)-gps(ji,jj))/e2v(ji,jj) * vmask(ji,jj,1)
334         END DO
335      END DO
336     
337      ! compute the max and min error
338     
339      zemax1 = 0.e0
340      zemin1 = 0.e0
341      zemax2 = 0.e0
342      zemin2 = 0.e0
343      DO jj = 2,jpj-1
344         DO ji = 2,jpi-1
345            z1 = ABS( spgum(ji,jj)-gpsuu(ji,jj) )*umask(ji,jj,1)
346            z2 = ABS( spgvm(ji,jj)-gpsvv(ji,jj) )*vmask(ji,jj,1)
347            z3 = MAX ( ABS( spgum(ji,jj) ), ABS( spgvm(ji,jj) ) )
348            z4 = MAX ( ABS( gpsuu(ji,jj) ), ABS( gpsvv(ji,jj) ) )
349            zemax1 = MAX(z1,zemax1)
350            zemax2 = MAX(z2,zemax2)
351            zemin1 = MAX(z3,zemin1)
352            zemin2 = MAX(z4,zemin2)
353         END DO
354      END DO
355      IF( lk_mpp )   CALL  mpp_sum( zemax1 )   ! sum over the global domain
356      IF( lk_mpp )   CALL  mpp_sum( zemax2 )   ! sum over the global domain
357      IF( lk_mpp )   CALL  mpp_sum( zemin1 )   ! sum over the global domain
358      IF( lk_mpp )   CALL  mpp_sum( zemin2 )   ! sum over the global domain
359
360      IF(lwp) THEN
361         WRITE(numout,*)
362         WRITE(numout,*) 'pserro : time step = ', kt
363         WRITE(numout,*) '******** ------------------'
364         WRITE(numout,*)
365         WRITE(numout,*) '         gpsx error  max=',zemax1
366         WRITE(numout,*) '         gpsy error  max=',zemax2
367         WRITE(numout,*) '         gps max =',zemin1
368         WRITE(numout,*) '         gpsc max =',zemin2
369         WRITE(numout,*)
370      ENDIF
371
372      ! compute the norme and variance of this error
373
374      zcompt = 0.e0
375      zdif1  = 0.e0
376      zdif2  = 0.e0
377      zvar1  = 0.e0
378      zvar2  = 0.e0
379      DO jj = 2, jpj-1
380         DO ji = 2, jpi-1
381            z1 = ( spgum(ji,jj)-gpsuu(ji,jj) ) * umask(ji,jj,1)
382            z2 = ( spgvm(ji,jj)-gpsvv(ji,jj) ) * vmask(ji,jj,1)
383            zcompt=zcompt+tmask(ji,jj,1)
384            zdif1=zdif1+z1
385            zdif2=zdif2+z2
386            zvar1=zvar1+z1*z1
387            zvar2=zvar2+z2*z2
388         END DO
389      END DO
390      IF( lk_mpp )   CALL  mpp_sum( zcompt )   ! sum over the global domain
391      IF( lk_mpp )   CALL  mpp_sum( zdif1  )   ! sum over the global domain
392      IF( lk_mpp )   CALL  mpp_sum( zdif2  )   ! sum over the global domain
393      IF( lk_mpp )   CALL  mpp_sum( zvar1  )   ! sum over the global domain
394      IF( lk_mpp )   CALL  mpp_sum( zvar2  )   ! sum over the global domain
395
396      IF(lwp) WRITE(numout,*) '        zcompt = ',zcompt
397      zdif1=zdif1/zcompt
398      zdif2=zdif2/zcompt
399      IF( zvar1 < 0.) THEN
400         IF(lwp) THEN
401            WRITE(numout,*)'rrrrrrrrrrrrrrr'
402            WRITE(numout,*)'dia_spr-2:',zvar1
403            WRITE(numout,*)'rrrrrrrrrrrrrrr'
404         ENDIF
405      ENDIF
406      zvar1 = SQRT(zvar1)/zcompt
407      IF( zvar2 < 0. ) THEN
408         IF(lwp)THEN
409            WRITE(numout,*)'rrrrrrrrrrrrrrr'
410            WRITE(numout,*)'dia_spr-3:',zvar2
411            WRITE(numout,*)'rrrrrrrrrrrrrrr'
412         ENDIF
413      ENDIF
414      zvar2 = SQRT(zvar2)/zcompt
415     
416      IF(lwp) THEN
417         WRITE(numout,*)
418         WRITE(numout,*) '         gpsx mean error = ',zdif1
419         WRITE(numout,*) '         gpsy mean error = ',zdif2
420         WRITE(numout,*)
421         WRITE(numout,*) '         gpsx var. error = ',zvar1
422         WRITE(numout,*) '         gpsy var. error = ',zvar2
423         WRITE(numout,*)
424         WRITE(numout,*)
425      ENDIF
426     
427      ! reset to zero nmoyps and the mean surface pressure gradient
428      nmoyps = 0
429      spgum(:,:) = 0.e0
430      spgvm(:,:) = 0.e0
431     
432   END SUBROUTINE dia_spr
433
434
435   SUBROUTINE sprmat
436      !!---------------------------------------------------------------------
437      !!                  ***  ROUTINE sprmat  ***
438      !!               
439      !! ** Purpose :   construction of the matrix of the surface pressure
440      !!      system and the diagonal preconditioning matrix.
441      !!
442      !! ** Method :
443      !!
444      !! History :
445      !!        !  98-01  (G. Madec & M. Ioualalen)  Original code
446      !!   8.5  !  02-08  (G. Madec)  F90: Free form
447      !!----------------------------------------------------------------------
448      !! * Local declarations
449      INTEGER ::   ji, jj, jl                ! dummy loop indices
450      REAL(wp) ::   zcoefs, zcoefw, zcoefe, zcoefn
451      !!----------------------------------------------------------------------
452     
453     
454      ! 0. ocean/land mask at ps-point (computed from surface tmask): spmsk
455      ! --------------------------------------------------------------------
456     
457      ! computation
458      spmsk(:,:) = tmask(:,:,1)
459     
460      ! boundary conditions
461      ! south symmetry: psmsk must be set to 0. on 1
462      IF( nperio == 2 ) THEN
463         spmsk(:, 1 ) = 0.e0
464      ENDIF
465     
466      ! east-west cyclic: spmsk must be set to 0. on 1 and jpi
467      IF( nperio == 1 .OR. nperio == 4 .OR.nperio == 6) THEN
468         spmsk( 1 ,:) = 0.e0
469         spmsk(jpi,:) = 0.e0
470      ENDIF
471     
472      ! north fold: spmsk must be set to 0. on ligne jpj and on half
473      !                   ligne jpj-1
474      ! T-point pivot
475      IF( nperio == 3 .OR. nperio == 4 ) THEN
476         spmsk(:,jpj) = 0.e0
477         DO ji = jpi/2+1, jpi
478            spmsk(ji,jpjm1) = 0.e0
479         END DO
480      ENDIF
481      ! F-point pivot
482      IF( nperio == 5 .OR. nperio == 6 ) THEN
483         spmsk(:,jpj) = 0.e0
484      ENDIF
485     
486      ! mpp boundary cond.: spmsk is initialized at zero on the overlap
487      ! region for both the preconjugate gradient and the sor algorithms
488
489      IF( nbondi /= -1 .AND. nbondi /= 2 ) THEN
490         DO jl = 1, jpreci
491            spmsk(jl,:) = 0.e0
492         END DO
493      ENDIF
494      IF( nbondi /= 1 .AND. nbondi /= 2 ) THEN
495         DO ji = nlci, jpi
496            spmsk(ji,:) = 0.e0
497         END DO
498      ENDIF
499      IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
500         DO jl=1,jprecj
501            spmsk(:,jl) = 0.e0
502         END DO
503      ENDIF
504      IF( nbondj /= 1 .AND. nbondj /= 2 ) THEN
505         DO jj = nlcj, jpj
506            spmsk(:,jj) = 0.e0
507         END DO
508      ENDIF
509     
510      ! 1. construction of the matrix
511      ! -----------------------------
512     
513      DO jj = 1, jpj
514         DO ji = 1, jpi
515           
516            IF( spmsk(ji,jj) == 0. ) THEN
517               ! land points
518               gcps  (ji,jj,1) = 0.e0
519               gcps  (ji,jj,2) = 0.e0
520               gcps  (ji,jj,3) = 0.e0
521               gcps  (ji,jj,4) = 0.e0
522               gcdpsc(ji,jj  ) = 0.e0
523               gcsmat(ji,jj  ) = 0.e0
524            ELSE
525               ! south coefficient
526               zcoefs = -e1v(ji,jj-1) / e2v(ji,jj-1) * vmask(ji,jj-1,1)
527               gcps(ji,jj,1) = zcoefs
528               ! west coefficient
529               zcoefw = -e2u(ji-1,jj) / e1u(ji-1,jj) * umask(ji-1,jj,1)
530               gcps(ji,jj,2) = zcoefw
531               ! east coefficient
532               zcoefe = -e2u(ji  ,jj) / e1u(ji  ,jj) * umask(ji  ,jj,1)
533               gcps(ji,jj,3) = zcoefe
534               ! north coefficient
535               zcoefn = -e1v(ji,jj  ) / e2v(ji,jj  ) * vmask(ji,jj  ,1)
536               gcps(ji,jj,4) = zcoefn
537               
538               ! diagonal coefficient
539               gcsmat(ji,jj) = -zcoefs-zcoefw-zcoefe-zcoefn
540            ENDIF
541         END DO
542      END DO
543     
544
545      ! 2. boundary conditions
546      ! ----------------------
547     
548      ! cyclic east-west boundary conditions
549      ! ji=2 is the column east of ji=jpim1 and reciprocally,
550      ! ji=jpim1 is the column west of ji=2
551      ! all the coef are already set to zero as spmask is initialized to
552      ! zero for ji=1 and ji=jpj.
553     
554      ! symetrical conditions
555      ! the diagonal coefficient of the southern grid points must be modify to
556      ! account for the existence of the south symmetric bassin.
557      IF( nperio == 2 ) THEN
558         DO ji = 1, jpi
559            IF( spmsk(ji,2) /= 0 ) THEN
560               zcoefs = e1v(ji,1) / e2v(ji,1)
561               gcsmat(ji,2) = gcsmat(ji,2) - zcoefs
562            ENDIF
563         END DO
564      ENDIF
565     
566      ! North fold boundary condition
567      ! all the coef are already set to zero as bmask is initialized to
568      ! zero on duplicated lignes and portion of lignes
569     
570     
571      ! 3. preconditioned matrix
572      ! ------------------------
573     
574      DO jj = 1, jpj
575         DO ji = 1, jpi
576            IF( spmsk(ji,jj) /= 0. ) gcdpsc(ji,jj) = 1.e0 / gcsmat(ji,jj)
577         END DO
578      END DO
579     
580      gcps(:,:,1) = gcps(:,:,1) * gcdpsc(:,:)
581      gcps(:,:,2) = gcps(:,:,2) * gcdpsc(:,:)
582      gcps(:,:,3) = gcps(:,:,3) * gcdpsc(:,:)
583      gcps(:,:,4) = gcps(:,:,4) * gcdpsc(:,:)
584     
585     
586      ! 3. initialization the arrays used in sp solver
587      ! ----------------------------------------------
588     
589      gps  (:,:) = 0.e0
590      gpsuu(:,:) = 0.e0
591      gpsvv(:,:) = 0.e0
592     
593   END SUBROUTINE sprmat
594
595#else
596   !!----------------------------------------------------------------------
597   !!   Default option :                    NO surface pressure diagnostics
598   !!----------------------------------------------------------------------
599   LOGICAL, PUBLIC, PARAMETER ::   lk_diaspr = .FALSE.   !: surface pressure diag. flag
600CONTAINS
601   SUBROUTINE dia_spr( kt )      ! Empty routine
602      WRITE(*,*) 'dia_spr: You should not have seen this print! error?', kt
603   END SUBROUTINE dia_spr
604#endif
605
606   !!======================================================================
607END MODULE diaspr
Note: See TracBrowser for help on using the repository browser.