MODULE PHY
         USE dimphys_mod
         LOGICAL:: firstcall,lastcall
	contains      

	SUBROUTINE phyparam_lmd(it,ngrid,nlayer,nq,
     s            ptimestep,lati,long,rjourvrai,gmtime,
     s            pplev,pplay,pphi,pphis,
     s            pu,pv,pt,pq,
     s            pdu,pdv,pdt,pdq,pdpsrf)

	USE ICOSA
	USE dimphys_mod
	USE RADIATION
	USE SURFACE_PROCESS
c
      IMPLICIT NONE
c=======================================================================
c
c   subject:
c   --------
c
c   Organisation of the physical parametrisations of the LMD 
c   20 parameters GCM for planetary atmospheres.
c   It includes:
c   raditive transfer (long and shortwave) for CO2 and dust.
c   vertical turbulent mixing
c   convective adjsutment
c
c   author: Frederic Hourdin 15 / 10 /93
c   -------
c
c   arguments:
c   ----------
c
c   input:
c   ------
c
c    ngrid                 Size of the horizontal grid.
c                          All internal loops are performed on that grid.
c    nlayer                Number of vertical layers.
c    nq                    Number of advected fields
c    firstcall             True at the first call
c    lastcall              True at the last call
c    rjourvrai                  Number of days counted from the North. Spring
c                          equinoxe.
c    gmtime                 hour (s)
c    ptimestep             timestep (s)
c    pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa)
c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
c    pu(ngrid,nlayer)      u component of the wind (ms-1)
c    pv(ngrid,nlayer)      v component of the wind (ms-1)
c    pt(ngrid,nlayer)      Temperature (K)
c    pq(ngrid,nlayer,nq)   Advected fields
c    pudyn(ngrid,nlayer)    \ 
c    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
c    ptdyn(ngrid,nlayer)     / corresponding variables
c    pqdyn(ngrid,nlayer,nq) /
c    pw(ngrid,?)           vertical velocity
c
c   output:
c   -------
c
c    pdu(ngrid,nlayer)        \
c    pdv(ngrid,nlayer)         \  Temporal derivative of the corresponding
c    pdt(ngrid,nlayer)         /  variables due to physical processes.
c    pdq(ngrid,nlayer)      /
c    pdpsrf(ngrid)        /
c
c=======================================================================
c
!c-----------------------------------------------------------------------
!c
!c    0.  Declarations :
!c    ------------------
!c    Arguments :
!c    -----------

!c    inputs:
!c    -------
	INTEGER ngrid,nlayer,nq,it,ij,i
      REAL ptimestep
      real zdtime
      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
      REAL pphi(ngrid,nlayer)
      REAL pphis(ngrid)
      REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
      REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
      REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)

!c   dynamial tendencies
      REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq)
      REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer)
      REAL pw(ngrid,nlayer)

!c   Time
      real rjourvrai
      REAL gmtime

!c     outputs:
!c     --------

!c   physical tendencies
      REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
      REAL pdpsrf(ngrid)


!c    Local variables :
!c    -----------------

      INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,unit,isoil
      REAL zps_av
      REAL zday
      REAL zh(ngrid,nlayer),z1,z2
      REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
      REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
      REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
      REAL zflubid(ngrid),zpmer(ngrid)
      REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
      REAL zdum1(ngrid,nlayer)
      REAL zdum2(ngrid,nlayer)
      REAL zdum3(ngrid,nlayer)
      REAL ztim1,ztim2,ztim3
      REAL zls,zinsol
      REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
      REAL zfluxsw(ngrid),zfluxlw(ngrid)
      character*2 str2

!c   Local saved variables:
!c   ----------------------
	REAL(rstd)::long(ngrid),lati(ngrid) 
	REAL(rstd)::mu0(ngrid),fract(ngrid),coslat(ngrid)
        REAL(rstd)::sinlon(ngrid),coslon(ngrid),sinlat(ngrid)
	REAL(rstd)::dist_sol,declin
        REAL::totarea !sarvesh  
