source: trunk/src/zenithangle.f90 @ 28

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