source: trunk/LMDZ.MARS/libf/phymars/suaer.F90 @ 2316

Last change on this file since 2316 was 2311, checked in by emillour, 5 years ago

Mars GCM:
Code tidying: use getin_p() instead of getin() and use "call abort_physic"
instead of "stop" or "call abort".
EM

File size: 19.8 KB
Line 
1SUBROUTINE suaer
2use dimradmars_mod, only: longrefvis, longrefir, nsizemax, long1vis, &
3                    long2vis, long3vis, long1ir, long2ir, long1co2, &
4                    long2co2, nsun, nir,&
5                    naerkind, name_iaer, &
6                    iaer_dust_conrath,iaer_dust_doubleq,&
7                    iaer_dust_submicron,iaer_h2o_ice,&
8                    iaer_stormdust_doubleq,iaer_topdust_doubleq,&
9                    file_id,radiustab, gvis, omegavis, &
10                    QVISsQREF, gIR, omegaIR, &
11                    QIRsQREF, QREFvis, QREFir, &
12                    omegaREFvis, omegaREFir, &
13                    nsize
14use datafile_mod, only: datadir
15IMPLICIT NONE
16!==================================================================
17!     Purpose.
18!     --------
19!     initialize yomaer, the common that contains the
20!     radiative characteristics of the aerosols
21!     
22!     AUTHORS.
23!     --------
24!     Richard Fournier (1996) Francois Forget (1996)
25!     Frederic Hourdin
26!     Jean-jacques morcrette *ECMWF*
27!     MODIF Francois Forget (2000)
28!     MODIF Franck Montmessin (add water ice)
29!     MODIF J.-B. Madeleine 2008W27
30!       - Optical properties read in ASCII files
31!       - Add varying radius of the particules
32!
33!    Summary.
34!    --------
35!
36!    Read the optical properties -> Mean  -> Variable assignment
37!                  (ASCII files)                  (see yomaer.h)
38!    wvl(nwvl)                      longsun
39!    ep(nwvl)                       epav     QVISsQREF(nsun)
40!    omeg(nwvl)                     omegav   omegavis(nsun)
41!    gfactor(nwvl)                  gav      gvis(nsun)
42!   
43!==================================================================
44
45! Includes:
46
47include "callkeys.h"
48
49! Optical properties (read in external ASCII files)
50INTEGER          :: nwvl     ! Number of wavelengths in
51                             ! the domain (VIS or IR)
52REAL, DIMENSION(:),&
53  ALLOCATABLE, SAVE :: wvl       ! Wavelength axis
54REAL, DIMENSION(:),&
55  ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis
56
57REAL, DIMENSION(:,:),&
58  ALLOCATABLE, SAVE :: ep,&    ! Extinction coefficient Qext
59                       omeg,&  ! Single Scattering Albedo
60                       gfactor ! Assymetry Factor
61
62! Local variables:
63
64INTEGER :: iaer                ! Aerosol index
65INTEGER :: idomain             ! Domain index (1=VIS,2=IR)
66INTEGER :: iir                 ! IR channel index
67                               ! iir=1: 15um CO2 bands
68                               ! iir=2 : CO2 band wings
69                               ! iir=3 : 9 um band
70                               ! iir=4 : Far IR
71INTEGER :: isun                ! Solar band index
72INTEGER :: isize               ! Particle size index
73INTEGER :: jfile               ! ASCII file scan index
74INTEGER :: file_unit = 60
75LOGICAL :: file_ok, endwhile
76CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
77INTEGER :: read_ok
78
79! I/O  of "aerave" (subroutine averaging spectrally
80!   sing.scat.parameters)
81
82REAL tsun            ! Sun brightness temperature (for SW)
83REAL tsol            ! Surface reference brightness temp (LW)
84REAL longref         ! reference wavelengths
85REAL longsun(nsun+1) ! solar band boundaries
86REAL longir(nir+1)   ! IR band boundaries
87REAL epref           ! reference extinction ep
88                     ! at wavelength "longref"
89REAL epav(nir)       ! average ep
90                     ! (= <Qext>/Qext(longref) if epref=1)
91REAL omegav(nir)     ! Average sing.scat.albedo
92REAL gav(nir)        ! Average assymetry parameter
93
94!==================================================================
95!---- Please indicate the names of the optical property files below
96!     Please also choose the reference wavelengths of each aerosol
97
98      DO iaer = 1, naerkind ! Loop on aerosol kind
99        aerkind: SELECT CASE (name_iaer(iaer))
100!==================================================================
101        CASE("dust_conrath") aerkind      ! Typical dust profile
102!==================================================================
103!       Visible domain:
104        file_id(iaer,1) = 'optprop_dustvis_TM.dat'     !M.Wolff TM
105!       file_id(iaer,1) = 'optprop_dustvis_clancy.dat' !Clancy-Lee
106!       file_id(iaer,1) = 'optprop_dustvis_ockert.dat' !Ockert-Bell
107!       Infrared domain:
108        file_id(iaer,2) = 'optprop_dustir_TM.dat'      !M.Wolff TM
109!       Toon-Forget + solsir=2 using Clancy-Lee
110!       file_id(iaer,2) = 'optprop_dustir_clancy.dat'
111!       Toon-Forget + solsir=2 using Ockert-Bell
112!       file_id(iaer,2) = 'optprop_dustir_ockert.dat'
113!       Reference wavelength in the visible:
114        longrefvis(iaer)=0.67E-6
115!                     For dust: change readtesassim accordingly;
116!       Reference wavelength in the infrared:
117        longrefir(iaer)=dustrefir
118!==================================================================
119        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
120!==================================================================
121!       Visible domain:
122        file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
123!       file_id(iaer,1) = 'optprop_dustvis_n50.dat'    !Mie
124!       Infrared domain:
125        file_id(iaer,2) = 'optprop_dustir_n50.dat'     !Mie
126!       Reference wavelength in the visible:
127        longrefvis(iaer)=0.67E-6
128!       If not equal to 0.67e-6 -> change readtesassim accordingly;
129!       Reference wavelength in the infrared:
130        longrefir(iaer)=dustrefir
131!==================================================================
132        CASE("dust_submicron") aerkind   ! Small dust population
133!==================================================================
134!       Visible domain:
135        file_id(iaer,1) = 'optprop_dustvis_01um_TM.dat' !M.Wolff
136!       Infrared domain:
137        file_id(iaer,2) = 'optprop_dustir_01um_TM.dat'  !M.Wolff
138!       Reference wavelength in the visible:
139        longrefvis(iaer)=0.67E-6
140!       If not equal to 0.67e-6 -> change readtesassim accordingly;
141!       Reference wavelength in the infrared:
142        longrefir(iaer)=dustrefir
143!==================================================================
144        CASE("h2o_ice") aerkind             ! Water ice crystals
145!==================================================================
146!       Visible domain:
147        file_id(iaer,1) = 'optprop_icevis_n30.dat' !Warren
148!       file_id(iaer,1) = 'optprop_icevis.dat'     !Warren
149!       Infrared domain:
150        file_id(iaer,2) = 'optprop_iceir_n30.dat'  !Warren
151!       file_id(iaer,2) = 'optprop_iceir.dat'      !Warren
152!       Reference wavelength in the visible:
153        longrefvis(iaer)=0.67E-6  ! 1.5um OMEGA/MEx
154!       Reference wavelength in the infrared:
155        longrefir(iaer)=12.1E-6  ! 825cm-1 TES/MGS
156!==================================================================
157        CASE("stormdust_doubleq") aerkind   ! Two-moment scheme for stormdust - radiative properties
158!==================================================================
159!       Visible domain:
160        file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
161!       Infrared domain:
162        file_id(iaer,2) = 'optprop_dustir_n50.dat'     !Mie
163!       Reference wavelength in the visible:
164        longrefvis(iaer)=0.67E-6
165!       If not equal to 0.67e-6 -> change readtesassim accordingly;
166!       Reference wavelength in the infrared:
167        longrefir(iaer)=dustrefir
168!==================================================================
169        CASE("topdust_doubleq") aerkind   ! Two-moment scheme for topdust - radiative properties
170!==================================================================
171!       Visible domain:
172        file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
173!       Infrared domain:
174        file_id(iaer,2) = 'optprop_dustir_n50.dat'     !Mie
175!       Reference wavelength in the visible:
176        longrefvis(iaer)=0.67E-6
177!       If not equal to 0.67e-6 -> change readtesassim accordingly;
178!       Reference wavelength in the infrared:
179        longrefir(iaer)=dustrefir
180!==================================================================
181        END SELECT aerkind
182!==================================================================
183        WRITE(*,*) "Scatterer: ",trim(name_iaer(iaer))
184        WRITE(*,*) "  corresponding files: "
185        WRITE(*,*) "VIS: ",trim(file_id(iaer,1))
186        WRITE(*,*) "IR : ",trim(file_id(iaer,2))
187!==================================================================
188      ENDDO ! iaer (loop on aerosol kind)
189
190! Initializations:
191
192radiustab(1:naerkind,1:2,1:nsizemax)=0
193
194gvis(1:nsun,1:naerkind,1:nsizemax)=0
195omegavis(1:nsun,1:naerkind,1:nsizemax)=0
196QVISsQREF(1:nsun,1:naerkind,1:nsizemax)=0
197
198gIR(1:nir,1:naerkind,1:nsizemax)=0
199omegaIR(1:nir,1:naerkind,1:nsizemax)=0
200QIRsQREF(1:nir,1:naerkind,1:nsizemax)=0
201
202QREFvis(1:naerkind,1:nsizemax)=0
203QREFir(1:naerkind,1:nsizemax)=0
204omegaREFvis(1:naerkind,1:nsizemax)=0
205omegaREFir(1:naerkind,1:nsizemax)=0
206
207DO iaer = 1, naerkind ! Loop on aerosol kind
208  DO idomain = 1, 2   ! Loop on radiation domain (VIS or IR)
209!==================================================================
210! 1. READ OPTICAL PROPERTIES
211!==================================================================
212
213!       1.1 Open the ASCII file
214
215INQUIRE(FILE=TRIM(datadir)//&
216  '/'//TRIM(file_id(iaer,idomain)),&
217  EXIST=file_ok)
218IF(.NOT.file_ok) THEN
219  write(*,*)'Problem opening ',&
220    TRIM(file_id(iaer,idomain))
221  write(*,*)'It should be in: ',&
222    TRIM(datadir)
223  write(*,*)'1) You can change this directory address in callfis.def with'
224  write(*,*)'   datadir=/path/to/datafiles'
225  write(*,*)'2) If ',&
226    TRIM(file_id(iaer,idomain)),&
227    ' is a LMD reference datafile, it'
228  write(*,*)' can be obtained online on:'
229  write(*,*)' http://www.lmd.jussieu.fr/',&
230    '~lmdz/planets/mars/datadir'
231  write(*,*)'3) If the name of the file is wrong, you can'
232  write(*,*)' change it in file phymars/suaer.F90. Just'
233  write(*,*)' modify the variable called file_id.'
234  CALL abort_physic("suaer","missing file "//trim(file_id(iaer,idomain)),1)
235ENDIF
236OPEN(UNIT=file_unit,&
237  FILE=TRIM(datadir)//&
238  '/'//TRIM(file_id(iaer,idomain)),&
239  FORM='formatted')
240
241!       1.2 Allocate the optical property table
242
243jfile = 1
244endwhile = .false.
245DO WHILE (.NOT.endwhile)
246  READ(file_unit,*,iostat=read_ok) scanline
247  if (read_ok.ne.0) then
248    write(*,*)' readoptprop: Error reading file',&
249    TRIM(datadir)//&
250    '/'//TRIM(file_id(iaer,idomain))
251    call abort_physic("suaer","problem reading "//trim(file_id(iaer,idomain)),1)
252  endif
253  IF ((scanline(1:1) .ne. '#').and.&
254    (scanline(1:1) .ne. ' ')) THEN
255    BACKSPACE(file_unit)
256    reading1_seq: SELECT CASE (jfile) ! ====================
257    CASE(1) reading1_seq ! nwvl ----------------------------
258        read(file_unit,*,iostat=read_ok) nwvl
259        if (read_ok.ne.0) then
260          write(*,*)' readoptprop: Error while reading line:',&
261          trim(scanline)
262          write(*,*)'   of file',&
263          TRIM(datadir)//&
264          '/'//TRIM(file_id(iaer,idomain))
265          call abort_physic("suaer",&
266               "problem reading "//trim(file_id(iaer,idomain)),1)
267        endif
268        jfile = jfile+1
269    CASE(2) reading1_seq ! nsize ---------------------------
270        read(file_unit,*,iostat=read_ok) nsize(iaer,idomain)
271        if (read_ok.ne.0) then
272          write(*,*)' readoptprop: Error while reading line:',&
273          trim(scanline)
274          write(*,*)'   of file',&
275          TRIM(datadir)//&
276          '/'//TRIM(file_id(iaer,idomain))
277          call abort_physic("suaer",&
278               "problem reading "//trim(file_id(iaer,idomain)),1)
279        endif
280        endwhile = .true.
281    CASE DEFAULT reading1_seq ! ----------------------------
282        WRITE(*,*) 'readoptprop: ',&
283          'Error while loading optical properties.'
284        call abort_physic("suaer",&
285               "problem loading optical properties",1)
286    END SELECT reading1_seq ! ==============================
287  ENDIF
288ENDDO
289
290ALLOCATE(wvl(nwvl))                            ! wvl
291ALLOCATE(radiusdyn(nsize(iaer,idomain)))       ! radiusdyn
292ALLOCATE(ep(nwvl,nsize(iaer,idomain)))         ! ep
293ALLOCATE(omeg(nwvl,nsize(iaer,idomain)))       ! omeg
294ALLOCATE(gfactor(nwvl,nsize(iaer,idomain)))    ! g
295
296!       1.3 Read the data
297
298jfile = 1
299endwhile = .false.
300DO WHILE (.NOT.endwhile)
301   READ(file_unit,*) scanline
302  IF ((scanline(1:1) .ne. '#').and.&
303    (scanline(1:1) .ne. ' ')) THEN
304    BACKSPACE(file_unit)
305    reading2_seq: SELECT CASE (jfile) ! ====================
306    CASE(1) reading2_seq ! wvl -----------------------------
307        read(file_unit,*) wvl
308        jfile = jfile+1
309    CASE(2) reading2_seq ! radiusdyn -----------------------
310        read(file_unit,*) radiusdyn
311        jfile = jfile+1
312    CASE(3) reading2_seq ! ep ------------------------------
313        isize = 1
314        DO WHILE (isize .le. nsize(iaer,idomain))
315          READ(file_unit,*) scanline
316          IF ((scanline(1:1) .ne. '#').and.&
317            (scanline(1:1) .ne. ' ')) THEN
318          BACKSPACE(file_unit)
319          read(file_unit,*) ep(:,isize)
320          isize = isize + 1
321          ENDIF
322        ENDDO
323        jfile = jfile+1
324    CASE(4) reading2_seq ! omeg ----------------------------
325        isize = 1
326        DO WHILE (isize .le. nsize(iaer,idomain))
327          READ(file_unit,*) scanline
328          IF ((scanline(1:1) .ne. '#').and.&
329            (scanline(1:1) .ne. ' ')) THEN
330          BACKSPACE(file_unit)
331          read(file_unit,*) omeg(:,isize)
332          isize = isize + 1
333          ENDIF
334        ENDDO
335        jfile = jfile+1
336    CASE(5) reading2_seq ! gfactor -------------------------
337        isize = 1
338        DO WHILE (isize .le. nsize(iaer,idomain))
339          READ(file_unit,*) scanline
340          IF ((scanline(1:1) .ne. '#').and.&
341            (scanline(1:1) .ne. ' ')) THEN
342          BACKSPACE(file_unit)
343          read(file_unit,*) gfactor(:,isize)
344          isize = isize + 1
345          ENDIF
346        ENDDO
347        endwhile = .true.
348    CASE DEFAULT reading2_seq ! ----------------------------
349        WRITE(*,*) 'suaer.F90: ',&
350          'Error while loading optical properties.'
351        call abort_physic("suaer",&
352               "problem loading optical properties",1)
353    END SELECT reading2_seq ! ==============================
354  ENDIF
355ENDDO
356
357!       1.4 Close the file
358
359CLOSE(file_unit)
360
361!==================================================================
362! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
363!==================================================================
364domain: SELECT CASE (idomain)
365!==================================================================
366CASE(1) domain !                   VISIBLE DOMAIN (idomain=1)
367!==================================================================
368
369! 2.1 Parameters
370  tsun=6000.E+0
371  longsun(1)=long1vis
372  longsun(2)=long2vis
373  longsun(3)=long3vis
374  longref=longrefvis(iaer)
375  epref=1.E+0
376
377DO isize=1,nsize(iaer,idomain)
378! test that there is enough room to store the data
379 if (isize.gt.nsizemax) then
380   write(*,*) "suaer: Error ! nsizemax is too small!"
381   write(*,*) "       nsizemax=",nsizemax
382   write(*,*) "       you must increase the value of nsizemax"
383   write(*,*) "       in dimradmars_mod !"
384   call abort_physic("suaer","nsizemax too small",1)
385 endif
386! ------------------------------------------------
387! 2.2 Save the particle sizes
388  radiustab(iaer,idomain,isize)=radiusdyn(isize)
389! 2.3 Averaged optical properties (GCM channels)
390! Notice: Aerave also computes the extinction coefficient and
391!   single scattering albedo at reference wavelength
392!   (called QREFvis and OMEGAREFvis, same in the IR,
393!   and not epref, which is a different parameter);
394!   Reference wavelengths are defined for each aerosol in
395!   dimradmars_mod.
396
397  CALL aerave ( nwvl,&
398       wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
399       longref,epref,tsun,&
400       nsun,longsun, epav,omegav,gav,&
401       QREFvis(iaer,isize),omegaREFvis(iaer,isize) )
402! 2.4 Variable assignements (declared by yomaer.h)
403  DO isun=1,nsun
404    QVISsQREF(isun,iaer,isize)=epav(isun)
405    gvis(isun,iaer,isize)=gav(isun)
406    omegavis(isun,iaer,isize)=omegav(isun)
407  END DO
408! 2.5 Output display
409!  WRITE(*,*) 'Les donnees spectrales :'
410!  WRITE(*,*) 'Solaire (SW) ---->'
411!  WRITE(*,*) 'Aerosol number: ', iaer
412!  WRITE(*,*) 'Rayon aerosol: ', radiustab(iaer,idomain,isize)
413!  WRITE(*,*) '<Qext>/Qext(longrefvis) ; omega ; g'
414!  DO isun=1,nsun
415!    WRITE(*,*) QVISsQREF(isun,iaer,isize),&
416!         omegavis(isun,iaer,isize),&
417!         gvis(isun,iaer,isize)
418!  ENDDO
419!  WRITE(*,*) 'QREFvis(',iaer,isize,') = ',QREFvis(iaer,isize)
420!  WRITE(*,*) 'omegaREFvis(',iaer,isize,') = ',&
421!                                      omegaREFvis(iaer,isize)
422! ------------------------------------------------
423ENDDO
424
425!==================================================================
426CASE(2) domain !                  INFRARED DOMAIN (idomain=2)
427!==================================================================
428
429DO isize=1,nsize(iaer,idomain) ! ----------------------------------
430
431! 2.1 solsir is not used anymore; division of Qext(IR) by solsir
432!     has to be done in the input ASCII files (if necessary).
433
434! 2.2 Save the particle sizes
435  radiustab(iaer,idomain,isize)=radiusdyn(isize)
436
437! 2.3 Parameters
438
439  tsol=215.D+0
440  longir(1)=long1ir
441  longir(2)=long1co2
442  longir(3)=long2co2
443  longir(4)=long2ir
444  longref=longrefir(iaer)
445  epref=1.E+0
446
447! 2.4 Averaged optical properties (GCM channels)
448!           epav is <QIR>/Qext(longrefir) since epref=1
449! Notice: Aerave also Computes the extinction coefficient at
450!   reference wavelength (called QREFvis or QREFir,
451!   and not epref, which is a different parameter);
452!   Reference wavelengths are defined for each aerosol in
453!   dimradmar_mod.
454
455  CALL aerave ( nwvl,&
456       wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
457       longref,epref,tsol,&
458       nir-1,longir,epav,omegav,gav,&
459       QREFir(iaer,isize),omegaREFir(iaer,isize) )
460!  WRITE(*,*) 'QREFir(',iaer,isize,') = ',QREFir(iaer,isize)
461!  WRITE(*,*) 'omegaREFir(',iaer,isize,') = ',&
462!                                      omegaREFir(iaer,isize)
463
464! 2.5 Computing  <QIR>/Qext(longrefvis)
465
466  DO iir=1,nir-1
467!    WRITE(*,*) 'QIRsQREFir Channel ',iir,': ',epav(iir)
468    epav(iir)=  epav(iir) * QREFir(iaer,isize) / &
469                            QREFvis(iaer,isize)
470  ENDDO
471!  WRITE(*,*) 'Aerosol number', iaer
472!  WRITE(*,*) 'Particle size: ',radiustab(iaer,idomain,isize)
473!  WRITE(*,*) 'Rapport Solaire/IR:',&
474!             QREFvis(iaer,isize) / QREFir(iaer,isize)
475
476! 2.6 Variable assignements
477!           (variables are declared by yomaer.h)
478
479!         Single scattering properties
480!           in each of the "nir" bands
481!           (cf. dimradmars_mod)
482
483! iir=1 : central 15um CO2 bands   
484  QIRsQREF(1,iaer,isize)=epav(2)
485  omegaIR(1,iaer,isize)=omegav(2)
486  gIR(1,iaer,isize)=gav(2)
487
488! iir=2 : CO2 band wings
489!           (same properties than for central part)
490  QIRsQREF(2,iaer,isize)=epav(2)
491  omegaIR(2,iaer,isize)=omegav(2)
492  gIR(2,iaer,isize)=gav(2)
493
494! iir=3 : 9 um band [long1ir - long1co2]
495  QIRsQREF(3,iaer,isize)=epav(1)
496  omegaIR(3,iaer,isize)=omegav(1)
497  gIR(3,iaer,isize)=gav(1)
498
499! iir=4 : Far IR    [long2co2 - long2ir]
500  QIRsQREF(4,iaer,isize)=epav(3)
501  omegaIR(4,iaer,isize)=omegav(3)
502  gIR(4,iaer,isize)=gav(3)
503
504! 2.7 Output display
505
506!  WRITE(*,*) 'AEROSOL PROPERTIES: Number ',iaer
507!  WRITE(*,*) 'Thermal IR (LW) ---->'
508!  WRITE(*,*) 'Particle size: ',radiustab(iaer,idomain,isize)
509!  WRITE(*,*) '<Qext>/Qext(longrefvis) ; omega ; g'
510!  DO iir=1,nir
511!    WRITE(*,*) QIRsQREF(iir,iaer,isize),omegaIR(iir,iaer,isize),&
512!         gIR(iir,iaer,isize)
513!  ENDDO
514!  WRITE(*,*) 'CO2: <Qabs>/Qext(longrefvis) = ',&
515!       QIRsQREF(1,iaer,isize)*(1-omegaIR(1,iaer,isize))
516!  WRITE(*,*) ''
517
518ENDDO ! isize (particle size) -------------------------------------
519
520END SELECT domain
521!==================================================================
522! 3. Deallocate temporary variables (read in ASCII files)
523!==================================================================
524
525DEALLOCATE(wvl)        ! wvl
526DEALLOCATE(radiusdyn)  ! radiusdyn
527DEALLOCATE(ep)         ! ep
528DEALLOCATE(omeg)       ! omeg
529DEALLOCATE(gfactor)    ! g
530
531  END DO ! Loop on iaer
532END DO   ! Loop on idomain
533!==================================================================
534RETURN
535END
Note: See TracBrowser for help on using the repository browser.