Context Navigation

simu_airs.F90 @ 5006

Last change on this file since 5006 was 3531, checked in by Laurent Fairhead, 6 years ago
Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86 Still some work to be done in phylmd subdirectories
File size: 36.7 KB

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

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

Download in other formats:

Original Format