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

Last change on this file since 1365 was 1315, checked in by milmd, 11 years ago

LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.

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