Context Navigation

simu_airs.F90 @ 3709

Last change on this file since 3709 was 3531, checked in by Laurent Fairhead, 5 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

Rev	Line
[2580]	1
	2	module m_simu_airs
	3
[3531]	4	USE print_control_mod, ONLY: prt_level,lunout
	5
[2580]	6	implicit none
	7
[3531]	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.
[2580]	13
	14	contains
	15
[3531]	16	REAL function search_tropopause(P,T,alt,N) result(P_tropo)
[2580]	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
[3531]	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, &
[2580]	25	& delta_alt_limit,tmp,delta_alt
[3531]	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
[2580]	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
[3531]	97	INTEGER :: i, n, nss
[2580]	98
[3531]	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
[2580]	103
[3531]	104	REAL, intent(out) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, &
[2580]	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
[3531]	113	REAL, DIMENSION(len_cs) :: rneb_ord
	114	REAL :: rneb_min
[2580]	115
[3531]	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,&
[2580]	119	& emis_ss
[3531]	120	REAL, DIMENSION(:), allocatable :: deltaz_ss, rad_ss
[2580]	121
[3531]	122	CHARACTER (len = 50) :: modname = 'simu_airs.cloud_structure'
	123	CHARACTER (len = 160) :: abort_message
	124
[2580]	125
[3531]	126	write(lunout,*) 'dans cloud_structure'
[2580]	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
[3531]	307	WRITE(abort_message,*) 'cc_tot_cs > max rneb_cs', cc_tot_cs, maxval(rneb_cs)
	308	CALL abort_physic(modname,abort_message,1)
[2580]	309	endif
	310
	311	if (iwp_hc_cs .lt. 0.) then
[3531]	312	abort_message= 'cloud_structure: iwp_hc_cs < 0'
	313	CALL abort_physic(modname,abort_message,1)
[2580]	314	endif
	315
	316	end subroutine cloud_structure
	317
	318
	319	subroutine normal_undef(num, den)
	320
[3531]	321	REAL, intent(in) :: den
	322	REAL, intent(inout) :: num
[2580]	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
[3531]	335	REAL, intent(in) :: den
	336	REAL, intent(in) :: num
	337	REAL, intent(out) :: res
[2580]	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
[3531]	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, &
[2580]	358	& rneb_ord
[3531]	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
[2580]	366
[3531]	367	INTEGER :: i, ideb, ibeg, iend, nuage, visible
	368	REAL :: som, som_tau, som_iwc, som_dz, som_rad, tau
[2580]	369
[3531]	370	CHARACTER (len = 50) :: modname = 'simu_airs.sous_section'
	371	CHARACTER (len = 160) :: abort_message
[2580]	372
[3531]	373
[2580]	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
[3531]	499	WRITE(abort_message,*) 'ideb iwp =', ideb, iwp
	500	CALL abort_physic(modname,abort_message,1)
[2580]	501	endif
	502
	503	if (deltaz .lt. 0.) then
[3531]	504	WRITE(abort_message,*)'ideb deltaz =', ideb, deltaz
	505	CALL abort_physic(modname,abort_message,1)
[2580]	506	endif
	507
	508	if (emis .lt. 0.048 .and. emis .ne. 0.) then
[3531]	509	WRITE(abort_message,*) 'ideb emis =', ideb, emis
	510	CALL abort_physic(modname,abort_message,1)
[2580]	511	endif
	512
	513	end subroutine sous_section
	514
	515
	516	subroutine masque (ibeg, iend, som_tau, &
	517	& visible, w)
	518
[3531]	519	INTEGER, intent(in) :: ibeg, iend
	520	REAL, intent(in) :: som_tau
[2580]	521
[3531]	522	INTEGER, intent(inout) :: visible
	523	REAL, DIMENSION(:), intent(inout) :: w
[2580]	524
[3531]	525	INTEGER :: i
[2580]	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
[3531]	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
[2580]	567
[3531]	568	INTEGER :: i, ibase, imid
[2580]	569
[3531]	570	CHARACTER (len = 50) :: modname = 'simu_airs.caract'
	571	CHARACTER (len = 160) :: abort_message
	572
[2580]	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
[3531]	596	som_rad = som_rad/som_dz
[2580]	597	else
[3531]	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)
[2580]	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
[3531]	623	INTEGER, intent(in) :: ideb
	624	REAL, DIMENSION(:), intent(in) :: w
	625	INTEGER, intent(out) :: ibeg, iend
[2580]	626
[3531]	627	INTEGER :: i, itest
[2580]	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
[3531]	654	INTEGER, intent(in) :: len_cs
	655	REAL, DIMENSION(:), intent(in) :: rneb_cs
	656	REAL, DIMENSION(:), intent(out) :: rneb_ord
[2580]	657
[3531]	658	INTEGER :: i, j, ind_min
[2580]	659
[3531]	660	REAL, DIMENSION(len_cs) :: rneb
	661	REAL :: rneb_max
