- Timestamp:
- Jun 11, 2020, 11:09:42 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90
r3713 r3714 126 126 REAL :: cpm !Mixed draught heat capacity 127 127 REAL :: Tm !Mixed draught temperature 128 LOGICAL, DIMENSION(nloc) :: lwork129 128 130 129 REAL amxupcrit, df, ff … … 188 187 ! ===================================================================== 189 188 190 191 IF (ok_entrain) THEN 192 ! NOTE : (04/2020 AD) this loop order gave good results in cv3a_tracer 193 ! if( icb(il) <= i <= inb(il) ) is probably more expensive than reading non-contiguously 194 DO il = 1, ncum 195 DO i = icb(il), inb(il) 196 DO j = (icb(il) - 1), inb(il) 197 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 198 bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 199 IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN 189 IF (ok_entrain) THEN 190 ! NOTE : (04/2020 AD) this loop order gave good results in cv3a_tracer 191 ! if( icb(il) <= i <= inb(il) ) is probably more expensive than reading non-contiguously 192 DO il = 1, ncum 193 DO i = icb(il), inb(il) 194 DO j = (icb(il) - 1), inb(il) 195 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 196 bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 197 IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN 200 198 bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & 201 199 lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 202 END IF203 anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j))204 denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j)205 dei = denom206 IF (abs(dei) < 0.01) dei = 0.01207 Sij(il, i, j) = anum/dei208 Sij(il, i, i) = 1.0209 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j)210 altem = altem/bf2211 cwat = clw(il, j)*(1.-ep(il, j))212 stemp = Sij(il, i, j)213 IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN214 IF (cvflag_ice) THEN215 anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2)216 denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti)217 ELSE218 anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2)219 denom = denom + lv(il, j)*(rr(il, i) - rti)220 200 END IF 221 IF (abs(denom) < 0.01) denom = 0.01 222 Sij(il, i, j) = anum/denom 201 anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j)) 202 denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j) 203 dei = denom 204 IF (abs(dei) < 0.01) dei = 0.01 205 Sij(il, i, j) = anum/dei 206 Sij(il, i, i) = 1.0 223 207 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 224 altem = altem - (bf2 - 1.)*cwat 225 END IF 226 IF (Sij(il, i, j) > 0.0) THEN 227 Ment(il, i, j) = 1. 228 elij(il, i, j) = altem 229 elij(il, i, j) = amax1(0.0, elij(il, i, j)) 230 nent(il, i) = nent(il, i) + 1 231 END IF 232 233 Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) 234 Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) 208 altem = altem/bf2 209 cwat = clw(il, j)*(1.-ep(il, j)) 210 stemp = Sij(il, i, j) 211 IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN 212 IF (cvflag_ice) THEN 213 anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2) 214 denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti) 215 ELSE 216 anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) 217 denom = denom + lv(il, j)*(rr(il, i) - rti) 218 END IF 219 IF (abs(denom) < 0.01) denom = 0.01 220 Sij(il, i, j) = anum/denom 221 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 222 altem = altem - (bf2 - 1.)*cwat 223 END IF 224 IF (Sij(il, i, j) > 0.0) THEN 225 Ment(il, i, j) = 1. 226 elij(il, i, j) = altem 227 elij(il, i, j) = amax1(0.0, elij(il, i, j)) 228 nent(il, i) = nent(il, i) + 1 229 END IF 230 231 Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) 232 Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) 233 END DO 234 IF (nent(il, i) == 0) THEN 235 Ment(il, i, i) = 1. 236 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 237 uent(il, i, i) = unk(il) 238 vent(il, i, i) = vnk(il) 239 IF (fl_cor_ebil .GE. 2) THEN 240 hent(il, i, i) = hp(il, i) 241 ENDIF 242 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 243 Sij(il, i, i) = 0.0 244 END IF 235 245 END DO 236 IF ( nent(il, i) == 0) THEN 237 Ment(il, i, i) = 1. 238 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 239 uent(il, i, i) = unk(il) 240 vent(il, i, i) = vnk(il) 241 IF (fl_cor_ebil .GE. 2) THEN 242 hent(il, i, i) = hp(il, i) 243 ENDIF 244 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 245 Sij(il, i, i) = 0.0 246 END IF 246 END DO 247 ELSE 248 DO i = minorig + 1, nl 249 DO il = 1, ncum 250 IF ((i >= icb(il)) .AND. (i <= inb(il))) THEN 251 Ment(il, i, i) = 1. 252 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 253 uent(il, i, i) = unk(il) 254 vent(il, i, i) = vnk(il) 255 IF (fl_cor_ebil .GE. 2) THEN 256 hent(il, i, i) = hp(il, i) 257 ENDIF 258 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 259 Sij(il, i, i) = 0.0 260 END IF 261 END DO 262 END DO 263 ENDIF 264 265 DO i = minorig + 1, nl 266 DO j = minorig, nl 267 IF (prt_level >= 10) THEN 268 print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout) 269 !IF (nent(igout, i) .gt. 0) THEN 270 ! print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout)) 271 !ENDIF 272 ENDIF 247 273 END DO 248 274 END DO 249 ELSE 250 DO i = minorig + 1, nl 251 DO il = 1, ncum 252 IF ((i >= icb(il)) .AND. (i <= inb(il))) THEN 253 Ment(il, i, i) = 1. 254 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 255 uent(il, i, i) = unk(il) 256 vent(il, i, i) = vnk(il) 257 IF (fl_cor_ebil .GE. 2) THEN 258 hent(il, i, i) = hp(il, i) 259 ENDIF 260 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 261 Sij(il, i, i) = 0.0 262 END IF 275 276 Sigij(1:ncum, 1:nd, 1:nd) = 0.0 277 DO il = 1, ncum 278 DO i = icb(il), inb(il) 279 DO j = (icb(il) - 1), inb(il) 280 Sigij(il, i, j) = Sij(il, i, j) 281 END DO 263 282 END DO 264 283 END DO 265 ENDIF266 267 DO i = minorig + 1, nl268 DO j = minorig, nl269 IF (prt_level >= 10) THEN270 print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout)271 !IF (nent(igout, i) .gt. 0) THEN272 ! print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout))273 !ENDIF274 ENDIF275 END DO276 END DO277 278 Sigij(1:ncum, 1:nd, 1:nd) = 0.0279 DO il = 1, ncum280 DO i = icb(il), inb(il)281 DO j = (icb(il) - 1), inb(il)282 Sigij(il, i, j) = Sij(il, i, j)283 END DO284 END DO285 END DO286 287 284 288 285 ! ===================================================================== … … 292 289 293 290 DO il = 1, ncum 294 lwork(il) = .FALSE.295 END DO296 297 DO il = 1, ncum298 291 DO i = icb(il), inb(il) !Loop on origin level "i" 299 292 signhpmh(il) = sign(1., hp(il, i) - h(il, i)) … … 301 294 Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il)) 302 295 303 lwork(il) = (nent(il, i) /= 0)304 296 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 305 297 IF (cvflag_ice) THEN … … 341 333 342 334 DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" 343 IF ( lwork(il) .AND.Sij(il, i, j) > 0.0) THEN335 IF (Sij(il, i, j) > 0.0) THEN 344 336 ! ----------------------------------------------- 345 337 IF (j > i) THEN … … 467 459 DO il = 1, ncum 468 460 DO i = icb(il), inb(il) !Loop on origin level "i" 469 IF (lwork(il)) THEN 470 csum(il, i) = 0 471 DO j = (icb(il) - 1), inb(il) 472 ASij(il, i) = amax1(1.0E-16, ASij(il, i)) 473 ASij_inv(il) = 1.0/ASij(il, i) 474 ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs 475 IF (ASij_inv(il) > 100.) ASij_inv(il) = 0. 476 Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) 477 csum(il, i) = csum(il, i) + Ment(il, i, j) 478 END DO 479 IF (csum(il, i) < 1.) THEN 480 nent(il, i) = 0 481 Ment(il, i, i) = 1. 482 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) 483 uent(il, i, i) = unk(il) 484 vent(il, i, i) = vnk(il) 485 elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) 486 IF (fl_cor_ebil .GE. 2) THEN 487 hent(il, i, i) = hp(il, i) 488 Sigij(il, i, i) = 0.0 489 ELSE 490 Sij(il, i, i) = 0.0 491 ENDIF 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 492 482 ENDIF 493 END 483 ENDIF 494 484 END DO ! Loop il = 1, ncum 495 485 END DO ! End loop on origin level "i"
Note: See TracChangeset
for help on using the changeset viewer.
