Changeset 3716 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
r3715 r3716 108 108 109 109 !local variables: 110 INTEGER i, j, k, il, im, jm 111 REAL Smid 112 REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp 113 REAL :: alt, delp, delm 110 INTEGER i, j, k, il 114 111 REAL, DIMENSION(nloc, nd) :: ASij 115 112 REAL, DIMENSION(nloc, nd, nd) :: Sij 116 REAL :: awat 117 REAL :: cpm !Mixed draught heat capacity 118 REAL :: Tm !Mixed draught temperature 119 120 REAL amxupcrit, df, ff 121 INTEGER nstep 122 123 INTEGER, SAVE :: igout = 1 113 114 INTEGER, SAVE :: igout = 1 124 115 !$omp THREADPRIVATE(igout) 125 116 … … 183 174 DO i = icb(il), inb(il) 184 175 DO j = (icb(il) - 1), inb(il) 185 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 186 bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 187 IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN 188 bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & 189 lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 190 END IF 191 anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j)) 176 block 177 real :: rti, bf2, anum, denom, dei, altem, cwat, stemp 178 179 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 180 bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 181 IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN 182 bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & 183 lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 184 END IF 185 anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j)) 192 186 denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j) 193 187 dei = denom … … 198 192 altem = altem/bf2 199 193 cwat = clw(il, j)*(1.-ep(il, j)) 200 stemp = Sij(il, i, j) 201 IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN 202 IF (cvflag_ice) THEN 203 anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2) 204 denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti) 205 ELSE 206 anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) 207 denom = denom + lv(il, j)*(rr(il, i) - rti) 194 stemp = Sij(il, i, j) 195 IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN 196 IF (cvflag_ice) THEN 197 anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2) 198 denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti) 199 ELSE 200 anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) 201 denom = denom + lv(il, j)*(rr(il, i) - rti) 202 END IF 203 IF (abs(denom) < 0.01) denom = 0.01 204 Sij(il, i, j) = anum/denom 205 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 206 altem = altem - (bf2 - 1.)*cwat 208 207 END IF 209 IF (abs(denom) < 0.01) denom = 0.01 210 Sij(il, i, j) = anum/denom 211 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 212 altem = altem - (bf2 - 1.)*cwat 213 END IF 214 IF (Sij(il, i, j) > 0.0) THEN 215 Ment(il, i, j) = 1. 216 elij(il, i, j) = altem 217 elij(il, i, j) = amax1(0.0, elij(il, i, j)) 218 nent(il, i) = nent(il, i) + 1 219 END IF 220 221 Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) 222 Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) 208 IF (Sij(il, i, j) > 0.0) THEN 209 Ment(il, i, j) = 1. 210 elij(il, i, j) = altem 211 elij(il, i, j) = amax1(0.0, elij(il, i, j)) 212 nent(il, i) = nent(il, i) + 1 213 END IF 214 215 Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) 216 Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) 217 end block 223 218 END DO 224 219 IF (nent(il, i) == 0) THEN … … 280 275 DO i = icb(il), inb(il) !Loop on origin level "i" 281 276 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 277 real :: signhpmh, Scrit 278 real :: Sjmin(icb(il)-1:inb(il)), Sjmax(icb(il)-1:inb(il)), Qmixmax(icb(il)-1:inb(il)), Qmixmin(icb(il)-1:inb(il)), Rmixmax(icb(il)-1:inb(il)), Rmixmin(icb(il)-1:inb(il)) 279 280 block 281 real :: rti, anum, denom, alt, Sx 282 283 signhpmh = sign(1., hp(il, i) - h(il, i)) 284 285 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 286 IF (cvflag_ice) THEN 287 anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & 288 (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 289 denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & 290 (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 291 ELSE 292 anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & 293 (cpv - cpd)*t(il, i)*(rti - rr(il, i)) 294 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & 295 (cpd - cpv)*t(il, i)*(rr(il, i) - rti) 296 END IF 297 IF (abs(denom) < 0.01) denom = 0.01 298 Scrit = min(anum/denom, 1.) 299 alt = rti - rs(il, i) + Scrit*(rr(il, i) - rti) 300 301 if (alt <= 0) then 302 Scrit = 1.0 303 else 304 ! Get max of Sij 305 Sx = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh) 306 ! Find new critical value Scrit 307 ! such that : Sij > Scrit => mixed draught will detrain at J<I 308 ! Sij < Scrit => mixed draught will detrain at J>I 309 Scrit = max(0., min(Scrit, Sx*max(0., -signhpmh)) + Scrit*max(0., signhpmh)) 310 endif 311 312 end block 313 316 314 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 315 316 ! Compute all Sjmax(j) and Sjmin(j) 317 block 318 real :: smin, smax, sup, Sbef, Smid 319 smin = 1 320 smax = 0.0 321 sup = 0. ! upper S-value reached by descending draughts 322 ! Glitchy : why so complicated? 323 if (i < inb(il)) then 324 Sbef = Sij(il, i, inb(il)) 325 else 326 Sbef = max(0., signhpmh) 327 endif 328 329 DO j = (icb(il) - 1), i - 1 330 IF (Sij(il, i, j) > 0.0) THEN 331 Smid = max(Sij(il, i, j), Scrit) 332 Sjmax(j) = Smid 333 Sjmin(j) = Smid 334 IF (Smid > smax .AND. Sij(il, i, j + 1) > Smid) THEN 335 smax = Smid 336 Sjmax(j) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) 337 Sjmax(j) = max(Sjmax(j), Scrit) 338 Sjmin(j) = min((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) 339 Sjmin(j) = max(Sjmin(j), Scrit) 340 Sbef = Sij(il, i, j) 341 END IF 342 IF (abs(Sjmin(j) - Sjmax(j)) > 1.E-10) & 343 sup = max(Sjmin(j), Sjmax(j), sup) 344 ENDIF 345 END DO 346 j = i 347 IF (Sij(il, i, j) > 0.0) THEN 348 Smid = 1. 349 Sjmin(j) = max((Sij(il, i, j - 1) + Smid)/2., Scrit)*max(0., -signhpmh) + & 350 min((Sij(il, i, j + 1) + Smid)/2., Scrit)*max(0., signhpmh) 351 Sjmin(j) = max(Sjmin(j), sup) 352 Sjmax(j) = 1. 353 ! preparation des variables Scrit, Smin et Sbef pour la partie j>i 354 Scrit = min(Sjmin(j), Sjmax(j), Scrit) 355 Sbef = max(0., signhpmh) 356 supmax(il, i) = sign(Scrit, -signhpmh) 357 ENDIF 358 DO j = i + 1, inb(il) 359 IF (Sij(il, i, j) > 0.0) THEN 360 Smid = min(Sij(il, i, j), Scrit) 361 Sjmax(j) = Smid 362 Sjmin(j) = Smid 363 IF (Smid < smin .AND. Sij(il, i, j + 1) < Smid) THEN 364 smin = min(smin, Smid) 365 Sjmax(j) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit) 366 Sjmin(j) = max((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) 367 Sjmin(j) = min(Sjmin(j), Scrit) 368 Sbef = Sij(il, i, j) 369 END IF 370 ENDIF 371 END DO 372 end block 373 374 DO j = (icb(il) - 1), inb(il) 375 Qmixmax(j) = Qmix(Sjmax(j)) 376 Qmixmin(j) = Qmix(Sjmin(j)) 377 Rmixmax(j) = Rmix(Sjmax(j)) 378 Rmixmin(j) = Rmix(Sjmin(j)) 379 END DO 325 380 326 381 DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" 327 382 IF (Sij(il, i, j) > 0.0) THEN 328 383 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 ! ----------------------------------------------- 384 real :: sqmrmax, sqmrmin, cpm, awat, Tm, rti 370 385 371 386 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) 387 388 sqmrmax = Sjmax(j)*Qmixmax(j) - Rmixmax(j) 389 sqmrmin = Sjmin(j)*Qmixmin(j) - Rmixmin(j) 390 391 Ment(il, i, j) = abs(Qmixmax(j) - Qmixmin(j))*Ment(il, i, j) 392 IF (abs(Qmixmax(j) - Qmixmin(j)) > 1.E-10) THEN 393 Sigij(il, i, j) = (sqmrmax - sqmrmin)/(Qmixmax(j) - Qmixmin(j)) 382 394 ELSE 383 395 Sigij(il, i, j) = 0. … … 427 439 428 440 ! 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)441 ASij(il, i) = ASij(il, i) + abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - & 442 Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j)) 431 443 432 444 ! If I=J (detrainement and entrainement at the same level), then only the … … 434 446 IF (i == j) THEN 435 447 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)448 Ment(il, i, i) = abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - & 449 Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j)) 438 450 Qent(il, i, i) = rti 439 451 uent(il, i, i) = unk(il)
Note: See TracChangeset
for help on using the changeset viewer.
