! =================================================================================================================================
! MODULE       : solar
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2011)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF        This module define solar 
!!
!!\n DESCRIPTION: 
!!
!! RECENT CHANGE(S): None
!!
!! REFERENCE(S)	: None
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_global/solar.f90 $
!! $Date: 2017-10-18 11:15:06 +0200 (mer. 18 oct. 2017) $
!! $Revision: 4693 $
!! \n
!_ ================================================================================================================================

MODULE solar

  USE defprec
  USE constantes
  USE ioipsl_para 

  IMPLICIT NONE



CONTAINS


!! ================================================================================================================================
!! SUBROUTINE   : solarang
!!
!>\BRIEF         This subroutine computes the solar angle according to the method used by GSWP and developed by J.C. Morill.
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::csang
!!
!! REFERENCE(S) : See for further details :
!! http://www.atmo.arizona.edu/~morrill/swrad.html
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE solarang (julian, julian0, iim, jjm, lon, lat, csang)

    USE calendar
    USE time, ONLY : one_year, calendar_str

    IMPLICIT NONE

    !! 0. Parameters and variables declaration
    
    !! 0.1 Input variables

    INTEGER, INTENT(in)                      :: iim, jjm        !!
    REAL, INTENT(in)                         :: julian          !!
    REAL, INTENT(in)                         :: julian0         !!
    REAL, DIMENSION(iim,jjm), INTENT(in)     :: lon, lat        !!

    !! 0.2 Output variables

    REAL, DIMENSION(iim,jjm), INTENT(out)    :: csang           !!

    !! 0.4 Local variables

    REAL                                     :: gamma           !!
    REAL                                     :: dec             !!
    REAL                                     :: decd            !!
    REAL                                     :: et              !!
    REAL                                     :: gmt             !!
    REAL                                     :: le              !!
    REAL                                     :: ls              !!
    REAL                                     :: lcorr           !!
    REAL                                     :: latime          !!
    REAL                                     :: omega           !!
    REAL                                     :: omegad          !!
    REAL                                     :: llatd, llat     !!
    INTEGER                                  :: igmt            !!
    INTEGER                                  :: ilon, ilat      !!
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: zone            !!
!$OMP THREADPRIVATE(zone)
    REAL, SAVE, ALLOCATABLE, DIMENSION(:)    :: lhour           !!
!$OMP THREADPRIVATE(lhour)

!_ ================================================================================================================================

    IF (printlev >= 4) WRITE(numout,*) 'Calendar information', julian, julian0
    !-
    !- 1) Day angle gamma
    !-
    !   gamma = 2.*pi*MOD(julian,one_year)/one_year
    gamma = 2.*pi*(julian-julian0)/one_year
    !-
    !- 2) Solar declination (assumed constant for a 24 hour period)  page 7
    !-    in radians
    !-
    IF (printlev >= 4) WRITE(numout,*) 'Solar declination', gamma, one_year
    !
    dec = ( 0.006918-0.399912*COS(gamma)+0.070257*SIN(gamma) &
         &       -0.006758*COS(2*gamma)+0.000907*SIN(2*gamma)      &
         &       -0.002697*COS(3*gamma)+0.00148*SIN(3*gamma))
    decd = dec*(180/pi)
    !-
    !- maximum error 0.0006 rad (<3'), leads to error
    !- of less than 1/2 degree in zenith angle
    !-
    IF (printlev >= 4) WRITE(numout,*) 'Equation of time', decd
    !- 3)  Equation of time  page 11
    !-
    et = ( 0.000075+0.001868*COS(gamma)-0.032077*SIN(gamma)&
         &      -0.014615*COS(2*gamma)-0.04089*SIN(2*gamma))*229.18
    !-
    !- Get the time zones for the current time
    !-
    gmt = 24.*(julian-INT(julian))
    IF (.NOT.ALLOCATED(zone))  ALLOCATE(zone(iim))
    IF (.NOT.ALLOCATED(lhour)) ALLOCATE(lhour(iim))
    !-
    !igmt = NINT(gmt)
    IF (printlev >= 4) WRITE(numout,*) 'Get time zone', et, gmt
    CALL time_zone(gmt, iim, jjm, lon, zone, lhour)
    !-
    !- Start a loop through the grid
    !-
    IF (printlev >= 4) WRITE(numout,*) 'Start a loop through the grid'
    DO ilon=1,iim
       !---
       !--- 4) Local apparent time  page 13
       !---
       !--- ls     standard longitude (nearest 15 degree meridian)
       !--- le     local longtitude
       !--- lhour  local standard time
       !--- latime local apparent time
       !--- lcorr  longitudunal correction (minutes)
       !---
       le = lon(ilon,1)
       ls = ((zone(ilon)-1)*15)-180.
       lcorr = 4.*(ls-le)*(-1)
       latime = lhour(ilon)+lcorr/60.+et/60.
       IF (latime <  0.) latime = latime+24
       IF (latime > 24.) latime = latime-24
       !---
       !--- 5) Hour angle omega  page 15
       !---
       !--- hour angle is zero at noon, positive in the morning
       !--- It ranges from 180 to -180
       !---
       !--- omegad is omega in degrees, omega is in radians
       !---
       omegad = (latime-12.)*(-15.)
       omega  = omegad*pi/180.
       !---
       DO ilat=1,jjm
          !-----
          !----- 6)  Zenith angle  page 15
          !-----
          !----- csang cosine of zenith angle (radians)
          !----- llatd =  local latitude in degrees
          !----- llat  =  local latitude in radians
          !-----
          llatd = lat(1,ilat)
          llat  = llatd*pi/180.
          csang(ilon,ilat) = &
               &  MAX(zero,SIN(dec)*SIN(llat)+COS(dec)*COS(llat)*COS(omega))
       ENDDO
    ENDDO
    !----------------------
  END SUBROUTINE solarang


