Changeset 495 for LMDZ.3.3/branches/rel-LF/libf
- Timestamp:
- Mar 4, 2004, 4:11:16 PM (21 years ago)
- Location:
- LMDZ.3.3/branches/rel-LF/libf
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F
r477 r495 153 153 LOGICAL ok_nudge 154 154 PARAMETER (ok_nudge = .false.) 155 c 155 156 156 157 c----------------------------------------------------------------------- … … 367 368 c ---------------------------------- 368 369 370 369 371 1 CONTINUE 370 372 … … 449 451 ELSE IF( iq.EQ. nqmx ) THEN 450 452 c 451 iapp_tracvl = 5 453 c iapp_tracvl = 5 454 iapp_tracvl = iperiod 455 print*,'***WARNING***: iapp_tracvl = iperiod' 452 456 c 453 457 cccc iapp_tracvl est la frequence en pas du groupement des flux 454 458 cccc de masse pour Van-Leer dans la routine tracvl . 455 459 c 456 CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv, 460 c CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv, 461 CALL vanleer(numvanle,iapp_tracvl,2,q,pbaru,pbarv, 457 462 * p, masse, dq, iadv(1), teta, pk ) 458 463 c 464 print*,'***WARNING***: Appel vanleer avec 2 traceurs' 465 print*,' et non nqmx' 459 466 c ... Modif F.Codron .... 460 467 c … … 483 490 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis , 484 491 $ finvmaold ) 492 485 493 486 494 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd) … … 660 668 IF( itau. EQ. itaufinp1 ) then 661 669 670 662 671 abort_message = 'Simulation finished' 663 672 call abort_gcm(modname,abort_message,0) -
LMDZ.3.3/branches/rel-LF/libf/filtrez/filtreg.F
r231 r495 49 49 INTEGER i,j,l,k 50 50 INTEGER iim2,immjm 51 INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil ,jffil51 INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil(2),jffil(2) 52 52 53 53 REAL champ( iip1,nlat,nbniv) … … 56 56 , , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs) 57 57 , , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus) 58 REAL eignq(iim), sdd1(iim),sdd2(iim) 58 cIM REAL eignq(iim), sdd1(iim),sdd2(iim) 59 REAL eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim) 59 60 LOGICAL griscal 60 61 INTEGER hemisph, iaire … … 127 128 END IF 128 129 c 130 jdfil(1) = jdfil1 131 jffil(1) = jffil1 132 jdfil(2) = jdfil2 133 jffil(2) = jffil2 129 134 c 130 135 DO 100 hemisph = 1, 2 131 136 c 132 133 134 135 136 137 138 137 c IF ( hemisph.EQ.1 ) THEN 138 c jdfil = jdfil1 139 c jffil = jffil1 140 c ELSE 141 c jdfil = jdfil2 142 c jffil = jffil2 143 c END IF 139 144 140 145 141 146 DO 50 l = 1, nbniv 142 DO 30 j = jdfil ,jffil147 DO 30 j = jdfil(hemisph),jffil(hemisph) 143 148 144 149 … … 147 152 5 CONTINUE 148 153 c 154 30 CONTINUE 155 50 CONTINUE 149 156 150 157 IF( hemisph. EQ. 1 ) THEN … … 152 159 IF( ifiltre. EQ. -2 ) THEN 153 160 #ifdef CRAY 161 DO l = 1, nbniv 162 DO j = jdfil(hemisph),jffil(hemisph) 154 163 CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 155 164 * 1, iim, iim ) 156 #else 157 #ifdef BLAS 165 ENDDO 166 ENDDO 167 #else 168 #ifdef BLAS 169 DO l = 1, nbniv 170 DO j = jdfil(hemisph),jffil(hemisph) 158 171 CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim, 159 172 . champ(1,j,l), 1, 0.0, eignq, 1) 160 #else 161 DO k = 1, iim 162 eignq(k) = 0.0 163 ENDDO 164 DO k = 1, iim 165 DO i = 1, iim 166 eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l) 167 ENDDO 168 ENDDO 173 ENDDO 174 ENDDO 175 #else 176 DO l = 1, nbniv 177 DO j = jdfil(hemisph),jffil(hemisph) 178 DO k = 1, iim 179 c eignq(k) = 0.0 180 eignq(k,j,l) = 0.0 181 ENDDO 182 DO k = 1, iim 183 DO i = 1, iim 184 c eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l) 185 eignq(k,j,l) = eignq(k,j,l) + matrinvn(k,i,j)*champ(i,j,l) 186 ENDDO 187 ENDDO 188 ENDDO 189 ENDDO 169 190 #endif 170 191 #endif 171 192 ELSE IF ( griscal ) THEN 172 193 #ifdef CRAY 194 DO l = 1, nbniv 195 DO j = jdfil(hemisph),jffil(hemisph) 173 196 CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 174 197 * 1, iim, iim ) 175 #else 176 #ifdef BLAS 198 ENDDO 199 ENDDO 200 #else 201 #ifdef BLAS 202 DO l = 1, nbniv 203 DO j = jdfil(hemisph),jffil(hemisph) 177 204 CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim, 178 205 . champ(1,j,l), 1, 0.0, eignq, 1) 179 #else 180 DO k = 1, iim 181 eignq(k) = 0.0 182 ENDDO 183 DO i = 1, iim 184 DO k = 1, iim 185 eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l) 186 ENDDO 187 ENDDO 206 ENDDO 207 ENDDO 208 #else 209 DO l = 1, nbniv 210 DO j = jdfil(hemisph),jffil(hemisph) 211 DO k = 1, iim 212 c eignq(k) = 0.0 213 eignq(k,j,l) = 0.0 214 ENDDO 215 DO i = 1, iim 216 DO k = 1, iim 217 c eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l) 218 eignq(k,j,l) = eignq(k,j,l) + matriceun(k,i,j)*champ(i,j,l) 219 ENDDO 220 ENDDO 221 ENDDO 222 ENDDO 188 223 #endif 189 224 #endif 190 225 ELSE 191 226 #ifdef CRAY 227 DO l = 1, nbniv 228 DO j = jdfil(hemisph),jffil(hemisph) 192 229 CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 193 230 * 1, iim, iim ) 194 #else 195 #ifdef BLAS 231 ENDDO 232 ENDDO 233 #else 234 #ifdef BLAS 235 DO l = 1, nbniv 236 DO j = jdfil(hemisph),jffil(hemisph) 196 237 CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim, 197 238 . champ(1,j,l), 1, 0.0, eignq, 1) 198 #else 199 DO k = 1, iim 200 eignq(k) = 0.0 201 ENDDO 202 DO i = 1, iim 203 DO k = 1, iim 204 eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l) 205 ENDDO 206 ENDDO 239 ENDDO 240 ENDDO 241 #else 242 DO l = 1, nbniv 243 DO j = jdfil(hemisph),jffil(hemisph) 244 DO k = 1, iim 245 c eignq(k) = 0.0 246 eignq(k,j,l) = 0.0 247 ENDDO 248 DO i = 1, iim 249 DO k = 1, iim 250 c eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l) 251 eignq(k,j,l) = eignq(k,j,l) + matricevn(k,i,j)*champ(i,j,l) 252 ENDDO 253 ENDDO 254 ENDDO 255 ENDDO 207 256 #endif 208 257 #endif … … 213 262 IF( ifiltre. EQ. -2 ) THEN 214 263 #ifdef CRAY 264 DO l = 1, nbniv 265 DO j = jdfil(hemisph),jffil(hemisph) 215 266 CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , 216 267 * eignq, 1, iim, iim ) 217 #else 218 #ifdef BLAS 268 ENDDO 269 ENDDO 270 #else 271 #ifdef BLAS 272 DO l = 1, nbniv 273 DO j = jdfil(hemisph),jffil(hemisph) 219 274 CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim, 220 275 . champ(1,j,l), 1, 0.0, eignq, 1) 221 #else 222 DO k = 1, iim 223 eignq(k) = 0.0 224 ENDDO 225 DO i = 1, iim 226 DO k = 1, iim 227 eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) 228 ENDDO 229 ENDDO 276 ENDDO 277 ENDDO 278 #else 279 DO l = 1, nbniv 280 DO j = jdfil(hemisph),jffil(hemisph) 281 DO k = 1, iim 282 c eignq(k) = 0.0 283 eignq(k,j,l) = 0.0 284 ENDDO 285 DO i = 1, iim 286 DO k = 1, iim 287 c eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) 288 eignq(k,j,l) = eignq(k,j,l) + 289 .matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) 290 ENDDO 291 ENDDO 292 ENDDO 293 ENDDO 230 294 #endif 231 295 #endif 232 296 ELSE IF ( griscal ) THEN 233 297 #ifdef CRAY 298 DO l = 1, nbniv 299 DO j = jdfil(hemisph),jffil(hemisph) 234 300 CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , 235 301 * eignq, 1, iim, iim ) 236 #else 237 #ifdef BLAS 302 ENDDO 303 ENDDO 304 #else 305 #ifdef BLAS 306 DO l = 1, nbniv 307 DO j = jdfil(hemisph),jffil(hemisph) 238 308 CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim, 239 309 . champ(1,j,l), 1, 0.0, eignq, 1) 240 #else 241 DO k = 1, iim 242 eignq(k) = 0.0 243 ENDDO 244 DO i = 1, iim 245 DO k = 1, iim 246 eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) 247 ENDDO 248 ENDDO 310 ENDDO 311 ENDDO 312 #else 313 DO l = 1, nbniv 314 DO j = jdfil(hemisph),jffil(hemisph) 315 DO k = 1, iim 316 c eignq(k) = 0.0 317 eignq(k,j,l) = 0.0 318 ENDDO 319 DO i = 1, iim 320 DO k = 1, iim 321 c eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) 322 eignq(k,j,l) = eignq(k,j,l) + 323 .matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) 324 ENDDO 325 ENDDO 326 ENDDO 327 ENDDO 249 328 #endif 250 329 #endif 251 330 ELSE 252 331 #ifdef CRAY 332 DO l = 1, nbniv 333 DO j = jdfil(hemisph),jffil(hemisph) 253 334 CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 , 254 335 * eignq, 1, iim, iim ) 255 #else 256 #ifdef BLAS 336 ENDDO 337 ENDDO 338 #else 339 #ifdef BLAS 340 DO l = 1, nbniv 341 DO j = jdfil(hemisph),jffil(hemisph) 257 342 CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim, 258 343 . champ(1,j,l), 1, 0.0, eignq, 1) 259 #else 260 DO k = 1, iim 261 eignq(k) = 0.0 262 ENDDO 263 DO i = 1, iim 264 DO k = 1, iim 265 eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) 266 ENDDO 267 ENDDO 344 ENDDO 345 ENDDO 346 #else 347 DO l = 1, nbniv 348 DO j = jdfil(hemisph),jffil(hemisph) 349 DO k = 1, iim 350 c eignq(k) = 0.0 351 eignq(k,j,l) = 0.0 352 ENDDO 353 DO i = 1, iim 354 DO k = 1, iim 355 c eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) 356 eignq(k,j,l) = eignq(k,j,l) + 357 .matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) 358 ENDDO 359 ENDDO 360 ENDDO 361 ENDDO 268 362 #endif 269 363 #endif … … 273 367 c 274 368 IF( ifiltre.EQ. 2 ) THEN 369 DO l = 1, nbniv 370 DO j = jdfil(hemisph),jffil(hemisph) 275 371 DO 15 i = 1, iim 276 champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) 372 c champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) 373 champ( i,j,l ) = ( champ(i,j,l) + eignq(i,j,l) ) * sdd2(i) 277 374 15 CONTINUE 375 ENDDO 376 ENDDO 278 377 ELSE 378 DO l = 1, nbniv 379 DO j = jdfil(hemisph),jffil(hemisph) 279 380 DO 16 i=1,iim 280 champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) 381 c champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) 382 champ( i,j,l ) = ( champ(i,j,l) - eignq(i,j,l) ) * sdd2(i) 281 383 16 CONTINUE 384 ENDDO 385 ENDDO 282 386 ENDIF 283 387 c 388 DO l = 1, nbniv 389 DO j = jdfil(hemisph),jffil(hemisph) 284 390 champ( iip1,j,l ) = champ( 1,j,l ) 285 c 286 30 CONTINUE 287 c 288 50 CONTINUE 391 ENDDO 392 ENDDO 393 c 394 c 30 CONTINUE 395 c 396 c 50 CONTINUE 289 397 c 290 398 100 CONTINUE -
LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F
r418 r495 48 48 REAL xx(klon), aux(klon), coeff(klon), block(klon) 49 49 REAL dist(klon), fprime(klon), det(klon) 50 REAL pi, u(klon), v(klon), erf u(klon), erfv(klon)50 REAL pi, u(klon), v(klon), erfcu(klon), erfcv(klon) 51 51 REAL xx1(klon), xx2(klon) 52 52 real erf,kkk 53 real sqrtpi,sqrt2,zx1,zx2,exdel 53 54 c lconv = true si le calcul a converge (entre autre si qsub < min_q) 54 55 LOGICAL lconv(klon) … … 56 57 57 58 pi = ACOS(-1.) 59 sqrtpi=sqrt(pi) 60 sqrt2=sqrt(2.) 58 61 59 62 ptconv=.false. … … 123 126 xx(i) = -0.0001 124 127 else 125 xx1(i)=-SQRT(2.)*vmax(i)*(1.0-SQRT(1.0+delta(i)/(vmax(i)**2.))) 126 xx2(i)=-SQRT(2.)*vmax(i)*(1.0+SQRT(1.0+delta(i)/(vmax(i)**2.))) 128 zx1=-sqrt2*vmax(i) 129 zx2=SQRT(1.0+delta(i)/(vmax(i)**2.)) 130 xx1(i)=zx1*(1.0-zx2) 131 xx2(i)=zx1*(1.0+zx2) 127 132 xx(i) = 1.01 * xx1(i) 128 133 if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i) … … 143 148 if (.not.lconv(i)) then 144 149 145 u(i) = delta(i)/(xx(i)*sqrt (2.)) + xx(i)/(2.*sqrt(2.))146 v(i) = delta(i)/(xx(i)*sqrt (2.)) - xx(i)/(2.*sqrt(2.))150 u(i) = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2) 151 v(i) = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2) 147 152 148 153 IF ( v(i) .GT. vmax(i) ) THEN … … 153 158 c -- use asymptotic expression of erf for u and v large: 154 159 c ( -> analytic solution for xx ) 155 156 aux(i) = 2.0*delta(i)*(1.- beta(i)*EXP(delta(i)))157 : /(1.+ beta(i)*EXP(delta(i)))160 exdel=beta(i)*EXP(delta(i)) 161 aux(i) = 2.0*delta(i)*(1.-exdel) 162 : /(1.+exdel) 158 163 if (aux(i).lt.0.) then 159 print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i)164 c print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i) 160 165 aux(i)=0. 161 166 endif 162 167 xx(i) = -SQRT(aux(i)) 163 block(i) = EXP(-v(i)*v(i)) / v(i) / SQRT(pi)168 block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi 164 169 dist(i) = 0.0 165 170 fprime(i) = 1.0 … … 169 174 c -- erfv -> 1.0, use an asymptotic expression of erfv for v large: 170 175 171 erf u(i) =ERF(u(i))176 erfcu(i) = 1.0-ERF(u(i)) 172 177 c !!! ATTENTION : rajout d'un seuil pour l'exponentiel 173 aux(i) = SQRT(pi)*(1.0-erfu(i))*EXP(min(v(i)*v(i),100.))178 aux(i) = sqrtpi*erfcu(i)*EXP(min(v(i)*v(i),100.)) 174 179 coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.) 175 block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / SQRT(pi)180 block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi 176 181 dist(i) = v(i) * aux(i) / coeff(i) - beta(i) 177 182 fprime(i) = 2.0 / xx(i) * (v(i)**2.) … … 185 190 c -- general case: 186 191 187 erf u(i) =ERF(u(i))188 erf v(i) =ERF(v(i))189 block(i) = 1.0-erfv(i)190 dist(i) = (1.0 - erfu(i)) / (1.0 - erfv(i)) - beta(i)192 erfcu(i) = 1.0-ERF(u(i)) 193 erfcv(i) = 1.0-ERF(v(i)) 194 block(i) = erfcv(i) 195 dist(i) = erfcu(i) / erfcv(i) - beta(i) 191 196 zu2(i)=u(i)*u(i) 192 197 zv2(i)=v(i)*v(i) 193 198 if(zu2(i).gt.20..or. zv2(i).gt.20.) then 194 print*,'ATTENTION !!! xx(',i,') =', xx(i)195 print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',196 .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),197 .CLDF(i,k)198 print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i)199 c print*,'ATTENTION !!! xx(',i,') =', xx(i) 200 c print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF', 201 c .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k), 202 c .CLDF(i,k) 203 c print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i) 199 204 zu2(i)=20. 200 205 zv2(i)=20. 201 206 fprime(i) = 0. 202 207 else 203 fprime(i) = 2. / SQRT(pi) /xx(i) /(1.0-erfv(i))**2.204 : * ( (1.0-erfv(i))*v(i)*EXP(-zu2(i))205 : - (1.0-erfu(i))*u(i)*EXP(-zv2(i)) )208 fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2. 209 : * ( erfcv(i)*v(i)*EXP(-zu2(i)) 210 : - erfcu(i)*u(i)*EXP(-zv2(i)) ) 206 211 endif 207 212 ENDIF ! x … … 219 224 ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.) 220 225 CLDF(i,K) = 0.5 * block(i) 221 cccccccccccccccccccccccc222 c kkk=-sqrt(log(1.+ratqsc(i,k)**2))223 c u(i)=delta(i)/(kkk*sqrt(2.))-kkk/(2.*sqrt(2.))224 c v(i)=delta(i)/(kkk*sqrt(2.))+kkk/(2.*sqrt(2.))225 c erfu(i)=erf(u(i))226 c erfv(i)=erf(v(i))227 c print*,'SIG ',k,qsub(i,k)228 c s ,mu(i)*((1.-erfv(i))/(1.-erfu(i))-qsat(i)/mu(i))229 c s ,0.5*erfu(i)230 cccccccccccccccccccccccc231 226 else 232 227 xx(i) = xx(i) - dist(i)/fprime(i) -
LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F
r486 r495 804 804 110 continue 805 805 806 do 121 j=1,ntra 807 ccccc do 111 k=1,nl+1 808 do 111 k=1,nd 809 nn=0 810 do 101 i=1,len 811 if(iflag1(i).eq.0)then 812 nn=nn+1 813 tra(nn,k,j)=tra1(i,k,j) 814 endif 815 101 continue 816 111 continue 817 121 continue 806 c do 121 j=1,ntra 807 c do 111 k=1,nd 808 c nn=0 809 c do 101 i=1,len 810 c if(iflag1(i).eq.0)then 811 c nn=nn+1 812 c tra(nn,k,j)=tra1(i,k,j) 813 c endif 814 c 101 continue 815 c 111 continue 816 c 121 continue 818 817 819 818 if (nn.ne.ncum) then … … 1493 1492 400 continue 1494 1493 1495 do k=1,ntra1496 do j=1,nd ! instead nlp1497 do i=1,nd ! instead nlp1498 do il=1,ncum1499 traent(il,i,j,k)=tra(il,j,k)1500 enddo1501 enddo1502 enddo1503 enddo1494 c do k=1,ntra 1495 c do j=1,nd ! instead nlp 1496 c do i=1,nd ! instead nlp 1497 c do il=1,ncum 1498 c traent(il,i,j,k)=tra(il,j,k) 1499 c enddo 1500 c enddo 1501 c enddo 1502 c enddo 1504 1503 zm(:,:)=0. 1505 1504 … … 1557 1556 710 continue 1558 1557 1559 do k=1,ntra1560 do j=minorig,nl1561 do il=1,ncum1562 if( (i.ge.icb(il)).and.(i.le.inb(il)).and.1563 : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then1564 traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)1565 : +(1.-sij(il,i,j))*tra(il,nk(il),k)1566 endif1567 enddo1568 enddo1569 enddo1558 c do k=1,ntra 1559 c do j=minorig,nl 1560 c do il=1,ncum 1561 c if( (i.ge.icb(il)).and.(i.le.inb(il)).and. 1562 c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then 1563 c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) 1564 c : +(1.-sij(il,i,j))*tra(il,nk(il),k) 1565 c endif 1566 c enddo 1567 c enddo 1568 c enddo 1570 1569 1571 1570 c … … 1590 1589 750 continue 1591 1590 1592 do j=1,ntra1593 do i=minorig+1,nl1594 do il=1,ncum1595 if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then1596 traent(il,i,i,j)=tra(il,nk(il),j)1597 endif1598 enddo1599 enddo1600 enddo1591 c do j=1,ntra 1592 c do i=minorig+1,nl 1593 c do il=1,ncum 1594 c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then 1595 c traent(il,i,i,j)=tra(il,nk(il),j) 1596 c endif 1597 c enddo 1598 c enddo 1599 c enddo 1601 1600 1602 1601 do 100 j=minorig,nl … … 1764 1763 enddo ! il 1765 1764 1766 do j=1,ntra 1767 do il=1,ncum 1768 if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) 1769 : .and. csum(il,i).lt.m(il,i) ) then 1770 traent(il,i,i,j)=tra(il,nk(il),j) 1771 endif 1772 enddo 1773 enddo 1774 1765 c do j=1,ntra 1766 c do il=1,ncum 1767 c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) 1768 c : .and. csum(il,i).lt.m(il,i) ) then 1769 c traent(il,i,i,j)=tra(il,nk(il),j) 1770 c endif 1771 c enddo 1772 c enddo 1775 1773 789 continue 1776 1774 c … … 1869 1867 enddo 1870 1868 1871 do k=1,ntra1872 do i=1,nd1873 do il=1,ncum1874 trap(il,i,k)=tra(il,i,k)1875 enddo1876 enddo1877 enddo1869 c do k=1,ntra 1870 c do i=1,nd 1871 c do il=1,ncum 1872 c trap(il,i,k)=tra(il,i,k) 1873 c enddo 1874 c enddo 1875 c enddo 1878 1876 1879 1877 c … … 2103 2101 vp(il,i)=vp(il,i)/mp(il,i) 2104 2102 2105 do j=1,ntra2106 trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)2103 c do j=1,ntra 2104 c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) 2107 2105 ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) 2108 : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))2109 trap(il,i,j)=trap(il,i,j)/mp(il,i)2110 end do2106 c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1)) 2107 c trap(il,i,j)=trap(il,i,j)/mp(il,i) 2108 c end do 2111 2109 2112 2110 else … … 2125 2123 vp(il,i)=vp(il,i+1) 2126 2124 2127 do j=1,ntra2128 trap(il,i,j)=trap(il,i+1,j)2129 end do2125 c do j=1,ntra 2126 c trap(il,i,j)=trap(il,i+1,j) 2127 c end do 2130 2128 2131 2129 endif … … 2226 2224 enddo 2227 2225 2228 do j=1,ntra2229 do i=1,nd2230 do il=1,ncum2231 ftra(il,i,j)=0.02232 enddo2233 enddo2234 enddo2226 c do j=1,ntra 2227 c do i=1,nd 2228 c do il=1,ncum 2229 c ftra(il,i,j)=0.0 2230 c enddo 2231 c enddo 2232 c enddo 2235 2233 2236 2234 do i=1,nl … … 2330 2328 enddo ! il 2331 2329 2332 do j=1,ntra2333 do il=1,ncum2334 if (cvflag_grav) then2335 ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)2336 : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))2337 : +am(il)*(tra(il,2,j)-tra(il,1,j)))2338 else2339 ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)2340 : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))2341 : +am(il)*(tra(il,2,j)-tra(il,1,j)))2342 endif2343 enddo2344 enddo2330 c do j=1,ntra 2331 c do il=1,ncum 2332 c if (cvflag_grav) then 2333 c ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il) 2334 c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) 2335 c : +am(il)*(tra(il,2,j)-tra(il,1,j))) 2336 c else 2337 c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il) 2338 c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) 2339 c : +am(il)*(tra(il,2,j)-tra(il,1,j))) 2340 c endif 2341 c enddo 2342 c enddo 2345 2343 2346 2344 do j=2,nl … … 2366 2364 enddo 2367 2365 2368 do k=1,ntra2369 do j=2,nl2370 do il=1,ncum2371 if (j.le.inb(il)) then2372 2373 if (cvflag_grav) then2374 ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)2375 : *(traent(il,j,1,k)-tra(il,1,k))2376 else2377 ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)2378 : *(traent(il,j,1,k)-tra(il,1,k))2379 endif2380 2381 endif2382 enddo2383 enddo2384 enddo2366 c do k=1,ntra 2367 c do j=2,nl 2368 c do il=1,ncum 2369 c if (j.le.inb(il)) then 2370 2371 c if (cvflag_grav) then 2372 c ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1) 2373 c : *(traent(il,j,1,k)-tra(il,1,k)) 2374 c else 2375 c ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1) 2376 c : *(traent(il,j,1,k)-tra(il,1,k)) 2377 c endif 2378 2379 c endif 2380 c enddo 2381 c enddo 2382 c enddo 2385 2383 2386 2384 c … … 2488 2486 1350 continue 2489 2487 2490 do k=1,ntra2491 do il=1,ncum2492 if (i.le.inb(il)) then2493 dpinv=1.0/(ph(il,i)-ph(il,i+1))2494 cpinv=1.0/cpn(il,i)2495 if (cvflag_grav) then2496 ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv2497 : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))2498 : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))2499 else2500 ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv2501 : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))2502 : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))2503 endif2504 endif2505 enddo2506 enddo2488 c do k=1,ntra 2489 c do il=1,ncum 2490 c if (i.le.inb(il)) then 2491 c dpinv=1.0/(ph(il,i)-ph(il,i+1)) 2492 c cpinv=1.0/cpn(il,i) 2493 c if (cvflag_grav) then 2494 c ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv 2495 c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) 2496 c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) 2497 c else 2498 c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv 2499 c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) 2500 c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) 2501 c endif 2502 c endif 2503 c enddo 2504 c enddo 2507 2505 2508 2506 do 480 k=1,i-1 … … 2538 2536 480 continue 2539 2537 2540 do j=1,ntra2541 do k=1,i-12542 do il=1,ncum2543 if (i.le.inb(il)) then2544 dpinv=1.0/(ph(il,i)-ph(il,i+1))2545 cpinv=1.0/cpn(il,i)2546 if (cvflag_grav) then2547 ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)2548 : *(traent(il,k,i,j)-tra(il,i,j))2549 else2550 ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)2551 : *(traent(il,k,i,j)-tra(il,i,j))2552 endif2553 endif2554 enddo2555 enddo2556 enddo2538 c do j=1,ntra 2539 c do k=1,i-1 2540 c do il=1,ncum 2541 c if (i.le.inb(il)) then 2542 c dpinv=1.0/(ph(il,i)-ph(il,i+1)) 2543 c cpinv=1.0/cpn(il,i) 2544 c if (cvflag_grav) then 2545 c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) 2546 c : *(traent(il,k,i,j)-tra(il,i,j)) 2547 c else 2548 c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) 2549 c : *(traent(il,k,i,j)-tra(il,i,j)) 2550 c endif 2551 c endif 2552 c enddo 2553 c enddo 2554 c enddo 2557 2555 2558 2556 do 490 k=i,nl+1 … … 2581 2579 490 continue 2582 2580 2583 do j=1,ntra2584 do k=i,nl+12585 do il=1,ncum2586 if (i.le.inb(il) .and. k.le.inb(il)) then2587 dpinv=1.0/(ph(il,i)-ph(il,i+1))2588 cpinv=1.0/cpn(il,i)2589 if (cvflag_grav) then2590 ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)2591 : *(traent(il,k,i,j)-tra(il,i,j))2592 else2593 ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)2594 : *(traent(il,k,i,j)-tra(il,i,j))2595 endif2596 endif ! i and k2597 enddo2598 enddo2599 enddo2581 c do j=1,ntra 2582 c do k=i,nl+1 2583 c do il=1,ncum 2584 c if (i.le.inb(il) .and. k.le.inb(il)) then 2585 c dpinv=1.0/(ph(il,i)-ph(il,i+1)) 2586 c cpinv=1.0/cpn(il,i) 2587 c if (cvflag_grav) then 2588 c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) 2589 c : *(traent(il,k,i,j)-tra(il,i,j)) 2590 c else 2591 c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) 2592 c : *(traent(il,k,i,j)-tra(il,i,j)) 2593 c endif 2594 c endif ! i and k 2595 c enddo 2596 c enddo 2597 c enddo 2600 2598 2601 2599 do 1400 il=1,ncum … … 2654 2652 enddo 2655 2653 2656 do j=1,ntra 2657 do il=1,ncum 2658 if (i.le.inb(il)) then 2659 dpinv=1.0/(ph(il,i)-ph(il,i+1)) 2660 cpinv=1.0/cpn(il,i) 2661 2662 if (cvflag_grav) then 2663 ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv 2664 : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) 2665 : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) 2666 else 2667 ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv 2668 : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) 2669 : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) 2670 endif 2671 endif ! i 2672 enddo 2673 enddo 2674 2654 c do j=1,ntra 2655 c do il=1,ncum 2656 c if (i.le.inb(il)) then 2657 c dpinv=1.0/(ph(il,i)-ph(il,i+1)) 2658 c cpinv=1.0/cpn(il,i) 2659 2660 c if (cvflag_grav) then 2661 c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv 2662 c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) 2663 c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) 2664 c else 2665 c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv 2666 c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) 2667 c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) 2668 c endif 2669 c endif ! i 2670 c enddo 2671 c enddo 2675 2672 2676 2673 500 continue … … 2715 2712 503 continue 2716 2713 2717 do j=1,ntra2718 do il=1,ncum2719 ex=0.1*ment(il,inb(il),inb(il))2720 : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))2721 : /(ph(il,inb(il))-ph(il,inb(il)+1))2722 ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex2723 ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)2724 : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))2725 : /(ph(il,inb(il)-1)-ph(il,inb(il)))2726 enddo2727 enddo2714 c do j=1,ntra 2715 c do il=1,ncum 2716 c ex=0.1*ment(il,inb(il),inb(il)) 2717 c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) 2718 c : /(ph(il,inb(il))-ph(il,inb(il)+1)) 2719 c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex 2720 c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) 2721 c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) 2722 c : /(ph(il,inb(il)-1)-ph(il,inb(il))) 2723 c enddo 2724 c enddo 2728 2725 2729 2726 c … … 2981 2978 end 2982 2979 2980 SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na, 2981 & ment,sij,da,phi) 2982 implicit none 2983 c inputs: 2984 integer ncum, nd, na, nloc,len 2985 real ment(nloc,na,na),sij(nloc,na,na) 2986 c ouputs: 2987 real da(nloc,na),phi(nloc,na,na) 2988 c local variables: 2989 integer i,j,k 2990 c 2991 da(:,:)=0. 2992 c 2993 do j=1,na 2994 do k=1,na 2995 do i=1,ncum 2996 da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j) 2997 phi(i,j,k)=sij(i,k,j)*ment(i,k,j) 2998 c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j) 2999 end do 3000 end do 3001 end do 3002 3003 return 3004 end 3005 2983 3006 2984 3007 SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum … … 3051 3074 3052 3075 3053 do 2100 j=1,ntra3054 do 2110 k=1,nd ! oct33055 do 2120 i=1,ncum3056 ftra1(idcum(i),k,j)=ftra(i,k,j)3057 2120 continue3058 2110 continue3059 2100 continue3076 c do 2100 j=1,ntra 3077 c do 2110 k=1,nd ! oct3 3078 c do 2120 i=1,ncum 3079 c ftra1(idcum(i),k,j)=ftra(i,k,j) 3080 c 2120 continue 3081 c 2110 continue 3082 c 2100 continue 3060 3083 3061 3084 return -
LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F
r493 r495 315 315 call gwprofil 316 316 * ( nlon , nlev 317 * , kgwd , kdx 317 * , kgwd , kdx , ktest 318 318 * , ikcrith, icrit 319 319 * , paphm1, zrho , zstab , zvph … … 343 343 c 344 344 c 345 do 523 jl=1,kgwd 346 ji=kdx(jl) 345 c do 523 jl=1,kgwd 346 c ji=kdx(jl) 347 c Modif vectorisation 02/04/2004 348 do 523 ji=kidia,kfdia 349 if(ktest(ji).eq.1) then 350 347 351 zdelp=paphm1(ji,jk+1)-paphm1(ji,jk) 348 352 ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp) … … 401 405 pte(ji,jk)=0.0 402 406 407 endif 403 408 523 continue 404 409 … … 1007 1012 SUBROUTINE GWPROFIL 1008 1013 * ( NLON, NLEV 1009 * , kgwd, kdx 1014 * , kgwd, kdx , ktest 1010 1015 * , KKCRITH, KCRIT 1011 1016 * , PAPHM1, PRHO , PSTAB , PVPH , PRI , PTAU … … 1075 1080 integer nlon,nlev 1076 1081 INTEGER KKCRITH(NLON),KCRIT(NLON) 1077 * ,kdx(nlon) 1082 * ,kdx(nlon) , ktest(nlon) 1083 1078 1084 C 1079 1085 REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1), … … 1109 1115 ilevh=KLEV/3 1110 1116 C 1111 DO 400 ji=1,kgwd 1112 jl=kdx(ji) 1117 c DO 400 ji=1,kgwd 1118 c jl=kdx(ji) 1119 c Modif vectorisation 02/04/2004 1120 DO 400 jl=kidia,kfdia 1121 if (ktest(jl).eq.1) then 1113 1122 Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0) 1114 1123 ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1) 1124 endif 1115 1125 400 CONTINUE 1116 1126 … … 1123 1133 410 CONTINUE 1124 1134 C 1125 DO 411 ji=1,kgwd 1126 jl=kdx(ji) 1135 c DO 411 ji=1,kgwd 1136 c jl=kdx(ji) 1137 c Modif vectorisation 02/04/2004 1138 do 411 jl=kidia,kfdia 1139 if (ktest(jl).eq.1) then 1127 1140 IF(JK.GT.KKCRITH(JL)) THEN 1128 1141 PTAU(JL,JK)=ZTAU(JL,KLEV+1) … … 1132 1145 PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1) 1133 1146 ENDIF 1147 endif 1134 1148 411 CONTINUE 1135 1149 C … … 1143 1157 420 CONTINUE 1144 1158 C 1145 DO 421 ji=1,kgwd 1146 jl=kdx(ji) 1159 c DO 421 ji=1,kgwd 1160 c jl=kdx(ji) 1161 c Modif vectorisation 02/04/2004 1162 do 421 jl=kidia,kfdia 1163 if(ktest(jl).eq.1) then 1147 1164 IF(JK.LT.KKCRITH(JL)) THEN 1148 1165 ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK) … … 1150 1167 ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec) 1151 1168 ENDIF 1169 endif 1152 1170 421 CONTINUE 1153 1171 C … … 1157 1175 C 1158 1176 1159 DO 431 ji=1,kgwd 1160 jl=kdx(ji) 1177 c DO 431 ji=1,kgwd 1178 c jl=Kdx(ji) 1179 c Modif vectorisation 02/04/2004 1180 do 431 jl=kidia,kfdia 1181 if(ktest(jl).eq.1) then 1182 1161 1183 IF(JK.LT.KKCRITH(JL)) THEN 1162 1184 IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN … … 1178 1200 ENDIF 1179 1201 ENDIF 1202 endif 1180 1203 431 CONTINUE 1181 1204 … … 1185 1208 C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL 1186 1209 1187 DO 530 ji=1,kgwd 1188 jl=kdx(ji) 1210 c DO 530 ji=1,kgwd 1211 c jl=kdx(ji) 1212 c Modif vectorisation 02/04/2004 1213 do 530 jl=kidia,kfdia 1214 if(ktest(jl).eq.1) then 1189 1215 ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL)) 1190 1216 ZTAU(JL,NSTRA)=PTAU(JL,NSTRA) 1217 endif 1191 1218 530 CONTINUE 1192 1219 1193 1220 DO 531 JK=1,KLEV 1194 1221 1195 DO 532 ji=1,kgwd 1196 jl=kdx(ji) 1222 c DO 532 ji=1,kgwd 1223 c jl=kdx(ji) 1224 c Modif vectorisation 02/04/2004 1225 do 532 jl=kidia,kfdia 1226 if(ktest(jl).eq.1) then 1227 1197 1228 1198 1229 IF(JK.GT.KKCRITH(JL))THEN … … 1206 1237 ENDIF 1207 1238 1239 endif 1208 1240 532 CONTINUE 1209 1241 1210 1242 C REORGANISATION IN THE STRATOSPHERE 1211 1243 1212 DO 533 ji=1,kgwd 1213 jl=kdx(ji) 1244 c DO 533 ji=1,kgwd 1245 c jl=kdx(ji) 1246 c Modif vectorisation 02/04/2004 1247 do 533 jl=kidia,kfdia 1248 if(ktest(jl).eq.1) then 1249 1214 1250 1215 1251 IF(JK.LT.NSTRA)THEN … … 1221 1257 ENDIF 1222 1258 1259 endif 1223 1260 533 CONTINUE 1224 1261 1225 1262 C REORGANISATION IN THE TROPOSPHERE 1226 1263 1227 DO 534 ji=1,kgwd 1228 jl=kdx(ji) 1264 c DO 534 ji=1,kgwd 1265 c jl=kdx(ji) 1266 c Modif vectorisation 02/04/2004 1267 do 534 jl=kidia,kfdia 1268 if(ktest(jl).eq.1) then 1269 1229 1270 1230 1271 IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN … … 1236 1277 1237 1278 ENDIF 1279 endif 1238 1280 534 CONTINUE 1239 1281
Note: See TracChangeset
for help on using the changeset viewer.