ZENSUN.f

Go to the documentation of this file.

 
!-------------------------------------------------------------------------------------
      subroutine zensun(day,time,lat,lon,pi,sun_zenith,sun_azimuth)
 
!
  use kinds, only: r_kind,i_kind
 
  implicit none
 
  integer(i_kind),intent(in)::   day
  real(r_kind),intent(in) ::  pi,time,lat,lon
  integer(i_kind)   di
  real(r_kind)      z,gaz
  real(r_kind)      ut,noon
  real(r_kind)      y(5),y2(5),x(2,5),Tx(5,2),xTx(2,2),aTx(5,2),det
  real(r_kind)      tt,eqtime,decang,latsun,lonsun
  real(r_kind)      nday(74),eqt(74),dec(74)
  real(r_kind)      beta(2), beta2(2), a(2,2)
  real(r_kind)      t0,t1,p0,p1,zz,xx,yy
  real(r_kind),intent(out) ::     sun_zenith,sun_azimuth
 
  data   nday/1.0,   6.0,  11.0,  16.0,  21.0,  26.0,  31.0,  36.0,  41.0,  46.0,&
             51.0,  56.0,  61.0,  66.0,  71.0,  76.0,  81.0,  86.0,  91.0,  96.0,&
             101.0, 106.0, 111.0, 116.0, 121.0, 126.0, 131.0, 136.0, 141.0, 146.0,&
             151.0, 156.0, 161.0, 166.0, 171.0, 176.0, 181.0, 186.0, 191.0, 196.0,&
             201.0, 206.0, 211.0, 216.0, 221.0, 226.0, 231.0, 236.0, 241.0, 246.0,&
             251.0, 256.0, 261.0, 266.0, 271.0, 276.0, 281.0, 286.0, 291.0, 296.0,&
             301.0, 306.0, 311.0, 316.0, 321.0, 326.0, 331.0, 336.0, 341.0, 346.0,&
             351.0, 356.0, 361.0, 366.0/
 
  data  eqt/ -3.23, -5.49, -7.60, -9.48,-11.09,-12.39,-13.34,-13.95,-14.23,-14.19,&
            -13.85,-13.22,-12.35,-11.26,-10.01, -8.64, -7.18, -5.67, -4.16, -2.69,&
             -1.29, -0.02,  1.10,  2.05,  2.80,  3.33,  3.63,  3.68,  3.49,  3.09,&
              2.48,  1.71,  0.79, -0.24, -1.33, -2.41, -3.45, -4.39, -5.20, -5.84,&
             -6.28, -6.49, -6.44, -6.15, -5.60, -4.82, -3.81, -2.60, -1.19,  0.36,&
              2.03,  3.76,  5.54,  7.31,  9.04, 10.69, 12.20, 13.53, 14.65, 15.52,&
             16.12, 16.41, 16.36, 15.95, 15.19, 14.09, 12.67, 10.93,  8.93,  6.70,&
              4.32,  1.86, -0.62, -3.23/
 
  data dec/ -23.06,-22.57,-21.91,-21.06,-20.05,-18.88,-17.57,-16.13,-14.57,-12.91,&
            -11.16, -9.34, -7.46, -5.54, -3.59, -1.62,  0.36,  2.33,  4.28,  6.19,&
              8.06,  9.88, 11.62, 13.29, 14.87, 16.34, 17.70, 18.94, 20.04, 21.00,&
             21.81, 22.47, 22.95, 23.28, 23.43, 23.40, 23.21, 22.85, 22.32, 21.63,&
             20.79, 19.80, 18.67, 17.42, 16.05, 14.57, 13.00, 11.33,  9.60,  7.80,&
              5.95,  4.06,  2.13,  0.19, -1.75, -3.69, -5.62, -7.51, -9.36,-11.16,&
            -12.88,-14.53,-16.07,-17.50,-18.81,-19.98,-20.99,-21.85,-22.52,-23.02,&
            -23.33,-23.44,-23.35,-23.06/
 
!
! compute the subsolar coordinates
!
 
  tt= mod(real((int(day)+time/24.-1.)),365.25) +1.  ! fractional day number
                                                    ! with 12am 1jan = 1.
  do di = 1, 73
    if ((tt >= nday(di)) .and. (tt <= nday(di+1))) exit
  end do
 
!============== Perform a least squares regression on doy**3 ==============
 
  x(1,:) = 1.0
 
  if ((di >= 3) .and. (di <= 72)) then
    y(:) = eqt(di-2:di+2)
    y2(:) = dec(di-2:di+2)
 
    x(2,:) = nday(di-2:di+2)**3
  end if
  if (di == 2) then
    y(1) = eqt(73)
    y(2:5) = eqt(di-1:di+2)
    y2(1) = dec(73)
    y2(2:5) = dec(di-1:di+2)
 
    x(2,1) = nday(73)**3
    x(2,2:5) = (365.+nday(di-1:di+2))**3
  end if
  if (di == 1) then
    y(1:2) = eqt(72:73)
    y(3:5) = eqt(di:di+2)
    y2(1:2) = dec(72:73)
    y2(3:5) = dec(di:di+2)
 
    x(2,1:2) = nday(72:73)**3
    x(2,3:5) = (365.+nday(di:di+2))**3
  end if
  if (di == 73) then
    y(1:4) = eqt(di-2:di+1)
    y(5) = eqt(2)
    y2(1:4) = dec(di-2:di+1)
    y2(5) = dec(2)
 
    x(2,1:4) = nday(di-2:di+1)**3
    x(2,5) = (365.+nday(2))**3
  end if
  if (di == 74) then
    y(1:3) = eqt(di-2:di)
    y(4:5) = eqt(2:3)
    y2(1:3) = dec(di-2:di)
    y2(4:5) = dec(2:3)
 
    x(2,1:3) = nday(di-2:di)**3
    x(2,4:5) = (365.+nday(2:3))**3
  end if
 
