Context Navigation

simu_airs.F90 @ 4947

Last change on this file since 4947 was 2585, checked in by idelkadi, 8 years ago
Corrections des erreurs pour l'implementation du simulateur AIRS
File size: 35.3 KB

Rev	Line
[2580]	1
	2	module m_simu_airs
	3
	4	implicit none
	5
	6	real, parameter :: tau_thresh = 0.05 ! seuil nuages detectables
	7	real, parameter :: p_thresh = 445. ! seuil nuages hauts
	8	real, parameter :: em_min = 0.2 ! seuils nuages semi-transp
	9	real, parameter :: em_max = 0.85
	10	real, parameter :: undef = 999.
	11
	12	contains
	13
	14	real function search_tropopause(P,T,alt,N) result(P_tropo)
	15	! this function searches for the tropopause pressure in [hPa].
	16	! The search is based on ideology described in
	17	! Reichler et al., Determining the tropopause height from gridded data,
	18	! GRL, 30(20) 2042, doi:10.1029/2003GL018240, 2003
	19
	20	integer N,i,i_lev,first_point,exit_flag,i_dir
	21	real P(N),T(N),alt(N),slope(N)
	22	real P_min, P_max, slope_limit,slope_2km, &
	23	& delta_alt_limit,tmp,delta_alt
	24	parameter(P_min=75.0, P_max=470.0) ! hPa
	25	parameter(slope_limit=0.002) ! 2 K/km converted to K/m
	26	parameter(delta_alt_limit=2000.0) ! 2000 meters
	27
	28	do i=1,N-1
	29	slope(i)=-(T(i+1)-T(i))/(alt(i+1)-alt(i))
	30	end do
	31	slope(N)=slope(N-1)
	32
	33	if (P(1).gt.P(N)) then
	34	i_dir= 1
	35	i=1
	36	else
	37	i_dir=-1
	38	i=N
	39	end if
	40
	41	first_point=0
	42	exit_flag=0
	43	do while(exit_flag.eq.0)
	44	if (P(i).ge.P_min.and.P(i).le.P_max) then
	45	if (first_point.gt.0) then
	46	delta_alt=alt(i)-alt(first_point)
	47	if (delta_alt.ge.delta_alt_limit) then
	48	slope_2km=(T(first_point)-T(i))/delta_alt
	49	if (slope_2km.lt.slope_limit) then
	50	exit_flag=1
	51	else
	52	first_point=0
	53	end if
	54	end if
	55	end if
	56	if (first_point.eq.0.and.slope(i).lt.slope_limit) first_point=i
	57	end if
	58	i=i+i_dir
	59	if (i.le.1.or.i.ge.N) exit_flag=1
	60	end do
	61
	62	if (first_point.le.0) P_tropo=65.4321
	63	if (first_point.eq.1) P_tropo=64.5432
	64	if (first_point.gt.1) then
	65	tmp=(slope_limit-slope(first_point))/(slope(first_point+1)- &
	66	& slope(first_point))*(P(first_point+1)-P(first_point))
	67	P_tropo=P(first_point)+tmp
	68	! print*, 'P_tropo= ', tmp, P(first_point), P_tropo
	69	end if
	70
	71	! Ajout Marine
	72	if (P_tropo .lt. 60. .or. P_tropo .gt. 470.) then
	73	P_tropo = 999.
	74	endif
	75
	76	return
	77	end function search_tropopause
	78
	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
	94
	95	integer :: i, n, nss
	96
	97	integer, intent(in) :: len_cs
	98	real, dimension(:), intent(in) :: rneb_cs, temp_cs
	99	real, dimension(:), intent(in) :: emis_cs, iwco_cs, rad_cs
	100	real, dimension(:), intent(in) :: pres_cs, dz_cs, rhodz_cs
	101
	102	real, intent(out) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, &
	103	& cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
	104	& pcld_hc_cs, tcld_hc_cs, em_hc_cs, iwp_hc_cs, &
	105	& pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
	106	& pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
	107	& pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
	108	& em_hist_cs, iwp_hist_cs, &
	109	& deltaz_hc_cs, deltaz_hist_cs, rad_hist_cs
	110
	111	real, dimension(len_cs) :: rneb_ord
	112	real :: rneb_min
	113
	114	real, dimension(:), allocatable :: s, s_hc, s_hist, rneb_max
	115	real, dimension(:), allocatable :: sCb, sThCi, sAnv
	116	real, dimension(:), allocatable :: iwp_ss, pcld_ss, tcld_ss,&
	117	& emis_ss
	118	real, dimension(:), allocatable :: deltaz_ss, rad_ss
	119
	120
	121	write(,) '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) .gt. 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) .gt. 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 .gt. maxval(rneb_cs) .and. &
	301	& abs(cc_tot_cs-maxval(rneb_cs)) .gt. 1.e-4 ) then
	302	write(,) 'cc_tot_cs > max rneb_cs'
	303	write(,) cc_tot_cs, maxval(rneb_cs)
	304	STOP
	305	endif
	306
	307	if (iwp_hc_cs .lt. 0.) then
	308	write(,) 'cloud_structure:'
	309	write(,) 'iwp_hc_cs < 0'
	310	STOP
	311	endif
	312
	313	end subroutine cloud_structure
	314
	315
	316	subroutine normal_undef(num, den)
	317
	318	real, intent(in) :: den
	319	real, intent(inout) :: num
	320
	321	if (den .ne. 0) then
	322	num = num/den
	323	else
	324	num = undef
	325	endif
	326
	327	end subroutine normal_undef
	328
	329
	330	subroutine normal2_undef(res,num,den)
	331
	332	real, intent(in) :: den
	333	real, intent(in) :: num
	334	real, intent(out) :: res
	335
	336	if (den .ne. 0.) then
	337	res = num/den
	338	else
	339	res = undef
	340	endif
	341
	342	end subroutine normal2_undef
	343
	344
	345	subroutine sous_section(len_cs, rneb_cs, temp_cs, &
	346	& emis_cs, iwco_cs, &
	347	& pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
	348	& rnebmax, stot, shc, shist, &
	349	& sCb, sThCi, sAnv, &
	350	& emis, pcld, tcld, iwp, deltaz, rad)
	351
	352	integer, intent(in) :: len_cs
	353	real, dimension(len_cs), intent(in) :: rneb_cs, temp_cs
	354	real, dimension(len_cs), intent(in) :: emis_cs, iwco_cs, &
	355	& rneb_ord
	356	real, dimension(len_cs), intent(in) :: pres_cs, dz_cs, rad_cs
	357	real, dimension(len_cs), intent(in) :: rhodz_cs
	358	real, dimension(len_cs) :: tau_cs, w
	359	real, intent(in) :: rnebmax
	360	real, intent(inout) :: stot, shc, shist
	361	real, intent(inout) :: sCb, sThCi, sAnv
	362	real, intent(out) :: emis, pcld, tcld, iwp, deltaz, rad
	363
	364	integer :: i, ideb, ibeg, iend, nuage, visible
	365	real :: som, som_tau, som_iwc, som_dz, som_rad, tau
	366
	367
	368	! Ponderation: 1 pour les nuages, 0 pour les trous
	369
	370	do i = 1, len_cs
	371	if (rneb_cs(i) .ge. rnebmax) then
	372	w(i) = 1.
	373	else
	374	w(i) = 0.
	375	endif
	376	enddo
	377
	378	! Calcul des epaisseurs optiques a partir des emissivites
	379
	380	som = 0.
	381	do i = 1, len_cs
	382	if (emis_cs(i) .eq. 1.) then
	383	tau_cs(i) = 10.
	384	else
	385	tau_cs(i) = -log(1.-emis_cs(i))
	386	endif
	387	som = som+tau_cs(i)
	388	enddo
	389
	390
	391	ideb = 1
	392	nuage = 0
	393	visible = 0
	394
	395
	396	! Boucle sur les nuages
	397	do while (ideb .ne. 0 .and. ideb .le. 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 .gt. 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 .gt. 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 .gt. 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 .lt. em_min .or. emis .gt. em_max &
	453	& .or. tcld .gt. 230.) then
	454	shist = 0.
	455	endif
	456
	457	if (emis .lt. 0.98) then
	458	sCb = 0.
	459	endif
	460
	461	if (emis .gt. 0.5 .or. emis .lt. 0.1) then
	462	sThCi = 0.
	463	endif
	464
	465	if (emis .le. 0.5 .or. emis .ge. 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 .lt. 0.) then
	493	write(,) 'ideb iwp =', ideb, iwp
	494	STOP
	495	endif
	496
	497	if (deltaz .lt. 0.) then
	498	write(,) 'ideb deltaz =', ideb, deltaz
	499	STOP
	500	endif
	501
	502	if (emis .lt. 0.048 .and. emis .ne. 0.) then
	503	write(,) 'ideb emis =', ideb, emis
	504	STOP
	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 .eq. 0) then
	528
	529	if (som_tau .lt. 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
	548	end subroutine masque
	549
	550
	551	subroutine caract (ibeg, iend, temp_cs, tau_cs, iwco_cs, &
	552	& pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
	553	& som_tau, som_iwc, som_dz, som_rad)
	554
	555	integer, intent(in) :: ibeg, iend
	556	real, dimension(:), intent(in) :: tau_cs, iwco_cs, temp_cs
	557	real, dimension(:), intent(in) :: pres_cs, dz_cs, rad_cs
	558	real, dimension(:), intent(in) :: rhodz_cs
	559	real, intent(out) :: som_tau, som_iwc, som_dz, som_rad
	560	real , intent(out) :: pcld, tcld
	561
	562	integer :: i, ibase, imid
	563
	564	! Somme des epaisseurs optiques et des contenus en glace sur le nuage
	565
	566	som_tau = 0.
	567	som_iwc = 0.
	568	som_dz = 0.
	569	som_rad = 0.
	570	ibase = -100
	571
	572	do i = ibeg, iend
	573
	574	som_tau = som_tau + tau_cs(i)
	575
	576	som_dz = som_dz + dz_cs(i)
	577	som_iwc = som_iwc + iwco_cs(i)1000rhodz_cs(i) ! en g/m2
	578	som_rad = som_rad + rad_cs(i)*dz_cs(i)
	579
	580	if (som_tau .gt. 3. .and. ibase .eq. -100) then
	581	ibase = i-1
	582	endif
	583
	584	enddo
	585
	586	if (som_dz .ne. 0.) then
	587	som_rad = som_rad/som_dz
	588	else
	589	write(,) 'som_dez = 0 STOP'
	590	write(,) 'ibeg, iend =', ibeg, iend
	591	do i = ibeg, iend
	592	write(,) dz_cs(i), rhodz_cs(i)
	593	enddo
	594	STOP
	595	endif
	596
	597	! Determination de Pcld et Tcld
	598
	599	if (ibase .lt. ibeg) then
	600	ibase = ibeg
	601	endif
	602
	603	imid = (ibeg+ibase)/2
	604
	605	pcld = pres_cs(imid)/100. ! pcld en hPa
	606	tcld = temp_cs(imid)
	607
	608
	609	end subroutine caract
	610
	611	subroutine topbot(ideb,w,ibeg,iend)
	612
	613	integer, intent(in) :: ideb
	614	real, dimension(:), intent(in) :: w
	615	integer, intent(out) :: ibeg, iend
	616
	617	integer :: i, itest
	618
	619	itest = 0
	620	ibeg = 0
	621	iend = 0
	622
	623	do i = ideb, size(w)
	624
	625	if (w(i) .eq. 1. .and. itest .eq. 0) then
	626	ibeg = i
	627	itest = 1
	628	endif
	629
	630	enddo
	631
	632
	633	i = ibeg
	634	do while (w(i) .eq. 1. .and. i .le. size(w))
	635	i = i+1
	636	enddo
	637	iend = i-1
	638
	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
	654	do i = 1, size(rneb_cs)
	655	rneb(i) = rneb_cs(i)
	656	enddo
	657
	658	do j = 1, size(rneb_cs)
	659
	660	rneb_max = 100.
	661
	662	do i = 1, size(rneb_cs)
	663	if (rneb(i) .lt. rneb_max) then
	664	rneb_max = rneb(i)
	665	ind_min = i
	666	endif
	667	enddo
	668
	669	rneb_ord(j) = rneb_max
	670	rneb(ind_min) = 100.
	671
	672	enddo
	673
	674	end subroutine ordonne
	675
	676
	677	subroutine sim_mesh(rneb_1D, temp_1D, emis_1D, &
	678	& iwcon_1D, rad_1D, &
	679	& pres, dz, &
	680	& rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh,&
	681	& tcld_hc_mesh, &
	682	& em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, &
	683	& cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
	684	& pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
	685	& pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
	686	& pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
	687	& em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
	688
	689	USE dimphy
	690
	691	real, dimension(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, &
	692	& iwcon_1D, rad_1D
	693	real, dimension(klev), intent(in) :: pres, dz, rhodz_1D
	694	real, intent(out) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
	695	real, intent(out) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
	696	real, intent(out) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, &
	697	& iwp_hc_mesh
	698
	699	real, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
	700	real, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
	701	& em_ThCi_mesh
	702	real, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
	703
	704	real, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
	705	real, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh
	706
	707	real, dimension(:), allocatable :: rneb_cs, temp_cs, emis_cs, &
	708	& iwco_cs
	709	real, dimension(:), allocatable :: pres_cs, dz_cs, rad_cs, &
	710	& rhodz_cs
	711
	712	integer :: i,j,l
	713	integer :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics
	714
	715	real :: som_emi_hc,som_pcl_hc,som_tcl_hc,som_iwp_hc,som_hc,&
	716	& som_hist
	717	real :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
	718	& som_deltaz_hist
	719	real :: som_rad_hist
	720	real :: som_Cb, som_ThCi, som_Anv
	721	real :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb
	722	real :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv
	723	real :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi
	724	real :: tsom_tot, tsom_hc, tsom_hist
	725	real :: prod, prod_hh
	726
	727	real :: cc_tot_cs, cc_hc_cs, cc_hist_cs
	728	real :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs
	729	real :: pcld_hc_cs, tcld_hc_cs
	730	real :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs
	731	real :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs
	732	real :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs
	733	real :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs
	734	real :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs
	735
	736	real, dimension(klev) :: test_tot, test_hc, test_hist
	737	real, dimension(klev) :: test_pcld, test_tcld, test_em, test_iwp
	738
	739
	740	do j = 1, klev
	741	write(,) 'simu_airs, j, rneb_1D =', rneb_1D(j)
	742	enddo
	743
	744	! Definition des structures nuageuses, de la plus haute a la plus basse
	745
	746	num_cs = 0
	747	ltop = klev-1
	748
	749	prod = 1.
	750
	751	som_emi_hc = 0.
	752	som_emi_hist = 0.
	753	som_pcl_hc = 0.
	754	som_tcl_hc = 0.
	755	som_iwp_hc = 0.
	756	som_iwp_hist = 0.
	757	som_deltaz_hc = 0.
	758	som_deltaz_hist = 0.
	759	som_rad_hist = 0.
	760	som_hc = 0.
	761	som_hist = 0.
	762
	763	som_Cb = 0.
	764	som_ThCi = 0.
	765	som_Anv = 0.
	766
	767	som_emi_Cb = 0.
	768	som_pcld_Cb = 0.
	769	som_tcld_Cb = 0.
	770
	771	som_emi_ThCi = 0.
	772	som_pcld_ThCi = 0.
	773	som_tcld_ThCi = 0.
	774
	775	som_emi_Anv = 0.
	776	som_pcld_Anv = 0.
	777	som_tcld_Anv = 0.
	778
	779	tsom_tot = 0.
	780	tsom_hc = 0.
	781	tsom_hist = 0.
	782
	783	do while (ltop .ge. 1) ! Boucle sur toute la colonne
	784
	785	itop = 0
	786
	787	do l = ltop,1,-1
	788
	789	if (itop .eq. 0 .and. rneb_1D(l) .gt. 0.001 ) then
	790	itop = l
	791	endif
	792
	793	enddo
	794
	795	ibot = itop
	796
	797	do while (rneb_1D(ibot) .gt. 0.001 .and. ibot .ge. 1)
	798	ibot = ibot -1
	799	enddo
	800
	801
	802	ibot = ibot+1
	803
	804	if (itop .gt. 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
	849	deallocate (rneb_cs)
	850	deallocate (temp_cs)
	851	deallocate (emis_cs)
	852	deallocate (iwco_cs)
	853	deallocate (pres_cs)
	854	deallocate (dz_cs)
	855	deallocate (rad_cs)
	856	deallocate (rhodz_cs)
	857
	858
	859	! Pour la couverture nuageuse sur la maille
	860
	861	prod_hh = prod
	862
	863	prod = prod*(1.-cc_tot_cs)
	864
	865	! Pour les autres variables definies sur la maille
	866
	867	som_emi_hc = som_emi_hc + em_hc_cscc_hc_csprod_hh
	868	som_iwp_hc = som_iwp_hc + iwp_hc_cscc_hc_csprod_hh
	869	som_deltaz_hc = som_deltaz_hc + deltaz_hc_cscc_hc_csprod_hh
	870
	871	som_emi_Cb = som_emi_Cb + em_Cb_cscc_Cb_csprod_hh
	872	som_tcld_Cb = som_tcld_Cb + tcld_Cb_cscc_Cb_csprod_hh
	873	som_pcld_Cb = som_pcld_Cb + pcld_Cb_cscc_Cb_csprod_hh
	874
	875	som_emi_ThCi = som_emi_ThCi + em_ThCi_cscc_ThCi_csprod_hh
	876	som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cscc_ThCi_csprod_hh
	877	som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cscc_ThCi_csprod_hh
	878
	879	som_emi_Anv = som_emi_Anv + em_Anv_cscc_Anv_csprod_hh
	880	som_tcld_Anv = som_tcld_Anv + tcld_Anv_cscc_Anv_csprod_hh
	881	som_pcld_Anv = som_pcld_Anv + pcld_Anv_cscc_Anv_csprod_hh
	882
	883	som_emi_hist = som_emi_hist + em_hist_cscc_hist_csprod_hh
	884	som_iwp_hist = som_iwp_hist + iwp_hist_cscc_hist_csprod_hh
	885	som_deltaz_hist = som_deltaz_hist + &
	886	& deltaz_hist_cscc_hist_csprod_hh
	887	som_rad_hist = som_rad_hist + rad_hist_cscc_hist_csprod_hh
	888
	889	som_pcl_hc = som_pcl_hc + pcld_hc_cscc_hc_csprod_hh
	890	som_tcl_hc = som_tcl_hc + tcld_hc_cscc_hc_csprod_hh
	891
	892	som_hc = som_hc + cc_hc_cs*prod_hh
	893	som_hist = som_hist + cc_hist_cs*prod_hh
	894
	895	som_Cb = som_Cb + cc_Cb_cs*prod_hh
	896	som_ThCi = som_ThCi + cc_ThCi_cs*prod_hh
	897	som_Anv = som_Anv + cc_Anv_cs*prod_hh
	898
	899
	900	! Pour test
	901
	902	call test_bornes('cc_tot_cs ',cc_tot_cs,1.,0.)
	903	call test_bornes('cc_hc_cs ',cc_hc_cs,1.,0.)
	904	call test_bornes('cc_hist_cs ',cc_hist_cs,1.,0.)
	905	call test_bornes('pcld_hc_cs ',pcld_hc_cs,1200.,0.)
	906	call test_bornes('tcld_hc_cs ',tcld_hc_cs,1000.,100.)
	907	call test_bornes('em_hc_cs ',em_hc_cs,1000.,0.048)
	908
	909	test_tot(num_cs) = cc_tot_cs
	910	test_hc(num_cs) = cc_hc_cs
	911	test_hist(num_cs) = cc_hist_cs
	912	test_pcld(num_cs) = pcld_hc_cs
	913	test_tcld(num_cs) = tcld_hc_cs
	914	test_em(num_cs) = em_hc_cs
	915	test_iwp(num_cs) = iwp_hc_cs
	916
	917	tsom_tot = tsom_tot + cc_tot_cs
	918	tsom_hc = tsom_hc + cc_hc_cs
	919	tsom_hist = tsom_hist + cc_hist_cs
	920
	921
	922	endif ! itop > 0
	923
	924	ltop = ibot -1
	925
	926	enddo ! fin de la boucle sur la colonne
	927
	928	N_CS = num_cs
	929
	930
	931	! Determination des variables de sortie
	932
	933	if (N_CS .gt. 0) then ! if N_CS>0
	934
	935	cc_tot_mesh = 1. - prod
	936
	937	cc_hc_mesh = som_hc
	938	cc_hist_mesh = som_hist
	939
	940	cc_Cb_mesh = som_Cb
	941	cc_ThCi_mesh = som_ThCi
	942	cc_Anv_mesh = som_Anv
	943
	944	call normal2_undef(pcld_hc_mesh,som_pcl_hc, &
	945	& cc_hc_mesh)
	946	call normal2_undef(tcld_hc_mesh,som_tcl_hc, &
	947	& cc_hc_mesh)
	948	call normal2_undef(em_hc_mesh,som_emi_hc, &
	949	& cc_hc_mesh)
	950	call normal2_undef(iwp_hc_mesh,som_iwp_hc, &
	951	& cc_hc_mesh)
	952	call normal2_undef(deltaz_hc_mesh,som_deltaz_hc, &
	953	& cc_hc_mesh)
	954
	955	call normal2_undef(em_Cb_mesh,som_emi_Cb, &
	956	& cc_Cb_mesh)
	957	call normal2_undef(tcld_Cb_mesh,som_tcld_Cb, &
	958	& cc_Cb_mesh)
	959	call normal2_undef(pcld_Cb_mesh,som_pcld_Cb, &
	960	& cc_Cb_mesh)
	961
	962	call normal2_undef(em_ThCi_mesh,som_emi_ThCi, &
	963	& cc_ThCi_mesh)
	964	call normal2_undef(tcld_ThCi_mesh,som_tcld_ThCi, &
	965	& cc_ThCi_mesh)
	966	call normal2_undef(pcld_ThCi_mesh,som_pcld_ThCi, &
	967	& cc_ThCi_mesh)
	968
	969	call normal2_undef(em_Anv_mesh,som_emi_Anv, &
	970	& cc_Anv_mesh)
	971	call normal2_undef(tcld_Anv_mesh,som_tcld_Anv, &
	972	& cc_Anv_mesh)
	973	call normal2_undef(pcld_Anv_mesh,som_pcld_Anv, &
	974	& cc_Anv_mesh)
	975
	976
	977	call normal2_undef(em_hist_mesh,som_emi_hist, &
	978	& cc_hist_mesh)
	979	call normal2_undef(iwp_hist_mesh,som_iwp_hist, &
	980	& cc_hist_mesh)
	981	call normal2_undef(deltaz_hist_mesh,som_deltaz_hist, &
	982	& cc_hist_mesh)
	983	call normal2_undef(rad_hist_mesh,som_rad_hist, &
	984	& cc_hist_mesh)
	985
	986
	987	! Tests
	988
	989	! Tests
	990
	991	if (cc_tot_mesh .gt. tsom_tot .and. &
	992	& abs(cc_tot_mesh-tsom_tot) .gt. 1.e-4) then
	993	write(,) 'cc_tot_mesh > tsom_tot'
	994	write(,) cc_tot_mesh, tsom_tot
	995	STOP
	996	endif
	997
	998	if (cc_tot_mesh .lt. maxval(test_tot(1:N_CS)) .and. &
	999	& abs(cc_tot_mesh-maxval(test_tot(1:N_CS))) .gt. 1.e-4) then
	1000	write(,) 'cc_tot_mesh < max'
	1001	write(,) cc_tot_mesh, maxval(test_tot(1:N_CS))
	1002	STOP
	1003	endif
	1004
	1005	if (cc_hc_mesh .gt. tsom_hc .and. &
	1006	& abs(cc_hc_mesh-tsom_hc) .gt. 1.e-4) then
	1007	write(,) 'cc_hc_mesh > tsom_hc'
	1008	write(,) cc_hc_mesh, tsom_hc
	1009	STOP
	1010	endif
	1011
	1012	if (cc_hc_mesh .lt. maxval(test_hc(1:N_CS)) .and. &
	1013	& abs(cc_hc_mesh-maxval(test_hc(1:N_CS))) .gt. 1.e-4) then
	1014	write(,) 'cc_hc_mesh < max'
	1015	write(,) cc_hc_mesh, maxval(test_hc(1:N_CS))
	1016	STOP
	1017	endif
	1018
	1019	if (cc_hist_mesh .gt. tsom_hist .and. &
	1020	& abs(cc_hist_mesh-tsom_hist) .gt. 1.e-4) then
	1021	write(,) 'cc_hist_mesh > tsom_hist'
	1022	write(,) cc_hist_mesh, tsom_hist
	1023	STOP
	1024	endif
	1025
	1026	if (cc_hist_mesh .lt. 0.) then
	1027	write(,) 'cc_hist_mesh < 0'
	1028	write(,) cc_hist_mesh
	1029	STOP
	1030	endif
	1031
	1032	if ((pcld_hc_mesh .gt. maxval(test_pcld(1:N_CS)) .or. &
	1033	& pcld_hc_mesh .lt. minval(test_pcld(1:N_CS))) .and. &
	1034	& abs(pcld_hc_mesh-maxval(test_pcld(1:N_CS))) .gt. 1. .and. &
	1035	& maxval(test_pcld(1:N_CS)) .ne. 999. &
	1036	& .and. minval(test_pcld(1:N_CS)) .ne. 999.) then
	1037	write(,) 'pcld_hc_mesh est faux'
	1038	write(,) pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
	1039	& minval(test_pcld(1:N_CS))
	1040	STOP
	1041	endif
	1042
	1043	if ((tcld_hc_mesh .gt. maxval(test_tcld(1:N_CS)) .or. &
	1044	& tcld_hc_mesh .lt. minval(test_tcld(1:N_CS))) .and. &
	1045	& abs(tcld_hc_mesh-maxval(test_tcld(1:N_CS))) .gt. 0.1 .and. &
	1046	& maxval(test_tcld(1:N_CS)) .ne. 999. &
	1047	& .and. minval(test_tcld(1:N_CS)) .ne. 999.) then
	1048	write(,) 'tcld_hc_mesh est faux'
	1049	write(,) tcld_hc_mesh, maxval(test_tcld(1:N_CS)), &
	1050	& minval(test_tcld(1:N_CS))
	1051	endif
	1052
	1053	if ((em_hc_mesh .gt. maxval(test_em(1:N_CS)) .or. &
	1054	& em_hc_mesh .lt. minval(test_em(1:N_CS))) .and. &
	1055	& abs(em_hc_mesh-maxval(test_em(1:N_CS))) .gt. 1.e-4 .and. &
	1056	& minval(test_em(1:N_CS)) .ne. 999. .and. &
	1057	& maxval(test_em(1:N_CS)) .ne. 999. ) then
	1058	write(,) 'em_hc_mesh est faux'
	1059	write(,) em_hc_mesh, maxval(test_em(1:N_CS)), &
	1060	& minval(test_em(1:N_CS))
	1061	STOP
	1062	endif
	1063
	1064	else ! if N_CS>0
	1065
	1066	cc_tot_mesh = 0.
	1067	cc_hc_mesh = 0.
	1068	cc_hist_mesh = 0.
	1069
	1070	cc_Cb_mesh = 0.
	1071	cc_ThCi_mesh = 0.
	1072	cc_Anv_mesh = 0.
	1073
	1074	iwp_hc_mesh = undef
	1075	deltaz_hc_mesh = undef
	1076	em_hc_mesh = undef
	1077	iwp_hist_mesh = undef
	1078	deltaz_hist_mesh = undef
	1079	rad_hist_mesh = undef
	1080	em_hist_mesh = undef
	1081	pcld_hc_mesh = undef
	1082	tcld_hc_mesh = undef
	1083
	1084	pcld_Cb_mesh = undef
	1085	tcld_Cb_mesh = undef
	1086	em_Cb_mesh = undef
	1087
	1088	pcld_ThCi_mesh = undef
	1089	tcld_ThCi_mesh = undef
	1090	em_ThCi_mesh = undef
	1091
	1092	pcld_Anv_mesh = undef
	1093	tcld_Anv_mesh = undef
	1094	em_Anv_mesh = undef
	1095
	1096
	1097	endif ! if N_CS>0
	1098
	1099	end subroutine sim_mesh
	1100
	1101	subroutine test_bornes(sx,x,bsup,binf)
	1102
	1103	real, intent(in) :: x, bsup, binf
	1104	character*14, intent(in) :: sx
	1105
	1106	if (x .gt. bsup .or. x .lt. binf) then
	1107	write(,) sx, 'est faux'
	1108	write(,) sx, x
	1109	STOP
	1110	endif
	1111
	1112	end subroutine test_bornes
	1113
	1114	end module m_simu_airs
	1115
	1116
	1117	subroutine simu_airs &
	1118	& (itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, &
	1119	& geop_airs, pplay_airs, paprs_airs, &
	1120	& map_prop_hc,map_prop_hist,&
	1121	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
	1122	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
	1123	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
	1124	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
	1125	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
	1126	& map_ntot,map_hc,map_hist,&
	1127	& map_Cb,map_ThCi,map_Anv,alt_tropo )
	1128
	1129	USE dimphy
	1130	USE m_simu_airs
	1131
	1132	IMPLICIT NONE
	1133
	1134	include "YOMCST.h"
	1135
	1136	integer,intent(in) :: itap
	1137
	1138	real, dimension(klon,klev), intent(in) :: &
	1139	& rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
	1140	& rad_airs, geop_airs, pplay_airs, paprs_airs
	1141
	1142	real, dimension(klon,klev) :: &
	1143	& rhodz_airs, rho_airs, iwcon_airs
	1144
	1145	real, dimension(klon),intent(out) :: alt_tropo
	1146
	1147	real, dimension(klev) :: rneb_1D, temp_1D, &
	1148	& emis_1D, rad_1D, pres_1D, alt_1D, &
	1149	& rhodz_1D, dz_1D, iwcon_1D
	1150
	1151	integer :: i, j
	1152
	1153	real :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
	1154	real :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
	1155	real :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
	1156	real :: em_hist_mesh, iwp_hist_mesh
	1157	real :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh
	1158	real :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
	1159	real :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh
	1160	real :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
	1161
	1162	real, dimension(klon),intent(out) :: map_prop_hc, map_prop_hist
	1163	real, dimension(klon),intent(out) :: map_emis_hc, map_iwp_hc
	1164	real, dimension(klon),intent(out) :: map_deltaz_hc, map_pcld_hc
	1165	real, dimension(klon),intent(out) :: map_tcld_hc
[2585]	1166	real, dimension(klon),intent(out) :: map_emis_Cb,map_pcld_Cb,map_tcld_Cb
[2580]	1167	real, dimension(klon),intent(out) :: &
	1168	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi
	1169	real, dimension(klon),intent(out) :: &
	1170	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv
	1171	real, dimension(klon),intent(out) :: &
	1172	& map_emis_hist,map_iwp_hist,map_deltaz_hist,&
	1173	& map_rad_hist
[2585]	1174	real, dimension(klon),intent(out) :: map_ntot,map_hc,map_hist
	1175	real, dimension(klon),intent(out) :: map_Cb,map_ThCi,map_Anv
[2580]	1176
	1177
	1178	write(,) 'simu_airs'
	1179	write(,) 'itap, klon, klev', itap, klon, klev
	1180	write(,) 'RG, RD =', RG, RD
	1181
	1182
	1183	! Definition des variables 1D
	1184
	1185	do i = 1, klon
	1186	do j = 1, klev-1
	1187	rhodz_airs(i,j) = &
	1188	& (paprs_airs(i,j)-paprs_airs(i,j+1))/RG
	1189	enddo
	1190	rhodz_airs(i,klev) = 0.
	1191	enddo
	1192
	1193	do i = 1, klon
	1194	do j = 1,klev
	1195	rho_airs(i,j) = &
	1196	& pplay_airs(i,j)/(temp_airs(i,j)*RD)
	1197
	1198	if (rneb_airs(i,j) .gt. 0.001) then
	1199	iwcon_airs(i,j) = iwcon0_airs(i,j)/rneb_airs(i,j)
	1200	else
	1201	iwcon_airs(i,j) = 0.
	1202	endif
	1203
	1204	enddo
	1205	enddo
	1206
	1207	!=============================================================================
	1208
	1209	do i = 1, klon ! boucle sur les points de grille
	1210
	1211	!=============================================================================
	1212
	1213	do j = 1,klev
	1214
	1215	rneb_1D(j) = rneb_airs(i,j)
	1216	temp_1D(j) = temp_airs(i,j)
	1217	emis_1D(j) = cldemi_airs(i,j)
	1218	iwcon_1D(j) = iwcon_airs(i,j)
	1219	rad_1D(j) = rad_airs(i,j)
	1220	pres_1D(j) = pplay_airs(i,j)
	1221	alt_1D(j) = geop_airs(i,j)/RG
	1222	rhodz_1D(j) = rhodz_airs(i,j)
	1223	dz_1D(j) = rhodz_airs(i,j)/rho_airs(i,j)
	1224
	1225	enddo
	1226
	1227	alt_tropo(i) = &
	1228	& search_tropopause(pres_1D/100.,temp_1D,alt_1D,klev)
	1229
	1230
	1231	! Appel du ss-programme sim_mesh
	1232
	1233	! if (itap .eq. 1 ) then
	1234
	1235	call sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, &
	1236	& pres_1D, dz_1D, rhodz_1D, &
	1237	& cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, &
	1238	& pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, &
	1239	& deltaz_hc_mesh,&
	1240	& cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
	1241	& pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
	1242	& pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
	1243	& pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
	1244	& em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
	1245
	1246	write(,) '===================================='
	1247	write(,) 'itap, i:', itap, i
	1248	write(,) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, &
	1249	& iwp_hc, em_hist, iwp_hist ='
	1250	write(,) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
	1251	write(,) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
	1252	write(,) em_hist_mesh, iwp_hist_mesh
	1253
	1254	! endif
	1255
	1256	! Definition des variables a ecrire dans le fichier de sortie
	1257
	1258	call normal2_undef(map_prop_hc(i),cc_hc_mesh, &
	1259	& cc_tot_mesh)
	1260	call normal2_undef(map_prop_hist(i),cc_hist_mesh, &
	1261	& cc_tot_mesh)
	1262
	1263	map_emis_hc(i) = em_hc_mesh
	1264	map_iwp_hc(i) = iwp_hc_mesh
	1265	map_deltaz_hc(i) = deltaz_hc_mesh
	1266	map_pcld_hc(i) = pcld_hc_mesh
	1267	map_tcld_hc(i) = tcld_hc_mesh
	1268
	1269	map_emis_Cb(i) = em_Cb_mesh
	1270	map_pcld_Cb(i) = pcld_Cb_mesh
	1271	map_tcld_Cb(i) = tcld_Cb_mesh
	1272
	1273	map_emis_ThCi(i) = em_ThCi_mesh
	1274	map_pcld_ThCi(i) = pcld_ThCi_mesh
	1275	map_tcld_ThCi(i) = tcld_ThCi_mesh
	1276
	1277	map_emis_Anv(i) = em_Anv_mesh
	1278	map_pcld_Anv(i) = pcld_Anv_mesh
	1279	map_tcld_Anv(i) = tcld_Anv_mesh
	1280
	1281	map_emis_hist(i) = em_hist_mesh
	1282	map_iwp_hist(i) = iwp_hist_mesh
	1283	map_deltaz_hist(i) = deltaz_hist_mesh
	1284	map_rad_hist(i) = rad_hist_mesh
	1285
	1286	map_ntot(i) = cc_tot_mesh
	1287	map_hc(i) = cc_hc_mesh
	1288	map_hist(i) = cc_hist_mesh
	1289
	1290	map_Cb(i) = cc_Cb_mesh
	1291	map_ThCi(i) = cc_ThCi_mesh
	1292	map_Anv(i) = cc_Anv_mesh
	1293
	1294
	1295	enddo ! fin boucle sur les points de grille
	1296
	1297
	1298
	1299	end subroutine simu_airs
	1300

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

Download in other formats:

Original Format