Context Navigation

simu_airs.F90 @ 5112

Last change on this file since 5112 was 5112, checked in by abarral, 3 months ago
Rename modules in phy_common from _mod > lmdz_
File size: 35.6 KB

Line
1	module m_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
407	nuage = nuage + 1
408
409	! On determine les caracteristiques du nuage
410	! (ep. optique, ice water path, pression, temperature)
411
412	CALL caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
413	pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
414	som_tau, som_iwc, som_dz, som_rad)
415
416	! On masque le nuage s'il n'est pas detectable
417
418	CALL masque(ibeg, iend, som_tau, visible, w)
419
420	endif
421
422	! Fin boucle nuages
423	enddo
424
425
426	! Analyse du nuage detecte
427
428	CALL topbot(1, w, ibeg, iend)
429
430	if (ibeg > 0) then
431
432	CALL caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
433	pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
434	som_tau, som_iwc, som_dz, som_rad)
435
436	tau = som_tau
437	emis = 1. - exp(-tau)
438	iwp = som_iwc
439	deltaz = som_dz
440	rad = som_rad
441
442	if (pcld > p_thresh) then
443
444	shc = 0.
445	shist = 0.
446	sCb = 0.
447	sThCi = 0.
448	sAnv = 0.
449
450	else
451
452	if (emis < em_min .or. emis > em_max &
453	.or. tcld > 230.) then
454	shist = 0.
455	endif
456
457	if (emis < 0.98) then
458	sCb = 0.
459	endif
460
461	if (emis > 0.5 .or. emis < 0.1) then
462	sThCi = 0.
463	endif
464
465	if (emis <= 0.5 .or. emis >= 0.98) then
466	sAnv = 0.
467	endif
468
469	endif
470
471	else
472
473	tau = 0.
474	emis = 0.
475	iwp = 0.
476	deltaz = 0.
477	pcld = 0.
478	tcld = 0.
479	stot = 0.
480	shc = 0.
481	shist = 0.
482	rad = 0.
483	sCb = 0.
484	sThCi = 0.
485	sAnv = 0.
486
487	endif
488
489
490	! Tests
491
492	if (iwp < 0.) then
493	WRITE(abort_message, *) 'ideb iwp =', ideb, iwp
494	CALL abort_physic(modname, abort_message, 1)
495	endif
496
497	if (deltaz < 0.) then
498	WRITE(abort_message, *)'ideb deltaz =', ideb, deltaz
499	CALL abort_physic(modname, abort_message, 1)
500	endif
501
502	if (emis < 0.048 .and. emis /= 0.) then
503	WRITE(abort_message, *) 'ideb emis =', ideb, emis
504	CALL abort_physic(modname, abort_message, 1)
505	endif
506
507	END SUBROUTINE sous_section
508
509
510	SUBROUTINE masque(ibeg, iend, som_tau, &
511	visible, w)
512
513	INTEGER, intent(in) :: ibeg, iend
514	REAL, intent(in) :: som_tau
515
516	INTEGER, intent(inout) :: visible
517	REAL, DIMENSION(:), intent(inout) :: w
518
519	INTEGER :: i
520
521
522
523	! Masque
524
525	! Cas ou il n'y a pas de nuage visible au-dessus
526
527	if (visible == 0) then
528
529	if (som_tau < tau_thresh) then
530	do i = ibeg, iend
531	w(i) = 0.
532	enddo
533	else
534	visible = 1
535	endif
536
537	! Cas ou il y a un nuage visible au-dessus
538
539	else
540
541	do i = ibeg, iend
542	w(i) = 0.
543	enddo
544
545	endif
546
547	END SUBROUTINE masque
548
549
550	SUBROUTINE caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
551	pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
552	som_tau, som_iwc, som_dz, som_rad)
553
554	INTEGER, intent(in) :: ibeg, iend
555	REAL, DIMENSION(:), intent(in) :: tau_cs, iwco_cs, temp_cs
556	REAL, DIMENSION(:), intent(in) :: pres_cs, dz_cs, rad_cs
557	REAL, DIMENSION(:), intent(in) :: rhodz_cs
558	REAL, intent(out) :: som_tau, som_iwc, som_dz, som_rad
559	REAL, intent(out) :: pcld, tcld
560
561	INTEGER :: i, ibase, imid
562
563	CHARACTER (len = 50) :: modname = 'simu_airs.caract'
564	CHARACTER (len = 160) :: abort_message
565
566	! Somme des epaisseurs optiques et des contenus en glace sur le nuage
567
568	som_tau = 0.
569	som_iwc = 0.
570	som_dz = 0.
571	som_rad = 0.
572	ibase = -100
573
574	do i = ibeg, iend
575
576	som_tau = som_tau + tau_cs(i)
577
578	som_dz = som_dz + dz_cs(i)
579	som_iwc = som_iwc + iwco_cs(i) * 1000 * rhodz_cs(i) ! en g/m2
580	som_rad = som_rad + rad_cs(i) * dz_cs(i)
581
582	if (som_tau > 3. .and. ibase == -100) then
583	ibase = i - 1
584	endif
585
586	enddo
587
588	if (som_dz /= 0.) then
589	som_rad = som_rad / som_dz
590	else
591	write(, ) 'som_dez = 0 stop '
592	write(, ) 'ibeg, iend =', ibeg, iend
593	do i = ibeg, iend
594	write(, ) dz_cs(i), rhodz_cs(i)
595	enddo
596	abort_message = 'see above'
597	CALL abort_physic(modname, abort_message, 1)
598	endif
599
600	! Determination de Pcld et Tcld
601
602	if (ibase < ibeg) then
603	ibase = ibeg
604	endif
605
606	imid = (ibeg + ibase) / 2
607
608	pcld = pres_cs(imid) / 100. ! pcld en hPa
609	tcld = temp_cs(imid)
610
611	END SUBROUTINE caract
612
613	SUBROUTINE topbot(ideb, w, ibeg, iend)
614
615	INTEGER, intent(in) :: ideb
616	REAL, DIMENSION(:), intent(in) :: w
617	INTEGER, intent(out) :: ibeg, iend
618
619	INTEGER :: i, itest
620
621	itest = 0
622	ibeg = 0
623	iend = 0
624
625	do i = ideb, size(w)
626
627	if (w(i) == 1. .and. itest == 0) then
628	ibeg = i
629	itest = 1
630	endif
631
632	enddo
633
634	i = ibeg
635	do while (w(i) == 1. .and. i <= size(w))
636	i = i + 1
637	enddo
638	iend = i - 1
639
640	END SUBROUTINE topbot
641
642	SUBROUTINE ordonne(len_cs, rneb_cs, rneb_ord)
643
644	INTEGER, intent(in) :: len_cs
645	REAL, DIMENSION(:), intent(in) :: rneb_cs
646	REAL, DIMENSION(:), intent(out) :: rneb_ord
647
648	INTEGER :: i, j, ind_min
649
650	REAL, DIMENSION(len_cs) :: rneb
651	REAL :: rneb_max
652
653	do i = 1, size(rneb_cs)
654	rneb(i) = rneb_cs(i)
655	enddo
656
657	do j = 1, size(rneb_cs)
658
659	rneb_max = 100.
660
661	do i = 1, size(rneb_cs)
662	if (rneb(i) < rneb_max) then
663	rneb_max = rneb(i)
664	ind_min = i
665	endif
666	enddo
667
668	rneb_ord(j) = rneb_max
669	rneb(ind_min) = 100.
670
671	enddo
672
673	END SUBROUTINE ordonne
674
675
676	SUBROUTINE sim_mesh(rneb_1D, temp_1D, emis_1D, &
677	iwcon_1D, rad_1D, &
678	pres, dz, &
679	rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh, &
680	tcld_hc_mesh, &
681	em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, &
682	cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
683	pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
684	pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
685	pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
686	em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
687
688	USE dimphy
689
690	REAL, DIMENSION(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, &
691	iwcon_1D, rad_1D
692	REAL, DIMENSION(klev), intent(in) :: pres, dz, rhodz_1D
693	REAL, intent(out) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
694	REAL, intent(out) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
695	REAL, intent(out) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, &
696	iwp_hc_mesh
697
698	REAL, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
699	REAL, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
700	em_ThCi_mesh
701	REAL, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
702
703	REAL, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
704	REAL, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh
705
706	REAL, DIMENSION(:), allocatable :: rneb_cs, temp_cs, emis_cs, &
707	iwco_cs
708	REAL, DIMENSION(:), allocatable :: pres_cs, dz_cs, rad_cs, &
709	rhodz_cs
710
711	INTEGER :: i, j, l
712	INTEGER :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics
713
714	REAL :: som_emi_hc, som_pcl_hc, som_tcl_hc, som_iwp_hc, som_hc, &
715	som_hist
716	REAL :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
717	som_deltaz_hist
718	REAL :: som_rad_hist
719	REAL :: som_Cb, som_ThCi, som_Anv
720	REAL :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb
721	REAL :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv
722	REAL :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi
723	REAL :: tsom_tot, tsom_hc, tsom_hist
724	REAL :: prod, prod_hh
725
726	REAL :: cc_tot_cs, cc_hc_cs, cc_hist_cs
727	REAL :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs
728	REAL :: pcld_hc_cs, tcld_hc_cs
729	REAL :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs
730	REAL :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs
731	REAL :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs
732	REAL :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs
733	REAL :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs
734
735	REAL, DIMENSION(klev) :: test_tot, test_hc, test_hist
736	REAL, DIMENSION(klev) :: test_pcld, test_tcld, test_em, test_iwp
737
738	CHARACTER (len = 50) :: modname = 'simu_airs.sim_mesh'
739	CHARACTER (len = 160) :: abort_message
740
741	do j = 1, klev
742	WRITE(lunout, *) 'simu_airs, j, rneb_1D =', rneb_1D(j)
743	enddo
744
745	! Definition des structures nuageuses, de la plus haute a la plus basse
746
747	num_cs = 0
748	ltop = klev - 1
749
750	prod = 1.
751
752	som_emi_hc = 0.
753	som_emi_hist = 0.
754	som_pcl_hc = 0.
755	som_tcl_hc = 0.
756	som_iwp_hc = 0.
757	som_iwp_hist = 0.
758	som_deltaz_hc = 0.
759	som_deltaz_hist = 0.
760	som_rad_hist = 0.
761	som_hc = 0.
762	som_hist = 0.
763
764	som_Cb = 0.
765	som_ThCi = 0.
766	som_Anv = 0.
767
768	som_emi_Cb = 0.
769	som_pcld_Cb = 0.
770	som_tcld_Cb = 0.
771
772	som_emi_ThCi = 0.
773	som_pcld_ThCi = 0.
774	som_tcld_ThCi = 0.
775
776	som_emi_Anv = 0.
777	som_pcld_Anv = 0.
778	som_tcld_Anv = 0.
779
780	tsom_tot = 0.
781	tsom_hc = 0.
782	tsom_hist = 0.
783
784	do while (ltop >= 1) ! Boucle sur toute la colonne
785
786	itop = 0
787
788	do l = ltop, 1, -1
789
790	if (itop == 0 .and. rneb_1D(l) > 0.001) then
791	itop = l
792	endif
793
794	enddo
795
796	ibot = itop
797
798	do while (rneb_1D(ibot) > 0.001 .and. ibot >= 1)
799	ibot = ibot - 1
800	enddo
801
802	ibot = ibot + 1
803
804	if (itop > 0) then ! itop > 0
805
806	num_cs = num_cs + 1
807	len_cs = itop - ibot + 1
808
809	! Allocation et definition des variables de la structure nuageuse
810	! le premier indice denote ce qui est le plus haut
811
812	allocate (rneb_cs(len_cs))
813	allocate (temp_cs(len_cs))
814	allocate (emis_cs(len_cs))
815	allocate (iwco_cs(len_cs))
816	allocate (pres_cs(len_cs))
817	allocate (dz_cs(len_cs))
818	allocate (rad_cs(len_cs))
819	allocate (rhodz_cs(len_cs))
820
821	ics = 0
822
823	do i = itop, ibot, -1
824	ics = ics + 1
825	rneb_cs(ics) = rneb_1D(i)
826	temp_cs(ics) = temp_1D(i)
827	emis_cs(ics) = emis_1D(i)
828	iwco_cs(ics) = iwcon_1D(i)
829	rad_cs(ics) = rad_1D(i)
830	pres_cs(ics) = pres(i)
831	dz_cs(ics) = dz(i)
832	rhodz_cs(ics) = rhodz_1D(i)
833	enddo
834
835	! Appel du sous_programme cloud_structure
836
837	CALL cloud_structure(len_cs, rneb_cs, temp_cs, emis_cs, iwco_cs, &
838	pres_cs, dz_cs, rhodz_cs, rad_cs, &
839	cc_tot_cs, cc_hc_cs, cc_hist_cs, &
840	cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
841	pcld_hc_cs, tcld_hc_cs, &
842	em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
843	pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
844	pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
845	pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
846	em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
847
848	deallocate (rneb_cs)
849	deallocate (temp_cs)
850	deallocate (emis_cs)
851	deallocate (iwco_cs)
852	deallocate (pres_cs)
853	deallocate (dz_cs)
854	deallocate (rad_cs)
855	deallocate (rhodz_cs)
856
857
858	! Pour la couverture nuageuse sur la maille
859
860	prod_hh = prod
861
862	prod = prod * (1. - cc_tot_cs)
863
864	! Pour les autres variables definies sur la maille
865
866	som_emi_hc = som_emi_hc + em_hc_cs * cc_hc_cs * prod_hh
867	som_iwp_hc = som_iwp_hc + iwp_hc_cs * cc_hc_cs * prod_hh
868	som_deltaz_hc = som_deltaz_hc + deltaz_hc_cs * cc_hc_cs * prod_hh
869
870	som_emi_Cb = som_emi_Cb + em_Cb_cs * cc_Cb_cs * prod_hh
871	som_tcld_Cb = som_tcld_Cb + tcld_Cb_cs * cc_Cb_cs * prod_hh
872	som_pcld_Cb = som_pcld_Cb + pcld_Cb_cs * cc_Cb_cs * prod_hh
873
874	som_emi_ThCi = som_emi_ThCi + em_ThCi_cs * cc_ThCi_cs * prod_hh
875	som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cs * cc_ThCi_cs * prod_hh
876	som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cs * cc_ThCi_cs * prod_hh
877
878	som_emi_Anv = som_emi_Anv + em_Anv_cs * cc_Anv_cs * prod_hh
879	som_tcld_Anv = som_tcld_Anv + tcld_Anv_cs * cc_Anv_cs * prod_hh
880	som_pcld_Anv = som_pcld_Anv + pcld_Anv_cs * cc_Anv_cs * prod_hh
881
882	som_emi_hist = som_emi_hist + em_hist_cs * cc_hist_cs * prod_hh
883	som_iwp_hist = som_iwp_hist + iwp_hist_cs * cc_hist_cs * prod_hh
884	som_deltaz_hist = som_deltaz_hist + &
885	deltaz_hist_cs * cc_hist_cs * prod_hh
886	som_rad_hist = som_rad_hist + rad_hist_cs * cc_hist_cs * prod_hh
887
888	som_pcl_hc = som_pcl_hc + pcld_hc_cs * cc_hc_cs * prod_hh
889	som_tcl_hc = som_tcl_hc + tcld_hc_cs * cc_hc_cs * prod_hh
890
891	som_hc = som_hc + cc_hc_cs * prod_hh
892	som_hist = som_hist + cc_hist_cs * prod_hh
893
894	som_Cb = som_Cb + cc_Cb_cs * prod_hh
895	som_ThCi = som_ThCi + cc_ThCi_cs * prod_hh
896	som_Anv = som_Anv + cc_Anv_cs * prod_hh
897
898
899	! Pour test
900
901	CALL test_bornes('cc_tot_cs ', cc_tot_cs, 1., 0.)
902	CALL test_bornes('cc_hc_cs ', cc_hc_cs, 1., 0.)
903	CALL test_bornes('cc_hist_cs ', cc_hist_cs, 1., 0.)
904	CALL test_bornes('pcld_hc_cs ', pcld_hc_cs, 1200., 0.)
905	CALL test_bornes('tcld_hc_cs ', tcld_hc_cs, 1000., 100.)
906	CALL test_bornes('em_hc_cs ', em_hc_cs, 1000., 0.048)
907
908	test_tot(num_cs) = cc_tot_cs
909	test_hc(num_cs) = cc_hc_cs
910	test_hist(num_cs) = cc_hist_cs
911	test_pcld(num_cs) = pcld_hc_cs
912	test_tcld(num_cs) = tcld_hc_cs
913	test_em(num_cs) = em_hc_cs
914	test_iwp(num_cs) = iwp_hc_cs
915
916	tsom_tot = tsom_tot + cc_tot_cs
917	tsom_hc = tsom_hc + cc_hc_cs
918	tsom_hist = tsom_hist + cc_hist_cs
919
920	endif ! itop > 0
921
922	ltop = ibot - 1
923
924	enddo ! fin de la boucle sur la colonne
925
926	N_CS = num_cs
927
928
929	! Determination des variables de sortie
930
931	if (N_CS > 0) then ! if N_CS>0
932
933	cc_tot_mesh = 1. - prod
934
935	cc_hc_mesh = som_hc
936	cc_hist_mesh = som_hist
937
938	cc_Cb_mesh = som_Cb
939	cc_ThCi_mesh = som_ThCi
940	cc_Anv_mesh = som_Anv
941
942	CALL normal2_undef(pcld_hc_mesh, som_pcl_hc, &
943	cc_hc_mesh)
944	CALL normal2_undef(tcld_hc_mesh, som_tcl_hc, &
945	cc_hc_mesh)
946	CALL normal2_undef(em_hc_mesh, som_emi_hc, &
947	cc_hc_mesh)
948	CALL normal2_undef(iwp_hc_mesh, som_iwp_hc, &
949	cc_hc_mesh)
950	CALL normal2_undef(deltaz_hc_mesh, som_deltaz_hc, &
951	cc_hc_mesh)
952
953	CALL normal2_undef(em_Cb_mesh, som_emi_Cb, &
954	cc_Cb_mesh)
955	CALL normal2_undef(tcld_Cb_mesh, som_tcld_Cb, &
956	cc_Cb_mesh)
957	CALL normal2_undef(pcld_Cb_mesh, som_pcld_Cb, &
958	cc_Cb_mesh)
959
960	CALL normal2_undef(em_ThCi_mesh, som_emi_ThCi, &
961	cc_ThCi_mesh)
962	CALL normal2_undef(tcld_ThCi_mesh, som_tcld_ThCi, &
963	cc_ThCi_mesh)
964	CALL normal2_undef(pcld_ThCi_mesh, som_pcld_ThCi, &
965	cc_ThCi_mesh)
966
967	CALL normal2_undef(em_Anv_mesh, som_emi_Anv, &
968	cc_Anv_mesh)
969	CALL normal2_undef(tcld_Anv_mesh, som_tcld_Anv, &
970	cc_Anv_mesh)
971	CALL normal2_undef(pcld_Anv_mesh, som_pcld_Anv, &
972	cc_Anv_mesh)
973
974	CALL normal2_undef(em_hist_mesh, som_emi_hist, &
975	cc_hist_mesh)
976	CALL normal2_undef(iwp_hist_mesh, som_iwp_hist, &
977	cc_hist_mesh)
978	CALL normal2_undef(deltaz_hist_mesh, som_deltaz_hist, &
979	cc_hist_mesh)
980	CALL normal2_undef(rad_hist_mesh, som_rad_hist, &
981	cc_hist_mesh)
982
983
984	! Tests
985
986	! Tests
987
988	if (cc_tot_mesh > tsom_tot .and. &
989	abs(cc_tot_mesh - tsom_tot) > 1.e-4) then
990	WRITE(abort_message, *)'cc_tot_mesh > tsom_tot', cc_tot_mesh, tsom_tot
991	CALL abort_physic(modname, abort_message, 1)
992	endif
993
994	if (cc_tot_mesh < maxval(test_tot(1:N_CS)) .and. &
995	abs(cc_tot_mesh - maxval(test_tot(1:N_CS))) > 1.e-4) then
996	WRITE(abort_message, *) 'cc_tot_mesh < max', cc_tot_mesh, maxval(test_tot(1:N_CS))
997	CALL abort_physic(modname, abort_message, 1)
998	endif
999
1000	if (cc_hc_mesh > tsom_hc .and. &
1001	abs(cc_hc_mesh - tsom_hc) > 1.e-4) then
1002	WRITE(abort_message, *) 'cc_hc_mesh > tsom_hc', cc_hc_mesh, tsom_hc
1003	CALL abort_physic(modname, abort_message, 1)
1004	endif
1005
1006	if (cc_hc_mesh < maxval(test_hc(1:N_CS)) .and. &
1007	abs(cc_hc_mesh - maxval(test_hc(1:N_CS))) > 1.e-4) then
1008	WRITE(abort_message, *) 'cc_hc_mesh < max', cc_hc_mesh, maxval(test_hc(1:N_CS))
1009	CALL abort_physic(modname, abort_message, 1)
1010	endif
1011
1012	if (cc_hist_mesh > tsom_hist .and. &
1013	abs(cc_hist_mesh - tsom_hist) > 1.e-4) then
1014	WRITE(abort_message, *) 'cc_hist_mesh > tsom_hist', cc_hist_mesh, tsom_hist
1015	CALL abort_physic(modname, abort_message, 1)
1016	endif
1017
1018	if (cc_hist_mesh < 0.) then
1019	WRITE(abort_message, *) 'cc_hist_mesh < 0', cc_hist_mesh
1020	CALL abort_physic(modname, abort_message, 1)
1021	endif
1022
1023	if ((pcld_hc_mesh > maxval(test_pcld(1:N_CS)) .or. &
1024	pcld_hc_mesh < minval(test_pcld(1:N_CS))) .and. &
1025	abs(pcld_hc_mesh - maxval(test_pcld(1:N_CS))) > 1. .and. &
1026	maxval(test_pcld(1:N_CS)) /= 999. &
1027	.and. minval(test_pcld(1:N_CS)) /= 999.) then
1028	WRITE(abort_message, *) 'pcld_hc_mesh est faux', pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
1029	minval(test_pcld(1:N_CS))
1030	CALL abort_physic(modname, abort_message, 1)
1031	endif
1032
1033	if ((tcld_hc_mesh > maxval(test_tcld(1:N_CS)) .or. &
1034	tcld_hc_mesh < minval(test_tcld(1:N_CS))) .and. &
1035	abs(tcld_hc_mesh - maxval(test_tcld(1:N_CS))) > 0.1 .and. &
1036	maxval(test_tcld(1:N_CS)) /= 999. &
1037	.and. minval(test_tcld(1:N_CS)) /= 999.) then
1038	WRITE(abort_message, *) 'tcld_hc_mesh est faux', tcld_hc_mesh, maxval(test_tcld(1:N_CS)), &
1039	minval(test_tcld(1:N_CS))
1040	CALL abort_physic(modname, abort_message, 1)
1041	endif
1042
1043	if ((em_hc_mesh > maxval(test_em(1:N_CS)) .or. &
1044	em_hc_mesh < minval(test_em(1:N_CS))) .and. &
1045	abs(em_hc_mesh - maxval(test_em(1:N_CS))) > 1.e-4 .and. &
1046	minval(test_em(1:N_CS)) /= 999. .and. &
1047	maxval(test_em(1:N_CS)) /= 999.) then
1048	WRITE(abort_message, *) 'em_hc_mesh est faux', em_hc_mesh, maxval(test_em(1:N_CS)), &
1049	minval(test_em(1:N_CS))
1050	CALL abort_physic(modname, abort_message, 1)
1051	endif
1052
1053	else ! if N_CS>0
1054
1055	cc_tot_mesh = 0.
1056	cc_hc_mesh = 0.
1057	cc_hist_mesh = 0.
1058
1059	cc_Cb_mesh = 0.
1060	cc_ThCi_mesh = 0.
1061	cc_Anv_mesh = 0.
1062
1063	iwp_hc_mesh = undef
1064	deltaz_hc_mesh = undef
1065	em_hc_mesh = undef
1066	iwp_hist_mesh = undef
1067	deltaz_hist_mesh = undef
1068	rad_hist_mesh = undef
1069	em_hist_mesh = undef
1070	pcld_hc_mesh = undef
1071	tcld_hc_mesh = undef
1072
1073	pcld_Cb_mesh = undef
1074	tcld_Cb_mesh = undef
1075	em_Cb_mesh = undef
1076
1077	pcld_ThCi_mesh = undef
1078	tcld_ThCi_mesh = undef
1079	em_ThCi_mesh = undef
1080
1081	pcld_Anv_mesh = undef
1082	tcld_Anv_mesh = undef
1083	em_Anv_mesh = undef
1084
1085	endif ! if N_CS>0
1086
1087	END SUBROUTINE sim_mesh
1088
1089	SUBROUTINE test_bornes(sx, x, bsup, binf)
1090
1091	REAL, intent(in) :: x, bsup, binf
1092	character*14, intent(in) :: sx
1093	CHARACTER (len = 50) :: modname = 'simu_airs.test_bornes'
1094	CHARACTER (len = 160) :: abort_message
1095
1096	if (x > bsup .or. x < binf) then
1097	WRITE(abort_message, *) sx, 'est faux', sx, x
1098	CALL abort_physic(modname, abort_message, 1)
1099	endif
1100
1101	END SUBROUTINE test_bornes
1102
1103	end module m_simu_airs
1104
1105
1106	SUBROUTINE simu_airs &
1107	(itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, &
1108	geop_airs, pplay_airs, paprs_airs, &
1109	map_prop_hc, map_prop_hist, &
1110	map_emis_hc, map_iwp_hc, map_deltaz_hc, map_pcld_hc, map_tcld_hc, &
1111	map_emis_Cb, map_pcld_Cb, map_tcld_Cb, &
1112	map_emis_ThCi, map_pcld_ThCi, map_tcld_ThCi, &
1113	map_emis_Anv, map_pcld_Anv, map_tcld_Anv, &
1114	map_emis_hist, map_iwp_hist, map_deltaz_hist, map_rad_hist, &
1115	map_ntot, map_hc, map_hist, &
1116	map_Cb, map_ThCi, map_Anv, alt_tropo)
1117
1118	USE dimphy
1119	USE m_simu_airs
1120
1121	IMPLICIT NONE
1122
1123	include "YOMCST.h"
1124
1125	INTEGER, intent(in) :: itap
1126
1127	REAL, DIMENSION(klon, klev), intent(in) :: &
1128	rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
1129	rad_airs, geop_airs, pplay_airs, paprs_airs
1130
1131	REAL, DIMENSION(klon, klev) :: &
1132	rhodz_airs, rho_airs, iwcon_airs
1133
1134	REAL, DIMENSION(klon), intent(out) :: alt_tropo
1135
1136	REAL, DIMENSION(klev) :: rneb_1D, temp_1D, &
1137	emis_1D, rad_1D, pres_1D, alt_1D, &
1138	rhodz_1D, dz_1D, iwcon_1D
1139
1140	INTEGER :: i, j
1141
1142	REAL :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
1143	REAL :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
1144	REAL :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
1145	REAL :: em_hist_mesh, iwp_hist_mesh
1146	REAL :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh
1147	REAL :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
1148	REAL :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh
1149	REAL :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
1150
1151	REAL, DIMENSION(klon), intent(out) :: map_prop_hc, map_prop_hist
1152	REAL, DIMENSION(klon), intent(out) :: map_emis_hc, map_iwp_hc
1153	REAL, DIMENSION(klon), intent(out) :: map_deltaz_hc, map_pcld_hc
1154	REAL, DIMENSION(klon), intent(out) :: map_tcld_hc
1155	REAL, DIMENSION(klon), intent(out) :: map_emis_Cb, map_pcld_Cb, map_tcld_Cb
1156	REAL, DIMENSION(klon), intent(out) :: &
1157	map_emis_ThCi, map_pcld_ThCi, map_tcld_ThCi
1158	REAL, DIMENSION(klon), intent(out) :: &
1159	map_emis_Anv, map_pcld_Anv, map_tcld_Anv
1160	REAL, DIMENSION(klon), intent(out) :: &
1161	map_emis_hist, map_iwp_hist, map_deltaz_hist, &
1162	map_rad_hist
1163	REAL, DIMENSION(klon), intent(out) :: map_ntot, map_hc, map_hist
1164	REAL, DIMENSION(klon), intent(out) :: map_Cb, map_ThCi, map_Anv
1165
1166	write(, ) 'simu_airs'
1167	write(, ) 'itap, klon, klev', itap, klon, klev
1168	write(, ) 'RG, RD =', RG, RD
1169
1170
1171	! Definition des variables 1D
1172
1173	do i = 1, klon
1174	do j = 1, klev - 1
1175	rhodz_airs(i, j) = &
1176	(paprs_airs(i, j) - paprs_airs(i, j + 1)) / RG
1177	enddo
1178	rhodz_airs(i, klev) = 0.
1179	enddo
1180
1181	do i = 1, klon
1182	do j = 1, klev
1183	rho_airs(i, j) = &
1184	pplay_airs(i, j) / (temp_airs(i, j) * RD)
1185
1186	if (rneb_airs(i, j) > 0.001) then
1187	iwcon_airs(i, j) = iwcon0_airs(i, j) / rneb_airs(i, j)
1188	else
1189	iwcon_airs(i, j) = 0.
1190	endif
1191
1192	enddo
1193	enddo
1194
1195	!=============================================================================
1196
1197	do i = 1, klon ! boucle sur les points de grille
1198
1199	!=============================================================================
1200
1201	do j = 1, klev
1202
1203	rneb_1D(j) = rneb_airs(i, j)
1204	temp_1D(j) = temp_airs(i, j)
1205	emis_1D(j) = cldemi_airs(i, j)
1206	iwcon_1D(j) = iwcon_airs(i, j)
1207	rad_1D(j) = rad_airs(i, j)
1208	pres_1D(j) = pplay_airs(i, j)
1209	alt_1D(j) = geop_airs(i, j) / RG
1210	rhodz_1D(j) = rhodz_airs(i, j)
1211	dz_1D(j) = rhodz_airs(i, j) / rho_airs(i, j)
1212
1213	enddo
1214
1215	alt_tropo(i) = &
1216	search_tropopause(pres_1D / 100., temp_1D, alt_1D, klev)
1217
1218
1219	! Appel du ss-programme sim_mesh
1220
1221	! if (itap .eq. 1 ) then
1222
1223	CALL sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, &
1224	pres_1D, dz_1D, rhodz_1D, &
1225	cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, &
1226	pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, &
1227	deltaz_hc_mesh, &
1228	cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
1229	pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
1230	pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
1231	pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
1232	em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
1233
1234	write(, ) '===================================='
1235	write(, ) 'itap, i:', itap, i
1236	write(, ) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, &
1237	iwp_hc, em_hist, iwp_hist ='
1238	write(, ) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
1239	write(, ) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
1240	write(, ) em_hist_mesh, iwp_hist_mesh
1241
1242	! endif
1243
1244	! Definition des variables a ecrire dans le fichier de sortie
1245
1246	CALL normal2_undef(map_prop_hc(i), cc_hc_mesh, &
1247	cc_tot_mesh)
1248	CALL normal2_undef(map_prop_hist(i), cc_hist_mesh, &
1249	cc_tot_mesh)
1250
1251	map_emis_hc(i) = em_hc_mesh
1252	map_iwp_hc(i) = iwp_hc_mesh
1253	map_deltaz_hc(i) = deltaz_hc_mesh
1254	map_pcld_hc(i) = pcld_hc_mesh
1255	map_tcld_hc(i) = tcld_hc_mesh
1256
1257	map_emis_Cb(i) = em_Cb_mesh
1258	map_pcld_Cb(i) = pcld_Cb_mesh
1259	map_tcld_Cb(i) = tcld_Cb_mesh
1260
1261	map_emis_ThCi(i) = em_ThCi_mesh
1262	map_pcld_ThCi(i) = pcld_ThCi_mesh
1263	map_tcld_ThCi(i) = tcld_ThCi_mesh
1264
1265	map_emis_Anv(i) = em_Anv_mesh
1266	map_pcld_Anv(i) = pcld_Anv_mesh
1267	map_tcld_Anv(i) = tcld_Anv_mesh
1268
1269	map_emis_hist(i) = em_hist_mesh
1270	map_iwp_hist(i) = iwp_hist_mesh
1271	map_deltaz_hist(i) = deltaz_hist_mesh
1272	map_rad_hist(i) = rad_hist_mesh
1273
1274	map_ntot(i) = cc_tot_mesh
1275	map_hc(i) = cc_hc_mesh
1276	map_hist(i) = cc_hist_mesh
1277
1278	map_Cb(i) = cc_Cb_mesh
1279	map_ThCi(i) = cc_ThCi_mesh
1280	map_Anv(i) = cc_Anv_mesh
1281
1282	enddo ! fin boucle sur les points de grille
1283
1284	END SUBROUTINE simu_airs
1285

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

Download in other formats:

Original Format