Changeset 2817 for trunk/LMDZ.MARS/util
- Timestamp:
- Nov 14, 2022, 6:39:11 PM (2 years ago)
- Location:
- trunk/LMDZ.MARS/util
- Files:
-
- 1 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/util/README
r2571 r2817 192 192 -------------------------------------------------------------------- 193 193 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. 194 Program to compute opacities [1/km] of aerosols at a wavelength given by the user. 196 195 Computation is made from the aerosol mass mixing ratios and effective radius, 197 196 associated with air density (rho) and files containing the aerosol … … 199 198 The user have to precise the type of opacity : extinction or absorption. 200 199 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. 200 NB : this program requires to compile the module aeropt_mod.F90 201 at the same time 203 202 204 203 input : diagfi.nc / concat.nc / stats.nc kind of files -
trunk/LMDZ.MARS/util/aeroptical.F90
r2567 r2817 2 2 ! program for computation of aerosol opacities 3 3 ! author : Antoine Bierjon, April-May 2020 4 ! Rebranding in November 2022 4 5 ! 5 6 !=============================================================================== 6 7 ! PREFACE 7 8 !=============================================================================== 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 9 12 ! 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 16 18 ! containing the ASCII files of the aerosols' optical properties. 17 19 ! 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> 21 22 ! 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. 25 25 ! 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 32 27 !=============================================================================== 33 28 34 29 35 30 use netcdf 31 use aeropt_mod 36 32 37 33 !=============================================================================== … … 41 37 implicit none ! for no implicitly typed variables 42 38 39 ! User inputs 40 character(len=256) :: datadir ! directory containing the ASCII optical properties files 41 real :: wvl_val ! reference wavelength of the output opacities [m] 42 character(len=3) :: opatype ! opacity type : extinction (ext) or absorption (abs) 43 character(len=50),dimension(20) :: name_aer,mmrname_aer,reffname_aer ! aerosols' name, MMR varname, reff varname (or value) 44 integer,dimension(20) :: reffval_ok ! to know if the user gave a variable name or a value for reff_aer 45 real,dimension(20) :: reffval_aer ! value given by the user for reff_aer [m] 46 real,dimension(20) :: rho_aer ! Density of the aerosols [kg.m-3] 47 character(len=128),dimension(20) :: optpropfile_aer ! name of the ASCII optical properties file 48 integer :: read_ok ! to know if the user input is well read 49 integer :: naerkind ! nb of aerosols to consider 43 50 44 51 ! GCM file … … 57 64 real,dimension(:,:,:,:),allocatable :: mmr ! aerosols mass mixing ratios [kg/kg] 58 65 real,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]60 66 real,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? 67 integer :: iaer ! Loop index for aerosols 63 68 integer :: ilon,ilat,ialt,it ! Loop indices for coordinates 64 69 70 ! Opacity computation 71 real :: Qext_val ! Qext after interpolations 72 real :: omeg_val ! omeg after interpolations 65 73 66 74 ! Output file … … 73 81 ! timedimout: stores time dimension ID number 74 82 integer :: 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)78 83 integer :: varshape(4) ! stores a variable shape (of dimensions' IDs) 79 84 character(len=16) :: tmpvarname ! temporarily stores a variable name 85 character(len=128) :: tmplongname ! temporarily stores a variable long_name attribute 80 86 real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km] 81 87 real,dimension(:),allocatable :: missval ! Value to put in outfile when the reff is out of bounds … … 83 89 84 90 85 ! Optical properties (read in external ASCII files)86 character(len=256) :: datadir ! directory containing the ASCII files87 character(len=128) :: optpropfile ! name of the ASCII optical properties file88 logical :: file_ok ! to know if the file can be opened89 integer :: file_unit = 60 ! ASCII file ID90 91 integer :: jfile ! ASCII file scan index92 logical :: endwhile ! Reading loop boolean93 character(len=132) :: scanline ! ASCII file scanning line94 integer :: read_ok ! to know if the line is well read95 96 integer :: nwvl ! Number of wavelengths in the domain (VIS or IR)97 integer :: nsize ! Number of particle sizes stored in the file98 integer :: iwvl,isize ! Wavelength and Particle size loop indices99 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 computation106 integer :: iwvl1,iwvl2,isize1,isize2 ! Wavelength and Particle size indices for the interpolation107 real :: coeff ! Interpolation coefficient108 real,dimension(2) :: reffQext ! Qext after reff interpolation109 real :: Qext_val ! Qext after both interpolations110 real,dimension(2) :: reffomeg ! omeg after reff interpolation111 real :: omeg_val ! omeg after both interpolations112 real,dimension(:),allocatable :: rho_aer ! Density of the aerosols [kg.m-3]113 114 91 !=============================================================================== 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 116 103 !=============================================================================== 117 104 write(*,*) "Welcome in the aerosol opacities' computation program !" 118 105 write(*,*) 119 106 120 ! Ask the GCM file name107 ! 1) Ask the GCM file name 121 108 WRITE(*,*) "Enter an input file name (GCM diagfi/stats/concat) :" 122 109 READ(*,*) gcmfile 123 110 124 ! Ask the reffwice to the user111 ! 2) Ask the directory containing the optical properties files 125 112 write(*,*)"" 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 113 WRITE(*,*) "In which directory should we look for the optical properties' files ?" 114 READ(*,'(a)') datadir 115 116 ! 3) Ask the wavelength to the user 132 117 write(*,*)"" 133 118 WRITE(*,*) "Value of the reference wavelength for the opacities' computation (in meters) ?" 134 119 READ(*,*) wvl_val 135 120 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 142 122 write(*,*)"" 143 123 WRITE(*,*) "Do you want to compute extinction or absorption opacities ? (ext/abs)" … … 149 129 endif 150 130 131 132 133 ! 5-N) Ask for aerosols to compute to the user 134 write(*,*)"" 135 write(*,*) "For which aerosols do you want to compute the opacity?" 136 write(*,*) "List of aerosols (up to 20), separated by <Enter>s, with each line containing :" 137 write(*,*) "<aerosol_name>,<aerosol_mmr>,<aerosol_reff>,<aerosol_rho>,<optpropfile_name>" 138 write(*,*) "(aerosol_mmr is a variable name, aerosol_reff can be variable name or value in meter ; aerosol_rho is a value in kg/m3)" 139 write(*,*)"" 140 write(*,*) "An empty line , i.e: just <Enter>, implies end of list." 141 142 reffval_ok(:) = 0 ! initialize reffval_ok to 0 for all aerosols 143 144 naerkind=0 145 read(*,*,iostat=read_ok) name_aer(1),mmrname_aer(1),reffname_aer(1),rho_aer(1),optpropfile_aer(1) 146 147 do 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) 158 enddo 159 160 if(naerkind==0) then 161 write(*,*) "no aerosol... better stop now." 162 stop 163 endif 164 165 write(*,*) "naerkind=",naerkind 166 write(*,*) "name_aer=",name_aer(1:naerkind) 167 write(*,*) "mmrname_aer=",mmrname_aer(1:naerkind) 168 write(*,*) "reffname_aer=",reffname_aer(1:naerkind) 169 write(*,*) "rho_aer=",rho_aer(1:naerkind) 170 write(*,*) "optpropfile_aer=",optpropfile_aer(1:naerkind) 171 172 151 173 !=============================================================================== 152 174 ! 1. GCM FILE & OUTPUT FILE INITIALIZATIONS … … 195 217 196 218 !========================================================================== 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 !========================================================================== 202 221 203 222 ! Variables length allocation … … 213 232 aerok(:)=.true. 214 233 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 235 DO iaer = 1, naerkind 236 ! MASS MIXING RATIO 230 237 ! 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." 234 242 if (ierr.ne.nf90_noerr) then 235 243 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 261 ENDDO !iaer 320 262 321 263 ! 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..."264 IF (.NOT.ANY(aerok)) THEN 265 write(*,*) "No aerosol among those listed was found in the file. Better stop now..." 324 266 stop 325 267 ENDIF 326 268 327 269 !========================================================================== 328 ! 1.3 G CM ATMOSPHERIC DENSITY270 ! 1.3 Get GCM atmospheric density 329 271 !========================================================================== 330 272 … … 343 285 344 286 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 347 292 !========================================================================== 348 293 ! Define the ID shape of the output variables … … 352 297 varshape(4)=timedimout 353 298 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 301 DO iaer = 1, naerkind 363 302 if (aerok(iaer)) then ! check if this aerosol can be computed 364 303 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..." 376 305 377 306 ! Length allocation of the input variables … … 380 309 381 310 ! Datasets loading 382 !MMR311 ! ->MMR 383 312 ierr=NF90_GET_VAR(gcmfid,mmrvarid(iaer),mmr(:,:,:,:)) 384 313 error_text="Failed to load mass mixing ratio" … … 386 315 write(*,*) "Mass mixing ratio loaded from "//trim(gcmfile) 387 316 ! 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)) 389 318 if (ierr.ne.nf90_noerr) then 390 missval( 1) = 1.e+20319 missval(iaer) = 1.e+20 391 320 endif 392 321 393 !REFF394 if (iaer.ne.2) then322 ! ->REFF 323 IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value 395 324 ! Load dataset 396 325 ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,:,:)) … … 398 327 call status_check(ierr,error_text) 399 328 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 418 333 419 334 … … 421 336 ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen)) 422 337 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 641 346 !========================================================================== 642 347 ! Switch to netcdf define mode … … 645 350 call status_check(ierr,error_text) 646 351 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 710 368 711 369 ierr=nf90_put_att(outfid,tmpvaridout,"units","opacity/km") … … 718 376 error_text="ERROR: Couldn't end netcdf define mode" 719 377 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 724 390 do ilon=1,lonlen 725 391 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 740 402 opa_aer(ilon,ilat,ialt,it)=missval(iaer) 403 write(*,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)" 404 write(*,*)"" 741 405 742 else if ( mmr(ilon,ilat,ialt,it).eq.missval(iaer)) then406 else if ((mmr(ilon,ilat,ialt,it).eq.missval(iaer))) then 743 407 ! if there is a missing value in input file 744 408 opa_aer(ilon,ilat,ialt,it)=missval(iaer) 745 409 746 410 else 747 if (radiusdyn(isize).eq.reff(ilon,ilat,ialt,it)) then748 ! if the GCM reff is exactly in the optpropfile749 isize1 = isize750 isize2 = isize751 else ! radius(isize)>reff and radiusdyn(isize-1)<reff752 isize1 = isize-1753 isize2 = isize754 endif755 756 ! Qext COMPUTATION757 ! Linear interpolation in effective radius758 if (isize2.ne.isize1) then759 coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))760 else761 coeff = 0.762 endif763 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 wavelength766 if (iwvl2.ne.iwvl1) then767 coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))768 else769 coeff = 0.770 endif771 Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1))772 773 411 ! COMPUTATION OF THE EXTINCTION OPACITY [1/km] 774 412 opa_aer(ilon,ilat,ialt,it)= 750.*Qext_val*mmr(ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)& 775 413 / ( 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 776 421 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 803 424 enddo ! ilat 804 425 enddo ! ilon … … 810 431 error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile" 811 432 call status_check(ierr,error_text) 812 433 813 434 !========================================================================== 814 435 ! 3.4 Prepare next loop … … 817 438 DEALLOCATE(reff) 818 439 DEALLOCATE(opa_aer) 819 DEALLOCATE(wvl) 820 DEALLOCATE(radiusdyn) 821 DEALLOCATE(Qext) 822 DEALLOCATE(omeg) 440 CALL end_aeropt_mod 823 441 824 442 endif ! if aerok(iaer) … … 840 458 841 459 end program aeroptical 460 461 462 463 842 464 843 465 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -
trunk/LMDZ.MARS/util/aeroptical.def
r2443 r2817 1 1 diagfi.nc 2 0 2 ~/datadir 3 3 5.e-7 4 ~/datadir5 4 ext 6 5 dust,dustq,reffdust,2500,optprop_dustvis_TM_n50.dat 6 h2o_ice,h2o_ice,5.e-6,920,optprop_icevis_n30.dat 7 stormdust,rdsdustq,reffstormdust,2500,optprop_dustvis_TM_n50.dat 8 topdust,topdustq,refftopdust,2500,optprop_dustvis_TM_n50.dat 7 9 8 10 9 11 ##---------------------------------------------------------- 10 12 1) 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 ) 13 2) Directory where the optical properties files are 14 14 3) 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) 15 4) Opacity type : extinction (ext) or absorption (abs) 16 5) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name 17 6) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name 18 ... 19 N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name 17 20 21 ##---------------------------------------------------------- 22 USUAL 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 24 24 # in the C library and "-lnetcdff" should be replaced with "-lnetcdf") 25 25 26 $COMPILER $COMPILER_OPTIONS $1.F90 \ 26 code=$1.F90 # default 27 # some programs need to compile external modules 28 if [[ "$code" == aeroptical.F90 ]] 29 then 30 code="aeropt_mod.F90 $code" 31 fi 32 33 $COMPILER $COMPILER_OPTIONS $code \ 27 34 -I$NETCDF_HOME/include \ 28 35 -L$NETCDF_HOME/lib -lnetcdff \
Note: See TracChangeset
for help on using the changeset viewer.