!!!!!!!!sarvesh !!!!!!! CHECK SAVE ATTRIBUTE 
      INTEGER:: icount
      real:: zday_last
      REAL:: solarcst

      SAVE icount,zday_last
      SAVE solarcst
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      REAL stephan
      SAVE stephan
      DATA stephan/5.67e-08/
      DATA solarcst/1370./
      REAL presnivs(nlayer)
	INTEGER:: nn1,nn2
c-----------------------------------------------------------------------
c    1. Initialisations :
c    --------------------
c     call initial0(ngrid*nlayer*nqmx,pdq)

!	nn1=(jj_begin -1)*iim+ii_begin
!	nn2=(jj_end -1)*iim+ii_end

      nlevel=nlayer+1
      igout=ngrid/2+1

	 DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
           ig=(j-1)*iim+i
	   sinlat(ig) = sin(lati(ig))
	   coslat(ig) = cos(lati(ig))
	   sinlon(ig) = sin(long(ig))
	   coslon(ig) = cos(long(ig))
	  END DO
	 ENDDO 
        zday=rjourvrai+gmtime
	IF ( it == 0 ) then 
	 firstcall=.TRUE.
	ELSE
	 firstcall=.FALSE. 
	ENDIF 
	IF ( it == ndays*day_step ) Then
	 lastcall = .True.
        END IF 
	 
        IF(firstcall) THEN
        PRINT*,'FIRSTCALL  ',ngridmx,nlayermx,nsoilmx
         zday_last=rjourvrai
         inertie=2000
         albedo=0.2
         emissiv=1.
         z0=0.1
         rnatur=1.
         q2=1.e-10
         q2l=1.e-10
         tsurf(:)=300.
         tsoil(:,:)=300.
!         print*,tsoil(ngrid/2+1,nsoilmx/2+2)
!         print*,'TS ',tsurf(igout),tsoil(igout,5)
 
          IF (.not.callrad) then
	    DO j=jj_begin-offset,jj_end+offset
             DO i=ii_begin-offset,ii_end+offset
              ig=(j-1)*iim+i
              fluxrad(ig)=0.
             enddo
            enddo
           ENDIF

         IF(callsoil) THEN
             CALL soil(ngrid,nsoilmx,firstcall,inertie,
     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
          ELSE
            PRINT*,'WARNING!!! Thermal conduction in the soil
     s      turned off'
	    DO j=jj_begin-offset,jj_end+offset
             DO i=ii_begin-offset,ii_end+offset
               ig=(j-1)*iim+i
               capcal(ig)=1.e5
               fluxgrd(ig)=0.
             ENDDO
           ENDDO
          ENDIF
            icount=0
         ENDIF

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

       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
            pdv(ig,l)=0.
            pdu(ig,l)=0.
            pdt(ig,l)=0.
          ENDDO
         ENDDO
        ENDDO

         DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
          ig=(j-1)*iim+i
          pdpsrf(ig)=0.
          zflubid(ig)=0.
          zdtsrf(ig)=0.
         ENDDO
	ENDDO

        zps_av=0.0
         DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
           ig=(j-1)*iim+i
           zps_av=zps_av+pplev(ig,1)*Ai(ig)
           totarea=totarea+Ai(ig)
          END DO
         END DO
	zps_av=zps_av/totarea


	!print*,"maxplev",maxval(pplev(:,1)),minval(pplev(:,1))
c-----------------------------------------------------------------------
c   calcul du geopotentiel aux niveaux intercouches
c   ponderation des altitudes au niveau des couches en dp/p

       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
          zzlay(ig,l)=pphi(ig,l)/g
         ENDDO
        ENDDO
       ENDDO

	!print*,"zzlay",maxval(zzlay(:,1)),minval(zzlay(:,1))


        DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
          ig=(j-1)*iim+i
          zzlev(ig,1)=0.
         ENDDO
        ENDDO

      DO l=2,nlayer
         DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
          ig=(j-1)*iim+i
          z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
          z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
         ENDDO
        ENDDO
       ENDDO
       !print*,"zzlev",maxval(zzlev(:,1)),minval(zzlev(:,1))
c-----------------------------------------------------------------------
c   Transformation de la temperature en temperature potentielle
        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
           zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**kappa
          ENDDO
         ENDDO
        ENDDO

        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
           zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
          ENDDO
         ENDDO
	ENDDO
	!print*,"ph pot",maxval(zh(:,1)),minval(zh(:,1))
!	go to 101
c-----------------------------------------------------------------------
c    2. Calcul of the radiative tendencies :
c    ---------------------------------------
!      print*,'callrad0'
      IF(callrad) THEN
!         print*,'callrad'
!   WARNING !!! on calcule le ray a chaque appel
!        IF( MOD(icount,iradia).EQ.0) THEN

           CALL solarlong(zday,zls)
           CALL orbite(zls,dist_sol,declin)
!   
!      print*,'ATTENTIOn : pdeclin = 0','  L_s=',zls
!      print*,'diurnal=',diurnal
       IF(diurnal) THEN
         if ( 1.eq.1 ) then
               ztim1=SIN(declin)
               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))

               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
     s         ztim1,ztim2,ztim3,
     s         mu0,fract)
          else
               zdtime=ptimestep*float(iradia)
