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

Last change on this file was 3959, checked in by debatzbr, 2 months ago

Pluto PCM: Take clouds into account in radiative transfer.
BBT

File size: 17.3 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, clddprop_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, callmuclouds
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 (VIS & IR)
130      allocate(file_id(naerkind,2))
131
132      if (callmufi) then
133         if (is_master) then
134            write(*,*)'Suaer spherical aerosols optical properties, using: ', &
135                       TRIM(aersprop_file)
136            write(*,*)'Suaer fractal aerosols optical properties, using: ', &
137                       TRIM(aerfprop_file)
138            if (callmuclouds) then
139               write(*,*)'Suaer cloud drop optical properties, using: ', &
140                          TRIM(clddprop_file)
141            endif
142         endif
143         ! Visible
144         file_id(1,1)=TRIM(aersprop_file)
145         file_id(2,1)=TRIM(aerfprop_file)
146         if (callmuclouds) then
147            file_id(3,1)=TRIM(clddprop_file)
148         endif
149         ! Infrared
150         file_id(1,2)=file_id(1,1)
151         file_id(2,2)=file_id(2,1)
152         if (callmuclouds) then
153            file_id(3,2)=file_id(3,1)
154         endif
155
156         do iaer=1,naerkind
157            lamrefvis(iaer) = 0.607E-6 ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
158            lamrefir(iaer)  = 2.E-6    ! Reference wavelength for opacity IR (in the LEISA range)
159         enddo ! end iaer=1,naerkind
160
161      else
162         do iaer=1,naerkind
163            if (iaer.eq.1) then
164               ! Only one table of optical properties :
165               if (is_master) write(*,*)'Suaer haze optical properties, using: ', &
166                                       TRIM(hazeprop_file)
167               ! Visible
168               file_id(iaer,1)=TRIM(hazeprop_file)
169               ! Infrared
170               file_id(iaer,2)=file_id(iaer,1)
171
172               lamrefvis(iaer) = 0.607E-6 ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
173               lamrefir(iaer)  = 2.E-6    ! Reference wavelength for opacity IR (in the LEISA range)
174            endif ! end iaer.eq.1
175
176            if (naerkind.gt.1) then
177               ! Visible
178               file_id(iaer,1)=file_id(1,1)
179               ! Infrared
180               file_id(iaer,2)=file_id(1,1)
181
182               lamrefvis(iaer) = lamrefvis(1) ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
183               lamrefir(iaer)  = lamrefir(1)  ! Reference wavelength for opacity IR (in the LEISA range)
184            endif ! end naerkind.gt.1
185         enddo ! end iaer=1,naerkind
186      endif ! end callmufi
187
188!------------------------------------------------------------------
189
190!     Initializations:
191
192      radiustab(:,:,:) = 0.0
193      gvis(:,:,:)      = 0.0
194      omegavis(:,:,:)  = 0.0
195      QVISsQREF(:,:,:) = 0.0
196      gir(:,:,:)       = 0.0
197      omegair(:,:,:)   = 0.0
198      QIRsQREF(:,:,:)  = 0.0
199
200  DO iaer = 1, naerkind     ! Loop on aerosol kind
201    DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
202!==================================================================
203!     1. READ OPTICAL PROPERTIES
204!==================================================================
205
206!     1.1 Open the ASCII file
207
208!!!!$OMP MASTER
209      if (is_master) then
210
211            INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
212                    '/'//TRIM(file_id(iaer,idomain)),&
213                    EXIST=file_ok)
214            IF (file_ok) THEN
215              OPEN(UNIT=file_unit,&
216                   FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
217                        '/'//TRIM(file_id(iaer,idomain)),&
218                   FORM='formatted')
219            ELSE
220             ! In ye old days these files were stored in datadir;
221             ! let's be retro-compatible
222              INQUIRE(FILE=TRIM(datadir)//&
223                      '/'//TRIM(file_id(iaer,idomain)),&
224                      EXIST=file_ok)
225              IF (file_ok) THEN
226                OPEN(UNIT=file_unit,&
227                     FILE=TRIM(datadir)//&
228                          '/'//TRIM(file_id(iaer,idomain)),&
229                     FORM='formatted')
230              ENDIF
231            ENDIF
232            IF(.NOT.file_ok) THEN
233               write(*,*)'suaer_corrk: Problem opening ',&
234               TRIM(file_id(iaer,idomain))
235               write(*,*)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir)
236               write(*,*)'1) You can set this directory address ',&
237               'in callphys.def with:'
238               write(*,*)' datadir = /absolute/path/to/datagcm'
239               write(*,*)'2) If ',&
240              TRIM(file_id(iaer,idomain)),&
241               ' is a LMD reference datafile, it'
242               write(*,*)' can be obtained online at:'
243               write(*,*)' http://www.lmd.jussieu.fr/',&
244               '~lmdz/planets/LMDZ.GENERIC/datagcm/'
245               call abort_physic("suaer_corrk",'Problem with optical properties file',1)
246            ENDIF
247
248!     1.2 Allocate the optical property table
249
250            jfile = 1
251            endwhile = .false.
252            DO WHILE (.NOT.endwhile)
253               READ(file_unit,*) scanline
254               IF ((scanline(1:1) .ne. '#').and.&
255               (scanline(1:1) .ne. ' ')) THEN
256               BACKSPACE(file_unit)
257               reading1_seq: SELECT CASE (jfile) ! ====================
258            CASE(1) reading1_seq ! nwvl ----------------------------
259               read(file_unit,*) nwvl
260               jfile = jfile+1
261            CASE(2) reading1_seq ! nsize ---------------------------
262               read(file_unit,*) nsize(iaer,idomain)
263               endwhile = .true.
264            CASE DEFAULT reading1_seq ! ----------------------------
265               WRITE(*,*) 'readoptprop: ',&
266               'Error while loading optical properties.'
267               CALL ABORT
268            END SELECT reading1_seq ! ==============================
269         ENDIF
270      ENDDO
271
272      endif ! of if (is_master)
273
274      ! broadcast nwvl and nsize to all cores
275      call bcast(nwvl)
276      call bcast(nsize)
277
278      ALLOCATE(wvl(nwvl))       ! wvl
279      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
280      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
281      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
282      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
283
284
285!     1.3 Read the data
286
287      if (is_master) then
288      jfile = 1
289      endwhile = .false.
290      DO WHILE (.NOT.endwhile)
291         READ(file_unit,*) scanline
292         IF ((scanline(1:1) .ne. '#').and.&
293         (scanline(1:1) .ne. ' ')) THEN
294         BACKSPACE(file_unit)
295         reading2_seq: SELECT CASE (jfile) ! ====================
296      CASE(1) reading2_seq      ! wvl -----------------------------
297         read(file_unit,*) wvl
298         jfile = jfile+1
299      CASE(2) reading2_seq      ! radiusdyn -----------------------
300         read(file_unit,*) radiusdyn
301         jfile = jfile+1
302      CASE(3) reading2_seq      ! ep ------------------------------
303         isize = 1
304         DO WHILE (isize .le. nsize(iaer,idomain))
305            READ(file_unit,*) scanline
306            IF ((scanline(1:1) .ne. '#').and.&
307            (scanline(1:1) .ne. ' ')) THEN
308            BACKSPACE(file_unit)
309            read(file_unit,*) ep(:,isize)
310            isize = isize + 1
311         ENDIF
312      ENDDO
313
314      jfile = jfile+1
315      CASE(4) reading2_seq      ! omeg ----------------------------
316         isize = 1
317         DO WHILE (isize .le. nsize(iaer,idomain))
318            READ(file_unit,*) scanline
319            IF ((scanline(1:1) .ne. '#').and.&
320            (scanline(1:1) .ne. ' ')) THEN
321            BACKSPACE(file_unit)
322            read(file_unit,*) omeg(:,isize)
323            isize = isize + 1
324         ENDIF
325      ENDDO
326
327      jfile = jfile+1
328      CASE(5) reading2_seq      ! gfactor -------------------------
329         isize = 1
330         DO WHILE (isize .le. nsize(iaer,idomain))
331            READ(file_unit,*) scanline
332            IF ((scanline(1:1) .ne. '#').and.&
333            (scanline(1:1) .ne. ' ')) THEN
334            BACKSPACE(file_unit)
335            read(file_unit,*) gfactor(:,isize)
336            isize = isize + 1
337         ENDIF
338      ENDDO
339
340      jfile = jfile+1
341      print *, idomain, " ", iaer, " ", iaero_haze
342      IF ((idomain.NE.1).OR.(iaer.NE.1).OR.(jfile.EQ.6)) THEN
343         ! AF24: what does it mean? :-O
344         endwhile = .true.
345      ENDIF
346      CASE(6) reading2_seq
347         endwhile = .true.
348      CASE DEFAULT reading2_seq ! ----------------------------
349         WRITE(*,*) 'readoptprop: ',&
350         'Error while loading optical properties.'
351         CALL ABORT
352      END SELECT reading2_seq   ! ==============================
353      ENDIF
354      ENDDO
355
356!     1.4 Close the file
357
358      CLOSE(file_unit)
359
360!     1.5 If Francois' data, convert wvl to metres
361      !  if(iaer.eq.iaero_haze)then
362      !    if(forgetit)then
363      !    !   print*,'please sort out forgetit for naerkind>1'
364      !       do iwvl=1,nwvl
365      !          wvl(iwvl)=wvl(iwvl)*1.e-6
366      !       enddo
367      !    endif
368      ! endif
369
370      endif ! of if (is_master)
371
372      ! broadcast arrays to all cores
373      call bcast(wvl)
374      call bcast(radiusdyn)
375      call bcast(ep)
376      call bcast(omeg)
377      call bcast(gfactor)
378
379!==================================================================
380!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
381!==================================================================
382      domain: SELECT CASE (idomain)
383!==================================================================
384      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
385!==================================================================
386
387         lamref=lamrefvis(iaer)
388         epref=1.E+0
389
390         DO isize=1,nsize(iaer,idomain)
391
392!     Save the particle sizes
393            radiustab(iaer,idomain,isize)=radiusdyn(isize)
394
395!     Averaged optical properties (GCM channels)
396
397!            CALL aerave ( nwvl,&
398!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
399!            lamref,epref,tstellar,&
400!            L_NSPECTV,blamv,epavVI,&
401!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
402            CALL aerave_new ( nwvl,&
403            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
404            lamref,epref,tstellar,&
405            L_NSPECTV,blamv,epavVI,&
406            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
407
408!     Variable assignments (declared in radcommon)
409            DO isw=1,L_NSPECTV
410               QVISsQREF(isw,iaer,isize)=epavVI(isw)
411               gvis(isw,iaer,isize)=gavVI(isw)
412               omegavis(isw,iaer,isize)=omegavVI(isw)
413            END DO
414
415         ENDDO
416!==================================================================
417      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
418!==================================================================
419
420         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
421
422            lamref=lamrefir(iaer)
423            epref=1.E+0
424
425!     Save the particle sizes
426            radiustab(iaer,idomain,isize)=radiusdyn(isize)
427
428!     Averaged optical properties (GCM channels)
429
430!     epav is <QIR>/Qext(lamrefir) since epref=1
431!     Note: aerave also computes the extinction coefficient at
432!     the reference wavelength. This is called QREFvis or QREFir
433!     (not epref, which is a different parameter).
434!     Reference wavelengths SHOULD BE defined for each aerosol in
435!     radcommon_h.F90
436
437!            CALL aerave ( nwvl,&
438!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
439!            lamref,epref,tplanet,&
440!            L_NSPECTI,blami,epavIR,&
441!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
442            CALL aerave_new ( nwvl,&
443            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
444            lamref,epref,tplanet,&
445            L_NSPECTI,blami,epavIR,&
446            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
447
448
449!     Variable assignments (declared in radcommon)
450            DO ilw=1,L_NSPECTI
451               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
452               gir(ilw,iaer,isize)=gavIR(ilw)
453               omegair(ilw,iaer,isize)=omegavIR(ilw)
454            END DO
455
456
457         ENDDO ! isize (particle size) -------------------------------------
458
459      END SELECT domain
460
461!========================================================================
462!     3. Deallocate temporary variables that were read in the ASCII files
463!========================================================================
464
465      DEALLOCATE(wvl)             ! wvl
466      DEALLOCATE(radiusdyn)       ! radiusdyn
467      DEALLOCATE(ep)              ! ep
468      DEALLOCATE(omeg)            ! omeg
469      DEALLOCATE(gfactor)         ! g
470
471    END DO                    ! Loop on iaer
472  END DO                    ! Loop on idomain
473!========================================================================
474
475  ! cleanup
476  deallocate(file_id)
477
478end subroutine suaer_corrk
479
480end module suaer_corrk_mod
Note: See TracBrowser for help on using the repository browser.