oi_tsl.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 oi_tsl(KDAT,KSSSSS,PLAT,PLON,PMU0,PMU0M,KH)
!-----------------------------------------------------------------------
!
!     Computation of solar zenith angle 
!     ---------------------------------
!
!     INPUT PARAMETERS :
!
!     IDAT    : DATE in the following form => YYYYMMDD
!     NSSSSS  : TIME of the day in seconds
!     PLAT    : LATITUDE  in Degrees
!     PLON    : LONGITUDE in Degrees
!
!     OUTPUT PARAMETERS :
!
!     PMU0    : Cosine of solar zenith angle
!     PMU0M   : Cosine of solar zenith angle (mean value)
!     IH      : local time (hour)
!
!
!     J.F. Mahfouf (4/12/97) from IFS/ARPEGE routines
!                 
!
!     23/05/2009 : Fortran 90 recoding (IMPLICIT NONE + FUNCTIONS)
!
!-----------------------------------------------------------------------
!
! - Astronomical functions
!   you will find the description in the annex 1 of the documentation
!   RRS is the distance Sun-Earth
!   RDS is the declination of the Earth
!   RET is the equation of time
! 
USE modd_csts,  ONLY : xpi
USE modd_assim, ONLY : xrepsm, xrcdtr, nitrad
!
!
USE yomhook,  ONLY : lhook,   dr_hook
USE parkind1, ONLY : jprb
!
IMPLICIT NONE
!
INTEGER, INTENT(IN)  :: kdat
INTEGER, INTENT(IN)  :: ksssss 
REAL, INTENT(IN)     :: plat, plon
REAL, INTENT(OUT)    :: pmu0, pmu0m
INTEGER, INTENT(OUT) :: kh
!
REAL :: zgemu, zgelam, ztime, zteta, zrdecli, zreqtim, zrhgmt, zrsovr, &
        zrwsovr, zrcodec, zrsidec, zrcovsr, zrsivsr, zrtimtrm, ztetam,   &
        zrdeclim, zreqtimm, zrhgmtm, zrsovrm, zrwsovrm, zrcodecm,        &
        zrsidecm, zrcovsrm, zrsivsrm, zt  
INTEGER :: id, im, ia, insssss
REAL(KIND=JPRB) :: zhook_handle
!
!  Angle conversions 
!      
IF (lhook) CALL dr_hook('OI_TSL',0,zhook_handle)
!
zgemu  = sin(plat*xpi/180.)      ! sinus of latitude
zgelam = plon*xpi/180.          ! longitude
!
id = mod(kdat,100)
im = mod((kdat-id)/100,100)
ia = kdat/10000

ztime = rtime(ia,im,id,ksssss)
!
 CALL get_mu0(zgemu,zgelam,ztime,ksssss,pmu0)
!
!
!     Mean angle over the previous 6 hours
!     ------------------------------------
!  
ztime = rtime(ia,im,id,ksssss-nitrad)
!
IF ( (ksssss-nitrad).LT.0 ) THEN
  insssss = ksssss + 86400
ELSE
  insssss = ksssss
ENDIF  
!
 CALL get_mu0(zgemu,zgelam,ztime,insssss-nitrad,pmu0m)
!
! 
!  Local time in hours
!  Should be inside [1,24]
!
zt = (ksssss + plon*xrcdtr*3600.)/3600.
IF (zt < 0.0) THEN
  zt = zt + 24.
ELSEIF (zt > 24.0) THEN
  zt = zt - 24.
ENDIF
!
kh = int(zt)
IF ( kh==0 ) kh = 24 
!
IF (lhook) CALL dr_hook('OI_TSL',1,zhook_handle)
!
 CONTAINS
!
SUBROUTINE get_mu0(PGEMU,PGELAM,PTIME,KSSSSS,PMU)
!
USE modd_csts,  ONLY : xday
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: pgemu
REAL, INTENT(IN) :: pgelam
REAL, INTENT(IN) :: ptime
INTEGER, INTENT(IN) :: ksssss
REAL, INTENT(OUT) :: pmu
!
REAL :: zteta, zdecli, zeqtim, zhgmt, zwsovr
REAL :: zinter
!
zteta   = ptime/(xday*365.25)
!
zdecli = rds(zteta)             ! declinaison
!
zeqtim = ret(zteta)
zhgmt  = REAL( MOD(KSSSSS,NINT(XDAY)) )
zwsovr = (zeqtim + zhgmt) * 2. * xpi/xday      ! hour angle
!
zinter = sqrt(1.-pgemu**2)
pmu = max( sin(zdecli) * pgemu + cos(zdecli) * zinter * &
           ( sin(zwsovr)*sin(pgelam) - cos(zwsovr)*cos(pgelam) ) , 0. ) 
!
IF (pmu.GT.0.) pmu = sqrt(1224.*pmu*pmu +1.)/35. ! Magnification factor 
!
END SUBROUTINE get_mu0
!
FUNCTION rds(PT)
 
USE modd_assim, ONLY : xrepsm

IMPLICIT NONE

REAL, INTENT(IN) :: pt
REAL             :: rds, zls, zel
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('RDS',0,zhook_handle)
!
zel = 1.7535 + 6.283076 * pt
zls = 4.8952 + 6.283320 * pt &
     - 0.0075*sin(zel) - 0.0326*cos(zel) &
     - 0.0003*sin(2.*zel) + 0.0002*cos(2.*zel)  
!
rds   = asin(sin(xrepsm)*sin(zls))
!
IF (lhook) CALL dr_hook('RDS',1,zhook_handle)

END FUNCTION rds
! 
FUNCTION ret(PT)
 
IMPLICIT NONE

REAL, INTENT(IN) :: pt
REAL             :: ret, zem, zls
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('RET',0,zhook_handle)

zem  = 6.240075 + 6.283020 * pt
zls  = 4.8951   + 6.283076 * pt
ret  =  591.8*sin(2.*zls) - 12.7 * sin(4.*zls) &
      - 459.4*sin(zem) - 4.8 *sin(2.*zem) &
      + 39.5 * sin(zem) * cos(2.*zls) 
          
IF (lhook) CALL dr_hook('RET',1,zhook_handle)

END FUNCTION ret
!
FUNCTION rtime(KA,KM,KD,KS)

USE modd_csts, ONLY : xday

IMPLICIT NONE 

INTEGER, INTENT(IN) :: ka, km, kd, ks
REAL                :: rtime, zj
INTEGER             :: ia, im
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('RTIME',0,zhook_handle)

im  = km + 6*(1-isign(1,km-3))
ia  = ka - ( (1-isign(1,km-3))/2 )
!
zj = 1720994.5 + &
     float( 2 - ia/100 + (ia/100)/4 + &
           int(365.25*float(ia)) + int(30.601*float(im+1)) + kd )  
!
rtime  = (zj-2451545.)*xday + float(ks)
!
IF (lhook) CALL dr_hook('RTIME',1,zhook_handle)

END FUNCTION rtime

END SUBROUTINE oi_tsl