source: dynamico_lmdz/simple_physics/phyparam/physics/solar.F90 @ 4245

MODULE solar

#include "use_logging.h"

  IMPLICIT NONE
  PRIVATE

  REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local

  PUBLIC :: solang, zenang, mucorr

CONTAINS

  PURE SUBROUTINE solang ( kgrid,psilon,pcolon,psilat,pcolat, &
       ptim1,ptim2,ptim3,pmu0,pfract )

    !-----------------------------------------------------------------------
    !          CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID
    !
    !     ==== INPUTS  ===
    !
    ! PSILON(KGRID)   : SINUS OF THE LONGITUDE
    ! PCOLON(KGRID)   : COSINUS OF THE LONGITUDE
    ! PSILAT(KGRID)   : SINUS OF THE LATITUDE
    ! PCOLAT(KGRID)   : COSINUS OF THE LATITUDE
    ! PTIM1           : SIN(DECLI)
    ! PTIM2           : COS(DECLI)*COS(TIME)
    ! PTIM3           : SIN(DECLI)*SIN(TIME)
    !
    !     ==== OUTPUTS ===
    !
    ! PMU0 (KGRID)    : SOLAR ANGLE
    ! PFRACT(KGRID)   : DAY FRACTION OF THE TIME INTERVAL
    !
    !     REFERENCE.
    !     ----------
    !
    !         RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE
    !         PALTRIDGE AND PLATT
    !
    !     AUTHOR.
    !     -------
    !        FREDERIC HOURDIN
    !
    !     MODIFICATIONS.
    !     --------------
    !        ORIGINAL :90-01-14
    !                  92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher)
    !-----------------------------------------------------------------------

    INTEGER, INTENT(IN)  :: kgrid
    REAL,    INTENT(IN)  :: ptim1,ptim2,ptim3
    REAL,    INTENT(IN)  :: psilon(kgrid), pcolon(kgrid)
    REAL,    INTENT(IN)  :: psilat(kgrid), pcolat(kgrid)
    REAL,    INTENT(OUT) :: pmu0(kgrid), pfract(kgrid)

    INTEGER jl
    REAL ztim1,ztim2,ztim3
    !------------------------------------------------------------------------
    !------------------------------------------------------------------------
    !------------------------------------------------------------------------
    !
    !------------------------------------------------------------------------
    !
    !*     1.     INITIALISATION
    !             --------------
    !
    !
    DO jl=1,kgrid
       pmu0(jl)=0.
       pfract(jl)=0.
    ENDDO
    !
    !*     1.1     COMPUTATION OF THE SOLAR ANGLE
    !              ------------------------------
    !
    DO jl=1,kgrid
       ztim1=psilat(jl)*ptim1
       ztim2=pcolat(jl)*ptim2
       ztim3=pcolat(jl)*ptim3
       pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)
    ENDDO
    !
    !*     1.2      DISTINCTION BETWEEN DAY AND NIGHT
    !               ---------------------------------
    !
    DO jl=1,kgrid
       IF (pmu0(jl).gt.0.) THEN
          pfract(jl)=1.
       ELSE
          pmu0(jl)=0.
          pfract(jl)=0.
       ENDIF
    ENDDO

  END SUBROUTINE solang

  SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long,              &
       &                  pmu0,frac)
    USE astronomy

    !=============================================================
    ! Auteur : O. Boucher (LMD/CNRS)
    !          d apres les routines zenith et angle de Z.X. Li
    ! Objet  : calculer les valeurs moyennes du cos de l angle zenithal
    !          et l ensoleillement moyen entre gmtime1 et gmtime2
    !          connaissant la declinaison, la latitude et la longitude.
    ! Rque   : Different de la routine angle en ce sens que zenang
    !          fournit des moyennes de pmu0 et non des valeurs
    !          instantanees, du coup frac prend toutes les valeurs
    !          entre 0 et 1.
    ! Date   : premiere version le 13 decembre 1994
    !          revu pour  GCM  le 30 septembre 1996
    !===============================================================
    ! longi----INPUT : la longitude vraie de la terre dans son plan
    !                  solaire a partir de l equinoxe de printemps (degre)
    ! gmtime---INPUT : temps universel en fraction de jour
    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
    ! lat------INPUT : latitude en degres
    ! long-----INPUT : longitude en degres
    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
    !================================================================
    integer klon
    !================================================================
    real longi, gmtime, pdtrad
    real lat(klon), long(klon), pmu0(klon), frac(klon)
    !================================================================
    integer i
    real gmtime1, gmtime2
    real incl
    real omega1, omega2, omega
    ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.
    ! omega : heure en radian du coucher de soleil
    ! -omega est donc l heure en radian de lever du soleil
    real omegadeb, omegafin
    real zfrac1, zfrac2, z1_mu, z2_mu
    ! declinaison en radian
    real lat_sun
    ! longitude solaire en radian
    real lon_sun
    ! latitude du pt de grille en radian
    real latr
    !================================================================
    !
    !     incl=R_incl * pi_local / 180.
    WRITELOG(*,*) 'Obliquite =' ,obliquit
    LOG_INFO('solar')

    incl=obliquit * pi_local / 180.
    !
    !     lon_sun = longi * pi_local / 180.0
    lon_sun = longi
    lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )
    !
    gmtime1=gmtime*86400.
    gmtime2=gmtime*86400.+pdtrad
    !
    DO i = 1, klon
       !
       !     latr = lat(i) * pi_local / 180.
       latr = lat(i)
       !
       !--pose probleme quand lat=+/-90 degres
       !
       !      omega = -TAN(latr)*TAN(lat_sun)
       !      omega = ACOS(omega)
       !      IF (latr.GE.(pi_local/2.+lat_sun)
       !     .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN
       !         omega = 0.0       ! nuit polaire
       !      ENDIF
       !      IF (latr.GE.(pi_local/2.-lat_sun)
       !     .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
       !         omega = pi_local  ! journee polaire
       !      ENDIF
       !
       !--remplace par cela (le cas par defaut est different)
       !
       !--nuit polaire
       omega=0.0
       IF (latr.GE.(pi_local/2.-lat_sun)                                 &
            &          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
          ! journee polaire
          omega = pi_local
       ENDIF
       IF (latr.LT.(pi_local/2.+lat_sun).AND.                            &
            &    latr.GT.(-pi_local/2.+lat_sun).AND.                           &
            &    latr.LT.(pi_local/2.-lat_sun).AND.                            &
            &    latr.GT.(-pi_local/2.-lat_sun)) THEN
          omega = -TAN(latr)*TAN(lat_sun)
          omega = ACOS(omega)
       ENDIF
       !
       omega1 = gmtime1 + long(i)*86400.0/360.0
       omega1 = omega1 / 86400.0*deux_pi_local
       omega1 = MOD (omega1+deux_pi_local, deux_pi_local)
       omega1 = omega1 - pi_local
       !
       omega2 = gmtime2 + long(i)*86400.0/360.0
       omega2 = omega2 / 86400.0*deux_pi_local
       omega2 = MOD (omega2+deux_pi_local, deux_pi_local)
       omega2 = omega2 - pi_local
       !
       !--on est dans la meme journee locale
       IF (omega1.LE.omega2) THEN
          !
          IF (omega2.LE.-omega .OR. omega1.GE.omega .OR. omega.LT.1e-5)  &
               THEN
             !--nuit
             frac(i)=0.0
             pmu0(i)=0.0
             !--jour+nuit/jou
          ELSE
             omegadeb=MAX(-omega,omega1)
             omegafin=MIN(omega,omega2)
             frac(i)=(omegafin-omegadeb)/(omega2-omega1)
             pmu0(i)=SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*    &
                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
          ENDIF
          !
          !---omega1 GT omega2 -- a cheval sur deux journees
       ELSE
          !
          !-------------------entre omega1 et pi
          !--nuit
          IF (omega1.GE.omega) THEN
             zfrac1=0.0
             z1_mu =0.0
             !--jour+nuit
          ELSE
             omegadeb=MAX(-omega,omega1)
             omegafin=omega
             zfrac1=omegafin-omegadeb
             z1_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*     &
                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
          ENDIF
          !---------------------entre -pi et omega2
          !--nuit
          IF (omega2.LE.-omega) THEN
             zfrac2=0.0
             z2_mu =0.0
             !--jour+nuit
          ELSE
             omegadeb=-omega
             omegafin=MIN(omega,omega2)
             zfrac2=omegafin-omegadeb
             z2_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*     &
                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
             !
          ENDIF
          !-----------------------moyenne
          frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1)
          pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
          !
          !---comparaison omega1 et omega2
       ENDIF
       !
    ENDDO
    !
  END SUBROUTINE zenang

  PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)

    !=======================================================================
    !
    !   Calcul of equivalent solar angle and and fraction of day whithout
    !   diurnal cycle.
    !
    !   parmeters :
    !   -----------
    !
    !      Input :
    !      -------
    !         npts             number of points
    !         pdeclin          solar declinaison
    !         plat(npts)        latitude
    !         phaut            hauteur typique de l atmosphere
    !         prad             rayon planetaire
    !
    !      Output :
    !      --------
    !         pmu(npts)          equivalent cosinus of the solar angle
    !         pfract(npts)       fractionnal day
    !
    !=======================================================================

    !-----------------------------------------------------------------------
    !
    !    0. Declarations :
    !    -----------------

    !     Arguments :
    !     -----------
    INTEGER, INTENT(IN) :: npts
    REAL, INTENT(IN)    :: phaut, prad, pdeclin, plat(npts)
    REAL, INTENT(OUT)   :: pmu(npts), pfract(npts)
    !
    !     Local variables :
    !     -----------------
    INTEGER j
    REAL z,cz,sz,tz,phi,cphi,sphi,tphi
    REAL ap,a,t,b
    REAL alph

    REAL, PARAMETER :: pi=2.*asin(1.)

    !-----------------------------------------------------------------------

    z = pdeclin
    cz = cos (z)
    sz = sin (z)

    DO j = 1, npts

       phi = plat(j)
       cphi = cos(phi)
       if (cphi.le.1.e-9) cphi=1.e-9
       sphi = sin(phi)
       tphi = sphi / cphi
       b = cphi * cz
       t = -tphi * sz / cz
       a = 1.0 - t*t
       ap = a

       IF(t.eq.0.) then
          t=0.5*pi
       ELSE
          IF (a.lt.0.) a = 0.
          t = sqrt(a) / t
          IF (t.lt.0.) then
             t = -atan (-t) + pi
          ELSE
             t = atan(t)
          ENDIF
       ENDIF

       pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi
       pfract(j) = t / pi
       IF (ap .lt.0.) then
          pmu(j) = sphi * sz
          pfract(j) = 1.0
       ENDIF
       IF (pmu(j).le.0.0) pmu(j) = 0.0
       pmu(j) = pmu(j) / pfract(j)
       IF (pmu(j).eq.0.) pfract(j) = 0.

    END DO

    !-----------------------------------------------------------------------
    !   correction de rotondite:
    !   ------------------------

    alph=phaut/prad
    DO j=1,npts
       ! !!!!!!
       pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35.
       !    $          (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j))
    END DO

  END SUBROUTINE mucorr

END MODULE solar