Changeset 3715 for LMDZ6/branches/Optimisation_LMDZ/libf/phylmd
- Timestamp:
- Jun 11, 2020, 11:09:43 AM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90
r3714 r3715 112 112 REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp 113 113 REAL :: alt, delp, delm 114 REAL, DIMENSION(nloc) :: Qmixmax, Rmixmax, sqmrmax 115 REAL, DIMENSION(nloc) :: Qmixmin, Rmixmin, sqmrmin 116 REAL, DIMENSION(nloc) :: signhpmh 117 REAL, DIMENSION(nloc) :: Sx, smin 118 REAL :: Scrit2 119 REAL, DIMENSION(nloc) :: Sjmin, Sjmax 120 REAL, DIMENSION(nloc) :: Sbef, sup 121 REAL, DIMENSION(nloc, nd) :: ASij 122 REAL, DIMENSION(nloc) :: ASij_inv, smax, Scrit 123 REAL, DIMENSION(nloc, nd, nd) :: Sij 124 REAL, DIMENSION(nloc, nd) :: csum 114 REAL, DIMENSION(nloc, nd) :: ASij 115 REAL, DIMENSION(nloc, nd, nd) :: Sij 125 116 REAL :: awat 126 117 REAL :: cpm !Mixed draught heat capacity … … 179 170 Ment(1:ncum, 1:nd, 1:nd) = 0.0 180 171 Sij(1:ncum, 1:nd, 1:nd) = 0.0 181 Sigij(1:ncum, 1:nd, 1:nd) = 0.0182 172 183 173 ! ===================================================================== … … 287 277 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING 288 278 ! ===================================================================== 289 290 279 DO il = 1, ncum 291 280 DO i = icb(il), inb(il) !Loop on origin level "i" 292 signhpmh(il) = sign(1., hp(il, i) - h(il, i)) 293 294 Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il)) 295 296 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 297 IF (cvflag_ice) THEN 298 anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & 299 (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 300 denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & 301 (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 302 ELSE 303 anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & 304 (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 305 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & 306 (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 307 END IF 308 IF (abs(denom) < 0.01) denom = 0.01 309 Scrit(il) = min(anum/denom, 1.) 310 alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti) 311 312 ! Find new critical value Scrit2 313 ! such that : Sij > Scrit2 => mixed draught will detrain at J<I 314 ! Sij < Scrit2 => mixed draught will detrain at J>I 315 Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + & 316 Scrit(il)*max(0., signhpmh(il)) 317 Scrit(il) = Scrit2 318 319 ! Correction pour la nouvelle logique; la correction pour ALT 320 ! est un peu au hazard 321 IF (Scrit(il) <= 0.0) Scrit(il) = 0.0 322 IF (alt <= 0.0) Scrit(il) = 1.0 323 324 smax(il) = 0.0 325 ASij(il, i) = 0.0 326 sup(il) = 0. ! upper S-value reached by descending draughts 327 328 if (i < inb(il)) then 329 Sbef(il) = Sij(il, i, inb(il)) 330 else 331 Sbef(il) = max(0., signhpmh(il)) 332 endif 333 334 DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" 335 IF (Sij(il, i, j) > 0.0) THEN 336 ! ----------------------------------------------- 337 IF (j > i) THEN 338 Smid = min(Sij(il, i, j), Scrit(il)) 339 Sjmax(il) = Smid 340 Sjmin(il) = Smid 341 IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN 342 smin(il) = min(smin(il), Smid) 343 Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il)) 344 Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 345 Sjmin(il) = min(Sjmin(il), Scrit(il)) 346 Sbef(il) = Sij(il, i, j) 347 END IF 348 ELSE IF (j == i) THEN 349 Smid = 1. 350 Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + & 351 min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il)) 352 Sjmin(il) = max(Sjmin(il), sup(il)) 353 Sjmax(il) = 1. 354 ! preparation des variables Scrit, Smin et Sbef pour la partie j>i 355 Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il)) 356 357 smin(il) = 1. 358 Sbef(il) = max(0., signhpmh(il)) 359 supmax(il, i) = sign(Scrit(il), -signhpmh(il)) 360 ELSE IF (j < i) THEN 361 Smid = max(Sij(il, i, j), Scrit(il)) 362 Sjmax(il) = Smid 363 Sjmin(il) = Smid 364 IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN 365 smax(il) = Smid 366 Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) 367 Sjmax(il) = max(Sjmax(il), Scrit(il)) 368 Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 369 Sjmin(il) = max(Sjmin(il), Scrit(il)) 370 Sbef(il) = Sij(il, i, j) 371 END IF 372 IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) & 373 sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) 281 block 282 real :: signhpmh, Sx, Scrit 283 real :: smax, sup, Sbef, smin 284 285 signhpmh = sign(1., hp(il, i) - h(il, i)) 286 287 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 288 IF (cvflag_ice) THEN 289 anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & 290 (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 291 denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & 292 (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 293 ELSE 294 anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & 295 (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 296 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & 297 (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 298 END IF 299 IF (abs(denom) < 0.01) denom = 0.01 300 Scrit = min(anum/denom, 1.) 301 alt = rti - rs(il, i) + Scrit*(rr(il, i) - rti) 302 303 ! Get max of Sij 304 Sx = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh) 305 ! Find new critical value Scrit 306 ! such that : Sij > Scrit => mixed draught will detrain at J<I 307 ! Sij < Scrit => mixed draught will detrain at J>I 308 Scrit = min(Scrit, Sx*max(0., -signhpmh)) + Scrit*max(0., signhpmh) 309 310 ! Correction pour la nouvelle logique; la correction pour ALT 311 ! est un peu au hazard 312 IF (Scrit <= 0.0) Scrit = 0.0 313 IF (alt <= 0.0) Scrit = 1.0 314 315 smax = 0.0 316 ASij(il, i) = 0.0 317 sup = 0. ! upper S-value reached by descending draughts 318 319 ! Glitchy : why? 320 if (i < inb(il)) then 321 Sbef = Sij(il, i, inb(il)) 322 else 323 Sbef = max(0., signhpmh) 324 endif 325 326 DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" 327 IF (Sij(il, i, j) > 0.0) THEN 328 block 329 real :: Sjmax, Sjmin, Qmixmax, Qmixmin, Rmixmax, Rmixmin, sqmrmax, sqmrmin 330 ! ----------------------------------------------- 331 IF (j > i) THEN 332 Smid = min(Sij(il, i, j), Scrit) 333 Sjmax = Smid 334 Sjmin = Smid 335 IF (Smid < smin .AND. Sij(il, i, j + 1) < Smid) THEN 336 smin = min(smin, Smid) 337 Sjmax = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit) 338 Sjmin = max((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) 339 Sjmin = min(Sjmin, Scrit) 340 Sbef = Sij(il, i, j) 341 END IF 342 ELSE IF (j == i) THEN 343 Smid = 1. 344 Sjmin = max((Sij(il, i, j - 1) + Smid)/2., Scrit)*max(0., -signhpmh) + & 345 min((Sij(il, i, j + 1) + Smid)/2., Scrit)*max(0., signhpmh) 346 Sjmin = max(Sjmin, sup) 347 Sjmax = 1. 348 ! preparation des variables Scrit, Smin et Sbef pour la partie j>i 349 Scrit = min(Sjmin, Sjmax, Scrit) 350 351 smin = 1. 352 Sbef = max(0., signhpmh) 353 supmax(il, i) = sign(Scrit, -signhpmh) 354 ELSE IF (j < i) THEN 355 Smid = max(Sij(il, i, j), Scrit) 356 Sjmax = Smid 357 Sjmin = Smid 358 IF (Smid > smax .AND. Sij(il, i, j + 1) > Smid) THEN 359 smax = Smid 360 Sjmax = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) 361 Sjmax = max(Sjmax, Scrit) 362 Sjmin = min((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) 363 Sjmin = max(Sjmin, Scrit) 364 Sbef = Sij(il, i, j) 365 END IF 366 IF (abs(Sjmin - Sjmax) > 1.E-10) & 367 sup = max(Sjmin, Sjmax, sup) 368 END IF 369 ! ----------------------------------------------- 370 371 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 372 Qmixmax = Qmix(Sjmax) 373 Qmixmin = Qmix(Sjmin) 374 Rmixmax = Rmix(Sjmax) 375 Rmixmin = Rmix(Sjmin) 376 sqmrmax = Sjmax*Qmix(Sjmax) - Rmix(Sjmax) 377 sqmrmin = Sjmin*Qmix(Sjmin) - Rmix(Sjmin) 378 379 Ment(il, i, j) = abs(Qmixmax - Qmixmin)*Ment(il, i, j) 380 IF (abs(Qmixmax - Qmixmin) > 1.E-10) THEN 381 Sigij(il, i, j) = (sqmrmax - sqmrmin)/(Qmixmax - Qmixmin) 382 ELSE 383 Sigij(il, i, j) = 0. 384 END IF 385 386 ! Compute Qent, uent, vent according to the true mixing fraction 387 Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) 388 uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) 389 vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) 390 391 ! Compute liquid water static energy of mixed draughts 392 hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) 393 ! Heat capacity of mixed draught 394 cpm = cpd + Qent(il, i, j)*(cpv - cpd) 395 396 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 397 elij(il, i, j) = elij(il, i, j) + & 398 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 399 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 400 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 401 elij(il, i, j) = elij(il, i, j)/ & 402 (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & 403 (cpm*rrv*t(il, j)*t(il, j))) 404 ELSE 405 406 elij(il, i, j) = elij(il, i, j)/ & 407 (1.+lv(il, j)*lv(il, j)*rs(il, j)/ & 408 (cpm*rrv*t(il, j)*t(il, j))) 409 ENDIF 410 elij(il, i, j) = max(elij(il, i, j), 0.) 411 elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) 412 413 IF (j > i) THEN 414 awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) 415 awat = amax1(awat, 0.0) 416 ELSE 417 awat = 0. 418 END IF 419 420 ! Mixed draught temperature at level j 421 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)) 422 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 423 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat 424 ELSE 425 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat 426 ENDIF 427 428 ! ASij is the integral of P(F) over the relevant F interval 429 ASij(il, i) = ASij(il, i) + abs(Qmixmax*(1.-Sjmax) + Rmixmax - & 430 Qmixmin*(1.-Sjmin) - Rmixmin) 431 432 ! If I=J (detrainement and entrainement at the same level), then only the 433 ! adiabatic ascent part of the mixture is considered 434 IF (i == j) THEN 435 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 436 Ment(il, i, i) = abs(Qmixmax*(1.-Sjmax) + Rmixmax - & 437 Qmixmin*(1.-Sjmin) - Rmixmin) 438 Qent(il, i, i) = rti 439 uent(il, i, i) = unk(il) 440 vent(il, i, i) = vnk(il) 441 hent(il, i, i) = hp(il, i) 442 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 443 Sigij(il, i, i) = 0. 444 END IF 445 end block 374 446 END IF 375 ! ----------------------------------------------- 376 377 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 378 Qmixmax(il) = Qmix(Sjmax(il)) 379 Qmixmin(il) = Qmix(Sjmin(il)) 380 Rmixmax(il) = Rmix(Sjmax(il)) 381 Rmixmin(il) = Rmix(Sjmin(il)) 382 sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il)) 383 sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il)) 384 385 Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j) 386 IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN 387 Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il)) 388 ELSE 389 Sigij(il, i, j) = 0. 390 END IF 391 392 ! Compute Qent, uent, vent according to the true mixing fraction 393 Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) 394 uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) 395 vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) 396 397 ! Compute liquid water static energy of mixed draughts 398 hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) 399 ! Heat capacity of mixed draught 400 cpm = cpd + Qent(il, i, j)*(cpv - cpd) 401 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 402 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 403 elij(il, i, j) = elij(il, i, j) + & 404 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 405 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 406 elij(il, i, j) = elij(il, i, j)/ & 407 (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & 408 (cpm*rrv*t(il, j)*t(il, j))) 409 ELSE 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)*lv(il, j)*rs(il, j)/ & 416 (cpm*rrv*t(il, j)*t(il, j))) 417 ENDIF 418 elij(il, i, j) = max(elij(il, i, j), 0.) 419 420 elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) 421 422 IF (j > i) THEN 423 awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) 424 awat = amax1(awat, 0.0) 425 ELSE 426 awat = 0. 427 END IF 428 429 ! Mixed draught temperature at level j 430 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)) 431 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 432 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat 433 ELSE 434 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat 435 ENDIF 436 437 ! ASij is the integral of P(F) over the relevant F interval 438 ASij(il, i) = ASij(il, i) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 439 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 440 441 ! If I=J (detrainement and entrainement at the same level), then only the 442 ! adiabatic ascent part of the mixture is considered 443 IF (i == j) THEN 444 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 445 Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 446 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 447 Qent(il, i, i) = rti 448 uent(il, i, i) = unk(il) 449 vent(il, i, i) = vnk(il) 450 hent(il, i, i) = hp(il, i) 451 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 452 Sigij(il, i, i) = 0. 453 END IF 454 END IF 455 END DO ! Loop j = (icb(il) - 1), inb(il) 447 END DO ! Loop j = (icb(il) - 1), inb(il) 448 end block 456 449 END DO ! Loop i = icb(il), inb(il) 457 450 END DO ! Loop il = 1, ncum … … 459 452 DO il = 1, ncum 460 453 DO i = icb(il), inb(il) !Loop on origin level "i" 461 csum(il, i) = 0 462 DO j = (icb(il) - 1), inb(il) 463 ASij(il, i) = amax1(1.0E-16, ASij(il, i)) 464 ASij_inv(il) = 1.0/ASij(il, i) 465 ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs 466 IF (ASij_inv(il) > 100.) ASij_inv(il) = 0. 467 Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) 468 csum(il, i) = csum(il, i) + Ment(il, i, j) 469 END DO 470 IF (csum(il, i) < 1.) THEN 471 nent(il, i) = 0 472 Ment(il, i, i) = 1. 473 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 474 uent(il, i, i) = unk(il) 475 vent(il, i, i) = vnk(il) 476 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 477 IF (fl_cor_ebil .GE. 2) THEN 478 hent(il, i, i) = hp(il, i) 479 Sigij(il, i, i) = 0.0 480 ELSE 481 Sij(il, i, i) = 0.0 454 block 455 real :: csum, ASij_inv 456 csum = 0 457 DO j = (icb(il) - 1), inb(il) 458 ASij(il, i) = amax1(1.0E-16, ASij(il, i)) 459 ASij_inv = 1.0/ASij(il, i) 460 ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs 461 IF (ASij_inv > 100.) ASij_inv = 0. 462 Ment(il, i, j) = Ment(il, i, j)*ASij_inv 463 csum = csum + Ment(il, i, j) 464 END DO 465 IF (csum < 1.) THEN 466 nent(il, i) = 0 467 Ment(il, i, i) = 1. 468 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 469 uent(il, i, i) = unk(il) 470 vent(il, i, i) = vnk(il) 471 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 472 IF (fl_cor_ebil .GE. 2) THEN 473 hent(il, i, i) = hp(il, i) 474 Sigij(il, i, i) = 0.0 475 ELSE 476 Sij(il, i, i) = 0.0 477 ENDIF 482 478 ENDIF 483 ENDIF479 end block 484 480 END DO ! Loop il = 1, ncum 485 481 END DO ! End loop on origin level "i"
Note: See TracChangeset
for help on using the changeset viewer.
