C C================================================================================ C SUBROUTINE CLOUDS_GNO(klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF) IMPLICIT NONE C C-------------------------------------------------------------------------------- C C Inputs: C C ND----------: Number of vertical levels C R--------ND-: Domain-averaged mixing ratio of total water C RS-------ND-: Mean saturation humidity mixing ratio within the gridbox C QSUB-----ND-: Mixing ratio of condensed water within clouds associated C with SUBGRID-SCALE condensation processes (here, it is C predicted by the convection scheme) C Outputs: C C PTCONV-----ND-: Point convectif = TRUE C RATQSC-----ND-: Largeur normalisee de la distribution C CLDF-----ND-: Fraction nuageuse C C-------------------------------------------------------------------------------- INTEGER klon,ND REAL R(klon,ND), RS(klon,ND), QSUB(klon,ND) LOGICAL PTCONV(klon,ND) REAL RATQSC(klon,ND) REAL CLDF(klon,ND) c -- parameters controlling the iteration: c -- nmax : maximum nb of iterations (hopefully never reached) c -- epsilon : accuracy of the numerical resolution c -- vmax : v-value above which we use an asymptotic expression for ERF(v) INTEGER nmax PARAMETER ( nmax = 10) REAL epsilon, vmax0, vmax(klon) PARAMETER ( epsilon = 0.02, vmax0 = 2.0 ) REAL min_mu, min_Q PARAMETER ( min_mu = 1.e-12, min_Q=1.e-12 ) INTEGER i,K, n, m REAL mu(klon), qsat(klon), delta(klon), beta(klon) real zu2(klon),zv2(klon) REAL xx(klon), aux(klon), coeff(klon), block(klon) REAL dist(klon), fprime(klon), det(klon) REAL pi, u(klon), v(klon), erfu(klon), erfv(klon) REAL xx1(klon), xx2(klon) real erf,kkk c lconv = true si le calcul a converge (entre autre si qsub < min_q) LOGICAL lconv(klon) pi = ACOS(-1.) ptconv=.false. ratqsc=0. DO 500 K = 1, ND do i=1,klon ! vector mu(i) = R(i,K) mu(i) = MAX(mu(i),min_mu) qsat(i) = RS(i,K) qsat(i) = MAX(qsat(i),min_mu) delta(i) = log(mu(i)/qsat(i)) enddo ! vector C C *** There is no subgrid-scale condensation; *** C *** the scheme becomes equivalent to an "all-or-nothing" *** C *** large-scale condensation scheme. *** C C C *** Some condensation is produced at the subgrid-scale *** C *** *** C *** PDF = generalized log-normal distribution (GNO) *** C *** (k<0 because a lower bound is considered for the PDF) *** C *** *** C *** -> Determine x (the parameter k of the GNO PDF) such *** C *** that the contribution of subgrid-scale processes to *** C *** the in-cloud water content is equal to QSUB(K) *** C *** (equations (13), (14), (15) + Appendix B of the paper) *** C *** *** C *** Here, an iterative method is used for this purpose *** C *** (other numerical methods might be more efficient) *** C *** *** C *** NB: the "error function" is called ERF *** C *** (ERF in double precision) *** C c On commence par eliminer les cas pour lesquels on n'a pas c suffisamment d'eau nuageuse. do i=1,klon ! vector IF ( QSUB(i,K) .lt. min_Q ) THEN ptconv(i,k)=.false. ratqsc(i,k)=0. lconv(i) = .true. c Rien on a deja initialise ELSE lconv(i) = .FALSE. vmax(i) = vmax0 beta(i) = QSUB(i,K)/mu(i) + EXP( -MIN(0.0,delta(i)) ) c -- roots of equation v > vmax: det(i) = delta(i) + vmax(i)**2. if (det(i).LE.0.0) vmax(i) = vmax0 + 1.0 det(i) = delta(i) + vmax(i)**2. if (det(i).LE.0.) then xx(i) = -0.0001 else xx1(i)=-SQRT(2.)*vmax(i)*(1.0-SQRT(1.0+delta(i)/(vmax(i)**2.))) xx2(i)=-SQRT(2.)*vmax(i)*(1.0+SQRT(1.0+delta(i)/(vmax(i)**2.))) xx(i) = 1.01 * xx1(i) if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i) endif if (delta(i).LT.0.) xx(i) = -0.5*SQRT(log(2.)) ENDIF enddo ! vector c---------------------------------------------------------------------- c Debut des nmax iterations pour trouver la solution. c---------------------------------------------------------------------- DO n = 1, nmax do i=1,klon ! vector if (.not.lconv(i)) then u(i) = delta(i)/(xx(i)*sqrt(2.)) + xx(i)/(2.*sqrt(2.)) v(i) = delta(i)/(xx(i)*sqrt(2.)) - xx(i)/(2.*sqrt(2.)) IF ( v(i) .GT. vmax(i) ) THEN IF ( ABS(u(i)) .GT. vmax(i) : .AND. delta(i) .LT. 0. ) THEN c -- use asymptotic expression of erf for u and v large: c ( -> analytic solution for xx ) aux(i) = 2.0*delta(i)*(1.-beta(i)*EXP(delta(i))) : /(1.+beta(i)*EXP(delta(i))) if (aux(i).lt.0.) then print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i) aux(i)=0. endif xx(i) = -SQRT(aux(i)) block(i) = EXP(-v(i)*v(i)) / v(i) / SQRT(pi) dist(i) = 0.0 fprime(i) = 1.0 ELSE c -- erfv -> 1.0, use an asymptotic expression of erfv for v large: erfu(i) = ERF(u(i)) c !!! ATTENTION : rajout d'un seuil pour l'exponentiel aux(i) = SQRT(pi)*(1.0-erfu(i))*EXP(min(v(i)*v(i),100.)) coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.) block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / SQRT(pi) dist(i) = v(i) * aux(i) / coeff(i) - beta(i) fprime(i) = 2.0 / xx(i) * (v(i)**2.) : * ( coeff(i)*EXP(-delta(i)) - u(i) * aux(i) ) : / coeff(i) / coeff(i) ENDIF ! ABS(u) ELSE c -- general case: erfu(i) = ERF(u(i)) erfv(i) = ERF(v(i)) block(i) = 1.0-erfv(i) dist(i) = (1.0 - erfu(i)) / (1.0 - erfv(i)) - beta(i) zu2(i)=u(i)*u(i) zv2(i)=v(i)*v(i) if(zu2(i).gt.20..or. zv2(i).gt.20.) then print*,'ATTENTION !!! xx(',i,') =', xx(i) print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF', .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k), .CLDF(i,k) print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i) zu2(i)=20. zv2(i)=20. fprime(i) = 0. else fprime(i) = 2. /SQRT(pi) /xx(i) /(1.0-erfv(i))**2. : * ( (1.0-erfv(i))*v(i)*EXP(-zu2(i)) : - (1.0-erfu(i))*u(i)*EXP(-zv2(i)) ) endif ENDIF ! x c -- test numerical convergence: c print*,'avant test ',i,k,lconv(i),u(i),v(i) if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then c print*,'v-u **2',(v(i)-u(i))**2 c print*,'exp v-u **2',exp((v(i)-u(i))**2) ptconv(i,K) = .TRUE. lconv(i)=.true. c borne pour l'exponentielle ratqsc(i,k)=min(2.*(v(i)-u(i))**2,20.) ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.) CLDF(i,K) = 0.5 * block(i) cccccccccccccccccccccccc c kkk=-sqrt(log(1.+ratqsc(i,k)**2)) c u(i)=delta(i)/(kkk*sqrt(2.))-kkk/(2.*sqrt(2.)) c v(i)=delta(i)/(kkk*sqrt(2.))+kkk/(2.*sqrt(2.)) c erfu(i)=erf(u(i)) c erfv(i)=erf(v(i)) c print*,'SIG ',k,qsub(i,k) c s ,mu(i)*((1.-erfv(i))/(1.-erfu(i))-qsat(i)/mu(i)) c s ,0.5*erfu(i) cccccccccccccccccccccccc else xx(i) = xx(i) - dist(i)/fprime(i) endif c print*,'apres test ',i,k,lconv(i) endif ! lconv enddo ! vector c---------------------------------------------------------------------- c Fin des nmax iterations pour trouver la solution. ENDDO ! n c---------------------------------------------------------------------- 500 CONTINUE ! K RETURN END