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

Last change on this file since 317 was 222, checked in by aslmd, 13 years ago

MESSAGE DE JB MADELEINE

C'était juste pour vous dire que j'ai fait un peu de rangement dans le
DATAFILE concernant les propriétés de la poussière (directement dans le
répertoire /u/forget/WWW/datagcm/datafile/), car lorsque Ockert-Bell
dans le domaine solaire était utilisé avec Forget 1998 dans l'IR, le
rapport Solaire/IR "solsir" n'était pas rigoureusement égal à 2.

Du coup il faut changer les noms des fichiers dans suaer.F90

REMARQUE: Le rapport "solsir" est sorti dans les logs au firstcall, et
vous verrez qu'il n'est pas égal à 2 lorsque ces propriétés sont
utilisées. C'EST NORMAL! En effet, j'ai réglé la longueur d'onde de
référence de la poussière dans l'IR à 9.3um pour pouvoir comparer
l'opacité donnée en output avec TES. Or le fameux solsir=2 utilisé
partout date de l'époque où la longueur de référence était celle du
canal "T9" d'IRTM qui est 9um et non 9.3um!

En résumé, au firstcall, le GCM affiche Qext(0.67um)/Qext(9.3um), et non
l'ancien solsir qui lui vaut Qext(0.67um)/Qext(9um), d'où un résultat
qui n'est pas égal à deux. Mais rassurez-vous, Qext(0.67um)/Qext(9um)
dans les fichiers d'input vaut bien 2.

