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

Last change on this file since 1242 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)

EM

File size: 18.4 KB
Line 
1SUBROUTINE suaer
2use dimradmars_mod, only: longrefvis, longrefir, nsizemax, long1vis, &
3                         long2vis, long3vis, long1ir, long2ir, long1co2, &
4                         long2co2, nsun, nir
5use yomaer_h, only: radiustab, gvis, omegavis, QVISsQREF, gIR, omegaIR, &
6                    QIRsQREF, QREFvis, QREFir, omegaREFvis, omegaREFir, &
7                    nsize
8IMPLICIT 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)
51INTEGER          :: nwvl     ! Number of wavelengths in
52                             ! the domain (VIS or IR)
53REAL, DIMENSION(:),&
54  ALLOCATABLE, SAVE :: wvl       ! Wavelength axis
55REAL, DIMENSION(:),&
56  ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis
57
58REAL, DIMENSION(:,:),&
59  ALLOCATABLE, SAVE :: ep,&    ! Extinction coefficient Qext
60                       omeg,&  ! Single Scattering Albedo
61                       gfactor ! Assymetry Factor
62
63! Local variables:
64
65INTEGER :: iaer                ! Aerosol index
66INTEGER :: idomain             ! Domain index (1=VIS,2=IR)
67INTEGER :: 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
72INTEGER :: isun                ! Solar band index
73INTEGER :: isize               ! Particle size index
74INTEGER :: jfile               ! ASCII file scan index
75INTEGER :: file_unit = 60
76LOGICAL :: file_ok, endwhile
77CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
78INTEGER :: read_ok
79
80! I/O  of "aerave" (subroutine averaging spectrally
81!   sing.scat.parameters)
82
83REAL tsun            ! Sun brightness temperature (for SW)
84REAL tsol            ! Surface reference brightness temp (LW)
85REAL longref         ! reference wavelengths
86REAL longsun(nsun+1) ! solar band boundaries
87REAL longir(nir+1)   ! IR band boundaries
88REAL epref           ! reference extinction ep
89                     ! at wavelength "longref"
90REAL epav(nir)       ! average ep
91                     ! (= <Qext>/Qext(longref) if epref=1)
92REAL omegav(nir)     ! Average sing.scat.albedo
93REAL gav(nir)        ! Average assymetry parameter
94
95! Local saved variables:
96
97CHARACTER(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
173radiustab(1:naerkind,1:2,1:nsizemax)=0
174
175gvis(1:nsun,1:naerkind,1:nsizemax)=0
176omegavis(1:nsun,1:naerkind,1:nsizemax)=0
177QVISsQREF(1:nsun,1:naerkind,1:nsizemax)=0
178
179gIR(1:nir,1:naerkind,1:nsizemax)=0
180omegaIR(1:nir,1:naerkind,1:nsizemax)=0
181QIRsQREF(1:nir,1:naerkind,1:nsizemax)=0
182
183QREFvis(1:naerkind,1:nsizemax)=0
184QREFir(1:naerkind,1:nsizemax)=0
185omegaREFvis(1:naerkind,1:nsizemax)=0
186omegaREFir(1:naerkind,1:nsizemax)=0
187
188DO 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
196INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//&
197  '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
198  EXIST=file_ok)
199IF(.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
216ENDIF
217OPEN(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
224jfile = 1
225endwhile = .false.
226DO 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
266ENDDO
267
268ALLOCATE(wvl(nwvl))                            ! wvl
269ALLOCATE(radiusdyn(nsize(iaer,idomain)))       ! radiusdyn
270ALLOCATE(ep(nwvl,nsize(iaer,idomain)))         ! ep
271ALLOCATE(omeg(nwvl,nsize(iaer,idomain)))       ! omeg
272ALLOCATE(gfactor(nwvl,nsize(iaer,idomain)))    ! g
273
274!       1.3 Read the data
275
276jfile = 1
277endwhile = .false.
278DO 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
332ENDDO
333
334!       1.4 Close the file
335
336CLOSE(file_unit)
337
338!==================================================================
339! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
340!==================================================================
341domain: SELECT CASE (idomain)
342!==================================================================
343CASE(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
354DO 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! ------------------------------------------------
400ENDDO
401
402!==================================================================
403CASE(2) domain !                  INFRARED DOMAIN (idomain=2)
404!==================================================================
405
406DO 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
495ENDDO ! isize (particle size) -------------------------------------
496
497END SELECT domain
498!==================================================================
499! 3. Deallocate temporary variables (read in ASCII files)
500!==================================================================
501
502DEALLOCATE(wvl)        ! wvl
503DEALLOCATE(radiusdyn)  ! radiusdyn
504DEALLOCATE(ep)         ! ep
505DEALLOCATE(omeg)       ! omeg
506DEALLOCATE(gfactor)    ! g
507
508  END DO ! Loop on iaer
509END DO   ! Loop on idomain
510!==================================================================
511RETURN
512END
Note: See TracBrowser for help on using the repository browser.