Context Navigation

suaer_corrk.F90 @ 848

Last change on this file since 848 was 804, checked in by lkerber, 12 years ago

A bug was found in gfluxi.F. If co2_ice was the only aerosol, the single scattering albedo was occasionally equal to 1. This led to NaNs? in gfluxi.F because we divide by a (1-W0) statement. Now there is a test in gfluxi.F which puts any W0=1 to W0=99999. This change removed the necessity in aeropacity.F90 for aerosols which are supposed to be zero to be put at 1.E-9 (they are now set to zero). In the case of zero aerosols, a dummy co2 aerosol is created and set to zero abundance everywhere. An obsolete test was removed from inifis.F. Some extraneous print statements were removed from physiq.F90.

File size: 14.9 KB

Line
1	subroutine suaer_corrk
2
3	! inputs
4	use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,iim,jjm,naerkind
5	use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
6	use datafile_mod, only: datadir
7
8	! outputs
9	use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
10	use radcommon_h, only: radiustab,nsize,tstellar
11	use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
12	use aerosol_mod
13
14	implicit none
15
16	!==================================================================
17	! Purpose
18	! -------
19	! Initialize all radiative aerosol properties
20	!
21	! Notes
22	! -----
23	! Reads the optical properties -> Mean -> Variable assignment
24	! (ASCII files) (see radcommon_h.F90)
25	! wvl(nwvl) longsun
26	! ep(nwvl) epav QVISsQREF(L_NSPECTV)
27	! omeg(nwvl) omegav omegavis(L_NSPECTV)
28	! gfactor(nwvl) gav gvis(L_NSPECTV)
29	!
30	! Authors
31	! -------
32	! Richard Fournier (1996) Francois Forget (1996)
33	! Frederic Hourdin
34	! Jean-jacques morcrette ECMWF
35	! MODIF Francois Forget (2000)
36	! MODIF Franck Montmessin (add water ice)
37	! MODIF J.-B. Madeleine 2008W27
38	! - Optical properties read in ASCII files
39	! - Add varying radius of the particules
40	! - Add water ice clouds
41	! MODIF R. Wordsworth (2009)
42	! - generalisation to arbitrary spectral bands
43	!
44	!==================================================================
45
46	#include "callkeys.h"
47	! Optical properties (read in external ASCII files)
48	INTEGER :: nwvl ! Number of wavelengths in
49	! the domain (VIS or IR)
50
51	! REAL :: solsir ! visible to infrared ratio
52	! ! (tau_0.67um/tau_9um)
53
54	REAL, DIMENSION(:),&
55	ALLOCATABLE, SAVE :: wvl ! Wavelength axis
56	REAL, DIMENSION(:),&
57	ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis
58
59	REAL, DIMENSION(:,:),&
60	ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext
61	omeg,& ! Single Scattering Albedo
62	gfactor ! Assymetry Factor
63
64	! Local variables:
65
66	INTEGER :: iaer ! Aerosol index
67	INTEGER :: idomain ! Domain index (1=VIS,2=IR)
68	INTEGER :: ilw ! longwave index
69	INTEGER :: isw ! shortwave index
70	INTEGER :: isize ! Particle size index
71	INTEGER :: jfile ! ASCII file scan index
72	INTEGER :: file_unit = 60
73	LOGICAL :: file_ok, endwhile
74	CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
75
76	! I/O of "aerave" (subroutine that spectrally averages
77	! the single scattering parameters)
78
79	REAL lamref ! reference wavelengths
80	REAL epref ! reference extinction ep
81
82	! REAL epav(L_NSPECTI) ! Average ep (= <Qext>/Qext(lamref) if epref=1)
83	! REAL omegav(L_NSPECTI) ! Average single scattering albedo
84	! REAL gav(L_NSPECTI) ! Average assymetry parameter
85
86	REAL epavVI(L_NSPECTV) ! Average ep (= <Qext>/Qext(lamref) if epref=1)
87	REAL omegavVI(L_NSPECTV) ! Average single scattering albedo
88	REAL gavVI(L_NSPECTV) ! Average assymetry parameter
89
90	REAL epavIR(L_NSPECTI) ! Average ep (= <Qext>/Qext(lamref) if epref=1)
91	REAL omegavIR(L_NSPECTI) ! Average single scattering albedo
92	REAL gavIR(L_NSPECTI) ! Average assymetry parameter
93
94	logical forgetit ! use francois' data?
95	integer iwvl
96
97	! Local saved variables:
98
99	CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id
100	!---- Please indicate the names of the optical property files below
101	! Please also choose the reference wavelengths of each aerosol
102
103	if (noaero) then
104	print*, 'naerkind= 0'
105	endif
106	do iaer=1,naerkind
107	if (iaer.eq.iaero_co2) then
108	forgetit=.true.
109	if (.not.noaero) then
110	print*, 'naerkind= co2', iaer
111	end if
112	! visible
113	if(forgetit)then
114	file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
115	else
116	file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
117	endif
118	lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx ???
119
120	! infrared
121	if(forgetit)then
122	file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
123	else
124	file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
125	endif
126	lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS ???
127	endif ! CO2 aerosols
128	! NOTE: these lamref's are currently for dust, not CO2 ice.
129	! JB thinks this shouldn't matter too much, but it needs to be
130	! fixed / decided for the final version.
131
132	if (iaer.eq.iaero_h2o) then
133	print*, 'naerkind= h2o', iaer
134
135	! visible
136	file_id(iaer,1) = 'optprop_icevis_n50.dat'
137	lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx
138	! infrared
139	file_id(iaer,2) = 'optprop_iceir_n50.dat'
140	lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS
141	endif
142
143	if (iaer.eq.iaero_dust) then
144	print*, 'naerkind= dust', iaer
145
146	! visible
147	file_id(iaer,1) = 'optprop_dustvis_n50.dat'
148	!lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx
149	lamrefvis(iaer)=0.67e-6
150	! infrared
151	file_id(iaer,2) = 'optprop_dustir_n50.dat'
152	!lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS
153	lamrefir(iaer)=9.3E-6
154	endif
155
156	if (iaer.eq.iaero_h2so4) then
157	print*, 'naerkind= h2so4', iaer
158
159	! visible
160	file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
161	!lamrefir(iaer)= doesn't exist?
162	lamrefvis(iaer)=1.5E-6 ! no idea, must find
163	! infrared
164	file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
165	!lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS
166	lamrefir(iaer)=9.3E-6 ! no idea, must find
167	! added by LK
168	endif
169
170	enddo
171
172	!------------------------------------------------------------------
173
174	! Initializations:
175
176	radiustab(:,:,:) = 0.0
177	gvis(:,:,:) = 0.0
178	omegavis(:,:,:) = 0.0
179	QVISsQREF(:,:,:) = 0.0
180	gir(:,:,:) = 0.0
181	omegair(:,:,:) = 0.0
182	QIRsQREF(:,:,:) = 0.0
183
184	DO iaer = 1, naerkind ! Loop on aerosol kind
185	DO idomain = 1, 2 ! Loop on radiation domain (VIS or IR)
186	!==================================================================
187	! 1. READ OPTICAL PROPERTIES
188	!==================================================================
189
190	! 1.1 Open the ASCII file
191
192	INQUIRE(FILE=TRIM(datadir)//&
193	'/'//TRIM(file_id(iaer,idomain)),&
194	EXIST=file_ok)
195	IF(.NOT.file_ok) THEN
196	write(,)'suaer_corrk: Problem opening ',&
197	TRIM(file_id(iaer,idomain))
198	write(,)'It should be in: ',TRIM(datadir)
199	write(,)'1) You can set this directory address ',&
200	'in callphys.def with:'
201	write(,)' datadir = /absolute/path/to/datagcm'
202	write(,)'2) If ',&
203	TRIM(file_id(iaer,idomain)),&
204	' is a LMD reference datafile, it'
205	write(,)' can be obtained online at:'
206	write(,)' http://www.lmd.jussieu.fr/',&
207	'~forget/datagcm/datafile'
208	CALL ABORT
209	ENDIF
210	OPEN(UNIT=file_unit,&
211	FILE=TRIM(datadir)//&
212	'/'//TRIM(file_id(iaer,idomain)),&
213	FORM='formatted')
214
215	! 1.2 Allocate the optical property table
216
217	jfile = 1
218	endwhile = .false.
219	DO WHILE (.NOT.endwhile)
220	READ(file_unit,*) scanline
221	IF ((scanline(1:1) .ne. '#').and.&
222	(scanline(1:1) .ne. ' ')) THEN
223	BACKSPACE(file_unit)
224	reading1_seq: SELECT CASE (jfile) ! ====================
225	CASE(1) reading1_seq ! nwvl ----------------------------
226	read(file_unit,*) nwvl
227	jfile = jfile+1
228	CASE(2) reading1_seq ! nsize ---------------------------
229	read(file_unit,*) nsize(iaer,idomain)
230	endwhile = .true.
231	CASE DEFAULT reading1_seq ! ----------------------------
232	WRITE(,) 'readoptprop: ',&
233	'Error while loading optical properties.'
234	CALL ABORT
235	END SELECT reading1_seq ! ==============================
236	ENDIF
237	ENDDO
238
239	ALLOCATE(wvl(nwvl)) ! wvl
240	ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
241	ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
242	ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
243	ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
244
245
246	! 1.3 Read the data
247
248	jfile = 1
249	endwhile = .false.
250	DO WHILE (.NOT.endwhile)
251	READ(file_unit,*) scanline
252	IF ((scanline(1:1) .ne. '#').and.&
253	(scanline(1:1) .ne. ' ')) THEN
254	BACKSPACE(file_unit)
255	reading2_seq: SELECT CASE (jfile) ! ====================
256	CASE(1) reading2_seq ! wvl -----------------------------
257	read(file_unit,*) wvl
258	jfile = jfile+1
259	CASE(2) reading2_seq ! radiusdyn -----------------------
260	read(file_unit,*) radiusdyn
261	jfile = jfile+1
262	CASE(3) reading2_seq ! ep ------------------------------
263	isize = 1
264	DO WHILE (isize .le. nsize(iaer,idomain))
265	READ(file_unit,*) scanline
266	IF ((scanline(1:1) .ne. '#').and.&
267	(scanline(1:1) .ne. ' ')) THEN
268	BACKSPACE(file_unit)
269	read(file_unit,*) ep(:,isize)
270	isize = isize + 1
271	ENDIF
272	ENDDO
273
274	jfile = jfile+1
275	CASE(4) reading2_seq ! omeg ----------------------------
276	isize = 1
277	DO WHILE (isize .le. nsize(iaer,idomain))
278	READ(file_unit,*) scanline
279	IF ((scanline(1:1) .ne. '#').and.&
280	(scanline(1:1) .ne. ' ')) THEN
281	BACKSPACE(file_unit)
282	read(file_unit,*) omeg(:,isize)
283	isize = isize + 1
284	ENDIF
285	ENDDO
286
287	jfile = jfile+1
288	CASE(5) reading2_seq ! gfactor -------------------------
289	isize = 1
290	DO WHILE (isize .le. nsize(iaer,idomain))
291	READ(file_unit,*) scanline
292	IF ((scanline(1:1) .ne. '#').and.&
293	(scanline(1:1) .ne. ' ')) THEN
294	BACKSPACE(file_unit)
295	read(file_unit,*) gfactor(:,isize)
296	isize = isize + 1
297	ENDIF
298	ENDDO
299
300	jfile = jfile+1
301	IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
302	endwhile = .true.
303	ENDIF
304	CASE(6) reading2_seq
305	endwhile = .true.
306	CASE DEFAULT reading2_seq ! ----------------------------
307	WRITE(,) 'readoptprop: ',&
308	'Error while loading optical properties.'
309	CALL ABORT
310	END SELECT reading2_seq ! ==============================
311	ENDIF
312	ENDDO
313
314	! 1.4 Close the file
315
316	CLOSE(file_unit)
317
318	! 1.5 If Francois' data, convert wvl to metres
319	if(iaer.eq.iaero_co2)then
320	if(forgetit)then
321	! print*,'please sort out forgetit for naerkind>1'
322	do iwvl=1,nwvl
323	wvl(iwvl)=wvl(iwvl)*1.e-6
324	enddo
325	endif
326	endif
327
328
329
330
331
332
333	!==================================================================
334	! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
335	!==================================================================
336	domain: SELECT CASE (idomain)
337	!==================================================================
338	CASE(1) domain ! VISIBLE DOMAIN (idomain=1)
339	!==================================================================
340
341	lamref=lamrefvis(iaer)
342	epref=1.E+0
343
344	DO isize=1,nsize(iaer,idomain)
345
346	! Save the particle sizes
347	radiustab(iaer,idomain,isize)=radiusdyn(isize)
348
349	! Averaged optical properties (GCM channels)
350
351	! CALL aerave ( nwvl,&
352	! wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
353	! lamref,epref,tstellar,&
354	! L_NSPECTV,blamv,epavVI,&
355	! omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
356	CALL aerave_new ( nwvl,&
357	wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
358	lamref,epref,tstellar,&
359	L_NSPECTV,blamv,epavVI,&
360	omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
361
362	! Variable assignments (declared in radcommon)
363	DO isw=1,L_NSPECTV
364	QVISsQREF(isw,iaer,isize)=epavVI(isw)
365	gvis(isw,iaer,isize)=gavVI(isw)
366	omegavis(isw,iaer,isize)=omegavVI(isw)
367	END DO
368
369	ENDDO
370
371	!==================================================================
372	CASE(2) domain ! INFRARED DOMAIN (idomain=2)
373	!==================================================================
374
375
376	DO isize=1,nsize(iaer,idomain) ! ----------------------------------
377
378	lamref=lamrefir(iaer)
379	epref=1.E+0
380
381	! Save the particle sizes
382	radiustab(iaer,idomain,isize)=radiusdyn(isize)
383
384	! Averaged optical properties (GCM channels)
385
386	! epav is <QIR>/Qext(lamrefir) since epref=1
387	! Note: aerave also computes the extinction coefficient at
388	! the reference wavelength. This is called QREFvis or QREFir
389	! (not epref, which is a different parameter).
390	! Reference wavelengths SHOULD BE defined for each aerosol in
391	! radcommon_h.F90
392
393	! CALL aerave ( nwvl,&
394	! wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
395	! lamref,epref,tplanet,&
396	! L_NSPECTI,blami,epavIR,&
397	! omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
398	CALL aerave_new ( nwvl,&
399	wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
400	lamref,epref,tplanet,&
401	L_NSPECTI,blami,epavIR,&
402	omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
403
404
405	! Variable assignments (declared in radcommon)
406	DO ilw=1,L_NSPECTI
407	QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
408	gir(ilw,iaer,isize)=gavIR(ilw)
409	omegair(ilw,iaer,isize)=omegavIR(ilw)
410	END DO
411
412
413	ENDDO ! isize (particle size) -------------------------------------
414
415	END SELECT domain
416
417	!========================================================================
418	! 3. Deallocate temporary variables that were read in the ASCII files
419	!========================================================================
420
421	DEALLOCATE(wvl) ! wvl
422	DEALLOCATE(radiusdyn) ! radiusdyn
423	DEALLOCATE(ep) ! ep
424	DEALLOCATE(omeg) ! omeg
425	DEALLOCATE(gfactor) ! g
426
427	END DO ! Loop on iaer
428	END DO ! Loop on idomain
429	!========================================================================
430
431	RETURN
432
433
434
435	END subroutine suaer_corrk
436

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

Download in other formats:

Original Format