Également, j'ai un peu modifié aeropacity.F car il affichait le
"ccn_factor" à tout va dans les logs. Il est également attaché à ce mail.
File size: 17.9 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_clancy.dat' !Clancy-Lee
103!       file_id(iaer,1) = 'optprop_dustvis_ockert.dat' !Ockert-Bell
104!       Infrared domain:
105        file_id(iaer,2) = 'optprop_dustir_TM.dat'      !M.Wolff TM
106!       Toon-Forget + solsir=2 using Clancy-Lee
107!       file_id(iaer,2) = 'optprop_dustir_clancy.dat'
108!       Toon-Forget + solsir=2 using Ockert-Bell
109!       file_id(iaer,2) = 'optprop_dustir_ockert.dat'
110!       Reference wavelength in the visible:
111        longrefvis(iaer)=0.67E-6
112!                     For dust: change readtesassim accordingly;
113!       Reference wavelength in the infrared:
114        longrefir(iaer)=9.3E-6
115!==================================================================
116        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
117!==================================================================
118!       Visible domain:
119        file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
120!       file_id(iaer,1) = 'optprop_dustvis_n50.dat'    !Mie
121!       Infrared domain:
122        file_id(iaer,2) = 'optprop_dustir_n50.dat'     !Mie
123!       Reference wavelength in the visible:
124        longrefvis(iaer)=0.67E-6
125!       If not equal to 0.67e-6 -> change readtesassim accordingly;
126!       Reference wavelength in the infrared:
127        longrefir(iaer)=9.3E-6
128!==================================================================
129        CASE("dust_submicron") aerkind   ! Small dust population
130!==================================================================
131!       Visible domain:
132        file_id(iaer,1) = 'optprop_dustvis_01um_TM.dat' !M.Wolff
133!       Infrared domain:
134        file_id(iaer,2) = 'optprop_dustir_01um_TM.dat'  !M.Wolff
135!       Reference wavelength in the visible:
136        longrefvis(iaer)=0.67E-6
137!       If not equal to 0.67e-6 -> change readtesassim accordingly;
138!       Reference wavelength in the infrared:
139        longrefir(iaer)=9.3E-6
140!==================================================================
141        CASE("h2o_ice") aerkind             ! Water ice crystals
142!==================================================================
143!       Visible domain:
144        file_id(iaer,1) = 'optprop_icevis_n30.dat' !Warren
145!       file_id(iaer,1) = 'optprop_icevis.dat'     !Warren
146!       Infrared domain:
147        file_id(iaer,2) = 'optprop_iceir_n30.dat'  !Warren
148!       file_id(iaer,2) = 'optprop_iceir.dat'      !Warren
149!       Reference wavelength in the visible:
150        longrefvis(iaer)=0.67E-6  ! 1.5um OMEGA/MEx
151!       Reference wavelength in the infrared:
152        longrefir(iaer)=12.1E-6  ! 825cm-1 TES/MGS
153!==================================================================
154        END SELECT aerkind
155!==================================================================
156        WRITE(*,*) "Scatterer: ",trim(name_iaer(iaer))
157        WRITE(*,*) "  corresponding files: "
158        WRITE(*,*) "VIS: ",trim(file_id(iaer,1))
159        WRITE(*,*) "IR : ",trim(file_id(iaer,2))
160!==================================================================
161      ENDDO ! iaer (loop on aerosol kind)
162
163! Initializations:
164
165radiustab(1:naerkind,1:2,1:nsizemax)=0
166
167gvis(1:nsun,1:naerkind,1:nsizemax)=0
168omegavis(1:nsun,1:naerkind,1:nsizemax)=0
169QVISsQREF(1:nsun,1:naerkind,1:nsizemax)=0
170
171gIR(1:nir,1:naerkind,1:nsizemax)=0
172omegaIR(1:nir,1:naerkind,1:nsizemax)=0
173QIRsQREF(1:nir,1:naerkind,1:nsizemax)=0
174
175QREFvis(1:naerkind,1:nsizemax)=0
176QREFir(1:naerkind,1:nsizemax)=0
177omegaREFvis(1:naerkind,1:nsizemax)=0
178omegaREFir(1:naerkind,1:nsizemax)=0
179
180DO iaer = 1, naerkind ! Loop on aerosol kind
181  DO idomain = 1, 2   ! Loop on radiation domain (VIS or IR)
182!==================================================================
183! 1. READ OPTICAL PROPERTIES
184!==================================================================
185
186!       1.1 Open the ASCII file
187
188INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//&
189  '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
190  EXIST=file_ok)
191IF(.NOT.file_ok) THEN
192  write(*,*)'Problem opening ',&
193    file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))
194  write(*,*)'It should be in: ',&
195    datafile(1:LEN_TRIM(datafile))
196  write(*,*)'1) You can change this directory address in '
197  write(*,*)' file phymars/datafile.h'
198  write(*,*)'2) If ',&
199    file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
200    ' is a LMD reference datafile, it'
201  write(*,*)' can be obtained online on:'
202  write(*,*)' http://www.lmd.jussieu.fr/',&
203    '~forget/datagcm/datafile'
204  write(*,*)'3) If the name of the file is wrong, you can'
205  write(*,*)' change it in file phymars/suaer.F90. Just'
206  write(*,*)' modify the variable called file_id.'
207  CALL ABORT
208ENDIF
209OPEN(UNIT=file_unit,&
210  FILE=datafile(1:LEN_TRIM(datafile))//&
211  '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
212  FORM='formatted')
213
214!       1.2 Allocate the optical property table
215
216jfile = 1
217endwhile = .false.
218DO WHILE (.NOT.endwhile)
219  READ(file_unit,*,iostat=read_ok) scanline
220  if (read_ok.ne.0) then
221    write(*,*)' readoptprop: Error reading file',&
222    trim(datafile(1:LEN_TRIM(datafile))//&
223    '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))))
224    call abort
225  endif
226  IF ((scanline(1:1) .ne. '#').and.&
227    (scanline(1:1) .ne. ' ')) THEN
228    BACKSPACE(file_unit)
229    reading1_seq: SELECT CASE (jfile) ! ====================
230    CASE(1) reading1_seq ! nwvl ----------------------------
231        read(file_unit,*,iostat=read_ok) nwvl
232        if (read_ok.ne.0) then
233          write(*,*)' readoptprop: Error while reading line:',&
234          trim(scanline)
235          write(*,*)'   of file',&
236          trim(datafile(1:LEN_TRIM(datafile))//&
237          '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))))
238          call abort
239        endif
240        jfile = jfile+1
241    CASE(2) reading1_seq ! nsize ---------------------------
242        read(file_unit,*,iostat=read_ok) nsize(iaer,idomain)
243        if (read_ok.ne.0) then
244          write(*,*)' readoptprop: Error while reading line:',&
245          trim(scanline)
246          write(*,*)'   of file',&
247          trim(datafile(1:LEN_TRIM(datafile))//&
248          '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))))
249          call abort
250        endif
251        endwhile = .true.
252    CASE DEFAULT reading1_seq ! ----------------------------
253        WRITE(*,*) 'readoptprop: ',&
254          'Error while loading optical properties.'
255        CALL ABORT
256    END SELECT reading1_seq ! ==============================
257  ENDIF
258ENDDO
259
260ALLOCATE(wvl(nwvl))                            ! wvl
261ALLOCATE(radiusdyn(nsize(iaer,idomain)))       ! radiusdyn
262ALLOCATE(ep(nwvl,nsize(iaer,idomain)))         ! ep
263ALLOCATE(omeg(nwvl,nsize(iaer,idomain)))       ! omeg
264ALLOCATE(gfactor(nwvl,nsize(iaer,idomain)))    ! g
265
266!       1.3 Read the data
267
268jfile = 1
269endwhile = .false.
270DO WHILE (.NOT.endwhile)
271   READ(file_unit,*) scanline
272  IF ((scanline(1:1) .ne. '#').and.&
273    (scanline(1:1) .ne. ' ')) THEN
274    BACKSPACE(file_unit)
275    reading2_seq: SELECT CASE (jfile) ! ====================
276    CASE(1) reading2_seq ! wvl -----------------------------
277        read(file_unit,*) wvl
278        jfile = jfile+1
279    CASE(2) reading2_seq ! radiusdyn -----------------------
280        read(file_unit,*) radiusdyn
281        jfile = jfile+1
282    CASE(3) reading2_seq ! ep ------------------------------
283        isize = 1
284        DO WHILE (isize .le. nsize(iaer,idomain))
285          READ(file_unit,*) scanline
286          IF ((scanline(1:1) .ne. '#').and.&
287            (scanline(1:1) .ne. ' ')) THEN
288          BACKSPACE(file_unit)
289          read(file_unit,*) ep(:,isize)
290          isize = isize + 1
291          ENDIF
292        ENDDO
293        jfile = jfile+1
294    CASE(4) reading2_seq ! omeg ----------------------------
295        isize = 1
296        DO WHILE (isize .le. nsize(iaer,idomain))
297          READ(file_unit,*) scanline
298          IF ((scanline(1:1) .ne. '#').and.&
299            (scanline(1:1) .ne. ' ')) THEN
300          BACKSPACE(file_unit)
301          read(file_unit,*) omeg(:,isize)
302          isize = isize + 1
303          ENDIF
304        ENDDO
305        jfile = jfile+1
306    CASE(5) reading2_seq ! gfactor -------------------------
307        isize = 1
308        DO WHILE (isize .le. nsize(iaer,idomain))
309          READ(file_unit,*) scanline
310          IF ((scanline(1:1) .ne. '#').and.&
311            (scanline(1:1) .ne. ' ')) THEN
312          BACKSPACE(file_unit)
313          read(file_unit,*) gfactor(:,isize)
314          isize = isize + 1
315          ENDIF
316        ENDDO
317        endwhile = .true.
318    CASE DEFAULT reading2_seq ! ----------------------------
319        WRITE(*,*) 'suaer.F90: ',&
320          'Error while loading optical properties.'
321        CALL ABORT
322    END SELECT reading2_seq ! ==============================
323  ENDIF
324ENDDO
325
326!       1.4 Close the file
327
328CLOSE(file_unit)
329
330!==================================================================
331! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
332!==================================================================
333domain: SELECT CASE (idomain)
334!==================================================================
335CASE(1) domain !                   VISIBLE DOMAIN (idomain=1)
336!==================================================================
337
338! 2.1 Parameters
339  tsun=6000.E+0
340  longsun(1)=long1vis
341  longsun(2)=long2vis
342  longsun(3)=long3vis
343  longref=longrefvis(iaer)
344  epref=1.E+0
345
346DO isize=1,nsize(iaer,idomain)
347! test that there is enough room to store the data
348 if (isize.gt.nsizemax) then
349   write(*,*) "suaer: Error ! nsizemax is too small!"
350   write(*,*) "       nsizemax=",nsizemax
351   write(*,*) "       you must increase the value of nsizemax"
352   write(*,*) "       in dimradmars.h !"
353   stop
354 endif
355! ------------------------------------------------
356! 2.2 Save the particle sizes
357  radiustab(iaer,idomain,isize)=radiusdyn(isize)
358! 2.3 Averaged optical properties (GCM channels)
359! Notice: Aerave also computes the extinction coefficient and
360!   single scattering albedo at reference wavelength
361!   (called QREFvis and OMEGAREFvis, same in the IR,
362!   and not epref, which is a different parameter);
363!   Reference wavelengths are defined for each aerosol in
364!   dimradmars.h.
365
366  CALL aerave ( nwvl,&
367       wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
368       longref,epref,tsun,&
369       nsun,longsun, epav,omegav,gav,&
370       QREFvis(iaer,isize),omegaREFvis(iaer,isize) )
371! 2.4 Variable assignements (declared by yomaer.h)
372  DO isun=1,nsun
373    QVISsQREF(isun,iaer,isize)=epav(isun)
374    gvis(isun,iaer,isize)=gav(isun)
375    omegavis(isun,iaer,isize)=omegav(isun)
376  END DO
377! 2.5 Output display
378  WRITE(*,*) 'Les donnees spectrales :'
379  WRITE(*,*) 'Solaire (SW) ---->'
380  WRITE(*,*) 'Aerosol number: ', iaer
381  WRITE(*,*) 'Rayon aerosol: ', radiustab(iaer,idomain,isize)
382  WRITE(*,*) '<Qext>/Qext(longrefvis) ; omega ; g'
383  DO isun=1,nsun
384    WRITE(*,*) QVISsQREF(isun,iaer,isize),&
385         omegavis(isun,iaer,isize),&
386         gvis(isun,iaer,isize)
387  ENDDO
388  WRITE(*,*) 'QREFvis(',iaer,isize,') = ',QREFvis(iaer,isize)
389  WRITE(*,*) 'omegaREFvis(',iaer,isize,') = ',&
390                                      omegaREFvis(iaer,isize)
391! ------------------------------------------------
392ENDDO
393
394!==================================================================
395CASE(2) domain !                  INFRARED DOMAIN (idomain=2)
396!==================================================================
397
398DO isize=1,nsize(iaer,idomain) ! ----------------------------------
399
400! 2.1 solsir is not used anymore; division of Qext(IR) by solsir
401!     has to be done in the input ASCII files (if necessary).
402
403! 2.2 Save the particle sizes
404  radiustab(iaer,idomain,isize)=radiusdyn(isize)
405
406! 2.3 Parameters
407
408  tsol=215.D+0
409  longir(1)=long1ir
410  longir(2)=long1co2
411  longir(3)=long2co2
412  longir(4)=long2ir
413  longref=longrefir(iaer)
414  epref=1.E+0
415
416! 2.4 Averaged optical properties (GCM channels)
417!           epav is <QIR>/Qext(longrefir) since epref=1
418! Notice: Aerave also Computes the extinction coefficient at
419!   reference wavelength (called QREFvis or QREFir,
420!   and not epref, which is a different parameter);
421!   Reference wavelengths are defined for each aerosol in
422!   dimradmars.h.
423
424  CALL aerave ( nwvl,&
425       wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
426       longref,epref,tsol,&
427       nir-1,longir,epav,omegav,gav,&
428       QREFir(iaer,isize),omegaREFir(iaer,isize) )
429  WRITE(*,*) 'QREFir(',iaer,isize,') = ',QREFir(iaer,isize)
430  WRITE(*,*) 'omegaREFir(',iaer,isize,') = ',&
431                                      omegaREFir(iaer,isize)
432
433! 2.5 Computing  <QIR>/Qext(longrefvis)
434
435  DO iir=1,nir-1
436    WRITE(*,*) 'QIRsQREFir Channel ',iir,': ',epav(iir)
437    epav(iir)=  epav(iir) * QREFir(iaer,isize) / &
438                            QREFvis(iaer,isize)
439  ENDDO
440  WRITE(*,*) 'Aerosol number', iaer
441  WRITE(*,*) 'Particle size: ',radiustab(iaer,idomain,isize)
442  WRITE(*,*) 'Rapport Solaire/IR:',&
443             QREFvis(iaer,isize) / QREFir(iaer,isize)
444
445! 2.6 Variable assignements
446!           (variables are declared by yomaer.h)
447
448!         Single scattering properties
449!           in each of the "nir" bands
450!           (cf. dimradmars.h)
451
452! iir=1 : central 15um CO2 bands   
453  QIRsQREF(1,iaer,isize)=epav(2)
454  omegaIR(1,iaer,isize)=omegav(2)
455  gIR(1,iaer,isize)=gav(2)
456
457! iir=2 : CO2 band wings
458!           (same properties than for central part)
459  QIRsQREF(2,iaer,isize)=epav(2)
460  omegaIR(2,iaer,isize)=omegav(2)
461  gIR(2,iaer,isize)=gav(2)
462
463! iir=3 : 9 um band [long1ir - long1co2]
464  QIRsQREF(3,iaer,isize)=epav(1)
465  omegaIR(3,iaer,isize)=omegav(1)
466  gIR(3,iaer,isize)=gav(1)
467
468! iir=4 : Far IR    [long2co2 - long2ir]
469  QIRsQREF(4,iaer,isize)=epav(3)
470  omegaIR(4,iaer,isize)=omegav(3)
471  gIR(4,iaer,isize)=gav(3)
472
473! 2.7 Output display
474
475  WRITE(*,*) 'AEROSOL PROPERTIES: Number ',iaer
476  WRITE(*,*) 'Thermal IR (LW) ---->'
477  WRITE(*,*) 'Particle size: ',radiustab(iaer,idomain,isize)
478  WRITE(*,*) '<Qext>/Qext(longrefvis) ; omega ; g'
479  DO iir=1,nir
480    WRITE(*,*) QIRsQREF(iir,iaer,isize),omegaIR(iir,iaer,isize),&
481         gIR(iir,iaer,isize)
482  ENDDO
483  WRITE(*,*) 'CO2: <Qabs>/Qext(longrefvis) = ',&
484       QIRsQREF(1,iaer,isize)*(1-omegaIR(1,iaer,isize))
485  WRITE(*,*) ''
486
487ENDDO ! isize (particle size) -------------------------------------
488
489END SELECT domain
490!==================================================================
491! 3. Deallocate temporary variables (read in ASCII files)
492!==================================================================
493
494DEALLOCATE(wvl)        ! wvl
495DEALLOCATE(radiusdyn)  ! radiusdyn
496DEALLOCATE(ep)         ! ep
497DEALLOCATE(omeg)       ! omeg
498DEALLOCATE(gfactor)    ! g
499
500  END DO ! Loop on iaer
501END DO   ! Loop on idomain
502!==================================================================
503RETURN
504END
Note: See TracBrowser for help on using the repository browser.