Context Navigation

suaer_corrk.F90 @ 2613

Last change on this file since 2613 was 2297, checked in by jvatant, 5 years ago

Add a generic n-layer aerosol scheme to replace the former buggy 2-layer scheme as well as the hard-coded NH3 cloud.

It can be called using 'aeronlay=.true.' in callphys.def, and set the number of layers (up to 4) with 'nlayaero'.
Then, the following parameters are read as arrays of size nlayaero in callphys.def (separated by blank space)

*aeronlay_tauref (Optical depth of aerosol layer at ref wavelenght)
*aeronlay_lamref (Ref wavelenght (m))
*aeronlay_choice (Choice for vertical profile - 1:tau follows atm scale height btwn top and bottom - 2:tau follows it own scale height)
*aeronlay_pbot (Bottom pressure (Pa))
*aeronlay_ptop (Top pressure (Pa) - useful only if choice=1)
*aeronlay_sclhght (Ratio of aerosol layer scale height / atmospheric scale height - useful only if choice=2 )
*aeronlay_size (Particle size (m))
*optprop_aeronlay_vis File for VIS opt properties.
*optprop_aeronlay_ir File for IR opt properties.

+Extra info :

+ In addition of solving the bug from 2-layer it enables different optical properties.
+ The former scheme are left for retrocompatibility (for now) but you should use the new one.
+ See aeropacity.F90 for the calculations

+ Each layer can have different optical properties, size of particle ...
+ You have different choices for vertical profile of the aerosol layers :

aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).

In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.

+ The reference wavelenght for input optical depth is now read as input (aeronlay_lamref)
+ Following the last point some comment is added in suaer_corrk about the 'not-really-dummy'ness of IR lamref..

