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

Last change on this file since 3473 was 3377, checked in by afalco, 19 months ago

Pluto PCM:
initialize some variables to 0.
AF

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