!               CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract) ! FIXME
              !print*,'ZENANG '
         endif

          IF(lverbose) THEN
                  PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
                  PRINT*,zday, declin,
     s            sinlon(igout),coslon(igout),
     s            sinlat(igout),coslat(igout)
           ENDIF
        ELSE
             !print*,'declin,ngrid,radius',declin,ngrid,radius
            CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,radius)
!	      open(100,file="mu0.txt")
!	      write(100,*)(mu0(ij),ij=1,ngrid) 
        ENDIF
!	print*,"iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"
!c    2.2 Calcul of the radiative tendencies and fluxes:
!c    --------------------------------------------------
!c  2.1.2 levels

            zinsol=solarcst/(dist_sol*dist_sol)
!            print*,'iim,jjm,llm,ngrid,nlayer,ngrid,nlayer'
!            print*,iim,jjm,llm,ngrid,nlayer,ngrid,nlayer
!		print*,"zinsol sol_dist",zinsol,dist_sol
!	STOP 

            CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,
     $              pplev,zps_av,
     $              mu0,fract,zinsol,
     $              zfluxsw,zdtsw,
     $              lverbose)

	!print*,"sw",maxval(zfluxsw),minval(zfluxsw),
!	$            maxval(zdtsw),minval(zdtsw), "   it",it

!	print*,"lllllllllllllllllllllllllllllllllllllllll"
!	print*,"pplev",maxval(pplev),minval(pplev)
!	print*,"zps,tsurf",zps_av,maxval(tsurf),minval(tsurf)
!	print*,"pt",maxval(pt),minval(pt) 


            CALL lw(ngrid,nlayer,coefir,emissiv,
     $             pplev,zps_av,tsurf,pt,
     $             zfluxlw,zdtlw,
     $             lverbose)

	!print*,"lw",maxval(zfluxlw),minval(zfluxlw),
!	$            maxval(zdtlw),minval(zdtlw)

!	print*,"lw",maxval(

!	print*,"after lw"
c    2.4 total flux and tendencies:
c    ------------------------------

c    2.4.1 fluxes

          DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
               ig=(j-1)*iim+i
               fluxrad(ig)=emissiv(ig)*zfluxlw(ig)
     $         +zfluxsw(ig)*(1.-albedo(ig))

               zplanck(ig)=tsurf(ig)*tsurf(ig)

               zplanck(ig)=emissiv(ig)*
     $         stephan*zplanck(ig)*zplanck(ig)

               fluxrad(ig)=fluxrad(ig)-zplanck(ig)
            ENDDO
          ENDDO
c    2.4.2 temperature tendencies

            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
                  dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
                 ENDDO
               ENDDO
             ENDDO


c    2.5 Transformation of the radiative tendencies:
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
              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
            ENDDO
           ENDDO
          ENDDO

         IF(lverbose) THEN
            PRINT*,'Diagnotique for the radaition'
            PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
            PRINT*,albedo(igout),emissiv(igout),mu0(igout),
     s           fract(igout),
     s           fluxrad(igout),zplanck(igout)
            PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
!            PRINT*,'unjours',unjours
            DO l=1,nlayer
               WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),
     s         pplay(igout,l),pplev(igout,l),
     s         zdtsw(igout,l),zdtlw(igout,l)
            ENDDO
          ENDIF
        ENDIF   !( CALL RADIATION ) 
!	print*,"eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee"

