Changeset 3711
- Timestamp:
- Jun 11, 2020, 11:09:41 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90
r3710 r3711 131 131 132 132 INTEGER, SAVE :: igout = 1 133 134 135 133 !$omp THREADPRIVATE(igout) 134 135 ! -- Mixing probability distribution functions 136 136 137 137 REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F … … 149 149 150 150 INTEGER, SAVE :: ifrst = 0 151 152 153 154 155 156 151 !$omp THREADPRIVATE(ifrst) 152 153 ! ===================================================================== 154 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS 155 ! ===================================================================== 156 ! -- Initialize mixing PDF coefficients 157 157 IF (ifrst == 0) THEN 158 158 ifrst = 1 … … 181 181 Sigij(1:ncum, 1:nd, 1:nd) = 0.0 182 182 183 184 185 186 187 183 ! ===================================================================== 184 ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING 185 ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING 186 ! --- FRACTION (Sij) 187 ! ===================================================================== 188 188 189 189 DO i = minorig + 1, nl … … 253 253 ENDIF 254 254 255 256 255 ! *** if no air can entrain at level i assume that updraft detrains *** 256 ! *** at that level and calculate detrained air flux and properties *** 257 257 258 258 DO il = 1, ncum … … 282 282 END DO 283 283 284 ! ===================================================================== 285 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES 286 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING 287 ! ===================================================================== 288 289 CALL zilch(csum, nloc*nd) 284 ! ===================================================================== 285 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES 286 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING 287 ! ===================================================================== 290 288 291 289 DO il = 1, ncum … … 293 291 END DO 294 292 295 293 ! --------------------------------------------------------------- 296 294 DO i = minorig + 1, nl !Loop on origin level "i" 297 295 ! --------------------------------------------------------------- 298 296 299 297 DO il = 1, ncum … … 341 339 Sbef(il) = max(0., signhpmh(il)) 342 340 endif 343 END IF 344 END DO 345 346 DO il = 1, ncum 347 DO j = minorig, nl !Loop on destination level "j" 348 IF (i >= icb(il) .AND. i <= inb(il) .AND. & 349 j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & 350 lwork(il) .AND. Sij(il, i, j) > 0.0) THEN 351 ! ----------------------------------------------- 352 IF (j > i) THEN 353 Smid = min(Sij(il, i, j), Scrit(il)) 354 Sjmax(il) = Smid 355 Sjmin(il) = Smid 356 IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN 357 smin(il) = min(smin(il), Smid) 358 Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il)) 359 Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 360 Sjmin(il) = min(Sjmin(il), Scrit(il)) 361 Sbef(il) = Sij(il, i, j) 341 342 DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" 343 IF (lwork(il) .AND. Sij(il, i, j) > 0.0) THEN 344 ! ----------------------------------------------- 345 IF (j > i) THEN 346 Smid = min(Sij(il, i, j), Scrit(il)) 347 Sjmax(il) = Smid 348 Sjmin(il) = Smid 349 IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN 350 smin(il) = min(smin(il), Smid) 351 Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il)) 352 Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 353 Sjmin(il) = min(Sjmin(il), Scrit(il)) 354 Sbef(il) = Sij(il, i, j) 355 END IF 356 ELSE IF (j == i) THEN 357 Smid = 1. 358 Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + & 359 min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il)) 360 Sjmin(il) = max(Sjmin(il), sup(il)) 361 Sjmax(il) = 1. 362 ! preparation des variables Scrit, Smin et Sbef pour la partie j>i 363 Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il)) 364 365 smin(il) = 1. 366 Sbef(il) = max(0., signhpmh(il)) 367 supmax(il, i) = sign(Scrit(il), -signhpmh(il)) 368 ELSE IF (j < i) THEN 369 Smid = max(Sij(il, i, j), Scrit(il)) 370 Sjmax(il) = Smid 371 Sjmin(il) = Smid 372 IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN 373 smax(il) = Smid 374 Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) 375 Sjmax(il) = max(Sjmax(il), Scrit(il)) 376 Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 377 Sjmin(il) = max(Sjmin(il), Scrit(il)) 378 Sbef(il) = Sij(il, i, j) 379 END IF 380 IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) & 381 sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) 362 382 END IF 363 ELSE IF (j == i) THEN 364 Smid = 1. 365 Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + & 366 min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il)) 367 Sjmin(il) = max(Sjmin(il), sup(il)) 368 Sjmax(il) = 1. 369 ! preparation des variables Scrit, Smin et Sbef pour la partie j>i 370 Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il)) 371 372 smin(il) = 1. 373 Sbef(il) = max(0., signhpmh(il)) 374 supmax(il, i) = sign(Scrit(il), -signhpmh(il)) 375 ELSE IF (j < i) THEN 376 Smid = max(Sij(il, i, j), Scrit(il)) 377 Sjmax(il) = Smid 378 Sjmin(il) = Smid 379 IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN 380 smax(il) = Smid 381 Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) 382 Sjmax(il) = max(Sjmax(il), Scrit(il)) 383 Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 384 Sjmin(il) = max(Sjmin(il), Scrit(il)) 385 Sbef(il) = Sij(il, i, j) 383 ! ----------------------------------------------- 384 385 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 386 Qmixmax(il) = Qmix(Sjmax(il)) 387 Qmixmin(il) = Qmix(Sjmin(il)) 388 Rmixmax(il) = Rmix(Sjmax(il)) 389 Rmixmin(il) = Rmix(Sjmin(il)) 390 sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il)) 391 sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il)) 392 393 Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j) 394 IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN 395 Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il)) 396 ELSE 397 Sigij(il, i, j) = 0. 386 398 END IF 387 IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) & 388 sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) 399 400 ! Compute Qent, uent, vent according to the true mixing fraction 401 Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) 402 uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) 403 vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) 404 405 ! Compute liquid water static energy of mixed draughts 406 hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) 407 ! Heat capacity of mixed draught 408 cpm = cpd + Qent(il, i, j)*(cpv - cpd) 409 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 410 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 411 elij(il, i, j) = elij(il, i, j) + & 412 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 413 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 414 elij(il, i, j) = elij(il, i, j)/ & 415 (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & 416 (cpm*rrv*t(il, j)*t(il, j))) 417 ELSE 418 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 419 elij(il, i, j) = elij(il, i, j) + & 420 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 421 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 422 elij(il, i, j) = elij(il, i, j)/ & 423 (1.+lv(il, j)*lv(il, j)*rs(il, j)/ & 424 (cpm*rrv*t(il, j)*t(il, j))) 425 ENDIF 426 elij(il, i, j) = max(elij(il, i, j), 0.) 427 428 elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) 429 430 IF (j > i) THEN 431 awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) 432 awat = amax1(awat, 0.0) 433 ELSE 434 awat = 0. 435 END IF 436 437 ! Mixed draught temperature at level j 438 Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j)) 439 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 440 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat 441 ELSE 442 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat 443 ENDIF 444 445 ! ASij is the integral of P(F) over the relevant F interval 446 ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 447 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 448 449 ! If I=J (detrainement and entrainement at the same level), then only the 450 ! adiabatic ascent part of the mixture is considered 451 IF (i == j) THEN 452 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 453 Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 454 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 455 Qent(il, i, i) = rti 456 uent(il, i, i) = unk(il) 457 vent(il, i, i) = vnk(il) 458 hent(il, i, i) = hp(il, i) 459 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 460 Sigij(il, i, i) = 0. 461 END IF 389 462 END IF 390 ! ----------------------------------------------- 391 392 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 393 Qmixmax(il) = Qmix(Sjmax(il)) 394 Qmixmin(il) = Qmix(Sjmin(il)) 395 Rmixmax(il) = Rmix(Sjmax(il)) 396 Rmixmin(il) = Rmix(Sjmin(il)) 397 sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il)) 398 sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il)) 399 400 Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j) 401 IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN 402 Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il)) 403 ELSE 404 Sigij(il, i, j) = 0. 405 END IF 406 407 ! Compute Qent, uent, vent according to the true mixing fraction 408 Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) 409 uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) 410 vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) 411 412 ! Compute liquid water static energy of mixed draughts 413 hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) 414 ! Heat capacity of mixed draught 415 cpm = cpd + Qent(il, i, j)*(cpv - cpd) 416 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 417 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 418 elij(il, i, j) = elij(il, i, j) + & 419 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 420 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 421 elij(il, i, j) = elij(il, i, j)/ & 422 (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & 423 (cpm*rrv*t(il, j)*t(il, j))) 424 ELSE 425 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 426 elij(il, i, j) = elij(il, i, j) + & 427 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 428 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 429 elij(il, i, j) = elij(il, i, j)/ & 430 (1.+lv(il, j)*lv(il, j)*rs(il, j)/ & 431 (cpm*rrv*t(il, j)*t(il, j))) 432 ENDIF 433 elij(il, i, j) = max(elij(il, i, j), 0.) 434 435 elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) 436 437 IF (j > i) THEN 438 awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) 439 awat = amax1(awat, 0.0) 440 ELSE 441 awat = 0. 442 END IF 443 444 ! Mixed draught temperature at level j 445 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 446 Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j)) 447 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat 448 ELSE 449 Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j)) 450 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat 451 ENDIF 452 453 ! ASij is the integral of P(F) over the relevant F interval 454 ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 455 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 456 457 ! If I=J (detrainement and entrainement at the same level), then only the 458 ! adiabatic ascent part of the mixture is considered 459 IF (i == j) THEN 460 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 461 Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 462 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 463 Qent(il, i, i) = rti 463 END DO ! Loop j = (icb(il) - 1), inb(il) 464 465 IF (lwork(il)) THEN 466 csum(il, i) = 0 467 DO j = (icb(il) - 1), inb(il) 468 ASij(il) = amax1(1.0E-16, ASij(il)) 469 ASij_inv(il) = 1.0/ASij(il) 470 ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs 471 IF (ASij_inv(il) > 100.) ASij_inv(il) = 0. 472 Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) 473 csum(il, i) = csum(il, i) + Ment(il, i, j) 474 END DO 475 IF (csum(il, i) < 1.) THEN 476 nent(il, i) = 0 477 Ment(il, i, i) = 1. 478 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 464 479 uent(il, i, i) = unk(il) 465 480 vent(il, i, i) = vnk(il) 466 hent(il, i, i) = hp(il, i)467 481 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 468 Sigij(il, i, i) = 0. 469 END IF 470 END IF 471 END DO ! End loop on destination level "j" 472 ! --------------------------------------------------------------- 473 END DO 474 475 DO il = 1, ncum 476 IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN 477 ASij(il) = amax1(1.0E-16, ASij(il)) 478 ASij_inv(il) = 1.0/ASij(il) 479 ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs 480 IF (ASij_inv(il) > 100.) ASij_inv(il) = 0. 481 csum(il, i) = 0.0 482 END IF 483 END DO 484 485 DO j = minorig, nl 486 DO il = 1, ncum 487 IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. & 488 j >= (icb(il) - 1) .AND. j <= inb(il)) THEN 489 Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) 490 END IF 491 END DO 492 END DO 493 494 DO j = minorig, nl 495 DO il = 1, ncum 496 IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. & 497 j >= (icb(il) - 1) .AND. j <= inb(il)) THEN 498 csum(il, i) = csum(il, i) + Ment(il, i, j) 499 END IF 500 END DO 501 END DO 502 503 DO il = 1, ncum 504 IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. csum(il, i) < 1.) THEN 505 nent(il, i) = 0 506 Ment(il, i, i) = 1. 507 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 508 uent(il, i, i) = unk(il) 509 vent(il, i, i) = vnk(il) 510 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 511 IF (fl_cor_ebil .GE. 2) THEN 512 hent(il, i, i) = hp(il, i) 513 Sigij(il, i, i) = 0.0 514 ELSE 515 Sij(il, i, i) = 0.0 516 ENDIF 517 END IF 518 END DO ! il 519 520 ! --------------------------------------------------------------- 521 END DO ! End loop on origin level "i" 522 ! --------------------------------------------------------------- 523 RETURN 482 IF (fl_cor_ebil .GE. 2) THEN 483 hent(il, i, i) = hp(il, i) 484 Sigij(il, i, i) = 0.0 485 ELSE 486 Sij(il, i, i) = 0.0 487 ENDIF 488 ENDIF 489 END IF 490 END IF ! i >= icb(il) .AND. i <= inb(il) 491 END DO ! Loop il = 1, ncum 492 END DO ! End loop on origin level "i" 524 493 END SUBROUTINE 525 494
Note: See TracChangeset
for help on using the changeset viewer.
