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

Last change on this file since 412 was 359, checked in by opalod, 18 years ago

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