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

Last change on this file since 1266 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

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