Changeset 2056 for LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F90
- Timestamp:
- Jun 11, 2014, 3:46:46 PM (10 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/cv3p_mixing.F90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 1998,2000-2023,2025-2029,2032,2034,2036-2049,2051-2055
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F90
r1999 r2056 1 SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, & 2 u, v, tra, h, lv, qnk, unk, vnk, hp, tv, tvp, ep, clw, sig, ment, qent, & 3 hent, uent, vent, nent, sigij, elij, supmax, ments, qents, traent) 4 ! ************************************************************** 5 ! * 6 ! CV3P_MIXING : compute mixed draught properties and, * 7 ! within a scaling factor, mixed draught * 8 ! mass fluxes. * 9 ! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15* 10 ! modified by : * 11 ! ************************************************************** 1 SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, & 2 ph, t, rr, rs, u, v, tra, h, lv, qnk, & 3 unk, vnk, hp, tv, tvp, ep, clw, sig, & 4 Ment, Qent, hent, uent, vent, nent, & 5 Sigij, elij, supmax, Ments, Qents, traent) 6 ! ************************************************************** 7 ! * 8 ! CV3P_MIXING : compute mixed draught properties and, * 9 ! within a scaling factor, mixed draught * 10 ! mass fluxes. * 11 ! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15* 12 ! modified by : * 13 ! ************************************************************** 12 14 13 15 IMPLICIT NONE … … 17 19 include "YOMCST2.h" 18 20 19 ! inputs: 20 INTEGER ncum, nd, na, ntra, nloc 21 INTEGER icb(nloc), inb(nloc), nk(nloc) 22 REAL sig(nloc, nd) 23 REAL qnk(nloc), unk(nloc), vnk(nloc) 24 REAL ph(nloc, nd+1) 25 REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd) 26 REAL u(nloc, nd), v(nloc, nd) 27 REAL tra(nloc, nd, ntra) ! input of convect3 28 REAL lv(nloc, na) 29 REAL h(nloc, na) !liquid water static energy of environment 30 REAL hp(nloc, na) !liquid water static energy of air shed from adiab. asc. 31 REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na) 32 33 ! outputs: 34 REAL ment(nloc, na, na), qent(nloc, na, na) 35 REAL uent(nloc, na, na), vent(nloc, na, na) 36 REAL sigij(nloc, na, na), elij(nloc, na, na) 37 REAL supmax(nloc, na) ! Highest mixing fraction of mixed updraughts 38 ! with the sign of (h-hp) 39 REAL traent(nloc, nd, nd, ntra) 40 REAL ments(nloc, nd, nd), qents(nloc, nd, nd) 41 REAL hent(nloc, nd, nd) 42 INTEGER nent(nloc, nd) 43 44 ! local variables: 21 !inputs: 22 INTEGER, INTENT (IN) :: ncum, nd, na 23 INTEGER, INTENT (IN) :: ntra, nloc 24 INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk 25 REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig 26 REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk 27 REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph 28 REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs 29 REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v 30 REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3 31 REAL, DIMENSION (nloc, na), INTENT (IN) :: lv 32 REAL, DIMENSION (nloc, na), INTENT (IN) :: h !liquid water static energy of environMent 33 REAL, DIMENSION (nloc, na), INTENT (IN) :: hp !liquid water static energy of air shed from adiab. asc. 34 REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp 35 REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, clw 36 37 !outputs: 38 REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Ment, Qent 39 REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent 40 REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Sigij, elij 41 REAL, DIMENSION (nloc, na), INTENT (OUT) :: supmax(nloc, na) ! Highest mixing fraction of mixed 42 ! updraughts with the sign of (h-hp) 43 REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent 44 REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: Ments, Qents 45 REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: hent 46 INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent 47 48 !local variables: 45 49 INTEGER i, j, k, il, im, jm 46 50 INTEGER num1, num2 47 REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp 48 REAL alt, delp, delm 49 REAL qmixmax(nloc), rmixmax(nloc), sqmrmax(nloc) 50 REAL qmixmin(nloc), rmixmin(nloc), sqmrmin(nloc) 51 REAL signhpmh(nloc) 52 REAL sx(nloc), scrit2 53 REAL smid(nloc), sjmin(nloc), sjmax(nloc) 54 REAL sbef(nloc), sup(nloc), smin(nloc) 55 REAL asij(nloc), smax(nloc), scrit(nloc) 56 REAL sij(nloc, nd, nd) 57 REAL csum(nloc, nd) 58 REAL awat 59 LOGICAL lwork(nloc) 51 REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp, qp 52 REAL :: alt, delp, delm 53 REAL, DIMENSION (nloc) :: Qmixmax, Rmixmax, sqmrmax 54 REAL, DIMENSION (nloc) :: Qmixmin, Rmixmin, sqmrmin 55 REAL, DIMENSION (nloc) :: signhpmh 56 REAL, DIMENSION (nloc) :: Sx 57 REAL :: Scrit2 58 REAL, DIMENSION (nloc) :: Smid, Sjmin, Sjmax 59 REAL, DIMENSION (nloc) :: Sbef, sup, smin 60 REAL, DIMENSION (nloc) :: ASij, smax, Scrit 61 REAL, DIMENSION (nloc, nd, nd) :: Sij 62 REAL, DIMENSION (nloc, nd) :: csum 63 REAL :: awat 64 LOGICAL, DIMENSION (nloc) :: lwork 60 65 61 66 REAL amxupcrit, df, ff 62 67 INTEGER nstep 63 68 64 ! -- Mixing probability distribution functions 65 66 REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f 67 68 qcoef1(f) = tanh(f/gammas) 69 qcoef2(f) = (tanh(f/gammas)+gammas*log(cosh((1.-f)/gammas)/cosh(f/gammas))) 70 qff(f) = max(min(f,1.), 0.) 71 qfff(f) = min(qff(f), scut) 72 qmix1(f) = (tanh((qff(f)-fmax)/gammas)+qcoef1max)/qcoef2max 73 rmix1(f) = (gammas*log(cosh((qff(f)-fmax)/gammas))+qff(f)*qcoef1max)/ & 74 qcoef2max 75 qmix2(f) = -log(1.-qfff(f))/scut 76 rmix2(f) = (qfff(f)+(1.-qff(f))*log(1.-qfff(f)))/scut 77 qmix(f) = qqa1*qmix1(f) + qqa2*qmix2(f) 78 rmix(f) = qqa1*rmix1(f) + qqa2*rmix2(f) 69 ! -- Mixing probability distribution functions 70 71 REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F 72 73 Qcoef1(F) = tanh(F/gammas) 74 Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas))) 75 QFF(F) = max(min(F,1.), 0.) 76 QFFf(F) = min(QFF(F), scut) 77 Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max 78 Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max 79 Qmix2(F) = -log(1.-QFFf(F))/scut 80 Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut 81 Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F) 82 Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F) 79 83 80 84 INTEGER, SAVE :: ifrst 81 85 DATA ifrst/0/ 82 83 84 85 86 87 88 89 86 !$OMP THREADPRIVATE(ifrst) 87 88 89 ! ===================================================================== 90 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS 91 ! ===================================================================== 92 93 ! -- Initialize mixing PDF coefficients 90 94 IF (ifrst==0) THEN 91 95 ifrst = 1 92 qcoef1max = qcoef1(fmax)93 qcoef2max = qcoef2(fmax)96 Qcoef1max = Qcoef1(Fmax) 97 Qcoef2max = Qcoef2(Fmax) 94 98 95 99 END IF 96 100 97 101 98 102 ! ori do 360 i=1,ncum*nlp 99 103 DO j = 1, nl 100 104 DO i = 1, ncum 101 105 nent(i, j) = 0 102 103 104 END DO 105 END DO 106 107 108 106 ! in convect3, m is computed in cv3_closure 107 ! ori m(i,1)=0.0 108 END DO 109 END DO 110 111 ! ori do 400 k=1,nlp 112 ! ori do 390 j=1,nlp 109 113 DO j = 1, nl 110 114 DO k = 1, nl 111 115 DO i = 1, ncum 112 qent(i, k, j) = rr(i, j)116 Qent(i, k, j) = rr(i, j) 113 117 uent(i, k, j) = u(i, j) 114 118 vent(i, k, j) = v(i, j) 115 119 elij(i, k, j) = 0.0 116 120 hent(i, k, j) = 0.0 117 ! AC! ment(i,k,j)=0.0118 ! AC! sij(i,k,j)=0.0119 END DO 120 END DO 121 END DO 122 123 !AC!124 ment(1:ncum, 1:nd, 1:nd) = 0.0125 sij(1:ncum, 1:nd, 1:nd) = 0.0126 !AC!121 !AC! Ment(i,k,j)=0.0 122 !AC! Sij(i,k,j)=0.0 123 END DO 124 END DO 125 END DO 126 127 !AC! 128 Ment(1:ncum, 1:nd, 1:nd) = 0.0 129 Sij(1:ncum, 1:nd, 1:nd) = 0.0 130 !AC! 127 131 128 132 DO k = 1, ntra … … 136 140 END DO 137 141 138 139 ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING140 141 ! --- FRACTION (sij)142 142 ! ===================================================================== 143 ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING 144 ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING 145 ! --- FRACTION (Sij) 146 ! ===================================================================== 143 147 144 148 DO i = minorig + 1, nl … … 146 150 DO j = minorig, nl 147 151 DO il = 1, ncum 148 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &149 1)).AND. (j<=inb(il))) THEN152 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) & 153 .AND. (j<=inb(il))) THEN 150 154 151 155 rti = qnk(il) - ep(il, i)*clw(il, i) … … 155 159 dei = denom 156 160 IF (abs(dei)<0.01) dei = 0.01 157 sij(il, i, j) = anum/dei158 sij(il, i, i) = 1.0159 altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)161 Sij(il, i, j) = anum/dei 162 Sij(il, i, i) = 1.0 163 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j) 160 164 altem = altem/bf2 161 165 cwat = clw(il, j)*(1.-ep(il,j)) 162 stemp = sij(il, i, j)166 stemp = Sij(il, i, j) 163 167 IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN 164 168 anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2) 165 169 denom = denom + lv(il, j)*(rr(il,i)-rti) 166 170 IF (abs(denom)<0.01) denom = 0.01 167 sij(il, i, j) = anum/denom 168 altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - & 169 rs(il, j) 171 Sij(il, i, j) = anum/denom 172 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j) 170 173 altem = altem - (bf2-1.)*cwat 171 174 END IF 172 IF ( sij(il,i,j)>0.0) THEN173 ! cc ment(il,i,j)=m(il,i)174 ment(il, i, j) = 1.175 IF (Sij(il,i,j)>0.0) THEN 176 !!! Ment(il,i,j)=m(il,i) 177 Ment(il, i, j) = 1. 175 178 elij(il, i, j) = altem 176 179 elij(il, i, j) = amax1(0.0, elij(il,i,j)) … … 178 181 END IF 179 182 180 sij(il, i, j) = amax1(0.0, sij(il,i,j))181 sij(il, i, j) = amin1(1.0, sij(il,i,j))183 Sij(il, i, j) = amax1(0.0, Sij(il,i,j)) 184 Sij(il, i, j) = amin1(1.0, Sij(il,i,j)) 182 185 END IF ! new 183 186 END DO … … 185 188 186 189 187 ! *** if no air can entrain at level i assume that updraft detrains 188 ! *** 189 ! *** at that level and calculate detrained air flux and properties 190 ! *** 191 192 193 ! @ do 170 i=icb(il),inb(il) 190 ! *** if no air can entrain at level i assume that updraft detrains *** 191 ! *** at that level and calculate detrained air flux and properties *** 192 193 194 ! @ do 170 i=icb(il),inb(il) 194 195 195 196 DO il = 1, ncum 196 197 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN 197 198 ! cc ment(il,i,i)=m(il,i)199 ment(il, i, i) = 1.200 qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)198 ! @ if(nent(il,i).eq.0)then 199 !!! Ment(il,i,i)=m(il,i) 200 Ment(il, i, i) = 1. 201 Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i) 201 202 uent(il, i, i) = unk(il) 202 203 vent(il, i, i) = vnk(il) 203 204 elij(il, i, i) = clw(il, i)*(1.-ep(il,i)) 204 sij(il, i, i) = 0.0205 Sij(il, i, i) = 0.0 205 206 END IF 206 207 END DO … … 220 221 DO i = minorig, nl 221 222 DO il = 1, ncum 222 IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=&223 inb(il))) THEN224 sigij(il, i, j) = sij(il, i, j)225 END IF 226 END DO 227 END DO 228 END DO 229 230 231 232 233 234 235 236 223 IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. & 224 (i>=icb(il)) .AND. (i<=inb(il))) THEN 225 Sigij(il, i, j) = Sij(il, i, j) 226 END IF 227 END DO 228 END DO 229 END DO 230 ! @ enddo 231 232 ! @170 continue 233 234 ! ===================================================================== 235 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES 236 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING 237 ! ===================================================================== 237 238 238 239 CALL zilch(csum, nloc*nd) … … 242 243 END DO 243 244 244 245 DO i = minorig + 1, nl !Loop on origin level "i"246 245 ! --------------------------------------------------------------- 246 DO i = minorig + 1, nl !Loop on origin level "i" 247 ! --------------------------------------------------------------- 247 248 248 249 num1 = 0 … … 253 254 254 255 255 ! jyg1 Find maximum of SIJ for J>I, if any.256 257 sx(:) = 0.256 !JYG1 Find maximum of SIJ for J>I, if any. 257 258 Sx(:) = 0. 258 259 259 260 DO il = 1, ncum 260 261 IF (i>=icb(il) .AND. i<=inb(il)) THEN 261 262 signhpmh(il) = sign(1., hp(il,i)-h(il,i)) 262 sbef(il) = max(0., signhpmh(il))263 Sbef(il) = max(0., signhpmh(il)) 263 264 END IF 264 265 END DO … … 267 268 DO il = 1, ncum 268 269 IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN 269 IF ( sbef(il)<sij(il,i,j)) THEN270 sx(il) = max(sij(il,i,j), sx(il))271 END IF 272 sbef(il) = sij(il, i, j)270 IF (Sbef(il)<Sij(il,i,j)) THEN 271 Sx(il) = max(Sij(il,i,j), Sx(il)) 272 END IF 273 Sbef(il) = Sij(il, i, j) 273 274 END IF 274 275 END DO … … 279 280 IF (i>=icb(il) .AND. i<=inb(il)) THEN 280 281 lwork(il) = (nent(il,i)/=0) 281 qp= qnk(il) - ep(il, i)*clw(il, i)282 anum = h(il, i) - hp(il, i) - lv(il, i)*( qp-rs(il,i)) + &283 (cpv-cpd)*t(il, i)*(qp-rr(il,i))284 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)- qp) + &285 (cpd-cpv)*t(il, i)*(rr(il,i)-qp)282 rti = qnk(il) - ep(il, i)*clw(il, i) 283 anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + & 284 (cpv-cpd)*t(il, i)*(rti-rr(il,i)) 285 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + & 286 (cpd-cpv)*t(il, i)*(rr(il,i)-rti) 286 287 IF (abs(denom)<0.01) denom = 0.01 287 scrit(il) = min(anum/denom, 1.)288 alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)289 290 ! jyg1 Find new critical value Scrit2291 !such that : Sij > Scrit2 => mixed draught will detrain at J<I292 !Sij < Scrit2 => mixed draught will detrain at J>I293 294 scrit2 = min(scrit(il), sx(il))*max(0., -signhpmh(il)) + &295 scrit(il)*max(0., signhpmh(il))296 297 scrit(il) = scrit2298 299 ! jygCorrection pour la nouvelle logique; la correction pour ALT300 301 IF ( scrit(il)<=0.0) scrit(il) = 0.0302 IF (alt<=0.0) scrit(il) = 1.0288 Scrit(il) = min(anum/denom, 1.) 289 alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti) 290 291 !JYG1 Find new critical value Scrit2 292 ! such that : Sij > Scrit2 => mixed draught will detrain at J<I 293 ! Sij < Scrit2 => mixed draught will detrain at J>I 294 295 Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + & 296 Scrit(il)*max(0., signhpmh(il)) 297 298 Scrit(il) = Scrit2 299 300 !JYG Correction pour la nouvelle logique; la correction pour ALT 301 ! est un peu au hazard 302 IF (Scrit(il)<=0.0) Scrit(il) = 0.0 303 IF (alt<=0.0) Scrit(il) = 1.0 303 304 304 305 smax(il) = 0.0 305 asij(il) = 0.0306 sup(il) = 0. ! upper S-value reached by descending draughts307 END IF 308 END DO 309 310 311 DO j = minorig, nl !Loop on destination level "j"312 306 ASij(il) = 0.0 307 sup(il) = 0. ! upper S-value reached by descending draughts 308 END IF 309 END DO 310 311 ! --------------------------------------------------------------- 312 DO j = minorig, nl !Loop on destination level "j" 313 ! --------------------------------------------------------------- 313 314 314 315 num2 = 0 315 316 DO il = 1, ncum 316 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & 317 il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1 317 IF (i>=icb(il) .AND. i<=inb(il) .AND. & 318 j>=(icb(il)-1) .AND. j<=inb(il) .AND. & 319 lwork(il)) num2 = num2 + 1 318 320 END DO 319 321 IF (num2<=0) GO TO 175 320 322 321 323 ! ----------------------------------------------- 322 324 IF (j>i) THEN 323 325 ! ----------------------------------------------- 324 326 DO il = 1, ncum 325 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb(&326 il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN327 IF (sij(il,i,j)>0.0) THEN328 smid(il) = min(sij(il,i,j), scrit(il))329 sjmax(il) = smid(il)330 sjmin(il) = smid(il)331 IF (smid(il)<smin(il) .AND. sij(il,i,j+1)<smid(il)) THEN332 smin(il) = smid(il)333 s jmax(il) = min((sij(il,i,j+1)+sij(il,i, &334 j))/2., sij(il,i,j), scrit(il))335 sjmin(il) = max((sbef(il)+sij(il,i,j))/2., sij(il,i,j))336 sjmin(il) = min(sjmin(il), scrit(il))337 sbef(il) = sij(il, i, j)327 IF (i>=icb(il) .AND. i<=inb(il) .AND. & 328 j>=(icb(il)-1) .AND. j<=inb(il) .AND. & 329 lwork(il)) THEN 330 IF (Sij(il,i,j)>0.0) THEN 331 Smid(il) = min(Sij(il,i,j), Scrit(il)) 332 Sjmax(il) = Smid(il) 333 Sjmin(il) = Smid(il) 334 IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN 335 smin(il) = Smid(il) 336 Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il)) 337 Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j)) 338 Sjmin(il) = min(Sjmin(il), Scrit(il)) 339 Sbef(il) = Sij(il, i, j) 338 340 END IF 339 341 END IF 340 342 END IF 341 343 END DO 342 344 ! ----------------------------------------------- 343 345 ELSE IF (j==i) THEN 344 346 ! ----------------------------------------------- 345 347 DO il = 1, ncum 346 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & 347 il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN 348 IF (sij(il,i,j)>0.0) THEN 349 smid(il) = 1. 350 sjmin(il) = max((sij(il,i,j-1)+smid(il))/2., scrit(il))*max(0., & 351 -signhpmh(il)) + min((sij(il,i,j+1)+smid(il))/2., scrit(il))* & 352 max(0., signhpmh(il)) 353 sjmin(il) = max(sjmin(il), sup(il)) 354 sjmax(il) = 1. 355 356 ! - preparation des variables Scrit, Smin et Sbef 357 ! pour la partie j>i 358 scrit(il) = min(sjmin(il), sjmax(il), scrit(il)) 348 IF (i>=icb(il) .AND. i<=inb(il) .AND. & 349 j>=(icb(il)-1) .AND. j<=inb(il) .AND. & 350 lwork(il)) THEN 351 IF (Sij(il,i,j)>0.0) THEN 352 Smid(il) = 1. 353 Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + & 354 min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il)) 355 Sjmin(il) = max(Sjmin(il), sup(il)) 356 Sjmax(il) = 1. 357 358 ! - preparation des variables Scrit, Smin et Sbef pour la partie j>i 359 Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il)) 359 360 360 361 smin(il) = 1. 361 sbef(il) = max(0., signhpmh(il))362 supmax(il, i) = sign( scrit(il), -signhpmh(il))363 END IF 364 END IF 365 END DO 366 362 Sbef(il) = max(0., signhpmh(il)) 363 supmax(il, i) = sign(Scrit(il), -signhpmh(il)) 364 END IF 365 END IF 366 END DO 367 ! ----------------------------------------------- 367 368 ELSE IF (j<i) THEN 368 369 ! ----------------------------------------------- 369 370 DO il = 1, ncum 370 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & 371 il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN 372 IF (sij(il,i,j)>0.0) THEN 373 smid(il) = max(sij(il,i,j), scrit(il)) 374 sjmax(il) = smid(il) 375 sjmin(il) = smid(il) 376 IF (smid(il)>smax(il) .AND. sij(il,i,j+1)>smid(il)) THEN 377 smax(il) = smid(il) 378 sjmax(il) = max((sij(il,i,j+1)+sij(il,i,j))/2., sij(il,i,j)) 379 sjmax(il) = max(sjmax(il), scrit(il)) 380 sjmin(il) = min((sbef(il)+sij(il,i,j))/2., sij(il,i,j)) 381 sjmin(il) = max(sjmin(il), scrit(il)) 382 sbef(il) = sij(il, i, j) 371 IF (i>=icb(il) .AND. i<=inb(il) .AND. & 372 j>=(icb(il)-1) .AND. j<=inb(il) .AND. & 373 lwork(il)) THEN 374 IF (Sij(il,i,j)>0.0) THEN 375 Smid(il) = max(Sij(il,i,j), Scrit(il)) 376 Sjmax(il) = Smid(il) 377 Sjmin(il) = Smid(il) 378 IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN 379 smax(il) = Smid(il) 380 Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j)) 381 Sjmax(il) = max(Sjmax(il), Scrit(il)) 382 Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j)) 383 Sjmin(il) = max(Sjmin(il), Scrit(il)) 384 Sbef(il) = Sij(il, i, j) 383 385 END IF 384 IF (abs(sjmin(il)-sjmax(il))>1.E-10) sup(il) = max(sjmin(il), & 385 sjmax(il), sup(il)) 386 END IF 387 END IF 388 END DO 389 ! ----------------------------------------------- 390 END IF 391 ! ----------------------------------------------- 392 393 394 DO il = 1, ncum 395 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & 396 il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN 397 IF (sij(il,i,j)>0.0) THEN 386 IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) & 387 sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) 388 END IF 389 END IF 390 END DO 391 ! ----------------------------------------------- 392 END IF 393 ! ----------------------------------------------- 394 395 396 DO il = 1, ncum 397 IF (i>=icb(il) .AND. i<=inb(il) .AND. & 398 j>=(icb(il)-1) .AND. j<=inb(il) .AND. & 399 lwork(il)) THEN 400 IF (Sij(il,i,j)>0.0) THEN 398 401 rti = qnk(il) - ep(il, i)*clw(il, i) 399 qmixmax(il) = qmix(sjmax(il)) 400 qmixmin(il) = qmix(sjmin(il)) 401 rmixmax(il) = rmix(sjmax(il)) 402 rmixmin(il) = rmix(sjmin(il)) 403 sqmrmax(il) = sjmax(il)*qmix(sjmax(il)) - rmix(sjmax(il)) 404 sqmrmin(il) = sjmin(il)*qmix(sjmin(il)) - rmix(sjmin(il)) 405 406 ment(il, i, j) = abs(qmixmax(il)-qmixmin(il))*ment(il, i, j) 407 408 ! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j) 409 IF (abs(qmixmax(il)-qmixmin(il))>1.E-10) THEN 410 sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/ & 411 (qmixmax(il)-qmixmin(il)) 402 Qmixmax(il) = Qmix(Sjmax(il)) 403 Qmixmin(il) = Qmix(Sjmin(il)) 404 Rmixmax(il) = Rmix(Sjmax(il)) 405 Rmixmin(il) = Rmix(Sjmin(il)) 406 sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il)) 407 sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il)) 408 409 Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j) 410 411 ! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j) 412 IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN 413 Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il)) 412 414 ELSE 413 sigij(il, i, j) = 0. 414 END IF 415 416 ! -- Compute Qent, uent, vent according to the true mixing 417 ! fraction 418 qent(il, i, j) = (1.-sigij(il,i,j))*rti + & 419 sigij(il, i, j)*rr(il, i) 420 uent(il, i, j) = (1.-sigij(il,i,j))*unk(il) + & 421 sigij(il, i, j)*u(il, i) 422 vent(il, i, j) = (1.-sigij(il,i,j))*vnk(il) + & 423 sigij(il, i, j)*v(il, i) 424 425 ! -- Compute liquid water static energy of mixed draughts 426 ! IF (j .GT. i) THEN 427 ! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j) 428 ! awat=amax1(awat,0.0) 429 ! ELSE 430 ! awat = 0. 431 ! ENDIF 432 ! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i) 433 ! : + Sigij(il,i,j)*H(il,i) 434 ! : + (LV(il,j)+(cpd-cpv)*t(il,j))*awat 435 ! IM 301008 beg 436 hent(il, i, j) = (1.-sigij(il,i,j))*hp(il, i) + & 437 sigij(il, i, j)*h(il, i) 438 439 elij(il, i, j) = qent(il, i, j) - rs(il, j) 440 elij(il, i, j) = elij(il, i, j) + ((h(il,j)-hent(il,i, & 441 j))*rs(il,j)*lv(il,j)/((cpd*(1.-qent(il,i,j))+ & 442 qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j))) 443 elij(il, i, j) = elij(il, i, j)/(1.+lv(il,j)*lv(il,j)*rs(il,j)/(( & 444 cpd*(1.-qent(il,i,j))+qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j))) 415 Sigij(il, i, j) = 0. 416 END IF 417 418 ! -- Compute Qent, uent, vent according to the true mixing fraction 419 Qent(il, i, j) = (1.-Sigij(il,i,j))*rti + Sigij(il, i, j)*rr(il, i) 420 uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i) 421 vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i) 422 423 ! -- Compute liquid water static energy of mixed draughts 424 ! IF (j .GT. i) THEN 425 ! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j) 426 ! awat=amax1(awat,0.0) 427 ! ELSE 428 ! awat = 0. 429 ! ENDIF 430 ! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i) 431 ! : + Sigij(il,i,j)*H(il,i) 432 ! : + (LV(il,j)+(cpd-cpv)*t(il,j))*awat 433 !IM 301008 beg 434 hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i) 435 436 elij(il, i, j) = Qent(il, i, j) - rs(il, j) 437 elij(il, i, j) = elij(il, i, j) + & 438 ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / & 439 ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j))) 440 elij(il, i, j) = elij(il, i, j) / & 441 (1.+lv(il,j)*lv(il,j)*rs(il,j) / & 442 ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j))) 445 443 446 444 elij(il, i, j) = max(elij(il,i,j), 0.) 447 445 448 elij(il, i, j) = min(elij(il,i,j), qent(il,i,j))446 elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j)) 449 447 450 448 IF (j>i) THEN … … 455 453 END IF 456 454 457 ! print 458 ! *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)* 459 ! : t(il,j)) 460 461 hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))* & 462 awat 463 ! IM 301008 end 464 465 ! print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ', 466 ! : i,j,hent(il,i,j),sigij(il,i,j) 467 468 ! -- ASij is the integral of P(F) over the relevant F 469 ! interval 470 asij(il) = asij(il) + abs(qmixmax(il)*(1.-sjmax(il))+rmixmax(il)- & 471 qmixmin(il)*(1.-sjmin(il))-rmixmin(il)) 455 ! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)* 456 ! : t(il,j)) 457 458 hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat 459 !IM 301008 end 460 461 ! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ', 462 ! : i,j,hent(il,i,j),Sigij(il,i,j) 463 464 ! -- ASij is the integral of P(F) over the relevant F interval 465 ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - & 466 Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il)) 472 467 473 468 END IF … … 476 471 DO k = 1, ntra 477 472 DO il = 1, ncum 478 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-&479 1)) .AND. (j<=inb(il)) .AND. lwork(il)) THEN480 IF (sij(il,i,j)>0.0) THEN481 traent(il, i, j, k) = sigij(il, i, j)*tra(il, i, k) + &482 (1.-sigij(il,i,j))*tra(il, nk(il), k)483 END IF484 END IF485 END DO486 END DO487 488 ! -- If I=J (detrainement and entrainement at the same level), then 489 !only the490 473 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. & 474 (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. & 475 lwork(il)) THEN 476 IF (Sij(il,i,j)>0.0) THEN 477 traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + & 478 (1.-Sigij(il,i,j))*tra(il, nk(il), k) 479 END IF 480 END IF 481 END DO 482 END DO 483 484 ! -- If I=J (detrainement and entrainement at the same level), then only the 485 ! -- adiabatic ascent part of the mixture is considered 491 486 IF (i==j) THEN 492 487 DO il = 1, ncum 493 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & 494 il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN 495 IF (sij(il,i,j)>0.0) THEN 488 IF (i>=icb(il) .AND. i<=inb(il) .AND. & 489 j>=(icb(il)-1) .AND. j<=inb(il) .AND. & 490 lwork(il)) THEN 491 IF (Sij(il,i,j)>0.0) THEN 496 492 rti = qnk(il) - ep(il, i)*clw(il, i) 497 ! cc Ment(il,i,i) = 498 ! m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il)) 499 ment(il, i, i) = abs(qmixmax(il)*(1.-sjmax( & 500 il))+rmixmax(il)-qmixmin(il)*(1.-sjmin(il))-rmixmin(il)) 501 qent(il, i, i) = rti 493 !!! Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il)) 494 Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - & 495 Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il)) 496 Qent(il, i, i) = rti 502 497 uent(il, i, i) = unk(il) 503 498 vent(il, i, i) = vnk(il) 504 499 hent(il, i, i) = hp(il, i) 505 500 elij(il, i, i) = clw(il, i)*(1.-ep(il,i)) 506 sigij(il, i, i) = 0.501 Sigij(il, i, i) = 0. 507 502 END IF 508 503 END IF … … 510 505 DO k = 1, ntra 511 506 DO il = 1, ncum 512 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- & 513 1)) .AND. (j<=inb(il)) .AND. lwork(il)) THEN 514 IF (sij(il,i,j)>0.0) THEN 507 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. & 508 (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. & 509 lwork(il)) THEN 510 IF (Sij(il,i,j)>0.0) THEN 515 511 traent(il, i, i, k) = tra(il, nk(il), k) 516 512 END IF … … 521 517 END IF 522 518 523 175 END DO 519 ! --------------------------------------------------------------- 520 175 END DO ! End loop on destination level "j" 521 ! --------------------------------------------------------------- 524 522 525 523 DO il = 1, ncum 526 524 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN 527 asij(il) = amax1(1.0E-16, asij(il))528 asij(il) = 1.0/asij(il)525 ASij(il) = amax1(1.0E-16, ASij(il)) 526 ASij(il) = 1.0/ASij(il) 529 527 csum(il, i) = 0.0 530 528 END IF … … 533 531 DO j = minorig, nl 534 532 DO il = 1, ncum 535 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(&536 il)-1) .AND. j<=inb(il)) THEN537 ment(il, i, j) = ment(il, i, j)*asij(il)533 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. & 534 j>=(icb(il)-1) .AND. j<=inb(il)) THEN 535 Ment(il, i, j) = Ment(il, i, j)*ASij(il) 538 536 END IF 539 537 END DO … … 542 540 DO j = minorig, nl 543 541 DO il = 1, ncum 544 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(&545 il)-1) .AND. j<=inb(il)) THEN546 csum(il, i) = csum(il, i) + ment(il, i, j)542 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. & 543 j>=(icb(il)-1) .AND. j<=inb(il)) THEN 544 csum(il, i) = csum(il, i) + Ment(il, i, j) 547 545 END IF 548 546 END DO … … 550 548 551 549 DO il = 1, ncum 552 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) & 553 THEN 554 ! cc : .and. csum(il,i).lt.m(il,i) ) then 550 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN 551 ! cc : .and. csum(il,i).lt.m(il,i) ) then 555 552 nent(il, i) = 0 556 ! cc ment(il,i,i)=m(il,i)557 ment(il, i, i) = 1.558 qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)553 ! cc Ment(il,i,i)=m(il,i) 554 Ment(il, i, i) = 1. 555 Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i) 559 556 uent(il, i, i) = unk(il) 560 557 vent(il, i, i) = vnk(il) 561 558 elij(il, i, i) = clw(il, i)*(1.-ep(il,i)) 562 sij(il, i, i) = 0.0559 Sij(il, i, i) = 0.0 563 560 END IF 564 561 END DO ! il … … 566 563 DO j = 1, ntra 567 564 DO il = 1, ncum 568 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) & 569 THEN 570 ! cc : .and. csum(il,i).lt.m(il,i) ) then 565 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN 566 ! cc : .and. csum(il,i).lt.m(il,i) ) then 571 567 traent(il, i, i, j) = tra(il, nk(il), j) 572 568 END IF … … 574 570 END DO 575 571 576 789 END DO 572 ! --------------------------------------------------------------- 573 789 END DO ! End loop on origin level "i" 574 ! --------------------------------------------------------------- 575 577 576 578 577 RETURN
Note: See TracChangeset
for help on using the changeset viewer.
