Context Navigation

radtitan.F @ 1243

Last change on this file since 1243 was 1058, checked in by slebonnois, 11 years ago
SL: goes with 1056-1057: additional details and bug corrections
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
115	C IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE
116	C PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG
117	IPRINT=1
118
119	C&&
120	FHAZE=0.3
121	C&&
122	if(iprem.eq.0) then
123	TAUFAC=0
124	c xvis et xir lus dans physiq.def (ancien fichier initpar)
125	FHVIS= xvis
126	FHIR = xir
127	c on initialise le paquet optcld
128	if (clouds.eq.1) call iniqcld()
129	iprem=1
130	endif
131
132	c-----------------------------------------------------------------------
133	c 2. Calcul of the atmospheric profile:
134	c -------------------------------------
135
136	print*,'dans radtitan ',klon
137	print*,notfirstcall
138	IF(notfirstcall) GOTO 300 !F au premier appel!
139	print*,notfirstcall
140
141	c pression moyenne globale
142	c passage au pressions en bar avec indice 1 au sommet.
143	c (similaire zp dans radlwsw)
144	DO 210 J=2,NLEVEL
145	PRESS(J)=plevmoy(NLEVEL+1-j)*1.e-5
146	210 CONTINUE
147	PRESS(1) = PRESS(2)*0.001
148
149	c a cause du tableau predefini dans lell.F (et lell_light.F)
150	c IF(press(nlevel-1).GE.1.44) then
151	IF(press(nlevel-1).GE.1.48) then
152	STOP'pression au sol trop grande'
153	PRINT*,'pression au sol trop grande'
154	endif
155
156	c PRESS(nlevel)=1.48
157	c XCORR=1.48/PRESS(nlevel)
158	c DO 211 J=1,NLEVEL
159	c PRESS(J)=XCORR*PRESS(J)
160	c11 CONTINUE
161
162	c *********************************************************
163	c + 20/1/00: S.Lebonnois: model with chemistry
164	c ++ 22/07/02: ajout HCN ++
165	c *********************************************************
166	if (ylellouch) then
167	c------------------------------------------------------
168	c initialisation de l'atmosphere et de la composition
169	c------------------------------------------------------
170	CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU,
171	& CH4,H2,XN2,AR,IPRINT)
172
173	print*,'LELLOUCH'
174	do i=1,55
175	print*,z(i),PRESS(i)
176	enddo
177	C
178	C
179	C NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
180	CALL GASSES(IPRINT)
181
182	else
183	c------------------------------------------------------
184	c initialisation seulement de l'atmosphere
185	c------------------------------------------------------
186	CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU,
187	& CH4,H2,XN2,AR,IPRINT)
188
189	print*,'LELLOUCH LIGHT'
190	do i=1,55
191	print*,z(i),PRESS(i)
192	enddo
193
194	c ++ remplace gasses.F ++
195
196	do i=1,nq
197	if (tname(i).eq."CH4") then
198	iradch4=i
199	elseif (tname(i).eq."C2H2") then
200	iradc2h2=i
201	elseif (tname(i).eq."C2H6") then
202	iradc2h6=i
203	elseif (tname(i).eq."HCN") then
204	iradhcn=i
205	elseif (tname(i).eq."N2") then
206	iradn2=i
207	elseif (tname(i).eq."H2") then
208	iradh2=i
209	endif
210	enddo
211
212	c print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2
213
214	print*,' ALT CH4 mass mixing ratio '
215
216	somcoslat=0.
217	do j=1,klon
218	coslat(j) = cos(rlatd(j)*RPI/180.)
219	somcoslat=somcoslat+coslat(j)
220	enddo
221	do i=1,nlayer
222	c attention ici, Z en km doit etre passe en m
223	colden(i)=rhop(press(i+1)-press(i))/effg(z(i)1000.)
224	gas1(i)=0.
225	emu=(xmu(i+1)+xmu(i))/2.
226	do j=1,klon
227	gas1(i) = gas1(i) +
228	$ coslat(j)/somcoslatycomp(j,i,iradch4)(16./emu)
229	enddo
230	print*,z(i),gas1(i)
231	enddo
232
233	RHCH4=0.
234	do j=1,klon
235	RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4)
236	enddo
237	RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel))
238	print*,'RHCH4 = ',RHCH4
239
240	endif
241
242	c *********************************************************
243
244	C
245	C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
246	CALL SETSPI(IPRINT)
247	CALL SETSPV(IPRINT)
248	C SET UP PIA COEFFICIENTS
249	CALL SETPIA(IPRINT,1)
250
251	IF (TAUFAC .GT. 0.) CALL CLD(IPRINT)
252
253	C
254	C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
255	C INFRARED. AND THEN IN THE VISIBLE.
256
257	C NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
258	C AND AT EACH CALL OF THE PHYSICS
259
260	print*,'aerosol/gas/cloud properties'
261
262	CALL OPTCI(ycomp,qaer,nmicro,IPRINT) ! #1
263	print*,'On sort de optci'
264
265	C NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
266	C INFRARED. AND THEN IN THE VISIBLE.
267
268	CALL OPTCV(qaer,nmicro,IPRINT) ! #2
269
270	do j=1,NLAYER
271	DZED(j)=Z(J)-Z(J+1)
272	enddo
273
274	c print*,wlnv
275	c print*,""
276	c print*,wlni
277	c stop
278
279	300 CONTINUE ! fin notfirstcall
280
281
282	c -----------------------------
283	c on ne recalcule pas optci si microfi=0 et compo lellouch
284	c -----------------------------
285	IF ((MICROFI.ge.1).or.(.not.ylellouch)) THEN
286	IF(notfirstcall) THEN !F au 1er appel T aux autres appels!!
287	print*,'aerosol/gas/cloud properties'
288	CALL OPTCI(ycomp,qaer,nmicro,IPRINT) ! #1
289	ENDIF
290	ENDIF
291
292	c ni optcv si microfi=0
293
294	IF (MICROFI.ge.1) THEN
295	IF(notfirstcall) THEN !F au 1er appel T aux autres appels!!
296	print*,'aerosol/gas/cloud properties'
297	CALL OPTCV(qaer,nmicro,IPRINT) ! #2
298	ENDIF
299	ENDIF
300
301	c -----------------------------
302	if (klon.eq.1) then
303	ig=1
304	else
305	ig=klon/2
306	endif
307	c print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1)
308	c print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10)
309	c print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI)
310	c print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2)
311	c print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2)
312	c print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2)
313	c stop
314
315	notfirstcall=.true.
316
317	RETURN
318	191 FORMAT(F8.2,1P10E10.2)
319	192 FORMAT(a8,1P10E10.2)
320	END

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

Download in other formats:

Original Format