Context Navigation

phytrac_chimie.F @ 1395

Last change on this file since 1395 was 1305, checked in by slebonnois, 11 years ago
SL: VENUS PHOTOCHEMISTRY. Needs Lapack (see arch files...)
Property svn:executable set to ``*
File size: 11.8 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
3	!
4	c
5	c
6	SUBROUTINE phytrac_chimie (
7	I debutphy,
8	I gmtime,
9	I nqmax,
10	I n_lon,
11	I lat,
12	I lon,
13	I n_lev,
14	I pdtphys,
15	I temp,
16	I pplev,
17	O trac,
18	O NBRTOT_droplet,
19	O W_H2SO4,
20	O rho)
21
22	c======================================================================
23	c Auteur(s) FH
24	c Objet: Moniteur general des tendances traceurs
25	c
26	cAA Remarques en vrac:
27	cAA--------------------
28	cAA 1/ le call phytrac se fait avec nqmax
29	c======================================================================
30	USE chemparam_mod
31	IMPLICIT none
32
33	#include "clesphys.h"
34	#include "YOMCST.h"
35	c======================================================================
36	c Arguments:
37	c
38	c
39	c EN ENTREE:
40	c ==========
41	c
42	c divers:
43	c -------
44	c
45
46	REAL sza_local
47	REAL gmtime
48	c INTEGER, SAVE :: cpt_cloudIO !un compteur pour fichier sortie cloud_parameter en 1D
49	INTEGER iq
50	INTEGER i
51	INTEGER ilon, ilev
52	integer n_lon ! nombre de points horizontaux
53	INTEGER n_lev ! nombre de couches verticales
54	INTEGER nqmax ! nombre de traceurs auxquels on applique la physique
55	c INTEGER nbapp_cloud, i_app_cloud
56	real pdtphys ! pas d'integration pour la physique (seconde)
57	real lat(n_lon), lat_local(n_lon)
58	real lon(n_lon), lon_local(n_lon)
59	real temp(n_lon,n_lev) ! temp
60	real trac(n_lon,n_lev,nqmax) ! traceur
61	real pplev(n_lon,n_lev) ! pression pour le mileu de chaque couche (en Pa)
62	real lon_sun
63	logical debutphy ! le flag de l'initialisation de la physique
64	c character*7 modname
65
66	C
67	C----------------------------------------------------------------------------
68	C Model cloud:
69	C Aerosol and PSC variables:
70	real NBRTOT_droplet(n_lon,n_lev)
71	real W_H2SO4(n_lon,n_lev)
72	real W_H2O(n_lon,n_lev)
73	real rho(n_lon,n_lev)
74	C----------------------------------------------------------------------------
75	C----------------------------------------------------------------------------
76	C Time variables:
77	REAL, save :: dT_cloud
78	C----------------------------------------------------------------------------
79	C----------------------------------------------------------------------------
80	C Auxilary variables:
81
82	REAL mrtwv,mrtsa,mrwv,mrsa,
83	+ ppwv, psatwv,
84	+ ps_sa,satps_sa
85	C ps_sa: satur pressure pure SA
86	C satps_sa: satur pres over mixture in dyne/cm2=Pa/10
87	C----------------------------------------------------------------------------
88
89	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
90
91	logical ok_sync
92	parameter (ok_sync = .true.)
93
94	c modname = 'phytrac'
95
96	c PRINT*,'DEBUT subroutine PHYTRAC'
97
98	c----------------------------------
99	c debutphy: Initiation des traceurs
100	c----------------------------------
101
102	if (debutphy) then
103
104	PRINT*,'PRECISION REAL'
105	PRINT*,precision(NBRTOT_droplet(1,1)), range(NBRTOT_droplet(1,1))
106
107	if (n_lon .EQ. 1) then
108	PRINT*,'n_lon 1D: ',n_lon
109	end if
110
111	c if ((n_lon .GT. 1) .AND. ok_chem) then
112	c !!! DONC 3D !!!
113	c CALL chemparam_ini()
114	c endif
115
116	c if ((n_lon .GT. 1) .AND. ok_cloud) then
117	c !!! DONC 3D !!!
118	c CALL cloud_ini(n_lon,n_lev)
119	c endif
120
121	IF (reinit_trac) THEN
122	PRINT*,'REINIT MIXING RATIO TRACEURS'
123
124	c =============================================================
125	c Passage de Rm à Rv
126	c =============================================================
127	c Necessaire si on reprend les start.nc qui sont en MMR
128	DO iq=1,nqmax
129	trac(:,:,iq)=trac(:,:,iq)*RMD/M_tr(iq)
130	END DO
131	c =============================================================
132
133
134	c=============================================================
135	c Initialisation des profils traceurs en Rv
136	c=============================================================
137	trac(:,:,:)=1.0d-30
138
139	trac(:,:,i_co2)=0.965d0 * RMD / M_tr(i_co2)
140
141	trac(:,:,i_co)=25.0d-6
142
143
144	trac(:,:,i_h2so4)=1.0d-21
145	trac(:,:,i_h2o)=1.0d-21
146
147	c !!! SANS NUAGE !!!
148	c trac(:,1:29,i_ocs)=1.0d-6
149	c trac(:,29:40,i_ocs)=1.0e-9
150	c trac(:,:,i_so2)=1.d-6
151	c trac(:,:,i_h2o)=1.0d-6
152
153	c !!! AVEC NUAGE !!!
154	trac(:,1:20,i_ocs)=3.d-6
155
156	DO i=21,26
157	trac(:,i,i_ocs)=trac(:,i-1,i_ocs)-0.3d-6
158	END DO
159
160	c DO i=21,30
161	c trac(:,i,i_ocs)=trac(:,i-1,i_ocs)-0.3d-6
162	c END DO
163
164	trac(:,1:26,i_hcl)=0.2d-6
165
166	c trac(:,:,i_hcl)=0.2d-6
167
168	c Initialisation SO2 Bertaux et De Bergh 2007 JGR
169	c trac(:,1:26,i_so2)=20.d-6
170	c DO i=2,20
171	c trac(:,i,i_so2)=trac(:,i-1,i_so2)+(100./19.)*0.25d-6
172	c END DO
173	c DO i=21,22
174	c trac(:,i,i_so2)=trac(:,i-1,i_so2)-(100./9.)*0.25d-6
175	c END DO
176
177	c DO i=21,29
178	c trac(:,i,i_so2)=trac(:,i-1,i_so2)-(100./9.)*1d-6
179	c END DO
180
181	c trac(:,1:30,i_so2)= 100.0d-6
182	c trac(:,30,i_so2)=20.0d-6
183	c trac(:,31,i_so2)=10.0d-6
184	c trac(:,32,i_so2)=1.0d-6
185	c trac(:,33,i_so2)=0.1d-6
186	c trac(:,34:42,i_so2)=0.02d-6
187	c trac(:,43:46,i_so2)=0.07d-6
188	c trac(:,47:50,i_so2)=0.05d-6
189
190	c trac(:,1:28,i_h2o)=30.0d-6
191	c trac(:,29:50,i_h2o)=5.0d-6
192	trac(:,15:50,i_h2o)=10.0d-6
193	c trac(:,15:35,i_h2so4)=17.0d-6
194	c DO i=23,35
195	c trac(:,i,i_h2o)=(3.d-6-30.0d-6)/12.0*(-23.0+i)+trac(:,22,i_h2o)
196	c END DO
197	c trac(:,36:50,i_h2o)=3.0d-6
198
199	trac(:,15:50,i_h2so4)=20.0d-6
200	c trac(:,29:50,i_h2so4)=1.0d-9
201	c trac(:,1:10,i_h2)=1.0d-10
202	c trac(:,11:20,i_h2)=1.0d-9
203	c trac(:,21:35,i_h2)=1.0d-8
204	c trac(:,36:50,i_h2)=1.0d-7
205
206	c=============================================================
207
208	c =============================================================
209	c Passage de Rv à Rm
210	c =============================================================
211	DO iq=1,nqmax
212	trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/RMD
213	END DO
214	c =============================================================
215
216	c Ecriture fichier initialisation
217	c PRINT*,'Ecriture Initial_State.csv'
218	c OPEN(88,file='Initial_State.csv',
219	c & form='formatted')
220
221	c DO ilon=1,n_lon
222	c DO ilev=1,n_lev
223	c WRITE(88,"(36(e15.8,','))") R_MEDIAN(ilon,ilev),
224	c & STDDEV(ilon,ilev),trac(ilon,ilev,1:nqmax)
225	c ENDDO
226	c ENDDO
227	c PRINT*,'FIN Ecriture Initial_State.csv'
228
229	ENDIF !FIN REINIT TRAC
230
231	c-------------
232	c fin debutphy
233	c-------------
234	ENDIF ! fin debutphy
235
236	c =============================================================
237	c Passage de Rm à Rv
238	c =============================================================
239	DO iq=1,nqmax
240	trac(:,:,iq)=trac(:,:,iq)*RMD/M_tr(iq)
241	END DO
242	c =============================================================
243
244
245	c=============================================================
246	c Boucle sur les lon, lat (n_lon)
247	c=============================================================
248	c PRINT, 'gmtime', gmtimeRDAY
249	c PRINT*, 'RDAY', RDAY
250
251	lon_sun = (0.5 - gmtime) * 2.0 * RPI
252	lon_local = lon * RPI/180.0d0
253	lat_local = lat * RPI/180.0d0
254
255	DO ilon=1, n_lon
256
257	c calcul sza_local pour obtenir des sza_local > 90, utile pour la chimie
258	sza_local = acos(cos(lat_local(ilon))cos(lon_local(ilon))
259	& cos(lon_sun) + cos(lat_local(ilon))*sin(lon_local(ilon))
260	& sin(lon_sun)) 180.0d0/RPI
261
262	c PRINT*,'sza_local :', sza_local
263
264	IF (ok_cloud) THEN
265	c PRINT*,'DEBUT CLOUD'
266
267	dT_cloud=pdtphys
268
269
270	c nbapp_cloud=NINT(pdtphys/dT_cloud)
271	c PRINT*,'pdtphys',pdtphys
272	c PRINT*,'nbapp_cloud',nbapp_cloud
273	c =============================================================
274	c Appel Microphysique (sans nucleation)
275	c Volume Mixing Ratio
276	c =============================================================
277
278	c FIXE un profil de temperature def dans fichier temp
279	if (n_lon .EQ. 1) then
280	OPEN(13,file='temp',status='old',form='formatted')
281	DO ilev=1,n_lev
282	READ (13,*) temp(n_lon,ilev)
283	ENDDO
284	CLOSE(13)
285	endif
286
287	DO ilev=1, n_lev
288	c PRINT*,'DEBUT INIT CALL CLOUD'
289	c ppwv et pplev en Pa
290
291	c PRINT*,'@@@@ IN CLOUD @@@@'
292
293	c On remet tout le RM liq dans la partie gaz
294	c !!! On reforme un nuage à chaque fois !!!
295
296	mrtwv=trac(ilon,ilev,i_h2o) + trac(ilon,ilev,i_h2oliq)
297	mrtsa=trac(ilon,ilev,i_h2so4) + trac(ilon,ilev,i_h2so4liq)
298	mrwv=mrtwv
299	mrsa=mrtsa
300
301
302	c !!! Remise a zero !!!
303	W_H2SO4(ilon,ilev)=0.0d0
304	W_H2O(ilon,ilev)=0.d0
305	rho(ilon,ilev)=0.0d0
306	NBRTOT_droplet(ilon,ilev)=0.d0
307	satps_sa=0.d0
308	ps_sa=0.d0
309
310	c pression partielle H2O
311	ppwv=pplev(ilon,ilev) * mrwv
312
313	c Pression saturante de vapeur d'eau, tirée du code d'Anni
314	psatwv=EXP(77.344913 - 7235.4247/temp(ilon,ilev)
315	& - 8.2DLOG(temp(ilon,ilev)) + 0.0057113temp(ilon,ilev))
316
317	c PRINT*,'DEBUT CALL CLOUD'
318
319	c Ne pas passer par la routine des nuages si on a des valeurs proches de 0 ?
320	c Empeche de foirer en parallèle ?
321
322
323	CALL new_cloud_venus(dT_cloud,
324	e NBRTOT_droplet(ilon,ilev),
325	e R_MEDIAN(ilon,ilev),STDDEV(ilon,ilev),
326	e temp(ilon,ilev),pplev(ilon,ilev),
327	e ppwv,
328	e mrwv,mrsa,
329	e ilev,
330	e mrtwv,mrtsa,
331	e W_H2SO4(ilon,ilev),
332	e ps_sa,satps_sa,
333	e rho(ilon,ilev))
334
335	c END DO
336
337	c =========================================
338	c Actualisation des mixing ratio liq et gaz
339	c =========================================
340	c Si la routine new_cloud_venus n'a pas actualisé mrwv et mrsa
341	c on a alors bien mr=mrt pour sa et wv, donc les parties liq sont=0 hors du nuage
342	c ou si on ne condense pas
343
344	c PRINT*,'DEBUT ACTUALISATION OUTPUT CLOUD'
345	c si tout se passe bien, mrtwv et mrtsa ne changent pas
346
347	trac(ilon,ilev,i_h2o) = mrwv
348	trac(ilon,ilev,i_h2oliq) = mrtwv - trac(ilon,ilev,i_h2o)
349
350	trac(ilon,ilev,i_h2so4) = mrsa
351	trac(ilon,ilev,i_h2so4liq) = mrtsa - trac(ilon,ilev,i_h2so4)
352
353	c ENDIF
354
355
356	IF (n_lon .EQ. 1) THEN
357	WRITE(66,"(i4,','11(e15.8,','))") ilev,temp(ilon,ilev),
358	& pplev(ilon,ilev),ps_sa,satps_sa,NBRTOT_droplet(ilon,ilev),
359	& W_H2SO4(ilon,ilev),trac(ilon,ilev,i_h2oliq),
360	& trac(ilon,ilev,i_h2so4liq),mrwv,mrsa,trac(ilon,ilev,i_so2)
361	ENDIF
362
363	END DO
364
365	c =============================================================
366	c PRINT*,'FIN CLOUD'
367	ENDIF
368
369	IF (ok_chem) THEN
370	c PRINT*,"vmr SO2 ht atmo: ",trac(1,50,i_so2)
371	c PRINT*,'DEBUT CHEMISTRY'
372	c =============================================================
373	c Appel Photochimie
374	c =============================================================
375	c Pression en hPa => pplev/100.
376
377	CALL new_photochemistry_venus(n_lev, n_lon, pdtphys,
378	e pplev(ilon,:)/100.,
379	e temp(ilon,:),
380	e trac(ilon,:,:),
381	e sza_local, nqmax)
382	c =============================================================
383	c PRINT*,'FIN CHEMISTRY'
384	c PRINT*,"vmr SO2 ht atmo: ",trac(1,50,i_so2)
385
386	END IF
387
388	END DO
389	c =============================================================
390	c Passage de Rv à Rm
391	c =============================================================
392	DO iq=1,nqmax
393	trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/RMD
394	END DO
395	c =============================================================
396
397	c PRINT*,'FIN PHYTRAC'
398	RETURN
399	END

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

Download in other formats:

Original Format