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