Context Navigation

optcv_1pt.F @ 3303

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

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

Download in other formats:

Original Format