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

Last change on this file since 253 was 253, checked in by emillour, 13 years ago

Generic GCM

  • Massive update to version 0.7

EM+RW

File size: 14.5 KB
Line 
1      subroutine suaer_corrk
2
3      ! inputs
4      use radinc_h,    only: L_NSPECTI,L_NSPECTV,naerkind,nsizemax,iim,jjm
5      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
6
7      ! outputs
8      use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
9      use radcommon_h, only: radiustab,nsize,tstellar
10      use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
11
12      implicit none
13
14!==================================================================
15!     Purpose
16!     -------
17!     Initialize all radiative aerosol properties
18!
19!     Notes
20!     -----
21!     Reads the optical properties -> Mean  -> Variable assignment
22!     (ASCII files)                  (see radcommon_h.F90)
23!     wvl(nwvl)                      longsun
24!     ep(nwvl)                       epav     QVISsQREF(L_NSPECTV)
25!     omeg(nwvl)                     omegav   omegavis(L_NSPECTV)
26!     gfactor(nwvl)                  gav      gvis(L_NSPECTV)
27!     
28!     Authors
29!     -------
30!     Richard Fournier (1996) Francois Forget (1996)
31!     Frederic Hourdin
32!     Jean-jacques morcrette *ECMWF*
33!     MODIF Francois Forget (2000)
34!     MODIF Franck Montmessin (add water ice)
35!     MODIF J.-B. Madeleine 2008W27
36!     - Optical properties read in ASCII files
37!     - Add varying radius of the particules
38!     - Add water ice clouds
39!     MODIF R. Wordsworth (2009)
40!     - generalisation to arbitrary spectral bands
41!
42!==================================================================
43
44#include "callkeys.h"
45#include "datafile.h"
46
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      forgetit=.true.
104
105!     naerkind=1
106!     visible
107      if(forgetit)then
108         file_id(1,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
109      else
110         file_id(1,1) = 'optprop_co2ice_vis_n50.dat'
111      endif
112      lamrefvis(1)=1.5E-6   ! 1.5um OMEGA/MEx ???
113
114!     infrared
115      if(forgetit)then
116         file_id(1,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
117      else
118         file_id(1,2) = 'optprop_co2ice_ir_n50.dat'
119      endif
120      lamrefir(1)=12.1E-6   ! 825cm-1 TES/MGS ???
121         
122! NOTE: these lamref's are currently for dust, not CO2 ice.
123! JB thinks this shouldn't matter too much, but it needs to be
124! fixed / decided for the final version.
125
126      IF (naerkind .gt. 1) THEN
127!     naerkind=2
128!     visible
129         file_id(2,1) = 'optprop_icevis_n50.dat'
130         lamrefvis(2)=1.5E-6   ! 1.5um OMEGA/MEx
131!     infrared
132         file_id(2,2) = 'optprop_iceir_n50.dat'
133         lamrefir(2)=12.1E-6   ! 825cm-1 TES/MGS
134      ENDIF
135
136      IF (naerkind .eq. 3) THEN
137!     naerkind=3
138!     visible
139         file_id(naerkind,1) = 'optprop_dustvis_n50.dat'
140         !lamrefvis(3)=1.5E-6   ! 1.5um OMEGA/MEx
141         lamrefvis(naerkind)=0.67e-6
142!     infrared
143         file_id(naerkind,2) = 'optprop_dustir_n50.dat'
144         !lamrefir(3)=12.1E-6   ! 825cm-1 TES/MGS
145         lamrefir(naerkind)=9.3E-6
146
147      ENDIF
148
149      IF (naerkind .gt. 3) THEN
150         print*,'naerkind = ',naerkind
151         print*,'but we only have data for 3 types, exiting.'
152         call abort
153      ENDIF
154
155!------------------------------------------------------------------
156
157!     Initializations:
158
159      radiustab(:,:,:) = 0.0
160      gvis(:,:,:)      = 0.0
161      omegavis(:,:,:)  = 0.0
162      QVISsQREF(:,:,:) = 0.0
163      gir(:,:,:)       = 0.0
164      omegair(:,:,:)   = 0.0
165      QIRsQREF(:,:,:)  = 0.0
166
167      DO iaer = 1, naerkind     ! Loop on aerosol kind
168         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
169!==================================================================
170!     1. READ OPTICAL PROPERTIES
171!==================================================================
172
173!     1.1 Open the ASCII file
174
175            INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//&
176         '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
177            EXIST=file_ok)
178            IF(.NOT.file_ok) THEN
179               write(*,*)'Problem opening ',&
180               file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))
181               write(*,*)'It should be in: ',&
182               datafile(1:LEN_TRIM(datafile))
183               write(*,*)'1) You can change this directory address in '
184               write(*,*)' file phymars/datafile.h'
185               write(*,*)'2) If ',&
186              file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
187               ' is a LMD reference datafile, it'
188               write(*,*)' can be obtained online at:'
189               write(*,*)' http://www.lmd.jussieu.fr/',&
190               '~forget/datagcm/datafile'
191               write(*,*)'3) If the name of the file is wrong, you can'
192               write(*,*)' change it in file phyuniv/suaer_corrk.F90. Just'
193               write(*,*)' modify the variable called file_id.'
194               CALL ABORT
195            ENDIF
196            OPEN(UNIT=file_unit,&
197            FILE=datafile(1:LEN_TRIM(datafile))//&
198         '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
199            FORM='formatted')
200
201!     1.2 Allocate the optical property table
202
203            jfile = 1
204            endwhile = .false.
205            DO WHILE (.NOT.endwhile)
206               READ(file_unit,*) scanline
207               IF ((scanline(1:1) .ne. '#').and.&
208               (scanline(1:1) .ne. ' ')) THEN
209               BACKSPACE(file_unit)
210               reading1_seq: SELECT CASE (jfile) ! ====================
211            CASE(1) reading1_seq ! nwvl ----------------------------
212               read(file_unit,*) nwvl
213               jfile = jfile+1
214            CASE(2) reading1_seq ! nsize ---------------------------
215               read(file_unit,*) nsize(iaer,idomain)
216               endwhile = .true.
217            CASE DEFAULT reading1_seq ! ----------------------------
218               WRITE(*,*) 'readoptprop: ',&
219               'Error while loading optical properties.'
220               CALL ABORT
221            END SELECT reading1_seq ! ==============================
222         ENDIF
223      ENDDO
224
225      ALLOCATE(wvl(nwvl))       ! wvl
226      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
227      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
228      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
229      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
230
231
232!     1.3 Read the data
233
234      jfile = 1
235      endwhile = .false.
236      DO WHILE (.NOT.endwhile)
237         READ(file_unit,*) scanline
238         IF ((scanline(1:1) .ne. '#').and.&
239         (scanline(1:1) .ne. ' ')) THEN
240         BACKSPACE(file_unit)
241         reading2_seq: SELECT CASE (jfile) ! ====================
242      CASE(1) reading2_seq      ! wvl -----------------------------
243         read(file_unit,*) wvl
244         jfile = jfile+1
245      CASE(2) reading2_seq      ! radiusdyn -----------------------
246         read(file_unit,*) radiusdyn
247         jfile = jfile+1
248      CASE(3) reading2_seq      ! ep ------------------------------
249         isize = 1
250         DO WHILE (isize .le. nsize(iaer,idomain))
251            READ(file_unit,*) scanline
252            IF ((scanline(1:1) .ne. '#').and.&
253            (scanline(1:1) .ne. ' ')) THEN
254            BACKSPACE(file_unit)
255            read(file_unit,*) ep(:,isize)
256            isize = isize + 1
257         ENDIF
258      ENDDO
259
260      jfile = jfile+1
261      CASE(4) reading2_seq      ! omeg ----------------------------
262         isize = 1
263         DO WHILE (isize .le. nsize(iaer,idomain))
264            READ(file_unit,*) scanline
265            IF ((scanline(1:1) .ne. '#').and.&
266            (scanline(1:1) .ne. ' ')) THEN
267            BACKSPACE(file_unit)
268            read(file_unit,*) omeg(:,isize)
269            isize = isize + 1
270         ENDIF
271      ENDDO
272
273      jfile = jfile+1
274      CASE(5) reading2_seq      ! gfactor -------------------------
275         isize = 1
276         DO WHILE (isize .le. nsize(iaer,idomain))
277            READ(file_unit,*) scanline
278            IF ((scanline(1:1) .ne. '#').and.&
279            (scanline(1:1) .ne. ' ')) THEN
280            BACKSPACE(file_unit)
281            read(file_unit,*) gfactor(:,isize)
282            isize = isize + 1
283         ENDIF
284      ENDDO
285
286      jfile = jfile+1
287      IF ((idomain.NE.1).OR.(iaer.NE.1)) THEN
288         endwhile = .true.
289      ENDIF
290      CASE(6) reading2_seq
291         endwhile = .true.
292      CASE DEFAULT reading2_seq ! ----------------------------
293         WRITE(*,*) 'readoptprop: ',&
294         'Error while loading optical properties.'
295         CALL ABORT
296      END SELECT reading2_seq   ! ==============================
297      ENDIF
298      ENDDO
299
300!     1.4 Close the file
301
302      CLOSE(file_unit)
303
304!     1.5 If Francois' data, convert wvl to metres
305      if(forgetit)then
306         !   print*,'please sort out forgetit for naerkind>1'
307         if(iaer.eq.1)then
308            do iwvl=1,nwvl
309               wvl(iwvl)=wvl(iwvl)*1.e-6
310            enddo
311         endif
312      endif
313
314
315
316
317
318
319!==================================================================
320!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
321!==================================================================
322      domain: SELECT CASE (idomain)
323!==================================================================
324      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
325!==================================================================
326
327         lamref=lamrefvis(iaer)
328         epref=1.E+0
329
330         DO isize=1,nsize(iaer,idomain)
331
332!     Save the particle sizes
333            radiustab(iaer,idomain,isize)=radiusdyn(isize)
334
335!     Averaged optical properties (GCM channels)
336
337!            CALL aerave ( nwvl,&
338!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
339!            lamref,epref,tstellar,&
340!            L_NSPECTV,blamv,epavVI,&
341!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
342            CALL aerave_new ( nwvl,&
343            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
344            lamref,epref,tstellar,&
345            L_NSPECTV,blamv,epavVI,&
346            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
347
348!     Variable assignments (declared in radcommon)
349            DO isw=1,L_NSPECTV
350               QVISsQREF(isw,iaer,isize)=epavVI(isw)
351               gvis(isw,iaer,isize)=gavVI(isw)
352               omegavis(isw,iaer,isize)=omegavVI(isw)
353            END DO
354
355         ENDDO
356
357!==================================================================
358      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
359!==================================================================
360
361
362         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
363
364            lamref=lamrefir(iaer)
365            epref=1.E+0
366
367!     Save the particle sizes
368            radiustab(iaer,idomain,isize)=radiusdyn(isize)
369
370!     Averaged optical properties (GCM channels)
371
372!     epav is <QIR>/Qext(lamrefir) since epref=1
373!     Note: aerave also computes the extinction coefficient at
374!     the reference wavelength. This is called QREFvis or QREFir
375!     (not epref, which is a different parameter).
376!     Reference wavelengths SHOULD BE defined for each aerosol in
377!     radcommon_h.F90
378
379!            CALL aerave ( nwvl,&
380!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
381!            lamref,epref,tplanet,&
382!            L_NSPECTI,blami,epavIR,&
383!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
384            CALL aerave_new ( nwvl,&
385            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
386            lamref,epref,tplanet,&
387            L_NSPECTI,blami,epavIR,&
388            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
389
390
391!     Variable assignments (declared in radcommon)
392            DO ilw=1,L_NSPECTI
393               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
394               gir(ilw,iaer,isize)=gavIR(ilw)
395               omegair(ilw,iaer,isize)=omegavIR(ilw)
396            END DO
397
398
399      ENDDO                     ! isize (particle size) -------------------------------------
400
401      END SELECT domain
402
403!========================================================================
404!     3. Deallocate temporary variables that were read in the ASCII files
405!========================================================================
406
407      DEALLOCATE(wvl)           ! wvl
408      DEALLOCATE(radiusdyn)     ! radiusdyn
409      DEALLOCATE(ep)            ! ep
410      DEALLOCATE(omeg)          ! omeg
411      DEALLOCATE(gfactor)       ! g
412
413      END DO                    ! Loop on iaer
414      END DO                    ! Loop on idomain
415!========================================================================
416
417      RETURN
418
419
420
421    END subroutine suaer_corrk
422     
Note: See TracBrowser for help on using the repository browser.