program aeropt ! ! Read and map the tabulated aerosol optical spectral data ! onto corresponding sw/lw radiation spectral bands. ! integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_phys = 8 integer, parameter :: KAERBND=61 ! num of bands for aer data (gocart) integer, parameter :: KRHLEV =36 ! num of rh levels for rh-dep components integer, parameter :: KCM1 = 8 ! num of rh independent aer species integer, parameter :: KCM2 = 5 ! num of rh dependent aer species integer, parameter :: NBDSW = 7 ! total num of sw bands integer, parameter :: NWVTOT = 57600 ! total num of wvnum included integer, parameter :: NWVSOL = 151 ! num of wvnum regions where solar ! flux is constant integer, parameter :: NIAERCM = 102 integer, parameter :: NOAER = 20, NODMP = 21 INTEGER, PARAMETER :: NP=100, NP2=200, nWave=100, nAero=6, n_p=36 real (kind=kind_phys), parameter :: f_zero = 0.0 real (kind=kind_phys), parameter :: f_one = 1.0 real :: wvnum1(NBDSW), wvnum2(NBDSW) real (kind=kind_phys) :: sumsol real (kind=kind_phys) :: sumk, sumok, sumokg, sumreft, & & sp, refb, reft, rsolbd integer, dimension(NBDSW) :: nv1, nv2 INTEGER :: ib, ii, iw1, iw2, nc, ni real (kind=kind_phys), dimension(NBDSW,KAERBND) :: solwaer real (kind=kind_phys), dimension(NBDSW) :: solbnd real (kind=kind_phys), dimension(NWVTOT) :: solfwv ! --- solar flux (w/m**2) in each wvnumb region where it is constant real (kind=kind_phys), dimension(NWVSOL) :: s0intv ! --- number of wavenumbers in each region where the solar flux is constant integer, dimension(NWVSOL) :: nwvns0 ! - spectral band structure: integer :: iendwv (KAERBND) ! ending wavenumber (cm**-1) for each band ! --- tabulated aerosol optical properties (input dataset) ! - relative humidity independent aerosol optical properties: ! ===> species : dust (8 bins) real (kind=kind_phys) :: rhidext0(KAERBND,KCM1) ! extinction coefficient real (kind=kind_phys) :: rhidssa0(KAERBND,KCM1) ! single scattering albedo real (kind=kind_phys) :: rhidasy0(KAERBND,KCM1) ! asymmetry parameter ! - relative humidity dependent aerosol optical properties: ! ===> species : soot, suso, waso, ssam, sscm real (kind=kind_phys) :: rhdpext0(KAERBND,KRHLEV,KCM2) ! extinction coefficient real (kind=kind_phys) :: rhdpssa0(KAERBND,KRHLEV,KCM2) ! single scattering albedo real (kind=kind_phys) :: rhdpasy0(KAERBND,KRHLEV,KCM2) ! asymmetry parameter ! --- aerosol optical properties mapped onto specified spectral bands ! - relative humidity dependent aerosol optical properties: real (kind=kind_phys) :: extrhd(KRHLEV,KCM2,NBDSW) ! extinction coefficient for sw band real (kind=kind_phys) :: ssarhd(KRHLEV,KCM2,NBDSW) ! single scattering albedo for sw band real (kind=kind_phys) :: asyrhd(KRHLEV,KCM2,NBDSW) ! asymmetry parameter for sw band ! - relative humidity independent aerosol optical properties: real (kind=kind_phys) :: extrhi(KCM1,NBDSW) ! extinction coefficient for sw band real (kind=kind_phys) :: ssarhi(KCM1,NBDSW) ! single scattering albedo for sw band real (kind=kind_phys) :: asyrhi(KCM1,NBDSW) ! asymmetry parameter for sw band INTEGER :: NW, NS, nH, n_bin real (kind=kind_io8), Dimension(NP2) :: Angle, Cos_Angle, Cos_Weight real (kind=kind_io8), Dimension(n_p,nAero) :: RH, rm, reff real (kind=kind_io8), Dimension(nWave,n_p,nAero):: ext0, sca0, asy0 real (kind=kind_io8), Dimension(NP2,n_p,nWave,nAero) :: ph0 real (kind=kind_io8) :: wavelength(nWave), density(nAero), & & sigma(nAero), wave, n_fac, PI, t1, s1, g1 CHARACTER(len=80) :: AerosolName(nAero) INTEGER :: i, j, k, l, ij character :: aerosol_file*30 character :: aerluts_file*30 logical :: file_exist integer :: indx_dust(8) ! map 36 dust bins to gocart size bins real (kind=kind_phys), dimension(KRHLEV) :: rhlev data rhlev (:)/ .00, .05, .10, .15, .20, .25, .30, .35, & & .40, .45, .50, .55, .60, .65, .70, .75, & & .80, .81, .82, .83, .84, .85, .86, .87, & & .88, .89, .90, .91, .92, .93, .94, .95, & & .96, .97, .98, .99/ ! ! From MODIS Aerosol Collection 5 (MOD04/MYD04) ATBD ! Band Number Bandwidth CentralWavelength Resolution ! 1 0.620-0.670 0.646 250 ! 2 0.841-0.876 0.855 250 ! 3 0.459-0.479 0.466 500 ! 4 0.545-0.565 0.553 500 ! 5 1.230-1.250 1.243 500 ! 6 1.628-1.652 1.632 500 ! 7 2.105-2.155 2.119 500 ! From AERONET Spectral Corrections/Components document ! 340 (2nm), 380 (4nm), 440 (10nm), 500 (10nm), 675 (10mm), 870 (10nm), ! 940 (10nm), 1020 (10nm), 1640 (25nm) ! ! AVHRR channels: 0.63, 0.86, 1.61 ! UX index spectrum: 0.29 - 0.40 ! IR channels (for STAR/IOSSPDT group): 11.1 micron ! AERONET Angstrom: 440-870 nm ! Therefore, compute the following channels: ! AERONET: 0.34, 0.44 ! MODIS: 0.66, 0.86, 0.55, 1.63 ! data /0.34, 0.44, 0.55, 0.66, 0.86, 1.63, 11.1/ data wvnum1/0.338, 0.430, 0.545, 0.62, 0.841, 1.628, 11.0/ data wvnum2/0.342, 0.450, 0.565, 0.67, 0.876, 1.652, 11.2/ data nwvns0 / 100, 11, 14, 18, 24, 33, 50, 83, 12, 12, & & 13, 15, 15, 17, 18, 20, 21, 24, 26, 30, 32, 37, 42, & & 47, 55, 64, 76, 91, 111, 139, 179, 238, 333, 41, 42, 45, & & 46, 48, 51, 53, 55, 58, 61, 64, 68, 71, 75, 79, 84, & & 89, 95, 101, 107, 115, 123, 133, 142, 154, 167, 181, 197, 217, & & 238, 263, 293, 326, 368, 417, 476, 549, 641, 758, 909, 101, 103, & & 105, 108, 109, 112, 115, 117, 119, 122, 125, 128, 130, 134, 137, & & 140, 143, 147, 151, 154, 158, 163, 166, 171, 175, 181, 185, 190, & & 196, 201, 207, 213, 219, 227, 233, 240, 248, 256, 264, 274, 282, & & 292, 303, 313, 325, 337, 349, 363, 377, 392, 408, 425, 444, 462, & & 483, 505, 529, 554, 580, 610, 641, 675, 711, 751, 793, 841, 891, & & 947,1008,1075,1150,1231,1323,1425,1538,1667,1633,14300 / data s0intv( 1: 50) / & & 1.60000E-6, 2.88000E-5, 3.60000E-5, 4.59200E-5, 6.13200E-5, & & 8.55000E-5, 1.28600E-4, 2.16000E-4, 2.90580E-4, 3.10184E-4, & & 3.34152E-4, 3.58722E-4, 3.88050E-4, 4.20000E-4, 4.57056E-4, & & 4.96892E-4, 5.45160E-4, 6.00600E-4, 6.53600E-4, 7.25040E-4, & & 7.98660E-4, 9.11200E-4, 1.03680E-3, 1.18440E-3, 1.36682E-3, & & 1.57560E-3, 1.87440E-3, 2.25500E-3, 2.74500E-3, 3.39840E-3, & & 4.34000E-3, 5.75400E-3, 7.74000E-3, 9.53050E-3, 9.90192E-3, & & 1.02874E-2, 1.06803E-2, 1.11366E-2, 1.15830E-2, 1.21088E-2, & & 1.26420E-2, 1.32250E-2, 1.38088E-2, 1.44612E-2, 1.51164E-2, & & 1.58878E-2, 1.66500E-2, 1.75140E-2, 1.84450E-2, 1.94106E-2 / data s0intv( 51:100) / & & 2.04864E-2, 2.17248E-2, 2.30640E-2, 2.44470E-2, 2.59840E-2, & & 2.75940E-2, 2.94138E-2, 3.13950E-2, 3.34800E-2, 3.57696E-2, & & 3.84054E-2, 4.13490E-2, 4.46880E-2, 4.82220E-2, 5.22918E-2, & & 5.70078E-2, 6.19888E-2, 6.54720E-2, 6.69060E-2, 6.81226E-2, & & 6.97788E-2, 7.12668E-2, 7.27100E-2, 7.31610E-2, 7.33471E-2, & & 7.34814E-2, 7.34717E-2, 7.35072E-2, 7.34939E-2, 7.35202E-2, & & 7.33249E-2, 7.31713E-2, 7.35462E-2, 7.36920E-2, 7.23677E-2, & & 7.25023E-2, 7.24258E-2, 7.20766E-2, 7.18284E-2, 7.32757E-2, & & 7.31645E-2, 7.33277E-2, 7.36128E-2, 7.33752E-2, 7.28965E-2, & & 7.24924E-2, 7.23307E-2, 7.21050E-2, 7.12620E-2, 7.10903E-2 / data s0intv(101:151) / 7.12714E-2, & & 7.08012E-2, 7.03752E-2, 7.00350E-2, 6.98639E-2, 6.90690E-2, & & 6.87621E-2, 6.52080E-2, 6.65184E-2, 6.60038E-2, 6.47615E-2, & & 6.44831E-2, 6.37206E-2, 6.24102E-2, 6.18698E-2, 6.06320E-2, & & 5.83498E-2, 5.67028E-2, 5.51232E-2, 5.48645E-2, 5.12340E-2, & & 4.85581E-2, 4.85010E-2, 4.79220E-2, 4.44058E-2, 4.48718E-2, & & 4.29373E-2, 4.15242E-2, 3.81744E-2, 3.16342E-2, 2.99615E-2, & & 2.92740E-2, 2.67484E-2, 1.76904E-2, 1.40049E-2, 1.46224E-2, & & 1.39993E-2, 1.19574E-2, 1.06386E-2, 1.00980E-2, 8.63808E-3, & & 6.52736E-3, 4.99410E-3, 4.39350E-3, 2.21676E-3, 1.33812E-3, & & 1.12320E-3, 5.59000E-4, 3.60000E-4, 2.98080E-4, 7.46294E-5 / data aerosol_file /"NCEP_AEROSOL.bin"/ data aerluts_file /"AEROSOL_LUTS.dat"/ data AerosolName/ ' Dust ', ' Soot ', ' SUSO ', ' WASO ', & & ' SSAM ', ' SSCM '/ !! 8 dust bins !! 1 2 3 4 5 6 7 8 !! .1-.18, .18-.3, .3-.6, 0.6-1.0, 1.0-1.8, 1.8-3, 3-6, 6-10 <-- def !! 0.1399 0.2399 0.4499 0.8000 1.3994 2.3964 4.4964 7.9887 <-- reff data indx_dust/4, 8, 12, 18, 21, 24, 30, 36/ PI = acos(-1.d0) ! ! read input gocart aerosol optical data from CRTM Mie code calculations ! inquire (file = aerosol_file, exist = file_exist) if ( file_exist ) then print *,' Open aerosol data file ',aerosol_file close (NIAERCM) open(unit=NIAERCM,file=aerosol_file,status='OLD',form='UNFORMATTED') else print *,' Requested aerosol data file not found,"', aerosol_file print *,' Exit now' stop 1003 endif ! end if_file_exist_block ! READ(NIAERCM) (Cos_Angle(i),i=1,NP) READ(NIAERCM) (Cos_Weight(i),i=1,NP) READ(NIAERCM) READ(NIAERCM) READ(NIAERCM) NW,NS READ(NIAERCM) READ(NIAERCM) (wavelength(i),i=1,NW) write(400,*) 'NW, NS, NP = ', NW, NS, NP ! --- check nAero and NW if (NW /= KAERBND) then print *, " Incorrect spectral band, abort ", NW stop 1004 endif ! --- convert wavelength to wavenumber DO i = 1, KAERBND iendwv(i) = int(10000. / wavelength(i)) print *,' i,wn,lamda:', i,iendwv(i),wavelength(i) write(400,*) 'i,wn,lamda:', i,iendwv(i),wavelength(i) ENDDO DO j = 1, nAero print *,' Read LUTs:', j,AerosolName(j) READ(NIAERCM) READ(NIAERCM) READ(NIAERCM) n_bin, density(j), sigma(j) READ(NIAERCM) READ(NIAERCM) (RH(i,j),i=1, n_bin) READ(NIAERCM) READ(NIAERCM) (rm(i,j),i=1, n_bin) READ(NIAERCM) READ(NIAERCM) (reff(i,j),i=1, n_bin) ! --- check n_bin if (n_bin /= KRHLEV ) then print *, "Incorrect rh levels, abort ", n_bin stop 1005 endif ! --- read luts if( j == 3 ) then write(400,*) 'n_bin:', n_bin write(400,*) 'nw:', nw endif DO k = 1, NW READ(NIAERCM) wave,(ext0(k,L,j),L=1,n_bin) READ(NIAERCM) (sca0(k,L,j),L=1,n_bin) READ(NIAERCM) (asy0(k,L,j),L=1,n_bin) READ(NIAERCM) (ph0(1:NP2,L,k,j),L=1,n_bin) ! if ( j == 3 .and. k == 54) then if ( j == 3 ) then write(400,*) k, (ext0(k,L,j),L=1,n_bin) endif END DO ! --- map luts input to module variables IF (AerosolName(j) == ' Dust ' ) then do i = 1, KAERBND do k = 1, KCM1 ij = indx_dust(k) rhidext0(i,k)=ext0(i,ij,j) rhidssa0(i,k)=sca0(i,ij,j) rhidasy0(i,k)=asy0(i,ij,j) enddo enddo ELSE if (AerosolName(j) == ' Soot ') ij = 1 if (AerosolName(j) == ' SUSO ') ij = 2 if (AerosolName(j) == ' WASO ') ij = 3 if (AerosolName(j) == ' SSAM ') ij = 4 if (AerosolName(j) == ' SSCM ') ij = 5 do i = 1, KAERBND do k = 1, KRHLEV rhdpext0(i,k,ij)=ext0(i,k,j) rhdpssa0(i,k,ij)=sca0(i,k,j) rhdpasy0(i,k,ij)=asy0(i,k,j) enddo enddo ! write(400,*) k, (ext0(k,L,j),L=1,n_bin) ENDIF END DO ! ! --- ... define the one wavenumber solar fluxes based on toa solar ! spectral distribution do ib = 1, NWVSOL if ( ib == 1 ) then iw1 = 1 else iw1 = iw1 + nwvns0(ib-1) endif iw2 = iw1 + nwvns0(ib) - 1 do ii = iw1, iw2 solfwv(ii) = s0intv(ib) enddo enddo ! --- ... compute solar flux weights and interval indices for mapping ! spectral bands between sw radiation and aerosol data ! solbnd (:) = f_zero solwaer(:,:) = f_zero do ib = 1, NBDSW ii = 1 iw1 = nint(10000./wvnum2(ib)) iw2 = nint(10000./wvnum1(ib)) Lab_swdowhile : do while ( iw1 > iendwv(ii) ) if ( ii == KAERBND ) exit Lab_swdowhile ii = ii + 1 enddo Lab_swdowhile sumsol = f_zero nv1(ib) = ii do iw = iw1, iw2 solbnd(ib) = solbnd(ib) + solfwv(iw) sumsol = sumsol + solfwv(iw) if ( iw == iendwv(ii) ) then solwaer(ib,ii) = sumsol if ( ii < KAERBND ) then sumsol = f_zero ii = ii + 1 endif endif enddo if ( iw2 /= iendwv(ii) ) then solwaer(ib,ii) = sumsol endif nv2(ib) = ii enddo ! --- ... loop for each sw radiation spectral band do ib = 1, NBDSW rsolbd = f_one / solbnd(ib) ! --- for rh independent aerosol species do nc = 1, KCM1 sumk = f_zero sumok = f_zero sumokg = f_zero sumreft = f_zero do ni = nv1(ib), nv2(ib) sp = sqrt( (f_one - rhidssa0(ni,nc)) & & / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) ) reft = (f_one - sp) / (f_one + sp) sumreft = sumreft + reft*solwaer(ib,ni) sumk = sumk + rhidext0(ni,nc)*solwaer(ib,ni) sumok = sumok + rhidssa0(ni,nc)*solwaer(ib,ni) & & * rhidext0(ni,nc) sumokg = sumokg + rhidssa0(ni,nc)*solwaer(ib,ni) & & * rhidext0(ni,nc)*rhidasy0(ni,nc) enddo refb = sumreft * rsolbd extrhi(nc,ib) = sumk * rsolbd asyrhi(nc,ib) = sumokg / (sumok + 1.0e-10) ssarhi(nc,ib) = 4.0*refb & & / ( (f_one+refb)**2 - asyrhi(nc,ib)*(f_one-refb)**2 ) enddo ! end do_nc_block for rh-ind aerosols ! --- for rh dependent aerosols species do nc = 1, KCM2 do nh = 1, KRHLEV sumk = f_zero sumok = f_zero sumokg = f_zero sumreft = f_zero do ni = nv1(ib), nv2(ib) sp = sqrt( (f_one - rhdpssa0(ni,nh,nc)) & & /(f_one-rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc))) reft = (f_one - sp) / (f_one + sp) sumreft = sumreft + reft*solwaer(ib,ni) sumk = sumk + rhdpext0(ni,nh,nc)*solwaer(ib,ni) sumok = sumok + rhdpssa0(ni,nh,nc)*solwaer(ib,ni) & & * rhdpext0(ni,nh,nc) sumokg = sumokg + rhdpssa0(ni,nh,nc)*solwaer(ib,ni) & & * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc) enddo refb = sumreft * rsolbd extrhd(nh,nc,ib) = sumk * rsolbd asyrhd(nh,nc,ib) = sumokg / (sumok + 1.0e-10) ssarhd(nh,nc,ib) = 4.0*refb & & /((f_one+refb)**2 - asyrhd(nh,nc,ib)*(f_one-refb)**2) enddo ! end do_nh_block enddo ! end do_nc_block for rh-dep aerosols enddo ! end do_ib_block for sw !!!! open(unit=NOAER,file=aerluts_file,status='UNKNOWN') open(unit=NODMP,file='aeropt_luts.out',status='UNKNOWN') write(NODMP,*) 'Spectral bands = ', NBDSW do ib = 1, NBDSW write(NODMP,*) ib, wvnum1(ib), wvnum2(ib) enddo write(NODMP,*) ' ' write(NODMP,*) 'RH levels = ', KRHLEV write(NODMP,*) (rhlev(k), k = 1, KRHLEV) write(NODMP,*) ' ' write(NODMP,*) 'Aerosol Optical Properties (EXT, ASY, SSA)' write(NODMP,*) 'Aerosol Name = ', AerosolName(1) write(NOAER,'(2x,a6)') AerosolName(1) do ib = 1,NBDSW write(NODMP,'(8f10.5)') (extrhi(ii, ib), ii=1, KCM1) write(NODMP,'(8f10.5)') (asyrhi(ii, ib), ii=1, KCM1) write(NODMP,'(8f10.5)') (ssarhi(ii, ib), ii=1, KCM1) write(NOAER,'(8f10.5)') (extrhi(ii, ib), ii=1, KCM1) ! write(NOAER,'(8f10.5)') (asyrhi(ii, ib), ii=1, KCM1) ! write(NOAER,'(8f10.5)') (ssarhi(ii, ib), ii=1, KCM1) enddo write(NODMP,*) ' ' do i = 1, KCM2 ii = i + 1 write(NODMP,*) 'Aerosol Name = ', AerosolName(ii) write(NOAER,'(2x,a6)') AerosolName(ii) do ib = 1, NBDSW write(NODMP,*) (extrhd(ii,i,ib), ii=1, KRHLEV) write(NODMP,*) (asyrhd(ii,i,ib), ii=1, KRHLEV) write(NODMP,*) (ssarhd(ii,i,ib), ii=1, KRHLEV) write(NOAER,*) (extrhd(ii,i,ib), ii=1, KRHLEV) ! write(NOAER,*) (asyrhd(ii,i,ib), ii=1, KRHLEV) ! write(NOAER,*) (ssarhd(ii,i,ib), ii=1, KRHLEV) ! for soot at 550nm if(ib == 3 ) then if (i == 1) then write(400,*) 'soot' write(400,*) (extrhd(ii,i,ib), ii=1, KRHLEV) endif if (i == 2) then write(400,*) 'suso' write(400,*) (extrhd(ii,i,ib), ii=1, KRHLEV) endif if (i == 3) then write(400,*) 'waso' write(400,*) (extrhd(ii,i,ib), ii=1, KRHLEV) endif if( i == 4) then write(400,*) 'ssam' write(400,*) (extrhd(ii,i,ib), ii=1, KRHLEV) endif if( i == 5) then write(400,*) 'sscm' write(400,*) (extrhd(ii,i,ib), ii=1, KRHLEV) endif endif enddo enddo ! stop end