CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890

      Program CalcTau

      real mass, zenith, dsol, Dobson, lat, lambda, h, min
      real tu, xlon, xlat, asol, phi0, tau_m, p, lambda_cal(10)
      real ozon_ext, V, V0(10), tau_ozon, tau, aot, x, Vd
      integer julday, month, day, year, iday, imon, juldayx
      integer i, j, ndetect, nlambda
      Write(*,*)
      Write(*,*) '---------------------------------------------'
      Write(*,*) '        PROGRAM AOT1.1         '
      Write(*,*) ' Autor Krzysztof Markowicz'
      Write(*,*)
      Write(*,*) 'Program oblicza grubosc optyczna aerozolu'
      Write(*,*) 'Opis programu znajduje sie w pliku Readme.txt'
      Write(*,*) 'Program korzysta z kliku kalibracyjbego lang.dat'
      Write(*,*) '---------------------------------------------'
      Write(*,*)
      Write(*,*) 'Podaj szerokosc geograficzna stacji:'
      Read(*,*) xlat
      Write(*,*) 'Podaj dlugosc geograficzna stacji:'
      Read(*,*) xlon
      Write(*,*) 'Podaj rok pomiaru:'
      Read(*,*) year
      Write(*,*) 'Podaj miesiac pomiaru:'
      Read(*,*) month
      Write(*,*) 'Podaj dzien pomiaru:'
      Read(*,*) day
      Write(*,*) 'Podaj pelna godzine pomiaru w UTC:'
      Read(*,*) h
      Write(*,*) 'Podaj minute pomiaru:'
      Read(*,*) min
      tu=h+min/60
      Write(*,*) 'Podaj cisnienie powietrza na stacji w hpa:'
      Read(*,*) p
      Write(*,*) 'Podaj dlugosc fali w nano-metrach:'
      Read(*,*) lambda
      Write(*,*) 'Podaj napiecie wysciowe z Fotometru w V :'
      Read(*,*) V
      Write(*,*) 'Podaj napiecie zacienionego Fotometru  w V :'
      Read(*,*) Vd
      V=V-Vd
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC read calibraton file
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      open(21,file='lang.dat',status='unknown')
      do i=1,8
       read(21,*)
      end do
      read(21,*) ndetect
      read(21,*) (lambda_cal(j),j=1,ndetect)
      read(21,*) (V0(j),j=1,ndetect)
      close(21)
       nlambda=0
      do i=1,ndetect
       if (lambda.eq.lambda_cal(i)) then
        nlambda=i
       end if
      end do
      if (nlambda.eq.0) then
       write(*,*) 'Niezgodnosc dlugosci fali z plikiem kalibracyjnym'
       write(*,*) 'Sprawdz plik kalibracyjny lang.dat'
       write(*,*) 'Konie programu'
       stop
      end if
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      call possol (month,day,tu,xlon,xlat,asol,phi0)
      write (6,*) 'Kat zenitalny i azymutalny slonca',asol, phi0
      zenith=asol
      call airmass(zenith, mass)
      write (6,*) 'Masa optycza atmosfery',mass

      call varsol(day,month,dsol)

      tau=-(alog(V)-alog(V0(nlambda))-alog(dsol))/mass
      write(6,*) 'Calkowita grubosc optyczna',tau

      call ozon (lambda,ozon_ext)

      call EstimDobson( month, xlat, Dobson )
      write (6,*) 'Klimatologiczna wartosci ozonu w [Du]', Dobson

      tau_ozon=Dobson/1000*ozon_ext
      write(6,*) 'Grubosc optyczna ozonu', tau_ozon

      call rayleigh(lambda,p,tau_m)
      write(6,*) 'Molekularna grubosc optyczna', tau_m

      aot=tau-tau_ozon-tau_m
      write(6,*) '**********************************'
      write(6,*) 'Grubosc optyczna aerozolu',aot
      write(6,*) '**********************************'


      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Get rayleigh optical depth
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      Subroutine rayleigh(lambda,p,tau_m)
       real A, B, C, p0, p, tau_m, lambda,l
       A=8436.e-6
       B=-1225.e-7
       C=14.e-5
       p0=1013
       l=lambda*1.e-3
       tau_m=p/p0*(A*l**-4+B*l**-5+C*l**-6)
       return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Get airmass from solar zenith angle
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Subroutine airmass(zenith,mass)
       real a, b, c, mass, z
       a = 0.50572
       b = 6.07995
       c = 1.6364
       el=90-zenith
       z=3.14159*el/180
       mass = 1.0 /(sin(z) + a * (el+ b)**-c )
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine varsol (jday,month,dsol)
      real dsol,pi,om
      integer jday,month,j end

      if (month.le.2) goto 1
      if (month.gt.8) goto 2
      j=31*(month-1)-((month-1)/2)-2+jday
      goto 3
    1 j=31*(month-1)+jday
      goto 3
    2 j=31*(month-1)-((month-2)/2)-2+jday
    3 pi=2.*acos (0.)
      om=(.9856*float(j-4))*pi/180.
      dsol=1./((1.-.01673*cos(om))**2)
      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

        Subroutine Jul2Greg( julday, year, mon, day )

        Integer*4 julday, year, mon, day
        Integer*4 days(12)

        Logical*1 leap_year

        Data days / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /

        leap_year = .FALSE.
        If ( MOD( year, 4 ) .Eq. 0 )  leap_year = .TRUE.
        If ( leap_year )  days(2) = days(2) + 1
        day = julday
        mon = 1
