1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
real function zenithangle(ylat,xlon,jul)
!*********************************************************************
! *
! Author: G. WOTAWA *
! Date: 1993-11-17 *
! Project: POP-M *
! Last update: *
! *
!*********************************************************************
! *
! DESCRIPTION: This function returns the sinus of solar *
! elevation as a function of geographic longitude, *
! latitude and GMT-Time. *
! *
!*********************************************************************
! *
! INPUT: *
! *
! ylat geographical latitude [DEG] *
! xlon geographical longitude [DEG] *
! jjjj Year *
! mm Month *
! dd Day *
! hh Hour *
! minute Minute *
! *
!*********************************************************************
use par_mod, only: dp
implicit none
integer :: jjjj,mm,id,iu,minute,yyyymmdd,hhmmss
integer :: ndaynum
real :: sinsol,solelev,ylat,xlon
real :: rnum,rylat,ttime,dekl,rdekl,eq
real,parameter :: pi=3.1415927
real(kind=dp) :: jul
call caldate(jul,yyyymmdd,hhmmss)
jjjj=yyyymmdd/10000
mm=yyyymmdd/100-jjjj*100
id=yyyymmdd-jjjj*10000-mm*100
iu=hhmmss/10000
minute=hhmmss/100-100*iu
ndaynum=31*(mm-1)+id
if(mm.gt.2) ndaynum=ndaynum-int(0.4*mm+2.3)
if((mm.gt.2).and.(jjjj/4*4.eq.jjjj)) ndaynum=ndaynum+1
rnum=2.*pi*ndaynum/365.
rylat=pi*ylat/180.
ttime=real(iu)+real(minute)/60.
dekl=0.396+3.631*sin(rnum)+0.038*sin(2.*rnum)+0.077*sin(3.*rnum)- &
22.97*cos(rnum)-0.389*cos(2.*rnum)-0.158*cos(3.*rnum)
rdekl=pi*dekl/180.
eq=(0.003-7.343*sin(rnum)-9.47*sin(2.*rnum)- &
0.329*sin(3.*rnum)-0.196*sin(4.*rnum)+ &
0.552*cos(rnum)-3.020*cos(2.*rnum)- &
0.076*cos(3.*rnum)-0.125*cos(4.*rnum))/60.
sinsol=sin(rylat)*sin(rdekl)+cos(rylat)*cos(rdekl)* &
cos((ttime-12.+xlon/15.+eq)*pi/12.)
! Calculate the maximum solar elevation on that day
!sinsol=sin(rylat)*sin(rdekl)+cos(rylat)*cos(rdekl)*
! & cos((eq)*pi/12.)
solelev=asin(sinsol)*180./pi
zenithangle=90.-solelev
return
end function zenithangle
|