Changeset 2029 for trunk/LMDZ.MARS/libf
- Timestamp:
- Oct 27, 2018, 5:58:25 PM (6 years ago)
- Location:
- trunk/LMDZ.MARS/libf/aeronomars
- Files:
-
- 1 added
- 3 edited
-
calchim.F90 (modified) (3 diffs)
-
photochemistry.F90 (modified) (3 diffs)
-
photolysis_mod.F90 (added)
-
photolysis_online.F (modified) (52 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
r2027 r2029 15 15 16 16 use conc_mod, only: mmean ! mean molecular mass of the atmosphere 17 USE comcstfi_h 17 use comcstfi_h 18 use photolysis_mod 18 19 19 20 implicit none … … 191 192 if (photochem) then 192 193 print*,'calchim: Read photolysis lookup table' 193 call read_phototable 194 call read_phototable ! off-line photolysis 195 print*,'calchim: Read UV absorption cross-sections' 196 call init_photolysis ! on-line photolysis 194 197 end if 195 198 ! find index of chemical tracers to use … … 600 603 ! latvl1= 22.27 601 604 ! lonvl1= -47.94 602 ! ig_vl1= 1+ int( (1.5-(latvl1-90.)* jjm/180.) -2 )*iim+ &603 ! int(1.5+(lonvl1+180)* iim/360.)605 ! ig_vl1= 1+ int( (1.5-(latvl1-90.)*48./180.) -2 )*64. + & 606 ! int(1.5+(lonvl1+180)*64./360.) 604 607 605 608 !======================================================================= -
trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
r2028 r2029 60 60 !=================================================================== 61 61 62 integer :: phychemrat ! (physical timestep)/(nominal chemical timestep) 63 integer :: j_o3_o1d, ilev 64 integer :: iesp, nesp 62 integer, parameter :: nesp = 17 ! number of species in the chemical code 63 64 integer :: phychemrat ! (physical timestep)/(nominal chemical timestep) 65 integer :: j_o3_o1d, ilev, iesp 65 66 integer :: lswitch 66 67 67 logical, save :: firstcall = .true. 68 68 logical :: jonline 69 70 parameter (nesp = 17) ! number of species in the chemistry71 69 72 70 ! tracer indexes in the chemistry: … … 139 137 140 138 !=================================================================== 141 ! interpolation of photolysis rates in the lookup table139 ! photolysis rates 142 140 !=================================================================== 143 141 … … 145 143 146 144 if (jonline) then 147 tau = tau*press(1)/6.1 ! temporary 148 call photolysis_online(nlayer, alt, press, temp, zmmean, & 149 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 150 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 151 i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, & 152 tau, sza, dist_sol, v_phot) 145 if (sza <= 120.) then ! day 146 tau = tau*press(1)/6.1 ! temporary 147 call photolysis_online(nlayer, alt, press, temp, zmmean, & 148 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 149 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 150 i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, & 151 tau, sza, dist_sol, v_phot) 152 else ! night 153 v_phot(:,:) = 0. 154 end if 153 155 else 154 156 call photolysis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, & -
trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F
r2028 r2029 7 7 $ tau, sza, dist_sol, v_phot) 8 8 9 use photolysis_mod 10 9 11 implicit none 10 12 … … 13 15 ! input 14 16 15 integer :: nesp! total number of chemical species16 integer :: nlayer17 integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,18 $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o,19 $ i_n, i_n2d, i_no, i_no2, i_n220 21 real, dimension(nlayer) :: press, temp, zmmean! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)22 real, dimension(nlayer) :: alt! altitude (km)23 real, dimension(nlayer,nesp) :: rm! mixing ratios24 real :: tau! integrated aerosol optical depth at the surface25 real :: sza! solar zenith angle (degrees)26 real :: dist_sol! solar distance (au)17 integer, intent(in) :: nesp ! total number of chemical species 18 integer, intent(in) :: nlayer 19 integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, 20 $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, 21 $ i_n, i_n2d, i_no, i_no2, i_n2 22 23 real, dimension(nlayer), intent(in) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1) 24 real, dimension(nlayer), intent(in) :: alt ! altitude (km) 25 real, dimension(nlayer,nesp), intent(in) :: rm ! mixing ratios 26 real, intent(in) :: tau ! integrated aerosol optical depth at the surface 27 real, intent(in) :: sza ! solar zenith angle (degrees) 28 real, intent(in) :: dist_sol ! solar distance (au) 27 29 28 30 ! output … … 30 32 real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1) 31 33 32 ! spectral grid 33 34 ! integer, parameter :: nw = 3789 ! number of spectral intervals (high-res) 35 integer, parameter :: nw = 152 ! number of spectral intervals (low-res) 36 real, dimension(nw), save :: wl, wc, wu ! lower, center, upper wavelength for each interval 37 38 ! solar flux 39 40 real, dimension(nw), save :: f ! solar flux (w.m-2.nm-1) at 1 au 34 ! solar flux at mars 35 36 real, dimension(nw) :: fmars ! solar flux (w.m-2.nm-1) 41 37 real :: factor 42 38 43 ! cross-sections and quantum yields39 ! cross-sections 44 40 45 41 integer, parameter :: nabs = 10 ! number of absorbing gases 46 42 integer, parameter :: nphot = 13 ! number of photolysis 47 43 real, dimension(nlayer,nw,nphot) :: sj ! general cross-section array (cm2) 48 real, dimension(nw), save :: xsco2_195, xsco2_295, xsco2_370 ! co2 absorption cross-section at 195-295-370 k (cm2)49 real, dimension(nw), save :: yieldco2 ! co2 photodissociation yield50 real, dimension(nw), save :: xso2_150, xso2_200, ! o2 absorption cross-section at 150-200-250-300 k (cm2)51 $ xso2_250, xso2_30052 real, dimension(nw), save :: yieldo2 ! o2 photodissociation yield53 real, dimension(nw), save :: xso3_218, xso3_298 ! o3 absorption cross-section at 218-298 k (cm2)54 real, dimension(nw), save :: xsh2o ! h2o absorption cross-section (cm2)55 real, dimension(nw), save :: xsh2o2 ! h2o2 absorption cross-section (cm2)56 real, dimension(nw), save :: xsho2 ! ho2 absorption cross-section (cm2)57 real, dimension(nw), save :: xsh2 ! h2 absorption cross-section (cm2)58 real, dimension(nw), save :: yieldh2 ! h2 photodissociation yield59 real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)60 real, dimension(nw), save :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k61 real, dimension(nw), save :: xsno ! no absorption cross-section (cm2)62 real, dimension(nw), save :: yieldno ! no photodissociation yield63 real, dimension(nw), save :: xsn2 ! n2 absorption cross-section (cm2)64 real, dimension(nw), save :: yieldn2 ! n2 photodissociation yield65 real, dimension(nw), save :: albedo ! surface albedo66 44 67 45 ! atmosphere … … 90 68 $ a_no2, a_n2 91 69 integer :: i, ilay, iw, ialt 92 integer, save :: mopt93 70 real :: deltaj 94 95 logical, save :: firstcall = .true.96 71 97 72 ! absorbing gas numbering … … 125 100 j_n2 = 13 ! n2 + hv -> n + n 126 101 127 !===== day/night criterion128 129 if (sza <= 120.) then130 131 !===== read once absorption cross-sections, quantum yields, and albedo132 133 if (firstcall) then134 135 ! high-res/low-res switch136 !137 ! mopt = 1 high-resolution138 ! mopt = 2 low-resolution (recommended for gcm use)139 140 mopt = 2141 142 ! set wavelength grid143 144 call gridw(nw,wl,wc,wu,mopt)145 146 ! read and grid solar flux data147 148 call rdsolarflux(nw,wl,wc,f)149 150 ! read and grid o2 cross-sections151 152 call rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)153 154 ! read and grid co2 cross-sections155 156 call rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)157 158 ! read and grid o3 cross-sections159 160 call rdxso3(nw,wl,xso3_218,xso3_298)161 162 ! read and grid h2o cross-sections163 164 call rdxsh2o(nw,wl,xsh2o)165 166 ! read and grid h2o2 cross-sections167 168 call rdxsh2o2(nw,wl,xsh2o2)169 170 ! read and grid ho2 cross-sections171 172 call rdxsho2(nw,wl,xsho2)173 174 ! read and grid h2 cross-sections175 176 call rdxsh2(nw,wl,wc,xsh2,yieldh2)177 178 ! read and grid no cross-sections179 180 call rdxsno(nw,wl,xsno,yieldno)181 182 ! read and grid no2 cross-sections183 184 call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,185 $ yldno2_248, yldno2_298)186 187 ! read and grid n2 cross-sections188 189 call rdxsn2(nw,wl,xsn2,yieldn2)190 191 ! set surface albedo192 193 call setalb(nw,wl,albedo)194 195 firstcall = .false.196 197 end if198 199 102 !==== air column increments and rayleigh optical depth 200 103 … … 270 173 end do 271 174 272 ! ialt = 49273 ! print*, 'altitude = ', alt(ialt), ' km'274 ! do iw = 1,nw275 ! write(35,'(f8.4,12e13.5)') wc(iw), dtrl(ialt,iw),276 ! $ (max(dtgas(ialt,iw,i),1.e-30),i = 1,nabs)277 ! end do278 279 175 !==== set aerosol properties and optical depth 280 176 … … 287 183 !==== slant path lengths in spherical geometry 288 184 289 call sphers(nlayer, alt, sza, dsdh,nid)185 call sphers(nlayer,alt,sza,dsdh,nid) 290 186 291 187 !==== solar flux at mars … … 293 189 factor = (1./dist_sol)**2. 294 190 do iw = 1,nw-1 295 f (iw) = f(iw)*factor191 fmars(iw) = f(iw)*factor 296 192 end do 297 193 … … 318 214 319 215 do ilay = 1,nlayer 320 saflux(ilay) = f(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))216 saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay)) 321 217 end do 322 218 … … 338 234 end do ! iw 339 235 340 else ! night 341 342 v_phot(:,1:nphot) = 0. 343 344 end if ! sza 345 346 end subroutine photolysis_online 347 348 !============================================================================== 349 350 subroutine gridw(nw,wl,wc,wu,mopt) 351 352 ! Create the wavelength grid for all interpolations and radiative transfer 353 ! calculations. Grid may be irregularly spaced. Wavelengths are in nm. 354 ! No gaps are allowed within the wavelength grid. 355 356 implicit none 357 358 ! input 359 360 integer :: nw ! number of wavelength grid points 361 integer :: mopt ! high-res/low-res switch 362 363 ! output 364 365 real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval 366 367 ! local 368 369 real :: wincr ! wavelength increment 370 integer :: iw, kw 371 372 ! mopt = 1 high-resolution mode (3789 intervals) 373 ! 374 ! 0-108 nm : 1.0 nm 375 ! 108-124 nm : 0.1 nm 376 ! 124-175 nm : 0.5 nm 377 ! 175-205 nm : 0.01 nm 378 ! 205-365 nm : 0.5 nm 379 ! 365-850 nm : 5.0 nm 380 ! 381 ! mopt = 2 low-resolution mode 382 ! 383 ! 0-60 nm : 6.0 nm 384 ! 60-80 nm : 2.0 nm 385 ! 80-120 nm : 5.0 nm 386 ! 120-123 nm : 0.2 nm 387 ! 123-163 nm : 5.0 nm 388 ! 163-175 nm : 2.0 nm 389 ! 175-205 nm : 0.5 nm 390 ! 205-245 nm : 5.0 nm 391 ! 245-415 nm : 10.0 nm 392 ! 415-815 nm : 50.0 nm 393 394 if (mopt == 1) then ! high-res 395 396 ! define wavelength intervals of width 1.0 nm from 0 to 108 nm: 397 398 kw = 0 399 wincr = 1.0 400 do iw = 0, 107 401 kw = kw + 1 402 wl(kw) = real(iw) 403 wu(kw) = wl(kw) + wincr 404 wc(kw) = (wl(kw) + wu(kw))/2. 405 end do 406 407 ! define wavelength intervals of width 0.1 nm from 108 to 124 nm: 408 409 wincr = 0.1 410 do iw = 1080, 1239, 1 411 kw = kw + 1 412 wl(kw) = real(iw)/10. 413 wu(kw) = wl(kw) + wincr 414 wc(kw) = (wl(kw) + wu(kw))/2. 415 end do 416 417 ! define wavelength intervals of width 0.5 nm from 124 to 175 nm: 418 419 wincr = 0.5 420 do iw = 1240, 1745, 5 421 kw = kw + 1 422 wl(kw) = real(iw)/10. 423 wu(kw) = wl(kw) + wincr 424 wc(kw) = (wl(kw) + wu(kw))/2. 425 end do 426 427 ! define wavelength intervals of width 0.01 nm from 175 to 205 nm: 428 429 wincr = 0.01 430 do iw = 17500, 20499, 1 431 kw = kw + 1 432 wl(kw) = real(iw)/100. 433 wu(kw) = wl(kw) + wincr 434 wc(kw) = (wl(kw) + wu(kw))/2. 435 end do 436 437 ! define wavelength intervals of width 0.5 nm from 205 to 365 nm: 438 439 wincr = 0.5 440 do iw = 2050, 3645, 5 441 kw = kw + 1 442 wl(kw) = real(iw)/10. 443 wu(kw) = wl(kw) + wincr 444 wc(kw) = (wl(kw) + wu(kw))/2. 445 end do 446 447 ! define wavelength intervals of width 5.0 nm from 365 to 855 nm: 448 449 wincr = 5.0 450 do iw = 365, 850, 5 451 kw = kw + 1 452 wl(kw) = real(iw) 453 wu(kw) = wl(kw) + wincr 454 wc(kw) = (wl(kw) + wu(kw))/2. 455 end do 456 wl(kw+1) = wu(kw) 457 458 else if (mopt == 2) then ! low-res 459 460 * define wavelength intervals of width 6.0 nm from 0 to 60 nm: 461 462 kw = 0 463 wincr = 6.0 464 DO iw = 0, 54, 6 465 kw = kw + 1 466 wl(kw) = real(iw) 467 wu(kw) = wl(kw) + wincr 468 wc(kw) = (wl(kw) + wu(kw))/2. 469 END DO 470 471 * define wavelength intervals of width 2.0 nm from 60 to 80 nm: 472 473 wincr = 2.0 474 DO iw = 60, 78, 2 475 kw = kw + 1 476 wl(kw) = real(iw) 477 wu(kw) = wl(kw) + wincr 478 wc(kw) = (wl(kw) + wu(kw))/2. 479 END DO 480 481 * define wavelength intervals of width 5.0 nm from 80 to 120 nm: 482 483 wincr = 5.0 484 DO iw = 80, 115, 5 485 kw = kw + 1 486 wl(kw) = real(iw) 487 wu(kw) = wl(kw) + wincr 488 wc(kw) = (wl(kw) + wu(kw))/2. 489 END DO 490 491 492 * define wavelength intervals of width 0.2 nm from 120 to 123 nm: 493 494 wincr = 0.2 495 DO iw = 1200, 1228, 2 496 kw = kw + 1 497 wl(kw) = real(iw)/10. 498 wu(kw) = wl(kw) + wincr 499 wc(kw) = (wl(kw) + wu(kw))/2. 500 ENDDO 501 502 * define wavelength intervals of width 5.0 nm from 123 to 163 nm: 503 504 wincr = 5.0 505 DO iw = 123, 158, 5 506 kw = kw + 1 507 wl(kw) = real(iw) 508 wu(kw) = wl(kw) + wincr 509 wc(kw) = (wl(kw) + wu(kw))/2. 510 ENDDO 511 512 * define wavelength intervals of width 2.0 nm from 163 to 175 nm: 513 514 wincr = 2.0 515 DO iw = 163, 173, 2 516 kw = kw + 1 517 wl(kw) = real(iw) 518 wu(kw) = wl(kw) + wincr 519 wc(kw) = (wl(kw) + wu(kw))/2. 520 ENDDO 521 522 * define wavelength intervals of width 0.5 nm from 175 to 205 nm: 523 524 wincr = 0.5 525 DO iw = 1750, 2045, 5 526 kw = kw + 1 527 wl(kw) = real(iw)/10. 528 wu(kw) = wl(kw) + wincr 529 wc(kw) = (wl(kw) + wu(kw))/2. 530 ENDDO 531 532 * define wavelength intervals of width 5.0 nm from 205 to 245 nm: 533 534 wincr = 5. 535 DO iw = 205, 240, 5 536 kw = kw + 1 537 wl(kw) = real(iw) 538 wu(kw) = wl(kw) + wincr 539 wc(kw) = (wl(kw) + wu(kw))/2. 540 ENDDO 541 542 * define wavelength intervals of width 10.0 nm from 245 to 415 nm: 543 544 wincr = 10.0 545 DO iw = 245, 405, 10 546 kw = kw + 1 547 wl(kw) = real(iw) 548 wu(kw) = wl(kw) + wincr 549 wc(kw) = (wl(kw) + wu(kw))/2. 550 ENDDO 551 552 * define wavelength intervals of width 50.0 nm from 415 to 815 nm: 553 554 wincr = 50.0 555 DO iw = 415, 815, 50 556 kw = kw + 1 557 wl(kw) = real(iw) 558 wu(kw) = wl(kw) + wincr 559 wc(kw) = (wl(kw) + wu(kw))/2. 560 ENDDO 561 562 wl(kw+1) = wu(kw) 563 564 end if ! mopt 565 566 ! do iw = 1,nw 567 ! write(20,*) iw, wl(iw), wu(iw) 568 ! end do 569 print*, 'number of spectral intervals : ', kw+1 570 571 return 572 end 573 574 !============================================================================== 575 576 subroutine rdsolarflux(nw,wl,wc,f) 577 578 ! Read and re-grid solar flux data. 579 580 use datafile_mod, only: datadir 581 582 implicit none 583 584 ! input 585 586 integer :: nw ! number of wavelength grid points 587 real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval 588 589 ! output 590 591 real, dimension(nw) :: f ! solar flux (w.m-2.nm-1) 592 593 ! local 594 595 integer, parameter :: kdata = 20000 ! max dimension of input solar flux 596 integer :: msun ! choice of solar flux 597 integer :: iw, nhead, ihead, n, i, ierr, kin 598 599 real, parameter :: deltax = 1.e-4 600 real, dimension(kdata) :: x1, y1 ! input solar flux 601 real, dimension(nw) :: yg1 ! gridded solar flux 602 603 character(len=100) :: fil 604 605 kin = 10 ! input logical unit 606 607 ! select desired extra-terrestrial solar irradiance, using msun: 608 609 ! 18 = atlas3_thuillier_tuv.txt 0-900 nm November 1994 610 ! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004 611 612 msun = 18 613 614 if (msun == 18) THEN 615 616 fil = trim(datadir)//'/solar_fluxes/atlas3_thuillier_tuv.txt' 617 print*, 'solar flux : ', fil 618 OPEN(UNIT=kin,FILE=fil,STATUS='old') 619 nhead = 9 620 n = 19193 621 DO ihead = 1, nhead 622 READ(kin,*) 623 ENDDO 624 DO i = 1, n 625 READ(kin,*) x1(i), y1(i) 626 y1(i) = y1(i)*1.e-3 ! mw -> w 627 ENDDO 628 CLOSE (kin) 629 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 630 CALL addpnt(x1,y1,kdata,n, 0.,0.) 631 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 632 CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) 633 CALL inter2(nw,wl,yg1,n,x1,y1,ierr) 634 635 IF (ierr .NE. 0) THEN 636 WRITE(*,*) ierr, fil 637 STOP 638 ENDIF 639 640 ! convert to photon.s-1.nm-1.cm-2 641 ! 5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8) 642 643 DO iw = 1, nw-1 644 f(iw) = yg1(iw)*wc(iw)*5.039e11 645 ! write(25,*) iw, wc(iw), f(iw) 646 ENDDO 647 648 end if 649 650 end subroutine rdsolarflux 651 652 !============================================================================== 653 654 subroutine addpnt ( x, y, ld, n, xnew, ynew ) 655 656 *-----------------------------------------------------------------------------* 657 *= PURPOSE: =* 658 *= Add a point <xnew,ynew> to a set of data pairs <x,y>. x must be in =* 659 *= ascending order =* 660 *-----------------------------------------------------------------------------* 661 *= PARAMETERS: =* 662 *= X - REAL vector of length LD, x-coordinates (IO)=* 663 *= Y - REAL vector of length LD, y-values (IO)=* 664 *= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=* 665 *= program =* 666 *= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=* 667 *= N < LD. On exit, N is incremented by 1. =* 668 *= XNEW - REAL, x-coordinate at which point is to be added (I)=* 669 *= YNEW - REAL, y-value of point to be added (I)=* 670 *-----------------------------------------------------------------------------* 671 672 IMPLICIT NONE 673 674 C calling parameters 675 676 INTEGER ld, n 677 REAL x(ld), y(ld) 678 REAL xnew, ynew 679 INTEGER ierr 680 681 C local variables 682 683 INTEGER insert 684 INTEGER i 685 686 C----------------------------------------------------------------------- 687 688 * initialize error flag 689 690 ierr = 0 691 692 * check n<ld to make sure x will hold another point 693 694 IF (n .GE. ld) THEN 695 WRITE(0,*) '>>> ERROR (ADDPNT) <<< Cannot expand array ' 696 WRITE(0,*) ' All elements used.' 697 STOP 698 ENDIF 699 700 insert = 1 701 i = 2 702 703 * check, whether x is already sorted. 704 * also, use this loop to find the point at which xnew needs to be inserted 705 * into vector x, if x is sorted. 706 707 10 CONTINUE 708 IF (i .LT. n) THEN 709 IF (x(i) .LT. x(i-1)) THEN 710 print*, x(i-1), x(i) 711 WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be '// 712 > 'in ascending order!' 713 STOP 714 ELSE 715 IF (xnew .GT. x(i)) insert = i + 1 716 ENDIF 717 i = i+1 718 GOTO 10 719 ENDIF 720 721 * if <xnew,ynew> needs to be appended at the end, just do so, 722 * otherwise, insert <xnew,ynew> at position INSERT 723 724 IF ( xnew .GT. x(n) ) THEN 725 726 x(n+1) = xnew 727 y(n+1) = ynew 728 729 ELSE 730 731 * shift all existing points one index up 732 733 DO i = n, insert, -1 734 x(i+1) = x(i) 735 y(i+1) = y(i) 736 ENDDO 737 738 * insert new point 739 740 x(insert) = xnew 741 y(insert) = ynew 742 743 ENDIF 744 745 * increase total number of elements in x, y 746 747 n = n+1 748 749 END 750 c 751 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 752 c 753 subroutine inter2(ng,xg,yg,n,x,y,ierr) 754 755 *-----------------------------------------------------------------------------* 756 *= PURPOSE: =* 757 *= Map input data given on single, discrete points onto a set of target =* 758 *= bins. =* 759 *= The original input data are given on single, discrete points of an =* 760 *= arbitrary grid and are being linearly interpolated onto a specified set =* 761 *= of target bins. In general, this is the case for most of the weighting =* 762 *= functions (action spectra, molecular cross section, and quantum yield =* 763 *= data), which have to be matched onto the specified wavelength intervals. =* 764 *= The average value in each target bin is found by averaging the trapezoi- =* 765 *= dal area underneath the input data curve (constructed by linearly connec-=* 766 *= ting the discrete input values). =* 767 *= Some caution should be used near the endpoints of the grids. If the =* 768 *= input data set does not span the range of the target grid, an error =* 769 *= message is printed and the execution is stopped, as extrapolation of the =* 770 *= data is not permitted. =* 771 *= If the input data does not encompass the target grid, use ADDPNT to =* 772 *= expand the input array. =* 773 *-----------------------------------------------------------------------------* 774 *= PARAMETERS: =* 775 *= NG - INTEGER, number of bins + 1 in the target grid (I)=* 776 *= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=* 777 *= as [XG(i),XG(i+1)] (i = 1..NG-1) =* 778 *= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=* 779 *= bin i (i = 1..NG-1) =* 780 *= N - INTEGER, number of points in input grid (I)=* 781 *= X - REAL, grid on which input data are defined (I)=* 782 *= Y - REAL, input y-data (I)=* 783 *-----------------------------------------------------------------------------* 784 785 IMPLICIT NONE 786 787 * input: 788 INTEGER ng, n 789 REAL x(n), y(n), xg(ng) 790 791 * output: 792 REAL yg(ng) 793 794 * local: 795 REAL area, xgl, xgu 796 REAL darea, slope 797 REAL a1, a2, b1, b2 798 INTEGER ngintv 799 INTEGER i, k, jstart 800 INTEGER ierr 801 *_______________________________________________________________________ 802 803 ierr = 0 804 805 * test for correct ordering of data, by increasing value of x 806 807 DO 10, i = 2, n 808 IF (x(i) .LE. x(i-1)) THEN 809 ierr = 1 810 WRITE(*,*)'data not sorted' 811 WRITE(*,*) x(i), x(i-1) 812 RETURN 813 ENDIF 814 10 CONTINUE 815 816 DO i = 2, ng 817 IF (xg(i) .LE. xg(i-1)) THEN 818 ierr = 2 819 WRITE(0,*) '>>> ERROR (inter2) <<< xg-grid not sorted!' 820 RETURN 821 ENDIF 822 ENDDO 823 824 * check for xg-values outside the x-range 825 826 IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN 827 WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span '// 828 > 'grid. ' 829 WRITE(0,*) ' Use ADDPNT to '// 830 > 'expand data and re-run.' 831 STOP 832 ENDIF 833 834 * find the integral of each grid interval and use this to 835 * calculate the average y value for the interval 836 * xgl and xgu are the lower and upper limits of the grid interval 837 838 jstart = 1 839 ngintv = ng - 1 840 DO 50, i = 1,ngintv 841 842 * initalize: 843 844 area = 0.0 845 xgl = xg(i) 846 xgu = xg(i+1) 847 848 * discard data before the first grid interval and after the 849 * last grid interval 850 * for internal grid intervals, start calculating area by interpolating 851 * between the last point which lies in the previous interval and the 852 * first point inside the current interval 853 854 k = jstart 855 IF (k .LE. n-1) THEN 856 857 * if both points are before the first grid, go to the next point 858 30 CONTINUE 859 IF (x(k+1) .LE. xgl) THEN 860 jstart = k - 1 861 k = k+1 862 IF (k .LE. n-1) GO TO 30 863 ENDIF 864 865 866 * if the last point is beyond the end of the grid, complete and go to the next 867 * grid 868 40 CONTINUE 869 IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN 870 871 jstart = k-1 872 873 * compute x-coordinates of increment 874 875 a1 = MAX(x(k),xgl) 876 a2 = MIN(x(k+1),xgu) 877 878 * if points coincide, contribution is zero 879 880 IF (x(k+1).EQ.x(k)) THEN 881 darea = 0.e0 882 ELSE 883 slope = (y(k+1) - y(k))/(x(k+1) - x(k)) 884 b1 = y(k) + slope*(a1 - x(k)) 885 b2 = y(k) + slope*(a2 - x(k)) 886 darea = (a2 - a1)*(b2 + b1)/2. 887 ENDIF 888 889 890 * find the area under the trapezoid from a1 to a2 891 892 area = area + darea 893 894 * go to next point 895 896 k = k+1 897 GO TO 40 898 899 ENDIF 900 901 ENDIF 902 903 * calculate the average y after summing the areas in the interval 904 yg(i) = area/(xgu - xgl) 905 906 50 CONTINUE 907 908 RETURN 909 END 910 911 !============================================================================== 912 913 subroutine rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2) 914 915 *-----------------------------------------------------------------------------* 916 *= PURPOSE: =* 917 *= Read and grid CO2 absorption cross-sections and photodissociation yield =* 918 *-----------------------------------------------------------------------------* 919 *= PARAMETERS: =* 920 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 921 *= wavelength grid =* 922 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 923 *= working wavelength grid =* 924 *= XSCO2 - REAL, molecular absoprtion cross section (cm^2) of CO2 at (O)=* 925 *= each specified wavelength =* 926 *-----------------------------------------------------------------------------* 927 928 use datafile_mod, only: datadir 929 930 implicit none 931 932 ! input 933 934 integer :: nw ! number of wavelength grid points 935 real, dimension(nw) :: wl ! lower and central wavelength for each interval 936 937 ! output 938 939 real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2) 940 real, dimension(nw) :: yieldco2 ! co2 photodissociation yield 941 942 ! local 943 944 integer, parameter :: kdata = 42000 945 real, parameter :: deltax = 1.e-4 946 real, dimension(kdata) :: x1, y1, y2, y3, xion, ion 947 real, dimension(nw) :: yg 948 real :: xl, xu 949 integer :: ierr, i, l, n, n1, n2, n3, n4 950 CHARACTER*100 fil 951 952 integer :: kin, kout ! input/ouput logical units 953 954 kin = 10 955 kout = 30 956 957 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 958 c 959 c CO2 absorption cross-sections 960 c 961 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 962 c 963 c 195K: huestis and berkowitz (2010) + Starck et al. (2006) 964 c + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation 965 c 966 c 295K: huestis and berkowitz (2010) + Starck et al. (2006) 967 c + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation 968 c 969 c 370K: huestis and berkowitz (2010) + Starck et al. (2006) 970 c + Lewis and Carver (1983) + extrapolation 971 c 972 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 973 974 n1 = 40769 975 n2 = 41586 976 n3 = 10110 977 978 ! 195K: 979 980 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_195k.txt' 981 print*, 'section efficace CO2 195K: ', fil 982 983 OPEN(UNIT=kin,FILE=fil,STATUS='old') 984 DO i = 1,11 985 read(kin,*) 986 END DO 987 988 DO i = 1, n1 989 READ(kin,*) x1(i), y1(i) 990 END DO 991 CLOSE (kin) 992 993 CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) 994 CALL addpnt(x1,y1,kdata,n1, 0.,0.) 995 CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) 996 CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.) 997 CALL inter2(nw,wl,yg,n1,x1,y1,ierr) 998 IF (ierr .NE. 0) THEN 999 WRITE(*,*) ierr, fil 1000 STOP 1001 ENDIF 1002 1003 DO l = 1, nw-1 1004 xsco2_195(l) = yg(l) 1005 END DO 1006 1007 ! 295K: 1008 1009 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_295k.txt' 1010 print*, 'section efficace CO2 295K: ', fil 1011 1012 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1013 DO i = 1,11 1014 read(kin,*) 1015 END DO 1016 1017 DO i = 1, n2 1018 READ(kin,*) x1(i), y1(i) 1019 END DO 1020 CLOSE (kin) 1021 1022 CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.) 1023 CALL addpnt(x1,y1,kdata,n2, 0.,0.) 1024 CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.) 1025 CALL addpnt(x1,y1,kdata,n2, 1.e+38,0.) 1026 CALL inter2(nw,wl,yg,n2,x1,y1,ierr) 1027 IF (ierr .NE. 0) THEN 1028 WRITE(*,*) ierr, fil 1029 STOP 1030 ENDIF 1031 1032 DO l = 1, nw-1 1033 xsco2_295(l) = yg(l) 1034 END DO 1035 1036 ! 370K: 1037 1038 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_370k.txt' 1039 print*, 'section efficace CO2 370K: ', fil 1040 1041 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1042 DO i = 1,11 1043 read(kin,*) 1044 END DO 1045 1046 DO i = 1, n3 1047 READ(kin,*) x1(i), y1(i) 1048 END DO 1049 CLOSE (kin) 1050 1051 CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.) 1052 CALL addpnt(x1,y1,kdata,n3, 0.,0.) 1053 CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.) 1054 CALL addpnt(x1,y1,kdata,n3, 1.e+38,0.) 1055 CALL inter2(nw,wl,yg,n3,x1,y1,ierr) 1056 IF (ierr .NE. 0) THEN 1057 WRITE(*,*) ierr, fil 1058 STOP 1059 ENDIF 1060 1061 DO l = 1, nw-1 1062 xsco2_370(l) = yg(l) 1063 END DO 1064 1065 ! photodissociation yield: 1066 1067 fil = trim(datadir)// 1068 $'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt' 1069 print*, 'photodissociation yield CO2: ', fil 1070 1071 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1072 1073 do i = 1,3 1074 read(kin,*) 1075 end do 1076 1077 n4 = 17 1078 do i = 1, n4 1079 read(kin,*) xl, xu, ion(i) 1080 xion(i) = (xl + xu)/2. 1081 ion(i) = max(ion(i), 0.) 1082 end do 1083 close(kin) 1084 1085 CALL addpnt(xion,ion,kdata,n4,xion(1)*(1.-deltax),0.) 1086 CALL addpnt(xion,ion,kdata,n4, 0.,0.) 1087 CALL addpnt(xion,ion,kdata,n4,xion(n4)*(1.+deltax),1.) 1088 CALL addpnt(xion,ion,kdata,n4, 1.e+38,1.) 1089 CALL inter2(nw,wl,yieldco2,n4,xion,ion,ierr) 1090 IF (ierr .NE. 0) THEN 1091 WRITE(*,*) ierr, fil 1092 STOP 1093 ENDIF 1094 1095 ! DO l = 1, nw-1 1096 ! write(kout,*) wl(l), xsco2_195(l), 1097 ! $ xsco2_295(l), 1098 ! $ xsco2_370(l), 1099 ! $ yieldco2(l) 1100 ! END DO 1101 1102 end subroutine rdxsco2 1103 1104 !============================================================================== 1105 1106 subroutine rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300, 1107 $ yieldo2) 1108 1109 *-----------------------------------------------------------------------------* 1110 *= PURPOSE: =* 1111 *= Read and grid O2 cross-sections and photodissociation yield =* 1112 *-----------------------------------------------------------------------------* 1113 *= PARAMETERS: =* 1114 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1115 *= wavelength grid =* 1116 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 1117 *= working wavelength grid =* 1118 *= XSO2 - REAL, molecular absorption cross section =* 1119 *-----------------------------------------------------------------------------* 1120 1121 use datafile_mod, only: datadir 1122 1123 implicit none 1124 1125 ! input 1126 1127 integer :: nw ! number of wavelength grid points 1128 real, dimension(nw) :: wl ! lower and central wavelength for each interval 1129 1130 ! output 1131 1132 real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 cross-sections (cm2) 1133 real, dimension(nw) :: yieldo2 ! o2 photodissociation yield 1134 1135 ! local 1136 1137 integer, parameter :: kdata = 18000 1138 real, parameter :: deltax = 1.e-4 1139 real, dimension(kdata) :: x1, y1, x2, y2, x3, y3, x4, y4 1140 real, dimension(kdata) :: xion, ion 1141 real :: factor, xl, xu, dummy 1142 integer :: i, ierr, n, n1, n2, n3, n4, nhead 1143 integer :: kin, kout ! input/output logical units 1144 character*100 fil 1145 1146 kin = 10 1147 kout = 30 1148 1149 ! read o2 cross section data 1150 1151 nhead = 22 1152 n = 17434 1153 1154 fil = trim(datadir)//'/cross_sections/o2_composite_2018_150K.txt' 1155 print*, 'section efficace O2 150K: ', fil 1156 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1157 1158 DO i = 1,nhead 1159 read(kin,*) 1160 END DO 1161 1162 n1 = n 1163 DO i = 1, n1 1164 READ(kin,*) x1(i), y1(i) 1165 END DO 1166 CLOSE (kin) 1167 1168 CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) 1169 CALL addpnt(x1,y1,kdata,n1, 0.,0.) 1170 CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) 1171 CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.) 1172 CALL inter2(nw,wl,xso2_150,n1,x1,y1,ierr) 1173 IF (ierr .NE. 0) THEN 1174 WRITE(*,*) ierr, fil 1175 STOP 1176 ENDIF 1177 1178 fil = trim(datadir)//'/cross_sections/o2_composite_2018_200K.txt' 1179 print*, 'section efficace O2 200K: ', fil 1180 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1181 1182 DO i = 1,nhead 1183 read(kin,*) 1184 END DO 1185 1186 n2 = n 1187 DO i = 1, n2 1188 READ(kin,*) x2(i), y2(i) 1189 END DO 1190 CLOSE (kin) 1191 1192 CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.) 1193 CALL addpnt(x2,y2,kdata,n2, 0.,0.) 1194 CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.) 1195 CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.) 1196 CALL inter2(nw,wl,xso2_200,n2,x2,y2,ierr) 1197 IF (ierr .NE. 0) THEN 1198 WRITE(*,*) ierr, fil 1199 STOP 1200 ENDIF 1201 1202 fil = trim(datadir)//'/cross_sections/o2_composite_2018_250K.txt' 1203 print*, 'section efficace O2 250K: ', fil 1204 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1205 1206 DO i = 1,nhead 1207 read(kin,*) 1208 END DO 1209 1210 n3 = n 1211 DO i = 1, n3 1212 READ(kin,*) x3(i), y3(i) 1213 END DO 1214 CLOSE (kin) 1215 1216 CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.) 1217 CALL addpnt(x3,y3,kdata,n3, 0.,0.) 1218 CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.) 1219 CALL addpnt(x3,y3,kdata,n3, 1.e+38,0.) 1220 CALL inter2(nw,wl,xso2_250,n3,x3,y3,ierr) 1221 IF (ierr .NE. 0) THEN 1222 WRITE(*,*) ierr, fil 1223 STOP 1224 ENDIF 1225 1226 fil = trim(datadir)//'/cross_sections/o2_composite_2018_300K.txt' 1227 print*, 'section efficace O2 300K: ', fil 1228 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1229 1230 DO i = 1,nhead 1231 read(kin,*) 1232 END DO 1233 1234 n4 = n 1235 DO i = 1, n4 1236 READ(kin,*) x4(i), y4(i) 1237 END DO 1238 CLOSE (kin) 1239 1240 CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.) 1241 CALL addpnt(x4,y4,kdata,n4, 0.,0.) 1242 CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.) 1243 CALL addpnt(x4,y4,kdata,n4, 1.e+38,0.) 1244 CALL inter2(nw,wl,xso2_300,n4,x4,y4,ierr) 1245 IF (ierr .NE. 0) THEN 1246 WRITE(*,*) ierr, fil 1247 STOP 1248 ENDIF 1249 1250 ! photodissociation yield 1251 1252 fil = trim(datadir)// 1253 $ '/cross_sections/efdis_co2-o2_schunkandnagy2000.txt' 1254 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1255 1256 do i = 1,11 1257 read(kin,*) 1258 end do 1259 1260 n = 9 1261 DO i = 1, n 1262 READ(kin,*) xl, xu, dummy, ion(i) 1263 xion(i) = (xl + xu)/2. 1264 ion(i) = max(ion(i), 0.) 1265 END DO 1266 CLOSE (kin) 1267 1268 CALL addpnt(xion,ion,kdata,n,xion(1)*(1.-deltax),0.) 1269 CALL addpnt(xion,ion,kdata,n, 0.,0.) 1270 CALL addpnt(xion,ion,kdata,n,xion(n)*(1.+deltax),1.) 1271 CALL addpnt(xion,ion,kdata,n, 1.e+38,1.) 1272 CALL inter2(nw,wl,yieldo2,n,xion,ion,ierr) 1273 IF (ierr .NE. 0) THEN 1274 WRITE(*,*) ierr, fil 1275 STOP 1276 ENDIF 1277 1278 end subroutine rdxso2 1279 1280 !============================================================================== 1281 1282 subroutine rdxso3(nw,wl,xso3_218,xso3_298) 1283 1284 *-----------------------------------------------------------------------------* 1285 *= PURPOSE: =* 1286 *= Read ozone molecular absorption cross section. Re-grid data to match =* 1287 *= specified wavelength working grid. =* 1288 *-----------------------------------------------------------------------------* 1289 *= PARAMETERS: =* 1290 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1291 *= wavelength grid =* 1292 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 1293 *= working wavelength grid =* 1294 *= XSO3_218 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=* 1295 *= each specified wavelength (JPL 2006) 218 K =* 1296 *= XSO3_298 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=* 1297 *= each specified wavelength (JPL 2006) 298 K =* 1298 *-----------------------------------------------------------------------------* 1299 1300 use datafile_mod, only: datadir 1301 1302 implicit none 1303 1304 ! input 1305 1306 integer :: nw ! number of wavelength grid points 1307 real, dimension(nw) :: wl ! lower and central wavelength for each interval 1308 1309 ! output 1310 1311 real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2) 1312 1313 ! local 1314 1315 integer, parameter :: kdata = 200 1316 real, parameter :: deltax = 1.e-4 1317 real, dimension(kdata) :: x1, x2, y1, y2 1318 real, dimension(nw) :: yg 1319 real :: a1, a2 1320 1321 integer :: i, ierr, iw, n, n1, n2 1322 integer :: kin, kout ! input/output logical units 1323 1324 character*100 fil 1325 1326 ! JPL 2006 218 K 1327 1328 fil = trim(datadir)// 1329 $ '/cross_sections/o3_cross-sections_jpl_2006_218K.txt' 1330 print*, 'section efficace O3 218K: ', fil 1331 1332 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1333 n1 = 167 1334 DO i = 1, n1 1335 READ(kin,*) a1, a2, y1(i) 1336 x1(i) = (a1+a2)/2. 1337 END DO 1338 CLOSE (kin) 1339 1340 CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) 1341 CALL addpnt(x1,y1,kdata,n1, 0.,0.) 1342 CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) 1343 CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.) 1344 CALL inter2(nw,wl,yg,n1,x1,y1,ierr) 1345 IF (ierr .NE. 0) THEN 1346 WRITE(*,*) ierr, fil 1347 STOP 1348 ENDIF 1349 c 1350 DO iw = 1, nw-1 1351 xso3_218(iw) = yg(iw) 1352 END DO 1353 1354 ! JPL 2006 298 K 1355 1356 fil = trim(datadir)// 1357 $ '/cross_sections/o3_cross-sections_jpl_2006_298K.txt' 1358 print*, 'section efficace O3 298K: ', fil 1359 1360 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1361 n2 = 167 1362 DO i = 1, n2 1363 READ(kin,*) a1, a2, y2(i) 1364 x2(i) = (a1+a2)/2. 1365 END DO 1366 CLOSE (kin) 1367 1368 CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.) 1369 CALL addpnt(x2,y2,kdata,n2, 0.,0.) 1370 CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.) 1371 CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.) 1372 CALL inter2(nw,wl,yg,n2,x2,y2,ierr) 1373 IF (ierr .NE. 0) THEN 1374 WRITE(*,*) ierr, fil 1375 STOP 1376 ENDIF 1377 c 1378 DO iw = 1, nw-1 1379 xso3_298(iw) = yg(iw) 1380 END DO 1381 1382 end subroutine rdxso3 1383 1384 !============================================================================== 1385 1386 subroutine rdxsh2o(nw, wl, yg) 1387 1388 *-----------------------------------------------------------------------------* 1389 *= PURPOSE: =* 1390 *= Read H2O molecular absorption cross section. Re-grid data to match =* 1391 *= specified wavelength working grid. =* 1392 *-----------------------------------------------------------------------------* 1393 *= PARAMETERS: =* 1394 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1395 *= wavelength grid =* 1396 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 1397 *= working wavelength grid =* 1398 *= YG - REAL, molecular absoprtion cross section (cm^2) of H2O at (O)=* 1399 *= each specified wavelength =* 1400 *-----------------------------------------------------------------------------* 1401 1402 use datafile_mod, only: datadir 1403 1404 IMPLICIT NONE 1405 1406 ! input 1407 1408 integer :: nw ! number of wavelength grid points 1409 real, dimension(nw) :: wl ! lower and central wavelength for each interval 1410 1411 ! output 1412 1413 real, dimension(nw) :: yg ! h2o cross-sections (cm2) 1414 1415 ! local 1416 1417 integer, parameter :: kdata = 500 1418 real, parameter :: deltax = 1.e-4 1419 REAL x1(kdata) 1420 REAL y1(kdata) 1421 INTEGER ierr 1422 INTEGER i, n 1423 CHARACTER*100 fil 1424 integer :: kin, kout ! input/output logical units 1425 1426 kin = 10 1427 1428 fil = trim(datadir)//'/cross_sections/h2o_composite_250K.txt' 1429 print*, 'section efficace H2O: ', fil 1430 1431 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1432 1433 DO i = 1,26 1434 read(kin,*) 1435 END DO 1436 1437 n = 420 1438 DO i = 1, n 1439 READ(kin,*) x1(i), y1(i) 1440 END DO 1441 CLOSE (kin) 1442 1443 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1444 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1445 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1446 CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) 1447 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1448 IF (ierr .NE. 0) THEN 1449 WRITE(*,*) ierr, fil 1450 STOP 1451 ENDIF 1452 1453 end subroutine rdxsh2o 1454 1455 !============================================================================== 1456 1457 subroutine rdxsh2o2(nw, wl, xsh2o2) 1458 1459 *-----------------------------------------------------------------------------* 1460 *= PURPOSE: =* 1461 *= Read and grid H2O2 cross-sections 1462 *= H2O2 + hv -> 2 OH =* 1463 *= Cross section: Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508 =* 1464 *= from 125 to 185 nm, then JPL97 from 190 to 350 nm. =* 1465 *-----------------------------------------------------------------------------* 1466 *= PARAMETERS: =* 1467 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1468 *= wavelength grid =* 1469 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 1470 *= working wavelength grid =* 1471 *-----------------------------------------------------------------------------* 1472 1473 use datafile_mod, only: datadir 1474 1475 implicit none 1476 1477 ! input 1478 1479 integer :: nw ! number of wavelength grid points 1480 real, dimension(nw) :: wl ! lower wavelength for each interval 1481 1482 ! output 1483 1484 real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2) 1485 1486 ! local 1487 1488 real, parameter :: deltax = 1.e-4 1489 integer, parameter :: kdata = 100 1490 real, dimension(kdata) :: x1, y1 1491 real, dimension(nw) :: yg 1492 integer :: i, ierr, iw, n, idum 1493 integer :: kin, kout ! input/output logical units 1494 character*100 fil 1495 1496 kin = 10 1497 1498 ! read cross-sections 1499 1500 fil = trim(datadir)//'/cross_sections/h2o2_composite.txt' 1501 print*, 'section efficace H2O2: ', fil 1502 1503 OPEN(kin,FILE=fil,STATUS='OLD') 1504 READ(kin,*) idum,n 1505 DO i = 1, idum-2 1506 READ(kin,*) 1507 ENDDO 1508 DO i = 1, n 1509 READ(kin,*) x1(i), y1(i) 1510 ENDDO 1511 CLOSE (kin) 1512 1513 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1514 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1515 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1516 CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) 1517 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1518 IF (ierr .NE. 0) THEN 1519 WRITE(*,*) ierr, fil 1520 STOP 1521 ENDIF 1522 1523 DO iw = 1, nw - 1 1524 xsh2o2(iw) = yg(iw) 1525 END DO 1526 1527 end subroutine rdxsh2o2 1528 1529 !============================================================================== 1530 1531 subroutine rdxsho2(nw, wl, yg) 1532 1533 *-----------------------------------------------------------------------------* 1534 *= PURPOSE: =* 1535 *= Read ho2 cross-sections =* 1536 *= JPL 2006 recommendation =* 1537 *-----------------------------------------------------------------------------* 1538 *= PARAMETERS: =* 1539 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1540 *= wavelength grid =* 1541 *-----------------------------------------------------------------------------* 1542 1543 use datafile_mod, only: datadir 1544 1545 IMPLICIT NONE 1546 1547 ! input 1548 1549 integer :: nw ! number of wavelength grid points 1550 real, dimension(nw) :: wl ! lower wavelength for each interval 1551 1552 ! output 1553 1554 real, dimension(nw) :: yg ! ho2 cross-sections (cm2) 1555 1556 ! local 1557 1558 real, parameter :: deltax = 1.e-4 1559 integer, parameter :: kdata = 100 1560 real, dimension(kdata) :: x1, y1 1561 integer :: i, n, ierr 1562 character*100 fil 1563 integer :: kin, kout ! input/output logical units 1564 1565 kin = 10 1566 1567 **** cross sections from Sander et al. [2003] 1568 1569 fil = trim(datadir)//'/cross_sections/ho2_jpl2003.txt' 1570 print*, 'section efficace HO2: ', fil 1571 1572 OPEN(kin,FILE=fil,STATUS='OLD') 1573 READ(kin,*) n 1574 DO i = 1, n 1575 READ(kin,*) x1(i), y1(i) 1576 ENDDO 1577 CLOSE(kin) 1578 1579 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1580 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1581 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1582 CALL addpnt(x1,y1,kdata,n, 1E38,0.) 1583 1584 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1585 1586 IF (ierr .NE. 0) THEN 1587 WRITE(*,*) ierr, fil 1588 STOP 1589 ENDIF 1590 1591 end subroutine rdxsho2 1592 1593 !============================================================================== 1594 1595 subroutine rdxsh2(nw, wl, wc, yg, yieldh2) 1596 1597 *-----------------------------------------------------------------------------* 1598 *= PURPOSE: =* 1599 *= Read h2 cross-sections and photodissociation yield =* 1600 *-----------------------------------------------------------------------------* 1601 *= PARAMETERS: =* 1602 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1603 *= wavelength grid =* 1604 *-----------------------------------------------------------------------------* 1605 1606 use datafile_mod, only: datadir 1607 1608 implicit none 1609 1610 ! input 1611 1612 integer :: nw ! number of wavelength grid points 1613 real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval 1614 1615 ! output 1616 1617 real, dimension(nw) :: yg ! h2 cross-sections (cm2) 1618 real, dimension(nw) :: yieldh2 ! photodissociation yield 1619 1620 ! local 1621 1622 integer, parameter :: kdata = 1000 1623 real, parameter :: deltax = 1.e-4 1624 real, dimension(kdata) :: x1, y1, x2, y2 1625 real :: xl, xu 1626 integer :: i, iw, n, ierr 1627 integer :: kin, kout ! input/output logical units 1628 character*100 fil 1629 1630 kin = 10 1631 1632 ! h2 cross sections 1633 1634 fil = trim(datadir)//'/cross_sections/h2secef.txt' 1635 print*, 'section efficace H2: ', fil 1636 1637 OPEN(kin,FILE=fil,STATUS='OLD') 1638 1639 n = 792 1640 read(kin,*) ! avoid first line with wavelength = 0. 1641 DO i = 1, n 1642 READ(kin,*) x1(i), y1(i) 1643 ENDDO 1644 CLOSE(kin) 1645 1646 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1647 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1648 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1649 CALL addpnt(x1,y1,kdata,n, 1E38,0.) 1650 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1651 1652 IF (ierr .NE. 0) THEN 1653 WRITE(*,*) ierr, fil 1654 STOP 1655 ENDIF 1656 1657 ! photodissociation yield 1658 1659 fil = trim(datadir)//'/cross_sections/h2_ionef_schunknagy2000.txt' 1660 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1661 1662 n = 19 1663 read(kin,*) 1664 DO i = 1, n 1665 READ(kin,*) xl, xu, y2(i) 1666 x2(i) = (xl + xu)/2. 1667 y2(i) = max(1. - y2(i),0.) 1668 END DO 1669 CLOSE (kin) 1670 1671 CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) 1672 CALL addpnt(x2,y2,kdata,n, 0.,0.) 1673 CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.) 1674 CALL addpnt(x2,y2,kdata,n, 1.e+38,1.) 1675 CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr) 1676 IF (ierr .NE. 0) THEN 1677 WRITE(*,*) ierr, fil 1678 STOP 1679 ENDIF 1680 1681 end subroutine rdxsh2 1682 1683 !============================================================================== 1684 1685 subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294, 1686 $ yldno2_248, yldno2_298) 1687 1688 *-----------------------------------------------------------------------------* 1689 *= PURPOSE: =* 1690 *= read and grid cross section + quantum yield for NO2 =* 1691 *= photolysis =* 1692 *= Jenouvrier et al., 1996 200-238 nm 1693 *= Vandaele et al., 1998 238-666 nm 220K and 294K 1694 *= quantum yield from jpl 2006 1695 *-----------------------------------------------------------------------------* 1696 *= PARAMETERS: =* 1697 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1698 *= wavelength grid =* 1699 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 1700 *= working wavelength grid =* 1701 *= SQ - REAL, cross section x quantum yield (cm^2) for each (O)=* 1702 *= photolysis reaction defined, at each defined wavelength and =* 1703 *= at each defined altitude level =* 1704 *-----------------------------------------------------------------------------* 1705 1706 use datafile_mod, only: datadir 1707 1708 implicit none 1709 1710 ! input 1711 1712 integer :: nw ! number of wavelength grid points 1713 real, dimension(nw) :: wl ! lower and central wavelength for each interval 1714 1715 ! output 1716 1717 real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2) 1718 real, dimension(nw) :: yldno2_248, yldno2_298 ! quantum yields at 248-298 k 1719 1720 ! local 1721 1722 integer, parameter :: kdata = 28000 1723 real, parameter :: deltax = 1.e-4 1724 real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5 1725 real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5 1726 real :: dum, qy 1727 integer :: i, iw, n, n1, n2, n3, n4, n5, ierr 1728 character*100 fil 1729 integer :: kin, kout ! input/output logical units 1730 1731 kin = 10 1732 1733 !*************** NO2 photodissociation 1734 1735 * Jenouvrier 1996 + Vandaele 1998 (JPL 2006) 1736 1737 fil = trim(datadir)//'/cross_sections/no2_xs_jenouvrier.txt' 1738 print*, 'section efficace NO2: ', fil 1739 1740 OPEN(UNIT=kin,FILE=fil,status='old') 1741 DO i = 1, 3 1742 READ(kin,*) 1743 ENDDO 1744 n1 = 10001 1745 DO i = 1, n1 1746 READ(kin,*) x1(i), y1(i) 1747 end do 1748 1749 CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.) 1750 CALL addpnt(x1,y1,kdata,n1, 0., 0.) 1751 CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) 1752 CALL addpnt(x1,y1,kdata,n1, 1.e+38, 0.) 1753 CALL inter2(nw,wl,yg1,n1,x1,y1,ierr) 1754 1755 fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_294K.txt' 1756 print*, 'section efficace NO2: ', fil 1757 1758 OPEN(UNIT=kin,FILE=fil,status='old') 1759 DO i = 1, 3 1760 READ(kin,*) 1761 ENDDO 1762 n2 = 27993 1763 DO i = 1, n2 1764 READ(kin,*) x2(i), y2(i) 1765 end do 1766 1767 CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.) 1768 CALL addpnt(x2,y2,kdata,n2, 0., 0.) 1769 CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.) 1770 CALL addpnt(x2,y2,kdata,n2, 1.e+38, 0.) 1771 CALL inter2(nw,wl,yg2,n2,x2,y2,ierr) 1772 1773 fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_220K.txt' 1774 print*, 'section efficace NO2: ', fil 1775 1776 OPEN(UNIT=kin,FILE=fil,status='old') 1777 DO i = 1, 3 1778 READ(kin,*) 1779 ENDDO 1780 n3 = 27993 1781 DO i = 1, n3 1782 READ(kin,*) x3(i), y3(i) 1783 end do 1784 1785 CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.) 1786 CALL addpnt(x3,y3,kdata,n3, 0., 0.) 1787 CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.) 1788 CALL addpnt(x3,y3,kdata,n3, 1.e+38, 0.) 1789 CALL inter2(nw,wl,yg3,n3,x3,y3,ierr) 1790 1791 do iw = 1, nw - 1 1792 xsno2(iw) = yg1(iw) 1793 xsno2_294(iw) = yg2(iw) 1794 xsno2_220(iw) = yg3(iw) 1795 end do 1796 1797 ! photodissociation efficiency from jpl 2006 1798 1799 fil = trim(datadir)//'/cross_sections/no2_yield_jpl2006.txt' 1800 print*, 'quantum yield NO2: ', fil 1801 1802 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1803 DO i = 1, 5 1804 READ(kin,*) 1805 ENDDO 1806 n = 25 1807 n4 = n 1808 n5 = n 1809 DO i = 1, n 1810 READ(kin,*) x4(i), y4(i), y5(i) 1811 x5(i) = x4(i) 1812 ENDDO 1813 CLOSE(kin) 1814 1815 CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1)) 1816 CALL addpnt(x4,y4,kdata,n4, 0.,y4(1)) 1817 CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax), 0.) 1818 CALL addpnt(x4,y4,kdata,n4, 1.e+38, 0.) 1819 CALL inter2(nw,wl,yg4,n4,x4,y4,ierr) 1820 IF (ierr .NE. 0) THEN 1821 WRITE(*,*) ierr, fil 1822 STOP 1823 ENDIF 1824 1825 CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1)) 1826 CALL addpnt(x5,y5,kdata,n5, 0.,y5(1)) 1827 CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax), 0.) 1828 CALL addpnt(x5,y5,kdata,n5, 1.e+38, 0.) 1829 CALL inter2(nw,wl,yg5,n5,x5,y5,ierr) 1830 IF (ierr .NE. 0) THEN 1831 WRITE(*,*) ierr, fil 1832 STOP 1833 ENDIF 1834 1835 do iw = 1, nw - 1 1836 yldno2_298(iw) = yg4(iw) 1837 yldno2_248(iw) = yg5(iw) 1838 end do 1839 1840 end subroutine rdxsno2 1841 1842 !============================================================================== 1843 1844 subroutine rdxsno(nw, wl, yg, yieldno) 1845 1846 *-----------------------------------------------------------------------------* 1847 *= PURPOSE: =* 1848 *= Read NO cross-sections and photodissociation efficiency =* 1849 *= Lida et al 1986 (provided by Francisco Gonzalez-Galindo) =* 1850 *-----------------------------------------------------------------------------* 1851 *= PARAMETERS: =* 1852 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1853 *= wavelength grid =* 1854 *-----------------------------------------------------------------------------* 1855 1856 use datafile_mod, only: datadir 1857 1858 implicit none 1859 1860 ! input 1861 1862 integer :: nw ! number of wavelength grid points 1863 real, dimension(nw) :: wl ! lower wavelength for each interval 1864 1865 ! output 1866 1867 real, dimension(nw) :: yg ! no cross-sections (cm2) 1868 real, dimension(nw) :: yieldno ! no photodissociation efficiency 1869 1870 ! local 1871 1872 integer, parameter :: kdata = 110 1873 real, parameter :: deltax = 1.e-4 1874 real, dimension(kdata) :: x1, y1, x2, y2 1875 integer :: i, iw, n, ierr 1876 character*100 fil 1877 integer :: kin, kout ! input/output logical units 1878 1879 kin = 10 1880 1881 ! no cross-sections 1882 1883 fil = trim(datadir)//'/cross_sections/no_xs_francisco.txt' 1884 print*, 'section efficace NO: ', fil 1885 OPEN(kin,FILE=fil,STATUS='OLD') 1886 1887 n = 99 1888 DO i = 1, n 1889 READ(kin,*) x1(i), y1(i) 1890 ENDDO 1891 CLOSE(kin) 1892 1893 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1894 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1895 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1896 CALL addpnt(x1,y1,kdata,n, 1E38,0.) 1897 1898 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1899 IF (ierr .NE. 0) THEN 1900 WRITE(*,*) ierr, fil 1901 STOP 1902 ENDIF 1903 1904 ! photodissociation yield 1905 1906 fil = trim(datadir)//'/cross_sections/noefdis.txt' 1907 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1908 1909 n = 33 1910 DO i = 1, n 1911 READ(kin,*) x2(n-i+1), y2(n-i+1) 1912 END DO 1913 CLOSE (kin) 1914 1915 CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) 1916 CALL addpnt(x2,y2,kdata,n, 0.,0.) 1917 CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.) 1918 CALL addpnt(x2,y2,kdata,n, 1.e+38,1.) 1919 CALL inter2(nw,wl,yieldno,n,x2,y2,ierr) 1920 IF (ierr .NE. 0) THEN 1921 WRITE(*,*) ierr, fil 1922 STOP 1923 ENDIF 1924 1925 end subroutine rdxsno 1926 1927 !============================================================================== 1928 1929 subroutine rdxsn2(nw, wl, yg, yieldn2) 1930 1931 *-----------------------------------------------------------------------------* 1932 *= PURPOSE: =* 1933 *= Read n2 cross-sections and photodissociation yield =* 1934 *-----------------------------------------------------------------------------* 1935 *= PARAMETERS: =* 1936 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 1937 *= wavelength grid =* 1938 *-----------------------------------------------------------------------------* 1939 1940 use datafile_mod, only: datadir 1941 1942 implicit none 1943 1944 ! input 1945 1946 integer :: nw ! number of wavelength grid points 1947 real, dimension(nw) :: wl ! lower wavelength for each interval 1948 1949 ! output 1950 1951 real, dimension(nw) :: yg ! n2 cross-sections (cm2) 1952 real, dimension(nw) :: yieldn2 ! n2 photodissociation yield 1953 1954 ! local 1955 1956 integer, parameter :: kdata = 1100 1957 real, parameter :: deltax = 1.e-4 1958 real, dimension(kdata) :: x1, y1, x2, y2 1959 real :: xl, xu 1960 integer :: i, iw, n, ierr 1961 integer :: kin, kout ! input/output logical units 1962 character*100 fil 1963 1964 kin = 10 1965 1966 ! n2 cross sections 1967 1968 fil = trim(datadir)//'/cross_sections/n2secef_01nm.txt' 1969 print*, 'section efficace N2: ', fil 1970 OPEN(kin,FILE=fil,STATUS='OLD') 1971 1972 n = 1020 1973 DO i = 1, n 1974 READ(kin,*) x1(i), y1(i) 1975 ENDDO 1976 CLOSE(kin) 1977 1978 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1979 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1980 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1981 CALL addpnt(x1,y1,kdata,n, 1E38,0.) 1982 1983 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1984 1985 IF (ierr .NE. 0) THEN 1986 WRITE(*,*) ierr, fil 1987 STOP 1988 ENDIF 1989 1990 ! photodissociation yield 1991 1992 fil = trim(datadir)//'/cross_sections/n2_ionef_schunknagy2000.txt' 1993 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1994 1995 n = 19 1996 read(kin,*) 1997 DO i = 1, n 1998 READ(kin,*) xl, xu, y2(i) 1999 x2(i) = (xl + xu)/2. 2000 y2(i) = 1. - y2(i) 2001 END DO 2002 CLOSE (kin) 2003 2004 CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) 2005 CALL addpnt(x2,y2,kdata,n, 0.,0.) 2006 CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.) 2007 CALL addpnt(x2,y2,kdata,n, 1.e+38,1.) 2008 CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr) 2009 IF (ierr .NE. 0) THEN 2010 WRITE(*,*) ierr, fil 2011 STOP 2012 ENDIF 2013 2014 end subroutine rdxsn2 2015 2016 !============================================================================== 2017 2018 subroutine setalb(nw,wl,albedo) 2019 2020 *-----------------------------------------------------------------------------* 2021 *= PURPOSE: =* 2022 *= Set the albedo of the surface. The albedo is assumed to be Lambertian, =* 2023 *= i.e., the reflected light is isotropic, and idependt of the direction =* 2024 *= of incidence of light. Albedo can be chosen to be wavelength dependent. =* 2025 *-----------------------------------------------------------------------------* 2026 *= PARAMETERS: =* 2027 *= NW - INTEGER, number of specified intervals + 1 in working (I)=* 2028 *= wavelength grid =* 2029 *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* 2030 *= working wavelength grid =* 2031 *= ALBEDO - REAL, surface albedo at each specified wavelength (O)=* 2032 *-----------------------------------------------------------------------------* 2033 2034 implicit none 2035 2036 ! input: (wavelength working grid data) 2037 2038 INTEGER nw 2039 REAL wl(nw) 2040 2041 ! output: 2042 2043 REAL albedo(nw) 2044 2045 ! local: 2046 2047 INTEGER iw 2048 REAL alb 2049 2050 ! 0.015: mean value from clancy et al., icarus, 49-63, 1999. 2051 2052 alb = 0.015 2053 2054 do iw = 1, nw - 1 2055 albedo(iw) = alb 2056 end do 2057 2058 return 2059 end 2060 236 contains 237 2061 238 !============================================================================== 2062 239 … … 2092 269 integer :: ilev, iw 2093 270 2094 ccompute column increments271 ! compute column increments 2095 272 2096 273 do ilev = 1, nlev-1 … … 2100 277 colinc(nlev) = avo*0.1*press(nlev)*100./(zmmean(nlev)*g) 2101 278 2102 ! do ilev = nlev,1,-12103 ! print*, ilev, press(ilev), temp(ilev), zmmean(ilev),2104 ! $ colinc(ilev)2105 ! end do2106 2107 279 do iw = 1, nw - 1 2108 280 2109 ccalcul de section efficace Rayleigh: voir Atreya and Gu, JGR,2110 c13133-13145, 1994.2111 c 2112 crho : normal depolarization ratio = 0.0774 for CO22113 cdf : depolarization factor2114 calpha : average dielectric polarizability = 2.911e-24 cm3 for CO2281 ! calcul de section efficace Rayleigh: voir Atreya and Gu, JGR, 282 ! 13133-13145, 1994. 283 ! 284 ! rho : normal depolarization ratio = 0.0774 for CO2 285 ! df : depolarization factor 286 ! alpha : average dielectric polarizability = 2.911e-24 cm3 for CO2 2115 287 2116 288 pi = acos(-1.0) … … 2127 299 end do 2128 300 2129 return 2130 end 301 end subroutine setair 2131 302 2132 303 !============================================================================== … … 2136 307 $ colinc, rm, dt, sj) 2137 308 2138 *-----------------------------------------------------------------------------*2139 *= PURPOSE: =*2140 *= Set up the CO2 temperature-dependent cross-sections and optical depth =*2141 *-----------------------------------------------------------------------------*309 !-----------------------------------------------------------------------------* 310 != PURPOSE: =* 311 != Set up the CO2 temperature-dependent cross-sections and optical depth =* 312 !-----------------------------------------------------------------------------* 2142 313 2143 314 implicit none … … 2216 387 end do 2217 388 2218 ! do l = 1,nw-12219 ! write(35,*) wc(l), sj(1,l,j_co2_o1d),2220 ! $ sj(1,l,j_co2_o), tlay(1)2221 ! end do2222 2223 389 end subroutine setco2 2224 390 … … 2229 395 $ j_o2_o1d, colinc, rm, dt, sj) 2230 396 2231 *-----------------------------------------------------------------------------*2232 *= PURPOSE: =*2233 *= Set up the O2 temperature-dependent cross-sections and optical depth =*2234 *-----------------------------------------------------------------------------*397 !-----------------------------------------------------------------------------* 398 != PURPOSE: =* 399 != Set up the O2 temperature-dependent cross-sections and optical depth =* 400 !-----------------------------------------------------------------------------* 2235 401 2236 402 implicit none … … 2316 482 $ colinc, rm, dt, sj) 2317 483 2318 *-----------------------------------------------------------------------------*2319 *= PURPOSE: =*2320 *= Set up the O3 temperature dependent cross-sections and optical depth =*2321 *-----------------------------------------------------------------------------*484 !-----------------------------------------------------------------------------* 485 != PURPOSE: =* 486 != Set up the O3 temperature dependent cross-sections and optical depth =* 487 !-----------------------------------------------------------------------------* 2322 488 2323 489 implicit none … … 2398 564 $ colinc, rm, dt, sj) 2399 565 2400 *-----------------------------------------------------------------------------*2401 *= PURPOSE: =*2402 *= Set up the h2o2 temperature dependent cross-sections and optical depth =*2403 *-----------------------------------------------------------------------------*566 !-----------------------------------------------------------------------------* 567 != PURPOSE: =* 568 != Set up the h2o2 temperature dependent cross-sections and optical depth =* 569 !-----------------------------------------------------------------------------* 2404 570 2405 571 implicit none … … 2477 643 $ colinc, rm, dt, sj) 2478 644 2479 *-----------------------------------------------------------------------------*2480 *= PURPOSE: =*2481 *= Set up the no2 temperature-dependent cross-sections and optical depth =*2482 *-----------------------------------------------------------------------------*645 !-----------------------------------------------------------------------------* 646 != PURPOSE: =* 647 != Set up the no2 temperature-dependent cross-sections and optical depth =* 648 !-----------------------------------------------------------------------------* 2483 649 2484 650 implicit none … … 2537 703 end do 2538 704 2539 ! do iw = 1, nw-12540 ! write(35,*) wc(iw), sj(1,iw,j_no2), tlay(1)2541 ! end do2542 ! stop2543 2544 705 end subroutine setno2 2545 706 … … 2548 709 subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer) 2549 710 2550 *-----------------------------------------------------------------------------*2551 *= PURPOSE: =*2552 *= Set aerosol properties for each specified altitude layer. Properties =*2553 *= may be wavelength dependent. =*2554 *-----------------------------------------------------------------------------*711 !-----------------------------------------------------------------------------* 712 != PURPOSE: =* 713 != Set aerosol properties for each specified altitude layer. Properties =* 714 != may be wavelength dependent. =* 715 !-----------------------------------------------------------------------------* 2555 716 2556 717 implicit none … … 2593 754 end do 2594 755 2595 ! print*, 'tau = ', tau2596 ! print*, 'q0 = ', q02597 756 q0 = tau/tautot 2598 757 do ilay = 1, nlayer-1 … … 2603 762 end do 2604 763 2605 ! do ilay = nlayer,1,-12606 ! print*, alt(ilay), q0*exp(gamma*(1. - exp(alt(ilay)/scaleh))),2607 ! $ dtaer(ilay,1)2608 ! end do2609 2610 764 end subroutine setaer 2611 765 … … 2614 768 subroutine setcld(nlayer,nw,dtcld,omcld,gcld) 2615 769 2616 *-----------------------------------------------------------------------------*2617 *= PURPOSE: =*2618 *= Set cloud properties for each specified altitude layer. Properties =*2619 *= may be wavelength dependent. =*2620 *-----------------------------------------------------------------------------*770 !-----------------------------------------------------------------------------* 771 != PURPOSE: =* 772 != Set cloud properties for each specified altitude layer. Properties =* 773 != may be wavelength dependent. =* 774 !-----------------------------------------------------------------------------* 2621 775 2622 776 implicit none … … 2637 791 integer :: ilay, iw 2638 792 2639 cdtcld : optical depth2640 comcld : single scattering albedo2641 cgcld : asymmetry factor793 ! dtcld : optical depth 794 ! omcld : single scattering albedo 795 ! gcld : asymmetry factor 2642 796 2643 797 do ilay = 1, nlayer - 1 … … 2649 803 end do 2650 804 2651 end 805 end subroutine setcld 2652 806 2653 807 !============================================================================== … … 2655 809 subroutine sphers(nlev, z, zen, dsdh, nid) 2656 810 2657 *-----------------------------------------------------------------------------*2658 *= PURPOSE: =*2659 *= Calculate slant path over vertical depth ds/dh in spherical geometry. =*2660 *= Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =*2661 *= for computing the radiation field available for photolysis and heating =*2662 *= at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =*2663 *-----------------------------------------------------------------------------*2664 *= PARAMETERS: =*2665 *= NZ - INTEGER, number of specified altitude levels in the working (I)=*2666 *= grid =*2667 *= Z - REAL, specified altitude working grid (km) (I)=*2668 *= ZEN - REAL, solar zenith angle (degrees) (I)=*2669 *= DSDH - REAL, slant path of direct beam through each layer crossed (O)=*2670 *= when travelling from the top of the atmosphere to layer i; =*2671 *= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*2672 *= NID - INTEGER, number of layers crossed by the direct beam when (O)=*2673 *= travelling from the top of the atmosphere to layer i; =*2674 *= NID(i), i = 0..NZ-1 =*2675 *-----------------------------------------------------------------------------*2676 c 811 !-----------------------------------------------------------------------------* 812 != PURPOSE: =* 813 != Calculate slant path over vertical depth ds/dh in spherical geometry. =* 814 != Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =* 815 != for computing the radiation field available for photolysis and heating =* 816 != at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =* 817 !-----------------------------------------------------------------------------* 818 != PARAMETERS: =* 819 != NZ - INTEGER, number of specified altitude levels in the working (I)=* 820 != grid =* 821 != Z - REAL, specified altitude working grid (km) (I)=* 822 != ZEN - REAL, solar zenith angle (degrees) (I)=* 823 != DSDH - REAL, slant path of direct beam through each layer crossed (O)=* 824 != when travelling from the top of the atmosphere to layer i; =* 825 != DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* 826 != NID - INTEGER, number of layers crossed by the direct beam when (O)=* 827 != travelling from the top of the atmosphere to layer i; =* 828 != NID(i), i = 0..NZ-1 =* 829 !-----------------------------------------------------------------------------* 830 2677 831 implicit none 2678 832 2679 * input 2680 INTEGER nlev 2681 REAL zen, z(nlev) 2682 2683 * output 833 ! input 834 835 integer, intent(in) :: nlev 836 real, dimension(nlev), intent(in) :: z 837 real, intent(in) :: zen 838 839 ! output 840 2684 841 INTEGER nid(0:nlev) 2685 842 REAL dsdh(0:nlev,nlev) 2686 843 2687 * more program constants 844 ! more program constants 845 2688 846 REAL re, ze(nlev) 2689 847 REAL dr … … 2691 849 parameter (radius = 3393.) 2692 850 2693 * local 2694 2695 REAL pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm 2696 INTEGER i, j, k 2697 INTEGER id 2698 2699 INTEGER nlay 851 ! local 852 853 real :: pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm 854 integer :: i, j, k, id, nlay 855 2700 856 REAL zd(0:nlev-1) 2701 857 2702 *-----------------------------------------------------------------------------858 !----------------------------------------------------------------------------- 2703 859 2704 860 pi = acos(-1.0) … … 2706 862 zenrad = zen*dr 2707 863 2708 * number of layers: 864 ! number of layers: 865 2709 866 nlay = nlev - 1 2710 867 2711 * include the elevation above sea level to the radius of Mars: 868 ! include the elevation above sea level to the radius of Mars: 869 2712 870 re = radius + z(1) 2713 * correspondingly z changed to the elevation above Mars surface: 871 872 ! correspondingly z changed to the elevation above Mars surface: 873 2714 874 DO k = 1, nlev 2715 875 ze(k) = z(k) - z(1) 2716 876 END DO 2717 877 2718 * inverse coordinate of z 878 ! inverse coordinate of z 879 2719 880 zd(0) = ze(nlev) 2720 881 DO k = 1, nlay … … 2722 883 END DO 2723 884 2724 * initialize dsdh(i,j), nid(i) 2725 DO i = 0, nlev 2726 nid(i) = 0 2727 DO j = 1, nlev 2728 dsdh(i,j) = 0. 2729 END DO 2730 END DO 2731 2732 * calculate ds/dh of every layer 2733 DO 100 i = 0, nlay 2734 2735 rpsinz = (re + zd(i)) * SIN(zenrad) 885 ! initialise dsdh(i,j), nid(i) 886 887 nid(:) = 0. 888 dsdh(:,:) = 0. 889 890 ! calculate ds/dh of every layer 891 892 do i = 0,nlay 893 rpsinz = (re + zd(i))*sin(zenrad) 2736 894 2737 895 IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN … … 2739 897 ELSE 2740 898 2741 * 2742 * Find index of layer in which the screening height lies 2743 * 899 ! Find index of layer in which the screening height lies 900 2744 901 id = i 2745 IF( zen .GT. 90.0 ) THEN2746 DO 10 j = 1,nlay902 if (zen > 90.) then 903 do j = 1,nlay 2747 904 IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND. 2748 905 $ (rpsinz .GE. ( zd(j) + re )) ) id = j 2749 10 CONTINUE2750 END IF906 end do 907 end if 2751 908 2752 DO 20 j = 1, id 2753 2754 sm = 1.0 2755 IF(j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0) 2756 $ sm = -1.0 909 do j = 1,id 910 sm = 1.0 911 IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0) 912 $ sm = -1.0 2757 913 2758 rj = re + zd(j-1)2759 rjp1 = re + zd(j)914 rj = re + zd(j-1) 915 rjp1 = re + zd(j) 2760 916 2761 dhj = zd(j-1) - zd(j)917 dhj = zd(j-1) - zd(j) 2762 918 2763 ga = rj*rj - rpsinz*rpsinz 2764 gb = rjp1*rjp1 - rpsinz*rpsinz 2765 IF (ga .LT. 0.0) ga = 0.0 2766 IF (gb .LT. 0.0) gb = 0.0 919 ga = rj*rj - rpsinz*rpsinz 920 gb = rjp1*rjp1 - rpsinz*rpsinz 921 922 IF (ga .LT. 0.0) ga = 0.0 923 IF (gb .LT. 0.0) gb = 0.0 2767 924 2768 IF(id.GT.i .AND. j.EQ.id) THEN 2769 dsj = SQRT( ga ) 2770 ELSE 2771 dsj = SQRT( ga ) - sm*SQRT( gb ) 2772 END IF 2773 dsdh(i,j) = dsj / dhj 2774 2775 20 CONTINUE 2776 2777 nid(i) = id 2778 2779 END IF 2780 2781 100 CONTINUE 2782 c 2783 RETURN 2784 END 925 IF (id.GT.i .AND. j.EQ.id) THEN 926 dsj = sqrt(ga) 927 ELSE 928 dsj = sqrt(ga) - sm*sqrt(gb) 929 END IF 930 dsdh(i,j) = dsj/dhj 931 end do 932 nid(i) = id 933 end if 934 end do ! i = 0,nlay 935 936 end subroutine sphers 2785 937 2786 938 !============================================================================== … … 2792 944 implicit none 2793 945 2794 * input 2795 2796 integer :: nw ! number of wavelength grid points 2797 INTEGER nlev, iw 946 ! input 947 948 integer, intent(in) :: nlev, nw, iw ! number of wavelength grid points 2798 949 REAL ag 2799 950 REAL zen … … 2806 957 REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw) 2807 958 2808 *output959 ! output 2809 960 2810 961 REAL edir(nlev), edn(nlev), eup(nlev) 2811 962 REAL fdir(nlev), fdn(nlev), fup(nlev) 2812 963 2813 *local:964 ! local: 2814 965 2815 966 REAL dt(nlev), om(nlev), g(nlev) … … 2818 969 real, parameter :: largest = 1.e+36 2819 970 2820 *specific two ps2str971 ! specific two ps2str 2821 972 2822 973 REAL ediri(nlev), edni(nlev), eupi(nlev) … … 2825 976 logical, save :: delta = .true. 2826 977 2827 *_______________________________________________________________________2828 2829 *initialize:978 !_______________________________________________________________________ 979 980 ! initialize: 2830 981 2831 982 do i = 1, nlev … … 2837 988 edn(i) = 0. 2838 989 end do 2839 c 990 2840 991 do i = 1, nlev - 1 2841 992 dscld = dtcld(i,iw)*omcld(i,iw) 2842 993 dacld = dtcld(i,iw)*(1.-omcld(i,iw)) 2843 c 994 2844 995 dsaer = dtaer(i,iw)*omaer(i,iw) 2845 996 daaer = dtaer(i,iw)*(1.-omaer(i,iw)) 2846 c 997 2847 998 dtsct = dtrl(i,iw) + dscld + dsaer 2848 999 dtabs = dagas(i,iw) + dacld + daaer 2849 c 1000 2850 1001 dtabs = amax1(dtabs,1./largest) 2851 1002 dtsct = amax1(dtsct,1./largest) 2852 c 2853 cinvert z-coordinate:2854 c 1003 1004 ! invert z-coordinate: 1005 2855 1006 ii = nlev - i 2856 1007 dt(ii) = dtsct + dtabs … … 2860 1011 $ gaer(i,iw)*dsaer)/dtsct 2861 1012 end do 2862 c 2863 ccall rt routine:2864 c 1013 1014 ! call rt routine: 1015 2865 1016 call ps2str(nlev, zen, ag, dt, om, g, 2866 1017 $ dsdh, nid, delta, 2867 1018 $ fdiri, fupi, fdni, ediri, eupi, edni) 2868 c 2869 coutput (invert z-coordinate)2870 c 1019 1020 ! output (invert z-coordinate) 1021 2871 1022 do i = 1, nlev 2872 1023 ii = nlev - i + 1 … … 2878 1029 edn(i) = edni(ii) 2879 1030 end do 2880 c 2881 return 2882 end 1031 1032 end subroutine rtlink 2883 1033 2884 1034 *=============================================================================* 2885 1035 2886 SUBROUTINEps2str(nlev,zen,rsfc,tauu,omu,gu,1036 subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu, 2887 1037 $ dsdh, nid, delta, 2888 1038 $ fdr, fup, fdn, edr, eup, edn) 2889 1039 2890 *-----------------------------------------------------------------------------* 2891 *= PURPOSE: =* 2892 *= Solve two-stream equations for multiple layers. The subroutine is based =* 2893 *= on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=* 2894 *= It contains 9 two-stream methods to choose from. A pseudo-spherical =* 2895 *= correction has also been added. =* 2896 *-----------------------------------------------------------------------------* 2897 *= PARAMETERS: =* 2898 *= NLEVEL - INTEGER, number of specified altitude levels in the working (I)=* 2899 *= grid =* 2900 *= ZEN - REAL, solar zenith angle (degrees) (I)=* 2901 *= RSFC - REAL, surface albedo at current wavelength (I)=* 2902 *= TAUU - REAL, unscaled optical depth of each layer (I)=* 2903 *= OMU - REAL, unscaled single scattering albedo of each layer (I)=* 2904 *= GU - REAL, unscaled asymmetry parameter of each layer (I)=* 2905 *= DSDH - REAL, slant path of direct beam through each layer crossed (I)=* 2906 *= when travelling from the top of the atmosphere to layer i; =* 2907 *= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* 2908 *= NID - INTEGER, number of layers crossed by the direct beam when (I)=* 2909 *= travelling from the top of the atmosphere to layer i; =* 2910 *= NID(i), i = 0..NZ-1 =* 2911 *= DELTA - LOGICAL, switch to use delta-scaling (I)=* 2912 *= .TRUE. -> apply delta-scaling =* 2913 *= .FALSE.-> do not apply delta-scaling =* 2914 *= FDR - REAL, contribution of the direct component to the total (O)=* 2915 *= actinic flux at each altitude level =* 2916 *= FUP - REAL, contribution of the diffuse upwelling component to (O)=* 2917 *= the total actinic flux at each altitude level =* 2918 *= FDN - REAL, contribution of the diffuse downwelling component to (O)=* 2919 *= the total actinic flux at each altitude level =* 2920 *= EDR - REAL, contribution of the direct component to the total (O)=* 2921 *= spectral irradiance at each altitude level =* 2922 *= EUP - REAL, contribution of the diffuse upwelling component to (O)=* 1040 !-----------------------------------------------------------------------------* 1041 != PURPOSE: =* 1042 != Solve two-stream equations for multiple layers. The subroutine is based =* 1043 != on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=* 1044 != It contains 9 two-stream methods to choose from. A pseudo-spherical =* 1045 != correction has also been added. =* 1046 !-----------------------------------------------------------------------------* 1047 != PARAMETERS: =* 1048 != NLEVEL - INTEGER, number of specified altitude levels in the working (I)=* 1049 != grid =* 1050 != ZEN - REAL, solar zenith angle (degrees) (I)=* 1051 != RSFC - REAL, surface albedo at current wavelength (I)=* 1052 != TAUU - REAL, unscaled optical depth of each layer (I)=* 1053 != OMU - REAL, unscaled single scattering albedo of each layer (I)=* 1054 != GU - REAL, unscaled asymmetry parameter of each layer (I)=* 1055 != DSDH - REAL, slant path of direct beam through each layer crossed (I)=* 1056 != when travelling from the top of the atmosphere to layer i; =* 1057 != DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* 1058 != NID - INTEGER, number of layers crossed by the direct beam when (I)=* 1059 != travelling from the top of the atmosphere to layer i; =* 1060 != NID(i), i = 0..NZ-1 =* 1061 != DELTA - LOGICAL, switch to use delta-scaling (I)=* 1062 != .TRUE. -> apply delta-scaling =* 1063 != .FALSE.-> do not apply delta-scaling =* 1064 != FDR - REAL, contribution of the direct component to the total (O)=* 1065 != actinic flux at each altitude level =* 1066 != FUP - REAL, contribution of the diffuse upwelling component to (O)=* 1067 != the total actinic flux at each altitude level =* 1068 != FDN - REAL, contribution of the diffuse downwelling component to (O)=* 1069 != the total actinic flux at each altitude level =* 1070 != EDR - REAL, contribution of the direct component to the total (O)=* 1071 != spectral irradiance at each altitude level =* 1072 != EUP - REAL, contribution of the diffuse upwelling component to (O)=* 1073 != the total spectral irradiance at each altitude level =* 1074 != EDN - REAL, contribution of the diffuse downwelling component to (O)=* 2923 1075 *= the total spectral irradiance at each altitude level =* 2924 *= EDN - REAL, contribution of the diffuse downwelling component to (O)=* 2925 *= the total spectral irradiance at each altitude level =* 2926 *-----------------------------------------------------------------------------* 2927 c 1076 !-----------------------------------------------------------------------------* 1077 2928 1078 implicit none 2929 c 2930 ******* 2931 * input: 2932 ******* 1079 1080 ! input: 1081 2933 1082 INTEGER nlev 2934 1083 REAL zen, rsfc … … 2938 1087 LOGICAL delta 2939 1088 2940 ******* 2941 * output: 2942 ******* 1089 ! output: 1090 2943 1091 REAL fup(nlev),fdn(nlev),fdr(nlev) 2944 1092 REAL eup(nlev),edn(nlev),edr(nlev) 2945 1093 2946 ******* 2947 * local: 2948 ******* 1094 ! local: 1095 2949 1096 REAL tausla(0:nlev), tauc(0:nlev) 2950 1097 REAL mu2(0:nlev), mu, sum 2951 1098 2952 * internal coefficients and matrix 1099 ! internal coefficients and matrix 1100 2953 1101 REAL lam(nlev),taun(nlev),bgam(nlev) 2954 1102 REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev) … … 2958 1106 REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev) 2959 1107 2960 ******* 2961 * other: 2962 ******* 1108 ! other: 1109 2963 1110 REAL pifs, fdn0 2964 1111 REAL gi(nlev), omi(nlev), tempg … … 2968 1115 real, parameter :: precis = 1.e-7 2969 1116 2970 *For calculations of Associated Legendre Polynomials for GAMA1,2,3,42971 *in delta-function, modified quadrature, hemispheric constant,2972 *Hybrid modified Eddington-delta function metods, p633,Table1.2973 *W.E.Meador and W.R.Weaver, GAS,1980,v37,p.6302974 *W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,2975 *uncomment the following two lines and the appropriate statements further2976 *down.2977 CREAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,2978 C> BETA1, BETAn, amu1, subd1117 ! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4 1118 ! in delta-function, modified quadrature, hemispheric constant, 1119 ! Hybrid modified Eddington-delta function metods, p633,Table1. 1120 ! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 1121 ! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, 1122 ! uncomment the following two lines and the appropriate statements further 1123 ! down. 1124 ! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0, 1125 ! > BETA1, BETAn, amu1, subd 2979 1126 2980 1127 REAL expon, expon0, expon1, divisr, temp, up, dn … … 2984 1131 INTEGER i, j 2985 1132 2986 *Some additional program constants:1133 ! Some additional program constants: 2987 1134 2988 1135 real pi, dr 2989 1136 REAL eps 2990 1137 PARAMETER (eps = 1.E-3) 2991 *_______________________________________________________________________2992 2993 *MU = cosine of solar zenith angle2994 *RSFC = surface albedo2995 *TAUU = unscaled optical depth of each layer2996 *OMU = unscaled single scattering albedo2997 *GU = unscaled asymmetry factor2998 *KLEV = max dimension of number of layers in atmosphere2999 *NLAYER = number of layers in the atmosphere3000 *NLEVEL = nlayer + 1 = number of levels3001 3002 *initial conditions: pi*solar flux = 1; diffuse incidence = 01138 !_______________________________________________________________________ 1139 1140 ! MU = cosine of solar zenith angle 1141 ! RSFC = surface albedo 1142 ! TAUU = unscaled optical depth of each layer 1143 ! OMU = unscaled single scattering albedo 1144 ! GU = unscaled asymmetry factor 1145 ! KLEV = max dimension of number of layers in atmosphere 1146 ! NLAYER = number of layers in the atmosphere 1147 ! NLEVEL = nlayer + 1 = number of levels 1148 1149 ! initial conditions: pi*solar flux = 1; diffuse incidence = 0 3003 1150 3004 1151 pifs = 1. … … 3011 1158 mu = COS(zen*dr) 3012 1159 3013 ************** compute coefficients for each layer:3014 *GAM1 - GAM4 = 2-stream coefficients, different for different approximations3015 *EXPON0 = calculation of e when TAU is zero3016 *EXPON1 = calculation of e when TAU is TAUN3017 *CUP and CDN = calculation when TAU is zero3018 *CUPTN and CDNTN = calc. when TAU is TAUN3019 *DIVISR = prevents division by zero1160 !************* compute coefficients for each layer: 1161 ! GAM1 - GAM4 = 2-stream coefficients, different for different approximations 1162 ! EXPON0 = calculation of e when TAU is zero 1163 ! EXPON1 = calculation of e when TAU is TAUN 1164 ! CUP and CDN = calculation when TAU is zero 1165 ! CUPTN and CDNTN = calc. when TAU is TAUN 1166 ! DIVISR = prevents division by zero 3020 1167 3021 1168 do j = 0, nlev … … 3024 1171 mu2(j) = 1./SQRT(largest) 3025 1172 end do 3026 c 1173 3027 1174 IF (.NOT. delta) THEN 3028 1175 DO i = 1, nlayer … … 3033 1180 ELSE 3034 1181 3035 *delta-scaling. Have to be done for delta-Eddington approximation,3036 *delta discrete ordinate, Practical Improved Flux Method, delta function,3037 *and Hybrid modified Eddington-delta function methods approximations1182 ! delta-scaling. Have to be done for delta-Eddington approximation, 1183 ! delta discrete ordinate, Practical Improved Flux Method, delta function, 1184 ! and Hybrid modified Eddington-delta function methods approximations 3038 1185 3039 1186 DO i = 1, nlayer … … 3044 1191 END DO 3045 1192 END IF 3046 * 3047 * calculate slant optical depth at the top of the atmosphere when zen>90. 3048 * in this case, higher altitude of the top layer is recommended which can 3049 * be easily changed in gridz.f. 3050 * 1193 1194 ! calculate slant optical depth at the top of the atmosphere when zen>90. 1195 ! in this case, higher altitude of the top layer is recommended. 1196 3051 1197 IF (zen .GT. 90.0) THEN 3052 1198 IF (nid(0) .LT. 0) THEN … … 3060 1206 END IF 3061 1207 END IF 3062 * 1208 3063 1209 DO 11, i = 1, nlayer 3064 1210 g = gi(i) … … 3066 1212 tauc(i) = tauc(i-1) + taun(i) 3067 1213 3068 *stay away from 1 by precision. For g, also stay away from -11214 ! stay away from 1 by precision. For g, also stay away from -1 3069 1215 3070 1216 tempg = AMIN1(abs(g),1. - precis) … … 3072 1218 om = AMIN1(om,1.-precis) 3073 1219 3074 *calculate slant optical depth3075 *1220 ! calculate slant optical depth 1221 3076 1222 IF (nid(i) .LT. 0) THEN 3077 1223 tausla(i) = largest … … 3093 1239 END IF 3094 1240 END IF 3095 * 3096 *** the following gamma equations are from pg 16,289, Table 13097 *** save mu1 for each approx. for use in converting irradiance to actinic flux3098 3099 *Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):1241 1242 !** the following gamma equations are from pg 16,289, Table 1 1243 !** save mu1 for each approx. for use in converting irradiance to actinic flux 1244 1245 ! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452): 3100 1246 3101 1247 c gam1 = (7. - om*(4. + 3.*g))/4. … … 3335 1481 j = j + 1 3336 1482 60 CONTINUE 3337 c 3338 RETURN 3339 END 1483 1484 end subroutine ps2str 3340 1485 3341 1486 *=============================================================================* 3342 1487 3343 SUBROUTINEtridiag(a,b,c,r,u,n)3344 3345 *_______________________________________________________________________3346 *solves tridiagonal system. From Numerical Recipies, p. 403347 *_______________________________________________________________________1488 subroutine tridiag(a,b,c,r,u,n) 1489 1490 !_______________________________________________________________________ 1491 ! solves tridiagonal system. From Numerical Recipies, p. 40 1492 !_______________________________________________________________________ 3348 1493 3349 1494 IMPLICIT NONE 3350 1495 3351 * input: 1496 ! input: 1497 3352 1498 INTEGER n 3353 1499 REAL a, b, c, r 3354 1500 DIMENSION a(n),b(n),c(n),r(n) 3355 1501 3356 * output: 1502 ! output: 1503 3357 1504 REAL u 3358 1505 DIMENSION u(n) 3359 1506 3360 * local: 1507 ! local: 1508 3361 1509 INTEGER j 3362 1510 3363 1511 REAL bet, gam 3364 1512 DIMENSION gam(n) 3365 *_______________________________________________________________________1513 !_______________________________________________________________________ 3366 1514 3367 1515 IF (b(1) .EQ. 0.) STOP 1001 … … 3377 1525 u(j) = u(j) - gam(j + 1)*u(j + 1) 3378 1526 12 CONTINUE 3379 *_______________________________________________________________________ 3380 3381 RETURN 3382 END 1527 !_______________________________________________________________________ 1528 1529 end subroutine tridiag 1530 1531 end subroutine photolysis_online
Note: See TracChangeset
for help on using the changeset viewer.
