sunpos.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE sunpos (KYEAR, KMONTH, KDAY, PTIME, &
                         PLON, PLAT, PTSUN, PZENITH, PAZIMSOL)
!     ####################################################################################
!
!!****  *SUNPOS * - routine to compute the position of the sun
!!
!!    PURPOSE
!!    -------
!!      The purpose of this routine is to compute the cosine and sinus of the 
!!    solar zenithal angle (angle defined by the local vertical at the position
!!    XLAT, XLON and the direction of the sun) and the azimuthal solar
!!    angle (angle between an horizontal direction (south or north according
!!    to the terrestrial hemisphere) and the horizontal projection of the
!!    direction of the sun.
!!
!!**  METHOD
!!    ------
!!      The cosine and sinus of the zenithal solar angle  and the azimuthal 
!!    solar angle are computed from the true universal time, valid for the (XLAT,
!!    XLON) location, and from the solar declination aPD_ICE(:)ngle of the day. There
!!    is a special convention to define the azimuthal solar angle.
!!     
!!    EXTERNAL
!!    --------
!!      NONE
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!      "Radiative Processes in Meteorology and Climatology"  
!!                          (1976)   Paltridge and Platt 
!!
!!    AUTHOR
!!    ------
!!      J.-P. Pinty      * Laboratoire d'Aerologie*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original             16/10/94 
!!      Revised              12/09/95
!!      (J.Stein)            01:04/96  bug correction for ZZEANG     
!!      (K. Suhre)           14/02/97  bug correction for ZLON0     
!!      (V. Masson)          01/03/03  add zenithal angle output
!!      (V. Masson)          14/03/14  avoid discontinuous declination at 00UTC each day
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_csts,          ONLY : xpi, xday
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef AIX64
!$ USE OMP_LIB
#endif
!
IMPLICIT NONE
!
#ifndef AIX64
!$ INCLUDE 'omp_lib.h'
#endif
!
!*       0.1   Declarations of dummy arguments :
!
INTEGER,                      INTENT(IN)   :: KYEAR      ! current year                        
INTEGER,                      INTENT(IN)   :: KMONTH     ! current month                        
INTEGER,                      INTENT(IN)   :: KDAY       ! current day                        
REAL,                         INTENT(IN)   :: PTIME      ! current time                        
REAL, DIMENSION(:),           INTENT(IN)   :: PLON       ! longitude
REAL, DIMENSION(:),           INTENT(IN)   :: PLAT       ! latutude
!
REAL, DIMENSION(:),           INTENT(OUT)  :: PZENITH    ! Solar zenithal angle
REAL, DIMENSION(:),           INTENT(OUT)  :: PAZIMSOL   ! Solar azimuthal angle
REAL, DIMENSION(:),           INTENT(OUT)  :: PTSUN      ! Solar time
!
!*       0.2   declarations of local variables
!
!
REAL                                       :: ZUT        ! Universal time
!
REAL, DIMENSION(SIZE(PLON))                :: ZTUT    ,&! True (absolute) Universal Time
                                                ZSOLANG ,&! Hourly solar angle
                                                ZSINAZI ,&! Sine of the solar azimuthal angle
                                                ZCOSAZI ,&! Cosine of the solar azimuthal angle
                                                ZLAT,    &
                                                ZLON,    &! Array of latitudes and longitudes
                                                ZSINZEN, &!Sine of zenithal angle
                                                ZCOSZEN   !Cosine of zenithal angle  
INTEGER, DIMENSION(0:11)                   :: IBIS, INOBIS ! Cumulative number of days per month
                                                           ! for bissextile and regular years
REAL                                       :: ZDATE         ! Julian day of the year
REAL                                       :: ZAD           ! Angular Julian day of the year
REAL                                       :: ZDECSOL       ! Daily solar declination angle 
REAL                                       :: ZA1, ZA2      ! Ancillary variables
REAL                                       :: ZTSIDER, &
                                                ZSINDEL, &!azimuthal angle
                                                ZCOSDEL !azimuthal angle  
!                                            
INTEGER                                    :: JI, JJ, INKPROMA
INTEGER    :: IINDX1, IINDX2
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!
!-------------------------------------------------------------------------------
!
!*       1.    TO COMPUTE THE TRUE SOLAR TIME
!              -------------------------------
!
IF (lhook) CALL dr_hook('SUNPOS_1',0,zhook_handle)
!
zut  = mod( 24.0+mod(ptime/3600.,24.0),24.0 )

inobis(:) = (/0,31,59,90,120,151,181,212,243,273,304,334/)
ibis(0:1) = inobis(0:1)
DO ji=2,11
  ibis(ji) = inobis(ji)+1
END DO
IF( mod(kyear,4).EQ.0 .AND. (mod(kyear,100).NE.0 .OR. mod(kyear,400).EQ.0)) THEN
  zdate = float(kday +   ibis(kmonth-1)) - 1 + ptime/xday
  zad = 2.0*xpi*zdate/366.0
ELSE
  zdate = float(kday + inobis(kmonth-1)) - 1 + ptime/xday
  zad = 2.0*xpi*zdate/365.0
END IF

za1 = (1.00554*zdate- 6.28306)*(xpi/180.0)
za2 = (1.93946*zdate+23.35089)*(xpi/180.0)
ztsider = (7.67825*sin(za1)+10.09176*sin(za2)) / 60.0
!
!-------------------------------------------------------------------------------
!
!*       2.     COMPUTE THE SOLAR DECLINATION ANGLE
!               -----------------------------------
!
zdecsol = 0.006918-0.399912*cos(zad)   +0.070257*sin(zad)    &
           -0.006758*cos(2.*zad)+0.000907*sin(2.*zad) &
           -0.002697*cos(3.*zad)+0.00148 *sin(3.*zad)  
zsindel = sin(zdecsol)
zcosdel = cos(zdecsol)
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('SUNPOS_1',1,zhook_handle)
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP) 
IF (lhook) CALL dr_hook('SUNPOS_2',0,zhook_handle_omp)
!$OMP DO PRIVATE(JJ)
DO jj = 1,SIZE(plat)
!
!*       3.    LOADS THE ZLAT, ZLON ARRAYS
!              ---------------------------
!
  zlat(jj) = plat(jj)*(xpi/180.)
  zlon(jj) = plon(jj)*(xpi/180.)
