Context Navigation

suaer.F90 @ 1047

Last change on this file since 1047 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 18.4 KB

Line
1	SUBROUTINE suaer
2	use dimradmars_mod, only: longrefvis, longrefir, nsizemax, long1vis, &
3	long2vis, long3vis, long1ir, long2ir, long1co2, &
4	long2co2, nsun, nir
5	use yomaer_h, only: radiustab, gvis, omegavis, QVISsQREF, gIR, omegaIR, &
6	QIRsQREF, QREFvis, QREFir, omegaREFvis, omegaREFir, &
7	nsize
8	IMPLICIT NONE
9	!==================================================================
10	! Purpose.
11	! --------
12	! initialize yomaer, the common that contains the
13	! radiative characteristics of the aerosols
14	!
15	! AUTHORS.
16	! --------
17	! Richard Fournier (1996) Francois Forget (1996)
18	! Frederic Hourdin
19	! Jean-jacques morcrette ECMWF
20	! MODIF Francois Forget (2000)
21	! MODIF Franck Montmessin (add water ice)
22	! MODIF J.-B. Madeleine 2008W27
23	! - Optical properties read in ASCII files
24	! - Add varying radius of the particules
25	!
26	! Summary.
27	! --------
28	!
29	! Read the optical properties -> Mean -> Variable assignment
30	! (ASCII files) (see yomaer.h)
31	! wvl(nwvl) longsun
32	! ep(nwvl) epav QVISsQREF(nsun)
33	! omeg(nwvl) omegav omegavis(nsun)
34	! gfactor(nwvl) gav gvis(nsun)
35	!
36	!==================================================================
37
38	! Includes:
39
40	#include "callkeys.h"
41	#include "datafile.h"
42	!#include "dimensions.h"
43	!#include "dimphys.h"
44	!#include "dimradmars.h"
45	!#include "yomaer.h"
46	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
47	#include"scatterers.h"
48	#include "aerkind.h"
49
50	! Optical properties (read in external ASCII files)
51	INTEGER :: nwvl ! Number of wavelengths in
52	! the domain (VIS or IR)
53	REAL, DIMENSION(:),&
54	ALLOCATABLE, SAVE :: wvl ! Wavelength axis
55	REAL, DIMENSION(:),&
56	ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis
57
58	REAL, DIMENSION(:,:),&
59	ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext
60	omeg,& ! Single Scattering Albedo
61	gfactor ! Assymetry Factor
62
63	! Local variables:
64
65	INTEGER :: iaer ! Aerosol index
66	INTEGER :: idomain ! Domain index (1=VIS,2=IR)
67	INTEGER :: iir ! IR channel index
68	! iir=1: 15um CO2 bands
69	! iir=2 : CO2 band wings
70	! iir=3 : 9 um band
71	! iir=4 : Far IR
72	INTEGER :: isun ! Solar band index
73	INTEGER :: isize ! Particle size index
74	INTEGER :: jfile ! ASCII file scan index
75	INTEGER :: file_unit = 60
76	LOGICAL :: file_ok, endwhile
77	CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
78	INTEGER :: read_ok
79
80	! I/O of "aerave" (subroutine averaging spectrally
81	! sing.scat.parameters)
82
83	REAL tsun ! Sun brightness temperature (for SW)
84	REAL tsol ! Surface reference brightness temp (LW)
85	REAL longref ! reference wavelengths
86	REAL longsun(nsun+1) ! solar band boundaries
87	REAL longir(nir+1) ! IR band boundaries
88	REAL epref ! reference extinction ep
89	! at wavelength "longref"
90	REAL epav(nir) ! average ep
91	! (= <Qext>/Qext(longref) if epref=1)
92	REAL omegav(nir) ! Average sing.scat.albedo
93	REAL gav(nir) ! Average assymetry parameter
94
95	! Local saved variables:
96
97	CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id
98
99	!==================================================================
100	!---- Please indicate the names of the optical property files below
101	! Please also choose the reference wavelengths of each aerosol
102
103	DO iaer = 1, naerkind ! Loop on aerosol kind
104	aerkind: SELECT CASE (name_iaer(iaer))
105	!==================================================================
106	CASE("dust_conrath") aerkind ! Typical dust profile
107	!==================================================================
108	! Visible domain:
109	file_id(iaer,1) = 'optprop_dustvis_TM.dat' !M.Wolff TM
110	! file_id(iaer,1) = 'optprop_dustvis_clancy.dat' !Clancy-Lee
111	! file_id(iaer,1) = 'optprop_dustvis_ockert.dat' !Ockert-Bell
112	! Infrared domain:
113	file_id(iaer,2) = 'optprop_dustir_TM.dat' !M.Wolff TM
114	! Toon-Forget + solsir=2 using Clancy-Lee
115	! file_id(iaer,2) = 'optprop_dustir_clancy.dat'
116	! Toon-Forget + solsir=2 using Ockert-Bell
117	! file_id(iaer,2) = 'optprop_dustir_ockert.dat'
118	! Reference wavelength in the visible:
119	longrefvis(iaer)=0.67E-6
120	! For dust: change readtesassim accordingly;
121	! Reference wavelength in the infrared:
122	longrefir(iaer)=9.3E-6
123	!==================================================================
124	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
125	!==================================================================
126	! Visible domain:
127	file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
128	! file_id(iaer,1) = 'optprop_dustvis_n50.dat' !Mie
129	! Infrared domain:
130	file_id(iaer,2) = 'optprop_dustir_n50.dat' !Mie
131	! Reference wavelength in the visible:
132	longrefvis(iaer)=0.67E-6
133	! If not equal to 0.67e-6 -> change readtesassim accordingly;
134	! Reference wavelength in the infrared:
135	longrefir(iaer)=9.3E-6
136	!==================================================================
137	CASE("dust_submicron") aerkind ! Small dust population
138	!==================================================================
139	! Visible domain:
140	file_id(iaer,1) = 'optprop_dustvis_01um_TM.dat' !M.Wolff
141	! Infrared domain:
142	file_id(iaer,2) = 'optprop_dustir_01um_TM.dat' !M.Wolff
143	! Reference wavelength in the visible:
144	longrefvis(iaer)=0.67E-6
145	! If not equal to 0.67e-6 -> change readtesassim accordingly;
146	! Reference wavelength in the infrared:
147	longrefir(iaer)=9.3E-6
148	!==================================================================
149	CASE("h2o_ice") aerkind ! Water ice crystals
150	!==================================================================
151	! Visible domain:
152	file_id(iaer,1) = 'optprop_icevis_n30.dat' !Warren
153	! file_id(iaer,1) = 'optprop_icevis.dat' !Warren
154	! Infrared domain:
155	file_id(iaer,2) = 'optprop_iceir_n30.dat' !Warren
156	! file_id(iaer,2) = 'optprop_iceir.dat' !Warren
157	! Reference wavelength in the visible:
158	longrefvis(iaer)=0.67E-6 ! 1.5um OMEGA/MEx
159	! Reference wavelength in the infrared:
160	longrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS
161	!==================================================================
162	END SELECT aerkind
163	!==================================================================
164	WRITE(,) "Scatterer: ",trim(name_iaer(iaer))
165	WRITE(,) " corresponding files: "
166	WRITE(,) "VIS: ",trim(file_id(iaer,1))
167	WRITE(,) "IR : ",trim(file_id(iaer,2))
168	!==================================================================
169	ENDDO ! iaer (loop on aerosol kind)
170
171	! Initializations:
172
173	radiustab(1:naerkind,1:2,1:nsizemax)=0
174
175	gvis(1:nsun,1:naerkind,1:nsizemax)=0
176	omegavis(1:nsun,1:naerkind,1:nsizemax)=0
177	QVISsQREF(1:nsun,1:naerkind,1:nsizemax)=0
178
179	gIR(1:nir,1:naerkind,1:nsizemax)=0
180	omegaIR(1:nir,1:naerkind,1:nsizemax)=0
181	QIRsQREF(1:nir,1:naerkind,1:nsizemax)=0
182
183	QREFvis(1:naerkind,1:nsizemax)=0
184	QREFir(1:naerkind,1:nsizemax)=0
185	omegaREFvis(1:naerkind,1:nsizemax)=0
186	omegaREFir(1:naerkind,1:nsizemax)=0
187
188	DO iaer = 1, naerkind ! Loop on aerosol kind
189	DO idomain = 1, 2 ! Loop on radiation domain (VIS or IR)
190	!==================================================================
191	! 1. READ OPTICAL PROPERTIES
192	!==================================================================
193
194	! 1.1 Open the ASCII file
195
196	INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//&
197	'/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
198	EXIST=file_ok)
199	IF(.NOT.file_ok) THEN
200	write(,)'Problem opening ',&
201	file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))
202	write(,)'It should be in: ',&
203	datafile(1:LEN_TRIM(datafile))
204	write(,)'1) You can change this directory address in '
205	write(,)' file phymars/datafile.h'
206	write(,)'2) If ',&
207	file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
208	' is a LMD reference datafile, it'
209	write(,)' can be obtained online on:'
210	write(,)' http://www.lmd.jussieu.fr/',&
211	'~forget/datagcm/datafile'
212	write(,)'3) If the name of the file is wrong, you can'
213	write(,)' change it in file phymars/suaer.F90. Just'
214	write(,)' modify the variable called file_id.'
215	CALL ABORT
216	ENDIF
217	OPEN(UNIT=file_unit,&
218	FILE=datafile(1:LEN_TRIM(datafile))//&
219	'/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
220	FORM='formatted')
221
222	! 1.2 Allocate the optical property table
223
224	jfile = 1
225	endwhile = .false.
226	DO WHILE (.NOT.endwhile)
227	READ(file_unit,*,iostat=read_ok) scanline
228	if (read_ok.ne.0) then
229	write(,)' readoptprop: Error reading file',&
230	trim(datafile(1:LEN_TRIM(datafile))//&
231	'/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))))
232	call abort
233	endif
234	IF ((scanline(1:1) .ne. '#').and.&
235	(scanline(1:1) .ne. ' ')) THEN
236	BACKSPACE(file_unit)
237	reading1_seq: SELECT CASE (jfile) ! ====================
238	CASE(1) reading1_seq ! nwvl ----------------------------
239	read(file_unit,*,iostat=read_ok) nwvl
240	if (read_ok.ne.0) then
241	write(,)' readoptprop: Error while reading line:',&
242	trim(scanline)
243	write(,)' of file',&
244	trim(datafile(1:LEN_TRIM(datafile))//&
245	'/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))))
246	call abort
247	endif
248	jfile = jfile+1
249	CASE(2) reading1_seq ! nsize ---------------------------
250	read(file_unit,*,iostat=read_ok) nsize(iaer,idomain)
251	if (read_ok.ne.0) then
252	write(,)' readoptprop: Error while reading line:',&
253	trim(scanline)
254	write(,)' of file',&
255	trim(datafile(1:LEN_TRIM(datafile))//&
256	'/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))))
257	call abort
258	endif
259	endwhile = .true.
260	CASE DEFAULT reading1_seq ! ----------------------------
261	WRITE(,) 'readoptprop: ',&
262	'Error while loading optical properties.'
263	CALL ABORT
264	END SELECT reading1_seq ! ==============================
265	ENDIF
266	ENDDO
267
268	ALLOCATE(wvl(nwvl)) ! wvl
269	ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
270	ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
271	ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
272	ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
273
274	! 1.3 Read the data
275
276	jfile = 1
277	endwhile = .false.
278	DO WHILE (.NOT.endwhile)
279	READ(file_unit,*) scanline
280	IF ((scanline(1:1) .ne. '#').and.&
281	(scanline(1:1) .ne. ' ')) THEN
282	BACKSPACE(file_unit)
283	reading2_seq: SELECT CASE (jfile) ! ====================
284	CASE(1) reading2_seq ! wvl -----------------------------
285	read(file_unit,*) wvl
286	jfile = jfile+1
287	CASE(2) reading2_seq ! radiusdyn -----------------------
288	read(file_unit,*) radiusdyn
289	jfile = jfile+1
290	CASE(3) reading2_seq ! ep ------------------------------
291	isize = 1
292	DO WHILE (isize .le. nsize(iaer,idomain))
293	READ(file_unit,*) scanline
294	IF ((scanline(1:1) .ne. '#').and.&
295	(scanline(1:1) .ne. ' ')) THEN
296	BACKSPACE(file_unit)
297	read(file_unit,*) ep(:,isize)
298	isize = isize + 1
299	ENDIF
300	ENDDO
301	jfile = jfile+1
302	CASE(4) reading2_seq ! omeg ----------------------------
303	isize = 1
304	DO WHILE (isize .le. nsize(iaer,idomain))
305	READ(file_unit,*) scanline
306	IF ((scanline(1:1) .ne. '#').and.&
307	(scanline(1:1) .ne. ' ')) THEN
308	BACKSPACE(file_unit)
309	read(file_unit,*) omeg(:,isize)
310	isize = isize + 1
311	ENDIF
312	ENDDO
313	jfile = jfile+1
314	CASE(5) reading2_seq ! gfactor -------------------------
315	isize = 1
316	DO WHILE (isize .le. nsize(iaer,idomain))
317	READ(file_unit,*) scanline
318	IF ((scanline(1:1) .ne. '#').and.&
319	(scanline(1:1) .ne. ' ')) THEN
320	BACKSPACE(file_unit)
321	read(file_unit,*) gfactor(:,isize)
322	isize = isize + 1
323	ENDIF
324	ENDDO
325	endwhile = .true.
326	CASE DEFAULT reading2_seq ! ----------------------------
327	WRITE(,) 'suaer.F90: ',&
328	'Error while loading optical properties.'
329	CALL ABORT
330	END SELECT reading2_seq ! ==============================
331	ENDIF
332	ENDDO
333
334	! 1.4 Close the file
335
336	CLOSE(file_unit)
337
338	!==================================================================
339	! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
340	!==================================================================
341	domain: SELECT CASE (idomain)
342	!==================================================================
343	CASE(1) domain ! VISIBLE DOMAIN (idomain=1)
344	!==================================================================
345
346	! 2.1 Parameters
347	tsun=6000.E+0
348	longsun(1)=long1vis
349	longsun(2)=long2vis
350	longsun(3)=long3vis
351	longref=longrefvis(iaer)
352	epref=1.E+0
353
354	DO isize=1,nsize(iaer,idomain)
355	! test that there is enough room to store the data
356	if (isize.gt.nsizemax) then
357	write(,) "suaer: Error ! nsizemax is too small!"
358	write(,) " nsizemax=",nsizemax
359	write(,) " you must increase the value of nsizemax"
360	write(,) " in dimradmars_mod !"
361	stop
362	endif
363	! ------------------------------------------------
364	! 2.2 Save the particle sizes
365	radiustab(iaer,idomain,isize)=radiusdyn(isize)
366	! 2.3 Averaged optical properties (GCM channels)
367	! Notice: Aerave also computes the extinction coefficient and
368	! single scattering albedo at reference wavelength
369	! (called QREFvis and OMEGAREFvis, same in the IR,
370	! and not epref, which is a different parameter);
371	! Reference wavelengths are defined for each aerosol in
372	! dimradmars_mod.
373
374	CALL aerave ( nwvl,&
375	wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
376	longref,epref,tsun,&
377	nsun,longsun, epav,omegav,gav,&
378	QREFvis(iaer,isize),omegaREFvis(iaer,isize) )
379	! 2.4 Variable assignements (declared by yomaer.h)
380	DO isun=1,nsun
381	QVISsQREF(isun,iaer,isize)=epav(isun)
382	gvis(isun,iaer,isize)=gav(isun)
383	omegavis(isun,iaer,isize)=omegav(isun)
384	END DO
385	! 2.5 Output display
386	! WRITE(,) 'Les donnees spectrales :'
387	! WRITE(,) 'Solaire (SW) ---->'
388	! WRITE(,) 'Aerosol number: ', iaer
389	! WRITE(,) 'Rayon aerosol: ', radiustab(iaer,idomain,isize)
390	! WRITE(,) '<Qext>/Qext(longrefvis) ; omega ; g'
391	! DO isun=1,nsun
392	! WRITE(,) QVISsQREF(isun,iaer,isize),&
393	! omegavis(isun,iaer,isize),&
394	! gvis(isun,iaer,isize)
395	! ENDDO
396	! WRITE(,) 'QREFvis(',iaer,isize,') = ',QREFvis(iaer,isize)
397	! WRITE(,) 'omegaREFvis(',iaer,isize,') = ',&
398	! omegaREFvis(iaer,isize)
399	! ------------------------------------------------
400	ENDDO
401
402	!==================================================================
403	CASE(2) domain ! INFRARED DOMAIN (idomain=2)
404	!==================================================================
405
406	DO isize=1,nsize(iaer,idomain) ! ----------------------------------
407
408	! 2.1 solsir is not used anymore; division of Qext(IR) by solsir
409	! has to be done in the input ASCII files (if necessary).
410
411	! 2.2 Save the particle sizes
412	radiustab(iaer,idomain,isize)=radiusdyn(isize)
413
414	! 2.3 Parameters
415
416	tsol=215.D+0
417	longir(1)=long1ir
418	longir(2)=long1co2
419	longir(3)=long2co2
420	longir(4)=long2ir
421	longref=longrefir(iaer)
422	epref=1.E+0
423
424	! 2.4 Averaged optical properties (GCM channels)
425	! epav is <QIR>/Qext(longrefir) since epref=1
426	! Notice: Aerave also Computes the extinction coefficient at
427	! reference wavelength (called QREFvis or QREFir,
428	! and not epref, which is a different parameter);
429	! Reference wavelengths are defined for each aerosol in
430	! dimradmar_mod.
431
432	CALL aerave ( nwvl,&
433	wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
434	longref,epref,tsol,&
435	nir-1,longir,epav,omegav,gav,&
436	QREFir(iaer,isize),omegaREFir(iaer,isize) )
437	! WRITE(,) 'QREFir(',iaer,isize,') = ',QREFir(iaer,isize)
438	! WRITE(,) 'omegaREFir(',iaer,isize,') = ',&
439	! omegaREFir(iaer,isize)
440
441	! 2.5 Computing <QIR>/Qext(longrefvis)
442
443	DO iir=1,nir-1
444	! WRITE(,) 'QIRsQREFir Channel ',iir,': ',epav(iir)
445	epav(iir)= epav(iir) * QREFir(iaer,isize) / &
446	QREFvis(iaer,isize)
447	ENDDO
448	! WRITE(,) 'Aerosol number', iaer
449	! WRITE(,) 'Particle size: ',radiustab(iaer,idomain,isize)
450	! WRITE(,) 'Rapport Solaire/IR:',&
451	! QREFvis(iaer,isize) / QREFir(iaer,isize)
452
453	! 2.6 Variable assignements
454	! (variables are declared by yomaer.h)
455
456	! Single scattering properties
457	! in each of the "nir" bands
458	! (cf. dimradmars_mod)
459
460	! iir=1 : central 15um CO2 bands
461	QIRsQREF(1,iaer,isize)=epav(2)
462	omegaIR(1,iaer,isize)=omegav(2)
463	gIR(1,iaer,isize)=gav(2)
464
465	! iir=2 : CO2 band wings
466	! (same properties than for central part)
467	QIRsQREF(2,iaer,isize)=epav(2)
468	omegaIR(2,iaer,isize)=omegav(2)
469	gIR(2,iaer,isize)=gav(2)
470
471	! iir=3 : 9 um band [long1ir - long1co2]
472	QIRsQREF(3,iaer,isize)=epav(1)
473	omegaIR(3,iaer,isize)=omegav(1)
474	gIR(3,iaer,isize)=gav(1)
475
476	! iir=4 : Far IR [long2co2 - long2ir]
477	QIRsQREF(4,iaer,isize)=epav(3)
478	omegaIR(4,iaer,isize)=omegav(3)
479	gIR(4,iaer,isize)=gav(3)
480
481	! 2.7 Output display
482
483	! WRITE(,) 'AEROSOL PROPERTIES: Number ',iaer
484	! WRITE(,) 'Thermal IR (LW) ---->'
485	! WRITE(,) 'Particle size: ',radiustab(iaer,idomain,isize)
486	! WRITE(,) '<Qext>/Qext(longrefvis) ; omega ; g'
487	! DO iir=1,nir
488	! WRITE(,) QIRsQREF(iir,iaer,isize),omegaIR(iir,iaer,isize),&
489	! gIR(iir,iaer,isize)
490	! ENDDO
491	! WRITE(,) 'CO2: <Qabs>/Qext(longrefvis) = ',&
492	! QIRsQREF(1,iaer,isize)*(1-omegaIR(1,iaer,isize))
493	! WRITE(,) ''
494
495	ENDDO ! isize (particle size) -------------------------------------
496
497	END SELECT domain
498	!==================================================================
499	! 3. Deallocate temporary variables (read in ASCII files)
500	!==================================================================
501
502	DEALLOCATE(wvl) ! wvl
503	DEALLOCATE(radiusdyn) ! radiusdyn
504	DEALLOCATE(ep) ! ep
505	DEALLOCATE(omeg) ! omeg
506	DEALLOCATE(gfactor) ! g
507
508	END DO ! Loop on iaer
509	END DO ! Loop on idomain
510	!==================================================================
511	RETURN
512	END

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

Download in other formats:

Original Format