1 |
module CLOUDS_GNO_m |
2 |
|
3 |
IMPLICIT NONE |
4 |
|
5 |
contains |
6 |
|
7 |
SUBROUTINE CLOUDS_GNO(klon, ND, R, RS, QSUB, PTCONV, RATQSC, CLDF) |
8 |
|
9 |
! From LMDZ4/libf/phylmd/clouds_gno.F, version 1.2, 2004/11/09 16:55:40 |
10 |
|
11 |
use numer_rec_95, only: nr_erf |
12 |
|
13 |
INTEGER, intent(in):: klon |
14 |
INTEGER, intent(in):: ND ! number of vertical levels |
15 |
|
16 |
REAL, intent(in):: R(klon, ND) |
17 |
! domain-averaged mixing ratio of total water |
18 |
|
19 |
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 |
! with SUBGRID-SCALE condensation processes (here, it is |
25 |
! predicted by the convection scheme) |
26 |
|
27 |
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 |
|
31 |
! Local: |
32 |
|
33 |
! 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 |
|
38 |
INTEGER nmax |
39 |
PARAMETER ( nmax = 10) |
40 |
REAL epsilon, vmax0, vmax(klon) |
41 |
PARAMETER ( epsilon = 0.02, vmax0 = 2.0 ) |
42 |
|
43 |
REAL min_mu, min_Q |
44 |
PARAMETER ( min_mu = 1.e-12, min_Q=1.e-12 ) |
45 |
|
46 |
INTEGER i, K, n |
47 |
REAL mu(klon), qsat(klon), delta(klon), beta(klon) |
48 |
real zu2(klon), zv2(klon) |
49 |
REAL xx(klon), aux(klon), coeff(klon), my_block(klon) |
50 |
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 |
|
57 |
!-------------------------------------------------------------- |
58 |
|
59 |
cldf=0.0 |
60 |
|
61 |
pi = ACOS(-1.) |
62 |
sqrtpi=sqrt(pi) |
63 |
sqrt2=sqrt(2.) |
64 |
|
65 |
ptconv=.false. |
66 |
ratqsc=0. |
67 |
|
68 |
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 |
|
77 |
! There is no subgrid-scale condensation; the scheme becomes |
78 |
! equivalent to an "all-or-nothing" large-scale condensation |
79 |
! scheme. |
80 |
|
81 |
! 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 |
|
97 |
! On commence par eliminer les cas pour lesquels on n'a pas |
98 |
! suffisamment d'eau nuageuse. |
99 |
|
100 |
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 |
|
109 |
beta(i) = QSUB(i, K)/mu(i) + EXP( -MIN(0.0, delta(i)) ) |
110 |
|
111 |
! roots of equation v > vmax: |
112 |
|
113 |
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 |
|
117 |
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 |
endif |
127 |
if (delta(i).LT.0.) xx(i) = -0.5*SQRT(log(2.)) |
128 |
ENDIF |
129 |
enddo |
130 |
|
131 |
! Debut des nmax iterations pour trouver la solution. |
132 |
loop_n: DO n = 1, nmax |
133 |
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 |
|
138 |
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.0*delta(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 |
my_block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi |
149 |
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 |
|
155 |
erfcu(i) = 1.0-NR_ERF(u(i)) |
156 |
! Attention : ajout d'un seuil pour l'exponentielle |
157 |
aux(i) = sqrtpi*erfcu(i)*EXP(min(v(i)*v(i), 80.)) |
158 |
coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.) |
159 |
my_block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi |
160 |
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 |
|
168 |
erfcu(i) = 1.0-NR_ERF(u(i)) |
169 |
erfcv(i) = 1.0-NR_ERF(v(i)) |
170 |
my_block(i) = erfcv(i) |
171 |
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 |
|
185 |
! 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 |
CLDF(i, K) = 0.5 * my_block(i) |
193 |
else |
194 |
xx(i) = xx(i) - dist(i)/fprime(i) |
195 |
endif |
196 |
endif test_lconv |
197 |
enddo loop_horizontal |
198 |
ENDDO loop_n |
199 |
end DO loop_vertical |
200 |
|
201 |
END SUBROUTINE CLOUDS_GNO |
202 |
|
203 |
end module CLOUDS_GNO_m |