!! ================================================================================================================================
!! SUBROUTINE   : time_zone
!!
!>\BRIEF         
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::zone, ::lhour
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE time_zone (gmt, iim, jjm, lon, zone, lhour)

    IMPLICIT NONE

    !! 0. Parameters and variables declaration
    
    !! 0.1 Input variables

    INTEGER, INTENT(in)                   :: iim, jjm        !!
    REAL, DIMENSION(iim,jjm), INTENT(in)  :: lon             !!
    REAL, INTENT(in)                      :: gmt             !!

    !! 0.2 Output variables

    INTEGER, DIMENSION(iim), INTENT(out)  :: zone            !!
    REAL, DIMENSION(iim), INTENT(out)     :: lhour           !!

    !! 0.4 Local variables

    INTEGER                               :: deg             !!
    !!??   REAL :: deg
    INTEGER                               :: i, ilon, change !! Indices

!_ ================================================================================================================================

    DO ilon=1,iim
       !---
       !--- Convert longitude index to longtitude in degrees
       !---
       deg = lon(ilon,1)
       !---
       !--- Determine into which time zone (15 degree interval) the
       !--- longitude falls.
       !---
       DO i=1,25
          IF (deg < (-187.5+(15*i))) THEN
             zone(ilon) = i
             IF (zone(ilon) == 25)   zone(ilon) = 1
             EXIT
          ENDIF
       ENDDO
       !---
       !--- Calculate change (in number of hours) from GMT time to
       !--- local hour.  Change will be negative for zones < 13 and
       !--- positive for zones > 13.
       !---
       !--- There is also a correction for lhour < 0 and lhour > 23
       !--- to lhour between 0 and 23.
       !---
       IF (zone(ilon) < 13) THEN
          change = zone(ilon)-13
          lhour(ilon) = gmt+change
       ENDIF
       !---
       IF (zone(ilon) == 13) THEN
          lhour(ilon) = gmt
       ENDIF
       !---
       IF (zone(ilon) > 13) THEN
          change = zone(ilon)-13
          lhour(ilon) = gmt+change
       ENDIF
       IF (lhour(ilon) <  0) lhour(ilon) = lhour(ilon)+24
       IF (lhour(ilon) >= 24) lhour(ilon) = lhour(ilon)-24
       !---
    ENDDO
    !-----------------------
  END SUBROUTINE time_zone



