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