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