source: trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90 @ 3585

Last change on this file since 3585 was 3585, checked in by debatzbr, 9 hours ago

Connecting microphysics to radiative transfer + miscellaneous cleans

File size: 17.0 KB
Line 
1module suaer_corrk_mod
2
3implicit none
4
5contains
6
7subroutine suaer_corrk
8
9      ! inputs
10      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
11      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
12      use datafile_mod, only: datadir, aerdir, &
13                              hazeprop_file, aersprop_file, aerfprop_file
14
15      ! outputs
16      use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
17      use radcommon_h, only: radiustab,nsize,tstellar
18      use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
19      use aerosol_mod, only: iaero_haze
20      use callkeys_mod, only: tplanet, callmufi
21      use tracer_h, only: noms
22
23      use mod_phys_lmdz_para, only : is_master, bcast
24
25      implicit none
26
27!==================================================================
28!     Purpose
29!     -------
30!     Initialize all radiative aerosol properties
31!
32!     Notes
33!     -----
34!     Reads the optical properties -> Mean  -> Variable assignment
35!     (ASCII files)                  (see radcommon_h.F90)
36!     wvl(nwvl)                      longsun
37!     ep(nwvl)                       epav     QVISsQREF(L_NSPECTV)
38!     omeg(nwvl)                     omegav   omegavis(L_NSPECTV)
39!     gfactor(nwvl)                  gav      gvis(L_NSPECTV)
40!
41!     Authors
42!     -------
43!     Richard Fournier (1996) Francois Forget (1996)
44!     Frederic Hourdin
45!     Jean-jacques morcrette *ECMWF*
46!     MODIF Francois Forget (2000)
47!     MODIF Franck Montmessin (add water ice)
48!     MODIF J.-B. Madeleine 2008W27
49!     - Optical properties read in ASCII files
50!     - Add varying radius of the particules
51!     - Add water ice clouds
52!     MODIF R. Wordsworth (2009)
53!     - generalisation to arbitrary spectral bands
54!
55!==================================================================
56
57!     Optical properties (read in external ASCII files)
58      INTEGER      :: nwvl  ! Number of wavelengths in
59                                ! the domain (VIS or IR), read by master
60
61!      REAL             :: solsir ! visible to infrared ratio
62!                                ! (tau_0.67um/tau_9um)
63
64      REAL, DIMENSION(:),&
65      ALLOCATABLE :: wvl  ! Wavelength axis, read by master
66      REAL, DIMENSION(:),&
67      ALLOCATABLE :: radiusdyn ! Particle size axis, read by master
68
69      REAL, DIMENSION(:,:),&
70      ALLOCATABLE :: ep,& ! Extinction coefficient Qext, read by master
71      omeg,&                    ! Single Scattering Albedo, read by master
72      gfactor                   ! Assymetry Factor, read by master
73
74!     Local variables:
75
76      INTEGER :: iaer           ! Aerosol index
77      INTEGER :: idomain        ! Domain index (1=VIS,2=IR)
78      INTEGER :: ilw            ! longwave index
79      INTEGER :: isw            ! shortwave index
80      INTEGER :: isize          ! Particle size index
81      INTEGER :: jfile          ! ASCII file scan index
82      INTEGER :: file_unit = 60
83      LOGICAL :: file_ok, endwhile
84      CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
85
86!     I/O  of "aerave" (subroutine that spectrally averages
87!     the single scattering parameters)
88
89      REAL lamref                      ! reference wavelengths
90      REAL epref                       ! reference extinction ep
91
92      REAL epavVI(L_NSPECTV)            ! Average ep (= <Qext>/Qext(lamrefvis) if epref=1)
93      REAL omegavVI(L_NSPECTV)          ! Average single scattering albedo
94      REAL gavVI(L_NSPECTV)             ! Average assymetry parameter
95
96      REAL epavIR(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamrefir) if epref=1)
97      REAL omegavIR(L_NSPECTI)          ! Average single scattering albedo
98      REAL gavIR(L_NSPECTI)             ! Average assymetry parameter
99
100      logical forgetit                  ! use francois' data?
101      integer iwvl, ia
102
103!     Local saved variables:
104
105      CHARACTER(LEN=50),ALLOCATABLE :: file_id(:,:)
106
107!---- Please indicate the names of the optical property files below
108!     Please also choose the reference wavelengths of each aerosol
109
110!--------- README TO UNDERSTAND WHAT FOLLOWS (JVO 20) -------
111!     The lamref variables comes from the Martian model
112!     where the visible one is the one used for computing
113!     and the IR one is used to output scaled opacity to
114!     match instrumental data ... This is done (at least
115!     for now) in the generic, so lamrefir is dummy*!
116
117!     The important one is the VISIBLE one as it will be used
118!     to rescale the values in callcork.F90 assuming 'aerosol' is
119!     at this visible reference wavelenght.
120
121!     *Actually if you change lamrefir here there is a
122!     slight sensitvity in the outputs because of some
123!     machine precision differences that amplifys and lead
124!     up to 10-6 differences in the radiative balance...
125!     This could be good to clean the code but would require
126!     a lot of modifs and to take care !
127
128!--------------------------------------------------------------
129      ! allocate file_id, as naerkind is a variable
130      allocate(file_id(naerkind,2))
131
132
133      if (callmufi) then
134         if (is_master) then
135            write(*,*)'Suaer spherical aerosols optical properties, using: ', &
136                       TRIM(aersprop_file)
137            write(*,*)'Suaer fractal aerosols optical properties, using: ', &
138                       TRIM(aerfprop_file)
139         endif
140         ! Visible
141         file_id(1,1)=TRIM(aersprop_file)
142         file_id(2,1)=TRIM(aerfprop_file)
143         ! Infrared
144         file_id(1,2)=file_id(1,1)
145         file_id(2,2)=file_id(2,1)
146
147         do iaer=1,naerkind
148            lamrefvis(iaer) = 0.607E-6 ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
149            lamrefir(iaer)  = 2.E-6    ! Reference wavelength for opacity IR (in the LEISA range)
150         enddo ! end iaer=1,naerkind
151         
152      else
153         do iaer=1,naerkind
154            if (iaer.eq.1) then
155               ! Only one table of optical properties :
156               if (is_master) write(*,*)'Suaer haze optical properties, using: ', &
157                                       TRIM(hazeprop_file)
158               ! Visible
159               file_id(iaer,1)=TRIM(hazeprop_file)
160               ! Infrared
161               file_id(iaer,2)=file_id(iaer,1)
162
163               lamrefvis(iaer) = 0.607E-6 ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
164               lamrefir(iaer)  = 2.E-6    ! Reference wavelength for opacity IR (in the LEISA range)
165            endif ! end iaer.eq.1
166
167            if (naerkind.gt.1) then
168               ! Visible
169               file_id(iaer,1)=file_id(1,1)
170               ! Infrared
171               file_id(iaer,2)=file_id(1,1)
172
173               lamrefvis(iaer) = lamrefvis(1) ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
174               lamrefir(iaer)  = lamrefir(1)  ! Reference wavelength for opacity IR (in the LEISA range)
175            endif ! end naerkind.gt.1
176         enddo ! end iaer=1,naerkind
177      endif ! end callmufi
178
179!------------------------------------------------------------------
180
181!     Initializations:
182
183      radiustab(:,:,:) = 0.0
184      gvis(:,:,:)      = 0.0
185      omegavis(:,:,:)  = 0.0
186      QVISsQREF(:,:,:) = 0.0
187      gir(:,:,:)       = 0.0
188      omegair(:,:,:)   = 0.0
189      QIRsQREF(:,:,:)  = 0.0
190
191  DO iaer = 1, naerkind     ! Loop on aerosol kind
192    DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
193!==================================================================
194!     1. READ OPTICAL PROPERTIES
195!==================================================================
196
197!     1.1 Open the ASCII file
198
199!!!!$OMP MASTER
200      if (is_master) then
201
202            INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
203                    '/'//TRIM(file_id(iaer,idomain)),&
204                    EXIST=file_ok)
205            IF (file_ok) THEN
206              OPEN(UNIT=file_unit,&
207                   FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
208                        '/'//TRIM(file_id(iaer,idomain)),&
209                   FORM='formatted')
210            ELSE
211             ! In ye old days these files were stored in datadir;
212             ! let's be retro-compatible
213              INQUIRE(FILE=TRIM(datadir)//&
214                      '/'//TRIM(file_id(iaer,idomain)),&
215                      EXIST=file_ok)
216              IF (file_ok) THEN
217                OPEN(UNIT=file_unit,&
218                     FILE=TRIM(datadir)//&
219                          '/'//TRIM(file_id(iaer,idomain)),&
220                     FORM='formatted')
221              ENDIF
222            ENDIF
223            IF(.NOT.file_ok) THEN
224               write(*,*)'suaer_corrk: Problem opening ',&
225               TRIM(file_id(iaer,idomain))
226               write(*,*)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir)
227               write(*,*)'1) You can set this directory address ',&
228               'in callphys.def with:'
229               write(*,*)' datadir = /absolute/path/to/datagcm'
230               write(*,*)'2) If ',&
231              TRIM(file_id(iaer,idomain)),&
232               ' is a LMD reference datafile, it'
233               write(*,*)' can be obtained online at:'
234               write(*,*)' http://www.lmd.jussieu.fr/',&
235               '~lmdz/planets/LMDZ.GENERIC/datagcm/'
236               call abort_physic("suaer_corrk",'Problem with optical properties file',1)
237            ENDIF
238
239!     1.2 Allocate the optical property table
240
241            jfile = 1
242            endwhile = .false.
243            DO WHILE (.NOT.endwhile)
244               READ(file_unit,*) scanline
245               IF ((scanline(1:1) .ne. '#').and.&
246               (scanline(1:1) .ne. ' ')) THEN
247               BACKSPACE(file_unit)
248               reading1_seq: SELECT CASE (jfile) ! ====================
249            CASE(1) reading1_seq ! nwvl ----------------------------
250               read(file_unit,*) nwvl
251               jfile = jfile+1
252            CASE(2) reading1_seq ! nsize ---------------------------
253               read(file_unit,*) nsize(iaer,idomain)
254               endwhile = .true.
255            CASE DEFAULT reading1_seq ! ----------------------------
256               WRITE(*,*) 'readoptprop: ',&
257               'Error while loading optical properties.'
258               CALL ABORT
259            END SELECT reading1_seq ! ==============================
260         ENDIF
261      ENDDO
262
263      endif ! of if (is_master)
264
265      ! broadcast nwvl and nsize to all cores
266      call bcast(nwvl)
267      call bcast(nsize)
268
269      ALLOCATE(wvl(nwvl))       ! wvl
270      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
271      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
272      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
273      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
274
275
276!     1.3 Read the data
277
278      if (is_master) then
279      jfile = 1
280      endwhile = .false.
281      DO WHILE (.NOT.endwhile)
282         READ(file_unit,*) scanline
283         IF ((scanline(1:1) .ne. '#').and.&
284         (scanline(1:1) .ne. ' ')) THEN
285         BACKSPACE(file_unit)
286         reading2_seq: SELECT CASE (jfile) ! ====================
287      CASE(1) reading2_seq      ! wvl -----------------------------
288         read(file_unit,*) wvl
289         jfile = jfile+1
290      CASE(2) reading2_seq      ! radiusdyn -----------------------
291         read(file_unit,*) radiusdyn
292         jfile = jfile+1
293      CASE(3) reading2_seq      ! ep ------------------------------
294         isize = 1
295         DO WHILE (isize .le. nsize(iaer,idomain))
296            READ(file_unit,*) scanline
297            IF ((scanline(1:1) .ne. '#').and.&
298            (scanline(1:1) .ne. ' ')) THEN
299            BACKSPACE(file_unit)
300            read(file_unit,*) ep(:,isize)
301            isize = isize + 1
302         ENDIF
303      ENDDO
304
305      jfile = jfile+1
306      CASE(4) reading2_seq      ! omeg ----------------------------
307         isize = 1
308         DO WHILE (isize .le. nsize(iaer,idomain))
309            READ(file_unit,*) scanline
310            IF ((scanline(1:1) .ne. '#').and.&
311            (scanline(1:1) .ne. ' ')) THEN
312            BACKSPACE(file_unit)
313            read(file_unit,*) omeg(:,isize)
314            isize = isize + 1
315         ENDIF
316      ENDDO
317
318      jfile = jfile+1
319      CASE(5) reading2_seq      ! gfactor -------------------------
320         isize = 1
321         DO WHILE (isize .le. nsize(iaer,idomain))
322            READ(file_unit,*) scanline
323            IF ((scanline(1:1) .ne. '#').and.&
324            (scanline(1:1) .ne. ' ')) THEN
325            BACKSPACE(file_unit)
326            read(file_unit,*) gfactor(:,isize)
327            isize = isize + 1
328         ENDIF
329      ENDDO
330
331      jfile = jfile+1
332      print *, idomain, " ", iaer, " ", iaero_haze
333      IF ((idomain.NE.1).OR.(iaer.NE.1).OR.(jfile.EQ.6)) THEN
334         ! AF24: what does it mean? :-O
335         endwhile = .true.
336      ENDIF
337      CASE(6) reading2_seq
338         endwhile = .true.
339      CASE DEFAULT reading2_seq ! ----------------------------
340         WRITE(*,*) 'readoptprop: ',&
341         'Error while loading optical properties.'
342         CALL ABORT
343      END SELECT reading2_seq   ! ==============================
344      ENDIF
345      ENDDO
346
347!     1.4 Close the file
348
349      CLOSE(file_unit)
350
351!     1.5 If Francois' data, convert wvl to metres
352      !  if(iaer.eq.iaero_haze)then
353      !    if(forgetit)then
354      !    !   print*,'please sort out forgetit for naerkind>1'
355      !       do iwvl=1,nwvl
356      !          wvl(iwvl)=wvl(iwvl)*1.e-6
357      !       enddo
358      !    endif
359      ! endif
360
361      endif ! of if (is_master)
362
363      ! broadcast arrays to all cores
364      call bcast(wvl)
365      call bcast(radiusdyn)
366      call bcast(ep)
367      call bcast(omeg)
368      call bcast(gfactor)
369
370!==================================================================
371!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
372!==================================================================
373      domain: SELECT CASE (idomain)
374!==================================================================
375      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
376!==================================================================
377
378         lamref=lamrefvis(iaer)
379         epref=1.E+0
380
381         DO isize=1,nsize(iaer,idomain)
382
383!     Save the particle sizes
384            radiustab(iaer,idomain,isize)=radiusdyn(isize)
385
386!     Averaged optical properties (GCM channels)
387
388!            CALL aerave ( nwvl,&
389!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
390!            lamref,epref,tstellar,&
391!            L_NSPECTV,blamv,epavVI,&
392!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
393            CALL aerave_new ( nwvl,&
394            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
395            lamref,epref,tstellar,&
396            L_NSPECTV,blamv,epavVI,&
397            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
398
399!     Variable assignments (declared in radcommon)
400            DO isw=1,L_NSPECTV
401               QVISsQREF(isw,iaer,isize)=epavVI(isw)
402               gvis(isw,iaer,isize)=gavVI(isw)
403               omegavis(isw,iaer,isize)=omegavVI(isw)
404            END DO
405
406         ENDDO
407!==================================================================
408      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
409!==================================================================
410
411         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
412
413            lamref=lamrefir(iaer)
414            epref=1.E+0
415
416!     Save the particle sizes
417            radiustab(iaer,idomain,isize)=radiusdyn(isize)
418
419!     Averaged optical properties (GCM channels)
420
421!     epav is <QIR>/Qext(lamrefir) since epref=1
422!     Note: aerave also computes the extinction coefficient at
423!     the reference wavelength. This is called QREFvis or QREFir
424!     (not epref, which is a different parameter).
425!     Reference wavelengths SHOULD BE defined for each aerosol in
426!     radcommon_h.F90
427
428!            CALL aerave ( nwvl,&
429!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
430!            lamref,epref,tplanet,&
431!            L_NSPECTI,blami,epavIR,&
432!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
433            CALL aerave_new ( nwvl,&
434            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
435            lamref,epref,tplanet,&
436            L_NSPECTI,blami,epavIR,&
437            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
438
439
440!     Variable assignments (declared in radcommon)
441            DO ilw=1,L_NSPECTI
442               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
443               gir(ilw,iaer,isize)=gavIR(ilw)
444               omegair(ilw,iaer,isize)=omegavIR(ilw)
445            END DO
446
447
448         ENDDO ! isize (particle size) -------------------------------------
449
450      END SELECT domain
451
452!========================================================================
453!     3. Deallocate temporary variables that were read in the ASCII files
454!========================================================================
455
456      DEALLOCATE(wvl)             ! wvl
457      DEALLOCATE(radiusdyn)       ! radiusdyn
458      DEALLOCATE(ep)              ! ep
459      DEALLOCATE(omeg)            ! omeg
460      DEALLOCATE(gfactor)         ! g
461
462    END DO                    ! Loop on iaer
463  END DO                    ! Loop on idomain
464!========================================================================
465
466  ! cleanup
467  deallocate(file_id)
468
469end subroutine suaer_corrk
470
471end module suaer_corrk_mod
Note: See TracBrowser for help on using the repository browser.