source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/suaer_corrk.F90 @ 298

Last change on this file since 298 was 298, checked in by milmd, 10 years ago

Less output messages are written. On 20000 cores it is better. In LMDZ, only master of MPI and OpenMP can write.

File size: 15.6 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      use mod_phys_lmdz_para, only : is_master
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#include "callkeys.h"
48!     Optical properties (read in external ASCII files)
49      INTEGER,SAVE      :: nwvl  ! Number of wavelengths in
50                                ! the domain (VIS or IR), read by master
51
52!      REAL             :: solsir ! visible to infrared ratio
53!                                ! (tau_0.67um/tau_9um)
54
55      REAL, DIMENSION(:),&
56      ALLOCATABLE, SAVE :: wvl  ! Wavelength axis, read by master
57      REAL, DIMENSION(:),&
58      ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis, read by master
59
60      REAL, DIMENSION(:,:),&
61      ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext, read by master
62      omeg,&                    ! Single Scattering Albedo, read by master
63      gfactor                   ! Assymetry Factor, read by master
64
65!     Local variables:
66
67      INTEGER :: iaer           ! Aerosol index
68      INTEGER :: idomain        ! Domain index (1=VIS,2=IR)
69      INTEGER :: ilw            ! longwave index
70      INTEGER :: isw           ! shortwave index
71      INTEGER :: isize          ! Particle size index
72      INTEGER :: jfile          ! ASCII file scan index
73      INTEGER :: file_unit = 60
74      LOGICAL :: file_ok, endwhile
75      CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
76
77!     I/O  of "aerave" (subroutine that spectrally averages
78!     the single scattering parameters)
79
80      REAL lamref                      ! reference wavelengths
81      REAL epref                       ! reference extinction ep
82
83!      REAL epav(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
84!      REAL omegav(L_NSPECTI)          ! Average single scattering albedo
85!      REAL gav(L_NSPECTI)             ! Average assymetry parameter
86
87      REAL epavVI(L_NSPECTV)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
88      REAL omegavVI(L_NSPECTV)          ! Average single scattering albedo
89      REAL gavVI(L_NSPECTV)             ! Average assymetry parameter
90
91      REAL epavIR(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
92      REAL omegavIR(L_NSPECTI)          ! Average single scattering albedo
93      REAL gavIR(L_NSPECTI)             ! Average assymetry parameter
94     
95      logical forgetit                  ! use francois' data?
96      integer iwvl
97
98!     Local saved variables:
99
100      CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id
101!$OMP THREADPRIVATE(file_id)
102!---- Please indicate the names of the optical property files below
103!     Please also choose the reference wavelengths of each aerosol
104
105      if (noaero) then
106        print*, 'naerkind= 0'
107      endif
108      do iaer=1,naerkind
109       if (iaer.eq.iaero_co2) then
110        forgetit=.true.
111          if (.not.noaero) then
112              print*, 'naerkind= co2', iaer
113          end if
114!     visible
115        if(forgetit)then
116           file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
117        else
118           file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
119        endif
120        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
121
122!     infrared
123        if(forgetit)then
124           file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
125        else
126           file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
127        endif
128        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
129       endif ! CO2 aerosols 
130! NOTE: these lamref's are currently for dust, not CO2 ice.
131! JB thinks this shouldn't matter too much, but it needs to be
132! fixed / decided for the final version.
133
134       if (iaer.eq.iaero_h2o) then
135        print*, 'naerkind= h2o', iaer
136
137!     visible
138         file_id(iaer,1) = 'optprop_icevis_n50.dat'
139         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
140!     infrared
141         file_id(iaer,2) = 'optprop_iceir_n50.dat'
142         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
143       endif
144
145       if (iaer.eq.iaero_dust) then
146        print*, 'naerkind= dust', iaer
147
148!     visible
149         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
150         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
151         lamrefvis(iaer)=0.67e-6
152!     infrared
153         file_id(iaer,2) = 'optprop_dustir_n50.dat'
154         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
155         lamrefir(iaer)=9.3E-6
156       endif
157
158       if (iaer.eq.iaero_h2so4) then
159         print*, 'naerkind= h2so4', iaer
160
161!     visible
162         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
163         !lamrefir(iaer)=  doesn't exist?
164         lamrefvis(iaer)=1.5E-6   ! no idea, must find
165!     infrared
166         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
167         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
168         lamrefir(iaer)=9.3E-6 ! no idea, must find
169! added by LK
170       endif
171
172       if (iaer.eq.iaero_back2lay) then
173         if (is_master) print*, 'naerkind= back2lay', iaer
174
175!     visible
176         file_id(iaer,1) = 'optprop_saturn_vis_n20.dat'
177         lamrefvis(iaer)=0.8E-6  !
178!     infrared
179         file_id(iaer,2) = 'optprop_saturn_ir_n20.dat'
180         lamrefir(iaer)=6.E-6  !
181! added by SG
182       endif
183       
184       
185      enddo
186
187!------------------------------------------------------------------
188
189!     Initializations:
190
191      radiustab(:,:,:) = 0.0
192      gvis(:,:,:)      = 0.0
193      omegavis(:,:,:)  = 0.0
194      QVISsQREF(:,:,:) = 0.0
195      gir(:,:,:)       = 0.0
196      omegair(:,:,:)   = 0.0
197      QIRsQREF(:,:,:)  = 0.0
198
199      DO iaer = 1, naerkind     ! Loop on aerosol kind
200         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
201!==================================================================
202!     1. READ OPTICAL PROPERTIES
203!==================================================================
204
205!     1.1 Open the ASCII file
206
207!$OMP MASTER
208            INQUIRE(FILE=TRIM(datadir)//&
209         '/'//TRIM(file_id(iaer,idomain)),&
210            EXIST=file_ok)
211            IF(.NOT.file_ok) THEN
212               write(*,*)'suaer_corrk: Problem opening ',&
213               TRIM(file_id(iaer,idomain))
214               write(*,*)'It should be in: ',TRIM(datadir)
215               write(*,*)'1) You can set this directory address ',&
216               'in callphys.def with:'
217               write(*,*)' datadir = /absolute/path/to/datagcm'
218               write(*,*)'2) If ',&
219              TRIM(file_id(iaer,idomain)),&
220               ' is a LMD reference datafile, it'
221               write(*,*)' can be obtained online at:'
222               write(*,*)' http://www.lmd.jussieu.fr/',&
223               '~forget/datagcm/datafile'
224               CALL ABORT
225            ENDIF
226            OPEN(UNIT=file_unit,&
227            FILE=TRIM(datadir)//&
228         '/'//TRIM(file_id(iaer,idomain)),&
229            FORM='formatted')
230
231!     1.2 Allocate the optical property table
232
233            jfile = 1
234            endwhile = .false.
235            DO WHILE (.NOT.endwhile)
236               READ(file_unit,*) scanline
237               IF ((scanline(1:1) .ne. '#').and.&
238               (scanline(1:1) .ne. ' ')) THEN
239               BACKSPACE(file_unit)
240               reading1_seq: SELECT CASE (jfile) ! ====================
241            CASE(1) reading1_seq ! nwvl ----------------------------
242               read(file_unit,*) nwvl
243               jfile = jfile+1
244            CASE(2) reading1_seq ! nsize ---------------------------
245               read(file_unit,*) nsize(iaer,idomain)
246               endwhile = .true.
247            CASE DEFAULT reading1_seq ! ----------------------------
248               WRITE(*,*) 'readoptprop: ',& 
249               'Error while loading optical properties.' 
250               CALL ABORT
251            END SELECT reading1_seq ! ==============================
252         ENDIF
253      ENDDO
254
255      ALLOCATE(wvl(nwvl))       ! wvl
256      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
257      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
258      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
259      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
260
261
262!     1.3 Read the data
263
264      jfile = 1
265      endwhile = .false.
266      DO WHILE (.NOT.endwhile)
267         READ(file_unit,*) scanline
268         IF ((scanline(1:1) .ne. '#').and.&
269         (scanline(1:1) .ne. ' ')) THEN
270         BACKSPACE(file_unit)
271         reading2_seq: SELECT CASE (jfile) ! ====================
272      CASE(1) reading2_seq      ! wvl -----------------------------
273         read(file_unit,*) wvl
274         jfile = jfile+1
275      CASE(2) reading2_seq      ! radiusdyn -----------------------
276         read(file_unit,*) radiusdyn
277         jfile = jfile+1
278      CASE(3) reading2_seq      ! ep ------------------------------
279         isize = 1
280         DO WHILE (isize .le. nsize(iaer,idomain))
281            READ(file_unit,*) scanline
282            IF ((scanline(1:1) .ne. '#').and.&
283            (scanline(1:1) .ne. ' ')) THEN
284            BACKSPACE(file_unit)
285            read(file_unit,*) ep(:,isize)
286            isize = isize + 1
287         ENDIF
288      ENDDO
289
290      jfile = jfile+1
291      CASE(4) reading2_seq      ! omeg ----------------------------
292         isize = 1
293         DO WHILE (isize .le. nsize(iaer,idomain))
294            READ(file_unit,*) scanline
295            IF ((scanline(1:1) .ne. '#').and.&
296            (scanline(1:1) .ne. ' ')) THEN
297            BACKSPACE(file_unit)
298            read(file_unit,*) omeg(:,isize)
299            isize = isize + 1
300         ENDIF
301      ENDDO
302
303      jfile = jfile+1
304      CASE(5) reading2_seq      ! gfactor -------------------------
305         isize = 1
306         DO WHILE (isize .le. nsize(iaer,idomain))
307            READ(file_unit,*) scanline
308            IF ((scanline(1:1) .ne. '#').and.&
309            (scanline(1:1) .ne. ' ')) THEN
310            BACKSPACE(file_unit)
311            read(file_unit,*) gfactor(:,isize)
312            isize = isize + 1
313         ENDIF
314      ENDDO
315
316      jfile = jfile+1
317      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
318         endwhile = .true.
319      ENDIF
320      CASE(6) reading2_seq
321         endwhile = .true.
322      CASE DEFAULT reading2_seq ! ----------------------------
323         WRITE(*,*) 'readoptprop: ',&
324         'Error while loading optical properties.'
325         CALL ABORT
326      END SELECT reading2_seq   ! ==============================
327      ENDIF
328      ENDDO
329
330!     1.4 Close the file
331
332      CLOSE(file_unit)
333
334!     1.5 If Francois' data, convert wvl to metres
335       if(iaer.eq.iaero_co2)then
336         if(forgetit)then
337         !   print*,'please sort out forgetit for naerkind>1'
338            do iwvl=1,nwvl
339               wvl(iwvl)=wvl(iwvl)*1.e-6
340            enddo
341         endif
342      endif
343
344!$OMP END MASTER
345!$OMP BARRIER
346
347
348
349
350
351!==================================================================
352!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
353!==================================================================
354      domain: SELECT CASE (idomain)
355!==================================================================
356      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
357!==================================================================
358
359         lamref=lamrefvis(iaer)
360         epref=1.E+0
361
362         DO isize=1,nsize(iaer,idomain)
363
364!     Save the particle sizes
365            radiustab(iaer,idomain,isize)=radiusdyn(isize)
366
367!     Averaged optical properties (GCM channels)
368
369!            CALL aerave ( nwvl,&
370!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
371!            lamref,epref,tstellar,&
372!            L_NSPECTV,blamv,epavVI,&
373!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
374            CALL aerave_new ( nwvl,&
375            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
376            lamref,epref,tstellar,&
377            L_NSPECTV,blamv,epavVI,&
378            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
379
380!     Variable assignments (declared in radcommon)
381            DO isw=1,L_NSPECTV
382               QVISsQREF(isw,iaer,isize)=epavVI(isw)
383               gvis(isw,iaer,isize)=gavVI(isw)
384               omegavis(isw,iaer,isize)=omegavVI(isw)
385            END DO
386
387         ENDDO
388
389!==================================================================
390      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
391!==================================================================
392
393
394         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
395
396            lamref=lamrefir(iaer)
397            epref=1.E+0
398
399!     Save the particle sizes
400            radiustab(iaer,idomain,isize)=radiusdyn(isize)
401
402!     Averaged optical properties (GCM channels)
403
404!     epav is <QIR>/Qext(lamrefir) since epref=1
405!     Note: aerave also computes the extinction coefficient at
406!     the reference wavelength. This is called QREFvis or QREFir
407!     (not epref, which is a different parameter).
408!     Reference wavelengths SHOULD BE defined for each aerosol in
409!     radcommon_h.F90
410
411!            CALL aerave ( nwvl,&
412!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
413!            lamref,epref,tplanet,&
414!            L_NSPECTI,blami,epavIR,&
415!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
416            CALL aerave_new ( nwvl,&
417            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
418            lamref,epref,tplanet,&
419            L_NSPECTI,blami,epavIR,&
420            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
421
422
423!     Variable assignments (declared in radcommon)
424            DO ilw=1,L_NSPECTI
425               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
426               gir(ilw,iaer,isize)=gavIR(ilw)
427               omegair(ilw,iaer,isize)=omegavIR(ilw)
428            END DO
429
430
431      ENDDO                     ! isize (particle size) -------------------------------------
432
433      END SELECT domain
434
435!========================================================================
436!     3. Deallocate temporary variables that were read in the ASCII files
437!========================================================================
438
439!$OMP BARRIER
440!$OMP MASTER
441      IF (ALLOCATED(wvl)) DEALLOCATE(wvl)                 ! wvl
442      IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn)     ! radiusdyn
443      IF (ALLOCATED(ep)) DEALLOCATE(ep)                   ! ep
444      IF (ALLOCATED(omeg)) DEALLOCATE(omeg)               ! omeg
445      IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor)         ! g
446!$OMP END MASTER
447!$OMP BARRIER
448
449      END DO                    ! Loop on iaer
450      END DO                    ! Loop on idomain
451!========================================================================
452
453      RETURN
454
455
456
457    END subroutine suaer_corrk
458     
Note: See TracBrowser for help on using the repository browser.