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