Context Navigation

suaer_corrk.F90 @ 1529

Last change on this file since 1529 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

File size: 16.2 KB

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

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

Download in other formats:

Original Format