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

Last change on this file since 1057 was 313, checked in by ymipsl, 10 years ago
  • implement splitting of XIOS file for lmdz physics
  • Termination is done properly in parallel by calling MPI_ABORT instead of abort or stop

YM

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_physiq 
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_physiq 
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_physiq
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.