MODULE RADIATION
	 USE ICOSA
	 USE dimphys_mod
!	 USE PHY
	contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	 SUBROUTINE zerophys(n,x)
      IMPLICIT NONE
      INTEGER n
      REAL x(n)
      INTEGER i

      DO i=1,n
         x(i)=0.
      ENDDO
      RETURN
      END subroutine zerophys
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	 SUBROUTINE solarlong(pday,psollong)
      IMPLICIT NONE
c=======================================================================
c
c   Objet:
c   ------
c
c      Calcul de la distance soleil-planete et de la declinaison
c   en fonction du jour de l'annee.
c
c
c   Methode:
c   --------
c
c      Calcul complet de l'elipse
c
c   Interface:
c   ----------
c
c      Uncommon comprenant les parametres orbitaux.
c
c   Arguments:
c   ----------
c
c   Input:
c   ------
c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
c   lwrite        clef logique pour sorties de controle
c
c   Output:
c   -------
c   pdist_sol     distance entre le soleil et la planete
c                 ( en unite astronomique pour utiliser la constante 
c                  solaire terrestre 1370 Wm-2 )
c   pdecli        declinaison ( en radians )
c
c=======================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

c arguments:
c ----------

      REAL pday,pdist_sol,pdecli,psollong
      LOGICAL lwrite

c Local:
c ------

      REAL zanom,xref,zx0,zdx,zteta,zz
      INTEGER iter


c-----------------------------------------------------------------------
c calcul de l'angle polaire et de la distance au soleil :
c -------------------------------------------------------

c  calcul de l'zanomalie moyenne

      zz=(pday-peri_day)/year_day
      zanom=2.*pi*(zz-nint(zz))
      xref=abs(zanom)

c  resolution de lequation horaire  zx0 - e * sin (zx0) = xref
c  methode de Newton

      zx0=xref+e_elips*sin(xref)
      DO 110 iter=1,10
         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
         if(abs(zdx).le.(1.e-7)) goto 120
         zx0=zx0+zdx
110   continue
120   continue
      zx0=zx0+zdx
      if(zanom.lt.0.) zx0=-zx0

c zteta est la longitude solaire

      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))

      psollong=zteta-timeperi

      IF(psollong.LT.0.) psollong=psollong+2.*pi
      IF(psollong.GT.2.*pi) psollong=psollong-2.*pi

      RETURN
      END SUBROUTINE solarlong

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        SUBROUTINE orbite(pls,pdist_sol,pdecli)
      IMPLICIT NONE
c====================================================================
c
c   Objet:
c   ------
c
c   Distance from sun and declimation as a function of the solar
c   longitude Ls
c
c   Interface:
c   ----------
c   Arguments:
c   ----------
c
c   Input:
c   ------
c   pls          Ls
c
c   Output:
c   -------
c   pdist_sol     Distance Sun-Planet in UA
c   pdecli        declinaison ( en radians )
c
c====================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------
c arguments:
c ----------

      REAL pday,pdist_sol,pdecli,pls

c-----------------------------------------------------------------------

c Distance Sun-Planet

      pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))

c Solar declination

      pdecli= asin (sin(pls)*sin(obliquit*pi/180.))

c-----------------------------------------------------------------------
c   sorties eventuelles:
c   ---------------------

      END SUBROUTINE orbite
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        SUBROUTINE iniorbit
     $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
      IMPLICIT NONE
c=====================================================================
c
c   Auteur:
c   -------
c     Frederic Hourdin      22 Fevrier 1991
c
c   Objet:
c   ------
c    Initialisation du sous programme orbite qui calcule
c    a une date donnee de l'annee de duree year_day commencant
c    a l'equinoxe de printemps et dont le perihelie se situe
c    a la date peri_day, la distance au soleil et la declinaison.
c
c   Interface:
c   ----------
c   - Doit etre appele avant d'utiliser orbite.
c   - initialise le common comorbit
c
c   Arguments:
c   ----------
c
c   Input:
c   ------
c   aphelie       \   aphelie et perihelie de l'orbite
c   periheli      /   en millions de kilometres.
c
c=====================================================================
c   Declarations:
c   -------------
c   Arguments:
c   ----------

      REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq

c   Local:
c   ------

      REAL zxref,zanom,zz,zx0,zdx
      INTEGER iter

c'-----------------------------------------------------------------------

      aphelie =paphelie
      periheli=pperiheli
      year_day=pyear_day
      obliquit=pobliq
      peri_day=pperi_day

      PRINT*,'Perihelie en Mkm  ',periheli
      PRINT*,'Aphelise  en Mkm  ',aphelie 
      PRINT*,'obliquite en degres  :',obliquit
      unitastr=149.597927
      e_elips=(aphelie-periheli)/(periheli+aphelie)
      p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr

      print*,'e_elips',e_elips
      print*,'p_elips',p_elips

