Context Navigation

optci.F @ 1243

Last change on this file since 1243 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 10.4 KB

Line
1	SUBROUTINE OPTCI(ykim,qaer,nmicro,IPRINT)
2	use dimphy
3	use infotrac
4	use common_mod, only:rmcbar,xfbar,ncount,TauHID,TauCID,TauGID
5	#include "dimensions.h"
6	#include "microtab.h"
7	#include "numchimrad.h"
8	#include "clesphys.h"
9
10	c Arguments:
11	c ---------
12	REAL ykim(klon,klev,nqtot)
13	real qaer(klon,klev,nqtot)
14	integer nmicro,IPRINT
15	c ---------
16
17	c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
18	INTEGER ngrid
19	PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon
20	c
21	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
22	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
23
24	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
25
26	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
27	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
28
29	COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
30	COMMON /STRAT2/ HCN(NLAYER)
31
32	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
33	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
34
35	COMMON /CLOUD/
36	& RCLDI(NSPECI), XICLDI(NSPECI)
37	& , RCLDV(NSPECV), XICLDV(NSPECV)
38	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
39	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
40
41	COMMON /TAUS/ TAUHI(ngrid,NSPECI),TAUCI(ngrid,NSPECI),
42	& TAUGI(ngrid,NSPECI),TAURV(ngrid,NSPECV),
43	& TAUHV(ngrid,NSPECV),TAUCV(ngrid,NSPECV),
44	& TAUGV(ngrid,NSPECV)
45
46	COMMON /OPTICI/ DTAUI(ngrid,NLAYER,NSPECI)
47	& ,TAUI (ngrid,NLEVEL,NSPECI)
48	& ,WBARI(ngrid,NLAYER,NSPECI)
49	& ,COSBI(ngrid,NLAYER,NSPECI)
50	& ,DTAUIP(ngrid,NLAYER,NSPECI)
51	& ,TAUIP(ngrid,NLEVEL,NSPECI)
52	& ,WBARIP(ngrid,NLAYER,NSPECI)
53	& ,COSBIP(ngrid,NLAYER,NSPECI)
54
55	COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
56	& DWNI(NSPECI), WLNI(NSPECI)
57
58	REAL DTAUP(ngrid,NLAYER,NSPECI),DTAUPP(ngrid,NLAYER,NSPECI)
59	COMMON /IRTAUS/ DTAUP,DTAUPP
60
61	COMMON /PLANT/ CSUBP,F0PI
62	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
63	COMMON /CONST/RGAS,RHOP,PI,SIGMA
64	COMMON /part/v,rayon,vrat,dr,dv
65
66	DIMENSION PROD(NLEVEL)
67	* nrad dans microtab.h
68	real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
69	real xv1(klev,nspeci),xv2(klev,nspeci)
70	real xv3(klev,nspeci)
71	REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI)
72	REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI)
73	REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
74	REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
75	real emu
76	REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
77
78	save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
79
80	integer iopti,iwarning ! iopti: premier appel, une seule boucle sur les l.d'o.
81	integer ig,seulmtunpt,iout
82	save iopti,iwarning,seulmtunpt
83	data iopti,iwarning,seulmtunpt/0,0,0/
84
85	real zqaer_1pt(NLAYER,2*nrad)
86	#include "optci_1pt.h"
87
88	character*100 dummy
89	real dummy2,dummy3
90
91	C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
92	C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
93	C
94	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
95	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
96	C LAYER: WBAR, DTAU, COSBAR
97	C LEVEL: TAU
98	C
99	print*,'START OPTCI'
100
101	c Diagnostic eventuellement:
102	c if (nmicro.gt.0) then
103	c sum=0.
104	c do nng=2,klon
105	c do i=1,klev
106	c do j=1,nmicro
107	c print*,'j,rj',j,rayon(j)
108	c print*,'paer',qaer(nng,i,j)
109	c sum=sum+qaer(nng,i,j)rayon(j)3.1.33333.14151000.
110	c enddo
111	c enddo
112	c enddo
113	c print*,sum/(klon-1),'SOMME COLONNE/OPTCI'
114	c endif
115
116
117	c do inq=1,nrad
118	c print*,inq,rayon(inq),vrat,qaer(12,25,inq)
119	c enddo
120
121	C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122	c INITIALISATIONS UNE SEULE FOIS
123	C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124
125	if (iopti.eq.0) then
126
127	c verif pour taille zqaer_1pt, sachant que si microfi=0 et nqtot=1,
128	c il faut quand meme qu on lise la look-up table de dim nrad=10
129	c et si microfi=1, on doit avoir nmicro=nrad (dans microtab.h)
130	c
131	c Nouvelle verif pour nuages :
132	c La condition ci-dessus n'est plus realisable !
133	c nmicro comprend maintenant aussi des glaces
134	c Donc on teste juste que nmicro soit > 2*nrad (ou nrad si on ne fait pas de nuages)
135	if (microfi.ge.1) then
136	if ((clouds.eq.1).and.(nmicro.lt.2*nrad)) then
137	print*,"OPTCI :"
138	print,"clouds = 1 MAIS nmicro < 2nrad"
139	print*,"Probleme pour zqaer_1pt dans optci."
140	stop
141	endif
142	if ((clouds.eq.0).and.(nmicro.lt.nrad)) then
143	print*,"OPTCI :"
144	print*,"nmicro < nrad"
145	print*,"Probleme pour zqaer_1pt dans optci."
146	stop
147	endif
148	endif
149
150	DO 420 K=1,NSPECI
151	C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE.
152	c CALL THOLIN(WLNI(K),TNR,TNI)
153	CALL THOLIN_CVD(WLNI(K),TNR,TNI)
154	REALI(K)=TNR
155	XIMGI(K)=TNI*FHIR
156	C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD
157	CALL LIQCH4(WLNI(K),TNR,TNI)
158	RCLDI(K)=TNR
159	XICLDI(K)=TNI
160	CALL LIQC2H6(WLNI(K),TNR,TNI)
161	RCLDI2(K)=TNR
162	XICLDI2(K)=TNI
163	420 CONTINUE
164
165	c DEBUG
166	c print*,"wnoi=",WNOI
167
168	C
169	C ZERO ALL OPTICAL DEPTHS.
170	C ??FLAG? FOR SOME APPLCIATIONS THE TOP OPACITY MAY NOT VANISH
171
172	c open (unit=1,file='xsetupi')
173	c do i=1,klev
174	c read(1,*) a
175	c do j=1,nspeci
176	c read(1,*) xv1(i,j),xv2(i,j),xv3(i,j)
177	c enddo
178	c enddo
179	c close(1)
180
181	endif ! fin initialisations premier appel
182
183	c************************************************************************
184	c************************************************************************
185	DO 79 ig=1,klon ! BOUCLE SUR GRILLE HORIZONTALE
186	c print*,'ig NEW optci',ig
187	c************************************************************************
188	c************************************************************************
189
190	if (.not.ylellouch) then
191
192	XN2(1) = ykim(ig,1,iradn2)
193	CH4(1) = ykim(ig,1,iradch4)
194	H2(1) = ykim(ig,1,iradh2)
195	do j=2,nlayer
196	XN2(j) = (ykim(ig,j,iradn2)+ykim(ig,j-1,iradn2))/2.
197	CH4(j) = (ykim(ig,j,iradch4)+ykim(ig,j-1,iradch4))/2.
198	H2(j) = (ykim(ig,j,iradh2)+ykim(ig,j-1,iradh2))/2.
199	enddo
200	XN2(nlevel) = ykim(ig,nlayer,iradn2)
201	CH4(nlevel) = ykim(ig,nlayer,iradch4)
202	H2(nlevel) = ykim(ig,nlayer,iradh2)
203
204	do j=1,nlayer
205	emu = ( xmu(j) + xmu(j+1) )/2.
206	C2H2(j) = ykim(ig,j,iradc2h2) * 26./emu
207	C2H6(j) = ykim(ig,j,iradc2h6) * 30./emu
208	HCN(j) = ykim(ig,j,iradhcn ) * 27./emu
209	enddo
210
211	endif
212
213	c if ((.not.ylellouch).and.(ig.eq.klon/2)) then
214	c print*,' LAYER C2H2 C2H6 HCN masmix ratios'
215	c do j=1,nlayer
216	c print*,j,C2H2(j),C2H6(j),HCN(j)
217	c enddo
218	c endif
219
220	if (microfi.ge.1) then
221	do iq=1,2*nrad
222	c si on ne fait pas de nuages on ne copie que les nrad premieres valeurs.
223	if (clouds.eq.0.and.iq.gt.nrad) then
224	zqaer_1pt(:,iq)=0.
225	else
226	do j=1,NLAYER
227	zqaer_1pt(j,iq)=qaer(ig,j,iq)
228	enddo
229	endif
230	enddo
231	else
232	if (ig.eq.1) then
233	zqaer_1pt = 0.
234	c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
235	c boucle sur nrad=10 (dans microtab.h)
236	open(10,file="qaer_eq_1d.dat")
237	do iq=1,15
238	read(10,'(A100)') dummy
239	enddo
240	do j=NLAYER,1,-1
241	read(10,*) dummy2,dummy3,(zqaer_1pt(j,iq),iq=1,nrad)
242	enddo
243	close(10)
244	c ici, les tableaux definissant la structure des aerosols sont
245	c remplis: rf,df(nq),rayon(nq,)v(nq)......
246	call rdf()
247	endif
248	endif
249
250	c if ((ig.eq.klon/2).or.(microfi.eq.0)) then
251	c print*,"Q01=",zqaer_1pt(:,1)
252	c print*,"Q05=",zqaer_1pt(:,5)
253	c print*,"Q10=",zqaer_1pt(:,10)
254	c stop
255	c endif
256
257	iout=0
258	c if ((microfi.eq.0).or.(ig.eq.(klon/2+16))) iout=1
259	if (seulmtunpt.eq.0) then
260	call optci_1pt3(zqaer_1pt,rmcbar(ig,:),xfbar(ig,:,:),
261	& iopti,iout)
262	iopti = 1
263	endif
264
265	c Pas de microphysique, ni de composition variable: un seul passage
266	c dans optci_1pt.
267	if ((microfi.eq.0).and.(ylellouch)) then
268	seulmtunpt = 1
269	endif
270
271	COSBI(ig,:,:) = MAX(MIN(COSBI_1pt(:,:),0.999999),1e-6)
272	WBARI(ig,:,:) = MAX(MIN(WBARI_1pt(:,:),0.999999),1e-6)
273	DTAUI(ig,:,:) = DTAUI_1pt(:,:)
274	TAUI(ig,:,:) = TAUI_1pt(:,:)
275
276	COSBIP(ig,:,:) = MAX(MIN(COSBIP_1pt(:,:),0.999999),1e-6)
277	WBARIP(ig,:,:) = MAX(MIN(WBARIP_1pt(:,:),0.999999),1e-6)
278	DTAUIP(ig,:,:) = DTAUIP_1pt(:,:)
279	TAUIP(ig,:,:) = TAUIP_1pt(:,:)
280
281	TAUHI(ig,:) = TAUHI_1pt(:)
282	TAUCI(ig,:) = TAUCI_1pt(:)
283	TAUGI(ig,:) = TAUGI_1pt(:)
284
285	TauHID(ig,:,:) = TAUHID_1pt(:,:)
286	TauCID(ig,:,:) = TAUCID_1pt(:,:)
287	TauGID(ig,:,:) = TAUGID_1pt(:,:)
288
289	c DEBUG
290	c if(ig.eq.(ngrid/2+16)) then
291	c print*,ig,'/',KLON,':'
292	c print*,'TauHID_1',TAUHID(ig,1,:)
293	c print*,'TauGID_1',TAUGID(ig,1,:)
294	c print*,'TauHID_50',TAUHID(ig,50,:)
295	c print*,'TauGID_50',TAUGID(ig,50,:)
296	c print*,"DTAUI_1=",DTAUI(ig,1,:)
297	c print*,"DTAUI_50=",DTAUI(ig,50,:)
298	c print*,'cosBI_1',COSBI(ig,1,:)
299	c print*,'cosBI_50',COSBI(ig,50,:)
300	c print*,'WBARI_1',WBARI(ig,1,:)
301	c print*,'WBARI_50',WBARI(ig,50,:)
302	c stop
303	c endif
304
305	c************************************************************************
306	c************************************************************************
307	79 CONTINUE ! FIN BOUCLE GRILLE HORIZONTALE
308	c************************************************************************
309	c************************************************************************
310	C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED
311	DO 225 IG=1,klon
312	DO 220 J=1,NLAYER
313	DO 230 K=1,NSPECI
314	DTAUP(IG,J,K)=DTAUI(IG,J,K)
315	DTAUPP(IG,J,K)=DTAUIP(IG,J,K)
316	230 CONTINUE
317	220 CONTINUE
318	225 CONTINUE
319
320	print*, 'FIN OPTCI'
321
322	RETURN
323	END
324

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

Download in other formats:

Original Format