!  Tx = transpose(x)
  tx(1:5,1)=x(1,1:5)
  tx(1:5,2)=x(2,1:5)
!  xTx = MATMUL(x,Tx)
  xtx(1,1)=x(1,1)*tx(1,1)+x(1,2)*tx(2,1)+x(1,3)*tx(3,1)+x(1,4)*tx(4,1)+x(1,5)*tx(5,1)
  xtx(1,2)=x(1,1)*tx(1,2)+x(1,2)*tx(2,2)+x(1,3)*tx(3,2)+x(1,4)*tx(4,2)+x(1,5)*tx(5,2)
  xtx(2,1)=x(2,1)*tx(1,1)+x(2,2)*tx(2,1)+x(2,3)*tx(3,1)+x(2,4)*tx(4,1)+x(2,5)*tx(5,1)
  xtx(2,2)=x(2,1)*tx(1,2)+x(2,2)*tx(2,2)+x(2,3)*tx(3,2)+x(2,4)*tx(4,2)+x(2,5)*tx(5,2)
 
  det = xtx(1,1)*xtx(2,2) - xtx(1,2)*xtx(2,1)
  a(1,1) = xtx(2,2)/det
  a(1,2) = -xtx(1,2)/det
  a(2,1) = -xtx(2,1)/det
  a(2,2) = xtx(1,1)/det
 
!  aTx = MATMUL(Tx,a)
  atx(1,1)=tx(1,1)*a(1,1)+tx(1,2)*a(2,1) 
  atx(2,1)=tx(2,1)*a(1,1)+tx(2,2)*a(2,1) 
  atx(3,1)=tx(3,1)*a(1,1)+tx(3,2)*a(2,1) 
  atx(4,1)=tx(4,1)*a(1,1)+tx(4,2)*a(2,1) 
  atx(5,1)=tx(5,1)*a(1,1)+tx(5,2)*a(2,1) 
  atx(1,2)=tx(1,1)*a(1,2)+tx(1,2)*a(2,2) 
  atx(2,2)=tx(2,1)*a(1,2)+tx(2,2)*a(2,2) 
  atx(3,2)=tx(3,1)*a(1,2)+tx(3,2)*a(2,2) 
  atx(4,2)=tx(4,1)*a(1,2)+tx(4,2)*a(2,2) 
  atx(5,2)=tx(5,1)*a(1,2)+tx(5,2)*a(2,2) 
 
!  beta = MATMUL(y,aTx)
  beta(1) = y(1)*atx(1,1)+y(2)*atx(2,1)+y(3)*atx(3,1)+y(4)*atx(4,1)+y(5)*atx(5,1)
  beta(2) = y(1)*atx(1,2)+y(2)*atx(2,2)+y(3)*atx(3,2)+y(4)*atx(4,2)+y(5)*atx(5,2)
 
!  beta2 = MATMUL(y2,aTx)
  beta2(1) = y2(1)*atx(1,1)+y2(2)*atx(2,1)+y2(3)*atx(3,1)+y2(4)*atx(4,1)+y2(5)*atx(5,1)
  beta2(2) = y2(1)*atx(1,2)+y2(2)*atx(2,2)+y2(3)*atx(3,2)+y2(4)*atx(4,2)+y2(5)*atx(5,2)
 
!============== finished least squares regression on doy**3 ==============
 
  if ((di < 3) .or. (di > 72)) tt = tt + 365.
 
  eqtime=(beta(1) + beta(2)*tt**3)/60.
  decang=beta2(1) + beta2(2)*tt**3
  latsun=decang
 
  ut=time
  noon=12.-lon/15.                             ! universal time of noon
 
  lonsun=-15.*(ut-12.+eqtime)
 
  t0=(90.-lat)*pi/180.
  t1=(90.-latsun)*pi/180.
 
  p0=lon*pi/180.
  p1=lonsun*pi/180.
 
  zz=cos(t0)*cos(t1)+sin(t0)*sin(t1)*cos(p1-p0)
!  zz2=sin(t0)*sin(t1)+cos(t0)*cos(t1)*cos(p1-p0)
  xx=sin(t1)*sin(p1-p0)
  yy=sin(t0)*cos(t1)-cos(t0)*sin(t1)*cos(p1-p0)
 
  sun_zenith=90-acos(zz)/(pi/180)
  sun_azimuth=atan2(xx,yy)/(pi/180)
 
  return
end subroutine zensun