Context Navigation

optcv_1pt.F @ 1243

Last change on this file since 1243 was 495, checked in by slebonnois, 13 years ago
Mise a jour physique Titan, ajout des forces de marees (dans la dynamique, sous flag titan). SL.
File size: 13.5 KB

Line
1	SUBROUTINE optcv_1pt(zqaer_1pt,rcdb,xfrb,ioptv,IPRINT)
2
3
4	use dimphy
5	#include "dimensions.h"
6	#include "microtab.h"
7	#include "clesphys.h"
8
9	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
10	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
11
12	c Arguments:
13	c ---------
14	integer IPRINT,ioptv
15	C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
16	* nrad dans microtab.h
17	real zqaer_1pt(NLAYER,2*nrad)
18	#include "optcv_1pt.h"
19	c ---------
20
21	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
22
23	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
24	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
25
26	COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
27	& ATERM(4,NSPECV),BTERM(4,NSPECV)
28
29	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
30	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
31
32	COMMON /CLOUD/
33	& RCLDI(NSPECI), XICLDI(NSPECI)
34	& , RCLDV(NSPECV), XICLDV(NSPECV)
35	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
36	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
37
38	COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
39	& ,DWNV(NSPECV),WLNV(NSPECV)
40
41	COMMON /PLANT/ CSUBP,F0PI
42	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
43	COMMON /CONST/ RGAS,RHOP,PI,SIGMA
44	* nrad dans microtab.h
45	COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
46
47	REAL QF1(nrad,NSPECV),QF2(nrad,NSPECV)
48	REAL QF3(nrad,NSPECV),QF4(nrad,NSPECV)
49	REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
50	REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
51	REAL QC1(nrad,NSPECV),QC2(nrad,NSPECV)
52	REAL QC3(nrad,NSPECV),QC4(nrad,NSPECV)
53
54	c---- NUAGES
55	real TNUABS,TNUSCAT
56	real rcdb(NLAYER)
57	real xfrb(NLAYER,4)
58
59	save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
60
61	integer ilat,jalt
62	common/toto/ilat,jalt
63
64	C*
65	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
66	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
67	C LAYER: WBAR, DTAU, COSBAR
68	C LEVEL: TAU
69	C
70	C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
71	C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
72	C MAY NOT BE ZERO.
73
74	c***** DEBUT DES BOUCLES **********************
75	DO 100 K=1,NSPECV !b! BOUCLE SUR LAMBDA
76
77	TAURV_1pt(K)=0.
78	TAUHV_1pt(K)=0. ! INTEGRATED TAU.......INITIALIZATION.
79	TAUCV_1pt(K)=0. ! Rayleigh, Haze, Cloud, Gas
80	TAUGV_1pt(K)=0. ! sca, abs, abs , abs
81
82	DO 100 J=1,NLAYER !a! BOUCLE SUR L"ALTITUDE
83	jalt=j
84	C #1: HAZE
85	c---------------------------
86
87	c CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
88	c USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um
89
90
91
92
93	c /\
94	c / \
95	c / \
96	c / _O \
97	c / \|/ \
98	c / / \ \
99	c / \|\ \/\ \
100	c / \|\| / \ \
101	c ----------------
102	c \| WARNING \|
103	c \| SLOW DOWN \|
104	c ----------------
105
106
107
108
109	c********* EN TRAVAUX *************************
110
111	TAEROS=0.
112	TAEROSCAT=0.
113	CBAR=0.
114
115	c print*,"rayon=",rayon
116	c print*,"RF=",RF
117
118	DO inq=1,nrad !BOUCLE SUR LES TAILLE D"AEROSOLS
119
120
121	IF (rayon(inq).lt.RF(inq)) THEN ! aerosols spheriques
122
123
124	if(ioptv.eq.0.and.J.eq.1) then
125	c CALL XMIE(rayon(inq)*1.e6,REALV(K),XIMGV(K),
126	c & QEXT,QSCT,QABS,QBAR,WNOV(K))
127
128	CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),rayon(inq),
129	& QEXT,QSCT,QABS,QBAR)
130
131	c print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT
132
133	QM1(inq,K)=QEXT
134	QM2(inq,K)=QSCT
135	QM3(inq,K)=QABS
136	QM4(inq,K)=QBAR
137	endif
138
139	TAEROS=QM1(inq,K)zqaer_1pt(NLAYER+1-J,inq)1.e-4+TAEROS
140	TAEROSCAT=QM2(inq,K)zqaer_1pt(NLAYER+1-J,inq)1.e-4+TAEROSCAT
141	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(NLAYER+1-J,inq)*1.e-4
142
143	ELSE ! aerosols fractals
144
145	XMONO=(rayon(inq)/RF(inq))**3.
146	XRULE=1.
147
148	if(XMONO.gt.16384./1.5) then
149	XRULE=(XMONO/16384.)
150	XMONO=16384.
151	endif
152
153	if(ioptv.eq.0.and.J.eq.1) then
154
155	c CALL OPTFRAC(XMONO,10000./WNOV(K)
156	c & ,QEXT,QSCT,QABS,QBAR)
157
158	CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
159	& ,XMONO,QSCT,QEXT,QABS,QBAR)
160
161
162	QF1(inq,K)=QEXT*XRULE
163	QF2(inq,K)=QSCT*XRULE
164	QF3(inq,K)=QABS*XRULE
165	QF4(inq,K)=QBAR
166
167	c print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
168
169	endif
170
171	TAEROS=QF1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROS
172	TAEROSCAT=QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROSCAT
173	CBAR=CBAR+QF4(inq,K)QF2(inq,K)zqaer_1pt(NLAYER+1-J,inq)
174
175	ENDIF
176
177
178	ENDDO ! nrad
179
180
181	CBAR=CBAR/TAEROSCAT
182
183	DELTAZ=Z(J)-Z(J+1)
184
185	c --------------------------------------------------------------------
186	c profil brume Pascal: fit T (sauf tropopause) et albedo
187	c -------------------
188	if( cutoff.eq.1) then
189	IF(PRESS(J).gt.9.e-3) THEN
190	TAEROS=TAEROSM1DELTAZ/DELTAZM10.85
191	TAEROSCAT=TAEROSCATM1DELTAZ/DELTAZM10.85
192	c TAEROS=0.
193	c TAEROSCAT=0.
194	ENDIF
195
196	IF(PRESS(J).gt.1.e-1) THEN
197	TAEROS=TAEROSM1DELTAZ/DELTAZM11.15
198	TAEROSCAT=TAEROSCATM1DELTAZ/DELTAZM11.15
199	c TAEROS=0.
200	c TAEROSCAT=0.
201	ENDIF
202	endif !cutoff=1
203
204	c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
205	c -----------------------
206	if( cutoff.eq.2) then
207	IF(PRESS(J).gt.1.e-1) THEN
208	TAEROS=0.
209	TAEROSCAT=0.
210	ENDIF
211	endif !cutoff=2
212	c --------------------------------------------------------------------
213
214	TAEROSM1=TAEROS
215	TAEROSCATM1=TAEROSCAT
216	DELTAZM1=DELTAZ
217
218
219	IF (TAEROSCAT.le.0.) CBAR=0.
220
221	c if (IPRINT.eq.1) then
222	c if (k.eq.NSPECV/2) then
223	c write(*,1699) '@VI',K,J,TAEROS,TAEROSCAT,CBAR
224	c write(*,1699) '@ ',K,J,QF1(1,K),QF2(1,K),zqaer_1pt(NLAYER+1-J,1)
225	c write(*,1699) '@ ',K,J,QF1(3,K),QF2(3,K),zqaer_1pt(NLAYER+1-J,3)
226	c write(*,1699) '@ ',K,J,QF1(5,K),QF2(5,K),zqaer_1pt(NLAYER+1-J,5)
227	c write(*,1699) '@ ',K,J,QF1(7,K),QF2(7,K),zqaer_1pt(NLAYER+1-J,7)
228	c write(*,1699) '@ ',K,J,QF1(9,K),QF2(9,K),zqaer_1pt(NLAYER+1-J,9)
229	c print*
230	c endif
231	c endif
232
233	1699 FORMAT(a3,2I3,3(ES15.7,1X))
234
235	c********* EN TRAVAUX *************************
236
237	C #2: RAYLEIGH
238	c-------------------------------
239
240	C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
241	C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
242	C COLUMN NUMBER ON EARTH. CM-2
243	C THIS IS THE SCATTERING BY THE ATMOSPHERE
244
245	TAURAY=(COLDEN(J)28.9/(XMU(J)1013.25))*
246	&(.008569/WLNV(K)*4)(1.+.0113/WLNV(K)2+.00013/WLNV(K)4)
247
248	c PRINT*,WLNV(K)
249	c COLX=0.
250	c COLP=0.
251	c COLT=0.
252	c DO IU=1,NLAYER
253	c COLP=COLDEN(IU)1.e+11.35+COLP
254	c TAURAY=(COLDEN(IU)28.9/(XMU(IU)1013.25))*
255	c & (.008569/WLNV(K)*4)(1.+.0113/WLNV(K)**2
256	c & +.00013/WLNV(K)**4)
257	c COLT=COLT+TAURAY
258	c COLX=COLDEN(IU)1.e+1/(1.E528./22.4E3)1.e-10.0933e-1+COLX
259	c \|
260	c \|
261	c g/cm2->kg/m2 \| m2/kg
262	c Print*,IU, tauray,
263	c & COLDEN(IU)1.e+1/(1.E528./22.4E3)1.e-10.543e-1
264	c ENDDO
265	c PRINT*,COLP,' PRESSURE AT GROUND;'
266	c PRINT*,COLX,' TAU_GAS AT GROUND;'
267	c print*,colt,colx,' COLT, COLX'
268	c STOP
269
270	c DZ=Z(J)-Z(J+1)
271	c PRINT*, Z(J),WLNV(K),
272	c &(28.9/(XMU(J)1013.25))(.008569/WLNV(K)*4)
273	c &(1.+.0113/WLNV(K)2+.00013/WLNV(K)4)
274	c & ,COLDEN(J)/DZ/100000.,
275	c &(28.9/(XMU(J)1013.25))(.008569/WLNV(K)*4)
276	c &(1.+.0113/WLNV(K)2+.00013/WLNV(K)4)
277	c & *COLDEN(J)/DZ/100000.
278
279
280
281	C #3: CLOUD
282	c----------------------------
283
284	C NEXT COMPUTE TAU CLOUD
285
286	IF (clouds.eq.0) THEN
287	CNBAR=0.
288	TNUSCAT=0.
289	TNUABS=0.
290	TBNUABS=0.
291	ELSE
292	CNBAR=0.
293	TNUSCAT=0.
294	TNUABS=0.
295	TBNUABS=0.
296	QEXTC=0.
297	QSCTC=0.
298	QABSC=0.
299	CBARC=0.
300
301	do inq=1,nrad
302	QC1(INQ,k)=0.
303	QC2(INQ,k)=0.
304	QC3(INQ,k)=0.
305	QC4(INQ,k)=0.
306	enddo
307
308	IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
309
310	** OPTICAL CONSTANT : MIXING RULES
311
312	XNR=xfrb(nlayer+1-J,1)*REALV(K) !
313	& +xfrb(nlayer+1-J,2)*RCLDV(K) !
314	& +xfrb(nlayer+1-J,3)*RCLDV2(K) !
315	& +xfrb(nlayer+1-J,4)*RCLDV2(K) !
316
317	XNI=xfrb(nlayer+1-J,1)*XIMGV(K)
318	& +xfrb(nlayer+1-J,2)*XICLDV(K)
319	& +xfrb(nlayer+1-J,3)*XICLDV2(K)
320	& +xfrb(nlayer+1-J,4)*XICLDV2(K)
321
322	** OPTICAL CONSTANT : LIQUID DROP = THOLIN
323	IF(xfrb(nlayer+1-J,1).ge.0.01) THEN
324	XNI=XIMGV(K)
325	XNR=REALV(K)
326	ENDIF
327
328	IF (XNI.gt.1.e-10 .and. XNR.gt.1.00) THEN
329	CALL CMIE(1.E-2/WNOV(K),XNR,XNI,
330	& rcdb(nlayer+1-J),
331	& QEXTC,QSCTC,QABSC,CBARC)
332	ELSE
333	PRINT*,' WARNING XNR/XNI in optcv: ',XNR,XNI
334	QEXTC=0.
335	QSCTC=0.
336	QABSC=0.
337	CBARC=0.
338	STOP
339	ENDIF
340	ELSE
341	QEXTC=0.
342	QSCTC=0.
343	QABSC=0.
344	CBARC=0.
345	ENDIF
346
347	DO inq=1,nrad
348
349	QC1(INQ,k)=QEXTC/xnuf
350	QC2(INQ,k)=QSCTC/xnuf
351	QC3(INQ,k)=QABSC/xnuf
352	QC4(INQ,k)=CBARC
353
354	TNUABS=QC1(inq,K)zqaer_1pt(NLAYER+1-J,inq+nrad)1.e-4
355	& +TNUABS
356
357	TNUSCAT=QC2(inq,K)zqaer_1pt(NLAYER+1-J,inq+nrad)1.e-4
358	& +TNUSCAT
359
360	CNBAR=QC4(inq,K)QC2(inq,K)
361	& zqaer_1pt(NLAYER+1-J,inq+nrad)*1.e-4 + CNBAR
362
363	ENDDO
364
365	IF(TNUSCAT.EQ.0) CNBAR=0.
366	IF(TNUSCAT.NE.0.) CNBAR=CNBAR/TNUSCAT
367
368
369	ENDIF ! Cond. CLD
370
371	TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
372	TAUGVD_1pt(J,K)=TAURV_1pt(K)
373
374	TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
375	TAUHVD_1pt(J,K)=TAUHV_1pt(K)
376
377	TAUCV_1pt(K)=TAUCV_1pt(K)+TNUABS
378	TAUCVD_1pt(J,K)=TAUCV_1pt(K)
379
380
381	C #4: TAUGAS
382	C----------------------------
383
384	C LOOP OVER THE NTERMS
385	C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
386
387
388	DO 909 NT=1,NTERM(K)
389	TAUGAS=COLDEN(J)GAS1(J)BTERM(NT,K)*
390	& ( (PRESS(J+1) + PRESS(J)).5 )*PEXPON(K)
391
392
393	* COSBV ET COSBVP
394	*-----------------
395
396	IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
397	COSBV_1pt(J,K,NT)=(CBARTAEROSCAT + CNBARTNUSCAT)
398	& /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
399	ELSE
400	COSBV_1pt(J,K,NT)=0.
401	ENDIF
402
403	IF(TAEROSCAT+TAURAY .ne. 0.) THEN
404	COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
405	& /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
406	ELSE
407	COSBVP_1pt(J,K,NT)=0.
408	ENDIF
409
410	* DTAUV ET DTAUVP
411	*-----------------
412
413	DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUABS !TAU_ABS_METH
414	DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY !TAU_ABS_METH
415
416	TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
417
418	* WBARV ET WBARVP
419	*-----------------
420
421	IF(TAUGAS+TAEROS+TAURAY+TNUABS .ne. 0.) THEN
422	WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
423	& /(TAUGAS+TAEROS+TAURAY+TNUABS)
424	ELSE
425	WBARV_1pt(J,K,NT)=0.
426	ENDIF
427
428	IF(TAUGAS+TAEROS+TAURAY .ne. 0.) THEN
429	WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
430	& /(TAUGAS+TAEROS+TAURAY)
431	ELSE
432	WBARVP_1pt(J,K,NT)=0.
433	ENDIF
434
435	909 CONTINUE
436	TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
437	100 CONTINUE
438	ioptv=1
439
440	c HERE END OF THE LOOPS *******
441	c******************************
442
443	C TOTAL EXTINCTION OPTICAL DEPTHS
444	DO 119 K=1,NSPECV
445	C LOOP OVER NTERMS
446	DO 119 NT=1,NTERM(K)
447	TAUV_1pt(1,K,NT)=0.0
448	TAUVP_1pt(1,K,NT)=0.0
449	DO 119 J=1,NLAYER
450	TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
451	TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
452	119 CONTINUE
453
454	c print*,'SETUP'
455	c do i=1,NSPECV
456	c print*,WLNV(i)
457	c do j=1,NLAYER+1
458	c print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
459	c enddo
460	c enddo
461	c
462	c IF (IPRINT .GT. 1) THEN
463	c NT=1
464	c IF (2 .GT. 1) THEN
465	c WRITE (6,120)
466	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
467	c WRITE(6,*) 'latitude:',ig
468	c DO 200 K=1,NSPECV
469	c WRITE (6,190)
470	c WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
471	c & ,BWNV(K)+DWNV(K),DWNV(K)
472	c WRITE (6,230)REALV(K),XIMGV(K)
473	c DO 195 J=1,NLAYER,NLAYER
474	C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
475	c WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
476	c & ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
477	c 195 CONTINUE
478	c WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
479	c 200 CONTINUE
480	c END IF
481
482	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
483	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
484	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
485	c &' #AEROSOLS WBAR COSBAR DTAU TAU'
486	c & ,9X,'RAY GAS AEROSOL')
487	c 220 FORMAT(8(1X,F9.3))
488	c 240 FORMAT(41X,F9.3)
489
490	if (IPRINT.eq.1) stop
491
492	RETURN
493	END

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

Download in other formats:

Original Format