source: trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90 @ 799

Last change on this file since 799 was 726, checked in by jleconte, 12 years ago

17/07/2012 == JL for LK

  • Generalization of aerosol scheme:
    • any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
    • addition of a module with the id numbers for aerosols (aerosol_mod.F90).
    • initialization of aerosols id numbers in iniaerosol.F90
    • compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
  • For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

  • If starting from an old start file, recreate start file with the q=0 option in newstart.e.
  • update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

  • To add other aerosols, see Laura Kerber.
File size: 14.8 KB
Line 
1      subroutine suaer_corrk
2
3      ! inputs
4      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,iim,jjm,naerkind
5      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
6      use datafile_mod, only: datadir
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,omegarefvis,omegarefir
12      use aerosol_mod
13
14      implicit none
15
16!==================================================================
17!     Purpose
18!     -------
19!     Initialize all radiative aerosol properties
20!
21!     Notes
22!     -----
23!     Reads the optical properties -> Mean  -> Variable assignment
24!     (ASCII files)                  (see radcommon_h.F90)
25!     wvl(nwvl)                      longsun
26!     ep(nwvl)                       epav     QVISsQREF(L_NSPECTV)
27!     omeg(nwvl)                     omegav   omegavis(L_NSPECTV)
28!     gfactor(nwvl)                  gav      gvis(L_NSPECTV)
29!     
30!     Authors
31!     -------
32!     Richard Fournier (1996) Francois Forget (1996)
33!     Frederic Hourdin
34!     Jean-jacques morcrette *ECMWF*
35!     MODIF Francois Forget (2000)
36!     MODIF Franck Montmessin (add water ice)
37!     MODIF J.-B. Madeleine 2008W27
38!     - Optical properties read in ASCII files
39!     - Add varying radius of the particules
40!     - Add water ice clouds
41!     MODIF R. Wordsworth (2009)
42!     - generalisation to arbitrary spectral bands
43!
44!==================================================================
45
46#include "callkeys.h"
47!     Optical properties (read in external ASCII files)
48      INTEGER          :: nwvl  ! Number of wavelengths in
49                                ! the domain (VIS or IR)
50
51!      REAL             :: solsir ! visible to infrared ratio
52!                                ! (tau_0.67um/tau_9um)
53
54      REAL, DIMENSION(:),&
55      ALLOCATABLE, SAVE :: wvl  ! Wavelength axis
56      REAL, DIMENSION(:),&
57      ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis
58
59      REAL, DIMENSION(:,:),&
60      ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext
61      omeg,&                    ! Single Scattering Albedo
62      gfactor                   ! Assymetry Factor
63
64!     Local variables:
65
66      INTEGER :: iaer           ! Aerosol index
67      INTEGER :: idomain        ! Domain index (1=VIS,2=IR)
68      INTEGER :: ilw            ! longwave index
69      INTEGER :: isw           ! shortwave index
70      INTEGER :: isize          ! Particle size index
71      INTEGER :: jfile          ! ASCII file scan index
72      INTEGER :: file_unit = 60
73      LOGICAL :: file_ok, endwhile
74      CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
75
76!     I/O  of "aerave" (subroutine that spectrally averages
77!     the single scattering parameters)
78
79      REAL lamref                      ! reference wavelengths
80      REAL epref                       ! reference extinction ep
81
82!      REAL epav(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
83!      REAL omegav(L_NSPECTI)          ! Average single scattering albedo
84!      REAL gav(L_NSPECTI)             ! Average assymetry parameter
85
86      REAL epavVI(L_NSPECTV)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
87      REAL omegavVI(L_NSPECTV)          ! Average single scattering albedo
88      REAL gavVI(L_NSPECTV)             ! Average assymetry parameter
89
90      REAL epavIR(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
91      REAL omegavIR(L_NSPECTI)          ! Average single scattering albedo
92      REAL gavIR(L_NSPECTI)             ! Average assymetry parameter
93     
94      logical forgetit                  ! use francois' data?
95      integer iwvl
96
97!     Local saved variables:
98
99      CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id
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
104       if (iaer.eq.iaero_co2)then
105        forgetit=.true.
106        print*, 'naerkind= co2', iaer
107
108!     visible
109        if(forgetit)then
110           file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
111        else
112           file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
113        endif
114        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
115
116!     infrared
117        if(forgetit)then
118           file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
119        else
120           file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
121        endif
122        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
123       endif ! CO2 aerosols 
124! NOTE: these lamref's are currently for dust, not CO2 ice.
125! JB thinks this shouldn't matter too much, but it needs to be
126! fixed / decided for the final version.
127
128       if (iaer.eq.iaero_h2o) then
129        print*, 'naerkind= h2o', iaer
130
131!     visible
132         file_id(iaer,1) = 'optprop_icevis_n50.dat'
133         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
134!     infrared
135         file_id(iaer,2) = 'optprop_iceir_n50.dat'
136         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
137       endif
138
139       if (iaer.eq.iaero_dust) then
140        print*, 'naerkind= dust', iaer
141
142!     visible
143         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
144         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
145         lamrefvis(iaer)=0.67e-6
146!     infrared
147         file_id(iaer,2) = 'optprop_dustir_n50.dat'
148         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
149         lamrefir(iaer)=9.3E-6
150       endif
151
152       if (iaer.eq.iaero_h2so4) then
153         print*, 'naerkind= h2so4', iaer
154
155!     visible
156         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
157         !lamrefir(iaer)=  doesn't exist?
158         lamrefvis(iaer)=1.5E-6   ! no idea, must find
159!     infrared
160         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
161         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
162         lamrefir(iaer)=9.3E-6 ! no idea, must find
163! added by LK
164       endif
165       
166      enddo
167
168!------------------------------------------------------------------
169
170!     Initializations:
171
172      radiustab(:,:,:) = 0.0
173      gvis(:,:,:)      = 0.0
174      omegavis(:,:,:)  = 0.0
175      QVISsQREF(:,:,:) = 0.0
176      gir(:,:,:)       = 0.0
177      omegair(:,:,:)   = 0.0
178      QIRsQREF(:,:,:)  = 0.0
179
180      DO 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
188            INQUIRE(FILE=TRIM(datadir)//&
189         '/'//TRIM(file_id(iaer,idomain)),&
190            EXIST=file_ok)
191            IF(.NOT.file_ok) THEN
192               write(*,*)'suaer_corrk: Problem opening ',&
193               TRIM(file_id(iaer,idomain))
194               write(*,*)'It should be in: ',TRIM(datadir)
195               write(*,*)'1) You can set this directory address ',&
196               'in callphys.def with:'
197               write(*,*)' datadir = /absolute/path/to/datagcm'
198               write(*,*)'2) If ',&
199              TRIM(file_id(iaer,idomain)),&
200               ' is a LMD reference datafile, it'
201               write(*,*)' can be obtained online at:'
202               write(*,*)' http://www.lmd.jussieu.fr/',&
203               '~forget/datagcm/datafile'
204               CALL ABORT
205            ENDIF
206            OPEN(UNIT=file_unit,&
207            FILE=TRIM(datadir)//&
208         '/'//TRIM(file_id(iaer,idomain)),&
209            FORM='formatted')
210
211!     1.2 Allocate the optical property table
212
213            jfile = 1
214            endwhile = .false.
215            DO WHILE (.NOT.endwhile)
216               READ(file_unit,*) scanline
217               IF ((scanline(1:1) .ne. '#').and.&
218               (scanline(1:1) .ne. ' ')) THEN
219               BACKSPACE(file_unit)
220               reading1_seq: SELECT CASE (jfile) ! ====================
221            CASE(1) reading1_seq ! nwvl ----------------------------
222               read(file_unit,*) nwvl
223               jfile = jfile+1
224            CASE(2) reading1_seq ! nsize ---------------------------
225               read(file_unit,*) nsize(iaer,idomain)
226               endwhile = .true.
227            CASE DEFAULT reading1_seq ! ----------------------------
228               WRITE(*,*) 'readoptprop: ',&
229               'Error while loading optical properties.'
230               CALL ABORT
231            END SELECT reading1_seq ! ==============================
232         ENDIF
233      ENDDO
234
235      ALLOCATE(wvl(nwvl))       ! wvl
236      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
237      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
238      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
239      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
240
241
242!     1.3 Read the data
243
244      jfile = 1
245      endwhile = .false.
246      DO WHILE (.NOT.endwhile)
247         READ(file_unit,*) scanline
248         IF ((scanline(1:1) .ne. '#').and.&
249         (scanline(1:1) .ne. ' ')) THEN
250         BACKSPACE(file_unit)
251         reading2_seq: SELECT CASE (jfile) ! ====================
252      CASE(1) reading2_seq      ! wvl -----------------------------
253         read(file_unit,*) wvl
254         jfile = jfile+1
255      CASE(2) reading2_seq      ! radiusdyn -----------------------
256         read(file_unit,*) radiusdyn
257         jfile = jfile+1
258      CASE(3) reading2_seq      ! ep ------------------------------
259         isize = 1
260         DO WHILE (isize .le. nsize(iaer,idomain))
261            READ(file_unit,*) scanline
262            IF ((scanline(1:1) .ne. '#').and.&
263            (scanline(1:1) .ne. ' ')) THEN
264            BACKSPACE(file_unit)
265            read(file_unit,*) ep(:,isize)
266            isize = isize + 1
267         ENDIF
268      ENDDO
269
270      jfile = jfile+1
271      CASE(4) reading2_seq      ! omeg ----------------------------
272         isize = 1
273         DO WHILE (isize .le. nsize(iaer,idomain))
274            READ(file_unit,*) scanline
275            IF ((scanline(1:1) .ne. '#').and.&
276            (scanline(1:1) .ne. ' ')) THEN
277            BACKSPACE(file_unit)
278            read(file_unit,*) omeg(:,isize)
279            isize = isize + 1
280         ENDIF
281      ENDDO
282
283      jfile = jfile+1
284      CASE(5) reading2_seq      ! gfactor -------------------------
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,*) gfactor(:,isize)
292            isize = isize + 1
293         ENDIF
294      ENDDO
295
296      jfile = jfile+1
297      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
298         endwhile = .true.
299      ENDIF
300      CASE(6) reading2_seq
301         endwhile = .true.
302      CASE DEFAULT reading2_seq ! ----------------------------
303         WRITE(*,*) 'readoptprop: ',&
304         'Error while loading optical properties.'
305         CALL ABORT
306      END SELECT reading2_seq   ! ==============================
307      ENDIF
308      ENDDO
309
310!     1.4 Close the file
311
312      CLOSE(file_unit)
313
314!     1.5 If Francois' data, convert wvl to metres
315       if(iaer.eq.iaero_co2)then
316         if(forgetit)then
317         !   print*,'please sort out forgetit for naerkind>1'
318            do iwvl=1,nwvl
319               wvl(iwvl)=wvl(iwvl)*1.e-6
320            enddo
321         endif
322      endif
323
324
325
326
327
328
329!==================================================================
330!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
331!==================================================================
332      domain: SELECT CASE (idomain)
333!==================================================================
334      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
335!==================================================================
336
337         lamref=lamrefvis(iaer)
338         epref=1.E+0
339
340         DO isize=1,nsize(iaer,idomain)
341
342!     Save the particle sizes
343            radiustab(iaer,idomain,isize)=radiusdyn(isize)
344
345!     Averaged optical properties (GCM channels)
346
347!            CALL aerave ( nwvl,&
348!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
349!            lamref,epref,tstellar,&
350!            L_NSPECTV,blamv,epavVI,&
351!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
352            CALL aerave_new ( nwvl,&
353            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
354            lamref,epref,tstellar,&
355            L_NSPECTV,blamv,epavVI,&
356            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
357
358!     Variable assignments (declared in radcommon)
359            DO isw=1,L_NSPECTV
360               QVISsQREF(isw,iaer,isize)=epavVI(isw)
361               gvis(isw,iaer,isize)=gavVI(isw)
362               omegavis(isw,iaer,isize)=omegavVI(isw)
363            END DO
364
365         ENDDO
366
367!==================================================================
368      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
369!==================================================================
370
371
372         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
373
374            lamref=lamrefir(iaer)
375            epref=1.E+0
376
377!     Save the particle sizes
378            radiustab(iaer,idomain,isize)=radiusdyn(isize)
379
380!     Averaged optical properties (GCM channels)
381
382!     epav is <QIR>/Qext(lamrefir) since epref=1
383!     Note: aerave also computes the extinction coefficient at
384!     the reference wavelength. This is called QREFvis or QREFir
385!     (not epref, which is a different parameter).
386!     Reference wavelengths SHOULD BE defined for each aerosol in
387!     radcommon_h.F90
388
389!            CALL aerave ( nwvl,&
390!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
391!            lamref,epref,tplanet,&
392!            L_NSPECTI,blami,epavIR,&
393!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
394            CALL aerave_new ( nwvl,&
395            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
396            lamref,epref,tplanet,&
397            L_NSPECTI,blami,epavIR,&
398            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
399
400
401!     Variable assignments (declared in radcommon)
402            DO ilw=1,L_NSPECTI
403               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
404               gir(ilw,iaer,isize)=gavIR(ilw)
405               omegair(ilw,iaer,isize)=omegavIR(ilw)
406            END DO
407
408
409      ENDDO                     ! isize (particle size) -------------------------------------
410
411      END SELECT domain
412
413!========================================================================
414!     3. Deallocate temporary variables that were read in the ASCII files
415!========================================================================
416
417      DEALLOCATE(wvl)           ! wvl
418      DEALLOCATE(radiusdyn)     ! radiusdyn
419      DEALLOCATE(ep)            ! ep
420      DEALLOCATE(omeg)          ! omeg
421      DEALLOCATE(gfactor)       ! g
422
423      END DO                    ! Loop on iaer
424      END DO                    ! Loop on idomain
425!========================================================================
426
427      RETURN
428
429
430
431    END subroutine suaer_corrk
432     
Note: See TracBrowser for help on using the repository browser.