Context Navigation

newmicro.F90 @ 5872

Last change on this file since 5872 was 3281, checked in by musat, 8 years ago
Correction commission precedente
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.3 KB

Line
1	! $Id: newmicro.F90 3281 2018-03-16 18:26:14Z evignon $
2
3	SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
4	pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
5	mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, &
6	reliq_pi, reice_pi)
7
8	USE dimphy
9	USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
10	reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, &
11	zfice, dNovrN
12	USE phys_state_var_mod, ONLY: rnebcon, clwcon
13	USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
14	USE ioipsl_getin_p_mod, ONLY : getin_p
15	USE print_control_mod, ONLY: lunout
16
17
18	IMPLICIT NONE
19	! ======================================================================
20	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
21	! O. Boucher (LMD/CNRS) mise a jour en 201212
22	! I. Musat (LMD/CNRS) : prise en compte de la meme hypothese de recouvrement
23	! pour les nuages que pour le rayonnement rrtm via
24	! le parametre novlp de radopt.h : 20160721
25	! Objet: Calculer epaisseur optique et emmissivite des nuages
26	! ======================================================================
27	! Arguments:
28	! ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols
29
30	! t-------input-R-temperature
31	! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie
32	! nuageuse (kg/kg)
33	! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
34	! mass_solu_aero-----input-R-total mass concentration for all soluble
35	! aerosols[ug/m^3]
36	! mass_solu_aero_pi--input-R-ditto, pre-industrial value
37
38	! bl95_b0-input-R-a PARAMETER, may be varied for tests (s-sea, l-land)
39	! bl95_b1-input-R-a PARAMETER, may be varied for tests ( -"- )
40
41	! re------output-R-Cloud droplet effective radius multiplied by fl [um]
42	! fl------output-R-Denominator to re, introduced to avoid problems in
43	! the averaging of the output. fl is the fraction of liquid
44	! water clouds within a grid cell
45
46	! pcltau--output-R-epaisseur optique des nuages
47	! pclemi--output-R-emissivite des nuages (0 a 1)
48	! pcldtaupi-output-R-pre-industrial value of cloud optical thickness,
49
50	! pcl-output-R-2D low-level cloud cover
51	! pcm-output-R-2D mid-level cloud cover
52	! pch-output-R-2D high-level cloud cover
53	! pct-output-R-2D total cloud cover
54	! ======================================================================
55
56	include "YOMCST.h"
57	include "nuage.h"
58	include "radepsi.h"
59	include "radopt.h"
60
61	! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
62	! !novlp=1: max-random
63	! !novlp=2: maximum
64	! !novlp=3: random
65	! LOGICAL random, maximum_random, maximum
66	! PARAMETER (random=.FALSE., maximum_random=.TRUE., maximum=.FALSE.)
67
68	LOGICAL, SAVE :: first = .TRUE.
69	!$OMP THREADPRIVATE(FIRST)
70	INTEGER flag_max
71
72	! threshold PARAMETERs
73	REAL thres_tau, thres_neb
74	PARAMETER (thres_tau=0.3, thres_neb=0.001)
75
76	REAL phase3d(klon, klev)
77	REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
78
79	REAL paprs(klon, klev+1)
80	REAL pplay(klon, klev)
81	REAL t(klon, klev)
82	REAL pclc(klon, klev)
83	REAL pqlwp(klon, klev)
84	REAL pcltau(klon, klev)
85	REAL pclemi(klon, klev)
86	REAL pcldtaupi(klon, klev)
87
88	REAL pct(klon)
89	REAL pcl(klon)
90	REAL pcm(klon)
91	REAL pch(klon)
92	REAL pctlwp(klon)
93
94	LOGICAL lo
95
96	! !Abderr modif JL mail du 19.01.2011 18:31
97	! REAL cetahb, cetamb
98	! PARAMETER (cetahb = 0.45, cetamb = 0.80)
99	! Remplacer
100	! cetahb*paprs(i,1) par prmhc
101	! cetamb*paprs(i,1) par prlmc
102	REAL prmhc ! Pressure between medium and high level cloud in Pa
103	REAL prlmc ! Pressure between low and medium level cloud in Pa
104	PARAMETER (prmhc=440.100., prlmc=680.100.)
105
106	INTEGER i, k
107	REAL xflwp(klon), xfiwp(klon)
108	REAL xflwc(klon, klev), xfiwc(klon, klev)
109
110	REAL radius
111
112	REAL coef_froi, coef_chau
113	PARAMETER (coef_chau=0.13, coef_froi=0.09)
114
115	REAL seuil_neb
116	PARAMETER (seuil_neb=0.001)
117
118	! JBM (3/14) nexpo is replaced by exposant_glace
119	! INTEGER nexpo ! exponentiel pour glace/eau
120	! PARAMETER (nexpo=6)
121	! PARAMETER (nexpo=1)
122	! if iflag_t_glace=0, the old values are used:
123	REAL, PARAMETER :: t_glace_min_old = 258.
124	REAL, PARAMETER :: t_glace_max_old = 273.13
125
126	REAL rel, tc, rei
127	REAL k_ice0, k_ice, df
128	PARAMETER (k_ice0=0.005) ! units=m2/g
129	PARAMETER (df=1.66) ! diffusivity factor
130
131	! jq for the aerosol indirect effect
132	! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
133	! jq
134	REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3]
135	REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
136	REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
137	REAL re(klon, klev) ! cloud droplet effective radius [um]
138	REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
139	REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
140
141	REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
142	! within the grid cell)
143
144	INTEGER flag_aerosol
145	LOGICAL ok_cdnc
146	REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
147
148	! jq-end
149	! IM cf. CR:parametres supplementaires
150	REAL zclear(klon)
151	REAL zcloud(klon)
152	REAL zcloudh(klon)
153	REAL zcloudm(klon)
154	REAL zcloudl(klon)
155	REAL rhodz(klon, klev) !--rho*dz pour la couche
156	REAL zrho(klon, klev) !--rho pour la couche
157	REAL dh(klon, klev) !--dz pour la couche
158	REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
159	REAL rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
160	REAL zflwp_var, zfiwp_var
161	REAL d_rei_dt
162
163	! Abderrahmane oct 2009
164	REAL reliq(klon, klev), reice(klon, klev)
165	REAL reliq_pi(klon, klev), reice_pi(klon, klev)
166
167	REAL,SAVE :: cdnc_min=-1.
168	REAL,SAVE :: cdnc_min_m3
169	!$OMP THREADPRIVATE(cdnc_min,cdnc_min_m3)
170	REAL,SAVE :: cdnc_max=-1.
171	REAL,SAVE :: cdnc_max_m3
172	!$OMP THREADPRIVATE(cdnc_max,cdnc_max_m3)
173
174	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
175	! FH : 2011/05/24
176
177	! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
178	! to be used for a temperature in celcius T(°C) < 0
179	! rei=rei_min for T(°C) < -81.4
180
181	! Calcul de la pente de la relation entre rayon effective des cristaux
182	! et la température.
183	! Pour retrouver les résultats numériques de la version d'origine,
184	! on impose 0.71 quand on est proche de 0.71
185
186	if (first) THEN
187	call getin_p('cdnc_min',cdnc_min)
188	cdnc_min_m3=cdnc_min*1.E6
189	IF (cdnc_min_m3<0.) cdnc_min_m3=20.E6 ! astuce pour retrocompatibilite
190	write(lunout,*)'cdnc_min=', cdnc_min_m3/1.E6
191	call getin_p('cdnc_max',cdnc_max)
192	cdnc_max_m3=cdnc_max*1.E6
193	IF (cdnc_max_m3<0.) cdnc_max_m3=1000.E6 ! astuce pour retrocompatibilite
194	write(lunout,*)'cdnc_max=', cdnc_max_m3/1.E6
195	ENDIF
196
197	d_rei_dt = (rei_max-rei_min)/81.4
198	IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71
199	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200
201	! Calculer l'epaisseur optique et l'emmissivite des nuages
202	! IM inversion des DO
203
204	xflwp = 0.D0
205	xfiwp = 0.D0
206	xflwc = 0.D0
207	xfiwc = 0.D0
208
209	reliq = 0.
210	reice = 0.
211	reliq_pi = 0.
212	reice_pi = 0.
213
214	IF (iflag_t_glace.EQ.0) THEN
215	DO k = 1, klev
216	DO i = 1, klon
217	! -layer calculation
218	rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
219	zrho(i, k) = pplay(i, k)/t(i, k)/rd ! kg/m3
220	dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
221	! -Fraction of ice in cloud using a linear transition
222	zfice(i, k) = 1.0 - (t(i,k)-t_glace_min_old)/(t_glace_max_old-t_glace_min_old)
223	zfice(i, k) = min(max(zfice(i,k),0.0), 1.0)
224	! -IM Total Liquid/Ice water content
225	xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
226	xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
227	ENDDO
228	ENDDO
229	ELSE ! of IF (iflag_t_glace.EQ.0)
230	DO k = 1, klev
231	CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
232
233
234	! JBM: icefrac_lsc is now contained icefrac_lsc_mod
235	! zfice(i, k) = icefrac_lsc(t(i,k), t_glace_min, &
236	! t_glace_max, exposant_glace)
237	DO i = 1, klon
238	! -layer calculation
239	rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
240	zrho(i, k) = pplay(i, k)/t(i, k)/rd ! kg/m3
241	dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
242	! -IM Total Liquid/Ice water content
243	xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
244	xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
245	ENDDO
246	ENDDO
247	ENDIF
248
249	IF (ok_cdnc) THEN
250
251	! --we compute cloud properties as a function of the aerosol load
252
253	DO k = 1, klev
254	DO i = 1, klon
255	! Formula "D" of Boucher and Lohmann, Tellus, 1995
256	! Cloud droplet number concentration (CDNC) is restricted
257	! to be within [20, 1000 cm^3]
258
259	! --pre-industrial case
260	cdnc_pi(i, k) = 10.*(bl95_b0+bl95_b1log(max(mass_solu_aero_pi(i,k), &
261	1.E-4))/log(10.))*1.E6 !-m-3
262	cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
263
264	ENDDO
265	ENDDO
266
267	!--flag_aerosol=7 => MACv2SP climatology
268	!--in this case there is an enhancement factor
269	IF (flag_aerosol .EQ. 7) THEN
270
271	!--present-day
272	DO k = 1, klev
273	DO i = 1, klon
274	cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
275	ENDDO
276	ENDDO
277
278	!--standard case
279	ELSE
280
281	DO k = 1, klev
282	DO i = 1, klon
283
284	! Formula "D" of Boucher and Lohmann, Tellus, 1995
285	! Cloud droplet number concentration (CDNC) is restricted
286	! to be within [20, 1000 cm^3]
287
288	! --present-day case
289	cdnc(i, k) = 10.*(bl95_b0+bl95_b1log(max(mass_solu_aero(i,k), &
290	1.E-4))/log(10.))*1.E6 !-m-3
291	cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
292
293	ENDDO
294	ENDDO
295
296	ENDIF !--flag_aerosol
297
298	!--computing cloud droplet size
299	DO k = 1, klev
300	DO i = 1, klon
301
302	! --present-day case
303	rad_chaud(i, k) = 1.1((pqlwp(i,k)pplay(i, &
304	k)/(rdt(i,k)))/(4./3rpi1000.cdnc(i,k)))**(1./3.)
305	rad_chaud(i, k) = max(rad_chaud(i,k)*1.E6, 5.)
306
307	! --pre-industrial case
308	rad_chaud_pi(i, k) = 1.1((pqlwp(i,k)pplay(i, &
309	k)/(rdt(i,k)))/(4./3.rpi1000.cdnc_pi(i,k)))**(1./3.)
310	rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
311
312	! --pre-industrial case
313	! --liquid/ice cloud water paths:
314	IF (pclc(i,k)<=seuil_neb) THEN
315
316	pcldtaupi(i, k) = 0.0
317
318	ELSE
319
320	zflwp_var = 1000.(1.-zfice(i,k))pqlwp(i, k)/pclc(i, k)* &
321	rhodz(i, k)
322	zfiwp_var = 1000.zfice(i, k)pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
323	tc = t(i, k) - 273.15
324	rei = d_rei_dt*tc + rei_max
325	IF (tc<=-81.4) rei = rei_min
326
327	! -- cloud optical thickness :
328	! [for liquid clouds, traditional formula,
329	! for ice clouds, Ebert & Curry (1992)]
330
331	IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
332	pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
333	zfiwp_var*(3.448E-03+2.431/rei)
334
335	ENDIF
336
337	ENDDO
338	ENDDO
339
340	ELSE !--not ok_cdnc
341
342	! -prescribed cloud droplet radius
343
344	DO k = 1, min(3, klev)
345	DO i = 1, klon
346	rad_chaud(i, k) = rad_chau2
347	rad_chaud_pi(i, k) = rad_chau2
348	ENDDO
349	ENDDO
350	DO k = min(3, klev) + 1, klev
351	DO i = 1, klon
352	rad_chaud(i, k) = rad_chau1
353	rad_chaud_pi(i, k) = rad_chau1
354	ENDDO
355	ENDDO
356
357	ENDIF !--ok_cdnc
358
359	! --computation of cloud optical depth and emissivity
360	! --in the general case
361
362	DO k = 1, klev
363	DO i = 1, klon
364
365	IF (pclc(i,k)<=seuil_neb) THEN
366
367	! effective cloud droplet radius (microns) for liquid water clouds:
368	! For output diagnostics cloud droplet effective radius [um]
369	! we multiply here with f * xl (fraction of liquid water
370	! clouds in the grid cell) to avoid problems in the averaging of the
371	! output.
372	! In the output of IOIPSL, derive the REAL cloud droplet
373	! effective radius as re/fl
374
375	fl(i, k) = seuil_neb*(1.-zfice(i,k))
376	re(i, k) = rad_chaud(i, k)*fl(i, k)
377	rel = 0.
378	rei = 0.
379	pclc(i, k) = 0.0
380	pcltau(i, k) = 0.0
381	pclemi(i, k) = 0.0
382
383	ELSE
384
385	! -- liquid/ice cloud water paths:
386
387	zflwp_var = 1000.(1.-zfice(i,k))pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
388	zfiwp_var = 1000.zfice(i, k)pqlwp(i, k)/pclc(i, k)*rhodz(i, k)
389
390	! effective cloud droplet radius (microns) for liquid water clouds:
391	! For output diagnostics cloud droplet effective radius [um]
392	! we multiply here with f * xl (fraction of liquid water
393	! clouds in the grid cell) to avoid problems in the averaging of the
394	! output.
395	! In the output of IOIPSL, derive the REAL cloud droplet
396	! effective radius as re/fl
397
398	fl(i, k) = pclc(i, k)*(1.-zfice(i,k))
399	re(i, k) = rad_chaud(i, k)*fl(i, k)
400
401	rel = rad_chaud(i, k)
402
403	! for ice clouds: as a function of the ambiant temperature
404	! [formula used by Iacobellis and Somerville (2000), with an
405	! asymptotical value of 3.5 microns at T<-81.4 C added to be
406	! consistent with observations of Heymsfield et al. 1986]:
407	! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
408	! rei_max=61.29
409
410	tc = t(i, k) - 273.15
411	rei = d_rei_dt*tc + rei_max
412	IF (tc<=-81.4) rei = rei_min
413
414	! -- cloud optical thickness :
415	! [for liquid clouds, traditional formula,
416	! for ice clouds, Ebert & Curry (1992)]
417
418	IF (zflwp_var==0.) rel = 1.
419	IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
420	pcltau(i, k) = 3.0/2.0(zflwp_var/rel) + zfiwp_var(3.448E-03+2.431/ &
421	rei)
422
423	! -- cloud infrared emissivity:
424	! [the broadband infrared absorption coefficient is PARAMETERized
425	! as a function of the effective cld droplet radius]
426	! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
427
428	k_ice = k_ice0 + 1.0/rei
429
430	pclemi(i, k) = 1.0 - exp(-coef_chauzflwp_var-dfk_ice*zfiwp_var)
431
432	ENDIF
433
434	reice(i, k) = rei
435
436	xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
437	xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
438
439	ENDDO
440	ENDDO
441
442	! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
443
444	IF (.NOT. ok_cdnc) THEN
445	DO k = 1, klev
446	DO i = 1, klon
447	pcldtaupi(i, k) = pcltau(i, k)
448	reice_pi(i, k) = reice(i, k)
449	ENDDO
450	ENDDO
451	ENDIF
452
453	DO k = 1, klev
454	DO i = 1, klon
455	reliq(i, k) = rad_chaud(i, k)
456	reliq_pi(i, k) = rad_chaud_pi(i, k)
457	reice_pi(i, k) = reice(i, k)
458	ENDDO
459	ENDDO
460
461	! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
462	! IM cf. CR:test: calcul prenant ou non en compte le recouvrement
463	! initialisations
464
465	DO i = 1, klon
466	zclear(i) = 1.
467	zcloud(i) = 0.
468	zcloudh(i) = 0.
469	zcloudm(i) = 0.
470	zcloudl(i) = 0.
471	pch(i) = 1.0
472	pcm(i) = 1.0
473	pcl(i) = 1.0
474	pctlwp(i) = 0.0
475	ENDDO
476
477	! --calculation of liquid water path
478
479	DO k = klev, 1, -1
480	DO i = 1, klon
481	pctlwp(i) = pctlwp(i) + pqlwp(i, k)*rhodz(i, k)
482	ENDDO
483	ENDDO
484
485	! --calculation of cloud properties with cloud overlap
486
487	IF (novlp==1) THEN
488	DO k = klev, 1, -1
489	DO i = 1, klon
490	zclear(i) = zclear(i)*(1.-max(pclc(i,k),zcloud(i)))/(1.-min(real( &
491	zcloud(i),kind=8),1.-zepsec))
492	pct(i) = 1. - zclear(i)
493	IF (paprs(i,k)<prmhc) THEN
494	pch(i) = pch(i)*(1.-max(pclc(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
495	(i),kind=8),1.-zepsec))
496	zcloudh(i) = pclc(i, k)
497	ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
498	pcm(i) = pcm(i)*(1.-max(pclc(i,k),zcloudm(i)))/(1.-min(real(zcloudm &
499	(i),kind=8),1.-zepsec))
500	zcloudm(i) = pclc(i, k)
501	ELSE IF (paprs(i,k)>=prlmc) THEN
502	pcl(i) = pcl(i)*(1.-max(pclc(i,k),zcloudl(i)))/(1.-min(real(zcloudl &
503	(i),kind=8),1.-zepsec))
504	zcloudl(i) = pclc(i, k)
505	ENDIF
506	zcloud(i) = pclc(i, k)
507	ENDDO
508	ENDDO
509	ELSE IF (novlp==2) THEN
510	DO k = klev, 1, -1
511	DO i = 1, klon
512	zcloud(i) = max(pclc(i,k), zcloud(i))
513	pct(i) = zcloud(i)
514	IF (paprs(i,k)<prmhc) THEN
515	pch(i) = min(pclc(i,k), pch(i))
516	ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
517	pcm(i) = min(pclc(i,k), pcm(i))
518	ELSE IF (paprs(i,k)>=prlmc) THEN
519	pcl(i) = min(pclc(i,k), pcl(i))
520	ENDIF
521	ENDDO
522	ENDDO
523	ELSE IF (novlp==3) THEN
524	DO k = klev, 1, -1
525	DO i = 1, klon
526	zclear(i) = zclear(i)*(1.-pclc(i,k))
527	pct(i) = 1 - zclear(i)
528	IF (paprs(i,k)<prmhc) THEN
529	pch(i) = pch(i)*(1.0-pclc(i,k))
530	ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
531	pcm(i) = pcm(i)*(1.0-pclc(i,k))
532	ELSE IF (paprs(i,k)>=prlmc) THEN
533	pcl(i) = pcl(i)*(1.0-pclc(i,k))
534	ENDIF
535	ENDDO
536	ENDDO
537	ENDIF
538
539	DO i = 1, klon
540	pch(i) = 1. - pch(i)
541	pcm(i) = 1. - pcm(i)
542	pcl(i) = 1. - pcl(i)
543	ENDDO
544
545	! ========================================================
546	! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
547	! ========================================================
548	! change by Nicolas Yan (LSCE)
549	! Cloud Droplet Number Concentration (CDNC) : 3D variable
550	! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
551	! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
552	! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
553	! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
554
555	IF (ok_cdnc) THEN
556
557	DO k = 1, klev
558	DO i = 1, klon
559	phase3d(i, k) = 1 - zfice(i, k)
560	IF (pclc(i,k)<=seuil_neb) THEN
561	lcc3d(i, k) = seuil_neb*phase3d(i, k)
562	ELSE
563	lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
564	ENDIF
565	scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
566	ENDDO
567	ENDDO
568
569	DO i = 1, klon
570	lcc(i) = 0.
571	reffclwtop(i) = 0.
572	cldncl(i) = 0.
573	IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
574	IF (novlp.EQ.2) tcc(i) = 0.
575	ENDDO
576
577	DO i = 1, klon
578	DO k = klev - 1, 1, -1 !From TOA down
579
580	! Test, if the cloud optical depth exceeds the necessary
581	! threshold:
582
583	IF (pcltau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
584
585	IF (novlp.EQ.2) THEN
586	IF (first) THEN
587	WRITE (, ) 'Hypothese de recouvrement: MAXIMUM'
588	first = .FALSE.
589	ENDIF
590	flag_max = -1.
591	ftmp(i) = max(tcc(i), pclc(i,k))
592	ENDIF
593
594	IF (novlp.EQ.3) THEN
595	IF (first) THEN
596	WRITE (, ) 'Hypothese de recouvrement: RANDOM'
597	first = .FALSE.
598	ENDIF
599	flag_max = 1.
600	ftmp(i) = tcc(i)*(1-pclc(i,k))
601	ENDIF
602
603	IF (novlp.EQ.1) THEN
604	IF (first) THEN
605	WRITE (, ) 'Hypothese de recouvrement: MAXIMUM_ &
606	& &
607	& RANDOM'
608	first = .FALSE.
609	ENDIF
610	flag_max = 1.
611	ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
612	k+1),1.-thres_neb))
613	ENDIF
614	! Effective radius of cloud droplet at top of cloud (m)
615	reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)1.0E-06phase3d(i, &
616	k)(tcc(i)-ftmp(i))flag_max
617	! CDNC at top of cloud (m-3)
618	cldncl(i) = cldncl(i) + cdnc(i, k)phase3d(i, k)(tcc(i)-ftmp(i))* &
619	flag_max
620	! Liquid Cloud Content at top of cloud
621	lcc(i) = lcc(i) + phase3d(i, k)(tcc(i)-ftmp(i))flag_max
622	! Total Cloud Content at top of cloud
623	tcc(i) = ftmp(i)
624
625	ENDIF ! is there a visible, not-too-small cloud?
626	ENDDO ! loop over k
627
628	IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
629
630	ENDDO ! loop over i
631
632	! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
633	! REFFCLWS)
634	DO i = 1, klon
635	DO k = 1, klev
636	! Weight to be used for outputs: eau_liquide*couverture nuageuse
637	lcc3dcon(i, k) = rnebcon(i, k)phase3d(i, k)clwcon(i, k) ! eau liquide convective
638	lcc3dstra(i, k) = pclc(i, k)pqlwp(i, k)phase3d(i, k)
639	lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
640	lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
641	!FC pour la glace (CAUSES)
642	icc3dcon(i, k) = rnebcon(i, k)(1-phase3d(i, k))clwcon(i, k) ! glace convective
643	icc3dstra(i, k)= pclc(i, k)pqlwp(i, k)(1-phase3d(i, k))
644	icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
645	icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
646	!FC (CAUSES)
647
648	! Compute cloud droplet radius as above in meter
649	radius = 1.1((pqlwp(i,k)pplay(i,k)/(rdt(i,k)))/(4./3rpi1000. &
650	cdnc(i,k)))**(1./3.)
651	radius = max(radius, 5.E-6)
652	! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
653	reffclwc(i, k) = radius
654	reffclwc(i, k) = reffclwc(i, k)*lcc3dcon(i, k)
655	! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
656	reffclws(i, k) = radius
657	reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
658	ENDDO !klev
659	ENDDO !klon
660
661	! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
662
663	DO i = 1, klon
664	cldnvi(i) = 0.
665	lcc_integrat(i) = 0.
666	height(i) = 0.
667	DO k = 1, klev
668	cldnvi(i) = cldnvi(i) + cdnc(i, k)lcc3d(i, k)dh(i, k)
669	lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
670	height(i) = height(i) + dh(i, k)
671	ENDDO ! klev
672	lcc_integrat(i) = lcc_integrat(i)/height(i)
673	IF (lcc_integrat(i)<=1.0E-03) THEN
674	cldnvi(i) = cldnvi(i)*lcc(i)/seuil_neb
675	ELSE
676	cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
677	ENDIF
678	ENDDO ! klon
679
680	DO i = 1, klon
681	DO k = 1, klev
682	IF (scdnc(i,k)<=0.0) scdnc(i, k) = 0.0
683	IF (reffclws(i,k)<=0.0) reffclws(i, k) = 0.0
684	IF (reffclwc(i,k)<=0.0) reffclwc(i, k) = 0.0
685	IF (lcc3d(i,k)<=0.0) lcc3d(i, k) = 0.0
686	IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
687	IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
688	!FC (CAUSES)
689	IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
690	IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
691	!FC (CAUSES)
692	ENDDO
693	IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
694	IF (cldncl(i)<=0.0) cldncl(i) = 0.0
695	IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
696	IF (lcc(i)<=0.0) lcc(i) = 0.0
697	ENDDO
698
699	ENDIF !ok_cdnc
700
701	first=.false. !to be sure
702
703	RETURN
704
705	END SUBROUTINE newmicro

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format