1 | SUBROUTINE RADTITAN(p,nq,nmicro,ycomp) |
---|
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 | IMPLICIT NONE |
---|
37 | #include "dimensions.h" |
---|
38 | #include "clesphys.h" |
---|
39 | #include "microtab.h" |
---|
40 | #include "numchimrad.h" |
---|
41 | #include "YOMCST.h" |
---|
42 | |
---|
43 | c Pour le CRAY, les block data doivent etre declares external |
---|
44 | c pour etre pris en compte |
---|
45 | EXTERNAL TGMDAT |
---|
46 | |
---|
47 | INTEGER NLEVEL,NLAYER,NSPECI,NSPC1I,NSPECV,NSPC1V,NSPV |
---|
48 | PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1) |
---|
49 | PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25) |
---|
50 | PARAMETER (NSPV=21) ! LDO POUR CALCUL ALBEDO |
---|
51 | c |
---|
52 | |
---|
53 | c Arguments: |
---|
54 | c ---------- |
---|
55 | |
---|
56 | INTEGER nq,nmicro |
---|
57 | |
---|
58 | REAL p(klon,nlevel) |
---|
59 | REAL ycomp(klon,nlayer,nq) |
---|
60 | |
---|
61 | c Local: |
---|
62 | c ------ |
---|
63 | |
---|
64 | INTEGER I,J,IG,K,IPRINT |
---|
65 | INTEGER IPREM |
---|
66 | LOGICAL notfirstcall |
---|
67 | SAVE IPREM,notfirstcall |
---|
68 | data notfirstcall/.false./ |
---|
69 | |
---|
70 | REAL emu,somcoslat,coslat(klon) |
---|
71 | |
---|
72 | REAL PCH4, effg,FH2L,RHCH4L,SSUM ! effg est une fonction(z) |
---|
73 | |
---|
74 | c COMMONS for interface with local subroutines: |
---|
75 | c --------------------------------------------- |
---|
76 | |
---|
77 | REAL DTAUP(klon,NLAYER,NSPECI) |
---|
78 | REAL UBARI,UBARV,UBAR0 |
---|
79 | REAL DZED(NLAYER) |
---|
80 | REAL Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL) |
---|
81 | REAL DTDP(NLAYER),CONVEQ |
---|
82 | REAL CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL) |
---|
83 | REAL XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER) |
---|
84 | REAL C2H2(NLAYER),C2H6(NLAYER),HCN(NLAYER) |
---|
85 | REAL SOLARF(NSPECV),PEXPON(NSPECV),ATERM(4,NSPECV) |
---|
86 | INTEGER NTERM(NSPECV) |
---|
87 | REAL BTERM(4,NSPECV) |
---|
88 | REAL RADIUS(NLAYER), XNUMB(NLAYER),REALI(NSPECI), XIMGI(NSPECI) |
---|
89 | REAL REALV(NSPECV), XIMGV(NSPECV) |
---|
90 | REAL RADCLD(NLAYER), XNCLD(NLAYER),RCLDI(NSPECI), XICLDI(NSPECI) |
---|
91 | REAL RCLDV(NSPECV), XICLDV(NSPECV) |
---|
92 | REAL TAUHI(klon,NSPECI),TAUCI(klon,NSPECI) |
---|
93 | REAL TAUGI(klon,NSPECI), TAUGV(klon,NSPECV) |
---|
94 | REAL TAURV(klon,NSPECV),TAUHV(klon,NSPECV) |
---|
95 | REAL TAUCV(klon,NSPECV) |
---|
96 | c |
---|
97 | REAL DTAUI(klon,NLAYER,NSPECI) |
---|
98 | REAL TAUI(klon,NLEVEL,NSPECI) |
---|
99 | REAL WBARI(klon,NLAYER,NSPECI) |
---|
100 | REAL COSBI(klon,NLAYER,NSPECI) |
---|
101 | REAL BWNI(NSPC1I),WNOI(NSPECI) |
---|
102 | REAL WLNI(NSPECI),DWNI(NSPECI) |
---|
103 | c |
---|
104 | REAL DTAUV(klon,NLAYER,NSPECV,4) |
---|
105 | REAL TAUV(klon,NLEVEL,NSPECV,4) |
---|
106 | REAL WBARV(klon,NLAYER,NSPECV,4) |
---|
107 | REAL COSBV(klon,NLAYER,NSPECV,4) |
---|
108 | REAL BWNV(NSPC1V), WNOV(NSPECV),DWNV(NSPECV), WLNV(NSPECV) |
---|
109 | REAL FNETV(klon,NLEVEL),FUPV(klon,NLEVEL,NSPECV) |
---|
110 | REAL FDV(klon,NLEVEL,NSPECV),FMNETV(klon,NLEVEL) |
---|
111 | REAL FMUPV(NLEVEL),FMDV(NLEVEL) |
---|
112 | REAL FNET(klon,NLEVEL),FMNET(klon,NLEVEL) |
---|
113 | REAL THEAT(klon,NLAYER) |
---|
114 | REAL CSUBP,RSFI,RSFV,F0PI |
---|
115 | REAL RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON |
---|
116 | REAL RGAS,RHOP,PII,SIGMA |
---|
117 | REAL TIDAL |
---|
118 | |
---|
119 | COMMON /IRTAUS/ DTAUP |
---|
120 | |
---|
121 | COMMON /VERTICAL/ DZED |
---|
122 | |
---|
123 | COMMON /ATM/ Z,PRESS |
---|
124 | & ,DEN,TEMP |
---|
125 | |
---|
126 | COMMON /LAPSE/ DTDP,CONVEQ |
---|
127 | COMMON /UBARED/ UBARI,UBARV,UBAR0 |
---|
128 | |
---|
129 | |
---|
130 | |
---|
131 | COMMON /GASS/ CH4,XN2 |
---|
132 | & ,H2,AR |
---|
133 | & ,XMU,GAS1 |
---|
134 | & ,COLDEN |
---|
135 | |
---|
136 | COMMON /STRATO/ C2H2,C2H6 |
---|
137 | COMMON /STRAT2/ HCN |
---|
138 | |
---|
139 | COMMON /VISGAS/ SOLARF,NTERM |
---|
140 | & ,PEXPON,ATERM |
---|
141 | & ,BTERM |
---|
142 | |
---|
143 | COMMON /AERSOL/ RADIUS, XNUMB |
---|
144 | & ,REALI, XIMGI |
---|
145 | & ,REALV, XIMGV |
---|
146 | |
---|
147 | COMMON /CLOUD/ RADCLD, XNCLD |
---|
148 | & , RCLDI, XICLDI |
---|
149 | & , RCLDV, XICLDV |
---|
150 | |
---|
151 | COMMON /TAUS/ TAUHI,TAUCI, |
---|
152 | & TAUGI, TAUGV, |
---|
153 | & TAURV,TAUHV, |
---|
154 | & TAUCV |
---|
155 | |
---|
156 | * INFRARED CHARACTERISTICS |
---|
157 | *------------------------------ |
---|
158 | |
---|
159 | |
---|
160 | COMMON /SPECTI/ BWNI,WNOI |
---|
161 | & ,DWNI,WLNI |
---|
162 | |
---|
163 | COMMON /OPTICI/ DTAUI, |
---|
164 | & TAUI, |
---|
165 | & WBARI, |
---|
166 | & COSBI |
---|
167 | |
---|
168 | |
---|
169 | |
---|
170 | * VISIBLE CHARACTERISTICS |
---|
171 | *------------------------------ |
---|
172 | |
---|
173 | |
---|
174 | |
---|
175 | COMMON /OPTICV/ DTAUV |
---|
176 | & ,TAUV |
---|
177 | & ,WBARV |
---|
178 | & ,COSBV |
---|
179 | |
---|
180 | COMMON /SPECTV/ BWNV, WNOV |
---|
181 | & ,DWNV, WLNV |
---|
182 | |
---|
183 | COMMON /FLUXvV/ FNETV, |
---|
184 | & FUPV, |
---|
185 | & FDV, |
---|
186 | & FMNETV |
---|
187 | |
---|
188 | COMMON /FLUX/ FNET, FMNET |
---|
189 | & ,THEAT |
---|
190 | |
---|
191 | |
---|
192 | COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI |
---|
193 | COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON |
---|
194 | COMMON /CONST/RGAS,RHOP,PII,SIGMA |
---|
195 | COMMON /IO/ TIDAL |
---|
196 | |
---|
197 | * common relatifs aux aerosols |
---|
198 | * nrad dans microtab.h |
---|
199 | REAL qaer(klon,nlayer,nqtot),volume(nrad),rayon(nrad),vrat, |
---|
200 | & drayon(nrad),dvolume(nrad) |
---|
201 | |
---|
202 | common/traceurs/qaer |
---|
203 | common/part/volume,rayon,vrat, |
---|
204 | & drayon,dvolume |
---|
205 | c----------------------------------------------------------------------- |
---|
206 | c 1. Initialisations: |
---|
207 | c ------------------- |
---|
208 | |
---|
209 | REAL xpoub,kkk,xvis,xir |
---|
210 | |
---|
211 | |
---|
212 | C IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE |
---|
213 | C PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG |
---|
214 | IPRINT=1 |
---|
215 | |
---|
216 | C MODIFY ADJUSTABLE NUMBERS HERE -- NOT IN COMMON |
---|
217 | C&& |
---|
218 | FHAZE=0.3 |
---|
219 | C&& |
---|
220 | if(iprem.eq.0) then |
---|
221 | TAUFAC=0 |
---|
222 | FHVIS=2.0 |
---|
223 | FHIR=.2 |
---|
224 | print*,'ouverture du fichier initpar' |
---|
225 | open (unit=1,file='initpar') |
---|
226 | read(1,*) xpoub,kkk,xvis,xir |
---|
227 | close(1) |
---|
228 | FHVIS= xvis |
---|
229 | FHIR = xir |
---|
230 | print*,'ouverture du fichier initpar ok' |
---|
231 | print*,'DANS RADTITAN' |
---|
232 | print*,'-------------' |
---|
233 | print*,'FHVIS = ',FHVIS |
---|
234 | print*,'FHIR = ',FHIR |
---|
235 | iprem=1 |
---|
236 | endif |
---|
237 | |
---|
238 | c----------------------------------------------------------------------- |
---|
239 | c 2. Calcul of the atmospheric profile: |
---|
240 | c ------------------------------------- |
---|
241 | |
---|
242 | print*,'dans radtitan ',klon |
---|
243 | print*,notfirstcall |
---|
244 | IF(notfirstcall) GOTO 300 !F au premier appel! |
---|
245 | print*,notfirstcall |
---|
246 | |
---|
247 | DO 210 J=1,NLEVEL |
---|
248 | PRESS(J)=SSUM(klon,p(1,j),1)/FLOAT(klon) |
---|
249 | 210 CONTINUE |
---|
250 | |
---|
251 | IF(press(nlevel-1).GE.1.44) then |
---|
252 | STOP'pression au sol trop grande' |
---|
253 | PRINT*,'pression au sol trop grande' |
---|
254 | endif |
---|
255 | |
---|
256 | c PRESS(nlevel)=1.44 |
---|
257 | c XCORR=1.44/PRESS(nlevel) |
---|
258 | c DO 211 J=1,NLEVEL |
---|
259 | c PRESS(J)=XCORR*PRESS(J) |
---|
260 | c11 CONTINUE |
---|
261 | |
---|
262 | c ********************************************************* |
---|
263 | c + 20/1/00: S.Lebonnois: model with chemistry |
---|
264 | c ++ 22/07/02: ajout HCN ++ |
---|
265 | c ********************************************************* |
---|
266 | if (ylellouch) then |
---|
267 | c------------------------------------------------------ |
---|
268 | c initialisation de l'atmosphere et de la composition |
---|
269 | c------------------------------------------------------ |
---|
270 | CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU, |
---|
271 | & CH4,H2,XN2,AR,IPRINT) |
---|
272 | |
---|
273 | print*,'LELLOUCH' |
---|
274 | do i=1,55 |
---|
275 | print*,z(i),PRESS(i) |
---|
276 | enddo |
---|
277 | C |
---|
278 | C |
---|
279 | C NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS. |
---|
280 | CALL GASSES(IPRINT) |
---|
281 | |
---|
282 | else |
---|
283 | c------------------------------------------------------ |
---|
284 | c initialisation seulement de l'atmosphere |
---|
285 | c------------------------------------------------------ |
---|
286 | CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU, |
---|
287 | & CH4,H2,XN2,AR,IPRINT) |
---|
288 | |
---|
289 | print*,'LELLOUCH LIGHT' |
---|
290 | do i=1,55 |
---|
291 | print*,z(i),PRESS(i) |
---|
292 | enddo |
---|
293 | |
---|
294 | c ++ remplace gasses.F ++ |
---|
295 | |
---|
296 | do i=1,nq |
---|
297 | if (tname(i).eq."CH4") then |
---|
298 | iradch4=i |
---|
299 | elseif (tname(i).eq."C2H2") then |
---|
300 | iradc2h2=i |
---|
301 | elseif (tname(i).eq."C2H6") then |
---|
302 | iradc2h6=i |
---|
303 | elseif (tname(i).eq."HCN") then |
---|
304 | iradhcn=i |
---|
305 | elseif (tname(i).eq."N2") then |
---|
306 | iradn2=i |
---|
307 | elseif (tname(i).eq."H2") then |
---|
308 | iradh2=i |
---|
309 | endif |
---|
310 | enddo |
---|
311 | |
---|
312 | c print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2 |
---|
313 | |
---|
314 | print*,' ALT CH4 mass mixing ratio ' |
---|
315 | |
---|
316 | somcoslat=0. |
---|
317 | do j=1,klon |
---|
318 | coslat(j) = cos(rlatd(j)*RPI/180.) |
---|
319 | somcoslat=somcoslat+coslat(j) |
---|
320 | enddo |
---|
321 | do i=1,nlayer |
---|
322 | colden(i)=rhop*(press(i+1)-press(i))/effg(z(i)) |
---|
323 | gas1(i)=0. |
---|
324 | emu=(xmu(i+1)+xmu(i))/2. |
---|
325 | do j=1,klon |
---|
326 | gas1(i) = gas1(i) + |
---|
327 | $ coslat(j)/somcoslat*ycomp(j,i,iradch4)*(16./emu) |
---|
328 | enddo |
---|
329 | print*,z(i),gas1(i) |
---|
330 | enddo |
---|
331 | |
---|
332 | RHCH4=0. |
---|
333 | do j=1,klon |
---|
334 | RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4) |
---|
335 | enddo |
---|
336 | RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel)) |
---|
337 | print*,'RHCH4 = ',RHCH4 |
---|
338 | |
---|
339 | endif |
---|
340 | |
---|
341 | c ********************************************************* |
---|
342 | |
---|
343 | C |
---|
344 | C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS |
---|
345 | CALL SETSPI(IPRINT) |
---|
346 | CALL SETSPV(IPRINT) |
---|
347 | C SET UP PIA COEFFICIENTS |
---|
348 | CALL SETPIA(IPRINT,1) |
---|
349 | |
---|
350 | IF (TAUFAC .GT. 0.) CALL CLD(IPRINT) |
---|
351 | |
---|
352 | C |
---|
353 | C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE |
---|
354 | C INFRARED. AND THEN IN THE VISIBLE. |
---|
355 | |
---|
356 | C NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon |
---|
357 | C AND AT EACH CALL OF THE PHYSICS |
---|
358 | |
---|
359 | print*,'aerosol/gas/cloud properties' |
---|
360 | |
---|
361 | CALL OPTCI(ycomp,nmicro,IPRINT) ! #1 |
---|
362 | print*,'On sort de optci' |
---|
363 | |
---|
364 | C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED |
---|
365 | DO 225 IG=1,klon |
---|
366 | DO 220 J=1,NLAYER |
---|
367 | DO 230 K=1,NSPECI |
---|
368 | DTAUP(IG,J,K)=DTAUI(IG,J,K) |
---|
369 | 230 CONTINUE |
---|
370 | 220 CONTINUE |
---|
371 | 225 CONTINUE |
---|
372 | |
---|
373 | C NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon |
---|
374 | C INFRARED. AND THEN IN THE VISIBLE. |
---|
375 | |
---|
376 | CALL OPTCV(nmicro,IPRINT) ! #2 |
---|
377 | |
---|
378 | do j=1,NLAYER |
---|
379 | DZED(j)=Z(J)-Z(J+1) |
---|
380 | enddo |
---|
381 | |
---|
382 | c print*,wlnv |
---|
383 | c print*,"" |
---|
384 | c print*,wlni |
---|
385 | c stop |
---|
386 | |
---|
387 | 300 CONTINUE ! fin notfirstcall |
---|
388 | |
---|
389 | |
---|
390 | c ----------------------------- |
---|
391 | c on ne recalcule pas optci si microfi=0 et compo lellouch |
---|
392 | c ----------------------------- |
---|
393 | IF ((MICROFI.eq.1).or.(.not.ylellouch)) THEN |
---|
394 | IF(notfirstcall) THEN !F au 1er appel T aux autres appels!! |
---|
395 | print*,'aerosol/gas/cloud properties' |
---|
396 | CALL OPTCI(ycomp,nmicro,IPRINT) ! #1 |
---|
397 | DO IG=1,klon |
---|
398 | DO J=1,NLAYER |
---|
399 | DO K=1,NSPECI |
---|
400 | DTAUP(IG,J,K)=DTAUI(IG,J,K) |
---|
401 | ENDDO |
---|
402 | ENDDO |
---|
403 | ENDDO |
---|
404 | ENDIF |
---|
405 | ENDIF |
---|
406 | |
---|
407 | c ni optcv si microfi=0 |
---|
408 | |
---|
409 | IF (MICROFI.eq.1) THEN |
---|
410 | IF(notfirstcall) THEN !F au 1er appel T aux autres appels!! |
---|
411 | print*,'aerosol/gas/cloud properties' |
---|
412 | CALL OPTCV(nmicro,IPRINT) ! #2 |
---|
413 | ENDIF |
---|
414 | ENDIF |
---|
415 | |
---|
416 | c ----------------------------- |
---|
417 | if (klon.eq.1) then |
---|
418 | ig=1 |
---|
419 | else |
---|
420 | ig=klon/2 |
---|
421 | endif |
---|
422 | c print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1) |
---|
423 | c print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10) |
---|
424 | c print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI) |
---|
425 | c print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2) |
---|
426 | c print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2) |
---|
427 | c print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2) |
---|
428 | c stop |
---|
429 | |
---|
430 | notfirstcall=.true. |
---|
431 | |
---|
432 | RETURN |
---|
433 | 191 FORMAT(F8.2,1P10E10.2) |
---|
434 | 192 FORMAT(a8,1P10E10.2) |
---|
435 | END |
---|