source: trunk/LMDZ.MARS/util/aeropt_mod.F90 @ 2817

Last change on this file since 2817 was 2817, checked in by abierjon, 2 years ago

Mars GCM:
Non-retrocompatible changes of the program util/aeroptical.F90 :

  • make the program more generic, as the user now specify in aeroptical.def the aerosols to be computed, as well as the necessary information for each aerosol, that are : aerosol_name,aerosol_mmr(GCM variable name),aerosol_reff(GCM variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
  • creation of a separate module aeropt_mod.F90 that can be used by other programs (e.g. for the MCD), with one subroutine reading the ASCII optical properties file, one subroutine interpolating optical properties at a given (wvl,reff) tuple, and one subroutine deallocating the module variables.
  • modify the compile program to compile this module aeropt_mod.F90 along with the aeroptical.F90 program. There is no change for the user, which just has to run the command 'compile aeroptical' like before.
  • correct the handling of user's wavelength or particle sizes that do not lie within the boundaries of the optical properties file.

AB
File size: 13.5 KB
Line 
1! Fortran module for computation of aerosol opacities
2! author : Antoine Bierjon, November 2022
3
4MODULE aeropt_mod
5
6 IMPLICIT NONE
7     
8 integer,save :: nwvl                            ! Number of wavelengths in the domain (VIS or IR)
9 integer,save :: nsize                           ! Number of particle sizes stored in the file
10 real,dimension(:),save,allocatable :: wvl       ! Wavelength axis [m]
11 real,dimension(:),save,allocatable :: radiusdyn ! Particle size axis [m]
12 real,dimension(:,:),save,allocatable :: Qext    ! Extinction efficiency coefficient [/]
13 real,dimension(:,:),save,allocatable :: omeg    ! Single scattering albedo [/]
14
15 CONTAINS
16 
17!*******************************************************************************
18 SUBROUTINE read_optpropfile(datadir,optpropfile)
19  !==============================================================================
20  ! Purpose:
21  ! Read optical properties from ASCII file given as input
22  ! and store them in Fortran arrays
23  !==============================================================================
24
25  use netcdf
26
27  implicit none
28
29  !==============================================================================
30  ! Arguments:
31  !==============================================================================
32  character(len=256), intent(in) :: datadir     ! directory containing the ASCII file
33  character(len=128), intent(in) :: optpropfile ! name of the ASCII optical properties file
34 
35  ! No intent(out) arguments, since only aeropt_mod module variable are updated
36  ! by this subroutine : wvl,radiusdyn,Qext,omeg
37 
38  !==============================================================================
39  ! Local variables:
40  !==============================================================================
41  logical :: file_ok                         ! to know if the file can be opened
42  integer :: file_unit = 60                  ! ASCII file ID
43  integer :: jfile                           ! ASCII file scan index
44  logical :: endwhile                        ! Reading loop boolean
45  character(len=132) :: scanline             ! ASCII file scanning line
46  integer :: read_ok                         ! to know if the line is well read
47  integer :: iwvl,isize                      ! Wavelength and Particle size loop indices
48
49  !==========================================================================
50  ! 1. Open the ASCII file
51  !==========================================================================
52  INQUIRE(FILE=trim(datadir)//'/'//trim(optpropfile),EXIST=file_ok)
53  if(.not.file_ok) then
54    write(*,*)'Problem opening ',trim(optpropfile)
55    write(*,*)'It should be in: ',trim(datadir)
56    stop
57  endif
58  OPEN(UNIT=file_unit,FILE=trim(datadir)//'/'//trim(optpropfile),FORM='formatted')
59
60  !==========================================================================
61  ! 2. Retrieve the optical properties' dimensions
62  !==========================================================================
63  jfile = 1
64  endwhile = .false.
65  DO WHILE (.NOT.endwhile)
66    READ(file_unit,*,iostat=read_ok) scanline
67    if (read_ok.ne.0) then
68      write(*,*)'Error reading file ',&
69                trim(datadir)//'/'//trim(optpropfile)
70      stop
71    endif
72    IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
73      BACKSPACE(file_unit)
74reading1_seq: SELECT CASE (jfile) ! FIRST READING SEQUENCE
75      CASE(1) reading1_seq ! nwvl ----------------------------
76          read(file_unit,*,iostat=read_ok) nwvl
77          if (read_ok.ne.0) then
78            write(*,*)'reading1_seq: Error while reading line: ',&
79                      trim(scanline)
80            write(*,*)'   of file ',&
81                      trim(datadir)//'/'//trim(optpropfile)
82            stop
83          endif
84          jfile = jfile+1
85      CASE(2) reading1_seq ! nsize ---------------------------
86          read(file_unit,*,iostat=read_ok) nsize
87          if (read_ok.ne.0) then
88            write(*,*)'reading1_seq: Error while reading line: ',&
89                      trim(scanline)
90            write(*,*)'   of file ',&
91                      trim(datadir)//'/'//trim(optpropfile)
92            stop
93          endif
94          endwhile = .true.
95      CASE DEFAULT reading1_seq ! default --------------------
96          write(*,*) 'reading1_seq: ',&
97                     'Error while loading optical properties.'
98          stop
99      END SELECT reading1_seq ! ==============================
100    ENDIF
101  ENDDO ! DO WHILE(.NOT.endwhile)
102
103  write(*,*) "nwvl=",nwvl," ; nsize=",nsize
104
105  ALLOCATE(wvl(nwvl))             ! Wavelength axis
106  ALLOCATE(radiusdyn(nsize))      ! Aerosol radius axis
107  ALLOCATE(Qext(nwvl,nsize))      ! Extinction efficiency coeff
108  ALLOCATE(omeg(nwvl,nsize))      ! Single scattering albedo
109
110  !==========================================================================
111  ! 3. Retrieve the optical properties
112  !==========================================================================
113  jfile = 1
114  endwhile = .false.
115  DO WHILE (.NOT.endwhile)
116    READ(file_unit,*) scanline
117    IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
118      BACKSPACE(file_unit)
119reading2_seq: SELECT CASE (jfile) ! SECOND READING SEQUENCE
120      CASE(1) reading2_seq ! wvl -----------------------------
121          read(file_unit,*,iostat=read_ok) wvl
122          if (read_ok.ne.0) then
123            write(*,*)'reading2_seq: Error while reading line: ',&
124                      trim(scanline)
125            write(*,*)'   of file ',&
126                      trim(datadir)//'/'//trim(optpropfile)
127            stop
128          endif
129          jfile = jfile+1
130      CASE(2) reading2_seq ! radiusdyn -----------------------
131          read(file_unit,*,iostat=read_ok) radiusdyn
132          if (read_ok.ne.0) then
133            write(*,*)'reading2_seq: Error while reading line: ',&
134                      trim(scanline)
135            write(*,*)'   of file ',&
136                      trim(datadir)//'/'//trim(optpropfile)
137            stop
138          endif
139          jfile = jfile+1
140      CASE(3) reading2_seq ! Qext ----------------------------
141          isize = 1
142          DO WHILE (isize .le. nsize)
143            READ(file_unit,*,iostat=read_ok) scanline
144            if (read_ok.ne.0) then
145              write(*,*)'reading2_seq: Error while reading line: ',&
146                        trim(scanline)
147              write(*,*)'   of file ',&
148                        trim(datadir)//'/'//trim(optpropfile)
149              stop
150            endif
151            IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
152              BACKSPACE(file_unit)
153              read(file_unit,*) Qext(:,isize)
154              isize = isize + 1
155            ENDIF
156          ENDDO
157          jfile = jfile+1
158      CASE(4) reading2_seq ! omeg ----------------------------
159          isize = 1
160          DO WHILE (isize .le. nsize)
161            READ(file_unit,*,iostat=read_ok) scanline
162            if (read_ok.ne.0) then
163              write(*,*)'reading2_seq: Error while reading line: ',&
164                        trim(scanline)
165              write(*,*)'   of file ',&
166                        trim(datadir)//'/'//trim(optpropfile)
167              stop
168            endif
169            IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
170              BACKSPACE(file_unit)
171              read(file_unit,*) omeg(:,isize)
172              isize = isize + 1
173            ENDIF
174          ENDDO
175          endwhile = .true.
176      CASE DEFAULT reading2_seq ! default --------------------
177          write(*,*) 'reading2_seq: ',&
178                     'Error while loading optical properties.'
179          stop
180      END SELECT reading2_seq ! ==============================
181    ENDIF
182  ENDDO
183
184  ! Close the file
185  CLOSE(file_unit)
186
187  write(*,*)""
188  write(*,*) "Wavelength array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
189  write(*,*) "ranging from ",wvl(1)," to ",wvl(nwvl)," meters"
190  write(*,*) "Effective radius array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
191  write(*,*) "ranging from ",radiusdyn(1)," to ",radiusdyn(nsize)," meters"
192 
193 
194 END SUBROUTINE read_optpropfile
195 
196 
197!*******************************************************************************
198 
199 SUBROUTINE interp_wvl_reff(wvl_val,reff_val,missval,&
200                            Qext_val,omeg_val)
201  !==============================================================================
202  ! Purpose:
203  ! Interpolate the optical properties to a given wavelength and effective radius
204  ! from properties tables
205  !==============================================================================
206
207  use netcdf
208
209  implicit none
210
211  !==============================================================================
212  ! Arguments:
213  !==============================================================================
214  real, intent(in) :: wvl_val                    ! Reference wavelength [m]
215  real, intent(in) :: reff_val                   ! Input reff [m]
216  real, intent(in) :: missval                    ! Value to put in outfile when the reff is out of bounds
217 
218  real, intent(out) :: Qext_val                  ! Qext after both interpolations
219  real, intent(out) :: omeg_val                  ! omega after both interpolations
220 
221 
222  !==============================================================================
223  ! Local variables:
224  !==============================================================================
225  integer :: iwvl,isize                 ! Wavelength and Particle size loop indices
226  integer :: iwvl1,iwvl2,isize1,isize2  ! Wavelength and Particle size indices for the interpolation
227  real :: coeff                         ! Interpolation coefficient
228  real,dimension(2) :: reffQext         ! Qext after reff interpolation
229  real,dimension(2) :: reffomeg         ! omeg after reff interpolation
230
231  !==========================================================================
232  ! 1. Get the optpropfile wavelength values encompassing the ref wavelength
233  !==========================================================================
234  iwvl=1
235  DO WHILE ((iwvl.lt.nwvl).and.(wvl(iwvl).lt.wvl_val))
236    iwvl=iwvl+1
237  ENDDO
238  if ((wvl(nwvl).lt.wvl_val).or.(wvl(1).gt.wvl_val)) then
239    write(*,*) "ERROR: the reference wavelength (",wvl_val,"m) is out of the bounds"
240    write(*,*) "of the optical properties' tables"
241    write(*,*) "(wvl_inf=",wvl(1),"m ; wvl_sup=",wvl(nwvl),"m)"
242    write(*,*) "Ensure you wrote it well (in meters),"
243    write(*,*) "or supply new optical properties' tables."
244    stop
245  endif
246  if (wvl(iwvl).eq.wvl_val) then
247    ! if the ref wavelength is in the optpropfile
248    iwvl1 = iwvl
249    iwvl2 = iwvl
250  else ! wvl(iwvl)>wvl_val and wvl(iwvl-1)<wvl_val
251    iwvl1 = iwvl-1
252    iwvl2 = iwvl
253  endif
254   
255  !==========================================================================
256  ! 2. Get the optpropfile reff values encompassing the input reff
257  !==========================================================================
258     
259  isize=1
260  DO WHILE ((isize.lt.nsize).and.(radiusdyn(isize).lt.reff_val))
261    isize=isize+1
262  ENDDO
263
264  if ((radiusdyn(nsize).lt.reff_val).or.(radiusdyn(1).gt.reff_val)) then
265    write(*,*) "WARNING: the effective radius (",reff_val,"m) is out of the bounds"
266    write(*,*) "of the optical properties' tables"
267    write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
268    write(*,*) "A missing value will be written for opacity."
269
270    ! NB: this should especially handle cases when reff=0
271    Qext_val = missval
272    omeg_val = missval
273
274  else
275   
276  !==========================================================================
277  ! 3. Interpolate the properties
278  !==========================================================================
279    if (radiusdyn(isize).eq.reff_val) then
280      ! if the user reff is exactly in the optpropfile
281      isize1 = isize
282      isize2 = isize
283    else ! radius(isize)>reff and radiusdyn(isize-1)<reff
284      isize1 = isize-1
285      isize2 = isize
286    endif
287
288    ! Qext COMPUTATION
289    ! Linear interpolation in effective radius
290    if (isize2.ne.isize1) then
291      coeff = (reff_val-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
292    else
293      coeff = 0.
294    endif
295    reffQext(1) = Qext(iwvl1,isize1)+coeff*(Qext(iwvl1,isize2)-Qext(iwvl1,isize1))
296    reffQext(2) = Qext(iwvl2,isize1)+coeff*(Qext(iwvl2,isize2)-Qext(iwvl2,isize1))
297    ! Linear interpolation in wavelength
298    if (iwvl2.ne.iwvl1) then
299      coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
300    else
301      coeff = 0.
302    endif
303    Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1))             
304
305
306    ! omeg COMPUTATION
307    ! Linear interpolation in effective radius
308    if (isize2.ne.isize1) then
309      coeff = (reff_val-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
310    else
311      coeff = 0.
312    endif
313    reffomeg(1) = omeg(iwvl1,isize1)+coeff*(omeg(iwvl1,isize2)-omeg(iwvl1,isize1))
314    reffomeg(2) = omeg(iwvl2,isize1)+coeff*(omeg(iwvl2,isize2)-omeg(iwvl2,isize1))
315    ! Linear interpolation in wavelength
316    if (iwvl2.ne.iwvl1) then
317      coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
318    else
319      coeff = 0.
320    endif
321    omeg_val = reffomeg(1)+coeff*(reffomeg(2)-reffomeg(1))
322
323  endif
324 
325 END SUBROUTINE interp_wvl_reff
326   
327!*******************************************************************************
328 SUBROUTINE end_aeropt_mod
329  !==============================================================================
330  ! Purpose:
331  ! Deallocate module variables
332  !==============================================================================
333
334  implicit none
335
336  if (allocated(wvl)) deallocate(wvl)
337  if (allocated(radiusdyn)) deallocate(radiusdyn)
338  if (allocated(Qext)) deallocate(Qext)
339  if (allocated(omeg)) deallocate(omeg)
340
341 END SUBROUTINE end_aeropt_mod
342 
343!*******************************************************************************
344
345END MODULE aeropt_mod
Note: See TracBrowser for help on using the repository browser.