- 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
r3712 r3713 188 188 ! ===================================================================== 189 189 190 DO i = minorig + 1, nl 191 IF (ok_entrain) THEN 192 DO j = minorig, nl 193 DO il = 1, ncum 194 IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (j >= (icb(il) - 1)) & 195 .AND. (j <= inb(il))) THEN 196 rti = qta(il, i - 1) - ep(il, i)*clw(il, i) 197 bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 198 IF (cvflag_ice) THEN 199 IF (t(il, j) <= 263.15) THEN 200 bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & 201 lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 202 END IF 203 END IF 204 anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j)) 205 denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j) 206 dei = denom 207 IF (abs(dei) < 0.01) dei = 0.01 208 Sij(il, i, j) = anum/dei 209 Sij(il, i, i) = 1.0 210 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 211 altem = altem/bf2 212 cwat = clw(il, j)*(1.-ep(il, j)) 213 stemp = Sij(il, i, j) 214 IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN 215 IF (cvflag_ice) THEN 216 anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2) 217 denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti) 218 ELSE 219 anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) 220 denom = denom + lv(il, j)*(rr(il, i) - rti) 221 END IF 222 IF (abs(denom) < 0.01) denom = 0.01 223 Sij(il, i, j) = anum/denom 224 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 225 altem = altem - (bf2 - 1.)*cwat 226 END IF 227 IF (Sij(il, i, j) > 0.0) THEN 228 Ment(il, i, j) = 1. 229 elij(il, i, j) = altem 230 elij(il, i, j) = amax1(0.0, elij(il, i, j)) 231 nent(il, i) = nent(il, i) + 1 232 END IF 233 234 Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) 235 Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) 236 ELSE IF (j > i) THEN 237 IF (prt_level >= 10) THEN 238 print *, 'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il, i, j) 239 ENDIF 240 END IF ! new 241 END DO 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 200 bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & 201 lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) 202 END IF 203 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 = denom 206 IF (abs(dei) < 0.01) dei = 0.01 207 Sij(il, i, j) = anum/dei 208 Sij(il, i, i) = 1.0 209 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) 210 altem = altem/bf2 211 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) THEN 214 IF (cvflag_ice) THEN 215 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 ELSE 218 anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) 219 denom = denom + lv(il, j)*(rr(il, i) - rti) 220 END IF 221 IF (abs(denom) < 0.01) denom = 0.01 222 Sij(il, i, j) = anum/denom 223 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)) 242 235 END DO 243 ELSE ! (ok_entrain) 244 DO il = 1, ncum 245 nent(il, i) = 0 246 ENDDO 247 ENDIF ! (ok_entrain) 248 249 IF (prt_level >= 10) THEN 250 print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout) 251 IF (nent(igout, i) .gt. 0) THEN 252 print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout)) 253 ENDIF 254 ENDIF 255 256 ! *** if no air can entrain at level i assume that updraft detrains *** 257 ! *** at that level and calculate detrained air flux and properties *** 258 259 DO il = 1, ncum 260 IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (nent(il, i) == 0)) THEN 236 IF ( nent(il, i) == 0) THEN 261 237 Ment(il, i, i) = 1. 262 238 Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) … … 270 246 END IF 271 247 END DO 272 END DO ! i = minorig + 1, nl 273 248 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 263 END DO 264 END DO 265 ENDIF 266 267 DO i = minorig + 1, nl 274 268 DO j = minorig, nl 275 DO i = minorig, nl 276 DO il = 1, ncum 277 IF ((j >= (icb(il) - 1)) .AND. (j <= inb(il)) .AND. & 278 (i >= icb(il)) .AND. (i <= inb(il))) THEN 279 Sigij(il, i, j) = Sij(il, i, j) 280 END IF 281 END DO 282 END DO 283 END DO 269 IF (prt_level >= 10) THEN 270 print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout) 271 !IF (nent(igout, i) .gt. 0) THEN 272 ! print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout)) 273 !ENDIF 274 ENDIF 275 END DO 276 END DO 277 278 Sigij(1:ncum, 1:nd, 1:nd) = 0.0 279 DO il = 1, ncum 280 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 DO 284 END DO 285 END DO 286 284 287 285 288 ! =====================================================================
Note: See TracChangeset
for help on using the changeset viewer.
