Changeset 2817 for trunk/LMDZ.MARS/util


Ignore:
Timestamp:
Nov 14, 2022, 6:39:11 PM (2 years ago)
Author:
abierjon
Message:

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

Location:
trunk/LMDZ.MARS/util
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/util/README

    r2571 r2817  
    192192--------------------------------------------------------------------
    193193
    194 Program to compute opacities [1/km] of all aerosols from dust, water ice,
    195 stormdust and topdust present in the input file, at a wavelength given by the user.
     194Program to compute opacities [1/km] of aerosols at a wavelength given by the user.
    196195Computation is made from the aerosol mass mixing ratios and effective radius,
    197196associated with air density (rho) and files containing the aerosol
     
    199198The user have to precise the type of opacity : extinction or absorption.
    200199
    201 NB : water ice effective radius is known to be a bit inaccurate in the GCM,
    202 so reffwice can be prescribed by the user.
     200NB : this program requires to compile the module aeropt_mod.F90
     201at the same time
    203202
    204203input : diagfi.nc  / concat.nc / stats.nc kind of files
  • trunk/LMDZ.MARS/util/aeroptical.F90

    r2567 r2817  
    22! program for computation of aerosol opacities
    33! author : Antoine Bierjon, April-May 2020
     4! Rebranding in November 2022
    45!
    56!===============================================================================
    67!     PREFACE
    78!===============================================================================
    8 ! This program takes as inputs a GCM output file (diagfi,stats,concat) that
     9! This program can be used to compute the opacity of any aerosol.
     10!
     11! It takes as input a GCM output file (diagfi,stats,concat) that
    912! contains:
    10 !      - the Mass Mixing Ratios of dust (dustq), water ice (h2o_ice),
    11 !                               stormdust (rdsdustq) & topdust (topdustq)
    12 !      - their effective radius (reffdust, reffice(*),
    13 !                           reffstormdust, refftopdust)
    14 !      - the atmospheric density (rho)
    15 ! as well as a wavelength to calibrate the opacities, and the directory
     13!      - the Mass Mixing Ratios of the aerosols to be computed
     14!      - their effective radius (replaceable by a single value given by the user)
     15!      - the atmospheric density
     16! The user also inputs the wavelength to calibrate the opacities,
     17! the type of opacity (extinction or absorption) and the directory
    1618! containing the ASCII files of the aerosols' optical properties.
    1719!
    18 ! It outputs a netcdf file containing the opacities [1/km] of the dust,
    19 ! water ice and stormdust aerosols calibrated at the prescribed wavelength.
    20 ! The type of opacity (extinction or absorption) is given by the user.
     20! Finally, for each aerosol to be computed, the user inputs :
     21! <aerosol_name>,<aerosol_mmr_varname>,<aerosol_reff_varname/value>,<aerosol_rho_value>,<optpropfile_name>
    2122!
    22 ! (*) due to a high uncertainty of reffice in the GCM, the value is asked
    23 ! directly to the user (if the user returns 0, then the program reads the GCM
    24 ! file to get reffice)
     23! The program outputs a netcdf file containing the opacities [1/km] of
     24! the aerosols calibrated at the prescribed wavelength.
    2525!
    26 ! NOTA BENE:
    27 ! 1) if one wanted to add another aerosol to compute, one should look for
    28 ! the comments + NEW AER? that are disseminated all along the code to indicate
    29 ! the parts of the code that should be modified.
    30 ! 2) another enhancement of this program could be the possibility to read its
    31 ! own product files and recalibrate the opacities at another wavelength
     26! NB: this program requires to compile the module aeropt_mod.F90 at the same time
    3227!===============================================================================
    3328
    3429
    3530use netcdf
     31use aeropt_mod
    3632
    3733!===============================================================================
     
    4137implicit none ! for no implicitly typed variables
    4238
     39! User inputs
     40character(len=256) :: datadir                                        ! directory containing the ASCII optical properties files
     41real :: wvl_val                                                      ! reference wavelength of the output opacities [m]
     42character(len=3) :: opatype                                          ! opacity type : extinction (ext) or absorption (abs)
     43character(len=50),dimension(20) :: name_aer,mmrname_aer,reffname_aer ! aerosols' name, MMR varname, reff varname (or value)
     44integer,dimension(20) :: reffval_ok                                  ! to know if the user gave a variable name or a value for reff_aer
     45real,dimension(20) :: reffval_aer                                    ! value given by the user for reff_aer [m]
     46real,dimension(20) :: rho_aer                                        ! Density of the aerosols [kg.m-3]
     47character(len=128),dimension(20) :: optpropfile_aer                  ! name of the ASCII optical properties file
     48integer :: read_ok                                                   ! to know if the user input is well read
     49integer :: naerkind                                                  ! nb of aerosols to consider
    4350
    4451! GCM file
     
    5764real,dimension(:,:,:,:),allocatable :: mmr    ! aerosols mass mixing ratios [kg/kg]
    5865real,dimension(:,:,:,:),allocatable :: reff   ! aerosols effective radii [m]
    59 real :: reffwice_val                          ! value given by the user for reffwice (0 if read in the GCM file) [m]
    6066real,dimension(:,:,:,:),allocatable :: rho    ! atmospheric density [kg.m-3]
    61 integer :: naerkind                           ! nb of aerosols to consider
    62 integer :: iaer                               ! aerosol kind index (1=dust,2=water ice,3=stormdust,4=topdust) ! + NEW AER?
     67integer :: iaer                               ! Loop index for aerosols
    6368integer :: ilon,ilat,ialt,it                  ! Loop indices for coordinates
    6469
     70! Opacity computation
     71real :: Qext_val                               ! Qext after interpolations
     72real :: omeg_val                               ! omeg after interpolations
    6573
    6674! Output file
     
    7381                                               ! timedimout: stores time dimension ID number
    7482integer :: tmpvaridout                         ! temporarily stores a variable ID number
    75 
    76 real :: wvl_val                                ! reference wavelength of the output opacities (given by the user) [m]
    77 character(len=3) :: opatype                    ! opacity type : extinction (ext) or absorption (abs) (given by the user)
    7883integer :: varshape(4)                         ! stores a variable shape (of dimensions' IDs)
    7984character(len=16) :: tmpvarname                ! temporarily stores a variable name
     85character(len=128) :: tmplongname              ! temporarily stores a variable long_name attribute
    8086real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km]
    8187real,dimension(:),allocatable :: missval       ! Value to put in outfile when the reff is out of bounds
     
    8389
    8490
    85 ! Optical properties (read in external ASCII files)
    86 character(len=256) :: datadir              ! directory containing the ASCII files
    87 character(len=128) :: optpropfile          ! name of the ASCII optical properties file
    88 logical :: file_ok                         ! to know if the file can be opened
    89 integer :: file_unit = 60                  ! ASCII file ID
    90 
    91 integer :: jfile                           ! ASCII file scan index
    92 logical :: endwhile                        ! Reading loop boolean
    93 character(len=132) :: scanline             ! ASCII file scanning line
    94 integer :: read_ok                         ! to know if the line is well read
    95                    
    96 integer :: nwvl                            ! Number of wavelengths in the domain (VIS or IR)
    97 integer :: nsize                           ! Number of particle sizes stored in the file
    98 integer :: iwvl,isize                      ! Wavelength and Particle size loop indices
    99 real,dimension(:),allocatable :: wvl       ! Wavelength axis [m]
    100 real,dimension(:),allocatable :: radiusdyn ! Particle size axis [m]
    101 real,dimension(:,:),allocatable :: Qext    ! Extinction efficiency coefficient [/]
    102 real,dimension(:,:),allocatable :: omeg    ! Single scattering albedo [/]
    103 
    104 
    105 ! Opacity computation
    106 integer :: iwvl1,iwvl2,isize1,isize2     ! Wavelength and Particle size indices for the interpolation
    107 real :: coeff                            ! Interpolation coefficient
    108 real,dimension(2) :: reffQext            ! Qext after reff interpolation
    109 real :: Qext_val                         ! Qext after both interpolations
    110 real,dimension(2) :: reffomeg            ! omeg after reff interpolation
    111 real :: omeg_val                         ! omeg after both interpolations
    112 real,dimension(:),allocatable :: rho_aer ! Density of the aerosols [kg.m-3]
    113 
    11491!===============================================================================
    115 ! 0. USER INPUTS
     92! 0. USER INPUTS FOR NEW AEROPTICAL.DEF
     93
     94!1) Name of the GCM input file
     95!2) Directory where the optical properties files are
     96!3) The wavelength at which the output opacities are calibrated (in meters)
     97!4) Opacity type : extinction (ext) or absorption (abs)
     98!5) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     99!6) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     100!...
     101!N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     102
    116103!===============================================================================
    117104write(*,*) "Welcome in the aerosol opacities' computation program !"
    118105write(*,*)
    119106
    120 ! Ask the GCM file name
     107! 1) Ask the GCM file name
    121108WRITE(*,*) "Enter an input file name (GCM diagfi/stats/concat) :"
    122109READ(*,*) gcmfile
    123110
    124 ! Ask the reffwice to the user
     111! 2) Ask the directory containing the optical properties files
    125112write(*,*)""
    126 write(*,*) "The water ice effective radius computed by the GCM is known to be a bit inaccurate."
    127 write(*,*) "Which value do you want to use for it (in meters) ?"
    128 write(*,*) "(put 0 if you still want the program to read the value in "//trim(gcmfile)//")"
    129 READ(*,*) reffwice_val
    130 
    131 ! Ask the wavelength to the user
     113WRITE(*,*) "In which directory should we look for the optical properties' files ?"
     114READ(*,'(a)') datadir
     115
     116! 3) Ask the wavelength to the user
    132117write(*,*)""
    133118WRITE(*,*) "Value of the reference wavelength for the opacities' computation (in meters) ?"
    134119READ(*,*) wvl_val
    135120
    136 ! Ask the directory containing the optical properties files
    137 write(*,*)""
    138 WRITE(*,*) "In which directory should we look for the optical properties' files ?"
    139 READ(*,'(a)') datadir
    140 
    141 ! Ask the type of opacity that has to be computed
     121! 4) Ask the type of opacity that has to be computed
    142122write(*,*)""
    143123WRITE(*,*) "Do you want to compute extinction or absorption opacities ? (ext/abs)"
     
    149129endif
    150130
     131
     132
     133! 5-N) Ask for aerosols to compute to the user
     134write(*,*)""
     135write(*,*) "For which aerosols do you want to compute the opacity?"
     136write(*,*) "List of aerosols (up to 20), separated by <Enter>s, with each line containing :"
     137write(*,*) "<aerosol_name>,<aerosol_mmr>,<aerosol_reff>,<aerosol_rho>,<optpropfile_name>"
     138write(*,*) "(aerosol_mmr is a variable name, aerosol_reff can be variable name or value in meter ; aerosol_rho is a value in kg/m3)"
     139write(*,*)""
     140write(*,*) "An empty line , i.e: just <Enter>, implies end of list."
     141
     142reffval_ok(:) = 0 ! initialize reffval_ok to 0 for all aerosols
     143
     144naerkind=0
     145read(*,*,iostat=read_ok) name_aer(1),mmrname_aer(1),reffname_aer(1),rho_aer(1),optpropfile_aer(1)
     146
     147do while ((read_ok.eq.0).and.(name_aer(naerkind+1)/=""))
     148  naerkind=naerkind+1
     149 
     150  ! Check if reffname_aer is a variable name or a value (in meters)
     151  read(reffname_aer(naerkind),*,iostat=reffval_ok(naerkind)) reffval_aer(naerkind)
     152  if (reffval_ok(naerkind).eq.0) then
     153    write(*,*) "This value for reff of ",trim(name_aer(naerkind))," will be broadcasted to the whole grid :",reffval_aer(naerkind)
     154  endif
     155
     156  ! Read next line
     157  read(*,*,iostat=read_ok) name_aer(naerkind+1),mmrname_aer(naerkind+1),reffname_aer(naerkind+1),rho_aer(naerkind+1),optpropfile_aer(naerkind+1)
     158enddo
     159
     160if(naerkind==0) then
     161   write(*,*) "no aerosol... better stop now."
     162   stop
     163endif
     164
     165write(*,*) "naerkind=",naerkind
     166write(*,*) "name_aer=",name_aer(1:naerkind)
     167write(*,*) "mmrname_aer=",mmrname_aer(1:naerkind)
     168write(*,*) "reffname_aer=",reffname_aer(1:naerkind)
     169write(*,*) "rho_aer=",rho_aer(1:naerkind)
     170write(*,*) "optpropfile_aer=",optpropfile_aer(1:naerkind)
     171
     172
    151173!===============================================================================
    152174! 1. GCM FILE & OUTPUT FILE INITIALIZATIONS
     
    195217
    196218!==========================================================================
    197 ! 1.2 GCM aerosols
    198 !==========================================================================
    199 
    200 ! Number of aerosols
    201 naerkind = 4 ! dust, water ice, stormdust, topdust ! + NEW AER?
     219! 1.2 Check GCM aerosols
     220!==========================================================================
    202221
    203222! Variables length allocation
     
    213232aerok(:)=.true.
    214233
    215 !==========================================================================
    216 ! --> 1.2.1 DUST
    217 
    218 ! DUST MASS MIXING RATIO
    219 ! Check that the GCM file contains that variable
    220 ierr=nf90_inq_varid(gcmfid,"dustq",mmrvarid(1))
    221 error_text="Failed to find variable dustq in "//trim(gcmfile)&
    222             //" We'll skip the dust aerosol."
    223 if (ierr.ne.nf90_noerr) then
    224   write(*,*)trim(error_text)
    225   aerok(1)=.false.
    226 endif
    227 
    228 ! DUST EFFECTIVE RADIUS
    229 if (aerok(1)) then
     234! Loop on all aerosols
     235DO iaer = 1, naerkind
     236  ! MASS MIXING RATIO
    230237  ! Check that the GCM file contains that variable
    231   ierr=nf90_inq_varid(gcmfid,"reffdust",reffvarid(1))
    232   error_text="Failed to find variable reffdust in "//trim(gcmfile)&
    233               //" We'll skip the dust aerosol."
     238  ierr=nf90_inq_varid(gcmfid,mmrname_aer(iaer),mmrvarid(iaer))
     239  error_text="Failed to find variable "//trim(mmrname_aer(iaer))//&
     240             " in "//trim(gcmfile)//&
     241             " We'll skip the "//name_aer(iaer)//" aerosol."
    234242  if (ierr.ne.nf90_noerr) then
    235243    write(*,*)trim(error_text)
    236     aerok(1)=.false.
    237   endif
    238 endif
    239 
    240 !==========================================================================
    241 ! --> 1.2.2 WATER ICE
    242 
    243 ! WATER ICE MASS MIXING RATIO
    244 ! Check that the GCM file contains that variable
    245 ierr=nf90_inq_varid(gcmfid,"h2o_ice",mmrvarid(2))
    246 error_text="Failed to find variable h2o_ice in "//trim(gcmfile)&
    247             //" We'll skip the water ice aerosol."
    248 if (ierr.ne.nf90_noerr) then
    249   write(*,*)trim(error_text)
    250   aerok(2)=.false.
    251 endif
    252 
    253 ! WATER ICE EFFECTIVE RADIUS
    254 if (aerok(2)) then 
    255   IF (reffwice_val.eq.0) THEN
    256     ! Check that the GCM file contains that variable
    257     ierr=nf90_inq_varid(gcmfid,"reffice",reffvarid(2))
    258     error_text="Failed to find variable reffice in "//trim(gcmfile)&
    259                 //" We'll skip the water ice aerosol."
    260     if (ierr.ne.nf90_noerr) then
    261       write(*,*)trim(error_text)
    262       aerok(2)=.false.
    263     endif
    264   ENDIF
    265 endif
    266 
    267 !==========================================================================
    268 ! --> 1.2.3 STORMDUST
    269 
    270 ! STORMDUST MASS MIXING RATIO
    271 ! Check that the GCM file contains that variable
    272 ierr=nf90_inq_varid(gcmfid,"rdsdustq",mmrvarid(3))
    273 error_text="Failed to find variable rdsdustq in "//trim(gcmfile)&
    274             //" We'll skip the stormdust aerosol."
    275 if (ierr.ne.nf90_noerr) then
    276   write(*,*)trim(error_text)
    277   aerok(3)=.false.
    278 endif
    279 
    280 ! STORMDUST EFFECTIVE RADIUS
    281 if (aerok(3)) then
    282   ! Check that the GCM file contains that variable
    283   ierr=nf90_inq_varid(gcmfid,"reffstormdust",reffvarid(3))
    284   error_text="Failed to find variable reffstormdust in "//trim(gcmfile)&
    285               //" We'll skip the stormdust aerosol."
    286   if (ierr.ne.nf90_noerr) then
    287     write(*,*)trim(error_text)
    288     aerok(3)=.false.
    289   endif
    290 endif
    291 
    292 !==========================================================================
    293 ! --> 1.2.4 TOPDUST
    294 
    295 ! TOPDUST MASS MIXING RATIO
    296 ! Check that the GCM file contains that variable
    297 ierr=nf90_inq_varid(gcmfid,"topdustq",mmrvarid(4))
    298 error_text="Failed to find variable topdustq in "//trim(gcmfile)&
    299             //" We'll skip the topdust aerosol."
    300 if (ierr.ne.nf90_noerr) then
    301   write(*,*)trim(error_text)
    302   aerok(4)=.false.
    303 endif
    304 
    305 ! TOPDUST EFFECTIVE RADIUS
    306 if (aerok(4)) then
    307   ! Check that the GCM file contains that variable
    308   ierr=nf90_inq_varid(gcmfid,"refftopdust",reffvarid(4))
    309   error_text="Failed to find variable refftopdust in "//trim(gcmfile)&
    310               //" We'll skip the topdust aerosol."
    311   if (ierr.ne.nf90_noerr) then
    312     write(*,*)trim(error_text)
    313     aerok(4)=.false.
    314   endif
    315 endif
    316 
    317 !==========================================================================
    318 ! --> 1.2.5 + NEW AER?
    319 
     244    aerok(iaer)=.false.
     245  endif
     246
     247  ! EFFECTIVE RADIUS
     248  if (aerok(iaer)) then
     249    IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value
     250      ! Check that the GCM file contains that variable
     251      ierr=nf90_inq_varid(gcmfid,reffname_aer(iaer),reffvarid(iaer))
     252      error_text="Failed to find variable "//trim(reffname_aer(iaer))//&
     253                 " in "//trim(gcmfile)//&
     254                 " We'll skip the "//trim(name_aer(iaer))//" aerosol."
     255      if (ierr.ne.nf90_noerr) then
     256        write(*,*)trim(error_text)
     257        aerok(iaer)=.false.
     258      endif
     259    ENDIF
     260  endif
     261ENDDO !iaer
    320262
    321263! Check if there is still at least one aerosol to compute
    322 IF (.NOT.ANY(aerok)) THEN ! + NEW AER?
    323   write(*,*) "No aerosol among [dust/water ice/stormdust/topdust] was found in the file. Better stop now..."
     264IF (.NOT.ANY(aerok)) THEN
     265  write(*,*) "No aerosol among those listed was found in the file. Better stop now..."
    324266  stop
    325267ENDIF
    326268
    327269!==========================================================================
    328 ! 1.3 GCM ATMOSPHERIC DENSITY
     270! 1.3 Get GCM atmospheric density
    329271!==========================================================================
    330272
     
    343285
    344286
    345 !==========================================================================
    346 ! 1.4 Output variable's initializations and datasets loading
     287!===============================================================================
     288! 2. OUTPUT FILE VARIABLES
     289!===============================================================================
     290!==========================================================================
     291! 2.1 Output variable's initializations and datasets loading
    347292!==========================================================================
    348293! Define the ID shape of the output variables
     
    352297varshape(4)=timedimout
    353298
    354 ! Fill the aerosol density array
    355 allocate(rho_aer(naerkind))
    356 rho_aer(1)=2500. ! dust
    357 rho_aer(2)=920.  ! water ice
    358 rho_aer(3)=2500. ! stormdust
    359 rho_aer(4)=2500. ! topdust
    360 ! + NEW AER?
    361 
    362 DO iaer = 1, naerkind ! Loop on aerosol kind
     299
     300! Loop on all aerosols
     301DO iaer = 1, naerkind
    363302  if (aerok(iaer)) then ! check if this aerosol can be computed
    364303    write(*,*)""
    365     SELECT CASE(iaer)
    366     CASE(1) ! DUST
    367       write(*,*)"Computing dust opacity..."
    368     CASE(2) ! WATER ICE
    369       write(*,*)"Computing water ice opacity..."
    370     CASE(3) ! STORMDUST
    371       write(*,*)"Computing stormdust opacity..."
    372     CASE(4) ! TOPDUST
    373       write(*,*)"Computing topdust opacity..."
    374 ! + NEW AER?
    375     END SELECT ! iaer
     304    write(*,*)"Computing ",trim(name_aer(iaer))," opacity..."
    376305   
    377306    ! Length allocation of the input variables
     
    380309   
    381310    ! Datasets loading
    382      ! MMR
     311    ! ->MMR
    383312    ierr=NF90_GET_VAR(gcmfid,mmrvarid(iaer),mmr(:,:,:,:))
    384313    error_text="Failed to load mass mixing ratio"
     
    386315    write(*,*) "Mass mixing ratio loaded from "//trim(gcmfile)
    387316    ! Get missing value
    388     ierr=nf90_get_att(gcmfid,mmrvarid(1),"missing_value",missval(1))
     317    ierr=nf90_get_att(gcmfid,mmrvarid(iaer),"missing_value",missval(iaer))
    389318    if (ierr.ne.nf90_noerr) then
    390       missval(1) = 1.e+20
     319      missval(iaer) = 1.e+20
    391320    endif
    392321   
    393      ! REFF
    394     if (iaer.ne.2) then
     322    ! ->REFF
     323    IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value
    395324      ! Load dataset
    396325      ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,:,:))
     
    398327      call status_check(ierr,error_text)
    399328      write(*,*) "Effective radius loaded from "//trim(gcmfile)
    400    
    401     else ! water ice special case
    402       IF (reffwice_val.eq.0) THEN
    403         ! Load dataset
    404         ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,1,:)) ! reffice is actually a GCM
    405                                                               ! surface (column-averaged) variable
    406         error_text="Failed to load effective radius"
    407         call status_check(ierr,error_text)
    408         do ialt=2,altlen
    409           reff(:,:,ialt,:) = reff(:,:,1,:) ! build up along altitude axis
    410         enddo
    411         write(*,*) "Effective radius loaded from "//trim(gcmfile)
    412       ELSE ! if reffwice_val/=0
    413         reff(:,:,:,:) = reffwice_val
    414         write(*,*) "Effective radius loaded from the user input"
    415       ENDIF
    416      
    417     endif ! iaer.ne.2
     329    ELSE ! if reffval_ok(iaer)=0
     330      reff(:,:,:,:) = reffval_aer(iaer)
     331      write(*,*) "Effective radius loaded from the user input"
     332    ENDIF
    418333   
    419334 
     
    421336    ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen))
    422337   
    423  
    424 !===============================================================================
    425 ! 2. OPTICAL PROPERTIES
    426 !===============================================================================
    427 !==========================================================================
    428 ! 2.1 Open the ASCII file
    429 !==========================================================================
    430     IF (wvl_val.lt.5.e-6) THEN
    431       ! "VISIBLE" DOMAIN
    432      
    433       SELECT CASE(iaer)
    434       CASE(1) ! DUST
    435         ! Dust file
    436         optpropfile = "optprop_dustvis_TM_n50.dat"
    437       CASE(2) ! WATER ICE
    438         ! Water ice file
    439         optpropfile = "optprop_icevis_n30.dat"
    440       CASE(3) ! STORMDUST
    441         ! Dust file
    442         optpropfile = "optprop_dustvis_TM_n50.dat"
    443       CASE(4) ! TOPDUST
    444         ! Dust file
    445         optpropfile = "optprop_dustvis_TM_n50.dat"
    446 ! + NEW AER?
    447       END SELECT ! iaer
    448      
    449     ELSE ! wvl_val.ge.5.e-6
    450       ! "INFRARED" DOMAIN
    451      
    452       SELECT CASE(iaer)
    453       CASE(1) ! DUST
    454         ! Dust file
    455         optpropfile = "optprop_dustir_n50.dat"
    456       CASE(2) ! WATER ICE
    457         ! Water ice file
    458         optpropfile = "optprop_iceir_n30.dat"
    459       CASE(3) ! STORMDUST
    460         ! Dust file
    461         optpropfile = "optprop_dustir_n50.dat"
    462       CASE(4) ! TOPDUST
    463         ! Dust file
    464         optpropfile = "optprop_dustir_n50.dat"
    465 ! + NEW AER?
    466       END SELECT ! iaer
    467      
    468     ENDIF ! wvl_val
    469 
    470     INQUIRE(FILE=trim(datadir)//'/'//trim(optpropfile),EXIST=file_ok)
    471     if(.not.file_ok) then
    472       write(*,*)'Problem opening ',trim(optpropfile)
    473       write(*,*)'It should be in: ',trim(datadir)
    474       stop
    475     endif
    476     OPEN(UNIT=file_unit,FILE=trim(datadir)//'/'//trim(optpropfile),FORM='formatted')
    477 
    478 !==========================================================================
    479 ! 2.2 Allocate the optical property table
    480 !==========================================================================
    481     jfile = 1
    482     endwhile = .false.
    483     DO WHILE (.NOT.endwhile)
    484       READ(file_unit,*,iostat=read_ok) scanline
    485       if (read_ok.ne.0) then
    486         write(*,*)'Error reading file ',&
    487                   trim(datadir)//'/'//trim(optpropfile)
    488         stop
    489       endif
    490       IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
    491         BACKSPACE(file_unit)
    492 reading1_seq: SELECT CASE (jfile) ! FIRST READING SEQUENCE
    493         CASE(1) reading1_seq ! nwvl ----------------------------
    494             read(file_unit,*,iostat=read_ok) nwvl
    495             if (read_ok.ne.0) then
    496               write(*,*)'reading1_seq: Error while reading line: ',&
    497                         trim(scanline)
    498               write(*,*)'   of file ',&
    499                         trim(datadir)//'/'//trim(optpropfile)
    500               stop
    501             endif
    502             jfile = jfile+1
    503         CASE(2) reading1_seq ! nsize ---------------------------
    504             read(file_unit,*,iostat=read_ok) nsize
    505             if (read_ok.ne.0) then
    506               write(*,*)'reading1_seq: Error while reading line: ',&
    507                         trim(scanline)
    508               write(*,*)'   of file ',&
    509                         trim(datadir)//'/'//trim(optpropfile)
    510               stop
    511             endif
    512             endwhile = .true.
    513         CASE DEFAULT reading1_seq ! default --------------------
    514             write(*,*) 'reading1_seq: ',&
    515                        'Error while loading optical properties.'
    516             stop
    517         END SELECT reading1_seq ! ==============================
    518       ENDIF
    519     ENDDO ! DO WHILE(.NOT.endwhile)
    520 
    521     ALLOCATE(wvl(nwvl))             ! Wavelength axis
    522     ALLOCATE(radiusdyn(nsize))      ! Aerosol radius axis
    523     ALLOCATE(Qext(nwvl,nsize))      ! Extinction efficiency coeff
    524     ALLOCATE(omeg(nwvl,nsize))      ! Single scattering albedo
    525 
    526 !==========================================================================
    527 ! 2.3 Read the data
    528 !==========================================================================
    529     jfile = 1
    530     endwhile = .false.
    531     DO WHILE (.NOT.endwhile)
    532       READ(file_unit,*) scanline
    533       IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
    534         BACKSPACE(file_unit)
    535 reading2_seq: SELECT CASE (jfile) ! SECOND READING SEQUENCE
    536         CASE(1) reading2_seq ! wvl -----------------------------
    537             read(file_unit,*,iostat=read_ok) wvl
    538             if (read_ok.ne.0) then
    539               write(*,*)'reading2_seq: Error while reading line: ',&
    540                         trim(scanline)
    541               write(*,*)'   of file ',&
    542                         trim(datadir)//'/'//trim(optpropfile)
    543               stop
    544             endif
    545             jfile = jfile+1
    546         CASE(2) reading2_seq ! radiusdyn -----------------------
    547             read(file_unit,*,iostat=read_ok) radiusdyn
    548             if (read_ok.ne.0) then
    549               write(*,*)'reading2_seq: Error while reading line: ',&
    550                         trim(scanline)
    551               write(*,*)'   of file ',&
    552                         trim(datadir)//'/'//trim(optpropfile)
    553               stop
    554             endif
    555             jfile = jfile+1
    556         CASE(3) reading2_seq ! Qext ----------------------------
    557             isize = 1
    558             DO WHILE (isize .le. nsize)
    559               READ(file_unit,*,iostat=read_ok) scanline
    560               if (read_ok.ne.0) then
    561                 write(*,*)'reading2_seq: Error while reading line: ',&
    562                           trim(scanline)
    563                 write(*,*)'   of file ',&
    564                           trim(datadir)//'/'//trim(optpropfile)
    565                 stop
    566               endif
    567               IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
    568                 BACKSPACE(file_unit)
    569                 read(file_unit,*) Qext(:,isize)
    570                 isize = isize + 1
    571               ENDIF
    572             ENDDO
    573             jfile = jfile+1
    574         CASE(4) reading2_seq ! omeg ----------------------------
    575             isize = 1
    576             DO WHILE (isize .le. nsize)
    577               READ(file_unit,*,iostat=read_ok) scanline
    578               if (read_ok.ne.0) then
    579                 write(*,*)'reading2_seq: Error while reading line: ',&
    580                           trim(scanline)
    581                 write(*,*)'   of file ',&
    582                           trim(datadir)//'/'//trim(optpropfile)
    583                 stop
    584               endif
    585               IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
    586                 BACKSPACE(file_unit)
    587                 read(file_unit,*) omeg(:,isize)
    588                 isize = isize + 1
    589               ENDIF
    590             ENDDO
    591             endwhile = .true.
    592         CASE DEFAULT reading2_seq ! default --------------------
    593             write(*,*) 'reading2_seq: ',&
    594                        'Error while loading optical properties.'
    595             stop
    596         END SELECT reading2_seq ! ==============================
    597       ENDIF
    598     ENDDO
    599 
    600     ! Close the file
    601     CLOSE(file_unit)
    602    
    603     write(*,*)""
    604     write(*,*) "Wavelength array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
    605     write(*,*) "ranging from ",wvl(1)," to ",wvl(nwvl)," meters"
    606     write(*,*) "Effective radius array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
    607     write(*,*) "ranging from ",radiusdyn(1)," to ",radiusdyn(nsize)," meters"
    608 
    609 ! + NEW AER? : one may adapt this part to handle the properties' file of the new aerosol
    610 
    611 !==========================================================================
    612 ! 2.4 Get the optpropfile wavelength values encompassing the ref wavelength
    613 !==========================================================================
    614     iwvl=1
    615     DO WHILE ((iwvl.le.nwvl).and.(wvl(iwvl).lt.wvl_val))
    616       iwvl=iwvl+1
    617     ENDDO
    618     if ((iwvl.gt.nwvl).or.(wvl(1).gt.wvl_val)) then
    619       write(*,*) "ERROR: the reference wavelength is out of the bounds"
    620       write(*,*) "of the file ",trim(datadir)//'/'//trim(optpropfile)
    621       write(*,*) "(wvl_inf=",wvl(1),"m ; wvl_sup=",wvl(nwvl),"m)"
    622       write(*,*) "Ensure you wrote it well (in meters),"
    623       write(*,*) "or supply a new optical properties' file (change in aeroptical.F90 directly)"
    624       stop
    625     endif
    626     if (wvl(iwvl).eq.wvl_val) then
    627       ! if the ref wavelength is in the optpropfile
    628       iwvl1 = iwvl
    629       iwvl2 = iwvl
    630     else ! wvl(iwvl)>wvl_val and wvl(iwvl-1)<wvl_val
    631       iwvl1 = iwvl-1
    632       iwvl2 = iwvl
    633     endif
    634 
    635 
    636 !===============================================================================
    637 ! 3. OUTPUT FILE VARIABLES
    638 !===============================================================================
    639 !==========================================================================
    640 ! 3.1 Creation of the output individual aerosol opacities
     338!==========================================================================
     339! 2.2 Reading optical properties
     340!==========================================================================
     341
     342    CALL read_optpropfile(datadir,optpropfile_aer(iaer))
     343
     344!==========================================================================
     345! 2.3 Creation of the output aerosol opacity variable
    641346!==========================================================================
    642347    ! Switch to netcdf define mode
     
    645350    call status_check(ierr,error_text)
    646351    write(*,*)""
    647     SELECT CASE (iaer)
    648     CASE(1) ! DUST
    649       ! Insert the definition of the variable
    650       tmpvarname="opadust"
    651       ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
    652       error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
    653       call status_check(ierr,error_text)
    654       write(*,*) trim(tmpvarname)//" has been created in the outfile"
    655      
    656       ! Write the attributes
    657       if (opatype.eq."ext") then
    658         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust extinction opacity at reference wavelength")
    659       else ! opatype.eq."abs"
    660         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust absorption opacity at reference wavelength")
    661       endif
    662      
    663     CASE(2) ! WATER ICE
    664       ! Insert the definition of the variable
    665       tmpvarname="opawice"
    666       ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
    667       error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
    668       call status_check(ierr,error_text)
    669       write(*,*) trim(tmpvarname)//" has been created in the outfile"
    670      
    671       ! Write the attributes
    672       if (opatype.eq."ext") then
    673         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice extinction opacity at reference wavelength")
    674       else ! opatype.eq."abs"
    675         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice absorption opacity at reference wavelength")
    676       endif
    677      
    678     CASE(3) ! STORMDUST
    679       ! Insert the definition of the variable
    680       tmpvarname="opastormdust"
    681       ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
    682       error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
    683       call status_check(ierr,error_text)
    684       write(*,*) trim(tmpvarname)//" has been created in the outfile"
    685      
    686       ! Write the attributes
    687       if (opatype.eq."ext") then
    688         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust extinction opacity at reference wavelength")
    689       else ! opatype.eq."abs"
    690         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust absorption opacity at reference wavelength")
    691       endif
    692 
    693     CASE(4) ! TOPDUST
    694       ! Insert the definition of the variable
    695       tmpvarname="opatopdust"
    696       ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
    697       error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
    698       call status_check(ierr,error_text)
    699       write(*,*) trim(tmpvarname)//" has been created in the outfile"
    700      
    701       ! Write the attributes
    702       if (opatype.eq."ext") then
    703         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Topdust extinction opacity at reference wavelength")
    704       else ! opatype.eq."abs"
    705         ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Topdust absorption opacity at reference wavelength")
    706       endif
    707 
    708 ! + NEW AER?   
    709     END SELECT ! iaer
     352   
     353    ! Insert the definition of the variable
     354    tmpvarname="opa_"//trim(name_aer(iaer))
     355    ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
     356    error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
     357    call status_check(ierr,error_text)
     358    write(*,*) trim(tmpvarname)//" has been created in the outfile"
     359
     360    ! Write the attributes
     361    if (opatype.eq."ext") then
     362      tmplongname=trim(name_aer(iaer))//" extinction opacity at reference wavelength"
     363      ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
     364    else ! opatype.eq."abs"
     365      tmplongname=trim(name_aer(iaer))//" absorption opacity at reference wavelength"
     366      ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
     367    endif
    710368   
    711369    ierr=nf90_put_att(outfid,tmpvaridout,"units","opacity/km")
     
    718376    error_text="ERROR: Couldn't end netcdf define mode"
    719377    call status_check(ierr,error_text)
    720      
    721 !==========================================================================
    722 ! 3.2 Computation of the opacities
    723 !==========================================================================
     378
     379!==========================================================================
     380! 2.4 Computation of the opacities
     381!==========================================================================
     382    IF (reffval_ok(iaer).eq.0) THEN
     383    ! if the user gave a reff value,
     384    ! do the interpolation only once
     385      CALL interp_wvl_reff(wvl_val,reffval_aer(iaer),missval(iaer),&
     386                           Qext_val,omeg_val)
     387      write(*,*) "Qext(reff_val)=",Qext_val," ; omeg(reff_val)=",omeg_val
     388    ENDIF
     389
    724390    do ilon=1,lonlen
    725391     do ilat=1,latlen
    726       do ialt=1,altlen
    727        do it=1,timelen
    728          ! Get the optpropfile reff values encompassing the GCM reff
    729          isize=1
    730          do while((isize.le.nsize).and.(radiusdyn(isize).lt.reff(ilon,ilat,ialt,it)))
    731            isize=isize+1
    732          enddo
    733          if ((isize.gt.nsize).or.(radiusdyn(1).gt.reff(ilon,ilat,ialt,it))) then
    734 !           write(*,*) "WARNING: the GCM reff (",reff(ilon,ilat,ialt,it),"m) is out of the bounds"
    735 !           write(*,*) "of the file ",trim(datadir)//'/'//trim(optpropfile)
    736 !           write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
    737 !           write(*,*) "A missing value will be written for opacity"
    738 
    739            ! NB: this should especially handle cases when reff=0
     392      do it=1,timelen
     393       do ialt=1,altlen
     394         IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value
     395           CALL interp_wvl_reff(wvl_val,reff(ilon,ilat,ialt,it),missval(iaer),&
     396                                Qext_val,omeg_val)
     397         ENDIF
     398         
     399         if ((Qext_val.eq.missval(iaer))) then
     400           ! if the user's effective radius is
     401           ! out of the optpropfile boundaries
    740402           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
     403           write(*,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)"
     404           write(*,*)""
    741405           
    742          else if (mmr(ilon,ilat,ialt,it).eq.missval(iaer)) then
     406         else if ((mmr(ilon,ilat,ialt,it).eq.missval(iaer))) then
    743407           ! if there is a missing value in input file
    744408           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
    745409           
    746410         else
    747            if (radiusdyn(isize).eq.reff(ilon,ilat,ialt,it)) then
    748              ! if the GCM reff is exactly in the optpropfile
    749              isize1 = isize
    750              isize2 = isize
    751            else ! radius(isize)>reff and radiusdyn(isize-1)<reff
    752              isize1 = isize-1
    753              isize2 = isize
    754            endif
    755          
    756            ! Qext COMPUTATION
    757            ! Linear interpolation in effective radius
    758            if (isize2.ne.isize1) then
    759              coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
    760            else
    761              coeff = 0.
    762            endif
    763            reffQext(1) = Qext(iwvl1,isize1)+coeff*(Qext(iwvl1,isize2)-Qext(iwvl1,isize1))
    764            reffQext(2) = Qext(iwvl2,isize1)+coeff*(Qext(iwvl2,isize2)-Qext(iwvl2,isize1))
    765            ! Linear interpolation in wavelength
    766            if (iwvl2.ne.iwvl1) then
    767              coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
    768            else
    769              coeff = 0.
    770            endif
    771            Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1))
    772            
    773411           ! COMPUTATION OF THE EXTINCTION OPACITY [1/km]
    774412           opa_aer(ilon,ilat,ialt,it)= 750.*Qext_val*mmr(ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)&
    775413                                               / ( rho_aer(iaer) * reff(ilon,ilat,ialt,it) )
     414       
     415           if (opatype.eq."abs") then
     416              ! COMPUTATION OF THE ABSORPTION OPACITY [1/km]
     417              opa_aer(ilon,ilat,ialt,it)= (1 - omeg_val) * opa_aer(ilon,ilat,ialt,it)
     418           endif ! opatype.eq."abs"
     419
     420         endif
    776421             
    777              
    778            if (opatype.eq."abs") then
    779              ! omeg COMPUTATION
    780              ! Linear interpolation in effective radius
    781              if (isize2.ne.isize1) then
    782                coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
    783              else
    784                coeff = 0.
    785              endif
    786              reffomeg(1) = omeg(iwvl1,isize1)+coeff*(omeg(iwvl1,isize2)-omeg(iwvl1,isize1))
    787              reffomeg(2) = omeg(iwvl2,isize1)+coeff*(omeg(iwvl2,isize2)-omeg(iwvl2,isize1))
    788              ! Linear interpolation in wavelength
    789              if (iwvl2.ne.iwvl1) then
    790                coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
    791              else
    792                coeff = 0.
    793              endif
    794              omeg_val = reffomeg(1)+coeff*(reffomeg(2)-reffomeg(1))
    795 
    796              ! COMPUTATION OF THE ABSORPTION OPACITY [1/km]
    797              opa_aer(ilon,ilat,ialt,it)= (1 - omeg_val) * opa_aer(ilon,ilat,ialt,it)
    798            endif ! opatype.eq."abs"
    799                  
    800          endif
    801        enddo ! it
    802       enddo ! ialt
     422       enddo ! ialt
     423      enddo ! it
    803424     enddo ! ilat
    804425    enddo ! ilon
     
    810431    error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
    811432    call status_check(ierr,error_text)
    812 
     433   
    813434!==========================================================================
    814435! 3.4 Prepare next loop
     
    817438    DEALLOCATE(reff)
    818439    DEALLOCATE(opa_aer)
    819     DEALLOCATE(wvl)
    820     DEALLOCATE(radiusdyn)
    821     DEALLOCATE(Qext)
    822     DEALLOCATE(omeg)
     440    CALL end_aeropt_mod
    823441   
    824442  endif ! if aerok(iaer)
     
    840458
    841459end program aeroptical
     460
     461
     462
     463
    842464
    843465!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  • trunk/LMDZ.MARS/util/aeroptical.def

    r2443 r2817  
    11diagfi.nc
    2 0
     2~/datadir
    335.e-7
    4 ~/datadir
    54ext
    6 
     5dust,dustq,reffdust,2500,optprop_dustvis_TM_n50.dat
     6h2o_ice,h2o_ice,5.e-6,920,optprop_icevis_n30.dat
     7stormdust,rdsdustq,reffstormdust,2500,optprop_dustvis_TM_n50.dat
     8topdust,topdustq,refftopdust,2500,optprop_dustvis_TM_n50.dat
    79
    810
    911##----------------------------------------------------------
    10121) Name of the GCM input file
    11 2) Value of the water ice effective radius (in meters)
    12   (if 0, reffwice is read in the input file ;
    13    usually between 3.e-6m and 5.e-6m )
     132) Directory where the optical properties files are
    14143) The wavelength at which the output opacities are calibrated (in meters)
    15 4) Directory where the optical properties files are
    16 5) Opacity type : extinction (ext) or absorption (abs)
     154) Opacity type : extinction (ext) or absorption (abs)
     165) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     176) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     18...
     19N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
    1720
     21##----------------------------------------------------------
     22USUAL GCM OPTPROPFILES :
     23- for water ice : "optprop_icevis_n30.dat" for visible (wvl<=5.e-6) ; "optprop_iceir_n30.dat" for infrared (wvl>5.e-6)
     24- for dust (any type) : "optprop_dustvis_TM_n50.dat" for visible (wvl<=5.e-6) ; "optprop_dustir_n50.dat" for infrared (wvl>5.e-6)
  • trunk/LMDZ.MARS/util/compile

    r2451 r2817  
    2424#  in the C library and "-lnetcdff" should be replaced with "-lnetcdf")
    2525
    26 $COMPILER $COMPILER_OPTIONS $1.F90 \
     26code=$1.F90 # default
     27# some programs need to compile external modules
     28if [[ "$code" == aeroptical.F90 ]]
     29then
     30  code="aeropt_mod.F90 $code"
     31fi
     32
     33$COMPILER $COMPILER_OPTIONS $code \
    2734-I$NETCDF_HOME/include \
    2835-L$NETCDF_HOME/lib -lnetcdff \
Note: See TracChangeset for help on using the changeset viewer.