MODULE radiative
  IMPLICIT NONE
  SAVE
  
  PRIVATE

  PUBLIC :: mucorr, sw, lwtr

CONTAINS
  
  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.)

    !-----------------------------------------------------------------------
    
!    print*,'npts,pdeclin',npts,pdeclin
    z = pdeclin
    cz = cos (z)
    sz = sin (z)
!    print*,'cz,sz',cz,sz
    
    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
  
  PURE SUBROUTINE monGATHER(n,a,b,index)      
    INTEGER, INTENT(IN) ::  n,index(n)
    REAL, INTENT(IN) :: b(n)
    REAL, INTENT(OUT) :: a(n)
    INTEGER :: i

    DO i=1,n
       a(i)=b(index(i))
    END DO
  END SUBROUTINE monGATHER
  
  SUBROUTINE sw(ngrid,nlayer,ldiurn, &
       coefvis,albedo, &
       plevel,ps_rad,pmu,pfract,psolarf0, &
       fsrfvis,dtsw, &
       lwrite)
    USE phys_const
!=======================================================================
!
!   Rayonnement solaire en atmosphere non diffusante avec un 
!   coefficient d'absoprption gris.
!
!=======================================================================
!
!   declarations:
!   -------------
!
!
!   arguments:
!   ----------
!
      INTEGER ngrid,nlayer
      REAL albedo(ngrid),coefvis
      REAL pmu(ngrid),pfract(ngrid)
      REAL plevel(ngrid,nlayer+1),ps_rad
      REAL psolarf0
      REAL fsrfvis(ngrid),dtsw(ngrid,nlayer)
      LOGICAL lwrite,ldiurn
!
!   variables locales:
!   ------------------
!

      REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)
      REAL zplev(ngrid,nlayer+1)
      REAL zflux(ngrid),zdtsw(ngrid,nlayer)

      INTEGER ig,l,nlevel,index(ngrid),ncount,igout
      REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)
      REAL zfsrfref(ngrid)
      REAL z1(ngrid)
      REAL zu(ngrid,nlayer+1)
      REAL tau0

!-----------------------------------------------------------------------
!   1. initialisations:
!   -------------------

 
      nlevel=nlayer+1

!-----------------------------------------------------------------------
!   Definitions des tableaux locaux pour les points ensoleilles:
!   ------------------------------------------------------------

      IF (ldiurn) THEN
         ncount=0
         DO ig=1,ngrid
            index(ig)=0
         ENDDO
         DO ig=1,ngrid
            IF(pfract(ig).GT.1.e-6) THEN
               ncount=ncount+1
               index(ncount)=ig
            ENDIF
         ENDDO
         CALL monGATHER(ncount,zfract,pfract,index)
         CALL monGATHER(ncount,zmu,pmu,index)
         CALL monGATHER(ncount,zalb,albedo,index)
         DO l=1,nlevel
            CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index)
         ENDDO
      ELSE
         ncount=ngrid
         zfract(:)=pfract(:)
         zmu(:)=pmu(:)
         zalb(:)=albedo(:)
         zplev(:,:)=plevel(:,:)
      ENDIF

!-----------------------------------------------------------------------
!   calcul des profondeurs optiques integres depuis p=0:
!   ----------------------------------------------------

      tau0=-.5*log(coefvis)

! calcul de la partie homogene de l'opacite
      tau0=tau0/ps_rad
      DO l=1,nlayer+1
         DO ig=1,ncount
            zu(ig,l)=tau0*zplev(ig,l)
         ENDDO
      ENDDO

!-----------------------------------------------------------------------
!   2. calcul de la transmission depuis le sommet de l'atmosphere:
!   -----------------------------------------------------------

      DO l=1,nlevel
         DO ig=1,ncount
            ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig))
         ENDDO
      ENDDO

      IF (lwrite) THEN
         igout=ncount/2+1
         PRINT*
         PRINT*,'Diagnostique des transmission dans le spectre solaire'
         PRINT*,'zfract, zmu, zalb'
         PRINT*,zfract(igout), zmu(igout), zalb(igout)
         PRINT*,'Pression, quantite d abs, transmission'
         DO l=1,nlayer+1
            PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l)
         ENDDO
      ENDIF

