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

Last change on this file since 3572 was 3572, checked in by debatzbr, 3 weeks ago

Remove generic_aerosols and generic_condensation, along with their related variables (useless). RENAME THE VARIABLE AEROHAZE TO OPTICHAZE.

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