!
!-------------------------------------------------------------------------------
!
!*       4.    COMPUTE THE TRUE SOLAR TIME
!              ----------------------------
!
  ztut(jj) = zut - ztsider + zlon(jj)*((180./xpi)/15.0)
!
  ptsun(jj) = mod(ptime -ztsider*3600. +plon(jj)*240., xday)
!
!-------------------------------------------------------------------------------
!*       3.    COMPUTES THE COSINE AND SINUS OF THE ZENITHAL SOLAR ANGLE
!              ---------------------------------------------------------
!
  zsolang(jj) = (ztut(jj)-12.0)*15.0*(xpi/180.)          ! hour angle in radians
!
  zcoszen(jj) = sin(zlat(jj))*zsindel +                 &! Cosine of the zenithal
                 cos(zlat(jj))*zcosdel*cos(zsolang(jj))  !       solar angle  
!
  zsinzen(jj)  = sqrt( 1. - zcoszen(jj)*zcoszen(jj) )
!
!-------------------------------------------------------------------------------
!
!*       5.    ZENITHAL SOLAR ANGLE
!              --------------------
!
  pzenith(jj) = acos(zcoszen(jj))
!
!-------------------------------------------------------------------------------
!
!*       6.    COMPUTE THE AZIMUTHAL SOLAR ANGLE (PAZIMSOL)
!              --------------------------------------------
!
  IF (zsinzen(jj)/=0.) THEN
    !Azimuth is measured clockwise from north
    zsinazi(jj)  = - zcosdel * sin(zsolang(jj)) / zsinzen(jj)
    zcosazi(jj)  = (-sin(zlat(jj))*zcosdel*cos(zsolang(jj))      &
                       +cos(zlat(jj))*zsindel                       &
                      ) / zsinzen(jj)  
    pazimsol(jj) = atan2(zsinazi(jj),zcosazi(jj))
  ELSE
    pazimsol(jj) = xpi
  ENDIF
!
ENDDO
!$OMP END DO 
IF (lhook) CALL dr_hook('SUNPOS_2',1,zhook_handle_omp)
!$OMP END PARALLEL
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE sunpos