20      If ( day .GT. days(mon) ) Then
           day = day - days(mon)
           mon = mon + 1
           GoTo 20
        EndIf

        Return
        End
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Estimate Dobson units from climatology, given the month and
CC latitude.  Table has 12 columns, one per month, and 35 rows,
CC from 85 N to -85, in steps of 5 deg lat.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        Subroutine EstimDobson( month, lat, Dobson )

        Integer*4 month, i1, i2

        Real*4 lat, TabDobson(35,12), Dobson

        Data TabDobson /
     o 395, 395, 395, 395, 395, 392, 390, 387, 376,
     o 354, 322, 292, 269, 254, 248, 246, 247, 251,
     o 255, 260, 266, 271, 277, 286, 295, 306, 319,
     o 334, 344, 344, 338, 331, 324, 320, 316,
     o 433, 433, 433, 436, 432, 428, 426, 418, 402,
     o 374, 338, 303, 278, 261, 251, 246, 248, 250,
     o 254, 258, 262, 265, 270, 278, 286, 294, 303,
     o 313, 322, 325, 324, 317, 306, 299, 294,
     o 467, 470, 460, 459, 451, 441, 433, 420, 401,
     o 377, 347, 316, 291, 271, 260, 254, 254, 255,
     o 257, 259, 261, 264, 269, 277, 284, 289, 296,
     o 305, 312, 315, 317, 312, 305, 299, 295,
     o 467, 465, 462, 455, 444, 431, 421, 410, 395,
     o 373, 348, 325, 304, 287, 275, 267, 261, 259,
     o 258, 259, 260, 263, 271, 278, 284, 289, 297,
     o 306, 314, 318, 319, 313, 302, 302, 302,
     o 411, 414, 416, 415, 410, 406, 402, 394, 382,
     o 363, 342, 324, 307, 291, 279, 271, 264, 260,
     o 258, 257, 258, 264, 271, 281, 291, 303, 312,
     o 318, 322, 323, 322, 322, 322, 322, 322,
     o 371, 371, 370, 368, 367, 372, 375, 372, 360,
     o 341, 323, 311, 301, 290, 282, 275, 268, 263,
     o 259, 256, 258, 264, 273, 289, 306, 319, 327,
     o 328, 328, 337, 337, 337, 337, 337, 337,
     o 333, 332, 332, 334, 338, 346, 350, 346, 335,
     o 321, 310, 302, 396, 296, 284, 280, 274, 268,
     o 262, 259, 261, 268, 279, 295, 315, 331, 340,
     o 342, 338, 344, 340, 340, 340, 340, 340,
     o 311, 308, 308, 313, 320, 327, 330, 326, 319,
     o 310, 303, 298, 291, 286, 283, 281, 277, 273,
     o 268, 264, 266, 274, 288, 306, 327, 343, 353,
     o 355, 351, 339, 325, 307, 294, 294, 294,
     o 283, 291, 302, 308, 312, 317, 318, 313, 307,
     o 300, 295, 290, 284, 279, 279, 279, 278, 276,
     o 272, 270, 273, 282, 295, 313, 333, 348, 360,
     o 367, 368, 353, 324, 291, 267, 253, 230,
     o 299, 299, 299, 309, 315, 317, 317, 312, 302,
     o 291, 283, 280, 275, 270, 268, 267, 263, 263,
     o 265, 269, 277, 287, 301, 317, 336, 354, 371,
     o 387, 402, 402, 374, 333, 294, 274, 259,
     o 314, 314, 314, 314, 332, 332, 327, 322, 311,
     o 297, 284, 276, 270, 263, 261, 260, 258, 259,
     o 264, 270, 278, 286, 298, 311, 323, 335, 350,
     o 366, 381, 390, 388, 376, 357, 346, 341,
     o 358, 358, 358, 358, 358, 358, 353, 349, 338,
     o 320, 299, 281, 267, 256, 252, 251, 251, 253,
     o 257, 264, 272, 279, 287, 297, 307, 318, 332,
     o 347, 358, 365, 366, 364, 358, 356, 353 /

      If ( lat .GE. 85.0 ) Then
           Dobson = TabDobson(1,month)
        Else If ( lat .LE. -85.0 ) Then
           Dobson = TabDobson(35,month)
        Else If ( lat .GE. 0.0 ) Then
           i1 = 17 - INT( lat / 5.0 )
           i2 = i1 + 1
           fac = ( TabDobson(i2,month) - TabDobson(i1,month) )
     & / ( 5.0 )
           diffLat = lat - ( 90.0 - (i1 * 5.0) )
           Dobson = TabDobson(i1,month) + fac * diffLat
        Else
           i1 = 18 - INT( lat / 5.0 )
           i2 = i1 + 1
           fac = ( TabDobson(i2,month) - TabDobson(i1,month) )
     & / ( 5.0 )
           diffLat = lat - ( 90.0 - (i1 * 5.0) )
           Dobson = TabDobson(i1,month) + fac * diffLat
        End If

        Return
        End

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Subroutine from 6S package to get Sun position
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine possol (month,jday,tu,xlon,xlat,
     a asol,phi0)

      real tu,xlon,xlat,asol,phi0
      integer month,jday,ia,nojour