!! ================================================================================================================================
!! SUBROUTINE   : downward_solar_flux
!!
!>\BRIEF         
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::solad, ::solai
!!
!! REFERENCE(S) :See for further details :
!! http://www.atmo.arizona.edu/~morrill/swrad.html
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE downward_solar_flux (npoi,latitude,jday,rtime,cloud,nband,solad,solai)

    USE time, ONLY : one_year, calendar_str
    
    IMPLICIT NONE

    !! 0. Parameter and variables declaration
    
    !! 0.1 Input variables
    
    INTEGER, INTENT(in)                      :: npoi               !! number of grid points (unitless)
    REAL, DIMENSION(npoi), INTENT(in)        :: latitude           !! latitude in degrees
    REAL, INTENT(in)                         :: jday
    REAL,INTENT(in)                          :: rtime
    REAL, DIMENSION(npoi), INTENT(in)        :: cloud              !! cloud fraction (0-1, unitless)
    INTEGER,INTENT(in)                       :: nband              !! number of visible bands (unitless)

    !! 0.2 Output variables

    REAL, DIMENSION(npoi,nband), INTENT(out) :: solad              !! direct downward solar flux 
                                                                   !! @tex $(W.m^{-2})$ @endtex
    REAL, DIMENSION(npoi,nband), INTENT(out) :: solai              !! diffuse downward solar flux 
                                                                   !! @tex $(W.m^{-2})$ @endtex

    !! 0.4 Local variables

    ! Parametres orbitaux:

    REAL, SAVE                               :: ecc                !! Eccentricity
!$OMP THREADPRIVATE(ecc)
    REAL, SAVE                               :: perh               !! Longitude of perihelie
!$OMP THREADPRIVATE(perh)
    REAL, SAVE                               :: xob                !! obliquity
!$OMP THREADPRIVATE(xob)
    REAL, PARAMETER                          :: pir = pi/180.      !!
    REAL, SAVE                               :: step               !!
!$OMP THREADPRIVATE(step)
    REAL                                     :: xl                 !!
    REAL                                     :: so                 !!
    REAL                                     :: xllp               !!
    REAL                                     :: xee                !!
    REAL                                     :: xse                !!  
    REAL                                     :: xlam               !!
    REAL                                     :: dlamm              !!
    REAL                                     :: anm                !!
    REAL                                     :: ranm               !!
    REAL                                     :: ranv               !!
    REAL                                     :: anv                !!
    REAL                                     :: tls                !!
    REAL                                     :: rlam               !!
    REAL                                     :: sd                 !!
    REAL                                     :: cd                 !!
    REAL                                     :: deltar             !!
    REAL                                     :: delta              !!
    REAL                                     :: Dis_ST             !!  
    REAL                                     :: ddt                !!
    REAL                                     :: orbit              !!
    REAL                                     :: angle              !!
    REAL                                     :: xdecl              !!
    REAL                                     :: xlat               !!
    REAL, DIMENSION(npoi)                    :: trans              !!
    REAL, DIMENSION(npoi)                    :: fdiffuse           !!
    REAL, DIMENSION(npoi)                    :: coszen             !! cosine of solar zenith angle
    REAL                                     :: sw                 !!
    REAL                                     :: frac               !!
    INTEGER                                  :: i,ib               !!
    LOGICAL, SAVE                            :: lstep_init = .TRUE. !!
!$OMP THREADPRIVATE(lstep_init)

