- 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
r3711 r3712 119 119 REAL, DIMENSION(nloc) :: Sjmin, Sjmax 120 120 REAL, DIMENSION(nloc) :: Sbef, sup 121 REAL, DIMENSION(nloc) :: ASij, ASij_inv, smax, Scrit 121 REAL, DIMENSION(nloc, nd) :: ASij 122 REAL, DIMENSION(nloc) :: ASij_inv, smax, Scrit 122 123 REAL, DIMENSION(nloc, nd, nd) :: Sij 123 124 REAL, DIMENSION(nloc, nd) :: csum … … 131 132 132 133 INTEGER, SAVE :: igout = 1 133 !$omp THREADPRIVATE(igout)134 135 ! -- Mixing probability distribution functions134 !$omp THREADPRIVATE(igout) 135 136 ! -- Mixing probability distribution functions 136 137 137 138 REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F … … 149 150 150 151 INTEGER, SAVE :: ifrst = 0 151 !$omp THREADPRIVATE(ifrst)152 153 ! =====================================================================154 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS155 ! =====================================================================156 ! -- Initialize mixing PDF coefficients152 !$omp THREADPRIVATE(ifrst) 153 154 ! ===================================================================== 155 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS 156 ! ===================================================================== 157 ! -- Initialize mixing PDF coefficients 157 158 IF (ifrst == 0) THEN 158 159 ifrst = 1 … … 181 182 Sigij(1:ncum, 1:nd, 1:nd) = 0.0 182 183 183 ! =====================================================================184 ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING185 ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING186 ! --- FRACTION (Sij)187 ! =====================================================================184 ! ===================================================================== 185 ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING 186 ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING 187 ! --- FRACTION (Sij) 188 ! ===================================================================== 188 189 189 190 DO i = minorig + 1, nl … … 253 254 ENDIF 254 255 255 ! *** if no air can entrain at level i assume that updraft detrains ***256 ! *** at that level and calculate detrained air flux and properties ***256 ! *** if no air can entrain at level i assume that updraft detrains *** 257 ! *** at that level and calculate detrained air flux and properties *** 257 258 258 259 DO il = 1, ncum … … 282 283 END DO 283 284 284 ! =====================================================================285 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES286 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING287 ! =====================================================================285 ! ===================================================================== 286 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES 287 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING 288 ! ===================================================================== 288 289 289 290 DO il = 1, ncum … … 291 292 END DO 292 293 293 ! --------------------------------------------------------------- 294 DO i = minorig + 1, nl !Loop on origin level "i" 295 ! --------------------------------------------------------------- 296 297 DO il = 1, ncum 298 IF (i >= icb(il) .AND. i <= inb(il)) THEN 299 signhpmh(il) = sign(1., hp(il, i) - h(il, i)) 300 301 Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il)) 302 303 lwork(il) = (nent(il, i) /= 0) 304 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 305 IF (cvflag_ice) THEN 306 anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & 307 (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 308 denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & 309 (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 310 ELSE 311 anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & 312 (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 313 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & 314 (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 315 END IF 316 IF (abs(denom) < 0.01) denom = 0.01 317 Scrit(il) = min(anum/denom, 1.) 318 alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti) 319 320 ! Find new critical value Scrit2 321 ! such that : Sij > Scrit2 => mixed draught will detrain at J<I 322 ! Sij < Scrit2 => mixed draught will detrain at J>I 323 Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + & 324 Scrit(il)*max(0., signhpmh(il)) 325 Scrit(il) = Scrit2 326 327 ! Correction pour la nouvelle logique; la correction pour ALT 328 ! est un peu au hazard 329 IF (Scrit(il) <= 0.0) Scrit(il) = 0.0 330 IF (alt <= 0.0) Scrit(il) = 1.0 331 332 smax(il) = 0.0 333 ASij(il) = 0.0 334 sup(il) = 0. ! upper S-value reached by descending draughts 335 336 if (i < inb(il)) then 337 Sbef(il) = Sij(il, i, inb(il)) 338 else 339 Sbef(il) = max(0., signhpmh(il)) 340 endif 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)) 294 DO il = 1, ncum 295 DO i = icb(il), inb(il) !Loop on origin level "i" 296 signhpmh(il) = sign(1., hp(il, i) - h(il, i)) 297 298 Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il)) 299 300 lwork(il) = (nent(il, i) /= 0) 301 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 302 IF (cvflag_ice) THEN 303 anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & 304 (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 305 denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & 306 (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 307 ELSE 308 anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & 309 (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 310 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & 311 (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 312 END IF 313 IF (abs(denom) < 0.01) denom = 0.01 314 Scrit(il) = min(anum/denom, 1.) 315 alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti) 316 317 ! Find new critical value Scrit2 318 ! such that : Sij > Scrit2 => mixed draught will detrain at J<I 319 ! Sij < Scrit2 => mixed draught will detrain at J>I 320 Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + & 321 Scrit(il)*max(0., signhpmh(il)) 322 Scrit(il) = Scrit2 323 324 ! Correction pour la nouvelle logique; la correction pour ALT 325 ! est un peu au hazard 326 IF (Scrit(il) <= 0.0) Scrit(il) = 0.0 327 IF (alt <= 0.0) Scrit(il) = 1.0 328 329 smax(il) = 0.0 330 ASij(il, i) = 0.0 331 sup(il) = 0. ! upper S-value reached by descending draughts 332 333 if (i < inb(il)) then 334 Sbef(il) = Sij(il, i, inb(il)) 335 else 336 Sbef(il) = max(0., signhpmh(il)) 337 endif 338 339 DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" 340 IF (lwork(il) .AND. Sij(il, i, j) > 0.0) THEN 341 ! ----------------------------------------------- 342 IF (j > i) THEN 343 Smid = min(Sij(il, i, j), Scrit(il)) 344 Sjmax(il) = Smid 345 Sjmin(il) = Smid 346 IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN 347 smin(il) = min(smin(il), Smid) 348 Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il)) 349 Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 350 Sjmin(il) = min(Sjmin(il), Scrit(il)) 351 Sbef(il) = Sij(il, i, j) 382 352 END IF 383 ! ----------------------------------------------- 384 353 ELSE IF (j == i) THEN 354 Smid = 1. 355 Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + & 356 min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il)) 357 Sjmin(il) = max(Sjmin(il), sup(il)) 358 Sjmax(il) = 1. 359 ! preparation des variables Scrit, Smin et Sbef pour la partie j>i 360 Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il)) 361 362 smin(il) = 1. 363 Sbef(il) = max(0., signhpmh(il)) 364 supmax(il, i) = sign(Scrit(il), -signhpmh(il)) 365 ELSE IF (j < i) THEN 366 Smid = max(Sij(il, i, j), Scrit(il)) 367 Sjmax(il) = Smid 368 Sjmin(il) = Smid 369 IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN 370 smax(il) = Smid 371 Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) 372 Sjmax(il) = max(Sjmax(il), Scrit(il)) 373 Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) 374 Sjmin(il) = max(Sjmin(il), Scrit(il)) 375 Sbef(il) = Sij(il, i, j) 376 END IF 377 IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) & 378 sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) 379 END IF 380 ! ----------------------------------------------- 381 382 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 383 Qmixmax(il) = Qmix(Sjmax(il)) 384 Qmixmin(il) = Qmix(Sjmin(il)) 385 Rmixmax(il) = Rmix(Sjmax(il)) 386 Rmixmin(il) = Rmix(Sjmin(il)) 387 sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il)) 388 sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il)) 389 390 Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j) 391 IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN 392 Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il)) 393 ELSE 394 Sigij(il, i, j) = 0. 395 END IF 396 397 ! Compute Qent, uent, vent according to the true mixing fraction 398 Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) 399 uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) 400 vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) 401 402 ! Compute liquid water static energy of mixed draughts 403 hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) 404 ! Heat capacity of mixed draught 405 cpm = cpd + Qent(il, i, j)*(cpv - cpd) 406 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 407 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 408 elij(il, i, j) = elij(il, i, j) + & 409 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 410 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 411 elij(il, i, j) = elij(il, i, j)/ & 412 (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & 413 (cpm*rrv*t(il, j)*t(il, j))) 414 ELSE 415 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 416 elij(il, i, j) = elij(il, i, j) + & 417 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & 418 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) 419 elij(il, i, j) = elij(il, i, j)/ & 420 (1.+lv(il, j)*lv(il, j)*rs(il, j)/ & 421 (cpm*rrv*t(il, j)*t(il, j))) 422 ENDIF 423 elij(il, i, j) = max(elij(il, i, j), 0.) 424 425 elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) 426 427 IF (j > i) THEN 428 awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) 429 awat = amax1(awat, 0.0) 430 ELSE 431 awat = 0. 432 END IF 433 434 ! Mixed draught temperature at level j 435 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)) 436 IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN 437 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat 438 ELSE 439 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat 440 ENDIF 441 442 ! ASij is the integral of P(F) over the relevant F interval 443 ASij(il, i) = ASij(il, i) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 444 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 445 446 ! If I=J (detrainement and entrainement at the same level), then only the 447 ! adiabatic ascent part of the mixture is considered 448 IF (i == j) THEN 385 449 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. 398 END IF 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 462 END IF 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) 450 Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & 451 Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) 452 Qent(il, i, i) = rti 479 453 uent(il, i, i) = unk(il) 480 454 vent(il, i, i) = vnk(il) 455 hent(il, i, i) = hp(il, i) 481 456 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 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 457 Sigij(il, i, i) = 0. 458 END IF 459 END IF 460 END DO ! Loop j = (icb(il) - 1), inb(il) 461 END DO ! Loop i = icb(il), inb(il) 462 END DO ! Loop il = 1, ncum 463 464 DO il = 1, ncum 465 DO i = icb(il), inb(il) !Loop on origin level "i" 466 IF (lwork(il)) THEN 467 csum(il, i) = 0 468 DO j = (icb(il) - 1), inb(il) 469 ASij(il, i) = amax1(1.0E-16, ASij(il, i)) 470 ASij_inv(il) = 1.0/ASij(il, i) 471 ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs 472 IF (ASij_inv(il) > 100.) ASij_inv(il) = 0. 473 Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) 474 csum(il, i) = csum(il, i) + Ment(il, i, j) 475 END DO 476 IF (csum(il, i) < 1.) THEN 477 nent(il, i) = 0 478 Ment(il, i, i) = 1. 479 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 480 uent(il, i, i) = unk(il) 481 vent(il, i, i) = vnk(il) 482 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 483 IF (fl_cor_ebil .GE. 2) THEN 484 hent(il, i, i) = hp(il, i) 485 Sigij(il, i, i) = 0.0 486 ELSE 487 Sij(il, i, i) = 0.0 488 488 ENDIF 489 END 490 END IF ! i >= icb(il) .AND. i <= inb(il)489 ENDIF 490 END IF 491 491 END DO ! Loop il = 1, ncum 492 492 END DO ! End loop on origin level "i"
Note: See TracChangeset
for help on using the changeset viewer.