c solar position (zenithal angle asol,azimuthal angle phi0
c in degrees)
c jday is the number of the day in the month

      ia = 0
      call day_number(jday,month,ia,nojour)
      call pos_fft (nojour, tu, xlon, xlat, asol, phi0)
      if(asol.gt.90) call print_error(
     s 'The sun is not raised')
      return
      end

      subroutine day_number(jday,month,ia,j)
      integer jday, month, ia, j

      if (month.le.2) then
                      j=31*(month-1)+jday
                      return
                      endif
      if (month.gt.8) then
                      j=31*(month-1)-((month-2)/2)-2+jday
                      else
                      j=31*(month-1)-((month-1)/2)-2+jday
                      endif
      if(ia.ne.0 .and. mod(ia,4).eq.0) j=j+1
      return
      end

      subroutine pos_fft (j,tu,xlon,xlat,asol,phi0)
      real tu, xlat, asol,phi0, tsm, xlon,xla, xj, tet,
     a a1, a2, a3, a4, a5, et, tsv, ah, b1, b2, b3, b4,
     a b5, b6, b7, delta, amuzero, elev, az, caz, azim, pi2
      integer j
      parameter (pi=3.14159265,fac=pi/180.)
c solar position (zenithal angleasol,azimuthal angle phi0
c in degrees)
c j is the day number in the year
c
c mean solar time (heure decimale)

      tsm=tu+xlon/15.
      xla=xlat*fac
      xj=float(j)
      tet=2.*pi*xj/365.

c time equation (in mn.dec)
      a1=.000075
      a2=.001868
      a3=.032077
      a4=.014615
      a5=.040849
      et=a1+a2*cos(tet)-a3*sin(tet)-a4*cos(2.*tet)-a5*sin(2.*tet)
      et=et*12.*60./pi

c true solar time

      tsv=tsm+et/60.
      tsv=(tsv-12.)

c hour angle
      ah=tsv*15.*fac

c solar declination (in radian)
      b1=.006918
      b2=.399912
      b3=.070257
      b4=.006758
      b5=.000907
      b6=.002697
      b7=.001480
      delta=b1-b2*cos(tet)+b3*sin(tet)-b4*cos(2.*tet)+b5*sin(2.*tet)-
     &b6*cos(3.*tet)+b7*sin(3.*tet)

