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

Last change on this file since 2062 was 2062, checked in by aboissinot, 6 years ago

Add keys (optprop_back2lay_vis, optprop_back2lay_ir) to set aeroback2lay aerosol optprop ir and visible files names
and solve the name conflict between Saturn and Jupiter.
Default values are the previously hard coded values, respectively 'optprop_saturn_vis_n20.dat' and 'optprop_saturn_ir_n20'.

File size: 16.9 KB
Line 
1      subroutine suaer_corrk
2
3      ! inputs
4      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
5      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
6      use datafile_mod, only: datadir, aerdir
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,omegarefir !,omegarefvis
12      use aerosol_mod
13      use callkeys_mod, only: tplanet, optprop_back2lay_vis, optprop_back2lay_ir
14
15      implicit none
16
17!==================================================================
18!     Purpose
19!     -------
20!     Initialize all radiative aerosol properties
21!
22!     Notes
23!     -----
24!     Reads the optical properties -> Mean  -> Variable assignment
25!     (ASCII files)                  (see radcommon_h.F90)
26!     wvl(nwvl)                      longsun
27!     ep(nwvl)                       epav     QVISsQREF(L_NSPECTV)
28!     omeg(nwvl)                     omegav   omegavis(L_NSPECTV)
29!     gfactor(nwvl)                  gav      gvis(L_NSPECTV)
30!     
31!     Authors
32!     -------
33!     Richard Fournier (1996) Francois Forget (1996)
34!     Frederic Hourdin
35!     Jean-jacques morcrette *ECMWF*
36!     MODIF Francois Forget (2000)
37!     MODIF Franck Montmessin (add water ice)
38!     MODIF J.-B. Madeleine 2008W27
39!     - Optical properties read in ASCII files
40!     - Add varying radius of the particules
41!     - Add water ice clouds
42!     MODIF R. Wordsworth (2009)
43!     - generalisation to arbitrary spectral bands
44!
45!==================================================================
46
47!     Optical properties (read in external ASCII files)
48      INTEGER,SAVE      :: nwvl  ! Number of wavelengths in
49                                ! the domain (VIS or IR), read by master
50
51!      REAL             :: solsir ! visible to infrared ratio
52!                                ! (tau_0.67um/tau_9um)
53
54      REAL, DIMENSION(:),&
55      ALLOCATABLE, SAVE :: wvl  ! Wavelength axis, read by master
56      REAL, DIMENSION(:),&
57      ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis, read by master
58
59      REAL, DIMENSION(:,:),&
60      ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext, read by master
61      omeg,&                    ! Single Scattering Albedo, read by master
62      gfactor                   ! Assymetry Factor, read by master
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!$OMP THREADPRIVATE(file_id)
101!---- Please indicate the names of the optical property files below
102!     Please also choose the reference wavelengths of each aerosol     
103
104      if (noaero) then
105        print*, 'naerkind= 0'
106      endif
107      do iaer=1,naerkind
108       if (iaer.eq.iaero_co2) then
109        forgetit=.true.
110          if (.not.noaero) then
111              print*, 'naerkind= co2', iaer
112          end if
113!     visible
114        if(forgetit)then
115           file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
116        else
117           file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
118        endif
119        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
120
121!     infrared
122        if(forgetit)then
123           file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
124        else
125           file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
126        endif
127        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
128       endif ! CO2 aerosols 
129! NOTE: these lamref's are currently for dust, not CO2 ice.
130! JB thinks this shouldn't matter too much, but it needs to be
131! fixed / decided for the final version.
132
133       if (iaer.eq.iaero_h2o) then
134        print*, 'naerkind= h2o', iaer
135
136!     visible
137         file_id(iaer,1) = 'optprop_icevis_n50.dat'
138         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
139!     infrared
140         file_id(iaer,2) = 'optprop_iceir_n50.dat'
141         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
142       endif
143
144       if (iaer.eq.iaero_dust) then
145        print*, 'naerkind= dust', iaer
146
147!     visible
148         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
149         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
150         lamrefvis(iaer)=0.67e-6
151!     infrared
152         file_id(iaer,2) = 'optprop_dustir_n50.dat'
153         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
154         lamrefir(iaer)=9.3E-6
155       endif
156
157       if (iaer.eq.iaero_h2so4) then
158         print*, 'naerkind= h2so4', iaer
159
160!     visible
161         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
162         !lamrefir(iaer)=  doesn't exist?
163         lamrefvis(iaer)=1.5E-6   ! no idea, must find
164!     infrared
165         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
166         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
167         lamrefir(iaer)=9.3E-6 ! no idea, must find
168! added by LK
169       endif
170
171       if (iaer.eq.iaero_back2lay) then
172         print*, 'naerkind= back2lay', iaer
173         
174!     visible
175         file_id(iaer,1) = TRIM(optprop_back2lay_vis)
176         lamrefvis(iaer)=0.8E-6  !
177!     infrared
178         file_id(iaer,2) = TRIM(optprop_back2lay_ir)
179         lamrefir(iaer)=6.E-6  !
180! added by SG
181       endif
182     
183       if (iaer.eq.iaero_nh3) then
184         print*, 'naerkind= nh3', iaer
185
186!     visible
187         file_id(iaer,1) = 'optprop_nh3ice_vis.dat'
188         lamrefvis(iaer)=0.756E-6  !
189!     infrared
190         file_id(iaer,2) = 'optprop_nh3ice_ir.dat'
191         lamrefir(iaer)=6.E-6  !
192! added by SG
193       endif
194
195       if (iaer.eq.iaero_aurora) then
196         print*, 'naerkind= aurora', iaer
197
198!     visible
199         file_id(iaer,1) = 'optprop_aurora_vis.dat'
200         lamrefvis(iaer)=0.3E-6  !
201!     infrared
202         file_id(iaer,2) = 'optprop_aurora_ir.dat'
203         lamrefir(iaer)=6.E-6  !
204! added by SG
205       endif
206 
207       
208      enddo
209
210!------------------------------------------------------------------
211
212!     Initializations:
213
214      radiustab(:,:,:) = 0.0
215      gvis(:,:,:)      = 0.0
216      omegavis(:,:,:)  = 0.0
217      QVISsQREF(:,:,:) = 0.0
218      gir(:,:,:)       = 0.0
219      omegair(:,:,:)   = 0.0
220      QIRsQREF(:,:,:)  = 0.0
221
222      DO iaer = 1, naerkind     ! Loop on aerosol kind
223         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
224!==================================================================
225!     1. READ OPTICAL PROPERTIES
226!==================================================================
227
228!     1.1 Open the ASCII file
229
230!$OMP MASTER
231
232            INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
233                    '/'//TRIM(file_id(iaer,idomain)),&
234                    EXIST=file_ok)
235            IF (file_ok) THEN
236              OPEN(UNIT=file_unit,&
237                   FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
238                        '/'//TRIM(file_id(iaer,idomain)),&
239                   FORM='formatted')
240            ELSE
241             ! In ye old days these files were stored in datadir;
242             ! let's be retro-compatible
243              INQUIRE(FILE=TRIM(datadir)//&
244                      '/'//TRIM(file_id(iaer,idomain)),&
245                      EXIST=file_ok)
246              IF (file_ok) THEN
247                OPEN(UNIT=file_unit,&
248                     FILE=TRIM(datadir)//&
249                          '/'//TRIM(file_id(iaer,idomain)),&
250                     FORM='formatted')
251              ENDIF             
252            ENDIF
253            IF(.NOT.file_ok) THEN
254               write(*,*)'suaer_corrk: Problem opening ',&
255               TRIM(file_id(iaer,idomain))
256               write(*,*)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir)
257               write(*,*)'1) You can set this directory address ',&
258               'in callphys.def with:'
259               write(*,*)' datadir = /absolute/path/to/datagcm'
260               write(*,*)'2) If ',&
261              TRIM(file_id(iaer,idomain)),&
262               ' is a LMD reference datafile, it'
263               write(*,*)' can be obtained online at:'
264               write(*,*)' http://www.lmd.jussieu.fr/',&
265               '~lmdz/planets/LMDZ.GENERIC/datagcm/'
266               CALL ABORT
267            ENDIF
268
269!     1.2 Allocate the optical property table
270
271            jfile = 1
272            endwhile = .false.
273            DO WHILE (.NOT.endwhile)
274               READ(file_unit,*) scanline
275               IF ((scanline(1:1) .ne. '#').and.&
276               (scanline(1:1) .ne. ' ')) THEN
277               BACKSPACE(file_unit)
278               reading1_seq: SELECT CASE (jfile) ! ====================
279            CASE(1) reading1_seq ! nwvl ----------------------------
280               read(file_unit,*) nwvl
281               jfile = jfile+1
282            CASE(2) reading1_seq ! nsize ---------------------------
283               read(file_unit,*) nsize(iaer,idomain)
284               endwhile = .true.
285            CASE DEFAULT reading1_seq ! ----------------------------
286               WRITE(*,*) 'readoptprop: ',&
287               'Error while loading optical properties.'
288               CALL ABORT
289            END SELECT reading1_seq ! ==============================
290         ENDIF
291      ENDDO
292
293      ALLOCATE(wvl(nwvl))       ! wvl
294      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
295      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
296      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
297      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
298
299
300!     1.3 Read the data
301
302      jfile = 1
303      endwhile = .false.
304      DO WHILE (.NOT.endwhile)
305         READ(file_unit,*) scanline
306         IF ((scanline(1:1) .ne. '#').and.&
307         (scanline(1:1) .ne. ' ')) THEN
308         BACKSPACE(file_unit)
309         reading2_seq: SELECT CASE (jfile) ! ====================
310      CASE(1) reading2_seq      ! wvl -----------------------------
311         read(file_unit,*) wvl
312         jfile = jfile+1
313      CASE(2) reading2_seq      ! radiusdyn -----------------------
314         read(file_unit,*) radiusdyn
315         jfile = jfile+1
316      CASE(3) reading2_seq      ! ep ------------------------------
317         isize = 1
318         DO WHILE (isize .le. nsize(iaer,idomain))
319            READ(file_unit,*) scanline
320            IF ((scanline(1:1) .ne. '#').and.&
321            (scanline(1:1) .ne. ' ')) THEN
322            BACKSPACE(file_unit)
323            read(file_unit,*) ep(:,isize)
324            isize = isize + 1
325         ENDIF
326      ENDDO
327
328      jfile = jfile+1
329      CASE(4) reading2_seq      ! omeg ----------------------------
330         isize = 1
331         DO WHILE (isize .le. nsize(iaer,idomain))
332            READ(file_unit,*) scanline
333            IF ((scanline(1:1) .ne. '#').and.&
334            (scanline(1:1) .ne. ' ')) THEN
335            BACKSPACE(file_unit)
336            read(file_unit,*) omeg(:,isize)
337            isize = isize + 1
338         ENDIF
339      ENDDO
340
341      jfile = jfile+1
342      CASE(5) reading2_seq      ! gfactor -------------------------
343         isize = 1
344         DO WHILE (isize .le. nsize(iaer,idomain))
345            READ(file_unit,*) scanline
346            IF ((scanline(1:1) .ne. '#').and.&
347            (scanline(1:1) .ne. ' ')) THEN
348            BACKSPACE(file_unit)
349            read(file_unit,*) gfactor(:,isize)
350            isize = isize + 1
351         ENDIF
352      ENDDO
353
354      jfile = jfile+1
355      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
356         endwhile = .true.
357      ENDIF
358      CASE(6) reading2_seq
359         endwhile = .true.
360      CASE DEFAULT reading2_seq ! ----------------------------
361         WRITE(*,*) 'readoptprop: ',&
362         'Error while loading optical properties.'
363         CALL ABORT
364      END SELECT reading2_seq   ! ==============================
365      ENDIF
366      ENDDO
367
368!     1.4 Close the file
369
370      CLOSE(file_unit)
371
372!     1.5 If Francois' data, convert wvl to metres
373       if(iaer.eq.iaero_co2)then
374         if(forgetit)then
375         !   print*,'please sort out forgetit for naerkind>1'
376            do iwvl=1,nwvl
377               wvl(iwvl)=wvl(iwvl)*1.e-6
378            enddo
379         endif
380      endif
381
382!$OMP END MASTER
383!$OMP BARRIER
384
385
386
387
388
389!==================================================================
390!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
391!==================================================================
392      domain: SELECT CASE (idomain)
393!==================================================================
394      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
395!==================================================================
396
397         lamref=lamrefvis(iaer)
398         epref=1.E+0
399
400         DO isize=1,nsize(iaer,idomain)
401
402!     Save the particle sizes
403            radiustab(iaer,idomain,isize)=radiusdyn(isize)
404
405!     Averaged optical properties (GCM channels)
406
407!            CALL aerave ( nwvl,&
408!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
409!            lamref,epref,tstellar,&
410!            L_NSPECTV,blamv,epavVI,&
411!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
412            CALL aerave_new ( nwvl,&
413            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
414            lamref,epref,tstellar,&
415            L_NSPECTV,blamv,epavVI,&
416            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
417
418!     Variable assignments (declared in radcommon)
419            DO isw=1,L_NSPECTV
420               QVISsQREF(isw,iaer,isize)=epavVI(isw)
421               gvis(isw,iaer,isize)=gavVI(isw)
422               omegavis(isw,iaer,isize)=omegavVI(isw)
423            END DO
424
425         ENDDO
426
427!==================================================================
428      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
429!==================================================================
430
431
432         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
433
434            lamref=lamrefir(iaer)
435            epref=1.E+0
436
437!     Save the particle sizes
438            radiustab(iaer,idomain,isize)=radiusdyn(isize)
439
440!     Averaged optical properties (GCM channels)
441
442!     epav is <QIR>/Qext(lamrefir) since epref=1
443!     Note: aerave also computes the extinction coefficient at
444!     the reference wavelength. This is called QREFvis or QREFir
445!     (not epref, which is a different parameter).
446!     Reference wavelengths SHOULD BE defined for each aerosol in
447!     radcommon_h.F90
448
449!            CALL aerave ( nwvl,&
450!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
451!            lamref,epref,tplanet,&
452!            L_NSPECTI,blami,epavIR,&
453!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
454            CALL aerave_new ( nwvl,&
455            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
456            lamref,epref,tplanet,&
457            L_NSPECTI,blami,epavIR,&
458            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
459
460
461!     Variable assignments (declared in radcommon)
462            DO ilw=1,L_NSPECTI
463               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
464               gir(ilw,iaer,isize)=gavIR(ilw)
465               omegair(ilw,iaer,isize)=omegavIR(ilw)
466            END DO
467
468
469      ENDDO                     ! isize (particle size) -------------------------------------
470
471      END SELECT domain
472
473!========================================================================
474!     3. Deallocate temporary variables that were read in the ASCII files
475!========================================================================
476
477!$OMP BARRIER
478!$OMP MASTER
479      IF (ALLOCATED(wvl)) DEALLOCATE(wvl)                 ! wvl
480      IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn)     ! radiusdyn
481      IF (ALLOCATED(ep)) DEALLOCATE(ep)                   ! ep
482      IF (ALLOCATED(omeg)) DEALLOCATE(omeg)               ! omeg
483      IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor)         ! g
484!$OMP END MASTER
485!$OMP BARRIER
486
487      END DO                    ! Loop on iaer
488      END DO                    ! Loop on idomain
489!========================================================================
490
491      RETURN
492
493
494
495    END subroutine suaer_corrk
496     
Note: See TracBrowser for help on using the repository browser.