Sources/phylmd/clouds_gno.f

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

    ! Inputs:

    ! ND : Number of vertical levels
    ! R ND: Domain-averaged mixing ratio of total water 
    ! RS ND: Mean saturation humidity mixing ratio within the gridbox

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

    ! Outputs:

    ! PTCONV ND: Point convectif = TRUE
    ! RATQSC ND: Largeur normalisee de la distribution
    ! CLDF ND: Fraction nuageuse

    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)

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