c-----------------------------------------------------------------------
c calcul de l'angle polaire et de la distance au soleil :
c -------------------------------------------------------

c  calcul de l'zanomalie moyenne

      zz=(year_day-pperi_day)/year_day
      zanom=2.*pi*(zz-nint(zz))
      zxref=abs(zanom)
      PRINT*,'zanom  ',zanom

c  resolution de lequation horaire  zx0 - e * sin (zx0) = zxref
c  methode de Newton

      zx0=zxref+e_elips*sin(zxref)
      DO 110 iter=1,100
         zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
         if(abs(zdx).le.(1.e-12)) goto 120
         zx0=zx0+zdx
110   continue
120   continue
      zx0=zx0+zdx
      if(zanom.lt.0.) zx0=-zx0
      PRINT*,'zx0   ',zx0

c zteta est la longitude solaire

      timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
      PRINT*,'longitude solaire du perihelie timeperi = ',timeperi

      RETURN
      END SUBROUTINE iniorbit
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

c-----------------------------------------------------------------------
c
c    0. Declarations :
c    -----------------

c     Arguments :
c     -----------
      INTEGER npts
      REAL plat(npts),pmu(npts), pfract(npts)
      REAL phaut,prad,pdeclin
c
c     Local variables :
c     -----------------
      INTEGER j,i,ij,ig
      REAL pi
      REAL z,cz,sz,tz,phi,cphi,sphi,tphi
      REAL ap,a,t,b
      REAL alph

c-----------------------------------------------------------------------

      !print*,'npts,pdeclin'
      !print*,npts,pdeclin
      pi = 4. * atan(1.0)
      !print*,'PI=',pi
      pi=2.*asin(1.)
      z = pdeclin
      cz = cos (z)
      sz = sin (z)
      !print*,'cz,sz',cz,sz

!	DO j=jj_begin-offset,jj_end+offset
!         DO i=ii_begin-offset,ii_end+offset
!             ig=(j-1)*iim+i

	  DO ig=1,npts
             phi = plat(ig)
             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(ig) = (sphi*sz*t) / pi + b*sin(t)/pi
            pfract(ig) = t / pi
            IF (ap .lt.0.) then
             pmu(ig) = sphi * sz
             pfract(ig) = 1.0
            ENDIF
            IF (pmu(ig).le.0.0) pmu(ig) = 0.0
             pmu(ig) = pmu(ig) / pfract(ig)
            IF (pmu(ig).eq.0.) pfract(ig) = 0.
	  ENDDO
!	ENDDO
c-----------------------------------------------------------------------
c   correction de rotondite:
c   ------------------------

	      alph=phaut/prad
!	DO j=jj_begin-offset,jj_end+offset
!         DO i=ii_begin-offset,ii_end+offset
!		           ig=(j-1)*iim+i
	   DO ig = 1,npts
            pmu(ig)=sqrt(1224.*pmu(ig)*pmu(ig)+1.)/35.

c    $    (sqrt(alph*alph*pmu(ig)*pmu(ig)+2.*alph+1.)-alph*pmu(ig))
	 ENDDO
!	ENDDO

      RETURN
      END SUBROUTINE mucorr
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	      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   arguments:
c   ----------
c
      INTEGER ngrid,nlayer,i,j,ij
      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,ncount,igout
	INTEGER,DIMENSION(ngrid)::index
      REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)
      REAL zfsrfref(ngrid)
      REAL z1(ngrid)
      REAL zu(ngrid,nlayer+1)
      REAL tau0

      EXTERNAL SSUM
      EXTERNAL ismax
      REAL ismax

      LOGICAL firstcall
      SAVE firstcall
      DATA firstcall/.true./

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


 
      IF (firstcall) THEN
 
      IF (ngrid.NE.ngrid) THEN
         PRINT*,'STOP in inifis'
         PRINT*,'Probleme de dimenesions :'
         PRINT*,'ngrid     = ',ngrid
         PRINT*,'ngrid   = ',ngrid
         STOP
      ENDIF
 
 
      IF (nlayer.NE.llm) THEN
         PRINT*,'STOP in inifis'
         PRINT*,'Probleme de dimenesions :'
         PRINT*,'nlayer     = ',nlayer
         PRINT*,'llm   = ',llm
         STOP
      ENDIF
 
      ENDIF

      nlevel=nlayer+1
c-----------------------------------------------------------------------
c   Definitions des tableaux locaux pour les points ensoleilles:
c   ------------------------------------------------------------

      IF (ldiurn) THEN
         ncount=0
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
            ig=(j-1)*iim+i
            index(ig)=0
           ENDDO
	 ENDDO
         DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
             ig=(j-1)*iim+i
           IF(pfract(ig).GT.1.e-6) THEN
             ncount=ncount+1
             index(ncount)=ig
           ENDIF
          ENDDO
         ENDDO
