Context Navigation

lmdz_simu_airs.f90 @ 5151

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago
Put YOMCST.h into modules
File size: 35.8 KB

Line
1	module lmdz_simu_airs
2
3	USE lmdz_print_control, ONLY: prt_level, lunout
4	USE lmdz_abort_physic, ONLY: abort_physic
5
6	IMPLICIT NONE
7
8	REAL, PARAMETER :: tau_thresh = 0.05 ! seuil nuages detectables
9	REAL, PARAMETER :: p_thresh = 445. ! seuil nuages hauts
10	REAL, PARAMETER :: em_min = 0.2 ! seuils nuages semi-transp
11	REAL, PARAMETER :: em_max = 0.85
12	REAL, PARAMETER :: undef = 999.
13
14	CONTAINS
15
16	REAL function search_tropopause(P, T, alt, N) result(P_tropo)
17	! this function searches for the tropopause pressure in [hPa].
18	! The search is based on ideology described in
19	! Reichler et al., Determining the tropopause height from gridded data,
20	! GRL, 30(20) 2042, doi:10.1029/2003GL018240, 2003
21
22	INTEGER N, i, i_lev, first_point, exit_flag, i_dir
23	REAL P(N), T(N), alt(N), slope(N)
24	REAL P_min, P_max, slope_limit, slope_2km, &
25	delta_alt_limit, tmp, delta_alt
26	PARAMETER(P_min = 75.0, P_max = 470.0) ! hPa
27	PARAMETER(slope_limit = 0.002) ! 2 K/km converted to K/m
28	PARAMETER(delta_alt_limit = 2000.0) ! 2000 meters
29
30	do i = 1, N - 1
31	slope(i) = -(T(i + 1) - T(i)) / (alt(i + 1) - alt(i))
32	END DO
33	slope(N) = slope(N - 1)
34
35	IF (P(1)>P(N)) THEN
36	i_dir = 1
37	i = 1
38	else
39	i_dir = -1
40	i = N
41	end if
42
43	first_point = 0
44	exit_flag = 0
45	do while(exit_flag==0)
46	IF (P(i)>=P_min.AND.P(i)<=P_max) THEN
47	IF (first_point>0) THEN
48	delta_alt = alt(i) - alt(first_point)
49	IF (delta_alt>=delta_alt_limit) THEN
50	slope_2km = (T(first_point) - T(i)) / delta_alt
51	IF (slope_2km<slope_limit) THEN
52	exit_flag = 1
53	else
54	first_point = 0
55	end if
56	end if
57	end if
58	IF (first_point==0.AND.slope(i)<slope_limit) first_point = i
59	end if
60	i = i + i_dir
61	IF (i<=1.OR.i>=N) exit_flag = 1
62	END DO
63
64	IF (first_point<=0) P_tropo = 65.4321
65	IF (first_point==1) P_tropo = 64.5432
66	IF (first_point>1) THEN
67	tmp = (slope_limit - slope(first_point)) / (slope(first_point + 1) - &
68	slope(first_point)) * (P(first_point + 1) - P(first_point))
69	P_tropo = P(first_point) + tmp
70	! PRINT*, 'P_tropo= ', tmp, P(first_point), P_tropo
71	end if
72
73	! Ajout Marine
74	IF (P_tropo < 60. .OR. P_tropo > 470.) THEN
75	P_tropo = 999.
76	endif
77
78	END FUNCTION search_tropopause
79
80
81	SUBROUTINE cloud_structure(len_cs, rneb_cs, temp_cs, &
82	emis_cs, iwco_cs, &
83	pres_cs, dz_cs, rhodz_cs, rad_cs, &
84	cc_tot_cs, cc_hc_cs, cc_hist_cs, &
85	cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
86	pcld_hc_cs, tcld_hc_cs, &
87	em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
88	pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
89	pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
90	pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
91	em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
92
93	INTEGER :: i, n, nss
94
95	INTEGER, INTENT(IN) :: len_cs
96	REAL, DIMENSION(:), INTENT(IN) :: rneb_cs, temp_cs
97	REAL, DIMENSION(:), INTENT(IN) :: emis_cs, iwco_cs, rad_cs
98	REAL, DIMENSION(:), INTENT(IN) :: pres_cs, dz_cs, rhodz_cs
99
100	REAL, INTENT(OUT) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, &
101	cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
102	pcld_hc_cs, tcld_hc_cs, em_hc_cs, iwp_hc_cs, &
103	pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
104	pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
105	pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
106	em_hist_cs, iwp_hist_cs, &
107	deltaz_hc_cs, deltaz_hist_cs, rad_hist_cs
108
109	REAL, DIMENSION(len_cs) :: rneb_ord
110	REAL :: rneb_min
111
112	REAL, DIMENSION(:), ALLOCATABLE :: s, s_hc, s_hist, rneb_max
113	REAL, DIMENSION(:), ALLOCATABLE :: sCb, sThCi, sAnv
114	REAL, DIMENSION(:), ALLOCATABLE :: iwp_ss, pcld_ss, tcld_ss, &
115	emis_ss
116	REAL, DIMENSION(:), ALLOCATABLE :: deltaz_ss, rad_ss
117
118	CHARACTER (len = 50) :: modname = 'simu_airs.cloud_structure'
119	CHARACTER (len = 160) :: abort_message
120
121	WRITE(lunout, *) 'dans cloud_structure'
122
123	CALL ordonne(len_cs, rneb_cs, rneb_ord)
124
125
126	! Definition des sous_sections
127
128	rneb_min = rneb_ord(1)
129	nss = 1
130
131	do i = 2, size(rneb_cs)
132	IF (rneb_ord(i) > rneb_min) THEN
133	nss = nss + 1
134	rneb_min = rneb_ord(i)
135	endif
136	enddo
137
138	allocate (s(nss))
139	allocate (s_hc(nss))
140	allocate (s_hist(nss))
141	allocate (rneb_max(nss))
142	allocate (emis_ss(nss))
143	allocate (pcld_ss(nss))
144	allocate (tcld_ss(nss))
145	allocate (iwp_ss(nss))
146	allocate (deltaz_ss(nss))
147	allocate (rad_ss(nss))
148	allocate (sCb(nss))
149	allocate (sThCi(nss))
150	allocate (sAnv(nss))
151
152	rneb_min = rneb_ord(1)
153	n = 1
154	s(1) = rneb_ord(1)
155	s_hc(1) = rneb_ord(1)
156	s_hist(1) = rneb_ord(1)
157	sCb(1) = rneb_ord(1)
158	sThCi(1) = rneb_ord(1)
159	sAnv(1) = rneb_ord(1)
160
161	rneb_max(1) = rneb_ord(1)
162
163	do i = 2, size(rneb_cs)
164	IF (rneb_ord(i) > rneb_min) THEN
165	n = n + 1
166	s(n) = rneb_ord(i) - rneb_min
167	s_hc(n) = rneb_ord(i) - rneb_min
168	s_hist(n) = rneb_ord(i) - rneb_min
169	sCb(n) = rneb_ord(i) - rneb_min
170	sThCi(n) = rneb_ord(i) - rneb_min
171	sAnv(n) = rneb_ord(i) - rneb_min
172
173	rneb_max(n) = rneb_ord(i)
174	rneb_min = rneb_ord(i)
175	endif
176	enddo
177
178	! Appel de sous_section
179
180	do i = 1, nss
181	CALL sous_section(len_cs, rneb_cs, temp_cs, &
182	emis_cs, iwco_cs, &
183	pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
184	rneb_max(i), s(i), s_hc(i), s_hist(i), &
185	sCb(i), sThCi(i), sAnv(i), &
186	emis_ss(i), &
187	pcld_ss(i), tcld_ss(i), iwp_ss(i), deltaz_ss(i), rad_ss(i))
188	enddo
189
190	! Caracteristiques de la structure nuageuse
191
192	cc_tot_cs = 0.
193	cc_hc_cs = 0.
194	cc_hist_cs = 0.
195
196	cc_Cb_cs = 0.
197	cc_ThCi_cs = 0.
198	cc_Anv_cs = 0.
199
200	em_hc_cs = 0.
201	iwp_hc_cs = 0.
202	deltaz_hc_cs = 0.
203
204	em_hist_cs = 0.
205	iwp_hist_cs = 0.
206	deltaz_hist_cs = 0.
207	rad_hist_cs = 0.
208
209	pcld_hc_cs = 0.
210	tcld_hc_cs = 0.
211
212	pcld_Cb_cs = 0.
213	tcld_Cb_cs = 0.
214	em_Cb_cs = 0.
215
216	pcld_ThCi_cs = 0.
217	tcld_ThCi_cs = 0.
218	em_ThCi_cs = 0.
219
220	pcld_Anv_cs = 0.
221	tcld_Anv_cs = 0.
222	em_Anv_cs = 0.
223
224	do i = 1, nss
225
226	cc_tot_cs = cc_tot_cs + s(i)
227	cc_hc_cs = cc_hc_cs + s_hc(i)
228	cc_hist_cs = cc_hist_cs + s_hist(i)
229
230	cc_Cb_cs = cc_Cb_cs + sCb(i)
231	cc_ThCi_cs = cc_ThCi_cs + sThCi(i)
232	cc_Anv_cs = cc_Anv_cs + sAnv(i)
233
234	iwp_hc_cs = iwp_hc_cs + s_hc(i) * iwp_ss(i)
235	em_hc_cs = em_hc_cs + s_hc(i) * emis_ss(i)
236	deltaz_hc_cs = deltaz_hc_cs + s_hc(i) * deltaz_ss(i)
237
238	iwp_hist_cs = iwp_hist_cs + s_hist(i) * iwp_ss(i)
239	em_hist_cs = em_hist_cs + s_hist(i) * emis_ss(i)
240	deltaz_hist_cs = deltaz_hist_cs + s_hist(i) * deltaz_ss(i)
241	rad_hist_cs = rad_hist_cs + s_hist(i) * rad_ss(i)
242
243	pcld_hc_cs = pcld_hc_cs + s_hc(i) * pcld_ss(i)
244	tcld_hc_cs = tcld_hc_cs + s_hc(i) * tcld_ss(i)
245
246	pcld_Cb_cs = pcld_Cb_cs + sCb(i) * pcld_ss(i)
247	tcld_Cb_cs = tcld_Cb_cs + sCb(i) * tcld_ss(i)
248	em_Cb_cs = em_Cb_cs + sCb(i) * emis_ss(i)
249
250	pcld_ThCi_cs = pcld_ThCi_cs + sThCi(i) * pcld_ss(i)
251	tcld_ThCi_cs = tcld_ThCi_cs + sThCi(i) * tcld_ss(i)
252	em_ThCi_cs = em_ThCi_cs + sThCi(i) * emis_ss(i)
253
254	pcld_Anv_cs = pcld_Anv_cs + sAnv(i) * pcld_ss(i)
255	tcld_Anv_cs = tcld_Anv_cs + sAnv(i) * tcld_ss(i)
256	em_Anv_cs = em_Anv_cs + sAnv(i) * emis_ss(i)
257
258	enddo
259
260	deallocate(s)
261	deallocate (s_hc)
262	deallocate (s_hist)
263	deallocate (rneb_max)
264	deallocate (emis_ss)
265	deallocate (pcld_ss)
266	deallocate (tcld_ss)
267	deallocate (iwp_ss)
268	deallocate (deltaz_ss)
269	deallocate (rad_ss)
270	deallocate (sCb)
271	deallocate (sThCi)
272	deallocate (sAnv)
273
274	CALL normal_undef(pcld_hc_cs, cc_hc_cs)
275	CALL normal_undef(tcld_hc_cs, cc_hc_cs)
276	CALL normal_undef(iwp_hc_cs, cc_hc_cs)
277	CALL normal_undef(em_hc_cs, cc_hc_cs)
278	CALL normal_undef(deltaz_hc_cs, cc_hc_cs)
279
280	CALL normal_undef(iwp_hist_cs, cc_hist_cs)
281	CALL normal_undef(em_hist_cs, cc_hist_cs)
282	CALL normal_undef(deltaz_hist_cs, cc_hist_cs)
283	CALL normal_undef(rad_hist_cs, cc_hist_cs)
284
285	CALL normal_undef(pcld_Cb_cs, cc_Cb_cs)
286	CALL normal_undef(tcld_Cb_cs, cc_Cb_cs)
287	CALL normal_undef(em_Cb_cs, cc_Cb_cs)
288
289	CALL normal_undef(pcld_ThCi_cs, cc_ThCi_cs)
290	CALL normal_undef(tcld_ThCi_cs, cc_ThCi_cs)
291	CALL normal_undef(em_ThCi_cs, cc_ThCi_cs)
292
293	CALL normal_undef(pcld_Anv_cs, cc_Anv_cs)
294	CALL normal_undef(tcld_Anv_cs, cc_Anv_cs)
295	CALL normal_undef(em_Anv_cs, cc_Anv_cs)
296
297
298	! Tests
299
300	IF (cc_tot_cs > maxval(rneb_cs) .AND. &
301	abs(cc_tot_cs - maxval(rneb_cs)) > 1.e-4) THEN
302	WRITE(abort_message, *) 'cc_tot_cs > max rneb_cs', cc_tot_cs, maxval(rneb_cs)
303	CALL abort_physic(modname, abort_message, 1)
304	endif
305
306	IF (iwp_hc_cs < 0.) THEN
307	abort_message = 'cloud_structure: iwp_hc_cs < 0'
308	CALL abort_physic(modname, abort_message, 1)
309	endif
310
311	END SUBROUTINE cloud_structure
312
313
314	SUBROUTINE normal_undef(num, den)
315
316	REAL, INTENT(IN) :: den
317	REAL, INTENT(INOUT) :: num
318
319	IF (den /= 0) THEN
320	num = num / den
321	else
322	num = undef
323	endif
324
325	END SUBROUTINE normal_undef
326
327
328	SUBROUTINE normal2_undef(res, num, den)
329
330	REAL, INTENT(IN) :: den
331	REAL, INTENT(IN) :: num
332	REAL, INTENT(OUT) :: res
333
334	IF (den /= 0.) THEN
335	res = num / den
336	else
337	res = undef
338	endif
339
340	END SUBROUTINE normal2_undef
341
342
343	SUBROUTINE sous_section(len_cs, rneb_cs, temp_cs, &
344	emis_cs, iwco_cs, &
345	pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
346	rnebmax, stot, shc, shist, &
347	sCb, sThCi, sAnv, &
348	emis, pcld, tcld, iwp, deltaz, rad)
349
350	INTEGER, INTENT(IN) :: len_cs
351	REAL, DIMENSION(len_cs), INTENT(IN) :: rneb_cs, temp_cs
352	REAL, DIMENSION(len_cs), INTENT(IN) :: emis_cs, iwco_cs, &
353	rneb_ord
354	REAL, DIMENSION(len_cs), INTENT(IN) :: pres_cs, dz_cs, rad_cs
355	REAL, DIMENSION(len_cs), INTENT(IN) :: rhodz_cs
356	REAL, DIMENSION(len_cs) :: tau_cs, w
357	REAL, INTENT(IN) :: rnebmax
358	REAL, INTENT(INOUT) :: stot, shc, shist
359	REAL, INTENT(INOUT) :: sCb, sThCi, sAnv
360	REAL, INTENT(OUT) :: emis, pcld, tcld, iwp, deltaz, rad
361
362	INTEGER :: i, ideb, ibeg, iend, nuage, visible
363	REAL :: som, som_tau, som_iwc, som_dz, som_rad, tau
364
365	CHARACTER (len = 50) :: modname = 'simu_airs.sous_section'
366	CHARACTER (len = 160) :: abort_message
367
368
369	! Ponderation: 1 pour les nuages, 0 pour les trous
370
371	do i = 1, len_cs
372	IF (rneb_cs(i) >= rnebmax) THEN
373	w(i) = 1.
374	else
375	w(i) = 0.
376	endif
377	enddo
378
379	! Calcul des epaisseurs optiques a partir des emissivites
380
381	som = 0.
382	do i = 1, len_cs
383	IF (emis_cs(i) == 1.) THEN
384	tau_cs(i) = 10.
385	else
386	tau_cs(i) = -log(1. - emis_cs(i))
387	endif
388	som = som + tau_cs(i)
389	enddo
390
391	ideb = 1
392	nuage = 0
393	visible = 0
394
395
396	! Boucle sur les nuages
397	do while (ideb /= 0 .AND. ideb <= len_cs)
398
399
400	! Definition d'un nuage de la sous-section
401
402	CALL topbot(ideb, w, ibeg, iend)
403	ideb = iend + 1
404
405	IF (ibeg > 0) THEN
406	nuage = nuage + 1
407
408	! On determine les caracteristiques du nuage
409	! (ep. optique, ice water path, pression, temperature)
410
411	CALL caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
412	pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
413	som_tau, som_iwc, som_dz, som_rad)
414
415	! On masque le nuage s'il n'est pas detectable
416
417	CALL masque(ibeg, iend, som_tau, visible, w)
418
419	endif
420
421	! Fin boucle nuages
422	enddo
423
424
425	! Analyse du nuage detecte
426
427	CALL topbot(1, w, ibeg, iend)
428
429	IF (ibeg > 0) THEN
430	CALL caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
431	pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
432	som_tau, som_iwc, som_dz, som_rad)
433
434	tau = som_tau
435	emis = 1. - exp(-tau)
436	iwp = som_iwc
437	deltaz = som_dz
438	rad = som_rad
439
440	IF (pcld > p_thresh) THEN
441	shc = 0.
442	shist = 0.
443	sCb = 0.
444	sThCi = 0.
445	sAnv = 0.
446
447	else
448
449	IF (emis < em_min .OR. emis > em_max &
450	.OR. tcld > 230.) THEN
451	shist = 0.
452	endif
453
454	IF (emis < 0.98) THEN
455	sCb = 0.
456	endif
457
458	IF (emis > 0.5 .OR. emis < 0.1) THEN
459	sThCi = 0.
460	endif
461
462	IF (emis <= 0.5 .OR. emis >= 0.98) THEN
463	sAnv = 0.
464	endif
465
466	endif
467
468	else
469
470	tau = 0.
471	emis = 0.
472	iwp = 0.
473	deltaz = 0.
474	pcld = 0.
475	tcld = 0.
476	stot = 0.
477	shc = 0.
478	shist = 0.
479	rad = 0.
480	sCb = 0.
481	sThCi = 0.
482	sAnv = 0.
483
484	endif
485
486
487	! Tests
488
489	IF (iwp < 0.) THEN
490	WRITE(abort_message, *) 'ideb iwp =', ideb, iwp
491	CALL abort_physic(modname, abort_message, 1)
492	endif
493
494	IF (deltaz < 0.) THEN
495	WRITE(abort_message, *)'ideb deltaz =', ideb, deltaz
496	CALL abort_physic(modname, abort_message, 1)
497	endif
498
499	IF (emis < 0.048 .AND. emis /= 0.) THEN
500	WRITE(abort_message, *) 'ideb emis =', ideb, emis
501	CALL abort_physic(modname, abort_message, 1)
502	endif
503
504	END SUBROUTINE sous_section
505
506
507	SUBROUTINE masque(ibeg, iend, som_tau, &
508	visible, w)
509
510	INTEGER, INTENT(IN) :: ibeg, iend
511	REAL, INTENT(IN) :: som_tau
512
513	INTEGER, INTENT(INOUT) :: visible
514	REAL, DIMENSION(:), INTENT(INOUT) :: w
515
516	INTEGER :: i
517
518
519
520	! Masque
521
522	! Cas ou il n'y a pas de nuage visible au-dessus
523
524	IF (visible == 0) THEN
525	IF (som_tau < tau_thresh) THEN
526	do i = ibeg, iend
527	w(i) = 0.
528	enddo
529	else
530	visible = 1
531	endif
532
533	! Cas ou il y a un nuage visible au-dessus
534
535	else
536
537	do i = ibeg, iend
538	w(i) = 0.
539	enddo
540
541	endif
542
543	END SUBROUTINE masque
544
545
546	SUBROUTINE caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
547	pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
548	som_tau, som_iwc, som_dz, som_rad)
549
550	INTEGER, INTENT(IN) :: ibeg, iend
551	REAL, DIMENSION(:), INTENT(IN) :: tau_cs, iwco_cs, temp_cs
552	REAL, DIMENSION(:), INTENT(IN) :: pres_cs, dz_cs, rad_cs
553	REAL, DIMENSION(:), INTENT(IN) :: rhodz_cs
554	REAL, INTENT(OUT) :: som_tau, som_iwc, som_dz, som_rad
555	REAL, INTENT(OUT) :: pcld, tcld
556
557	INTEGER :: i, ibase, imid
558
559	CHARACTER (len = 50) :: modname = 'simu_airs.caract'
560	CHARACTER (len = 160) :: abort_message
561
562	! Somme des epaisseurs optiques et des contenus en glace sur le nuage
563
564	som_tau = 0.
565	som_iwc = 0.
566	som_dz = 0.
567	som_rad = 0.
568	ibase = -100
569
570	do i = ibeg, iend
571
572	som_tau = som_tau + tau_cs(i)
573
574	som_dz = som_dz + dz_cs(i)
575	som_iwc = som_iwc + iwco_cs(i) * 1000 * rhodz_cs(i) ! en g/m2
576	som_rad = som_rad + rad_cs(i) * dz_cs(i)
577
578	IF (som_tau > 3. .AND. ibase == -100) THEN
579	ibase = i - 1
580	endif
581
582	enddo
583
584	IF (som_dz /= 0.) THEN
585	som_rad = som_rad / som_dz
586	else
587	WRITE(, ) 'som_dez = 0 stop '
588	WRITE(, ) 'ibeg, iend =', ibeg, iend
589	do i = ibeg, iend
590	WRITE(, ) dz_cs(i), rhodz_cs(i)
591	enddo
592	abort_message = 'see above'
593	CALL abort_physic(modname, abort_message, 1)
594	endif
595
596	! Determination de Pcld et Tcld
597
598	IF (ibase < ibeg) THEN
599	ibase = ibeg
600	endif
601
602	imid = (ibeg + ibase) / 2
603
604	pcld = pres_cs(imid) / 100. ! pcld en hPa
605	tcld = temp_cs(imid)
606
607	END SUBROUTINE caract
608
609	SUBROUTINE topbot(ideb, w, ibeg, iend)
610
611	INTEGER, INTENT(IN) :: ideb
612	REAL, DIMENSION(:), INTENT(IN) :: w
613	INTEGER, INTENT(OUT) :: ibeg, iend
614
615	INTEGER :: i, itest
616
617	itest = 0
618	ibeg = 0
619	iend = 0
620
621	do i = ideb, size(w)
622
623	IF (w(i) == 1. .AND. itest == 0) THEN
624	ibeg = i
625	itest = 1
626	endif
627
628	enddo
629
630	i = ibeg
631	do while (w(i) == 1. .AND. i <= size(w))
632	i = i + 1
633	enddo
634	iend = i - 1
635
636	END SUBROUTINE topbot
637
638	SUBROUTINE ordonne(len_cs, rneb_cs, rneb_ord)
639
640	INTEGER, INTENT(IN) :: len_cs
641	REAL, DIMENSION(:), INTENT(IN) :: rneb_cs
642	REAL, DIMENSION(:), INTENT(OUT) :: rneb_ord
643
644	INTEGER :: i, j, ind_min
645
646	REAL, DIMENSION(len_cs) :: rneb
647	REAL :: rneb_max
648
649	do i = 1, size(rneb_cs)
650	rneb(i) = rneb_cs(i)
651	enddo
652
653	do j = 1, size(rneb_cs)
654
655	rneb_max = 100.
656
657	do i = 1, size(rneb_cs)
658	IF (rneb(i) < rneb_max) THEN
659	rneb_max = rneb(i)
660	ind_min = i
661	endif
662	enddo
663
664	rneb_ord(j) = rneb_max
665	rneb(ind_min) = 100.
666
667	enddo
668
669	END SUBROUTINE ordonne
670
671
672	SUBROUTINE sim_mesh(rneb_1D, temp_1D, emis_1D, &
673	iwcon_1D, rad_1D, &
674	pres, dz, &
675	rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh, &
676	tcld_hc_mesh, &
677	em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, &
678	cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
679	pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
680	pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
681	pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
682	em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
683
684	USE dimphy
685
686	REAL, DIMENSION(klev), INTENT(IN) :: rneb_1D, temp_1D, emis_1D, &
687	iwcon_1D, rad_1D
688	REAL, DIMENSION(klev), INTENT(IN) :: pres, dz, rhodz_1D
689	REAL, INTENT(OUT) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
690	REAL, INTENT(OUT) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
691	REAL, INTENT(OUT) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, &
692	iwp_hc_mesh
693
694	REAL, INTENT(OUT) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
695	REAL, INTENT(OUT) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
696	em_ThCi_mesh
697	REAL, INTENT(OUT) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
698
699	REAL, INTENT(OUT) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
700	REAL, INTENT(OUT) :: deltaz_hc_mesh, deltaz_hist_mesh
701
702	REAL, DIMENSION(:), ALLOCATABLE :: rneb_cs, temp_cs, emis_cs, &
703	iwco_cs
704	REAL, DIMENSION(:), ALLOCATABLE :: pres_cs, dz_cs, rad_cs, &
705	rhodz_cs
706
707	INTEGER :: i, j, l
708	INTEGER :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics
709
710	REAL :: som_emi_hc, som_pcl_hc, som_tcl_hc, som_iwp_hc, som_hc, &
711	som_hist
712	REAL :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
713	som_deltaz_hist
714	REAL :: som_rad_hist
715	REAL :: som_Cb, som_ThCi, som_Anv
716	REAL :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb
717	REAL :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv
718	REAL :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi
719	REAL :: tsom_tot, tsom_hc, tsom_hist
720	REAL :: prod, prod_hh
721
722	REAL :: cc_tot_cs, cc_hc_cs, cc_hist_cs
723	REAL :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs
724	REAL :: pcld_hc_cs, tcld_hc_cs
725	REAL :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs
726	REAL :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs
727	REAL :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs
728	REAL :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs
729	REAL :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs
730
731	REAL, DIMENSION(klev) :: test_tot, test_hc, test_hist
732	REAL, DIMENSION(klev) :: test_pcld, test_tcld, test_em, test_iwp
733
734	CHARACTER (len = 50) :: modname = 'simu_airs.sim_mesh'
735	CHARACTER (len = 160) :: abort_message
736
737	do j = 1, klev
738	WRITE(lunout, *) 'simu_airs, j, rneb_1D =', rneb_1D(j)
739	enddo
740
741	! Definition des structures nuageuses, de la plus haute a la plus basse
742
743	num_cs = 0
744	ltop = klev - 1
745
746	prod = 1.
747
748	som_emi_hc = 0.
749	som_emi_hist = 0.
750	som_pcl_hc = 0.
751	som_tcl_hc = 0.
752	som_iwp_hc = 0.
753	som_iwp_hist = 0.
754	som_deltaz_hc = 0.
755	som_deltaz_hist = 0.
756	som_rad_hist = 0.
757	som_hc = 0.
758	som_hist = 0.
759
760	som_Cb = 0.
761	som_ThCi = 0.
762	som_Anv = 0.
763
764	som_emi_Cb = 0.
765	som_pcld_Cb = 0.
766	som_tcld_Cb = 0.
767
768	som_emi_ThCi = 0.
769	som_pcld_ThCi = 0.
770	som_tcld_ThCi = 0.
771
772	som_emi_Anv = 0.
773	som_pcld_Anv = 0.
774	som_tcld_Anv = 0.
775
776	tsom_tot = 0.
777	tsom_hc = 0.
778	tsom_hist = 0.
779
780	do while (ltop >= 1) ! Boucle sur toute la colonne
781
782	itop = 0
783
784	do l = ltop, 1, -1
785
786	IF (itop == 0 .AND. rneb_1D(l) > 0.001) THEN
787	itop = l
788	endif
789
790	enddo
791
792	ibot = itop
793
794	do while (rneb_1D(ibot) > 0.001 .AND. ibot >= 1)
795	ibot = ibot - 1
796	enddo
797
798	ibot = ibot + 1
799
800	IF (itop > 0) then ! itop > 0
801
802	num_cs = num_cs + 1
803	len_cs = itop - ibot + 1
804
805	! Allocation et definition des variables de la structure nuageuse
806	! le premier indice denote ce qui est le plus haut
807
808	allocate (rneb_cs(len_cs))
809	allocate (temp_cs(len_cs))
810	allocate (emis_cs(len_cs))
811	allocate (iwco_cs(len_cs))
812	allocate (pres_cs(len_cs))
813	allocate (dz_cs(len_cs))
814	allocate (rad_cs(len_cs))
815	allocate (rhodz_cs(len_cs))
816
817	ics = 0
818
819	do i = itop, ibot, -1
820	ics = ics + 1
821	rneb_cs(ics) = rneb_1D(i)
822	temp_cs(ics) = temp_1D(i)
823	emis_cs(ics) = emis_1D(i)
824	iwco_cs(ics) = iwcon_1D(i)
825	rad_cs(ics) = rad_1D(i)
826	pres_cs(ics) = pres(i)
827	dz_cs(ics) = dz(i)
828	rhodz_cs(ics) = rhodz_1D(i)
829	enddo
830
831	! Appel du sous_programme cloud_structure
832
833	CALL cloud_structure(len_cs, rneb_cs, temp_cs, emis_cs, iwco_cs, &
834	pres_cs, dz_cs, rhodz_cs, rad_cs, &
835	cc_tot_cs, cc_hc_cs, cc_hist_cs, &
836	cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
837	pcld_hc_cs, tcld_hc_cs, &
838	em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
839	pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
840	pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
841	pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
842	em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
843
844	deallocate (rneb_cs)
845	deallocate (temp_cs)
846	deallocate (emis_cs)
847	deallocate (iwco_cs)
848	deallocate (pres_cs)
849	deallocate (dz_cs)
850	deallocate (rad_cs)
851	deallocate (rhodz_cs)
852
853
854	! Pour la couverture nuageuse sur la maille
855
856	prod_hh = prod
857
858	prod = prod * (1. - cc_tot_cs)
859
860	! Pour les autres variables definies sur la maille
861
862	som_emi_hc = som_emi_hc + em_hc_cs * cc_hc_cs * prod_hh
863	som_iwp_hc = som_iwp_hc + iwp_hc_cs * cc_hc_cs * prod_hh
864	som_deltaz_hc = som_deltaz_hc + deltaz_hc_cs * cc_hc_cs * prod_hh
865
866	som_emi_Cb = som_emi_Cb + em_Cb_cs * cc_Cb_cs * prod_hh
867	som_tcld_Cb = som_tcld_Cb + tcld_Cb_cs * cc_Cb_cs * prod_hh
868	som_pcld_Cb = som_pcld_Cb + pcld_Cb_cs * cc_Cb_cs * prod_hh
869
870	som_emi_ThCi = som_emi_ThCi + em_ThCi_cs * cc_ThCi_cs * prod_hh
871	som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cs * cc_ThCi_cs * prod_hh
872	som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cs * cc_ThCi_cs * prod_hh
873
874	som_emi_Anv = som_emi_Anv + em_Anv_cs * cc_Anv_cs * prod_hh
875	som_tcld_Anv = som_tcld_Anv + tcld_Anv_cs * cc_Anv_cs * prod_hh
876	som_pcld_Anv = som_pcld_Anv + pcld_Anv_cs * cc_Anv_cs * prod_hh
877
878	som_emi_hist = som_emi_hist + em_hist_cs * cc_hist_cs * prod_hh
879	som_iwp_hist = som_iwp_hist + iwp_hist_cs * cc_hist_cs * prod_hh
880	som_deltaz_hist = som_deltaz_hist + &
881	deltaz_hist_cs * cc_hist_cs * prod_hh
882	som_rad_hist = som_rad_hist + rad_hist_cs * cc_hist_cs * prod_hh
883
884	som_pcl_hc = som_pcl_hc + pcld_hc_cs * cc_hc_cs * prod_hh
885	som_tcl_hc = som_tcl_hc + tcld_hc_cs * cc_hc_cs * prod_hh
886
887	som_hc = som_hc + cc_hc_cs * prod_hh
888	som_hist = som_hist + cc_hist_cs * prod_hh
889
890	som_Cb = som_Cb + cc_Cb_cs * prod_hh
891	som_ThCi = som_ThCi + cc_ThCi_cs * prod_hh
892	som_Anv = som_Anv + cc_Anv_cs * prod_hh
893
894
895	! Pour test
896
897	CALL test_bornes('cc_tot_cs ', cc_tot_cs, 1., 0.)
898	CALL test_bornes('cc_hc_cs ', cc_hc_cs, 1., 0.)
899	CALL test_bornes('cc_hist_cs ', cc_hist_cs, 1., 0.)
900	CALL test_bornes('pcld_hc_cs ', pcld_hc_cs, 1200., 0.)
901	CALL test_bornes('tcld_hc_cs ', tcld_hc_cs, 1000., 100.)
902	CALL test_bornes('em_hc_cs ', em_hc_cs, 1000., 0.048)
903
904	test_tot(num_cs) = cc_tot_cs
905	test_hc(num_cs) = cc_hc_cs
906	test_hist(num_cs) = cc_hist_cs
907	test_pcld(num_cs) = pcld_hc_cs
908	test_tcld(num_cs) = tcld_hc_cs
909	test_em(num_cs) = em_hc_cs
910	test_iwp(num_cs) = iwp_hc_cs
911
912	tsom_tot = tsom_tot + cc_tot_cs
913	tsom_hc = tsom_hc + cc_hc_cs
914	tsom_hist = tsom_hist + cc_hist_cs
915
916	endif ! itop > 0
917
918	ltop = ibot - 1
919
920	enddo ! fin de la boucle sur la colonne
921
922	N_CS = num_cs
923
924
925	! Determination des variables de sortie
926
927	IF (N_CS > 0) then ! if N_CS>0
928
929	cc_tot_mesh = 1. - prod
930
931	cc_hc_mesh = som_hc
932	cc_hist_mesh = som_hist
933
934	cc_Cb_mesh = som_Cb
935	cc_ThCi_mesh = som_ThCi
936	cc_Anv_mesh = som_Anv
937
938	CALL normal2_undef(pcld_hc_mesh, som_pcl_hc, &
939	cc_hc_mesh)
940	CALL normal2_undef(tcld_hc_mesh, som_tcl_hc, &
941	cc_hc_mesh)
942	CALL normal2_undef(em_hc_mesh, som_emi_hc, &
943	cc_hc_mesh)
944	CALL normal2_undef(iwp_hc_mesh, som_iwp_hc, &
945	cc_hc_mesh)
946	CALL normal2_undef(deltaz_hc_mesh, som_deltaz_hc, &
947	cc_hc_mesh)
948
949	CALL normal2_undef(em_Cb_mesh, som_emi_Cb, &
950	cc_Cb_mesh)
951	CALL normal2_undef(tcld_Cb_mesh, som_tcld_Cb, &
952	cc_Cb_mesh)
953	CALL normal2_undef(pcld_Cb_mesh, som_pcld_Cb, &
954	cc_Cb_mesh)
955
956	CALL normal2_undef(em_ThCi_mesh, som_emi_ThCi, &
957	cc_ThCi_mesh)
958	CALL normal2_undef(tcld_ThCi_mesh, som_tcld_ThCi, &
959	cc_ThCi_mesh)
960	CALL normal2_undef(pcld_ThCi_mesh, som_pcld_ThCi, &
961	cc_ThCi_mesh)
962
963	CALL normal2_undef(em_Anv_mesh, som_emi_Anv, &
964	cc_Anv_mesh)
965	CALL normal2_undef(tcld_Anv_mesh, som_tcld_Anv, &
966	cc_Anv_mesh)
967	CALL normal2_undef(pcld_Anv_mesh, som_pcld_Anv, &
968	cc_Anv_mesh)
969
970	CALL normal2_undef(em_hist_mesh, som_emi_hist, &
971	cc_hist_mesh)
972	CALL normal2_undef(iwp_hist_mesh, som_iwp_hist, &
973	cc_hist_mesh)
974	CALL normal2_undef(deltaz_hist_mesh, som_deltaz_hist, &
975	cc_hist_mesh)
976	CALL normal2_undef(rad_hist_mesh, som_rad_hist, &
977	cc_hist_mesh)
978
979
980	! Tests
981
982	! Tests
983
984	IF (cc_tot_mesh > tsom_tot .AND. &
985	abs(cc_tot_mesh - tsom_tot) > 1.e-4) THEN
986	WRITE(abort_message, *)'cc_tot_mesh > tsom_tot', cc_tot_mesh, tsom_tot
987	CALL abort_physic(modname, abort_message, 1)
988	endif
989
990	IF (cc_tot_mesh < maxval(test_tot(1:N_CS)) .AND. &
991	abs(cc_tot_mesh - maxval(test_tot(1:N_CS))) > 1.e-4) THEN
992	WRITE(abort_message, *) 'cc_tot_mesh < max', cc_tot_mesh, maxval(test_tot(1:N_CS))
993	CALL abort_physic(modname, abort_message, 1)
994	endif
995
996	IF (cc_hc_mesh > tsom_hc .AND. &
997	abs(cc_hc_mesh - tsom_hc) > 1.e-4) THEN
998	WRITE(abort_message, *) 'cc_hc_mesh > tsom_hc', cc_hc_mesh, tsom_hc
999	CALL abort_physic(modname, abort_message, 1)
1000	endif
1001
1002	IF (cc_hc_mesh < maxval(test_hc(1:N_CS)) .AND. &
1003	abs(cc_hc_mesh - maxval(test_hc(1:N_CS))) > 1.e-4) THEN
1004	WRITE(abort_message, *) 'cc_hc_mesh < max', cc_hc_mesh, maxval(test_hc(1:N_CS))
1005	CALL abort_physic(modname, abort_message, 1)
1006	endif
1007
1008	IF (cc_hist_mesh > tsom_hist .AND. &
1009	abs(cc_hist_mesh - tsom_hist) > 1.e-4) THEN
1010	WRITE(abort_message, *) 'cc_hist_mesh > tsom_hist', cc_hist_mesh, tsom_hist
1011	CALL abort_physic(modname, abort_message, 1)
1012	endif
1013
1014	IF (cc_hist_mesh < 0.) THEN
1015	WRITE(abort_message, *) 'cc_hist_mesh < 0', cc_hist_mesh
1016	CALL abort_physic(modname, abort_message, 1)
1017	endif
1018
1019	IF ((pcld_hc_mesh > maxval(test_pcld(1:N_CS)) .OR. &
1020	pcld_hc_mesh < minval(test_pcld(1:N_CS))) .AND. &
1021	abs(pcld_hc_mesh - maxval(test_pcld(1:N_CS))) > 1. .AND. &
1022	maxval(test_pcld(1:N_CS)) /= 999. &
1023	.AND. minval(test_pcld(1:N_CS)) /= 999.) THEN
1024	WRITE(abort_message, *) 'pcld_hc_mesh est faux', pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
1025	minval(test_pcld(1:N_CS))
1026	CALL abort_physic(modname, abort_message, 1)
1027	endif
1028
1029	IF ((tcld_hc_mesh > maxval(test_tcld(1:N_CS)) .OR. &
1030	tcld_hc_mesh < minval(test_tcld(1:N_CS))) .AND. &
1031	abs(tcld_hc_mesh - maxval(test_tcld(1:N_CS))) > 0.1 .AND. &
1032	maxval(test_tcld(1:N_CS)) /= 999. &
1033	.AND. minval(test_tcld(1:N_CS)) /= 999.) THEN
1034	WRITE(abort_message, *) 'tcld_hc_mesh est faux', tcld_hc_mesh, maxval(test_tcld(1:N_CS)), &
1035	minval(test_tcld(1:N_CS))
1036	CALL abort_physic(modname, abort_message, 1)
1037	endif
1038
1039	IF ((em_hc_mesh > maxval(test_em(1:N_CS)) .OR. &
1040	em_hc_mesh < minval(test_em(1:N_CS))) .AND. &
1041	abs(em_hc_mesh - maxval(test_em(1:N_CS))) > 1.e-4 .AND. &
1042	minval(test_em(1:N_CS)) /= 999. .AND. &
1043	maxval(test_em(1:N_CS)) /= 999.) THEN
1044	WRITE(abort_message, *) 'em_hc_mesh est faux', em_hc_mesh, maxval(test_em(1:N_CS)), &
1045	minval(test_em(1:N_CS))
1046	CALL abort_physic(modname, abort_message, 1)
1047	endif
1048
1049	else ! if N_CS>0
1050
1051	cc_tot_mesh = 0.
1052	cc_hc_mesh = 0.
1053	cc_hist_mesh = 0.
1054
1055	cc_Cb_mesh = 0.
1056	cc_ThCi_mesh = 0.
1057	cc_Anv_mesh = 0.
1058
1059	iwp_hc_mesh = undef
1060	deltaz_hc_mesh = undef
1061	em_hc_mesh = undef
1062	iwp_hist_mesh = undef
1063	deltaz_hist_mesh = undef
1064	rad_hist_mesh = undef
1065	em_hist_mesh = undef
1066	pcld_hc_mesh = undef
1067	tcld_hc_mesh = undef
1068
1069	pcld_Cb_mesh = undef
1070	tcld_Cb_mesh = undef
1071	em_Cb_mesh = undef
1072
1073	pcld_ThCi_mesh = undef
1074	tcld_ThCi_mesh = undef
1075	em_ThCi_mesh = undef
1076
1077	pcld_Anv_mesh = undef
1078	tcld_Anv_mesh = undef
1079	em_Anv_mesh = undef
1080
1081	endif ! if N_CS>0
1082
1083	END SUBROUTINE sim_mesh
1084
1085	SUBROUTINE test_bornes(sx, x, bsup, binf)
1086
1087	REAL, INTENT(IN) :: x, bsup, binf
1088	CHARACTER*14, INTENT(IN) :: sx
1089	CHARACTER (len = 50) :: modname = 'simu_airs.test_bornes'
1090	CHARACTER (len = 160) :: abort_message
1091
1092	IF (x > bsup .OR. x < binf) THEN
1093	WRITE(abort_message, *) sx, 'est faux', sx, x
1094	CALL abort_physic(modname, abort_message, 1)
1095	endif
1096
1097	END SUBROUTINE test_bornes
1098
1099
1100	SUBROUTINE simu_airs &
1101	(itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, &
1102	geop_airs, pplay_airs, paprs_airs, &
1103	map_prop_hc, map_prop_hist, &
1104	map_emis_hc, map_iwp_hc, map_deltaz_hc, map_pcld_hc, map_tcld_hc, &
1105	map_emis_Cb, map_pcld_Cb, map_tcld_Cb, &
1106	map_emis_ThCi, map_pcld_ThCi, map_tcld_ThCi, &
1107	map_emis_Anv, map_pcld_Anv, map_tcld_Anv, &
1108	map_emis_hist, map_iwp_hist, map_deltaz_hist, map_rad_hist, &
1109	map_ntot, map_hc, map_hist, &
1110	map_Cb, map_ThCi, map_Anv, alt_tropo)
1111
1112	USE dimphy
1113	USE lmdz_yomcst
1114
1115	IMPLICIT NONE
1116
1117	INTEGER, INTENT(IN) :: itap
1118
1119	REAL, DIMENSION(klon, klev), INTENT(IN) :: &
1120	rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
1121	rad_airs, geop_airs, pplay_airs, paprs_airs
1122
1123	REAL, DIMENSION(klon, klev) :: &
1124	rhodz_airs, rho_airs, iwcon_airs
1125
1126	REAL, DIMENSION(klon), INTENT(OUT) :: alt_tropo
1127
1128	REAL, DIMENSION(klev) :: rneb_1D, temp_1D, &
1129	emis_1D, rad_1D, pres_1D, alt_1D, &
1130	rhodz_1D, dz_1D, iwcon_1D
1131
1132	INTEGER :: i, j
1133
1134	REAL :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
1135	REAL :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
1136	REAL :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
1137	REAL :: em_hist_mesh, iwp_hist_mesh
1138	REAL :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh
1139	REAL :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
1140	REAL :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh
1141	REAL :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
1142
1143	REAL, DIMENSION(klon), INTENT(OUT) :: map_prop_hc, map_prop_hist
1144	REAL, DIMENSION(klon), INTENT(OUT) :: map_emis_hc, map_iwp_hc
1145	REAL, DIMENSION(klon), INTENT(OUT) :: map_deltaz_hc, map_pcld_hc
1146	REAL, DIMENSION(klon), INTENT(OUT) :: map_tcld_hc
1147	REAL, DIMENSION(klon), INTENT(OUT) :: map_emis_Cb, map_pcld_Cb, map_tcld_Cb
1148	REAL, DIMENSION(klon), INTENT(OUT) :: &
1149	map_emis_ThCi, map_pcld_ThCi, map_tcld_ThCi
1150	REAL, DIMENSION(klon), INTENT(OUT) :: &
1151	map_emis_Anv, map_pcld_Anv, map_tcld_Anv
1152	REAL, DIMENSION(klon), INTENT(OUT) :: &
1153	map_emis_hist, map_iwp_hist, map_deltaz_hist, &
1154	map_rad_hist
1155	REAL, DIMENSION(klon), INTENT(OUT) :: map_ntot, map_hc, map_hist
1156	REAL, DIMENSION(klon), INTENT(OUT) :: map_Cb, map_ThCi, map_Anv
1157
1158	WRITE(, ) 'simu_airs'
1159	WRITE(, ) 'itap, klon, klev', itap, klon, klev
1160	WRITE(, ) 'RG, RD =', RG, RD
1161
1162
1163	! Definition des variables 1D
1164
1165	do i = 1, klon
1166	do j = 1, klev - 1
1167	rhodz_airs(i, j) = &
1168	(paprs_airs(i, j) - paprs_airs(i, j + 1)) / RG
1169	enddo
1170	rhodz_airs(i, klev) = 0.
1171	enddo
1172
1173	do i = 1, klon
1174	do j = 1, klev
1175	rho_airs(i, j) = &
1176	pplay_airs(i, j) / (temp_airs(i, j) * RD)
1177
1178	IF (rneb_airs(i, j) > 0.001) THEN
1179	iwcon_airs(i, j) = iwcon0_airs(i, j) / rneb_airs(i, j)
1180	else
1181	iwcon_airs(i, j) = 0.
1182	endif
1183
1184	enddo
1185	enddo
1186
1187	!=============================================================================
1188
1189	do i = 1, klon ! boucle sur les points de grille
1190
1191	!=============================================================================
1192
1193	do j = 1, klev
1194
1195	rneb_1D(j) = rneb_airs(i, j)
1196	temp_1D(j) = temp_airs(i, j)
1197	emis_1D(j) = cldemi_airs(i, j)
1198	iwcon_1D(j) = iwcon_airs(i, j)
1199	rad_1D(j) = rad_airs(i, j)
1200	pres_1D(j) = pplay_airs(i, j)
1201	alt_1D(j) = geop_airs(i, j) / RG
1202	rhodz_1D(j) = rhodz_airs(i, j)
1203	dz_1D(j) = rhodz_airs(i, j) / rho_airs(i, j)
1204
1205	enddo
1206
1207	alt_tropo(i) = &
1208	search_tropopause(pres_1D / 100., temp_1D, alt_1D, klev)
1209
1210
1211	! Appel du ss-programme sim_mesh
1212
1213	! if (itap .EQ. 1 ) THEN
1214	CALL sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, &
1215	pres_1D, dz_1D, rhodz_1D, &
1216	cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, &
1217	pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, &
1218	deltaz_hc_mesh, &
1219	cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
1220	pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
1221	pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
1222	pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
1223	em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
1224
1225	WRITE(, ) '===================================='
1226	WRITE(, ) 'itap, i:', itap, i
1227	WRITE(, ) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, &
1228	iwp_hc, em_hist, iwp_hist ='
1229	WRITE(, ) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
1230	WRITE(, ) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
1231	WRITE(, ) em_hist_mesh, iwp_hist_mesh
1232
1233	! endif
1234
1235	! Definition des variables a ecrire dans le fichier de sortie
1236
1237	CALL normal2_undef(map_prop_hc(i), cc_hc_mesh, &
1238	cc_tot_mesh)
1239	CALL normal2_undef(map_prop_hist(i), cc_hist_mesh, &
1240	cc_tot_mesh)
1241
1242	map_emis_hc(i) = em_hc_mesh
1243	map_iwp_hc(i) = iwp_hc_mesh
1244	map_deltaz_hc(i) = deltaz_hc_mesh
1245	map_pcld_hc(i) = pcld_hc_mesh
1246	map_tcld_hc(i) = tcld_hc_mesh
1247
1248	map_emis_Cb(i) = em_Cb_mesh
1249	map_pcld_Cb(i) = pcld_Cb_mesh
1250	map_tcld_Cb(i) = tcld_Cb_mesh
1251
1252	map_emis_ThCi(i) = em_ThCi_mesh
1253	map_pcld_ThCi(i) = pcld_ThCi_mesh
1254	map_tcld_ThCi(i) = tcld_ThCi_mesh
1255
1256	map_emis_Anv(i) = em_Anv_mesh
1257	map_pcld_Anv(i) = pcld_Anv_mesh
1258	map_tcld_Anv(i) = tcld_Anv_mesh
1259
1260	map_emis_hist(i) = em_hist_mesh
1261	map_iwp_hist(i) = iwp_hist_mesh
1262	map_deltaz_hist(i) = deltaz_hist_mesh
1263	map_rad_hist(i) = rad_hist_mesh
1264
1265	map_ntot(i) = cc_tot_mesh
1266	map_hc(i) = cc_hc_mesh
1267	map_hist(i) = cc_hist_mesh
1268
1269	map_Cb(i) = cc_Cb_mesh
1270	map_ThCi(i) = cc_ThCi_mesh
1271	map_Anv(i) = cc_Anv_mesh
1272
1273	enddo ! fin boucle sur les points de grille
1274
1275	END SUBROUTINE simu_airs
1276
1277
1278	END MODULE lmdz_simu_airs
1279

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

Download in other formats:

Original Format