source: trunk/LMDZ.MARS/libf/phymars/suaer.F90 @ 171

Last change on this file since 171 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

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