MODULE radiative_sw

#include "use_logging.h"

  IMPLICIT NONE
  SAVE
  
  PRIVATE

  PUBLIC :: sw

CONTAINS
  
  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

  PURE subroutine monscatter(n,a,index,b)
    INTEGER, INTENT(IN) :: n,index(n)
    REAL, INTENT(IN) :: b(n)
    REAL, INTENT(OUT) :: a(n)
    INTEGER :: i
    DO i=1,n
       a(index(i))=b(i)
    END DO
  end subroutine monscatter
  
  SUBROUTINE sw(ngrid,nlayer,ldiurn, &
       coefvis,albedo, &
       plevel,ps_rad,pmu,pfract,psolarf0, &
       fsrfvis,dtsw, &
       lwrite)
    USE phys_const, ONLY : cpp, g
    !=======================================================================
    !
    !   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
       WRITELOG(*,*) 
       WRITELOG(*,*) 'Diagnostique des transmission dans le spectre solaire'
       WRITELOG(*,*) 'zfract, zmu, zalb'
       WRITELOG(*,*) zfract(igout), zmu(igout), zalb(igout)
       WRITELOG(*,*) 'Pression, quantite d abs, transmission'
       DO l=1,nlayer+1
          WRITELOG(*,*) 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
       WRITELOG(*,*) 
       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
       WRITELOG(*,*) ' 1 taux de chauffage lie au ray. solaire  direct'
       DO l=1,nlayer
          WRITELOG(*,*) 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
       WRITELOG(*,*) 
       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
       WRITELOG(*,*) ' 2 flux solaire net incident sur le sol'
       WRITELOG(*,*) 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
       WRITELOG(*,*) 
       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires'
       WRITELOG(*,*) ' 3 transmission avec les sol'
       WRITELOG(*,*) 'niveau     transmission'
       DO l=1,nlevel
          WRITELOG(*,*) 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
       WRITELOG(*,*) 
       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
       WRITELOG(*,*) ' 3 taux de chauffage total'
       DO l=1,nlayer
          WRITELOG(*,*) 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
       WRITELOG(*,*) '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     ')
    LOG_INFO('rad_sw')
    
  END SUBROUTINE sw
  
END MODULE radiative_sw