! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later

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