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