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

Last change on this file since 3506 was 3504, checked in by afalco, 2 weeks ago

Pluto: print only on master process.
AF

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