Context Navigation

radtitan.F @ 855

Last change on this file since 855 was 815, checked in by slebonnois, 12 years ago
SL: petites modifs Titan et Venus pour tableau controle dans la physique ; pour Titan, petits details lies a raz_date ; modif chemin ioipsl sur gnome ; + elimination d'un warning etrange dans gcm.F
File size: 8.9 KB

Line
1	SUBROUTINE RADTITAN(p,nq,nmicro,ycomp,qaer)
2
3	c=======================================================================
4	c
5	c Authors: C.P. Mc Kay 01/02/91
6	c -------
7	c
8	c Object: Computation of the solar and infra-red
9	c ------- Opacities (dans des common...)
10	c
11	c ON TITAN ADAPTED FROM BEST.FOR FEB 91
12	c C.P. McKAY
13	c
14	c Arguments:
15	c ----------
16	c
17	c Input:
18	c ------
19	c
20	c p(klon,nl) pressure (level)
21	c nq nombre de traceurs
22	c nmicro nombre de traceurs microphysiques
23	c ycomp(klon,nlayer,nq)
24	c
25	c Output:
26	c -------
27	c
28	c=======================================================================
29	c-----------------------------------------------------------------------
30	c Declarations:
31	c -------------
32
33	USE infotrac
34	use dimphy
35	USE comgeomphy
36	USE optcld, only : iniqcld
37	IMPLICIT NONE
38	#include "dimensions.h"
39	#include "clesphys.h"
40	#include "microtab.h"
41	#include "numchimrad.h"
42	#include "YOMCST.h"
43
44	c Pour le CRAY, les block data doivent etre declares external
45	c pour etre pris en compte
46	EXTERNAL TGMDAT
47
48	INTEGER NLEVEL,NLAYER,NSPECI,NSPC1I,NSPECV,NSPC1V,NSPV
49	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
50	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
51	PARAMETER (NSPV=21) ! LDO POUR CALCUL ALBEDO
52
53	c
54	c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
55	INTEGER ngrid
56	PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon
57	c
58
59	c Arguments:
60	c ----------
61
62	INTEGER nq,nmicro
63
64	REAL p(klon,nlevel)
65	REAL ycomp(klon,nlayer,nq)
66	REAL qaer(klon,klev,nq)
67
68	c Local:
69	c ------
70
71	INTEGER I,J,IG,K,IPRINT
72	INTEGER IPREM
73	LOGICAL notfirstcall
74	SAVE IPREM,notfirstcall
75	data notfirstcall/.false./
76
77	REAL emu,somcoslat,coslat(ngrid)
78
79	REAL PCH4, effg,FH2L,RHCH4L,SSUM ! effg est une fonction(z)
80
81	c COMMONS for interface with local subroutines:
82	c ---------------------------------------------
83
84	REAL DZED(NLAYER)
85	REAL Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
86	REAL CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
87	REAL XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
88	REAL C2H2(NLAYER),C2H6(NLAYER),HCN(NLAYER)
89	REAL RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
90	REAL RGAS,RHOP,PI,SIGMA
91
92	COMMON /VERTICAL/ DZED
93
94	COMMON /ATM/ Z,PRESS
95	& ,DEN,TEMP
96
97
98	COMMON /GASS/ CH4,XN2
99	& ,H2,AR
100	& ,XMU,GAS1
101	& ,COLDEN
102
103	COMMON /STRATO/ C2H2,C2H6
104	COMMON /STRAT2/ HCN
105
106	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
107	COMMON /CONST/RGAS,RHOP,PI,SIGMA
108
109	c-----------------------------------------------------------------------
110	c 1. Initialisations:
111	c -------------------
112
113	REAL xpoub,kkk,xvis,xir
114
115
116	C IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE
117	C PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG
118	IPRINT=1
119
120	C MODIFY ADJUSTABLE NUMBERS HERE -- NOT IN COMMON
121	C&&
122	FHAZE=0.3
123	C&&
124	if(iprem.eq.0) then
125	TAUFAC=0
126	FHVIS=2.0
127	FHIR=.2
128	print*,'ouverture du fichier initpar'
129	open (unit=1,file='initpar')
130	read(1,*) xpoub,kkk,xvis,xir
131	close(1)
132	FHVIS= xvis
133	FHIR = xir
134	print*,'ouverture du fichier initpar ok'
135	print*,'DANS RADTITAN'
136	print*,'-------------'
137	print*,'FHVIS = ',FHVIS
138	print*,'FHIR = ',FHIR
139	c on initialise le paquet optcld
140	if (clouds.eq.1) call iniqcld()
141	iprem=1
142	endif
143
144	c-----------------------------------------------------------------------
145	c 2. Calcul of the atmospheric profile:
146	c -------------------------------------
147
148	print*,'dans radtitan ',klon
149	print*,notfirstcall
150	IF(notfirstcall) GOTO 300 !F au premier appel!
151	print*,notfirstcall
152
153	DO 210 J=1,NLEVEL
154	PRESS(J)=SSUM(klon,p(1,j),1)/FLOAT(klon)
155	210 CONTINUE
156
157	c a cause du tableau predefini dans lell.F (et lell_light.F)
158	c IF(press(nlevel-1).GE.1.44) then
159	IF(press(nlevel-1).GE.1.48) then
160	STOP'pression au sol trop grande'
161	PRINT*,'pression au sol trop grande'
162	endif
163
164	c PRESS(nlevel)=1.48
165	c XCORR=1.48/PRESS(nlevel)
166	c DO 211 J=1,NLEVEL
167	c PRESS(J)=XCORR*PRESS(J)
168	c11 CONTINUE
169
170	c *********************************************************
171	c + 20/1/00: S.Lebonnois: model with chemistry
172	c ++ 22/07/02: ajout HCN ++
173	c *********************************************************
174	if (ylellouch) then
175	c------------------------------------------------------
176	c initialisation de l'atmosphere et de la composition
177	c------------------------------------------------------
178	CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU,
179	& CH4,H2,XN2,AR,IPRINT)
180
181	print*,'LELLOUCH'
182	do i=1,55
183	print*,z(i),PRESS(i)
184	enddo
185	C
186	C
187	C NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
188	CALL GASSES(IPRINT)
189
190	else
191	c------------------------------------------------------
192	c initialisation seulement de l'atmosphere
193	c------------------------------------------------------
194	CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU,
195	& CH4,H2,XN2,AR,IPRINT)
196
197	print*,'LELLOUCH LIGHT'
198	do i=1,55
199	print*,z(i),PRESS(i)
200	enddo
201
202	c ++ remplace gasses.F ++
203
204	do i=1,nq
205	if (tname(i).eq."CH4") then
206	iradch4=i
207	elseif (tname(i).eq."C2H2") then
208	iradc2h2=i
209	elseif (tname(i).eq."C2H6") then
210	iradc2h6=i
211	elseif (tname(i).eq."HCN") then
212	iradhcn=i
213	elseif (tname(i).eq."N2") then
214	iradn2=i
215	elseif (tname(i).eq."H2") then
216	iradh2=i
217	endif
218	enddo
219
220	c print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2
221
222	print*,' ALT CH4 mass mixing ratio '
223
224	somcoslat=0.
225	do j=1,klon
226	coslat(j) = cos(rlatd(j)*RPI/180.)
227	somcoslat=somcoslat+coslat(j)
228	enddo
229	do i=1,nlayer
230	c attention ici, Z en km doit etre passe en m
231	colden(i)=rhop(press(i+1)-press(i))/effg(z(i)1000.)
232	gas1(i)=0.
233	emu=(xmu(i+1)+xmu(i))/2.
234	do j=1,klon
235	gas1(i) = gas1(i) +
236	$ coslat(j)/somcoslatycomp(j,i,iradch4)(16./emu)
237	enddo
238	print*,z(i),gas1(i)
239	enddo
240
241	RHCH4=0.
242	do j=1,klon
243	RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4)
244	enddo
245	RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel))
246	print*,'RHCH4 = ',RHCH4
247
248	endif
249
250	c *********************************************************
251
252	C
253	C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
254	CALL SETSPI(IPRINT)
255	CALL SETSPV(IPRINT)
256	C SET UP PIA COEFFICIENTS
257	CALL SETPIA(IPRINT,1)
258
259	IF (TAUFAC .GT. 0.) CALL CLD(IPRINT)
260
261	C
262	C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
263	C INFRARED. AND THEN IN THE VISIBLE.
264
265	C NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
266	C AND AT EACH CALL OF THE PHYSICS
267
268	print*,'aerosol/gas/cloud properties'
269
270	CALL OPTCI(ycomp,qaer,nmicro,IPRINT) ! #1
271	print*,'On sort de optci'
272
273	C NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
274	C INFRARED. AND THEN IN THE VISIBLE.
275
276	CALL OPTCV(qaer,nmicro,IPRINT) ! #2
277
278	do j=1,NLAYER
279	DZED(j)=Z(J)-Z(J+1)
280	enddo
281
282	c print*,wlnv
283	c print*,""
284	c print*,wlni
285	c stop
286
287	300 CONTINUE ! fin notfirstcall
288
289
290	c -----------------------------
291	c on ne recalcule pas optci si microfi=0 et compo lellouch
292	c -----------------------------
293	IF ((MICROFI.ge.1).or.(.not.ylellouch)) THEN
294	IF(notfirstcall) THEN !F au 1er appel T aux autres appels!!
295	print*,'aerosol/gas/cloud properties'
296	CALL OPTCI(ycomp,qaer,nmicro,IPRINT) ! #1
297	ENDIF
298	ENDIF
299
300	c ni optcv si microfi=0
301
302	IF (MICROFI.ge.1) THEN
303	IF(notfirstcall) THEN !F au 1er appel T aux autres appels!!
304	print*,'aerosol/gas/cloud properties'
305	CALL OPTCV(qaer,nmicro,IPRINT) ! #2
306	ENDIF
307	ENDIF
308
309	c -----------------------------
310	if (klon.eq.1) then
311	ig=1
312	else
313	ig=klon/2
314	endif
315	c print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1)
316	c print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10)
317	c print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI)
318	c print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2)
319	c print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2)
320	c print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2)
321	c stop
322
323	notfirstcall=.true.
324
325	RETURN
326	191 FORMAT(F8.2,1P10E10.2)
327	192 FORMAT(a8,1P10E10.2)
328	END

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

Download in other formats:

Original Format