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

Last change on this file since 1243 was 1026, checked in by sglmd, 11 years ago

Added a flexible, 2-layer aerosol scenario (initially for Saturn, but can be simply tuned). Called aeroback2lay.

File size: 15.2 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      if (noaero) then
104        print*, 'naerkind= 0'
105      endif
106      do iaer=1,naerkind
107       if (iaer.eq.iaero_co2) then
108        forgetit=.true.
109          if (.not.noaero) then
110              print*, 'naerkind= co2', iaer
111          end if
112!     visible
113        if(forgetit)then
114           file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
115        else
116           file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
117        endif
118        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
119
120!     infrared
121        if(forgetit)then
122           file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
123        else
124           file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
125        endif
126        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
127       endif ! CO2 aerosols 
128! NOTE: these lamref's are currently for dust, not CO2 ice.
129! JB thinks this shouldn't matter too much, but it needs to be
130! fixed / decided for the final version.
131
132       if (iaer.eq.iaero_h2o) then
133        print*, 'naerkind= h2o', iaer
134
135!     visible
136         file_id(iaer,1) = 'optprop_icevis_n50.dat'
137         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
138!     infrared
139         file_id(iaer,2) = 'optprop_iceir_n50.dat'
140         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
141       endif
142
143       if (iaer.eq.iaero_dust) then
144        print*, 'naerkind= dust', iaer
145
146!     visible
147         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
148         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
149         lamrefvis(iaer)=0.67e-6
150!     infrared
151         file_id(iaer,2) = 'optprop_dustir_n50.dat'
152         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
153         lamrefir(iaer)=9.3E-6
154       endif
155
156       if (iaer.eq.iaero_h2so4) then
157         print*, 'naerkind= h2so4', iaer
158
159!     visible
160         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
161         !lamrefir(iaer)=  doesn't exist?
162         lamrefvis(iaer)=1.5E-6   ! no idea, must find
163!     infrared
164         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
165         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
166         lamrefir(iaer)=9.3E-6 ! no idea, must find
167! added by LK
168       endif
169
170       if (iaer.eq.iaero_back2lay) then
171         print*, 'naerkind= back2lay', iaer
172
173!     visible
174         file_id(iaer,1) = 'optprop_saturn_vis_n20.dat'
175         lamrefvis(iaer)=0.8E-6  !
176!     infrared
177         file_id(iaer,2) = 'optprop_saturn_ir_n20.dat'
178         lamrefir(iaer)=6.E-6  !
179! added by SG
180       endif
181       
182       
183      enddo
184
185!------------------------------------------------------------------
186
187!     Initializations:
188
189      radiustab(:,:,:) = 0.0
190      gvis(:,:,:)      = 0.0
191      omegavis(:,:,:)  = 0.0
192      QVISsQREF(:,:,:) = 0.0
193      gir(:,:,:)       = 0.0
194      omegair(:,:,:)   = 0.0
195      QIRsQREF(:,:,:)  = 0.0
196
197      DO iaer = 1, naerkind     ! Loop on aerosol kind
198         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
199!==================================================================
200!     1. READ OPTICAL PROPERTIES
201!==================================================================
202
203!     1.1 Open the ASCII file
204
205            INQUIRE(FILE=TRIM(datadir)//&
206         '/'//TRIM(file_id(iaer,idomain)),&
207            EXIST=file_ok)
208            IF(.NOT.file_ok) THEN
209               write(*,*)'suaer_corrk: Problem opening ',&
210               TRIM(file_id(iaer,idomain))
211               write(*,*)'It should be in: ',TRIM(datadir)
212               write(*,*)'1) You can set this directory address ',&
213               'in callphys.def with:'
214               write(*,*)' datadir = /absolute/path/to/datagcm'
215               write(*,*)'2) If ',&
216              TRIM(file_id(iaer,idomain)),&
217               ' is a LMD reference datafile, it'
218               write(*,*)' can be obtained online at:'
219               write(*,*)' http://www.lmd.jussieu.fr/',&
220               '~forget/datagcm/datafile'
221               CALL ABORT
222            ENDIF
223            OPEN(UNIT=file_unit,&
224            FILE=TRIM(datadir)//&
225         '/'//TRIM(file_id(iaer,idomain)),&
226            FORM='formatted')
227
228!     1.2 Allocate the optical property table
229
230            jfile = 1
231            endwhile = .false.
232            DO WHILE (.NOT.endwhile)
233               READ(file_unit,*) scanline
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,*) nwvl
240               jfile = jfile+1
241            CASE(2) reading1_seq ! nsize ---------------------------
242               read(file_unit,*) nsize(iaer,idomain)
243               endwhile = .true.
244            CASE DEFAULT reading1_seq ! ----------------------------
245               WRITE(*,*) 'readoptprop: ',&
246               'Error while loading optical properties.'
247               CALL ABORT
248            END SELECT reading1_seq ! ==============================
249         ENDIF
250      ENDDO
251
252      ALLOCATE(wvl(nwvl))       ! wvl
253      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
254      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
255      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
256      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
257
258
259!     1.3 Read the data
260
261      jfile = 1
262      endwhile = .false.
263      DO WHILE (.NOT.endwhile)
264         READ(file_unit,*) scanline
265         IF ((scanline(1:1) .ne. '#').and.&
266         (scanline(1:1) .ne. ' ')) THEN
267         BACKSPACE(file_unit)
268         reading2_seq: SELECT CASE (jfile) ! ====================
269      CASE(1) reading2_seq      ! wvl -----------------------------
270         read(file_unit,*) wvl
271         jfile = jfile+1
272      CASE(2) reading2_seq      ! radiusdyn -----------------------
273         read(file_unit,*) radiusdyn
274         jfile = jfile+1
275      CASE(3) reading2_seq      ! ep ------------------------------
276         isize = 1
277         DO WHILE (isize .le. nsize(iaer,idomain))
278            READ(file_unit,*) scanline
279            IF ((scanline(1:1) .ne. '#').and.&
280            (scanline(1:1) .ne. ' ')) THEN
281            BACKSPACE(file_unit)
282            read(file_unit,*) ep(:,isize)
283            isize = isize + 1
284         ENDIF
285      ENDDO
286
287      jfile = jfile+1
288      CASE(4) reading2_seq      ! omeg ----------------------------
289         isize = 1
290         DO WHILE (isize .le. nsize(iaer,idomain))
291            READ(file_unit,*) scanline
292            IF ((scanline(1:1) .ne. '#').and.&
293            (scanline(1:1) .ne. ' ')) THEN
294            BACKSPACE(file_unit)
295            read(file_unit,*) omeg(:,isize)
296            isize = isize + 1
297         ENDIF
298      ENDDO
299
300      jfile = jfile+1
301      CASE(5) reading2_seq      ! gfactor -------------------------
302         isize = 1
303         DO WHILE (isize .le. nsize(iaer,idomain))
304            READ(file_unit,*) scanline
305            IF ((scanline(1:1) .ne. '#').and.&
306            (scanline(1:1) .ne. ' ')) THEN
307            BACKSPACE(file_unit)
308            read(file_unit,*) gfactor(:,isize)
309            isize = isize + 1
310         ENDIF
311      ENDDO
312
313      jfile = jfile+1
314      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
315         endwhile = .true.
316      ENDIF
317      CASE(6) reading2_seq
318         endwhile = .true.
319      CASE DEFAULT reading2_seq ! ----------------------------
320         WRITE(*,*) 'readoptprop: ',&
321         'Error while loading optical properties.'
322         CALL ABORT
323      END SELECT reading2_seq   ! ==============================
324      ENDIF
325      ENDDO
326
327!     1.4 Close the file
328
329      CLOSE(file_unit)
330
331!     1.5 If Francois' data, convert wvl to metres
332       if(iaer.eq.iaero_co2)then
333         if(forgetit)then
334         !   print*,'please sort out forgetit for naerkind>1'
335            do iwvl=1,nwvl
336               wvl(iwvl)=wvl(iwvl)*1.e-6
337            enddo
338         endif
339      endif
340
341
342
343
344
345
346!==================================================================
347!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
348!==================================================================
349      domain: SELECT CASE (idomain)
350!==================================================================
351      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
352!==================================================================
353
354         lamref=lamrefvis(iaer)
355         epref=1.E+0
356
357         DO isize=1,nsize(iaer,idomain)
358
359!     Save the particle sizes
360            radiustab(iaer,idomain,isize)=radiusdyn(isize)
361
362!     Averaged optical properties (GCM channels)
363
364!            CALL aerave ( nwvl,&
365!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
366!            lamref,epref,tstellar,&
367!            L_NSPECTV,blamv,epavVI,&
368!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
369            CALL aerave_new ( nwvl,&
370            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
371            lamref,epref,tstellar,&
372            L_NSPECTV,blamv,epavVI,&
373            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
374
375!     Variable assignments (declared in radcommon)
376            DO isw=1,L_NSPECTV
377               QVISsQREF(isw,iaer,isize)=epavVI(isw)
378               gvis(isw,iaer,isize)=gavVI(isw)
379               omegavis(isw,iaer,isize)=omegavVI(isw)
380            END DO
381
382         ENDDO
383
384!==================================================================
385      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
386!==================================================================
387
388
389         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
390
391            lamref=lamrefir(iaer)
392            epref=1.E+0
393
394!     Save the particle sizes
395            radiustab(iaer,idomain,isize)=radiusdyn(isize)
396
397!     Averaged optical properties (GCM channels)
398
399!     epav is <QIR>/Qext(lamrefir) since epref=1
400!     Note: aerave also computes the extinction coefficient at
401!     the reference wavelength. This is called QREFvis or QREFir
402!     (not epref, which is a different parameter).
403!     Reference wavelengths SHOULD BE defined for each aerosol in
404!     radcommon_h.F90
405
406!            CALL aerave ( nwvl,&
407!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
408!            lamref,epref,tplanet,&
409!            L_NSPECTI,blami,epavIR,&
410!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
411            CALL aerave_new ( nwvl,&
412            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
413            lamref,epref,tplanet,&
414            L_NSPECTI,blami,epavIR,&
415            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
416
417
418!     Variable assignments (declared in radcommon)
419            DO ilw=1,L_NSPECTI
420               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
421               gir(ilw,iaer,isize)=gavIR(ilw)
422               omegair(ilw,iaer,isize)=omegavIR(ilw)
423            END DO
424
425
426      ENDDO                     ! isize (particle size) -------------------------------------
427
428      END SELECT domain
429
430!========================================================================
431!     3. Deallocate temporary variables that were read in the ASCII files
432!========================================================================
433
434      DEALLOCATE(wvl)           ! wvl
435      DEALLOCATE(radiusdyn)     ! radiusdyn
436      DEALLOCATE(ep)            ! ep
437      DEALLOCATE(omeg)          ! omeg
438      DEALLOCATE(gfactor)       ! g
439
440      END DO                    ! Loop on iaer
441      END DO                    ! Loop on idomain
442!========================================================================
443
444      RETURN
445
446
447
448    END subroutine suaer_corrk
449     
Note: See TracBrowser for help on using the repository browser.