!	SARVESH CHANGED NCOUNT TO NGRID TO BE CONSISTENT ???

         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
c'
      tau0=tau0/ps_rad


      DO l=1,nlayer+1
        DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
            ig=(j-1)*iim+i
            zu(ig,l)=tau0*zplev(ig,l)
           ENDDO
         ENDDO
        ENDDO

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

       DO l=1,nlevel
         DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
            ig=(j-1)*iim+i
            ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig))
           ENDDO
          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 j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
            ig=(j-1)*iim+i
            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
       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 j=jj_begin-offset,jj_end+offset
        DO i=ii_begin-offset,ii_end+offset
          ig=(j-1)*iim+i
          z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)
          zflux(ig)=(1.-zalb(ig))*z1(ig)
          zfsrfref(ig)=    zalb(ig)*z1(ig)
        ENDDO
       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 j=jj_begin-offset,jj_end+offset
        DO i=ii_begin-offset,ii_end+offset
          ig=(j-1)*iim+i
            ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66)
         ENDDO
        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 j=jj_begin-offset,jj_end+offset
        DO i=ii_begin-offset,ii_end+offset
          ig=(j-1)*iim+i
            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
       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

      RETURN
      END  SUBROUTINE sw
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,
     $             pp,ps_rad,ptsurf,pt,
     $             pfluxir,pdtlw,
     $             lwrite)

      IMPLICIT NONE
c=======================================================================
c
c   calcul de l'evolution de la temperature sous l'effet du rayonnement
c   infra-rouge.
c   Pour simplifier, les transmissions sont precalculees et ne
c   dependent que de l'altitude.
c
c   arguments:
c   ----------
c'
c   entree:
c   -------
c      ngrid             nombres de points de la grille horizontale
c      nlayer              nombre de couches
c      ptsurf(ngrid)     temperature de la surface
c      pt(ngrid,nlayer)    temperature des couches
c      pp(ngrid,nlayer+1)  pression entre les couches
c      lwrite            variable logique pour sorties
c
c   sortie:
c   -------
c      pdtlw(ngrid,nlayer) taux de refroidissement
c      pfluxir(ngrid)    flux infrarouge sur le sol
c
c=======================================================================

!c   declarations:
!c   -------------
!c   arguments:
!c'   ----------

      INTEGER ngrid,nlayer
      REAL coefir,emissiv(ngrid),ps_rad
      REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
      REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)
      LOGICAL lwrite

c   variables locales:
c   ------------------

      INTEGER nlevel,ilev,ig,i,il
      REAL zplanck(ngridmx,nlayermx+1),zcoef
      REAL zfluxup(ngridmx,nlayermx+1),zfluxdn(ngridmx,nlayermx+1)
      REAL zflux(ngridmx,nlayermx+1)
      REAL zlwtr1(ngridmx),zlwtr2(ngridmx)
      REAL zup(ngridmx,nlayermx+1),zdup(ngridmx)
      REAL stephan

      LOGICAL lstrong
      SAVE lstrong,stephan
      DATA lstrong/.true./
c-----------------------------------------------------------------------
c   initialisations:
c   ----------------

      nlevel=nlayer+1
      stephan=5.67e-08


	pfluxir=0.0
	pdtlw=0.0
	!print*,"ngr,nlay,coefi",ngrid,nlayer,coefir
c-----------------------------------------------------------------------
c   2. calcul des quantites d'absorbants:
c'   -------------------------------------

c   absorption forte
      IF(lstrong) THEN
         DO ilev=1,nlevel
            DO ig=1,ngrid
               zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
            ENDDO
         ENDDO
         IF(lwrite) THEN
            DO ilev=1,nlayer
             PRINT*,' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
            ENDDO
         ENDIF
         zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))

c   absorption faible
      ELSE
         DO ilev=1,nlevel
            DO ig=1,ngrid
               zup(ig,ilev)=pp(ig,ilev)
            ENDDO
         ENDDO
         zcoef=-log(coefir)/ps_rad
      ENDIF


c-----------------------------------------------------------------------
c   2. calcul de la fonction de corps noir:
c   ---------------------------------------

      DO ilev=1,nlayer
         DO ig=1,ngrid
            zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
            zplanck(ig,ilev)=stephan*
     $      zplanck(ig,ilev)*zplanck(ig,ilev)
	 ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   4. flux descendants:
c   --------------------

      DO ilev=1,nlayer
         DO ig=1,ngrid
            zfluxdn(ig,ilev)=0.
		 ENDDO
         DO ig=1,ngrid
            zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
         ENDDO
         CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)

         DO il=nlayer,ilev,-1
            zlwtr2(:)=zlwtr1(:)
            DO ig=1,ngrid
               zdup(ig)=zup(ig,ilev)-zup(ig,il)
            ENDDO
            CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
            DO ig=1,ngrid
               zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+
     $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
			ENDDO
		 ENDDO
      ENDDO

      DO ig=1,ngrid
         zfluxdn(ig,nlevel)=0.
         pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
      ENDDO

      DO ig=1,ngrid
         zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
         zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1)
     $   +(1.-emissiv(ig))*zfluxdn(ig,1)
      ENDDO

