trunk/phylmd/clouds_gno.f90

module CLOUDS_GNO_m

  IMPLICIT NONE

contains

  SUBROUTINE CLOUDS_GNO(klon, ND, R, RS, QSUB, PTCONV, RATQSC, CLDF)

    ! From LMDZ4/libf/phylmd/clouds_gno.F, version 1.2, 2004/11/09 16:55:40

    use numer_rec_95, only: nr_erf

    INTEGER, intent(in):: klon
    INTEGER, intent(in):: ND ! number of vertical levels

    REAL, intent(in):: R(klon, ND)
    ! domain-averaged mixing ratio of total water 

    REAL, intent(in):: RS(klon, ND)
    ! mean saturation humidity mixing ratio within the gridbox

    REAL, intent(in):: QSUB(klon, ND)
    ! mixing ratio of condensed water within clouds associated
    ! with SUBGRID-SCALE condensation processes (here, it is
    ! predicted by the convection scheme)

    LOGICAL, intent(out):: PTCONV(klon, ND) ! Point convectif = TRUE
    REAL, intent(out):: RATQSC(klon, ND) ! largeur normalisee de la distribution
    REAL, intent(out):: CLDF(klon, ND) ! fraction nuageuse

    ! Local:
    
    ! parameters controlling the iteration:
    ! nmax : maximum nb of iterations (hopefully never reached)
    ! epsilon : accuracy of the numerical resolution 
    ! 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
    REAL mu(klon), qsat(klon), delta(klon), beta(klon) 
    real zu2(klon), zv2(klon)
    REAL xx(klon), aux(klon), coeff(klon), my_block(klon)
    REAL dist(klon), fprime(klon), det(klon)
    REAL pi, u(klon), v(klon), erfcu(klon), erfcv(klon)
    REAL xx1(klon), xx2(klon)
    real sqrtpi, sqrt2, zx1, zx2, exdel
    ! lconv = true si le calcul a converge (entre autres si qsub < min_q)
    LOGICAL lconv(klon)

    !--------------------------------------------------------------

    cldf=0.0

    pi = ACOS(-1.)
    sqrtpi=sqrt(pi)
    sqrt2=sqrt(2.)

    ptconv=.false.
    ratqsc=0.

    loop_vertical: DO K = 1, ND
       do i=1, klon
          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

       ! There is no subgrid-scale condensation; the scheme becomes
       ! equivalent to an "all-or-nothing" large-scale condensation
       ! scheme.

       ! Some condensation is produced at the subgrid-scale 
       ! 
       ! PDF = generalized log-normal distribution (GNO) 
       ! (k<0 because a lower bound is considered for the PDF) 
       ! 
       ! -> Determine x (the parameter k of the GNO PDF) such that the
       ! contribution of subgrid-scale processes to the in-cloud water
       ! content is equal to QSUB(K) (equations (13), (14), (15) +
       ! Appendix B of the paper)
       ! 
       ! Here, an iterative method is used for this purpose (other
       ! numerical methods might be more efficient)
       ! 
       ! NB: the "error function" is called ERF (ERF in double
       ! precision)

       ! On commence par eliminer les cas pour lesquels on n'a pas
       ! suffisamment d'eau nuageuse.

       do i=1, klon
          IF ( QSUB(i, K) .lt. min_Q ) THEN
             ptconv(i, k)=.false.
             ratqsc(i, k)=0.
             lconv(i) = .true.
          ELSE 
             lconv(i) = .FALSE. 
             vmax(i) = vmax0

             beta(i) = QSUB(i, K)/mu(i) + EXP( -MIN(0.0, delta(i)) )

             ! 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 
                zx1=-sqrt2*vmax(i)
                zx2=SQRT(1.0+delta(i)/(vmax(i)**2.))
                xx1(i)=zx1*(1.0-zx2)
                xx2(i)=zx1*(1.0+zx2)
                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

       ! Debut des nmax iterations pour trouver la solution.
       loop_n: DO n = 1, nmax 
          loop_horizontal: do i = 1, klon
             test_lconv: if (.not.lconv(i)) then
                u(i) = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2)
                v(i) = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2)

                IF ( v(i) .GT. vmax(i) ) THEN 
                   IF ( ABS(u(i)) .GT. vmax(i) .AND. delta(i) .LT. 0. ) THEN
                      ! use asymptotic expression of erf for u and v large:
                      ! ( -> analytic solution for xx )
                      exdel=beta(i)*EXP(delta(i))
                      aux(i) = 2.0*delta(i)*(1.-exdel) /(1.+exdel)
                      if (aux(i).lt.0.) then
                         aux(i)=0.
                      endif
                      xx(i) = -SQRT(aux(i))
                      my_block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi
                      dist(i) = 0.0
                      fprime(i) = 1.0
                   ELSE
                      ! erfv -> 1.0, use an asymptotic expression of
                      ! erfv for v large:

                      erfcu(i) = 1.0-NR_ERF(u(i))
                      ! Attention : ajout d'un seuil pour l'exponentielle
                      aux(i) = sqrtpi*erfcu(i)*EXP(min(v(i)*v(i), 80.))
                      coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.)
                      my_block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi
                      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
                ELSE
                   ! general case:

                   erfcu(i) = 1.0-NR_ERF(u(i))
                   erfcv(i) = 1.0-NR_ERF(v(i))
                   my_block(i) = erfcv(i)
                   dist(i) = erfcu(i) / erfcv(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
                      zu2(i)=20.
                      zv2(i)=20.
                      fprime(i) = 0.
                   else
                      fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2. &
                           * ( erfcv(i)*v(i)*EXP(-zu2(i)) &
                           - erfcu(i)*u(i)*EXP(-zv2(i)) )
                   endif
                ENDIF

                ! test numerical convergence:
                if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then 
                   ptconv(i, K) = .TRUE. 
                   lconv(i)=.true.
                   ! 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 * my_block(i)
                else
                   xx(i) = xx(i) - dist(i)/fprime(i)
                endif
             endif test_lconv
          enddo loop_horizontal
       ENDDO loop_n
    end DO loop_vertical

  END SUBROUTINE CLOUDS_GNO

end module CLOUDS_GNO_m
1	guez	72	module CLOUDS_GNO_m
2	guez	3
3	guez	72	IMPLICIT NONE
4	guez	3
5	guez	72	contains
6	guez	3
7	guez	72	SUBROUTINE CLOUDS_GNO(klon, ND, R, RS, QSUB, PTCONV, RATQSC, CLDF)
8	guez	3
9	guez	97	! From LMDZ4/libf/phylmd/clouds_gno.F, version 1.2, 2004/11/09 16:55:40
10	guez	3
11	guez	72	use numer_rec_95, only: nr_erf
12	guez	3
13	guez	346	INTEGER, intent(in):: klon
14			INTEGER, intent(in):: ND ! number of vertical levels
15	guez	3
16	guez	346	REAL, intent(in):: R(klon, ND)
17			! domain-averaged mixing ratio of total water
18	guez	3
19	guez	346	REAL, intent(in):: RS(klon, ND)
20			! mean saturation humidity mixing ratio within the gridbox
21
22			REAL, intent(in):: QSUB(klon, ND)
23			! mixing ratio of condensed water within clouds associated
24	guez	72	! with SUBGRID-SCALE condensation processes (here, it is
25			! predicted by the convection scheme)
26	guez	3
27	guez	346	LOGICAL, intent(out):: PTCONV(klon, ND) ! Point convectif = TRUE
28			REAL, intent(out):: RATQSC(klon, ND) ! largeur normalisee de la distribution
29			REAL, intent(out):: CLDF(klon, ND) ! fraction nuageuse
30	guez	3
31	guez	346	! Local:
32
33	guez	72	! parameters controlling the iteration:
34			! nmax : maximum nb of iterations (hopefully never reached)
35			! epsilon : accuracy of the numerical resolution
36			! vmax : v-value above which we use an asymptotic expression for ERF(v)
37	guez	3
38	guez	72	INTEGER nmax
39			PARAMETER ( nmax = 10)
40			REAL epsilon, vmax0, vmax(klon)
41			PARAMETER ( epsilon = 0.02, vmax0 = 2.0 )
42	guez	3
43	guez	72	REAL min_mu, min_Q
44			PARAMETER ( min_mu = 1.e-12, min_Q=1.e-12 )
45	guez	3
46	guez	72	INTEGER i, K, n
47			REAL mu(klon), qsat(klon), delta(klon), beta(klon)
48			real zu2(klon), zv2(klon)
49	guez	346	REAL xx(klon), aux(klon), coeff(klon), my_block(klon)
50	guez	72	REAL dist(klon), fprime(klon), det(klon)
51			REAL pi, u(klon), v(klon), erfcu(klon), erfcv(klon)
52			REAL xx1(klon), xx2(klon)
53			real sqrtpi, sqrt2, zx1, zx2, exdel
54			! lconv = true si le calcul a converge (entre autres si qsub < min_q)
55			LOGICAL lconv(klon)
56	guez	3
57	guez	72	!--------------------------------------------------------------
58	guez	3
59	guez	346	cldf=0.0
60	guez	3
61	guez	72	pi = ACOS(-1.)
62			sqrtpi=sqrt(pi)
63			sqrt2=sqrt(2.)
64	guez	3
65	guez	72	ptconv=.false.
66			ratqsc=0.
67	guez	3
68	guez	72	loop_vertical: DO K = 1, ND
69			do i=1, klon
70			mu(i) = R(i, K)
71			mu(i) = MAX(mu(i), min_mu)
72			qsat(i) = RS(i, K)
73			qsat(i) = MAX(qsat(i), min_mu)
74			delta(i) = log(mu(i)/qsat(i))
75			enddo
76	guez	3
77	guez	72	! There is no subgrid-scale condensation; the scheme becomes
78			! equivalent to an "all-or-nothing" large-scale condensation
79			! scheme.
80	guez	3
81	guez	72	! Some condensation is produced at the subgrid-scale
82			!
83			! PDF = generalized log-normal distribution (GNO)
84			! (k<0 because a lower bound is considered for the PDF)
85			!
86			! -> Determine x (the parameter k of the GNO PDF) such that the
87			! contribution of subgrid-scale processes to the in-cloud water
88			! content is equal to QSUB(K) (equations (13), (14), (15) +
89			! Appendix B of the paper)
90			!
91			! Here, an iterative method is used for this purpose (other
92			! numerical methods might be more efficient)
93			!
94			! NB: the "error function" is called ERF (ERF in double
95			! precision)
96	guez	3
97	guez	72	! On commence par eliminer les cas pour lesquels on n'a pas
98			! suffisamment d'eau nuageuse.
99	guez	3
100	guez	72	do i=1, klon
101			IF ( QSUB(i, K) .lt. min_Q ) THEN
102			ptconv(i, k)=.false.
103			ratqsc(i, k)=0.
104			lconv(i) = .true.
105			ELSE
106			lconv(i) = .FALSE.
107			vmax(i) = vmax0
108	guez	3
109	guez	72	beta(i) = QSUB(i, K)/mu(i) + EXP( -MIN(0.0, delta(i)) )
110	guez	3
111	guez	72	! roots of equation v > vmax:
112	guez	3
113	guez	72	det(i) = delta(i) + vmax(i)**2.
114			if (det(i).LE.0.0) vmax(i) = vmax0 + 1.0
115			det(i) = delta(i) + vmax(i)**2.
116	guez	3
117	guez	72	if (det(i).LE.0.) then
118			xx(i) = -0.0001
119			else
120			zx1=-sqrt2*vmax(i)
121			zx2=SQRT(1.0+delta(i)/(vmax(i)**2.))
122			xx1(i)=zx1*(1.0-zx2)
123			xx2(i)=zx1*(1.0+zx2)
124			xx(i) = 1.01 * xx1(i)
125			if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i)
126	guez	3	endif
127	guez	72	if (delta(i).LT.0.) xx(i) = -0.5*SQRT(log(2.))
128			ENDIF
129			enddo
130	guez	3
131	guez	72	! Debut des nmax iterations pour trouver la solution.
132	guez	346	loop_n: DO n = 1, nmax
133	guez	72	loop_horizontal: do i = 1, klon
134			test_lconv: if (.not.lconv(i)) then
135			u(i) = delta(i)/(xx(i)sqrt2) + xx(i)/(2.sqrt2)
136			v(i) = delta(i)/(xx(i)sqrt2) - xx(i)/(2.sqrt2)
137	guez	3
138	guez	72	IF ( v(i) .GT. vmax(i) ) THEN
139			IF ( ABS(u(i)) .GT. vmax(i) .AND. delta(i) .LT. 0. ) THEN
140			! use asymptotic expression of erf for u and v large:
141			! ( -> analytic solution for xx )
142			exdel=beta(i)*EXP(delta(i))
143			aux(i) = 2.0delta(i)(1.-exdel) /(1.+exdel)
144			if (aux(i).lt.0.) then
145			aux(i)=0.
146			endif
147			xx(i) = -SQRT(aux(i))
148	guez	346	my_block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi
149	guez	72	dist(i) = 0.0
150			fprime(i) = 1.0
151			ELSE
152			! erfv -> 1.0, use an asymptotic expression of
153			! erfv for v large:
154	guez	3
155	guez	72	erfcu(i) = 1.0-NR_ERF(u(i))
156	guez	227	! Attention : ajout d'un seuil pour l'exponentielle
157			aux(i) = sqrtpierfcu(i)EXP(min(v(i)*v(i), 80.))
158	guez	72	coeff(i) = 1.0 - 1./2./(v(i)2.) + 3./4./(v(i)4.)
159	guez	346	my_block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi
160	guez	72	dist(i) = v(i) * aux(i) / coeff(i) - beta(i)
161			fprime(i) = 2.0 / xx(i) * (v(i)**2.) &
162			* ( coeff(i)EXP(-delta(i)) - u(i) aux(i) ) &
163			/ coeff(i) / coeff(i)
164			ENDIF
165			ELSE
166			! general case:
167	guez	3
168	guez	72	erfcu(i) = 1.0-NR_ERF(u(i))
169			erfcv(i) = 1.0-NR_ERF(v(i))
170	guez	346	my_block(i) = erfcv(i)
171	guez	72	dist(i) = erfcu(i) / erfcv(i) - beta(i)
172			zu2(i)=u(i)*u(i)
173			zv2(i)=v(i)*v(i)
174			if(zu2(i).gt.20..or. zv2(i).gt.20.) then
175			zu2(i)=20.
176			zv2(i)=20.
177			fprime(i) = 0.
178			else
179			fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2. &
180			* ( erfcv(i)v(i)EXP(-zu2(i)) &
181			- erfcu(i)u(i)EXP(-zv2(i)) )
182			endif
183			ENDIF
184	guez	3
185	guez	72	! test numerical convergence:
186			if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then
187			ptconv(i, K) = .TRUE.
188			lconv(i)=.true.
189			! borne pour l'exponentielle
190			ratqsc(i, k)=min(2.(v(i)-u(i))*2, 20.)
191			ratqsc(i, k)=sqrt(exp(ratqsc(i, k))-1.)
192	guez	346	CLDF(i, K) = 0.5 * my_block(i)
193	guez	72	else
194			xx(i) = xx(i) - dist(i)/fprime(i)
195			endif
196			endif test_lconv
197			enddo loop_horizontal
198	guez	346	ENDDO loop_n
199	guez	72	end DO loop_vertical
200	guez	3
201	guez	72	END SUBROUTINE CLOUDS_GNO
202	guez	3
203	guez	72	end module CLOUDS_GNO_m