Context Navigation

inicondens.F90 @ 1956

Last change on this file since 1956 was 1908, checked in by jvatant, 7 years ago

Making chemistry handling more flexible - Major and Final Step (hopefully) !
After preliminary commits r1871-86-87-91-94-98, here comes the major update of the interface
with photochemical module + fix how tendencies for chem and mufi tracers are managed in physiq_mod !
+ Major modifs in :
++ calchim.F90 to comply with flexible resolution, parallelism, upper pressure grid ...
++ physiq_mod.F90 where there was a bug on the update of the tracers and their tendencies for calchim

and calmufi ( we actually were sending non-updated fields to these processes )
We also now put the same tendency on all longitudes within a lat band and not
relative tendency if 2D chemistry ( and we set to zero if ever negs are created )

+ Also modifs to have chemistry in 1D in rcm1d ( and moved gr_kim_vervack in phytitan to be accessible for 1d )
+ In chemistry added a check.c to verify coherence of sizes between C and Fortran
--JVO

File size: 6.6 KB

Line
1	SUBROUTINE inicondens(press,temp,yc)
2
3	!==================================================================================
4	! Purpose
5	! -------
6	! Initialisation des profils de saturation des traceurs chimiques
7	!
8	! Authors
9	! -------
10	! S. Lebonnois
11	! J. Vatant d'Ollone (2017) : adaptation au nouveau titan et passage en F90
12	!==================================================================================
13
14
15	!-----------------------------------------------------------------------
16	! Declarations:
17	! -------------
18
19	use comchem_h, only: nkim, cnames
20	use dimphy
21	use mod_grid_phy_lmdz, only: nbp_lev
22	IMPLICIT NONE
23
24	! Arguments :
25	! -----------
26	REAL, INTENT(IN) :: press(nbp_lev),temp(nbp_lev) ! Pressure (mbar)
27	REAL, INTENT(OUT) :: yc(nbp_lev,nkim) ! Saturation profiles (mol/mol)
28
29	! Local variables :
30	! -----------------
31	INTEGER :: l,ic
32	REAL :: sy,x
33
34	do ic=1,nkim
35	print*, 'traceur CH(', ic, ')=', trim(cnames(ic)),'------------'
36	do l=1,nbp_lev
37
38	! Par defaut, yc est a 1 c'est a dire qu'on ne condense pas
39	yc(l,ic)=1.
40
41	if(trim(cnames(ic)).eq."CH4") then
42	if (temp(l).lt.90.65) then
43	yc(l,ic)= 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / TEMP(l) - &
44	115352.19e0 ) / TEMP(l) + 4055.6016e0 ) / TEMP(l) &
45	+ 453.92414e0 ) / TEMP(l) ) / PRESS(l) * 1013.25e0
46	else
47	yc(l,ic)= 10.0**(3.901408e0 - ( ( 154567.02e0 / TEMP(l) - &
48	1598.8512e0 ) / TEMP(l) + 437.54809e0 ) / TEMP(l)) &
49	/ PRESS(l) * 1013.25e0;
50	endif
51	! maintient a 1.4% minimum
52	if (yc(l,ic).lt.0.014) yc(l,ic)=0.014
53	endif
54
55	if(trim(cnames(ic)).eq."C2H2") then
56	yc(l,ic)= 10.0**(6.09748e0-1644.1e0/TEMP(l)+7.42346e0 &
57	* alog10(1.0e3/TEMP(l)) ) / PRESS(l)*1013.25e0/760.0
58	endif
59
60	if(trim(cnames(ic)).eq."C2H4") then
61	if (temp(l).lt.89.0) then
62	yc(l,ic)= 10.0**(1.5477e0 + (1.0e0/TEMP(l) - 0.011e0) &
63	(16537.0e0(1.0e0/TEMP(l) - 0.011e0) - 1038.1e0)) &
64	/ PRESS(l) * 1.01325e0 / 760.0
65	elseif (temp(l).lt.104.0) then
66	yc(l,ic)= 10.0**(8.724e0 - 901.6e0/(TEMP(l) - 2.555e0) ) &
67	/ PRESS(l) * 1013.25e0 / 760.0
68	elseif (temp(l).lt.120.0) then
69	yc(l,ic)= 10.0*(50.79e0 - 1703.0e0/TEMP(l) - 17.141e0 &
70	alog10(TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0
71	elseif (temp(l).lt.155.0) then
72	yc(l,ic)= 10.0**(6.74756e0 - 585.0e0/(TEMP(l) - 18.16e0) ) &
73	/ PRESS(l) * 1013.25e0 / 760.0
74	endif
75	endif
76
77	if(trim(cnames(ic)).eq."C2H6") then
78	if (temp(l).lt.90.) then
79	yc(l,ic)= 10.0**(10.01e0-1085.0e0/(TEMP(l)-0.561e0) ) &
80	/ PRESS(l) * 1013.25e0 / 760.0e0
81	else
82	yc(l,ic)= 10.0*(5.9366e0 - 1086.17e0/TEMP(l) + 3.83464e0 &
83	alog10(1.0e3/TEMP(l)) ) / PRESS(l)*1013.25e0/760.0
84	endif
85	endif
86
87	if((trim(cnames(ic)).eq."CH3CCH").or.(trim(cnames(ic)).eq."CH2CCH2")) then
88	yc(l,ic)= 10.0*(2.8808e0 - 4.5e0(249.9e0 - TEMP(l)) &
89	/(1.15e0TEMP(l) - 37.485e0) ) / PRESS(l) 1013.25e0 / 760.0e0
90	endif
91
92	if(trim(cnames(ic)).eq."C3H6") then
93	yc(l,ic)= 10.0*(7.4463e0 - 1028.5654e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
94	endif
95
96	if(trim(cnames(ic)).eq."C3H8") then
97	yc(l,ic)= 10.0*(7.217e0 - 994.30251e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
98	endif
99
100	if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
101	yc(l,ic)= 10.0**(96.26781e0 - 4651.872e0/TEMP(l) - 31.68595e0 &
102	alog10(TEMP(l)) ) / PRESS(l) 1013.25e0 / 760.0e0
103	endif
104
105	if(trim(cnames(ic)).eq."C4H4") then
106	yc(l,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(TEMP(l)-43.15e0) ) / PRESS(l)
107	endif
108
109	if(trim(cnames(ic)).eq."C4H6") then
110	yc(l,ic)= 10.0*(2.8808e0 - 4.6e0(262.3e0 - TEMP(l)) &
111	/(1.15e0TEMP(l) - 39.345e0) ) / PRESS(l) 1013.25e0 / 760.0e0
112	endif
113
114	if(trim(cnames(ic)).eq."C4H10") then
115	yc(l,ic)= 10.0*(8.446e0 - 1461.2e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
116	endif
117
118	if(trim(cnames(ic)).eq."C6H2") then
119	yc(l,ic)= 10.0*(4.666e0 - 4956e0/TEMP(l) + 25.845e0 &
120	alog10(1.0e3/TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
121	endif
122
123	if(trim(cnames(ic)).eq."C8H2") then
124	yc(l,ic)= 10.0*(3.95e0 - 6613e0/TEMP(l) + 35.055e0 &
125	alog10(1.0e3/TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
126	endif
127
128	if(trim(cnames(ic)).eq."AC6H6") then
129	x = 1.0e0 - TEMP(l) / 562.2e0
130	yc(l,ic)= 48.9e3 * exp( ( 1.33213 * x*1.5 - 6.98273 x &
131	- x*3 (2.62863 + 3.33399 * x*3) ) 562.2e0/TEMP(l) ) / PRESS(l)
132	endif
133
134	if(trim(cnames(ic)).eq."HCN") then
135	yc(l,ic)= 10.0*(8.6165e0 - 1516.5e0/(TEMP(l) - 26.2e0) ) / PRESS(l) 1013.25e0 / 760.0e0
136	endif
137
138	if(trim(cnames(ic)).eq."CH3CN") then
139	yc(l,ic)= 10.0*(8.458e0 - 1911.7e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
140	endif
141
142	if(trim(cnames(ic)).eq."C2H3CN") then
143	yc(l,ic)= 10.0*(9.3051e0 - 2782.21/(TEMP(l) - 51.15e0) ) / PRESS(l) 1013.25e0 / 760.0e0
144	endif
145
146	if(trim(cnames(ic)).eq."NCCN") then
147	yc(l,ic)= 10.0*(7.454e0 - 1832e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
148	endif
149
150	if(trim(cnames(ic)).eq."HC3N") then
151	yc(l,ic)= 10.0*(7.7446e0 - 1453.5609e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
152	endif
153
154	if(trim(cnames(ic)).eq."C4N2") then
155	yc(l,ic)= 10.0*(8.269e0 - 2155.0e0/TEMP(l) ) / PRESS(l) 1013.25e0 / 760.0e0
156	endif
157
158	enddo
159	enddo
160
161	end subroutine inicondens

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

Download in other formats:

Original Format