Context Navigation

optci.F @ 1525

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

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

Download in other formats:

Original Format