c-----------------------------------------------------------------------
c    3. Vertical diffusion (turbulent mixing):
c    -----------------------------------------
c
       IF(calldifv) THEN

         DO ig=1,ngrid
            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
         ENDDO

         CALL zerophys(ngrid*nlayer,zdum1)
         CALL zerophys(ngrid*nlayer,zdum2)
c        CALL zerophys(ngrid*nlayer,zdum3)
         do l=1,nlayer
            do ig=1,ngrid
               zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
            enddo
         enddo

         CALL vdif(ngrid,nlayer,zday,
     $        ptimestep,capcal,z0,
     $        pplay,pplev,zzlay,zzlev,
     $        pu,pv,zh,tsurf,emissiv,
     $        zdum1,zdum2,zdum3,zflubid,
     $        zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,
     $        lverbose)
c        igout=ngrid/2+1
c        PRINT*,'zdufr zdvfr zdhfr'
c        DO l=1,nlayer
c           PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
c        ENDDO
c        CALL difv  (ngrid,nlayer,
c    $                  area,lati,capcal,
c    $                  pplev,pphi,
c    $                  pu,pv,zh,tsurf,emissiv,
c    $                  zdum1,zdum2,zdum3,zflubid,
c    $                  zdufr,zdvfr,zdhfr,zdtsrf,
c    $                  .true.)
c        PRINT*,'zdufr zdvfr zdhfr'
c        DO l=1,nlayer
c           PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
c        ENDDO
c        STOP

         DO l=1,nlayer
            DO ig=1,ngrid
               pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
               pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
               pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
            ENDDO
         ENDDO

         DO ig=1,ngrid
            zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
         ENDDO

       ELSE

         DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
            ig=(j-1)*iim+i
            zdtsrf(ig)=zdtsrf(ig)+
     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
          ENDDO
         ENDDO

!       write(66,*)"tsrf",maxval(zdtsrf(:)),minval(zdtsrf(:))
!       write(66,*)"frd",maxval(fluxrad(:)),minval(fluxrad(:))
!       write(66,*)"fgd",maxval(fluxgrd(:)),minval(fluxgrd(:))
        ENDIF
c-----------------------------------------------------------------------
c   4. Dry convective adjustment:
c   -----------------------------

      IF(calladj) THEN

         DO l=1,nlayer
            DO ig=1,ngrid
               zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
            ENDDO
         ENDDO
         CALL zerophys(ngrid*nlayer,zdufr)
         CALL zerophys(ngrid*nlayer,zdvfr)
         CALL zerophys(ngrid*nlayer,zdhfr)
         CALL convadj(ngrid,nlayer,ptimestep,
     $                pplay,pplev,zpopsk,
     $                pu,pv,zh,
     $                pdu,pdv,zdum1,
     $                zdufr,zdvfr,zdhfr)

         DO l=1,nlayer
            DO ig=1,ngrid
               pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
               pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
               pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
            ENDDO
         ENDDO

        ENDIF

!101		continue
c-----------------------------------------------------------------------
c   On incremente les tendances physiques de la temperature du sol:
c   ---------------------------------------------------------------
!       WRITE(55,*)"tsurf",maxval(tsurf(:)),minval(tsurf(:)),it

	  DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
             ig=(j-1)*iim+i
             tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 
           ENDDO
         ENDDO
c-----------------------------------------------------------------------
c   soil temperatures:
c   --------------------

       IF (callsoil) THEN
         CALL soil(ngrid,nsoilmx,.false.,inertie,
     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
         IF(lverbose) THEN
            PRINT*,'Surface Heat capacity,conduction Flux, Ts,
     s          dTs, dt'
            PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout),
     s          zdtsrf(igout),ptimestep
         ENDIF
       ENDIF

c-----------------------------------------------------------------------
c   sorties:
c   --------
      if(zday.GT.zday_last+period_sort) then
       zday_last=zday
c  Ecriture/extension de la coordonnee temps

         do ig=1,ngrid
            zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(kappa*cpp*285.))
         enddo
       endif
c-----------------------------------------------------------------------
      IF(lastcall) THEN
       PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
      ENDIF

      icount=icount+1
!      RETURN
      END SUBROUTINE phyparam_lmd
	END MODULE PHY