[2580]	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
[3531]	701	REAL, DIMENSION(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, &
[2580]	702	& iwcon_1D, rad_1D
[3531]	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, &
[2580]	707	& iwp_hc_mesh
	708
[3531]	709	REAL, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
	710	REAL, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
[2580]	711	& em_ThCi_mesh
[3531]	712	REAL, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
[2580]	713
[3531]	714	REAL, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
	715	REAL, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh
[2580]	716
[3531]	717	REAL, DIMENSION(:), allocatable :: rneb_cs, temp_cs, emis_cs, &
[2580]	718	& iwco_cs
[3531]	719	REAL, DIMENSION(:), allocatable :: pres_cs, dz_cs, rad_cs, &
[2580]	720	& rhodz_cs
	721
[3531]	722	INTEGER :: i,j,l
	723	INTEGER :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics
[2580]	724
[3531]	725	REAL :: som_emi_hc,som_pcl_hc,som_tcl_hc,som_iwp_hc,som_hc,&
[2580]	726	& som_hist
[3531]	727	REAL :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
[2580]	728	& som_deltaz_hist
[3531]	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
[2580]	736
[3531]	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
[2580]	745
[3531]	746	REAL, DIMENSION(klev) :: test_tot, test_hc, test_hist
	747	REAL, DIMENSION(klev) :: test_pcld, test_tcld, test_em, test_iwp
[2580]	748
[3531]	749	CHARACTER (len = 50) :: modname = 'simu_airs.sim_mesh'
	750	CHARACTER (len = 160) :: abort_message
	751
[2580]	752
	753	do j = 1, klev
[3531]	754	WRITE(lunout,*) 'simu_airs, j, rneb_1D =', rneb_1D(j)
[2580]	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
[3531]	1006	WRITE(abort_message,*)'cc_tot_mesh > tsom_tot', cc_tot_mesh, tsom_tot
	1007	CALL abort_physic(modname,abort_message,1)
[2580]	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
[3531]	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)
[2580]	1014	endif
	1015
	1016	if (cc_hc_mesh .gt. tsom_hc .and. &
	1017	& abs(cc_hc_mesh-tsom_hc) .gt. 1.e-4) then
[3531]	1018	WRITE(abort_message,*) 'cc_hc_mesh > tsom_hc', cc_hc_mesh, tsom_hc
	1019	CALL abort_physic(modname,abort_message,1)
[2580]	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
[3531]	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)
[2580]	1026	endif
	1027
	1028	if (cc_hist_mesh .gt. tsom_hist .and. &
	1029	& abs(cc_hist_mesh-tsom_hist) .gt. 1.e-4) then
[3531]	1030	WRITE(abort_message,*) 'cc_hist_mesh > tsom_hist', cc_hist_mesh, tsom_hist
	1031	CALL abort_physic(modname,abort_message,1)
[2580]	1032	endif
	1033
	1034	if (cc_hist_mesh .lt. 0.) then
[3531]	1035	WRITE(abort_message,*) 'cc_hist_mesh < 0', cc_hist_mesh
	1036	CALL abort_physic(modname,abort_message,1)
[2580]	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
[3531]	1044	WRITE(abort_message,*) 'pcld_hc_mesh est faux', pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
[2580]	1045	& minval(test_pcld(1:N_CS))
[3531]	1046	CALL abort_physic(modname,abort_message,1)
[2580]	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
[3531]	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)
[2580]	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
[3531]	1064	WRITE(abort_message,*) 'em_hc_mesh est faux', em_hc_mesh, maxval(test_em(1:N_CS)), &
[2580]	1065	& minval(test_em(1:N_CS))
[3531]	1066	CALL abort_physic(modname,abort_message,1)
[2580]	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
[3531]	1108	REAL, intent(in) :: x, bsup, binf
[2580]	1109	character*14, intent(in) :: sx
[3531]	1110	CHARACTER (len = 50) :: modname = 'simu_airs.test_bornes'
	1111	CHARACTER (len = 160) :: abort_message
[2580]	1112
	1113	if (x .gt. bsup .or. x .lt. binf) then
[3531]	1114	WRITE(abort_message,*) sx, 'est faux', sx, x
	1115	CALL abort_physic(modname,abort_message,1)
[2580]	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
[3531]	1142	INTEGER,intent(in) :: itap
[2580]	1143
[3531]	1144	REAL, DIMENSION(klon,klev), intent(in) :: &
[2580]	1145	& rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
	1146	& rad_airs, geop_airs, pplay_airs, paprs_airs
	1147
[3531]	1148	REAL, DIMENSION(klon,klev) :: &
[2580]	1149	& rhodz_airs, rho_airs, iwcon_airs
	1150
[3531]	1151	REAL, DIMENSION(klon),intent(out) :: alt_tropo
[2580]	1152
[3531]	1153	REAL, DIMENSION(klev) :: rneb_1D, temp_1D, &
[2580]	1154	& emis_1D, rad_1D, pres_1D, alt_1D, &
	1155	& rhodz_1D, dz_1D, iwcon_1D
	1156
[3531]	1157	INTEGER :: i, j
[2580]	1158
[3531]	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
[2580]	1167
[3531]	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) :: &
[2580]	1174	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi
[3531]	1175	REAL, DIMENSION(klon),intent(out) :: &
[2580]	1176	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv
[3531]	1177	REAL, DIMENSION(klon),intent(out) :: &
[2580]	1178	& map_emis_hist,map_iwp_hist,map_deltaz_hist,&
	1179	& map_rad_hist
[3531]	1180	REAL, DIMENSION(klon),intent(out) :: map_ntot,map_hc,map_hist
	1181	REAL, DIMENSION(klon),intent(out) :: map_Cb,map_ThCi,map_Anv
[2580]	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