Context Navigation

radtitan.F @ 1056

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

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

Download in other formats:

Original Format