c-----------------------------------------------------------------------
c   3. flux montants:
c   ------------------

      DO ilev=1,nlayer
         DO ig=1,ngrid
            zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
         ENDDO
         CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
         DO ig=1,ngrid
            zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
		 ENDDO
         DO il=1,ilev
            zlwtr2(:)=zlwtr1(:)
            DO ig=1,ngrid
               zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
            ENDDO
            CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
            DO ig=1,ngrid
               zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+
     $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
			ENDDO
         ENDDO

      ENDDO

c-----------------------------------------------------------------------
c   5. calcul des flux nets:
c   ------------------------

      DO ilev=1,nlevel
         DO ig=1,ngrid
            zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
		 ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   6. Calcul des taux de refroidissement:
c   --------------------------------------
   
      DO ilev=1,nlayer
         DO ig=1,ngrid
            pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*
     $           g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
		 ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   10. sorties eventuelles:
c   ------------------------

      IF (lwrite) THEN

         PRINT*
         PRINT*,'Diagnostique rayonnement thermique'
         PRINT*
         PRINT*,'temperature     ',
     $   'flux montant    flux desc.     taux de refroid.'
         i=ngrid/2+1
         WRITE(6,9000) ptsurf(i)
         DO ilev=1,nlayer
            WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),
     $      zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
		 ENDDO
         WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)

      ENDIF

c-----------------------------------------------------------------------

      RETURN
9000  FORMAT(4e18.4)
      END SUBROUTINE lw 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      subroutine solang ( kgrid,psilon,pcolon,psilat,pcolat,
     &                    ptim1,ptim2,ptim3,pmu0,pfract )
      IMPLICIT NONE

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

C-----------------------------------------------------------------------
C
C*      0.1   ARGUMENTS
C             ---------
C
      INTEGER kgrid
      REAL ptim1,ptim2,ptim3
      REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)
      REAL psilat(kgrid), pcolat(kgrid)
C
      INTEGER jl,i,j
      REAL ztim1,ztim2,ztim3
C------------------------------------------------------------------------
C------------------------------------------------------------------------
C---------------------------------------------------------------------
C
C--------------------------------------------------------------------
C
C*     1.     INITIALISATION
C             --------------
C
!!!!!! 100  CONTINUE
C
      DO j=jj_begin-offset,jj_end+offset
       DO i=ii_begin-offset,ii_end+offset
         jl=(j-1)*iim+i
        pmu0(jl)=0.
        pfract(jl)=0.
       ENDDO
      ENDDO
!C
!C*     1.1     COMPUTATION OF THE SOLAR ANGLE
!C              ------------------------------
!C
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
          jl=(j-1)*iim+i
         ztim1=psilat(jl)*ptim1
         ztim2=pcolat(jl)*ptim2
         ztim3=pcolat(jl)*ptim3
         pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)
        ENDDO
      ENDDO
!C
!C*     1.2      DISTINCTION BETWEEN DAY AND NIGHT
!C               ---------------------------------
!C
       DO j=jj_begin-offset,jj_end+offset
        DO i=ii_begin-offset,ii_end+offset
          jl=(j-1)*iim+i      
          IF (pmu0(jl).gt.0.) THEN
            pfract(jl)=1.
	   ELSE
            pmu0(jl)=0.
	    pfract(jl)=0.
          ENDIF
         ENDDO
        ENDDO
      RETURN
      END SUBROUTINE solang
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      SUBROUTINE monGATHER(n,a,b,index)
c
      IMPLICIT NONE
C
      INTEGER n,ng,index(n),i,j,ij
      REAL a(n),b(n)
c
      DO 100 i=1,n
        a(i)=b(index(i))
100   CONTINUE

!       DO j=jj_begin-offset,jj_end+offset
!         DO i=ii_begin-offset,ii_end+offset
!          ij=(j-1)*iim+i
!          a(ij)=b(index(ij))
!c
      RETURN
      END  SUBROUTINE monGATHER
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

       subroutine monscatter(n,a,index,b)
C
       implicit none
       integer N,INDEX(n),I
       real A(n),B(n)
c
c
       DO 100 I=1,N
          A(INDEX(I))=B(I)
100    CONTINUE
c
       return
       end SUBROUTINE monscatter
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
      IMPLICIT NONE
      INTEGER ngrid
      REAL coef
      LOGICAL lstrong
      REAL dup(ngrid),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
      RETURN
      END subroutine lwtr 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       END MODULE RADIATION