File size: 18.1 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,omegarefir !,omegarefvis
12	use aerosol_mod
13	use callkeys_mod, only: tplanet, optprop_back2lay_vis, optprop_back2lay_ir, &
14	optprop_aeronlay_vis, optprop_aeronlay_ir, &
15	aeronlay_lamref, nlayaero
16
17	implicit none
18
19	!==================================================================
20	! Purpose
21	! -------
22	! Initialize all radiative aerosol properties
23	!
24	! Notes
25	! -----
26	! Reads the optical properties -> Mean -> Variable assignment
27	! (ASCII files) (see radcommon_h.F90)
28	! wvl(nwvl) longsun
29	! ep(nwvl) epav QVISsQREF(L_NSPECTV)
30	! omeg(nwvl) omegav omegavis(L_NSPECTV)
31	! gfactor(nwvl) gav gvis(L_NSPECTV)
32	!
33	! Authors
34	! -------
35	! Richard Fournier (1996) Francois Forget (1996)
36	! Frederic Hourdin
37	! Jean-jacques morcrette ECMWF
38	! MODIF Francois Forget (2000)
39	! MODIF Franck Montmessin (add water ice)
40	! MODIF J.-B. Madeleine 2008W27
41	! - Optical properties read in ASCII files
42	! - Add varying radius of the particules
43	! - Add water ice clouds
44	! MODIF R. Wordsworth (2009)
45	! - generalisation to arbitrary spectral bands
46	!
47	!==================================================================
48
49	! Optical properties (read in external ASCII files)
50	INTEGER,SAVE :: nwvl ! Number of wavelengths in
51	! the domain (VIS or IR), read by master
52
53	! REAL :: solsir ! visible to infrared ratio
54	! ! (tau_0.67um/tau_9um)
55
56	REAL, DIMENSION(:),&
57	ALLOCATABLE, SAVE :: wvl ! Wavelength axis, read by master
58	REAL, DIMENSION(:),&
59	ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis, read by master
60
61	REAL, DIMENSION(:,:),&
62	ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext, read by master
63	omeg,& ! Single Scattering Albedo, read by master
64	gfactor ! Assymetry Factor, read by master
65
66	! Local variables:
67
68	INTEGER :: iaer ! Aerosol index
69	INTEGER :: idomain ! Domain index (1=VIS,2=IR)
70	INTEGER :: ilw ! longwave index
71	INTEGER :: isw ! shortwave index
72	INTEGER :: isize ! Particle size index
73	INTEGER :: jfile ! ASCII file scan index
74	INTEGER :: file_unit = 60
75	LOGICAL :: file_ok, endwhile
76	CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
77
78	! I/O of "aerave" (subroutine that spectrally averages
79	! the single scattering parameters)
80
81	REAL lamref ! reference wavelengths
82	REAL epref ! reference extinction ep
83
84	REAL epavVI(L_NSPECTV) ! Average ep (= <Qext>/Qext(lamrefvis) if epref=1)
85	REAL omegavVI(L_NSPECTV) ! Average single scattering albedo
86	REAL gavVI(L_NSPECTV) ! Average assymetry parameter
87
88	REAL epavIR(L_NSPECTI) ! Average ep (= <Qext>/Qext(lamrefir) if epref=1)
89	REAL omegavIR(L_NSPECTI) ! Average single scattering albedo
90	REAL gavIR(L_NSPECTI) ! Average assymetry parameter
91
92	logical forgetit ! use francois' data?
93	integer iwvl, ia
94
95	! Local saved variables:
96
97	CHARACTER(LEN=50), DIMENSION(naerkind,2), SAVE :: file_id
98	!$OMP THREADPRIVATE(file_id)
99	!---- Please indicate the names of the optical property files below
100	! Please also choose the reference wavelengths of each aerosol
101
102	!--------- README TO UNDERSTAND WHAT FOLLOWS (JVO 20) -------
103	! The lamref variables comes from the Martian model
104	! where the visible one is the one used for computing
105	! and the IR one is used to output scaled opacity to
106	! match instrumental data ... This is done (at least
107	! for now) in the generic, so lamrefir is dummy*!
108
109	! The important one is the VISIBLE one as it will be used
110	! to rescale the values in callcork.F90 assuming 'aerosol' is
111	! at this visible reference wavelenght.
112
113	! *Actually if you change lamrefir here there is a
114	! slight sensitvity in the outputs because of some
115	! machine precision differences that amplifys and lead
116	! up to 10-6 differences in the radiative balance...
117	! This could be good to clean the code but would require
118	! a lot of modifs and to take care !
119
120	!--------------------------------------------------------------
121
122	if (noaero) then
123	print*, 'naerkind= 0'
124	endif
125	do iaer=1,naerkind
126	if (iaer.eq.iaero_co2) then
127	forgetit=.true.
128	if (.not.noaero) then
129	print*, 'naerkind= co2', iaer
130	end if
131	! visible
132	if(forgetit)then
133	file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
134	else
135	file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
136	endif
137	lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx ???
138
139	! infrared
140	if(forgetit)then
141	file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
142	else
143	file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
144	endif
145	lamrefir(iaer)=12.1E-6 ! Dummy in generic phys. (JVO 20)
146	endif ! CO2 aerosols
147	! NOTE: these lamref's are currently for dust, not CO2 ice.
148	! JB thinks this shouldn't matter too much, but it needs to be
149	! fixed / decided for the final version.
150
151	if (iaer.eq.iaero_h2o) then
152	print*, 'naerkind= h2o', iaer
153
154	! visible
155	file_id(iaer,1) = 'optprop_icevis_n50.dat'
156	lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx
157	! infrared
158	file_id(iaer,2) = 'optprop_iceir_n50.dat'
159	lamrefir(iaer)=12.1E-6 ! Dummy in generic phys. (JVO 20)
160	endif
161
162	if (iaer.eq.iaero_dust) then
163	print*, 'naerkind= dust', iaer
164
165	! visible
166	file_id(iaer,1) = 'optprop_dustvis_n50.dat'
167	!lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx
168	lamrefvis(iaer)=0.67e-6
169	! infrared
170	file_id(iaer,2) = 'optprop_dustir_n50.dat'
171	lamrefir(iaer)=9.3E-6 ! Dummy in generic phys. (JVO 20)
172	endif
173
174	if (iaer.eq.iaero_h2so4) then
175	print*, 'naerkind= h2so4', iaer
176
177	! visible
178	file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
179	lamrefvis(iaer)=1.5E-6 ! no idea, must find
180	! infrared
181	file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
182	lamrefir(iaer)=9.3E-6 ! ! Dummy in generic phys. (JVO 20)
183	! added by LK
184	endif
185
186	if (iaer.eq.iaero_back2lay) then
187	print*, 'naerkind= back2lay', iaer
188
189	! visible
190	file_id(iaer,1) = TRIM(optprop_back2lay_vis)
191	lamrefvis(iaer)=0.8E-6 ! This is the important one.
192	! infrared
193	file_id(iaer,2) = TRIM(optprop_back2lay_ir)
194	lamrefir(iaer)=6.E-6 ! This is dummy.
195	! added by SG
196	endif
197
198	if (iaer.eq.iaero_nh3) then
199	print*, 'naerkind= nh3', iaer
200
201	! visible
202	file_id(iaer,1) = 'optprop_nh3ice_vis.dat'
203	lamrefvis(iaer)=0.756E-6 !
204	! infrared
205	file_id(iaer,2) = 'optprop_nh3ice_ir.dat'
206	lamrefir(iaer)=6.E-6 ! dummy
207	! added by SG
208	endif
209
210	do ia=1,nlayaero
211	if (iaer.eq.iaero_nlay(ia)) then
212	print*, 'naerkind= nlay', iaer
213
214	! visible
215	file_id(iaer,1) = TRIM(optprop_aeronlay_vis(ia))
216	lamrefvis(iaer)=aeronlay_lamref(ia)
217	! infrared
218	file_id(iaer,2) = TRIM(optprop_aeronlay_ir(ia))
219	lamrefir(iaer)=6.E-6 ! Dummy value
220	endif
221	enddo
222	! added by JVO
223
224	if (iaer.eq.iaero_aurora) then
225	print*, 'naerkind= aurora', iaer
226
227	! visible
228	file_id(iaer,1) = 'optprop_aurora_vis.dat'
229	lamrefvis(iaer)=0.3E-6 !
230	! infrared
231	file_id(iaer,2) = 'optprop_aurora_ir.dat'
232	lamrefir(iaer)=6.E-6 ! dummy
233	! added by SG
234	endif
235
236
237	enddo
238
239	!------------------------------------------------------------------
240
241	! Initializations:
242
243	radiustab(:,:,:) = 0.0
244	gvis(:,:,:) = 0.0
245	omegavis(:,:,:) = 0.0
246	QVISsQREF(:,:,:) = 0.0
247	gir(:,:,:) = 0.0
248	omegair(:,:,:) = 0.0
249	QIRsQREF(:,:,:) = 0.0
250
251	DO iaer = 1, naerkind ! Loop on aerosol kind
252	DO idomain = 1, 2 ! Loop on radiation domain (VIS or IR)
253	!==================================================================
254	! 1. READ OPTICAL PROPERTIES
255	!==================================================================
256
257	! 1.1 Open the ASCII file
258
259	!$OMP MASTER
260
261	INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
262	'/'//TRIM(file_id(iaer,idomain)),&
263	EXIST=file_ok)
264	IF (file_ok) THEN
265	OPEN(UNIT=file_unit,&
266	FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
267	'/'//TRIM(file_id(iaer,idomain)),&
268	FORM='formatted')
269	ELSE
270	! In ye old days these files were stored in datadir;
271	! let's be retro-compatible
272	INQUIRE(FILE=TRIM(datadir)//&
273	'/'//TRIM(file_id(iaer,idomain)),&
274	EXIST=file_ok)
275	IF (file_ok) THEN
276	OPEN(UNIT=file_unit,&
277	FILE=TRIM(datadir)//&
278	'/'//TRIM(file_id(iaer,idomain)),&
279	FORM='formatted')
280	ENDIF
281	ENDIF
282	IF(.NOT.file_ok) THEN
283	write(,)'suaer_corrk: Problem opening ',&
284	TRIM(file_id(iaer,idomain))
285	write(,)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir)
286	write(,)'1) You can set this directory address ',&
287	'in callphys.def with:'
288	write(,)' datadir = /absolute/path/to/datagcm'
289	write(,)'2) If ',&
290	TRIM(file_id(iaer,idomain)),&
291	' is a LMD reference datafile, it'
292	write(,)' can be obtained online at:'
293	write(,)' http://www.lmd.jussieu.fr/',&
294	'~lmdz/planets/LMDZ.GENERIC/datagcm/'
295	CALL ABORT
296	ENDIF
297
298	! 1.2 Allocate the optical property table
299
300	jfile = 1
301	endwhile = .false.
302	DO WHILE (.NOT.endwhile)
303	READ(file_unit,*) scanline
304	IF ((scanline(1:1) .ne. '#').and.&
305	(scanline(1:1) .ne. ' ')) THEN
306	BACKSPACE(file_unit)
307	reading1_seq: SELECT CASE (jfile) ! ====================
308	CASE(1) reading1_seq ! nwvl ----------------------------
309	read(file_unit,*) nwvl
310	jfile = jfile+1
311	CASE(2) reading1_seq ! nsize ---------------------------
312	read(file_unit,*) nsize(iaer,idomain)
313	endwhile = .true.
314	CASE DEFAULT reading1_seq ! ----------------------------
315	WRITE(,) 'readoptprop: ',&
316	'Error while loading optical properties.'
317	CALL ABORT
318	END SELECT reading1_seq ! ==============================
319	ENDIF
320	ENDDO
321
322	ALLOCATE(wvl(nwvl)) ! wvl
323	ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
324	ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
325	ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
326	ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
327
328
329	! 1.3 Read the data
330
331	jfile = 1
332	endwhile = .false.
333	DO WHILE (.NOT.endwhile)
334	READ(file_unit,*) scanline
335	IF ((scanline(1:1) .ne. '#').and.&
336	(scanline(1:1) .ne. ' ')) THEN
337	BACKSPACE(file_unit)
338	reading2_seq: SELECT CASE (jfile) ! ====================
339	CASE(1) reading2_seq ! wvl -----------------------------
340	read(file_unit,*) wvl
341	jfile = jfile+1
342	CASE(2) reading2_seq ! radiusdyn -----------------------
343	read(file_unit,*) radiusdyn
344	jfile = jfile+1
345	CASE(3) reading2_seq ! ep ------------------------------
346	isize = 1
347	DO WHILE (isize .le. nsize(iaer,idomain))
348	READ(file_unit,*) scanline
349	IF ((scanline(1:1) .ne. '#').and.&
350	(scanline(1:1) .ne. ' ')) THEN
351	BACKSPACE(file_unit)
352	read(file_unit,*) ep(:,isize)
353	isize = isize + 1
354	ENDIF
355	ENDDO
356
357	jfile = jfile+1
358	CASE(4) reading2_seq ! omeg ----------------------------
359	isize = 1
360	DO WHILE (isize .le. nsize(iaer,idomain))
361	READ(file_unit,*) scanline
362	IF ((scanline(1:1) .ne. '#').and.&
363	(scanline(1:1) .ne. ' ')) THEN
364	BACKSPACE(file_unit)
365	read(file_unit,*) omeg(:,isize)
366	isize = isize + 1
367	ENDIF
368	ENDDO
369
370	jfile = jfile+1
371	CASE(5) reading2_seq ! gfactor -------------------------
372	isize = 1
373	DO WHILE (isize .le. nsize(iaer,idomain))
374	READ(file_unit,*) scanline
375	IF ((scanline(1:1) .ne. '#').and.&
376	(scanline(1:1) .ne. ' ')) THEN
377	BACKSPACE(file_unit)
378	read(file_unit,*) gfactor(:,isize)
379	isize = isize + 1
380	ENDIF
381	ENDDO
382
383	jfile = jfile+1
384	IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
385	endwhile = .true.
386	ENDIF
387	CASE(6) reading2_seq
388	endwhile = .true.
389	CASE DEFAULT reading2_seq ! ----------------------------
390	WRITE(,) 'readoptprop: ',&
391	'Error while loading optical properties.'
392	CALL ABORT
393	END SELECT reading2_seq ! ==============================
394	ENDIF
395	ENDDO
396
397	! 1.4 Close the file
398
399	CLOSE(file_unit)
400
401	! 1.5 If Francois' data, convert wvl to metres
402	if(iaer.eq.iaero_co2)then
403	if(forgetit)then
404	! print*,'please sort out forgetit for naerkind>1'
405	do iwvl=1,nwvl
406	wvl(iwvl)=wvl(iwvl)*1.e-6
407	enddo
408	endif
409	endif
410
411	!$OMP END MASTER
412	!$OMP BARRIER
413
414
415
416
417
418	!==================================================================
419	! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
420	!==================================================================
421	domain: SELECT CASE (idomain)
422	!==================================================================
423	CASE(1) domain ! VISIBLE DOMAIN (idomain=1)
424	!==================================================================
425
426	lamref=lamrefvis(iaer)
427	epref=1.E+0
428
429	DO isize=1,nsize(iaer,idomain)
430
431	! Save the particle sizes
432	radiustab(iaer,idomain,isize)=radiusdyn(isize)
433
434	! Averaged optical properties (GCM channels)
435
436	! CALL aerave ( nwvl,&
437	! wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
438	! lamref,epref,tstellar,&
439	! L_NSPECTV,blamv,epavVI,&
440	! omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
441	CALL aerave_new ( nwvl,&
442	wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
443	lamref,epref,tstellar,&
444	L_NSPECTV,blamv,epavVI,&
445	omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
446
447	! Variable assignments (declared in radcommon)
448	DO isw=1,L_NSPECTV
449	QVISsQREF(isw,iaer,isize)=epavVI(isw)
450	gvis(isw,iaer,isize)=gavVI(isw)
451	omegavis(isw,iaer,isize)=omegavVI(isw)
452	END DO
453
454	ENDDO
455
456	!==================================================================
457	CASE(2) domain ! INFRARED DOMAIN (idomain=2)
458	!==================================================================
459
460
461	DO isize=1,nsize(iaer,idomain) ! ----------------------------------
462
463	lamref=lamrefir(iaer)
464	epref=1.E+0
465
466	! Save the particle sizes
467	radiustab(iaer,idomain,isize)=radiusdyn(isize)
468
469	! Averaged optical properties (GCM channels)
470
471	! epav is <QIR>/Qext(lamrefir) since epref=1
472	! Note: aerave also computes the extinction coefficient at
473	! the reference wavelength. This is called QREFvis or QREFir
474	! (not epref, which is a different parameter).
475	! Reference wavelengths SHOULD BE defined for each aerosol in
476	! radcommon_h.F90
477
478	! CALL aerave ( nwvl,&
479	! wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
480	! lamref,epref,tplanet,&
481	! L_NSPECTI,blami,epavIR,&
482	! omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
483	CALL aerave_new ( nwvl,&
484	wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
485	lamref,epref,tplanet,&
486	L_NSPECTI,blami,epavIR,&
487	omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
488
489
490	! Variable assignments (declared in radcommon)
491	DO ilw=1,L_NSPECTI
492	QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
493	gir(ilw,iaer,isize)=gavIR(ilw)
494	omegair(ilw,iaer,isize)=omegavIR(ilw)
495	END DO
496
497
498	ENDDO ! isize (particle size) -------------------------------------
499
500	END SELECT domain
501
502	!========================================================================
503	! 3. Deallocate temporary variables that were read in the ASCII files
504	!========================================================================
505
506	!$OMP BARRIER
507	!$OMP MASTER
508	IF (ALLOCATED(wvl)) DEALLOCATE(wvl) ! wvl
509	IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn) ! radiusdyn
510	IF (ALLOCATED(ep)) DEALLOCATE(ep) ! ep
511	IF (ALLOCATED(omeg)) DEALLOCATE(omeg) ! omeg
512	IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor) ! g
513	!$OMP END MASTER
514	!$OMP BARRIER
515
516	END DO ! Loop on iaer
517	END DO ! Loop on idomain
518	!========================================================================
519
520	RETURN
521
522
523
524	END subroutine suaer_corrk
525

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

Download in other formats:

Original Format