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

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

Mars GCM:
The program util/aeroptical.F90 can now compute column-integrated optical depths
of the aerosols, in addition to the opacity profiles. The user has to specify it
('yes' or 'no') in aeroptical.def
Detailed changes :

  • update of the aeroptical.def file to ask the user if the column optical depth should be computed (non-retrocompatible change)
  • add in the init2 subroutine a computation of the layers' height delta_z, with adaptable method depending on the availability of some variables in the input file and on wether or not the input file has been zrecasted before. Preliminary validation shows that these different methods yield a +/-3% error on the output column optical depth tau_[aer]. The delta_z variable is also written in the output file.
  • add a log file for warnings in the interpolation subroutine in aeropt_mod.F90 + add some comments in the code
  • add the ouput "zzlev" (= interlayer altitudes) in libf/phymars/physiq_mod.F so that it can be directly used by aeroptical.F90 to compute delta_z

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