Changeset 2443 for trunk/LMDZ.MARS/util
- Timestamp:
- Jan 15, 2021, 5:41:46 PM (4 years ago)
- Location:
- trunk/LMDZ.MARS/util
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/util/README
r2435 r2443 191 191 -------------------------------------------------------------------- 192 192 193 Program to compute opacities [1/km] of all aerosols from dust, water ice 194 and stormdust present in the input file, at a wavelength given by the user.193 Program to compute opacities [1/km] of all aerosols from dust, water ice, 194 stormdust and topdust present in the input file, at a wavelength given by the user. 195 195 Computation is made from the aerosol mass mixing ratios and effective radius, 196 196 associated with air density (rho) and files containing the aerosol 197 197 optical properties (generally present in the GCM datadir/ directory). 198 The user have to precise the type of opacity : extinction or absorption. 198 199 199 200 NB : water ice effective radius is known to be a bit inaccurate in the GCM, … … 202 203 input : diagfi.nc / concat.nc / stats.nc kind of files 203 204 204 output file is : name_of_input_file_OPA.nc 205 output file is : 206 name_of_input_file_OPAext.nc for extinction opacities 207 name_of_input_file_OPAabs.nc for absorption opacities 208 205 209 aeroptical.e < aeroptical.def -
trunk/LMDZ.MARS/util/aeroptical.F90
r2435 r2443 8 8 ! This program takes as inputs a GCM output file (diagfi,stats,concat) that 9 9 ! contains: 10 ! - the Mass Mixing Ratios of dust (dustq), water ice (h2o_ice) 11 ! & stormdust (rdsdustq) 12 ! - their effective radius (reffdust, reffice(*), reffstormdust) 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) 13 14 ! - the atmospheric density (rho) 14 15 ! as well as a wavelength to calibrate the opacities, and the directory … … 17 18 ! It outputs a netcdf file containing the opacities [1/km] of the dust, 18 19 ! water ice and stormdust aerosols calibrated at the prescribed wavelength. 19 ! 20 ! The type of opacity (extinction or absorption) is given by the user. 20 21 ! 21 22 ! (*) due to a high uncertainty of reffice in the GCM, the value is asked … … 46 47 integer :: ierr ! [netcdf] subroutine returned error code 47 48 character(len=256) :: error_text ! text to display in case of error 48 integer :: tmpvarid ! temporarily stores a variable ID number49 49 integer :: lonlen,latlen,altlen,timelen ! nb of grid points along longitude,latitude,altitude,time 50 50 integer :: GCM_layers ! number of GCM layers … … 52 52 53 53 logical,dimension(:),allocatable :: aerok ! to know if the needed fields are in the GCM file 54 real,dimension(:,:,:,:,:),allocatable :: mmr ! aerosols mass mixing ratios [kg/kg] 55 real,dimension(:,:,:,:,:),allocatable :: reff ! aerosols effective radii [m] 54 integer,dimension(:),allocatable :: mmrvarid ! stores a MMR variable ID number 55 integer,dimension(:),allocatable :: reffvarid ! stores a reff variable ID number 56 integer :: tmpvarid ! temporarily stores a variable ID number 57 real,dimension(:,:,:,:),allocatable :: mmr ! aerosols mass mixing ratios [kg/kg] 58 real,dimension(:,:,:,:),allocatable :: reff ! aerosols effective radii [m] 56 59 real :: reffwice_val ! value given by the user for reffwice (0 if read in the GCM file) [m] 57 60 real,dimension(:,:,:,:),allocatable :: rho ! atmospheric density [kg.m-3] 58 61 integer :: naerkind ! nb of aerosols to consider 59 integer :: iaer ! aerosol kind index (1=dust, 2=water ice, 3=stormdust) ! + NEW AER?62 integer :: iaer ! aerosol kind index (1=dust,2=water ice,3=stormdust,4=topdust) ! + NEW AER? 60 63 integer :: ilon,ilat,ialt,it ! Loop indices for coordinates 61 64 62 65 63 66 ! Output file 64 character(len=128) :: outfile 65 integer :: outfid 67 character(len=128) :: outfile ! name of the netcdf output file 68 integer :: outfid ! [netcdf] file ID number 66 69 integer :: latdimout,londimout,altdimout,timedimout 67 ! latdimout: stores latitude dimension ID number 68 ! londimout: stores longitude dimension ID number 69 ! altdimout: stores altitude dimension ID number 70 ! timedimout: stores time dimension ID number 71 integer :: tmpvaridout ! temporarily stores a variable ID number 72 73 real :: wvl_val ! reference wavelength of the output opacities (given by the user) [m] 74 integer :: varshape(4) ! stores a variable shape (of dimensions' IDs) 75 character(len=16) :: tmpvarname ! temporarily stores a variable name 76 real,dimension(:,:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km] 77 real,dimension(:),allocatable :: missval ! Value to put in outfile when the reff is out of bounds 78 ! or when there is a mising value in input file 70 ! latdimout: stores latitude dimension ID number 71 ! londimout: stores longitude dimension ID number 72 ! altdimout: stores altitude dimension ID number 73 ! timedimout: stores time dimension ID number 74 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 integer :: varshape(4) ! stores a variable shape (of dimensions' IDs) 79 character(len=16) :: tmpvarname ! temporarily stores a variable name 80 real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km] 81 real,dimension(:),allocatable :: missval ! Value to put in outfile when the reff is out of bounds 82 ! or when there is a mising value in input file 79 83 80 84 … … 83 87 character(len=128) :: optpropfile ! name of the ASCII optical properties file 84 88 logical :: file_ok ! to know if the file can be opened 85 integer :: file_unit = 60 ! 89 integer :: file_unit = 60 ! ASCII file ID 86 90 87 91 integer :: jfile ! ASCII file scan index … … 96 100 real,dimension(:),allocatable :: radiusdyn ! Particle size axis [m] 97 101 real,dimension(:,:),allocatable :: Qext ! Extinction efficiency coefficient [/] 102 real,dimension(:,:),allocatable :: omeg ! Single scattering albedo [/] 98 103 99 104 … … 103 108 real,dimension(2) :: reffQext ! Qext after reff interpolation 104 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 105 112 real,dimension(:),allocatable :: rho_aer ! Density of the aerosols [kg.m-3] 106 113 … … 132 139 READ(*,'(a)') datadir 133 140 141 ! Ask the type of opacity that has to be computed 142 write(*,*)"" 143 WRITE(*,*) "Do you want to compute extinction or absorption opacities ? (ext/abs)" 144 READ(*,*) opatype 145 if ((trim(opatype).ne."ext").and.(trim(opatype).ne."abs")) then 146 write(*,*) "opatype = "//opatype 147 write(*,*) "Error: the opacity type must be ext or abs" 148 stop 149 endif 150 134 151 !=============================================================================== 135 152 ! 1. GCM FILE & OUTPUT FILE INITIALIZATIONS … … 149 166 !========================================================================== 150 167 ! --> 1.1.2 Creation of the outfile 151 outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPA.nc" 152 153 ierr=NF90_CREATE(outfile,nf90_clobber,outfid) ! NB: clobber mode=overwrite any dataset with the same file name ! 168 if (opatype.eq."ext") then 169 outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAext.nc" 170 else ! opatype.eq."abs" 171 outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAabs.nc" 172 endif 173 174 175 ierr=NF90_CREATE(outfile,IOR(nf90_clobber,nf90_64bit_offset),outfid) 176 ! NB: clobber mode=overwrite any dataset with the same file name ! 177 ! 64bit_offset enables the creation of large netcdf files with variables>2GB 154 178 error_text="Error: could not create file "//trim(outfile) 155 179 call status_check(ierr,error_text) … … 171 195 172 196 !========================================================================== 173 ! 1.2 GCM file variables197 ! 1.2 GCM aerosols 174 198 !========================================================================== 175 199 176 200 ! Number of aerosols 177 naerkind = 3 ! dust, water ice, stormdust ! + NEW AER?178 179 ! Length allocation of the variables180 allocate(mmr (naerkind,lonlen,latlen,altlen,timelen))181 allocate(reff (naerkind,lonlen,latlen,altlen,timelen))201 naerkind = 4 ! dust, water ice, stormdust, topdust ! + NEW AER? 202 203 ! Variables length allocation 204 allocate(mmrvarid(naerkind)) 205 allocate(reffvarid(naerkind)) 182 206 183 207 ! Initialize missing value … … 194 218 ! DUST MASS MIXING RATIO 195 219 ! Check that the GCM file contains that variable 196 ierr=nf90_inq_varid(gcmfid,"dustq", tmpvarid)220 ierr=nf90_inq_varid(gcmfid,"dustq",mmrvarid(1)) 197 221 error_text="Failed to find variable dustq in "//trim(gcmfile)& 198 222 //" We'll skip the dust aerosol." 199 223 if (ierr.ne.nf90_noerr) then 200 224 write(*,*)trim(error_text) 201 aerok(1)=.false. 202 else 203 ! Load datasets 204 ierr=NF90_GET_VAR(gcmfid,tmpvarid,mmr(1,:,:,:,:)) 205 error_text="Failed to load dust mass mixing ratio" 206 call status_check(ierr,error_text) 207 write(*,*) "Dust mass mixing ratio loaded from "//trim(gcmfile) 208 ! Get missing value 209 ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval(1)) 210 if (ierr.ne.nf90_noerr) then 211 missval(1) = 1.e+20 212 endif 225 aerok(1)=.false. 213 226 endif 214 227 … … 216 229 if (aerok(1)) then 217 230 ! Check that the GCM file contains that variable 218 ierr=nf90_inq_varid(gcmfid,"reffdust", tmpvarid)231 ierr=nf90_inq_varid(gcmfid,"reffdust",reffvarid(1)) 219 232 error_text="Failed to find variable reffdust in "//trim(gcmfile)& 220 233 //" We'll skip the dust aerosol." … … 222 235 write(*,*)trim(error_text) 223 236 aerok(1)=.false. 224 else225 ! Load datasets226 ierr=NF90_GET_VAR(gcmfid,tmpvarid,reff(1,:,:,:,:))227 error_text="Failed to load dust effective radius"228 call status_check(ierr,error_text)229 write(*,*) "Dust effective radius loaded from "//trim(gcmfile)230 237 endif 231 238 endif … … 236 243 ! WATER ICE MASS MIXING RATIO 237 244 ! Check that the GCM file contains that variable 238 ierr=nf90_inq_varid(gcmfid,"h2o_ice", tmpvarid)245 ierr=nf90_inq_varid(gcmfid,"h2o_ice",mmrvarid(2)) 239 246 error_text="Failed to find variable h2o_ice in "//trim(gcmfile)& 240 247 //" We'll skip the water ice aerosol." … … 242 249 write(*,*)trim(error_text) 243 250 aerok(2)=.false. 244 else245 ! Load datasets246 ierr=NF90_GET_VAR(gcmfid,tmpvarid,mmr(2,:,:,:,:))247 error_text="Failed to load water ice mass mixing ratio"248 call status_check(ierr,error_text)249 write(*,*) "Water ice mass mixing ratio loaded from "//trim(gcmfile)250 ! Get missing value251 ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval(2))252 if (ierr.ne.nf90_noerr) then253 missval(2) = 1.e+20254 endif255 251 endif 256 252 … … 259 255 IF (reffwice_val.eq.0) THEN 260 256 ! Check that the GCM file contains that variable 261 ierr=nf90_inq_varid(gcmfid,"reffice", tmpvarid)257 ierr=nf90_inq_varid(gcmfid,"reffice",reffvarid(2)) 262 258 error_text="Failed to find variable reffice in "//trim(gcmfile)& 263 259 //" We'll skip the water ice aerosol." … … 265 261 write(*,*)trim(error_text) 266 262 aerok(2)=.false. 267 else268 ! Load datasets269 ierr=NF90_GET_VAR(gcmfid,tmpvarid,reff(2,:,:,1,:)) ! reffice is actually a GCM270 ! surface (column-averaged) variable271 error_text="Failed to load water ice effective radius"272 call status_check(ierr,error_text)273 do ialt=2,altlen274 reff(2,:,:,ialt,:) = reff(2,:,:,1,:) ! build up along altitude axis275 enddo276 write(*,*) "Water ice effective radius loaded from "//trim(gcmfile)277 263 endif 278 ELSE ! if reffwice_val/=0279 reff(2,:,:,:,:) = reffwice_val280 write(*,*) "Water ice effective radius loaded from the user input"281 264 ENDIF 282 265 endif 283 284 ! Check if there is still at least one aerosol to compute285 IF (.NOT.ANY(aerok)) THEN286 write(*,*) "Neither dust nor water ice are found in the file. Better stop now..."287 stop288 ENDIF289 266 290 267 !========================================================================== … … 293 270 ! STORMDUST MASS MIXING RATIO 294 271 ! Check that the GCM file contains that variable 295 ierr=nf90_inq_varid(gcmfid,"rdsdustq", tmpvarid)272 ierr=nf90_inq_varid(gcmfid,"rdsdustq",mmrvarid(3)) 296 273 error_text="Failed to find variable rdsdustq in "//trim(gcmfile)& 297 274 //" We'll skip the stormdust aerosol." … … 299 276 write(*,*)trim(error_text) 300 277 aerok(3)=.false. 301 else302 ! Load datasets303 ierr=NF90_GET_VAR(gcmfid,tmpvarid,mmr(3,:,:,:,:))304 error_text="Failed to load stormdust mass mixing ratio"305 call status_check(ierr,error_text)306 write(*,*) "Stormdust mass mixing ratio loaded from "//trim(gcmfile)307 ! Get missing value308 ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval(3))309 if (ierr.ne.nf90_noerr) then310 missval(3) = 1.e+20311 endif312 278 endif 313 279 … … 315 281 if (aerok(3)) then 316 282 ! Check that the GCM file contains that variable 317 ierr=nf90_inq_varid(gcmfid,"reffstormdust", tmpvarid)283 ierr=nf90_inq_varid(gcmfid,"reffstormdust",reffvarid(3)) 318 284 error_text="Failed to find variable reffstormdust in "//trim(gcmfile)& 319 //" We'll skip the dust aerosol."285 //" We'll skip the stormdust aerosol." 320 286 if (ierr.ne.nf90_noerr) then 321 287 write(*,*)trim(error_text) 322 288 aerok(3)=.false. 323 else324 ! Load datasets325 ierr=NF90_GET_VAR(gcmfid,tmpvarid,reff(3,:,:,:,:))326 error_text="Failed to load stormdust effective radius"327 call status_check(ierr,error_text)328 write(*,*) "Stormdust effective radius loaded from "//trim(gcmfile)329 289 endif 330 290 endif 331 291 332 333 !========================================================================== 334 ! --> 1.2.4 + NEW AER? 335 336 337 !========================================================================== 338 ! --> 1.2.5 ATMOSPHERIC DENSITY 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 320 321 ! 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..." 324 stop 325 ENDIF 326 327 !========================================================================== 328 ! 1.3 GCM ATMOSPHERIC DENSITY 329 !========================================================================== 339 330 340 331 ! Check that the GCM file contains that variable … … 345 336 ! Length allocation for each dimension of the 4D variable 346 337 allocate(rho(lonlen,latlen,altlen,timelen)) 347 ! Load dataset s338 ! Load dataset 348 339 ierr=nf90_get_var(gcmfid,tmpvarid,rho) 349 340 error_text="Failed to load atmospheric density" 350 341 call status_check(ierr,error_text) 351 write(*,*) "Loaded atmospheric density rho from "//trim(gcmfile) 352 353 !========================================================================== 354 ! 1.3 Some output variable's initializations 342 write(*,*) "Atmospheric density rho loaded from "//trim(gcmfile) 343 344 345 !========================================================================== 346 ! 1.4 Output variable's initializations and datasets loading 355 347 !========================================================================== 356 348 ! Define the ID shape of the output variables … … 365 357 rho_aer(2)=920. ! water ice 366 358 rho_aer(3)=2500. ! stormdust 359 rho_aer(4)=2500. ! topdust 367 360 ! + NEW AER? 368 361 369 allocate(opa_aer(naerkind,lonlen,latlen,altlen,timelen)) 370 371 362 DO iaer = 1, naerkind ! Loop on aerosol kind 363 if (aerok(iaer)) then ! check if this aerosol can be computed 364 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 376 377 ! Length allocation of the input variables 378 ALLOCATE(mmr(lonlen,latlen,altlen,timelen)) 379 ALLOCATE(reff(lonlen,latlen,altlen,timelen)) 380 381 ! Datasets loading 382 ! MMR 383 ierr=NF90_GET_VAR(gcmfid,mmrvarid(iaer),mmr(:,:,:,:)) 384 error_text="Failed to load mass mixing ratio" 385 call status_check(ierr,error_text) 386 write(*,*) "Mass mixing ratio loaded from "//trim(gcmfile) 387 ! Get missing value 388 ierr=nf90_get_att(gcmfid,mmrvarid(1),"missing_value",missval(1)) 389 if (ierr.ne.nf90_noerr) then 390 missval(1) = 1.e+20 391 endif 392 393 ! REFF 394 if (iaer.ne.2) then 395 ! Load dataset 396 ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,:,:)) 397 error_text="Failed to load effective radius" 398 call status_check(ierr,error_text) 399 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 418 419 420 ! Length allocation of the output variable 421 ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen)) 422 423 372 424 !=============================================================================== 373 425 ! 2. OPTICAL PROPERTIES 374 426 !=============================================================================== 375 376 DO iaer = 1, naerkind ! Loop on aerosol kind377 if (aerok(iaer)) then ! check if this aerosol can be computed378 379 427 !========================================================================== 380 428 ! 2.1 Open the ASCII file … … 391 439 optpropfile = "optprop_icevis_n30.dat" 392 440 CASE(3) ! STORMDUST 441 ! Dust file 442 optpropfile = "optprop_dustvis_TM_n50.dat" 443 CASE(4) ! TOPDUST 393 444 ! Dust file 394 445 optpropfile = "optprop_dustvis_TM_n50.dat" … … 407 458 optpropfile = "optprop_iceir_n30.dat" 408 459 CASE(3) ! STORMDUST 460 ! Dust file 461 optpropfile = "optprop_dustir_n50.dat" 462 CASE(4) ! TOPDUST 409 463 ! Dust file 410 464 optpropfile = "optprop_dustir_n50.dat" … … 468 522 ALLOCATE(radiusdyn(nsize)) ! Aerosol radius axis 469 523 ALLOCATE(Qext(nwvl,nsize)) ! Extinction efficiency coeff 524 ALLOCATE(omeg(nwvl,nsize)) ! Single scattering albedo 470 525 471 526 !========================================================================== … … 516 571 ENDIF 517 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 518 591 endwhile = .true. 519 592 CASE DEFAULT reading2_seq ! default -------------------- … … 569 642 ! Switch to netcdf define mode 570 643 ierr=nf90_redef(outfid) 644 error_text="ERROR: Couldn't start netcdf define mode" 645 call status_check(ierr,error_text) 571 646 write(*,*)"" 572 647 SELECT CASE (iaer) … … 580 655 581 656 ! Write the attributes 582 ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust extinction opacity at reference wavelength") 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 583 662 584 663 CASE(2) ! WATER ICE … … 591 670 592 671 ! Write the attributes 593 ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice extinction opacity at reference wavelength") 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 594 677 595 678 CASE(3) ! STORMDUST … … 602 685 603 686 ! Write the attributes 604 ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust extinction opacity at reference wavelength") 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 605 707 606 708 ! + NEW AER? 607 709 END SELECT ! iaer 608 710 609 711 ierr=nf90_put_att(outfid,tmpvaridout,"units","opacity/km") 610 712 ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val) … … 614 716 ! End netcdf define mode 615 717 ierr=nf90_enddef(outfid) 616 718 error_text="ERROR: Couldn't end netcdf define mode" 719 call status_check(ierr,error_text) 720 617 721 !========================================================================== 618 722 ! 3.2 Computation of the opacities … … 624 728 ! Get the optpropfile reff values encompassing the GCM reff 625 729 isize=1 626 do while((isize.le.nsize).and.(radiusdyn(isize).lt.reff(i aer,ilon,ilat,ialt,it)))730 do while((isize.le.nsize).and.(radiusdyn(isize).lt.reff(ilon,ilat,ialt,it))) 627 731 isize=isize+1 628 732 enddo 629 if ((isize.gt.nsize).or.(radiusdyn(1).gt.reff(i aer,ilon,ilat,ialt,it))) then630 ! write(*,*) "WARNING: the GCM reff (",reff(i aer,ilon,ilat,ialt,it),"m) is out of the bounds"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" 631 735 ! write(*,*) "of the file ",trim(datadir)//'/'//trim(optpropfile) 632 736 ! write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)" … … 634 738 635 739 ! NB: this should especially handle cases when reff=0 636 opa_aer(i aer,ilon,ilat,ialt,it)=missval(iaer)740 opa_aer(ilon,ilat,ialt,it)=missval(iaer) 637 741 638 else if (mmr(i aer,ilon,ilat,ialt,it).eq.missval(iaer)) then742 else if (mmr(ilon,ilat,ialt,it).eq.missval(iaer)) then 639 743 ! if there is a missing value in input file 640 opa_aer(i aer,ilon,ilat,ialt,it)=missval(iaer)744 opa_aer(ilon,ilat,ialt,it)=missval(iaer) 641 745 642 746 else 643 if (radiusdyn(isize).eq.reff(i aer,ilon,ilat,ialt,it)) then644 ! if the GCM reff is in the optpropfile747 if (radiusdyn(isize).eq.reff(ilon,ilat,ialt,it)) then 748 ! if the GCM reff is exactly in the optpropfile 645 749 isize1 = isize 646 750 isize2 = isize … … 650 754 endif 651 755 756 ! Qext COMPUTATION 652 757 ! Linear interpolation in effective radius 653 758 if (isize2.ne.isize1) then 654 coeff = (reff(i aer,ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))759 coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1)) 655 760 else 656 761 coeff = 0. … … 658 763 reffQext(1) = Qext(iwvl1,isize1)+coeff*(Qext(iwvl1,isize2)-Qext(iwvl1,isize1)) 659 764 reffQext(2) = Qext(iwvl2,isize1)+coeff*(Qext(iwvl2,isize2)-Qext(iwvl2,isize1)) 660 661 765 ! Linear interpolation in wavelength 662 766 if (iwvl2.ne.iwvl1) then … … 667 771 Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1)) 668 772 669 ! Computation of the opacity [1/km] 670 opa_aer(iaer,ilon,ilat,ialt,it)= 750.*Qext_val*mmr(iaer,ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)& 671 / ( rho_aer(iaer) * reff(iaer,ilon,ilat,ialt,it) ) 672 773 ! COMPUTATION OF THE EXTINCTION OPACITY [1/km] 774 opa_aer(ilon,ilat,ialt,it)= 750.*Qext_val*mmr(ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)& 775 / ( rho_aer(iaer) * reff(ilon,ilat,ialt,it) ) 776 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 673 800 endif 674 801 enddo ! it … … 680 807 ! 3.3 Writing in the output file 681 808 !========================================================================== 682 ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer( iaer,:,:,:,:))809 ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(:,:,:,:)) 683 810 error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile" 684 811 call status_check(ierr,error_text) … … 687 814 ! 3.4 Prepare next loop 688 815 !========================================================================== 816 DEALLOCATE(mmr) 817 DEALLOCATE(reff) 818 DEALLOCATE(opa_aer) 689 819 DEALLOCATE(wvl) 690 820 DEALLOCATE(radiusdyn) 691 821 DEALLOCATE(Qext) 822 DEALLOCATE(omeg) 692 823 693 824 endif ! if aerok(iaer) -
trunk/LMDZ.MARS/util/aeroptical.def
r2317 r2443 3 3 5.e-7 4 4 ~/datadir 5 5 ext 6 6 7 7 … … 14 14 3) The wavelength at which the output opacities are calibrated (in meters) 15 15 4) Directory where the optical properties files are 16 5) Opacity type : extinction (ext) or absorption (abs) 16 17 -
trunk/LMDZ.MARS/util/simu_MCS.F90
r2404 r2443 1206 1206 !=============================================================================== 1207 1207 ! Observer variables to extract : 1208 IF (( trim(GCMvarname).eq."dsodust").or.(trim(GCMvarname).eq."dso") &1209 .or.(trim(GCMvarname).eq."dustq").or.(trim(GCMvarname).eq."dust") &1210 .or.(trim(GCMvarname).eq."reffdust")) THEN1208 IF ((index(GCMvarname,"dust").ne.0).or.(index(GCMvarname,"dso").ne.0)) THEN 1209 ! if the variable name contains "dust" or "dso". Especially for the targeted variables : 1210 ! dustq,dustN,dsodust,reffdust,opadust, and their equivalents for stormdust & topdust 1211 1211 OBSvarname = "dust" 1212 1212 numbin_name = "numbindust" 1213 1213 ELSE IF ((trim(GCMvarname).eq."h2o_ice").or.(trim(GCMvarname).eq."wice") & 1214 .or.(trim(GCMvarname).eq."reffice") ) THEN1214 .or.(trim(GCMvarname).eq."reffice").or.(trim(GCMvarname).eq."opawice")) THEN 1215 1215 OBSvarname = "wice" 1216 1216 numbin_name = "numbinwice"
Note: See TracChangeset
for help on using the changeset viewer.