c elevation,azimuth

      amuzero=sin(xla)*sin(delta)+cos(xla)*cos(delta)*cos(ah)
      elev=asin(amuzero)
      az=cos(delta)*sin(ah)/cos(elev)
      if ( (abs(az)-1.000).gt.0.00000) az = sign(1.,az)
      caz=(-cos(xla)*sin(delta)+sin(xla)*cos(delta)*cos(ah))/cos(elev)
      azim=asin(az)
      if(caz.le.0.) azim=pi-azim
      if(caz.gt.0.and.az.le.0) azim=2*pi+azim
      azim=azim+pi
      pi2=2*pi
      if(azim.gt.pi2) azim=azim-pi2
      elev=elev*180./pi
c conversion in degrees
      asol=90.-elev
      phi0=azim/fac
      return
      end

      subroutine print_error(tex)
      character *(*) tex
      logical ier
      integer iwr
c JM 6/99 common/sixs_ier/iwr,ier
      iwr = 6
      ier = .TRUE.
      write(iwr,'(a)')tex
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Subroutine julianday(year,month,day,jd)

      integer year, month, day, jd
      logical*1 leap_year
      leap_year=.FALSE.
      if (mod(year,4).eq.0) leap_year=.TRUE.

      if (month.eq.1) jd=day
      if (month.eq.2) jd=31+day

      if (leap_year) day=day+1
      if (month.eq.3) jd=31+28+day
      if (month.eq.4) jd=31+28+31+day
      if (month.eq.5) jd=31+28+31+30+day
      if (month.eq.6) jd=31+28+31+30+31+day
      if (month.eq.7) jd=31+28+31+30+31+30+day
      if (month.eq.8) jd=31+28+31+30+31+30+31+day
      if (month.eq.9) jd=31+28+31+30+31+30+31+31+day
      if (month.eq.10) jd=31+28+31+30+31+30+31+31+30+day
      if (month.eq.11) jd=31+28+31+30+31+30+31+31+30+31+day
      if (month.eq.12) jd=31+28+31+30+31+30+31+31+30+31+30+day
      if (leap_year) day=day-1
      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC get Ozon extinction coefition
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C234567890123456789012345678901234567890123456789012345678901234567890
      Subroutine  ozon(lambda, ozon_ext)

       real tabela(2,52), lam(52), ext(52), l, ozon_ext, lambda
       integer ni, i
       Data tabela /
     $ 0.445, 0.003, 0.450, 0.003, 0.455, 0.004, 0.460, 0.006, 0.465,
     $ 0.008, 0.470, 0.009, 0.475, 0.012, 0.480, 0.014, 0.485, 0.017,
     $ 0.490, 0.021, 0.495, 0.025, 0.500, 0.030, 0.505, 0.035, 0.510,
     $ 0.040, 0.515, 0.045, 0.520, 0.048, 0.525, 0.057, 0.530, 0.063,
     $ 0.535, 0.070, 0.540, 0.075, 0.545, 0.080, 0.550, 0.085, 0.555,
     $ 0.095, 0.560, 0.103, 0.565, 0.110, 0.570, 0.120, 0.575, 0.122,
     $ 0.580, 0.120, 0.585, 0.118, 0.590, 0.115, 0.595, 0.120, 0.600,
     $ 0.125, 0.605, 0.130, 0.610, 0.120, 0.620, 0.105, 0.630, 0.090,
     $ 0.640, 0.079, 0.650, 0.067, 0.660, 0.057, 0.670, 0.048, 0.680,
     $ 0.036, 0.690, 0.028, 0.700, 0.023, 0.710, 0.018, 0.720, 0.014,
     $ 0.730, 0.011, 0.740, 0.010, 0.750, 0.009, 0.760, 0.007, 0.770,
     $ 0.004, 0.780, 0.000, 0.790, 0.000 /

       l=lambda*1.e-3
       ni=52
       do i=1,ni
        lam(i)=tabela(1,i)
        ext(i)=tabela(2,i)
       end do
      do i=1,ni-1
       if ((lam(i).lt.l).and.(lam(i+1).ge.l)) then
        ozon_ext=ext(i)+(ext(i+1)-ext(i))/(lam(i+1)-lam(i))*
     $ (l-lam(i))
       end if
      end do
      if (l.lt.0.445) then
       ozon_ext=0
       Write(6,*) 'Dlugosc fali mniejsza od 445 nm'
       Write(6,*) 'przyjeto grubosc optyczna ozonu rowna O'
      end if
      if (l.gt.0.790) then
       ozon_ext=0
       Write(6,*) 'Dlugosc fali wieksza od 790nm'
       Write(6,*) 'przyjeto grubosc optyczna ozonu rowna O'
      end if
      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC