1 | MODULE 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 | !! $Id$ |
---|
65 | !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt |
---|
66 | !!---------------------------------------------------------------------- |
---|
67 | |
---|
68 | CONTAINS |
---|
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 |
---|
295 | 9100 FORMAT( /,' ===>>>> : w a r n i n g',/,' ===============',/ ) |
---|
296 | 9300 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 |
---|
600 | CONTAINS |
---|
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 | !!====================================================================== |
---|
607 | END MODULE diaspr |
---|