!-----------------------------------------------------------------------
!   3. taux de chauffage, ray. solaire direct:
!   ------------------------------------------

      DO l=1,nlayer
         DO ig=1,ncount
            zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)*      &
                          (ztrdir(ig,l+1)-ztrdir(ig,l))/    &
                          (cpp*(zplev(ig,l)-zplev(ig,l+1)))
         ENDDO
      ENDDO
      IF (lwrite) THEN
         PRINT*
         PRINT*,'Diagnostique des taux de chauffage solaires:'
         PRINT*,' 1 taux de chauffage lie au ray. solaire  direct'
         DO l=1,nlayer
            PRINT*,zdtsw(igout,l)
         ENDDO
      ENDIF


!-----------------------------------------------------------------------
!   4. calcul du flux solaire arrivant sur le sol:
!   ----------------------------------------------

      DO ig=1,ncount
         z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)
         zflux(ig)=(1.-zalb(ig))*z1(ig)
         zfsrfref(ig)=    zalb(ig)*z1(ig)
      ENDDO
      IF (lwrite) THEN
         PRINT*
         PRINT*,'Diagnostique des taux de chauffage solaires:'
         PRINT*,' 2 flux solaire net incident sur le sol'
         PRINT*,zflux(igout)
      ENDIF


!-----------------------------------------------------------------------
!   5.calcul des traansmissions depuis le sol, cas diffus:
!   ------------------------------------------------------

      DO l=1,nlevel
         DO ig=1,ncount
            ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66)
         ENDDO
      ENDDO

      IF (lwrite) THEN
         PRINT*
         PRINT*,'Diagnostique des taux de chauffage solaires'
         PRINT*,' 3 transmission avec les sol'
         PRINT*,'niveau     transmission'
         DO l=1,nlevel
            PRINT*,l,ztrref(igout,l)
         ENDDO
      ENDIF

!-----------------------------------------------------------------------
!   6.ajout a l'echauffement de la contribution du ray. sol. reflechit: 
!   -------------------------------------------------------------------

      DO l=1,nlayer
         DO ig=1,ncount
            zdtsw(ig,l)=zdtsw(ig,l)+ &
                 g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/ &
                 (cpp*(zplev(ig,l+1)-zplev(ig,l)))
         ENDDO
      ENDDO

!-----------------------------------------------------------------------
!   10. sorites eventuelles:
!   ------------------------

      IF (lwrite) THEN
         PRINT*
         PRINT*,'Diagnostique des taux de chauffage solaires:'
         PRINT*,' 3 taux de chauffage total'
         DO l=1,nlayer
            PRINT*,zdtsw(igout,l)
         ENDDO
      ENDIF

      IF (ldiurn) THEN
         fsrfvis(:)=0.
         CALL monscatter(ncount,fsrfvis,index,zflux)
         dtsw(:,:)=0.
         DO l=1,nlayer
            CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l))
         ENDDO
      ELSE
         print*,'NOT DIURNE'
         fsrfvis(:)=zflux(:)
         dtsw(:,:)=zdtsw(:,:)
      ENDIF
!        call dump2d(iim,jjm-1,zflux(2),'ZFLUX      ')
!        call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS     ')
!        call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir     ')
!        call dump2d(iim,jjm-1,pmu(2),'pmu        ')
!        call dump2d(iim,jjm-1,pfract(2),'pfract     ')
!        call dump2d(iim,jjm-1,albedo(2),'albedo     ')
!        call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir     ')


    END SUBROUTINE sw

    PURE SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
      INTEGER, INTENT(IN) :: ngrid
      REAL,    INTENT(IN) :: coef
      LOGICAL, INTENT(IN) :: lstrong
      REAL,    INTENT(IN) :: dup(ngrid)
      REAL,    INTENT(OUT) :: transm(ngrid)
      INTEGER ig
      IF(lstrong) THEN
         DO ig=1,ngrid
            transm(ig)=exp(-coef*sqrt(dup(ig)))
         ENDDO
      ELSE
         DO ig=1,ngrid
            transm(ig)=exp(-coef*dup(ig))
         ENDDO
      ENDIF
      
    END SUBROUTINE lwtr

END MODULE radiative