!_ ================================================================================================================================
    
    IF (lstep_init) THEN
       IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN  
          step = 1.0
       ELSE
          step = one_year/365.2425
       ENDIF
       !-
       ! Read Orbital Parameters
       !-

       !Config Key   = ECCENTRICITY
       !Config Desc  = Use prescribed values
       !Config If    = ALLOW_WEATHERGEN
       !Config Def   = 0.016724
       !Config Help  =
       !Config Units = [-]
       ecc = 0.016724
       CALL getin_p ('ECCENTRICITY',ecc)
       IF (printlev >= 1) WRITE(numout,*) 'ECCENTRICITY: ',ecc
       !
       !Config Key  = PERIHELIE
       !Config Desc = Use prescribed values
       !Config If   = ALLOW_WEATHERGEN
       !Config Def  = 102.04
       !Config Help  =
       !Config Units = [-]
       perh = 102.04
       CALL getin_p ('PERIHELIE',perh)
       IF (printlev >= 1) WRITE(numout,*) 'PERIHELIE: ',perh
       !
       !Config Key  = OBLIQUITY
       !Config Desc = Use prescribed values
       !Config If   = ALLOW_WEATHERGEN
       !Config Def  = 23.446
       !Config Help  =
       !Config Units = [Degrees]
       xob = 23.446
       CALL getin_p ('OBLIQUITY',xob)
       IF (printlev >= 1) WRITE(numout,*) 'OBLIQUITY: ',xob

       lstep_init = .FALSE.
    ENDIF

    !-
    ! calendar and orbital calculations
    !-

    !-
    ! calculate the earth's orbital angle (around the sun) in radians
    orbit = 2.0*pi*jday/365.2425
    !-
    ! calculate the solar hour angle in radians
    angle = 2.0*pi*(rtime-12.0)/24.0
    !-
    ! calculate the current solar declination angle
    ! ref: global physical climatology, hartmann, appendix a
    !
    ! xdecl = 0.006918 &
    !        -0.399912*cos(orbit) &
    !        +0.070257*sin(orbit) &
    !        -0.006758*cos(2.0*orbit) &
    !        +0.000907*sin(2.0*orbit) &
    !        -0.002697*cos(3.0*orbit) &
    !        +0.001480*sin(3.0*orbit)
    !
    ! calculate the effective solar constant,
    ! including effects of eccentricity
    ! ref: global physical climatology, hartmann, appendix a
    !
    ! sw = 1370.*( 1.000110 &
    !             +0.034221*cos(orbit) &
    !             +0.001280*sin(orbit) &
    !             +0.000719*cos(2.0*orbit) &
    !             +0.000077*sin(2.0*orbit))
    !
    ! correction Marie-France Loutre
    !
    !    orbital parameters
    !
    !    ecc = 0.016724
    !    perh = 102.04
    !    xob = 23.446
    !-
    xl = perh+180.0
    ! so : sinus of obliquity
    so = sin(xob*pir)
    !-
    xllp = xl*pir
    xee  = ecc*ecc
    xse  = sqrt(1.0d0-xee)
    ! xlam : true long. sun for mean long. = 0
    xlam = (ecc/2.0+ecc*xee/8.0d0)*(1.0+xse)*sin(xllp)-xee/4.0 &
         &      *(0.5+xse)*sin(2.0*xllp)+ecc*xee/8.0*(1.0/3.0+xse) &
         &      *sin(3.0*xllp)
    xlam  =2.0d0*xlam/pir
    ! dlamm : mean long. sun for ma-ja
    dlamm =xlam+(INT(jday)-79)*step
    anm   = dlamm-xl
    ranm  = anm*pir
    xee   = xee*ecc
    ! ranv : anomalie vraie    (radian)
    ranv = ranm+(2.0*ecc-xee/4.0)*sin(ranm)+5.0/4.0*ecc*ecc &
         &      *sin(2.0*ranm)+13.0/12.0*xee*sin(3.0*ranm)
    ! anv  : anomalie vraie    (degrees)
    anv  = ranv/pir
    ! tls  : longitude vraie   (degrees)
    tls  = anv+xl
    ! rlam : longitude vraie   (radian)
    rlam = tls*pir
    ! sd and cd: cosinus and sinus of solar declination angle (delta)
    ! sinus delta = sin (obl)*sin(lambda) with lambda = real longitude
    ! (Phd. thesis of Marie-France Loutre, ASTR-UCL, Belgium, 1993)
    sd    = so*sin(rlam)
    cd    = sqrt(1.0d0-sd*sd)
    ! delta : Solar Declination (degrees, angle sun at equator)
    deltar = atan(sd/cd)
    delta  = deltar/pir
    !-
    ! Eccentricity Effect
    !-
    Dis_ST  =(1.0-ecc*ecc)/(1.0+ecc*cos(ranv))
    ! ddt :    1 / normalized earth's sun distance
    ddt = 1.0/Dis_ST
    !-
    ! Insolation normal to the atmosphere (W/m2)
    !-
    sw  = ddt *ddt *1370.d0
    !-
    xdecl = deltar
    !-
    ! solar calculations
    !-
    do i=1,npoi
       !---
       !-- calculate the latitude in radians
       !---
       xlat = latitude(i)*pi/180.0
       !---
       !-- calculate the cosine of the solar zenith angle
       !---
       coszen(i) = MAX(0.0, (sin(xlat)*sin(xdecl) &
            &                      + cos(xlat)*cos(xdecl)*cos(angle)))
       !---
       !-- calculate the solar transmission through the atmosphere
       !-- using simple linear function of tranmission and cloud cover
       !---
       !-- note that the 'cloud cover' data is typically obtained from
       !-- sunshine hours -- not direct cloud observations
       !---
       !-- where, cloud cover = 1 - sunshine fraction
       !---
       !-- different authors present different values for the slope and
       !-- intercept terms of this equation
       !---
       !-- Friend, A: Parameterization of a global daily weather generator
       !-- for terrestrial ecosystem and biogeochemical modelling,
       !-- Ecological Modelling
       !---
       !-- Spitters et al., 1986: Separating the diffuse and direct component
       !-- of global radiation and its implications for modeling canopy
       !-- photosynthesis, Part I: Components of incoming radiation,
       !-- Agricultural and Forest Meteorology, 38, 217-229.
       !---
       !-- A. Friend       : trans = 0.251+0.509*(1.0-cloud(i))
       !-- Spitters et al. : trans = 0.200+0.560*(1.0-cloud(i))
       !---
       !-- we are using the values from A. Friend
       !---
       trans(i) = 0.251+0.509*(1.0-cloud(i))
       !---
       !-- calculate the fraction of indirect (diffuse) solar radiation
       !-- based upon the cloud cover
       !---
       !-- note that these relationships typically are measured for either
       !-- monthly or daily timescales, and may not be exactly appropriate
       !-- for hourly calculations -- however, in ibis, cloud cover is fixed
       !-- through the entire day so it may not make much difference
       !---
       !-- method i --
       !---
       !-- we use a simple empirical relationships
       !-- from Nikolov and Zeller (1992)
       !---
       !-- Nikolov, N. and K.F. Zeller, 1992:  A solar radiation algorithm
       !-- for ecosystem dynamics models, Ecological Modelling, 61, 149-168.
       !---
       fdiffuse(i) = 1.0045+trans(i)*( 0.0435+trans(i) &
            &                                 *(-3.5227+trans(i)*2.6313))
       !---
       IF (trans(i) > 0.75) fdiffuse(i) = 0.166
       !---
       !-- method ii --
       !---
       !-- another method was suggested by Spitters et al. (1986) based on
       !-- long-term data from the Netherlands
       !--
       !-- Spitters et al., 1986: Separating the diffuse and direct component
       !-- of global radiation and its implications for modeling canopy
       !-- photosynthesis, Part I: Components of incoming radiation,
       !-- Agricultural and Forest Meteorology, 38, 217-229.
       !---
       !--       if ((trans == 0.00).and.(trans < 0.07)) then
       !--         fdiffuse = 1.0
       !--       else if ((trans >= 0.07).and.(trans < 0.35)) then
       !--         fdiffuse = 1.0-2.3*(trans-0.07)**2
       !--       else if ((trans >= 0.35).and.(trans < 0.75)) then
       !--         fdiffuse = 1.33-1.46*trans
       !--       ELSE
       !--         fdiffuse = 0.23
       !--       endif
       !---
    ENDDO
    !-----------------------
    !-
    ! do for each waveband
    !-
    DO ib=1,nband
       !---
       !-- calculate the fraction in each waveband
       !---
       !-- GK010200
       IF (nband == 2) then
          frac = 0.46+0.08*REAL(ib-1)
       ELSE
          frac = 1./REAL(nband)
       ENDIF
       !---
       DO i=1,npoi
          !-----
          !---- calculate the direct and indirect solar radiation
          !-----
          solad(i,ib) = sw*coszen(i)*frac*trans(i)*(1.-fdiffuse(i))
          solai(i,ib) = sw*coszen(i)*frac*trans(i)*fdiffuse(i)
       ENDDO
    ENDDO


  END SUBROUTINE downward_solar_flux

 
END MODULE solar
