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

Last change on this file since 247 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

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: overrides dust
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 fixed for the
124! 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      ENDIF
136
137      IF (naerkind .gt. 2) THEN
138         print*,'naerkind = ',naerkind
139         print*,'but we only have data for 2 types, exiting.'
140         call abort
141      ENDIF
142
143!------------------------------------------------------------------
144
145!     Initializations:
146
147!      call zerophys(nsizemax*2*naerkind,radiustab)
148!      call zerophys(nsizemax*naerkind*L_NSPECTV,gvis)
149!      call zerophys(nsizemax*naerkind*L_NSPECTV,omegavis)
150!      call zerophys(nsizemax*naerkind*L_NSPECTV,QVISsQREF)
151!      call zerophys(nsizemax*naerkind*L_NSPECTI,gIR)
152!      call zerophys(nsizemax*naerkind*L_NSPECTI,omegaIR)
153!      call zerophys(nsizemax*naerkind*L_NSPECTI,QIRsQREF)
154
155      radiustab(:,:,:) = 0.0
156      gvis(:,:,:)      = 0.0
157      omegavis(:,:,:)  = 0.0
158      QVISsQREF(:,:,:) = 0.0
159      gir(:,:,:)       = 0.0
160      omegair(:,:,:)   = 0.0
161      QIRsQREF(:,:,:)  = 0.0
162
163      DO iaer = 1, naerkind     ! Loop on aerosol kind
164         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
165!==================================================================
166!     1. READ OPTICAL PROPERTIES
167!==================================================================
168
169!     1.1 Open the ASCII file
170
171            INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//&
172         '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
173            EXIST=file_ok)
174            IF(.NOT.file_ok) THEN
175               write(*,*)'Problem opening ',&
176               file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))
177               write(*,*)'It should be in: ',&
178               datafile(1:LEN_TRIM(datafile))
179               write(*,*)'1) You can change this directory address in '
180               write(*,*)' file phymars/datafile.h'
181               write(*,*)'2) If ',&
182              file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
183               ' is a LMD reference datafile, it'
184               write(*,*)' can be obtained online at:'
185               write(*,*)' http://www.lmd.jussieu.fr/',&
186               '~forget/datagcm/datafile'
187               write(*,*)'3) If the name of the file is wrong, you can'
188               write(*,*)' change it in file phyuniv/suaer_corrk.F90. Just'
189               write(*,*)' modify the variable called file_id.'
190               CALL ABORT
191            ENDIF
192            OPEN(UNIT=file_unit,&
193            FILE=datafile(1:LEN_TRIM(datafile))//&
194         '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),&
195            FORM='formatted')
196
197!     1.2 Allocate the optical property table
198
199            jfile = 1
200            endwhile = .false.
201            DO WHILE (.NOT.endwhile)
202               READ(file_unit,*) scanline
203               IF ((scanline(1:1) .ne. '#').and.&
204               (scanline(1:1) .ne. ' ')) THEN
205               BACKSPACE(file_unit)
206               reading1_seq: SELECT CASE (jfile) ! ====================
207            CASE(1) reading1_seq ! nwvl ----------------------------
208               read(file_unit,*) nwvl
209               jfile = jfile+1
210            CASE(2) reading1_seq ! nsize ---------------------------
211               read(file_unit,*) nsize(iaer,idomain)
212               endwhile = .true.
213            CASE DEFAULT reading1_seq ! ----------------------------
214               WRITE(*,*) 'readoptprop: ',&
215               'Error while loading optical properties.'
216               CALL ABORT
217            END SELECT reading1_seq ! ==============================
218         ENDIF
219      ENDDO
220
221      ALLOCATE(wvl(nwvl))       ! wvl
222      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
223      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
224      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
225      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
226
227
228!     1.3 Read the data
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         reading2_seq: SELECT CASE (jfile) ! ====================
238      CASE(1) reading2_seq      ! wvl -----------------------------
239         read(file_unit,*) wvl
240         jfile = jfile+1
241      CASE(2) reading2_seq      ! radiusdyn -----------------------
242         read(file_unit,*) radiusdyn
243         jfile = jfile+1
244      CASE(3) reading2_seq      ! ep ------------------------------
245         isize = 1
246         DO WHILE (isize .le. nsize(iaer,idomain))
247            READ(file_unit,*) scanline
248            IF ((scanline(1:1) .ne. '#').and.&
249            (scanline(1:1) .ne. ' ')) THEN
250            BACKSPACE(file_unit)
251            read(file_unit,*) ep(:,isize)
252            isize = isize + 1
253         ENDIF
254      ENDDO
255
256      jfile = jfile+1
257      CASE(4) reading2_seq      ! omeg ----------------------------
258         isize = 1
259         DO WHILE (isize .le. nsize(iaer,idomain))
260            READ(file_unit,*) scanline
261            IF ((scanline(1:1) .ne. '#').and.&
262            (scanline(1:1) .ne. ' ')) THEN
263            BACKSPACE(file_unit)
264            read(file_unit,*) omeg(:,isize)
265            isize = isize + 1
266         ENDIF
267      ENDDO
268
269      jfile = jfile+1
270      CASE(5) reading2_seq      ! gfactor -------------------------
271         isize = 1
272         DO WHILE (isize .le. nsize(iaer,idomain))
273            READ(file_unit,*) scanline
274            IF ((scanline(1:1) .ne. '#').and.&
275            (scanline(1:1) .ne. ' ')) THEN
276            BACKSPACE(file_unit)
277            read(file_unit,*) gfactor(:,isize)
278            isize = isize + 1
279         ENDIF
280      ENDDO
281
282      jfile = jfile+1
283      IF ((idomain.NE.1).OR.(iaer.NE.1)) THEN
284         endwhile = .true.
285      ENDIF
286      CASE(6) reading2_seq
287         endwhile = .true.
288      CASE DEFAULT reading2_seq ! ----------------------------
289         WRITE(*,*) 'readoptprop: ',&
290         'Error while loading optical properties.'
291         CALL ABORT
292      END SELECT reading2_seq   ! ==============================
293      ENDIF
294      ENDDO
295
296!     1.4 Close the file
297
298      CLOSE(file_unit)
299
300!     1.5 If Francois' data, convert wvl to metres
301      if(forgetit)then
302         !   print*,'please sort out forgetit for naerkind>1'
303         if(iaer.eq.1)then
304            do iwvl=1,nwvl
305               wvl(iwvl)=wvl(iwvl)*1.e-6
306            enddo
307         endif
308      endif
309
310!==================================================================
311!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
312!==================================================================
313      domain: SELECT CASE (idomain)
314!==================================================================
315      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
316!==================================================================
317
318         lamref=lamrefvis(iaer)
319         epref=1.E+0
320
321         DO isize=1,nsize(iaer,idomain)
322
323!     Save the particle sizes
324            radiustab(iaer,idomain,isize)=radiusdyn(isize)
325
326!     Averaged optical properties (GCM channels)
327
328!            CALL aerave ( nwvl,&
329!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
330!            lamref,epref,tstellar,&
331!            L_NSPECTV,blamv,epavVI,&
332!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
333            CALL aerave_new ( nwvl,&
334            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
335            lamref,epref,tstellar,&
336            L_NSPECTV,blamv,epavVI,&
337            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
338
339!     Variable assignments (declared in radcommon)
340            DO isw=1,L_NSPECTV
341               QVISsQREF(isw,iaer,isize)=epavVI(isw)
342               gvis(isw,iaer,isize)=gavVI(isw)
343               omegavis(isw,iaer,isize)=omegavVI(isw)
344            END DO
345
346         ENDDO
347
348!==================================================================
349      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
350!==================================================================
351
352
353         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
354
355            lamref=lamrefir(iaer)
356            epref=1.E+0
357
358!     Save the particle sizes
359            radiustab(iaer,idomain,isize)=radiusdyn(isize)
360
361!     Averaged optical properties (GCM channels)
362
363!     epav is <QIR>/Qext(lamrefir) since epref=1
364!     Note: aerave also computes the extinction coefficient at
365!     the reference wavelength. This is called QREFvis or QREFir
366!     (not epref, which is a different parameter).
367!     Reference wavelengths SHOULD BE defined for each aerosol in
368!     radcommon_h.F90
369
370!            CALL aerave ( nwvl,&
371!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
372!            lamref,epref,tplanet,&
373!            L_NSPECTI,blami,epavIR,&
374!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
375            CALL aerave_new ( nwvl,&
376            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
377            lamref,epref,tplanet,&
378            L_NSPECTI,blami,epavIR,&
379            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
380
381
382!     Variable assignments (declared in radcommon)
383            DO ilw=1,L_NSPECTI
384               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
385               gir(ilw,iaer,isize)=gavIR(ilw)
386               omegair(ilw,iaer,isize)=omegavIR(ilw)
387            END DO
388
389
390      ENDDO                     ! isize (particle size) -------------------------------------
391
392      END SELECT domain
393
394
395!========================================================================
396!     3. Deallocate temporary variables that were read in the ASCII files
397!========================================================================
398
399      DEALLOCATE(wvl)           ! wvl
400      DEALLOCATE(radiusdyn)     ! radiusdyn
401      DEALLOCATE(ep)            ! ep
402      DEALLOCATE(omeg)          ! omeg
403      DEALLOCATE(gfactor)       ! g
404
405      END DO                    ! Loop on iaer
406      END DO                    ! Loop on idomain
407!========================================================================
408      RETURN
409    END subroutine suaer_corrk
410     
Note: See TracBrowser for help on using the repository browser.