| 1 | |
|---|
| 2 | ! $Header$ |
|---|
| 3 | |
|---|
| 4 | |
|---|
| 5 | ! ================================================================================ |
|---|
| 6 | |
|---|
| 7 | SUBROUTINE clouds_gno(klon, nd, r, rs, qsub, ptconv, ratqsc, cldf) |
|---|
| 8 | IMPLICIT NONE |
|---|
| 9 | |
|---|
| 10 | ! -------------------------------------------------------------------------------- |
|---|
| 11 | |
|---|
| 12 | ! Inputs: |
|---|
| 13 | |
|---|
| 14 | ! ND----------: Number of vertical levels |
|---|
| 15 | ! R--------ND-: Domain-averaged mixing ratio of total water |
|---|
| 16 | ! RS-------ND-: Mean saturation humidity mixing ratio within the gridbox |
|---|
| 17 | ! QSUB-----ND-: Mixing ratio of condensed water within clouds associated |
|---|
| 18 | ! with SUBGRID-SCALE condensation processes (here, it is |
|---|
| 19 | ! predicted by the convection scheme) |
|---|
| 20 | ! Outputs: |
|---|
| 21 | |
|---|
| 22 | ! PTCONV-----ND-: Point convectif = TRUE |
|---|
| 23 | ! RATQSC-----ND-: Largeur normalisee de la distribution |
|---|
| 24 | ! CLDF-----ND-: Fraction nuageuse |
|---|
| 25 | |
|---|
| 26 | ! -------------------------------------------------------------------------------- |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | INTEGER klon, nd |
|---|
| 30 | REAL r(klon, nd), rs(klon, nd), qsub(klon, nd) |
|---|
| 31 | LOGICAL ptconv(klon, nd) |
|---|
| 32 | REAL ratqsc(klon, nd) |
|---|
| 33 | REAL cldf(klon, nd) |
|---|
| 34 | |
|---|
| 35 | ! -- parameters controlling the iteration: |
|---|
| 36 | ! -- nmax : maximum nb of iterations (hopefully never reached) |
|---|
| 37 | ! -- epsilon : accuracy of the numerical resolution |
|---|
| 38 | ! -- vmax : v-value above which we use an asymptotic expression for |
|---|
| 39 | ! ERF(v) |
|---|
| 40 | |
|---|
| 41 | INTEGER nmax |
|---|
| 42 | PARAMETER (nmax=10) |
|---|
| 43 | REAL epsilon, vmax0, vmax(klon) |
|---|
| 44 | PARAMETER (epsilon=0.02, vmax0=2.0) |
|---|
| 45 | |
|---|
| 46 | REAL min_mu, min_q |
|---|
| 47 | PARAMETER (min_mu=1.E-12, min_q=1.E-12) |
|---|
| 48 | |
|---|
| 49 | INTEGER i, k, n, m |
|---|
| 50 | REAL mu(klon), qsat, delta(klon), beta(klon) |
|---|
| 51 | REAL zu2, zv2 |
|---|
| 52 | REAL xx(klon), aux(klon), coeff, block |
|---|
| 53 | REAL dist, fprime, det |
|---|
| 54 | REAL pi, u, v, erfcu, erfcv |
|---|
| 55 | REAL xx1, xx2 |
|---|
| 56 | REAL erf, hsqrtlog_2, v2 |
|---|
| 57 | REAL sqrtpi, sqrt2, zx1, zx2, exdel |
|---|
| 58 | ! lconv = true si le calcul a converge (entre autre si qsub < min_q) |
|---|
| 59 | LOGICAL lconv(klon) |
|---|
| 60 | |
|---|
| 61 | ! cdir arraycomb |
|---|
| 62 | cldf(1:klon, 1:nd) = 0.0 ! cym |
|---|
| 63 | ratqsc(1:klon, 1:nd) = 0.0 |
|---|
| 64 | ptconv(1:klon, 1:nd) = .FALSE. |
|---|
| 65 | ! cdir end arraycomb |
|---|
| 66 | |
|---|
| 67 | pi = acos(-1.) |
|---|
| 68 | sqrtpi = sqrt(pi) |
|---|
| 69 | sqrt2 = sqrt(2.) |
|---|
| 70 | hsqrtlog_2 = 0.5*sqrt(log(2.)) |
|---|
| 71 | |
|---|
| 72 | DO k = 1, nd |
|---|
| 73 | |
|---|
| 74 | DO i = 1, klon ! vector |
|---|
| 75 | mu(i) = r(i, k) |
|---|
| 76 | mu(i) = max(mu(i), min_mu) |
|---|
| 77 | qsat = rs(i, k) |
|---|
| 78 | qsat = max(qsat, min_mu) |
|---|
| 79 | delta(i) = log(mu(i)/qsat) |
|---|
| 80 | ! enddo ! vector |
|---|
| 81 | |
|---|
| 82 | |
|---|
| 83 | ! *** There is no subgrid-scale condensation; *** |
|---|
| 84 | ! *** the scheme becomes equivalent to an "all-or-nothing" *** |
|---|
| 85 | ! *** large-scale condensation scheme. *** |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | |
|---|
| 89 | ! *** Some condensation is produced at the subgrid-scale *** |
|---|
| 90 | ! *** *** |
|---|
| 91 | ! *** PDF = generalized log-normal distribution (GNO) *** |
|---|
| 92 | ! *** (k<0 because a lower bound is considered for the PDF) *** |
|---|
| 93 | ! *** *** |
|---|
| 94 | ! *** -> Determine x (the parameter k of the GNO PDF) such *** |
|---|
| 95 | ! *** that the contribution of subgrid-scale processes to *** |
|---|
| 96 | ! *** the in-cloud water content is equal to QSUB(K) *** |
|---|
| 97 | ! *** (equations (13), (14), (15) + Appendix B of the paper) *** |
|---|
| 98 | ! *** *** |
|---|
| 99 | ! *** Here, an iterative method is used for this purpose *** |
|---|
| 100 | ! *** (other numerical methods might be more efficient) *** |
|---|
| 101 | ! *** *** |
|---|
| 102 | ! *** NB: the "error function" is called ERF *** |
|---|
| 103 | ! *** (ERF in double precision) *** |
|---|
| 104 | |
|---|
| 105 | |
|---|
| 106 | ! On commence par eliminer les cas pour lesquels on n'a pas |
|---|
| 107 | ! suffisamment d'eau nuageuse. |
|---|
| 108 | |
|---|
| 109 | ! do i=1,klon ! vector |
|---|
| 110 | |
|---|
| 111 | IF (qsub(i,k)<min_q) THEN |
|---|
| 112 | ptconv(i, k) = .FALSE. |
|---|
| 113 | ratqsc(i, k) = 0. |
|---|
| 114 | lconv(i) = .TRUE. |
|---|
| 115 | |
|---|
| 116 | ! Rien on a deja initialise |
|---|
| 117 | |
|---|
| 118 | ELSE |
|---|
| 119 | |
|---|
| 120 | lconv(i) = .FALSE. |
|---|
| 121 | vmax(i) = vmax0 |
|---|
| 122 | |
|---|
| 123 | beta(i) = qsub(i, k)/mu(i) + exp(-min(0.0,delta(i))) |
|---|
| 124 | |
|---|
| 125 | ! -- roots of equation v > vmax: |
|---|
| 126 | |
|---|
| 127 | det = delta(i) + vmax(i)*vmax(i) |
|---|
| 128 | IF (det<=0.0) vmax(i) = vmax0 + 1.0 |
|---|
| 129 | det = delta(i) + vmax(i)*vmax(i) |
|---|
| 130 | |
|---|
| 131 | IF (det<=0.) THEN |
|---|
| 132 | xx(i) = -0.0001 |
|---|
| 133 | ELSE |
|---|
| 134 | zx1 = -sqrt2*vmax(i) |
|---|
| 135 | zx2 = sqrt(1.0+delta(i)/(vmax(i)*vmax(i))) |
|---|
| 136 | xx1 = zx1*(1.0-zx2) |
|---|
| 137 | xx2 = zx1*(1.0+zx2) |
|---|
| 138 | xx(i) = 1.01*xx1 |
|---|
| 139 | IF (xx1>=0.0) xx(i) = 0.5*xx2 |
|---|
| 140 | END IF |
|---|
| 141 | IF (delta(i)<0.) xx(i) = -hsqrtlog_2 |
|---|
| 142 | |
|---|
| 143 | END IF |
|---|
| 144 | |
|---|
| 145 | END DO ! vector |
|---|
| 146 | |
|---|
| 147 | ! ---------------------------------------------------------------------- |
|---|
| 148 | ! Debut des nmax iterations pour trouver la solution. |
|---|
| 149 | ! ---------------------------------------------------------------------- |
|---|
| 150 | |
|---|
| 151 | DO n = 1, nmax |
|---|
| 152 | |
|---|
| 153 | DO i = 1, klon ! vector |
|---|
| 154 | IF (.NOT. lconv(i)) THEN |
|---|
| 155 | |
|---|
| 156 | u = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2) |
|---|
| 157 | v = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2) |
|---|
| 158 | v2 = v*v |
|---|
| 159 | |
|---|
| 160 | IF (v>vmax(i)) THEN |
|---|
| 161 | |
|---|
| 162 | IF (abs(u)>vmax(i) .AND. delta(i)<0.) THEN |
|---|
| 163 | |
|---|
| 164 | ! -- use asymptotic expression of erf for u and v large: |
|---|
| 165 | ! ( -> analytic solution for xx ) |
|---|
| 166 | exdel = beta(i)*exp(delta(i)) |
|---|
| 167 | aux(i) = 2.0*delta(i)*(1.-exdel)/(1.+exdel) |
|---|
| 168 | IF (aux(i)<0.) THEN |
|---|
| 169 | ! print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i) |
|---|
| 170 | aux(i) = 0. |
|---|
| 171 | END IF |
|---|
| 172 | xx(i) = -sqrt(aux(i)) |
|---|
| 173 | block = exp(-v*v)/v/sqrtpi |
|---|
| 174 | dist = 0.0 |
|---|
| 175 | fprime = 1.0 |
|---|
| 176 | |
|---|
| 177 | ELSE |
|---|
| 178 | |
|---|
| 179 | ! -- erfv -> 1.0, use an asymptotic expression of erfv for v |
|---|
| 180 | ! large: |
|---|
| 181 | |
|---|
| 182 | erfcu = 1.0 - erf(u) |
|---|
| 183 | ! !!! ATTENTION : rajout d'un seuil pour l'exponentiel |
|---|
| 184 | aux(i) = sqrtpi*erfcu*exp(min(v2,100.)) |
|---|
| 185 | coeff = 1.0 - 0.5/(v2) + 0.75/(v2*v2) |
|---|
| 186 | block = coeff*exp(-v2)/v/sqrtpi |
|---|
| 187 | dist = v*aux(i)/coeff - beta(i) |
|---|
| 188 | fprime = 2.0/xx(i)*(v2)*(exp(-delta(i))-u*aux(i)/coeff)/coeff |
|---|
| 189 | |
|---|
| 190 | END IF ! ABS(u) |
|---|
| 191 | |
|---|
| 192 | ELSE |
|---|
| 193 | |
|---|
| 194 | ! -- general case: |
|---|
| 195 | |
|---|
| 196 | erfcu = 1.0 - erf(u) |
|---|
| 197 | erfcv = 1.0 - erf(v) |
|---|
| 198 | block = erfcv |
|---|
| 199 | dist = erfcu/erfcv - beta(i) |
|---|
| 200 | zu2 = u*u |
|---|
| 201 | zv2 = v2 |
|---|
| 202 | IF (zu2>20. .OR. zv2>20.) THEN |
|---|
| 203 | ! print*,'ATTENTION !!! xx(',i,') =', xx(i) |
|---|
| 204 | ! print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF', |
|---|
| 205 | ! .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k), |
|---|
| 206 | ! .CLDF(i,k) |
|---|
| 207 | ! print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i) |
|---|
| 208 | zu2 = 20. |
|---|
| 209 | zv2 = 20. |
|---|
| 210 | fprime = 0. |
|---|
| 211 | ELSE |
|---|
| 212 | fprime = 2./sqrtpi/xx(i)/(erfcv*erfcv)* & |
|---|
| 213 | (erfcv*v*exp(-zu2)-erfcu*u*exp(-zv2)) |
|---|
| 214 | END IF |
|---|
| 215 | END IF ! x |
|---|
| 216 | |
|---|
| 217 | ! -- test numerical convergence: |
|---|
| 218 | |
|---|
| 219 | ! if (beta(i).lt.1.e-10) then |
|---|
| 220 | ! print*,'avant test ',i,k,lconv(i),u(i),v(i),beta(i) |
|---|
| 221 | ! stop |
|---|
| 222 | ! endif |
|---|
| 223 | IF (abs(fprime)<1.E-11) THEN |
|---|
| 224 | ! print*,'avant test fprime<.e-11 ' |
|---|
| 225 | ! s ,i,k,lconv(i),u(i),v(i),beta(i),fprime(i) |
|---|
| 226 | ! print*,'klon,ND,R,RS,QSUB', |
|---|
| 227 | ! s klon,ND,R(i,k),rs(i,k),qsub(i,k) |
|---|
| 228 | fprime = sign(1.E-11, fprime) |
|---|
| 229 | END IF |
|---|
| 230 | |
|---|
| 231 | |
|---|
| 232 | IF (abs(dist/beta(i))<epsilon) THEN |
|---|
| 233 | ! print*,'v-u **2',(v(i)-u(i))**2 |
|---|
| 234 | ! print*,'exp v-u **2',exp((v(i)-u(i))**2) |
|---|
| 235 | ptconv(i, k) = .TRUE. |
|---|
| 236 | lconv(i) = .TRUE. |
|---|
| 237 | ! borne pour l'exponentielle |
|---|
| 238 | ratqsc(i, k) = min(2.*(v-u)*(v-u), 20.) |
|---|
| 239 | ratqsc(i, k) = sqrt(exp(ratqsc(i,k))-1.) |
|---|
| 240 | cldf(i, k) = 0.5*block |
|---|
| 241 | ELSE |
|---|
| 242 | xx(i) = xx(i) - dist/fprime |
|---|
| 243 | END IF |
|---|
| 244 | ! print*,'apres test ',i,k,lconv(i) |
|---|
| 245 | |
|---|
| 246 | END IF ! lconv |
|---|
| 247 | END DO ! vector |
|---|
| 248 | |
|---|
| 249 | ! ---------------------------------------------------------------------- |
|---|
| 250 | ! Fin des nmax iterations pour trouver la solution. |
|---|
| 251 | END DO ! n |
|---|
| 252 | ! ---------------------------------------------------------------------- |
|---|
| 253 | |
|---|
| 254 | |
|---|
| 255 | END DO |
|---|
| 256 | ! K |
|---|
| 257 | RETURN |
|---|
| 258 | END SUBROUTINE clouds_gno |
|---|
| 259 | |
|---|
| 260 | |
|---|
| 261 | |
|---|
