[1] | 1 | ! |
---|
| 2 | ! $Header$ |
---|
| 3 | ! |
---|
| 4 | SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt) |
---|
[1422] | 5 | |
---|
| 6 | USE comconst_mod, ONLY: pi |
---|
| 7 | |
---|
[1] | 8 | IMPLICIT NONE |
---|
| 9 | |
---|
| 10 | c======================================================================= |
---|
| 11 | c Adaptation LMDZ: A.Armengaud (LGGE) |
---|
| 12 | c ---------------- |
---|
| 13 | c |
---|
| 14 | c ************************************************ |
---|
| 15 | c Transport des traceurs par la methode de prather |
---|
| 16 | c Ref : |
---|
| 17 | c |
---|
| 18 | c ************************************************ |
---|
| 19 | c q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg |
---|
| 20 | c |
---|
| 21 | c======================================================================= |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | #include "dimensions.h" |
---|
| 25 | #include "paramet.h" |
---|
| 26 | #include "comgeom2.h" |
---|
| 27 | |
---|
| 28 | c Arguments: |
---|
| 29 | c ---------- |
---|
| 30 | INTEGER iq,nt |
---|
| 31 | REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) |
---|
| 32 | REAL masse(iip1,jjp1,llm) |
---|
| 33 | REAL q( iip1,jjp1,llm,0:9) |
---|
| 34 | REAL w( ip1jmp1,llm ) |
---|
| 35 | integer ordre,ilim |
---|
| 36 | |
---|
| 37 | c Local: |
---|
| 38 | c ------ |
---|
| 39 | LOGICAL limit |
---|
| 40 | real zq(iip1,jjp1,llm) |
---|
| 41 | REAL sm ( iip1,jjp1, llm ) |
---|
| 42 | REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) |
---|
| 43 | REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm ) |
---|
| 44 | REAL sxx( iip1,jjp1,llm) |
---|
| 45 | REAL sxy( iip1,jjp1,llm) |
---|
| 46 | REAL sxz( iip1,jjp1,llm) |
---|
| 47 | REAL syy( iip1,jjp1,llm ) |
---|
| 48 | REAL syz( iip1,jjp1,llm ) |
---|
| 49 | REAL szz( iip1,jjp1,llm ),zz |
---|
| 50 | INTEGER i,j,l,indice |
---|
| 51 | real sxn(iip1),sxs(iip1) |
---|
| 52 | |
---|
| 53 | real sinlon(iip1),sinlondlon(iip1) |
---|
| 54 | real coslon(iip1),coslondlon(iip1) |
---|
| 55 | real qmin,qmax |
---|
| 56 | save qmin,qmax |
---|
| 57 | save sinlon,coslon,sinlondlon,coslondlon |
---|
| 58 | real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps |
---|
| 59 | real masn,mass |
---|
| 60 | c |
---|
| 61 | REAL SSUM |
---|
| 62 | integer ismax,ismin |
---|
| 63 | EXTERNAL SSUM, convflu,ismin,ismax |
---|
| 64 | logical first |
---|
| 65 | save first |
---|
| 66 | EXTERNAL advxp,advyp,advzp |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | data first/.true./ |
---|
| 70 | data qmin,qmax/-1.e33,1.e33/ |
---|
| 71 | |
---|
| 72 | |
---|
| 73 | c========================================================================== |
---|
| 74 | c========================================================================== |
---|
| 75 | c MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt |
---|
| 76 | c========================================================================== |
---|
| 77 | c========================================================================== |
---|
| 78 | REAL dt |
---|
| 79 | c========================================================================== |
---|
| 80 | limit = .TRUE. |
---|
| 81 | |
---|
| 82 | if(first) then |
---|
| 83 | print*,'SCHEMA PRATHER' |
---|
| 84 | first=.false. |
---|
| 85 | do i=2,iip1 |
---|
| 86 | coslon(i)=cos(rlonv(i)) |
---|
| 87 | sinlon(i)=sin(rlonv(i)) |
---|
| 88 | coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi |
---|
| 89 | sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi |
---|
| 90 | enddo |
---|
| 91 | coslon(1)=coslon(iip1) |
---|
| 92 | coslondlon(1)=coslondlon(iip1) |
---|
| 93 | sinlon(1)=sinlon(iip1) |
---|
| 94 | sinlondlon(1)=sinlondlon(iip1) |
---|
| 95 | |
---|
| 96 | DO l = 1,llm |
---|
| 97 | DO j = 1,jjp1 |
---|
| 98 | DO i = 1,iip1 |
---|
| 99 | q( i,j,l,1 )=0. |
---|
| 100 | q( i,j,l,2)=0. |
---|
| 101 | q( i,j,l,3)=0. |
---|
| 102 | q( i,j,l,4)=0. |
---|
| 103 | q( i,j,l,5)=0. |
---|
| 104 | q( i,j,l,6)=0. |
---|
| 105 | q( i,j,l,7)=0. |
---|
| 106 | q( i,j,l,8)=0. |
---|
| 107 | q( i,j,l,9)=0. |
---|
| 108 | ENDDO |
---|
| 109 | ENDDO |
---|
| 110 | ENDDO |
---|
| 111 | endif |
---|
| 112 | c Fin modif Fred |
---|
| 113 | |
---|
| 114 | c *** On calcule la masse d'air en kg |
---|
| 115 | |
---|
| 116 | DO l = 1,llm |
---|
| 117 | DO j = 1,jjp1 |
---|
| 118 | DO i = 1,iip1 |
---|
| 119 | sm( i,j,llm+1-l ) =masse(i,j,l) |
---|
| 120 | ENDDO |
---|
| 121 | ENDDO |
---|
| 122 | ENDDO |
---|
| 123 | |
---|
| 124 | c *** q contient les qqtes de traceur avant l'advection |
---|
| 125 | |
---|
| 126 | c *** Affectation des tableaux S a partir de Q |
---|
| 127 | |
---|
| 128 | DO l = 1,llm |
---|
| 129 | DO j = 1,jjp1 |
---|
| 130 | DO i = 1,iip1 |
---|
| 131 | s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l) |
---|
| 132 | sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l) |
---|
| 133 | sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l) |
---|
| 134 | sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l) |
---|
| 135 | sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l) |
---|
| 136 | sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l) |
---|
| 137 | sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l) |
---|
| 138 | syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l) |
---|
| 139 | syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l) |
---|
| 140 | szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l) |
---|
| 141 | ENDDO |
---|
| 142 | ENDDO |
---|
| 143 | ENDDO |
---|
| 144 | c *** Appel des subroutines d'advection en X, en Y et en Z |
---|
| 145 | c *** Advection avec "time-splitting" |
---|
| 146 | |
---|
| 147 | c----------------------------------------------------------- |
---|
| 148 | do indice =1,nt |
---|
| 149 | call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz |
---|
| 150 | . ,sxx,sxy,sxz,syy,syz,szz,1 ) |
---|
| 151 | end do |
---|
| 152 | do l=1,llm |
---|
| 153 | do i=1,iip1 |
---|
| 154 | sy(i,1,l)=0. |
---|
| 155 | sy(i,jjp1,l)=0. |
---|
| 156 | enddo |
---|
| 157 | enddo |
---|
| 158 | c--------------------------------------------------------- |
---|
| 159 | call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz |
---|
| 160 | . ,sxx,sxy,sxz,syy,syz,szz,1 ) |
---|
| 161 | c--------------------------------------------------------- |
---|
| 162 | |
---|
| 163 | c--------------------------------------------------------- |
---|
| 164 | do j=1,jjp1 |
---|
| 165 | do i=1,iip1 |
---|
| 166 | sz(i,j,1)=0. |
---|
| 167 | sz(i,j,llm)=0. |
---|
| 168 | sxz(i,j,1)=0. |
---|
| 169 | sxz(i,j,llm)=0. |
---|
| 170 | syz(i,j,1)=0. |
---|
| 171 | syz(i,j,llm)=0. |
---|
| 172 | szz(i,j,1)=0. |
---|
| 173 | szz(i,j,llm)=0. |
---|
| 174 | enddo |
---|
| 175 | enddo |
---|
| 176 | call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz |
---|
| 177 | . ,sxx,sxy,sxz,syy,syz,szz,1 ) |
---|
| 178 | do l=1,llm |
---|
| 179 | do i=1,iip1 |
---|
| 180 | sy(i,1,l)=0. |
---|
| 181 | sy(i,jjp1,l)=0. |
---|
| 182 | enddo |
---|
| 183 | enddo |
---|
| 184 | |
---|
| 185 | c--------------------------------------------------------- |
---|
| 186 | |
---|
| 187 | c--------------------------------------------------------- |
---|
| 188 | call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz |
---|
| 189 | . ,sxx,sxy,sxz,syy,syz,szz,1 ) |
---|
| 190 | c--------------------------------------------------------- |
---|
| 191 | DO l = 1,llm |
---|
| 192 | DO j = 1,jjp1 |
---|
| 193 | s0( iip1,j,l)=s0( 1,j,l ) |
---|
| 194 | sx( iip1,j,l)=sx( 1,j,l ) |
---|
| 195 | sy( iip1,j,l)=sy( 1,j,l ) |
---|
| 196 | sz( iip1,j,l)=sz( 1,j,l ) |
---|
| 197 | sxx( iip1,j,l)=sxx( 1,j,l ) |
---|
| 198 | sxy( iip1,j,l)=sxy( 1,j,l) |
---|
| 199 | sxz( iip1,j,l)=sxz( 1,j,l ) |
---|
| 200 | syy( iip1,j,l)=syy( 1,j,l ) |
---|
| 201 | syz( iip1,j,l)=syz( 1,j,l) |
---|
| 202 | szz( iip1,j,l)=szz( 1,j,l ) |
---|
| 203 | ENDDO |
---|
| 204 | ENDDO |
---|
| 205 | do indice=1,nt |
---|
| 206 | call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz |
---|
| 207 | . ,sxx,sxy,sxz,syy,syz,szz,1 ) |
---|
| 208 | end do |
---|
| 209 | c--------------------------------------------------------- |
---|
| 210 | c--------------------------------------------------------- |
---|
| 211 | c *** On repasse les S dans la variable qpr |
---|
| 212 | c *** On repasse les S dans la variable q directement 14/10/94 |
---|
| 213 | |
---|
| 214 | DO l = 1,llm |
---|
| 215 | DO j = 1,jjp1 |
---|
| 216 | DO i = 1,iip1 |
---|
| 217 | q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l) |
---|
| 218 | q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l) |
---|
| 219 | q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l) |
---|
| 220 | q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l) |
---|
| 221 | q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l) |
---|
| 222 | q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l) |
---|
| 223 | q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l) |
---|
| 224 | q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l) |
---|
| 225 | q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l) |
---|
| 226 | q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l) |
---|
| 227 | ENDDO |
---|
| 228 | ENDDO |
---|
| 229 | ENDDO |
---|
| 230 | |
---|
| 231 | c--------------------------------------------------------- |
---|
| 232 | c go to 777 |
---|
| 233 | c filtrages aux poles |
---|
| 234 | |
---|
| 235 | c Traitements specifiques au pole |
---|
| 236 | |
---|
| 237 | c filtrages aux poles |
---|
| 238 | DO l=1,llm |
---|
| 239 | c filtrages aux poles |
---|
| 240 | masn=ssum(iim,sm(1,1,l),1) |
---|
| 241 | mass=ssum(iim,sm(1,jjp1,l),1) |
---|
| 242 | qpn=ssum(iim,s0(1,1,l),1)/masn |
---|
| 243 | qps=ssum(iim,s0(1,jjp1,l),1)/mass |
---|
| 244 | dqzpn=ssum(iim,sz(1,1,l),1)/masn |
---|
| 245 | dqzps=ssum(iim,sz(1,jjp1,l),1)/mass |
---|
| 246 | do i=1,iip1 |
---|
| 247 | q( i,1,llm+1-l,3)=dqzpn |
---|
| 248 | q( i,jjp1,llm+1-l,3)=dqzps |
---|
| 249 | q( i,1,llm+1-l,0)=qpn |
---|
| 250 | q( i,jjp1,llm+1-l,0)=qps |
---|
| 251 | enddo |
---|
| 252 | c enddo |
---|
| 253 | c print*,'qpn',qpn,'qps',qps |
---|
| 254 | c print*,'dqzpn',dqzpn,'dqzps',dqzps |
---|
| 255 | c enddo |
---|
| 256 | dyn1=0. |
---|
| 257 | dys1=0. |
---|
| 258 | dyn2=0. |
---|
| 259 | dys2=0. |
---|
| 260 | do i=1,iim |
---|
| 261 | zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0) |
---|
| 262 | dyn1=dyn1+sinlondlon(i)*zz |
---|
| 263 | dyn2=dyn2+coslondlon(i)*zz |
---|
| 264 | zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l) |
---|
| 265 | dys1=dys1+sinlondlon(i)*zz |
---|
| 266 | dys2=dys2+coslondlon(i)*zz |
---|
| 267 | enddo |
---|
| 268 | do i=1,iim |
---|
| 269 | q(i,1,llm+1-l,2)= |
---|
| 270 | $ (sinlon(i)*dyn1+coslon(i)*dyn2)/2. |
---|
| 271 | q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0) |
---|
| 272 | $ +q(i,1,llm+1-l,2) |
---|
| 273 | q(i,jjp1,llm+1-l,2)= |
---|
| 274 | $ (sinlon(i)*dys1+coslon(i)*dys2)/2. |
---|
| 275 | q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) |
---|
| 276 | $ -q(i,jjp1,llm+1-l,2) |
---|
| 277 | enddo |
---|
| 278 | q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0) |
---|
| 279 | q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0) |
---|
| 280 | do i=1,iim |
---|
| 281 | sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0) |
---|
| 282 | sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0) |
---|
| 283 | enddo |
---|
| 284 | sxn(iip1)=sxn(1) |
---|
| 285 | sxs(iip1)=sxs(1) |
---|
| 286 | do i=1,iim |
---|
| 287 | q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1)) |
---|
| 288 | q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1)) |
---|
| 289 | END DO |
---|
| 290 | q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1) |
---|
| 291 | q(1,jjp1,llm+1-l,1)= |
---|
| 292 | $ q(iip1,jjp1,llm+1-l,1) |
---|
| 293 | enddo |
---|
| 294 | do l=1,llm |
---|
| 295 | do i=1,iim |
---|
| 296 | q( i,1,llm+1-l,4)=0. |
---|
| 297 | q( i,jjp1,llm+1-l,4)=0. |
---|
| 298 | q( i,1,llm+1-l,5)=0. |
---|
| 299 | q( i,jjp1,llm+1-l,5)=0. |
---|
| 300 | q( i,1,llm+1-l,6)=0. |
---|
| 301 | q( i,jjp1,llm+1-l,6)=0. |
---|
| 302 | q( i,1,llm+1-l,7)=0. |
---|
| 303 | q( i,jjp1,llm+1-l,7)=0. |
---|
| 304 | q( i,1,llm+1-l,8)=0. |
---|
| 305 | q( i,jjp1,llm+1-l,8)=0. |
---|
| 306 | q( i,1,llm+1-l,9)=0. |
---|
| 307 | q( i,jjp1,llm+1-l,9)=0. |
---|
| 308 | enddo |
---|
| 309 | ENDDO |
---|
| 310 | |
---|
| 311 | 777 continue |
---|
| 312 | c |
---|
| 313 | c bouclage en longitude |
---|
| 314 | do l=1,llm |
---|
| 315 | do j=1,jjp1 |
---|
| 316 | q(iip1,j,l,0)=q(1,j,l,0) |
---|
| 317 | q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0) |
---|
| 318 | q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1) |
---|
| 319 | q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2) |
---|
| 320 | q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3) |
---|
| 321 | q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4) |
---|
| 322 | q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5) |
---|
| 323 | q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6) |
---|
| 324 | q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7) |
---|
| 325 | q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8) |
---|
| 326 | q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9) |
---|
| 327 | enddo |
---|
| 328 | enddo |
---|
| 329 | DO l = 1,llm |
---|
| 330 | DO j = 2,jjm |
---|
| 331 | DO i = 1,iip1 |
---|
| 332 | IF (q(i,j,l,0).lt.0.) THEN |
---|
| 333 | PRINT*,'------------ BIP-----------' |
---|
| 334 | PRINT*,'S0(',i,j,l,')=',q(i,j,l,0), |
---|
| 335 | $ q(i,j-1,l,0) |
---|
| 336 | PRINT*,'SX(',i,j,l,')=',q(i,j,l,1) |
---|
| 337 | PRINT*,'SY(',i,j,l,')=',q(i,j,l,2), |
---|
| 338 | $ q(i,j-1,l,2) |
---|
| 339 | PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3) |
---|
| 340 | c PRINT*,' PBL EN SORTIE D'' ADVZP' |
---|
| 341 | q(i,j,l,0)=0. |
---|
| 342 | c STOP |
---|
| 343 | ENDIF |
---|
| 344 | ENDDO |
---|
| 345 | ENDDO |
---|
| 346 | do j=1,jjp1,jjm |
---|
| 347 | do i=1,iip1 |
---|
| 348 | IF (q(i,j,l,0).lt.0.) THEN |
---|
| 349 | PRINT*,'------------ BIP 2-----------' |
---|
| 350 | PRINT*,'S0(',i,j,l,')=',q(i,j,l,0) |
---|
| 351 | PRINT*,'SX(',i,j,l,')=',q(i,j,l,1) |
---|
| 352 | PRINT*,'SY(',i,j,l,')=',q(i,j,l,2) |
---|
| 353 | PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3) |
---|
| 354 | |
---|
| 355 | q(i,j,l,0)=0. |
---|
| 356 | c STOP |
---|
| 357 | ENDIF |
---|
| 358 | enddo |
---|
| 359 | enddo |
---|
| 360 | ENDDO |
---|
| 361 | RETURN |
---|
| 362 | END |
---|