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

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

CL : Add CVS Header and CeCILL licence information
  • 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      !!----------------------------------------------------------------------
99      !! * Arguments
100      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
101
102      !! * Local declarations
103      INTEGER  ::  ji, jj
104      INTEGER  ::   imax, ijt, iju
105      REAL(wp) ::   zpsmea, zeps, zmoyr
106      REAL(wp) ::   ztab(jpi,jpj,8)
107      REAL(wp) ::   zemin1, zemax1, zemin2, zemax2, zgwgt
108      REAL(wp) ::   z1, z2, zcompt,z3,z4
109      REAL(wp) ::   zdif1, zdif2, zvar1, zvar2
110      !!----------------------------------------------------------------------
111
112
113      ! 0. initialisation (the first time step)
114      ! ---------------------------------------
115     
116      IF( kt == nit000 ) THEN
117
118         ! Namelist namspr : surface pressure
119
120         nmaxp  = 2000
121         epsp   = 1.e-6
122         niterp = 16
123
124         ! Read Namelist namspr : surface pressure diagnostics
125         REWIND ( numnam )
126         READ(numnam,namspr)
127
128         IF(lwp) THEN
129            WRITE(numout,*) 'dia_spr : surface pressure diagnostic (rigid-lid case)'
130            WRITE(numout,*) '~~~~~~~'
131            WRITE(numout,*)
132            WRITE(numout,*) '          Namelist namspr : set solver parameters'
133            WRITE(numout,*)
134            WRITE(numout,*) '             maximum iterations for solver  nmaxp  = ', nmaxp
135            WRITE(numout,*) '             absolute precision of solver   epsp   = ', epsp
136            WRITE(numout,*) '             number of solver iterations    niterp = ', niterp
137            WRITE(numout,*) '             frequeny of averaged output    nwrite = ', nwrite
138            WRITE(numout,*)
139         ENDIF
140
141         ! control
142# if defined key_dynspg_fsc
143      IF(lwp) WRITE(numout,cform_err)
144      IF(lwp) WRITE(numout,*) '          surface pressure already explicitly computed !!'
145      nstop = nstop + 1
146# endif
147
148         ! compute the ocean surface
149         e1e2t = 0.e0
150         DO jj = 2, jpjm1
151            DO ji = 2, jpim1
152               e1e2t = e1e2t + e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,1)
153            END DO
154         END DO
155         IF( lk_mpp )   CALL  mpp_sum( e1e2t )   ! sum over the global domain
156         
157         ! build the matrix for the surface pressure
158         CALL sprmat
159         
160         ! set to zero the mean surface pressure gradient
161         nmoyps = 0
162         spgum(:,:) = 0.e0
163         spgvm(:,:) = 0.e0
164
165      ENDIF
166
167      ! 1. cumulate the surface pressure gradient (at each time step)
168      ! -----------------------------------------
169
170      nmoyps = nmoyps + 1
171      spgum(:,:) = spgum(:,:) + spgu(:,:)
172      spgvm(:,:) = spgvm(:,:) + spgv(:,:)
173     
174
175      ! 2. ps computation each nwrite time step
176      ! ---------------------------------------
177     
178      ! RETURN IF not the right time to compute ps
179      IF ( MOD(kt-nit000+1,nwrite) /= 0 ) RETURN
180     
181     
182      ! mean surface pressure gradient
183      !   averaging and mask
184      zmoyr = 1./float(nmoyps)
185      DO jj = 2, jpjm1
186         DO ji = 2, jpim1
187            spgum(ji,jj) = spgum(ji,jj) * zmoyr * umask(ji,jj,1)
188            spgvm(ji,jj) = spgvm(ji,jj) * zmoyr * vmask(ji,jj,1)
189         END DO
190      END DO
191
192      CALL  lbc_lnk(spgum, 'U', -1. )
193      CALL  lbc_lnk(spgvm, 'V', -1. )
194
195     
196      ! SAVE in local arrays and variables of solver informations
197      zeps   = eps
198      imax  = nmax 
199      ztab(:,:,1) = gcp   (:,:,1)
200      ztab(:,:,2) = gcp   (:,:,2)
201      ztab(:,:,3) = gcp   (:,:,3)
202      ztab(:,:,4) = gcp   (:,:,4)
203      ztab(:,:,5) = gcdprc(:,:  )
204      ztab(:,:,6) = gcdmat(:,:  )
205      ztab(:,:,7) = gcx   (:,:  )
206      ztab(:,:,8) = bmask (:,:  )
207
208      ! replace bsf solver informations by ps solver one
209      eps    = epsp
210      nmax   = nmaxp
211      gcp   (:,:,1) = gcps  (:,:,1)
212      gcp   (:,:,2) = gcps  (:,:,2)
213      gcp   (:,:,3) = gcps  (:,:,3)
214      gcp   (:,:,4) = gcps  (:,:,4)
215      gcdprc(:,:  ) = gcdpsc(:,:  )
216      gcdmat(:,:  ) = gcsmat(:,:  )
217      bmask (:,:  ) = spmsk (:,:  )
218      !    first guess: ps
219      gcx   (:,:  ) = gps   (:,:  )
220
221      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
222
223      ! right hand side: 2d div. of the surface pressure gradient
224      DO jj = 2, jpjm1
225         DO ji = 2, jpim1
226            gcb(ji,jj) = -gcdpsc(ji,jj)*   &
227               (  e2u(ji,jj)*spgum(ji,jj) - e2u(ji-1,jj)*spgum(ji-1,jj)   &
228               + e1v(ji,jj)*spgvm(ji,jj) - e1v(ji,jj-1)*spgvm(ji,jj-1) )
229         END DO
230      END DO
231     
232      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
233     
234      ! relative PRECISION
235      rnorme = 0.
236      DO jj = 1, jpj
237         DO ji = 1, jpi
238            rnorme = rnorme + gcb(ji,jj) * gcsmat(ji,jj) * gcb(ji,jj)
239         END DO
240      END DO
241      IF( lk_mpp )   CALL  mpp_sum( rnorme )   ! sum over the global domain
242
243      epsr=eps*eps*rnorme
244      ncut=0
245      !   IF the second member is 0 the solution is 0, solpcg isn't called
246      IF ( rnorme == 0.e0 ) THEN
247         gps(:,:) = 0.e0
248         res   = 0.e0
249         niter = 0
250         ncut  = 999
251      ENDIF
252     
253      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
254     
255      nindic = 0
256
257      ! iterarive solver of the spg system (except IF sol.=0)
258      !     (OUTPUT in gcx with boundary conditions applied)
259      IF ( ncut == 0 ) THEN
260         IF ( nsolv == 1 ) THEN
261            CALL sol_pcg( nindic )         !   diagonal preconditioned conjuguate gradient
262         ELSE IF ( nsolv == 2 ) THEN
263            CALL sol_sor( nindic )     !   successive-over-relaxation
264         ELSE IF(nsolv == 3) THEN
265            CALL sol_fet( nindic )         !   FETI solver
266         ELSE
267            !   e r r o r  in nsolv namelist PARAMETER
268            IF(lwp) THEN
269               WRITE(numout,*) ' dia_spr: e r r o r, nsolv = 1 or 2'
270               WRITE(numout,*) ' ******               not = ',nsolv
271            ENDIF
272            STOP 'dia_spr'
273         ENDIF
274      ENDIF
275     
276      !,,,,,,,,,,,,,,,,,,,,,,,,synchro IF macrotasking,,,,,,,,,,,,,,,,,,,,,,,
277
278     
279      ! sp solver statistics  (i.e. problem for the solver)
280      IF ( epsr < 0.) THEN
281         IF(lwp) THEN
282            WRITE(numout,*)'rrrrrrrrrrrrrrr'
283            IF(lwp)WRITE(numout,*)'dia_spr-1:',epsr
284            IF(lwp)WRITE(numout,*)'rrrrrrrrrrrrrrr'
285         ENDIF
286      ENDIF
287      IF(lwp)WRITE(numout,9300) kt, niter, res, SQRT(epsr)/eps
288      IF (nindic < 0) THEN
289         IF(lwp) THEN
290            WRITE(numout,9100)
291            WRITE(numout,*) ' dia_spr : the surface pressure solver DO not converge'
292            WRITE(numout,*) ' ====== ' 
293            WRITE(numout,*) 
294         ENDIF
295      ENDIF
2969100  FORMAT( /,' ===>>>> : w a r n i n g',/,'          ===============',/ )
2979300  FORMAT(' it :', i8, ' niter :', i4, ' res :',e20.10,' b :', e20.10)
298
299      ! recover bsf solver informations and SAVE ps for next computation
300      eps    = zeps
301      nmax   = imax 
302      gps   (:,:  ) = gcx (:,:)
303      gcp   (:,:,1) = ztab(:,:,1)
304      gcp   (:,:,2) = ztab(:,:,2)
305      gcp   (:,:,3) = ztab(:,:,3)
306      gcp   (:,:,4) = ztab(:,:,4)
307      gcdprc(:,:  ) = ztab(:,:,5)
308      gcdmat(:,:  ) = ztab(:,:,6)
309      gcx   (:,:  ) = ztab(:,:,7)
310      bmask (:,:  ) = ztab(:,:,8)
311     
312      ! compute and substract the mean value
313     
314      zpsmea = 0.e0
315      DO jj=2,jpjm1
316         DO ji=2,jpim1
317            zpsmea = zpsmea + gps(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1)
318         END DO
319      END DO
320      IF( lk_mpp )   CALL  mpp_sum( zpsmea )   ! sum over the global domain
321
322      zpsmea = zpsmea / e1e2t
323      gps(:,:) = ( gps(:,:) - zpsmea ) * tmask(:,:,1)
324 
325      IF(lwp)WRITE(numout,*) ' mean value of ps = ',zpsmea,' is substracted'
326      ! ----------------------------------------
327      ! i. compute the surface pressure gradient
328      !    from the computed surface pressure
329      ! ----------------------------------------
330
331      DO jj=2,jpjm1
332         DO ji=2,jpim1
333            gpsuu(ji,jj)=(gps(ji+1,jj)-gps(ji,jj))/e1u(ji,jj) * umask(ji,jj,1)
334            gpsvv(ji,jj)=(gps(ji,jj+1)-gps(ji,jj))/e2v(ji,jj) * vmask(ji,jj,1)
335         END DO
336      END DO
337     
338      ! compute the max and min error
339     
340      zemax1 = 0.e0
341      zemin1 = 0.e0
342      zemax2 = 0.e0
343      zemin2 = 0.e0
344      DO jj = 2,jpj-1
345         DO ji = 2,jpi-1
346            z1 = ABS( spgum(ji,jj)-gpsuu(ji,jj) )*umask(ji,jj,1)
347            z2 = ABS( spgvm(ji,jj)-gpsvv(ji,jj) )*vmask(ji,jj,1)
348            z3 = MAX ( ABS( spgum(ji,jj) ), ABS( spgvm(ji,jj) ) )
349            z4 = MAX ( ABS( gpsuu(ji,jj) ), ABS( gpsvv(ji,jj) ) )
350            zemax1 = MAX(z1,zemax1)
351            zemax2 = MAX(z2,zemax2)
352            zemin1 = MAX(z3,zemin1)
353            zemin2 = MAX(z4,zemin2)
354         END DO
355      END DO
356      IF( lk_mpp )   CALL  mpp_sum( zemax1 )   ! sum over the global domain
357      IF( lk_mpp )   CALL  mpp_sum( zemax2 )   ! sum over the global domain
358      IF( lk_mpp )   CALL  mpp_sum( zemin1 )   ! sum over the global domain
359      IF( lk_mpp )   CALL  mpp_sum( zemin2 )   ! sum over the global domain
360
361      IF(lwp) THEN
362         WRITE(numout,*)
363         WRITE(numout,*) 'pserro : time step = ', kt
364         WRITE(numout,*) '******** ------------------'
365         WRITE(numout,*)
366         WRITE(numout,*) '         gpsx error  max=',zemax1
367         WRITE(numout,*) '         gpsy error  max=',zemax2
368         WRITE(numout,*) '         gps max =',zemin1
369         WRITE(numout,*) '         gpsc max =',zemin2
370         WRITE(numout,*)
371      ENDIF
372
373      ! compute the norme and variance of this error
374
375      zcompt = 0.e0
376      zdif1  = 0.e0
377      zdif2  = 0.e0
378      zvar1  = 0.e0
379      zvar2  = 0.e0
380      DO jj = 2, jpj-1
381         DO ji = 2, jpi-1
382            z1 = ( spgum(ji,jj)-gpsuu(ji,jj) ) * umask(ji,jj,1)
383            z2 = ( spgvm(ji,jj)-gpsvv(ji,jj) ) * vmask(ji,jj,1)
384            zcompt=zcompt+tmask(ji,jj,1)
385            zdif1=zdif1+z1
386            zdif2=zdif2+z2
387            zvar1=zvar1+z1*z1
388            zvar2=zvar2+z2*z2
389         END DO
390      END DO
391      IF( lk_mpp )   CALL  mpp_sum( zcompt )   ! sum over the global domain
392      IF( lk_mpp )   CALL  mpp_sum( zdif1  )   ! sum over the global domain
393      IF( lk_mpp )   CALL  mpp_sum( zdif2  )   ! sum over the global domain
394      IF( lk_mpp )   CALL  mpp_sum( zvar1  )   ! sum over the global domain
395      IF( lk_mpp )   CALL  mpp_sum( zvar2  )   ! sum over the global domain
396
397      IF(lwp) WRITE(numout,*) '        zcompt = ',zcompt
398      zdif1=zdif1/zcompt
399      zdif2=zdif2/zcompt
400      IF( zvar1 < 0.) THEN
401         IF(lwp) THEN
402            WRITE(numout,*)'rrrrrrrrrrrrrrr'
403            WRITE(numout,*)'dia_spr-2:',zvar1
404            WRITE(numout,*)'rrrrrrrrrrrrrrr'
405         ENDIF
406      ENDIF
407      zvar1 = SQRT(zvar1)/zcompt
408      IF( zvar2 < 0. ) THEN
409         IF(lwp)THEN
410            WRITE(numout,*)'rrrrrrrrrrrrrrr'
411            WRITE(numout,*)'dia_spr-3:',zvar2
412            WRITE(numout,*)'rrrrrrrrrrrrrrr'
413         ENDIF
414      ENDIF
415      zvar2 = SQRT(zvar2)/zcompt
416     
417      IF(lwp) THEN
418         WRITE(numout,*)
419         WRITE(numout,*) '         gpsx mean error = ',zdif1
420         WRITE(numout,*) '         gpsy mean error = ',zdif2
421         WRITE(numout,*)
422         WRITE(numout,*) '         gpsx var. error = ',zvar1
423         WRITE(numout,*) '         gpsy var. error = ',zvar2
424         WRITE(numout,*)
425         WRITE(numout,*)
426      ENDIF
427     
428      ! reset to zero nmoyps and the mean surface pressure gradient
429      nmoyps = 0
430      spgum(:,:) = 0.e0
431      spgvm(:,:) = 0.e0
432     
433   END SUBROUTINE dia_spr
434
435
436   SUBROUTINE sprmat
437      !!---------------------------------------------------------------------
438      !!                  ***  ROUTINE sprmat  ***
439      !!               
440      !! ** Purpose :   construction of the matrix of the surface pressure
441      !!      system and the diagonal preconditioning matrix.
442      !!
443      !! ** Method :
444      !!
445      !! History :
446      !!        !  98-01  (G. Madec & M. Ioualalen)  Original code
447      !!   8.5  !  02-08  (G. Madec)  F90: Free form
448      !!----------------------------------------------------------------------
449      !! * Local declarations
450      INTEGER ::   ji, jj, jl                ! dummy loop indices
451      REAL(wp) ::   zcoefs, zcoefw, zcoefe, zcoefn
452      !!----------------------------------------------------------------------
453     
454     
455      ! 0. ocean/land mask at ps-point (computed from surface tmask): spmsk
456      ! --------------------------------------------------------------------
457     
458      ! computation
459      spmsk(:,:) = tmask(:,:,1)
460     
461      ! boundary conditions
462      ! south symmetry: psmsk must be set to 0. on 1
463      IF( nperio == 2 ) THEN
464         spmsk(:, 1 ) = 0.e0
465      ENDIF
466     
467      ! east-west cyclic: spmsk must be set to 0. on 1 and jpi
468      IF( nperio == 1 .OR. nperio == 4 .OR.nperio == 6) THEN
469         spmsk( 1 ,:) = 0.e0
470         spmsk(jpi,:) = 0.e0
471      ENDIF
472     
473      ! north fold: spmsk must be set to 0. on ligne jpj and on half
474      !                   ligne jpj-1
475      ! T-point pivot
476      IF( nperio == 3 .OR. nperio == 4 ) THEN
477         spmsk(:,jpj) = 0.e0
478         DO ji = jpi/2+1, jpi
479            spmsk(ji,jpjm1) = 0.e0
480         END DO
481      ENDIF
482      ! F-point pivot
483      IF( nperio == 5 .OR. nperio == 6 ) THEN
484         spmsk(:,jpj) = 0.e0
485      ENDIF
486     
487      ! mpp boundary cond.: spmsk is initialized at zero on the overlap
488      ! region for both the preconjugate gradient and the sor algorithms
489
490      IF( nbondi /= -1 .AND. nbondi /= 2 ) THEN
491         DO jl = 1, jpreci
492            spmsk(jl,:) = 0.e0
493         END DO
494      ENDIF
495      IF( nbondi /= 1 .AND. nbondi /= 2 ) THEN
496         DO ji = nlci, jpi
497            spmsk(ji,:) = 0.e0
498         END DO
499      ENDIF
500      IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
501         DO jl=1,jprecj
502            spmsk(:,jl) = 0.e0
503         END DO
504      ENDIF
505      IF( nbondj /= 1 .AND. nbondj /= 2 ) THEN
506         DO jj = nlcj, jpj
507            spmsk(:,jj) = 0.e0
508         END DO
509      ENDIF
510     
511      ! 1. construction of the matrix
512      ! -----------------------------
513     
514      DO jj = 1, jpj
515         DO ji = 1, jpi
516           
517            IF( spmsk(ji,jj) == 0. ) THEN
518               ! land points
519               gcps  (ji,jj,1) = 0.e0
520               gcps  (ji,jj,2) = 0.e0
521               gcps  (ji,jj,3) = 0.e0
522               gcps  (ji,jj,4) = 0.e0
523               gcdpsc(ji,jj  ) = 0.e0
524               gcsmat(ji,jj  ) = 0.e0
525            ELSE
526               ! south coefficient
527               zcoefs = -e1v(ji,jj-1) / e2v(ji,jj-1) * vmask(ji,jj-1,1)
528               gcps(ji,jj,1) = zcoefs
529               ! west coefficient
530               zcoefw = -e2u(ji-1,jj) / e1u(ji-1,jj) * umask(ji-1,jj,1)
531               gcps(ji,jj,2) = zcoefw
532               ! east coefficient
533               zcoefe = -e2u(ji  ,jj) / e1u(ji  ,jj) * umask(ji  ,jj,1)
534               gcps(ji,jj,3) = zcoefe
535               ! north coefficient
536               zcoefn = -e1v(ji,jj  ) / e2v(ji,jj  ) * vmask(ji,jj  ,1)
537               gcps(ji,jj,4) = zcoefn
538               
539               ! diagonal coefficient
540               gcsmat(ji,jj) = -zcoefs-zcoefw-zcoefe-zcoefn
541            ENDIF
542         END DO
543      END DO
544     
545
546      ! 2. boundary conditions
547      ! ----------------------
548     
549      ! cyclic east-west boundary conditions
550      ! ji=2 is the column east of ji=jpim1 and reciprocally,
551      ! ji=jpim1 is the column west of ji=2
552      ! all the coef are already set to zero as spmask is initialized to
553      ! zero for ji=1 and ji=jpj.
554     
555      ! symetrical conditions
556      ! the diagonal coefficient of the southern grid points must be modify to
557      ! account for the existence of the south symmetric bassin.
558      IF( nperio == 2 ) THEN
559         DO ji = 1, jpi
560            IF( spmsk(ji,2) /= 0 ) THEN
561               zcoefs = e1v(ji,1) / e2v(ji,1)
562               gcsmat(ji,2) = gcsmat(ji,2) - zcoefs
563            ENDIF
564         END DO
565      ENDIF
566     
567      ! North fold boundary condition
568      ! all the coef are already set to zero as bmask is initialized to
569      ! zero on duplicated lignes and portion of lignes
570     
571     
572      ! 3. preconditioned matrix
573      ! ------------------------
574     
575      DO jj = 1, jpj
576         DO ji = 1, jpi
577            IF( spmsk(ji,jj) /= 0. ) gcdpsc(ji,jj) = 1.e0 / gcsmat(ji,jj)
578         END DO
579      END DO
580     
581      gcps(:,:,1) = gcps(:,:,1) * gcdpsc(:,:)
582      gcps(:,:,2) = gcps(:,:,2) * gcdpsc(:,:)
583      gcps(:,:,3) = gcps(:,:,3) * gcdpsc(:,:)
584      gcps(:,:,4) = gcps(:,:,4) * gcdpsc(:,:)
585     
586     
587      ! 3. initialization the arrays used in sp solver
588      ! ----------------------------------------------
589     
590      gps  (:,:) = 0.e0
591      gpsuu(:,:) = 0.e0
592      gpsvv(:,:) = 0.e0
593     
594   END SUBROUTINE sprmat
595
596#else
597   !!----------------------------------------------------------------------
598   !!   Default option :                    NO surface pressure diagnostics
599   !!----------------------------------------------------------------------
600   LOGICAL, PUBLIC, PARAMETER ::   lk_diaspr = .FALSE.   !: surface pressure diag. flag
601CONTAINS
602   SUBROUTINE dia_spr( kt )      ! Empty routine
603      WRITE(*,*) 'dia_spr: You should not have seen this print! error?', kt
604   END SUBROUTINE dia_spr
605#endif
606
607   !!======================================================================
608END MODULE diaspr
Note: See TracBrowser for help on using the repository browser.