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

Last change on this file since 2033 was 2033, checked in by emillour, 6 years ago

Generic GCM:

  • Cleanup around aeroptproperties.90, remove arrays that were not initialized but still used, probably for nothing since removing them leads to no change in results. Possibly to revisit and further clean up later.

EM

File size: 16.8 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
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) = 'optprop_saturn_vis_n20.dat'
176         lamrefvis(iaer)=0.8E-6  !
177!     infrared
178         file_id(iaer,2) = 'optprop_saturn_ir_n20.dat'
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
196       if (iaer.eq.iaero_aurora) then
197         print*, 'naerkind= aurora', iaer
198
199!     visible
200         file_id(iaer,1) = 'optprop_aurora_vis.dat'
201         lamrefvis(iaer)=0.3E-6  !
202!     infrared
203         file_id(iaer,2) = 'optprop_aurora_ir.dat'
204         lamrefir(iaer)=6.E-6  !
205! added by SG
206       endif
207 
208       
209      enddo
210
211!------------------------------------------------------------------
212
213!     Initializations:
214
215      radiustab(:,:,:) = 0.0
216      gvis(:,:,:)      = 0.0
217      omegavis(:,:,:)  = 0.0
218      QVISsQREF(:,:,:) = 0.0
219      gir(:,:,:)       = 0.0
220      omegair(:,:,:)   = 0.0
221      QIRsQREF(:,:,:)  = 0.0
222
223      DO iaer = 1, naerkind     ! Loop on aerosol kind
224         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
225!==================================================================
226!     1. READ OPTICAL PROPERTIES
227!==================================================================
228
229!     1.1 Open the ASCII file
230
231!$OMP MASTER
232
233            INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
234                    '/'//TRIM(file_id(iaer,idomain)),&
235                    EXIST=file_ok)
236            IF (file_ok) THEN
237              OPEN(UNIT=file_unit,&
238                   FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
239                        '/'//TRIM(file_id(iaer,idomain)),&
240                   FORM='formatted')
241            ELSE
242             ! In ye old days these files were stored in datadir;
243             ! let's be retro-compatible
244              INQUIRE(FILE=TRIM(datadir)//&
245                      '/'//TRIM(file_id(iaer,idomain)),&
246                      EXIST=file_ok)
247              IF (file_ok) THEN
248                OPEN(UNIT=file_unit,&
249                     FILE=TRIM(datadir)//&
250                          '/'//TRIM(file_id(iaer,idomain)),&
251                     FORM='formatted')
252              ENDIF             
253            ENDIF
254            IF(.NOT.file_ok) THEN
255               write(*,*)'suaer_corrk: Problem opening ',&
256               TRIM(file_id(iaer,idomain))
257               write(*,*)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir)
258               write(*,*)'1) You can set this directory address ',&
259               'in callphys.def with:'
260               write(*,*)' datadir = /absolute/path/to/datagcm'
261               write(*,*)'2) If ',&
262              TRIM(file_id(iaer,idomain)),&
263               ' is a LMD reference datafile, it'
264               write(*,*)' can be obtained online at:'
265               write(*,*)' http://www.lmd.jussieu.fr/',&
266               '~lmdz/planets/LMDZ.GENERIC/datagcm/'
267               CALL ABORT
268            ENDIF
269
270!     1.2 Allocate the optical property table
271
272            jfile = 1
273            endwhile = .false.
274            DO WHILE (.NOT.endwhile)
275               READ(file_unit,*) scanline
276               IF ((scanline(1:1) .ne. '#').and.&
277               (scanline(1:1) .ne. ' ')) THEN
278               BACKSPACE(file_unit)
279               reading1_seq: SELECT CASE (jfile) ! ====================
280            CASE(1) reading1_seq ! nwvl ----------------------------
281               read(file_unit,*) nwvl
282               jfile = jfile+1
283            CASE(2) reading1_seq ! nsize ---------------------------
284               read(file_unit,*) nsize(iaer,idomain)
285               endwhile = .true.
286            CASE DEFAULT reading1_seq ! ----------------------------
287               WRITE(*,*) 'readoptprop: ',&
288               'Error while loading optical properties.'
289               CALL ABORT
290            END SELECT reading1_seq ! ==============================
291         ENDIF
292      ENDDO
293
294      ALLOCATE(wvl(nwvl))       ! wvl
295      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
296      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
297      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
298      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
299
300
301!     1.3 Read the data
302
303      jfile = 1
304      endwhile = .false.
305      DO WHILE (.NOT.endwhile)
306         READ(file_unit,*) scanline
307         IF ((scanline(1:1) .ne. '#').and.&
308         (scanline(1:1) .ne. ' ')) THEN
309         BACKSPACE(file_unit)
310         reading2_seq: SELECT CASE (jfile) ! ====================
311      CASE(1) reading2_seq      ! wvl -----------------------------
312         read(file_unit,*) wvl
313         jfile = jfile+1
314      CASE(2) reading2_seq      ! radiusdyn -----------------------
315         read(file_unit,*) radiusdyn
316         jfile = jfile+1
317      CASE(3) reading2_seq      ! ep ------------------------------
318         isize = 1
319         DO WHILE (isize .le. nsize(iaer,idomain))
320            READ(file_unit,*) scanline
321            IF ((scanline(1:1) .ne. '#').and.&
322            (scanline(1:1) .ne. ' ')) THEN
323            BACKSPACE(file_unit)
324            read(file_unit,*) ep(:,isize)
325            isize = isize + 1
326         ENDIF
327      ENDDO
328
329      jfile = jfile+1
330      CASE(4) reading2_seq      ! omeg ----------------------------
331         isize = 1
332         DO WHILE (isize .le. nsize(iaer,idomain))
333            READ(file_unit,*) scanline
334            IF ((scanline(1:1) .ne. '#').and.&
335            (scanline(1:1) .ne. ' ')) THEN
336            BACKSPACE(file_unit)
337            read(file_unit,*) omeg(:,isize)
338            isize = isize + 1
339         ENDIF
340      ENDDO
341
342      jfile = jfile+1
343      CASE(5) reading2_seq      ! gfactor -------------------------
344         isize = 1
345         DO WHILE (isize .le. nsize(iaer,idomain))
346            READ(file_unit,*) scanline
347            IF ((scanline(1:1) .ne. '#').and.&
348            (scanline(1:1) .ne. ' ')) THEN
349            BACKSPACE(file_unit)
350            read(file_unit,*) gfactor(:,isize)
351            isize = isize + 1
352         ENDIF
353      ENDDO
354
355      jfile = jfile+1
356      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
357         endwhile = .true.
358      ENDIF
359      CASE(6) reading2_seq
360         endwhile = .true.
361      CASE DEFAULT reading2_seq ! ----------------------------
362         WRITE(*,*) 'readoptprop: ',&
363         'Error while loading optical properties.'
364         CALL ABORT
365      END SELECT reading2_seq   ! ==============================
366      ENDIF
367      ENDDO
368
369!     1.4 Close the file
370
371      CLOSE(file_unit)
372
373!     1.5 If Francois' data, convert wvl to metres
374       if(iaer.eq.iaero_co2)then
375         if(forgetit)then
376         !   print*,'please sort out forgetit for naerkind>1'
377            do iwvl=1,nwvl
378               wvl(iwvl)=wvl(iwvl)*1.e-6
379            enddo
380         endif
381      endif
382
383!$OMP END MASTER
384!$OMP BARRIER
385
386
387
388
389
390!==================================================================
391!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
392!==================================================================
393      domain: SELECT CASE (idomain)
394!==================================================================
395      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
396!==================================================================
397
398         lamref=lamrefvis(iaer)
399         epref=1.E+0
400
401         DO isize=1,nsize(iaer,idomain)
402
403!     Save the particle sizes
404            radiustab(iaer,idomain,isize)=radiusdyn(isize)
405
406!     Averaged optical properties (GCM channels)
407
408!            CALL aerave ( nwvl,&
409!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
410!            lamref,epref,tstellar,&
411!            L_NSPECTV,blamv,epavVI,&
412!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
413            CALL aerave_new ( nwvl,&
414            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
415            lamref,epref,tstellar,&
416            L_NSPECTV,blamv,epavVI,&
417            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
418
419!     Variable assignments (declared in radcommon)
420            DO isw=1,L_NSPECTV
421               QVISsQREF(isw,iaer,isize)=epavVI(isw)
422               gvis(isw,iaer,isize)=gavVI(isw)
423               omegavis(isw,iaer,isize)=omegavVI(isw)
424            END DO
425
426         ENDDO
427
428!==================================================================
429      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
430!==================================================================
431
432
433         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
434
435            lamref=lamrefir(iaer)
436            epref=1.E+0
437
438!     Save the particle sizes
439            radiustab(iaer,idomain,isize)=radiusdyn(isize)
440
441!     Averaged optical properties (GCM channels)
442
443!     epav is <QIR>/Qext(lamrefir) since epref=1
444!     Note: aerave also computes the extinction coefficient at
445!     the reference wavelength. This is called QREFvis or QREFir
446!     (not epref, which is a different parameter).
447!     Reference wavelengths SHOULD BE defined for each aerosol in
448!     radcommon_h.F90
449
450!            CALL aerave ( nwvl,&
451!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
452!            lamref,epref,tplanet,&
453!            L_NSPECTI,blami,epavIR,&
454!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
455            CALL aerave_new ( nwvl,&
456            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
457            lamref,epref,tplanet,&
458            L_NSPECTI,blami,epavIR,&
459            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
460
461
462!     Variable assignments (declared in radcommon)
463            DO ilw=1,L_NSPECTI
464               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
465               gir(ilw,iaer,isize)=gavIR(ilw)
466               omegair(ilw,iaer,isize)=omegavIR(ilw)
467            END DO
468
469
470      ENDDO                     ! isize (particle size) -------------------------------------
471
472      END SELECT domain
473
474!========================================================================
475!     3. Deallocate temporary variables that were read in the ASCII files
476!========================================================================
477
478!$OMP BARRIER
479!$OMP MASTER
480      IF (ALLOCATED(wvl)) DEALLOCATE(wvl)                 ! wvl
481      IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn)     ! radiusdyn
482      IF (ALLOCATED(ep)) DEALLOCATE(ep)                   ! ep
483      IF (ALLOCATED(omeg)) DEALLOCATE(omeg)               ! omeg
484      IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor)         ! g
485!$OMP END MASTER
486!$OMP BARRIER
487
488      END DO                    ! Loop on iaer
489      END DO                    ! Loop on idomain
490!========================================================================
491
492      RETURN
493
494
495
496    END subroutine suaer_corrk
497     
Note: See TracBrowser for help on using the repository browser.