SUBROUTINE sw(ngrid,nlayer,ldiurn,
     $              coefvis,albedo,
     $              plevel,ps_rad,pmu,pfract,psolarf0,
     $              fsrfvis,dtsw,
     $              lwrite)
      IMPLICIT NONE
c=======================================================================
c
c   Rayonnement solaire en atmosphere non diffusante avec un 
c   coefficient d'absoprption gris.
c
c=======================================================================
c
c   declarations:
c   -------------
c
#include "comcstfi.h"
c
c   arguments:
c   ----------
c
      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
c
c   variables locales:
c   ------------------
c

      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

      LOGICAL firstcall
      SAVE firstcall
      DATA firstcall/.true./
!$OMP THREADPRIVATE(firstcall)

c-----------------------------------------------------------------------
c   1. initialisations:
c   -------------------

 
      nlevel=nlayer+1

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

      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

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

      tau0=-.5*log(coefvis)

c 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

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

      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

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

      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


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

      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


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

      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

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

      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

c-----------------------------------------------------------------------
c   10. sorites eventuelles:
c   ------------------------

      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
         CALL zerophys(ngrid,fsrfvis)
         CALL monscatter(ncount,fsrfvis,index,zflux)
         CALL zerophys(ngrid*nlayer,dtsw)
         DO l=1,nlayer
            CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l))
         ENDDO
      ELSE
         print*,'NOT DIURNE'
         fsrfvis(:)=zflux(:)
         dtsw(:,:)=zdtsw(:,:)
      ENDIF
c        call dump2d(iim,jjm-1,zflux(2),'ZFLUX      ')
c        call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS     ')
c        call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir     ')
c        call dump2d(iim,jjm-1,pmu(2),'pmu        ')
c        call dump2d(iim,jjm-1,pfract(2),'pfract     ')
c        call dump2d(iim,jjm-1,albedo(2),'albedo     ')
c        call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir     ')


      RETURN
      END