SUBROUTINE phyparam(ngrid,nlayer,nq,
     s            firstcall,lastcall,
     s            rjourvrai,gmtime,ptimestep,
     s            pplev,pplay,pphi,pphis,presnivs,
     s            pu,pv,pt,pq,
     s            pw,
     s            pdu,pdv,pdt,pdq,pdpsrf)

      USE comsaison
      USE dimphy
      USE comgeomfi
      USE phys_const
      USE planet
      USE astronomy
      USE vdif_mod, ONLY : vdif
      USE solar, ONLY : solang
      USE radiative, ONLY : mucorr, sw
      USE radiative_lw, ONLY : lw

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,nlayermx)        \
c    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
c    pdt(ngrid,nlayermx)         /  variables due to physical processes.
c    pdq(ngrid,nlayermx)      /
c    pdpsrf(ngrid)        /
c
c=======================================================================
c
c-----------------------------------------------------------------------
c
c    0.  Declarations :
c    ------------------

#include "dimensions.h"
#include "callkeys.h"
#include "surface.h"

c    Arguments :
c    -----------

c    inputs:
c    -------
      INTEGER ngrid,nlayer,nq

      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)
      LOGICAL firstcall,lastcall

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

      INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq
      INTEGER*4 day_ini
c
      REAL,DIMENSION(ngrid) :: mu0,fract
      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
      REAL factq(nq),tauq(nq)
      character*3 nomq

c   Local saved variables:
c   ----------------------

      INTEGER, SAVE :: icount
      real, SAVE :: zday_last
!$OMP THREADPRIVATE( icount,zday_last)

      REAL zps_av

      real,allocatable,SAVE :: tsurf(:),tsoil(:,:),rnatur(:)
      real,allocatable,SAVE :: capcal(:),fluxgrd(:)
      real,allocatable,SAVE :: dtrad(:,:),fluxrad(:)
      real,allocatable,SAVE ::  q2(:,:),q2l(:,:)
      real,allocatable,SAVE ::  albedo(:),emissiv(:),z0(:),inertie(:)
      real,SAVE :: solarcst=1370.
      real,SAVE :: stephan=5.67e-08

!$OMP THREADPRIVATE(tsurf,tsoil,rnatur)
!$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad)
!$OMP THREADPRIVATE( q2,q2l)
!$OMP THREADPRIVATE( albedo,emissiv,solarcst,z0,inertie)
!$OMP THREADPRIVATE( stephan)


      EXTERNAL convadj
      EXTERNAL ismin,ismax


      INTEGER        longcles
      PARAMETER    ( longcles = 20 )
      REAL clesphy0( longcles      )
      REAL presnivs(nlayer)



      print*,'OK DANS PHYPARAM'

c-----------------------------------------------------------------------
c    1. Initialisations :
c    --------------------

!     call initial0(ngrid*nlayermx*nqmx,pdq)
      nlevel=nlayer+1

!     print*,'OK ',nlevel

      igout=ngrid/2+1
!     print*,'OK PHYPARAM ',ngrid,igout


      zday=rjourvrai+gmtime

!     print*,'OK PHYPARAM 0A nq ',nq
      tauq(1)=1800.
      tauq(2)=10800.
      tauq(3)=86400.
      tauq(4)=864000.
!     print*,'OK PHYPARAM 0 B nq ',nq
      factq(1:4)=(1.-exp(-ptimestep/tauq(1:4)))/ptimestep

!     print*,'OK PHYPARAM 0 '
      print*,'nq ',nq
      print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
      print*,'nlayer',nlayer
      print*,'size pdq ',ngrid*nlayer*4,ngrid*nlayer*nq,
     s      size(pdq),size(lati),size(pq),size(factq)
      do iq=1,4
         do l=1,nlayer
             pdq(1:ngrid,l,iq)=
     &      (1.+sin(lati(1:ngrid))-pq(1:ngrid,l,iq))*factq(iq)
         enddo
      enddo

      IF(firstcall) THEN

      print*,'AKk',ngrid,nsoilmx
      allocate(tsurf(ngrid))
      print*,'AKa'
      allocate (tsoil(ngrid,nsoilmx))
      print*,'AKb'
      allocate (rnatur(ngrid))
      print*,'AK2'
      allocate(capcal(ngrid),fluxgrd(ngrid))
      print*,'AK3'
      allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
      print*,'AK4'
      allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
      print*,'AK5'
      allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
      print*,'AK6'


         do l=1,nlayer
            pdq(:,l,5)=1.+sin(lati(:))/ptimestep
         enddo
         PRINT*,'FIRSTCALL  '

!         zday_last=rjourvrai
         zday_last=zday-ptimestep/unjours
c        CALL readfi(ngrid,nlayer,nsoilmx,ldrs,
c    .      day_ini,gmtime,albedo,inertie,emissiv,z0,rnatur,
c    .      q2,q2l,tsurf,tsoil)
         rnatur=1.
         emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
         inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter
         albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
         z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
         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)
         CALL iniorbit

         if (.not.callrad) then
            do ig=1,ngrid
               fluxrad(ig)=0.
            enddo
         endif

!     print*,'OK PHYPARAM 1 '
         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 ig=1,ngrid
               capcal(ig)=1.e5
               fluxgrd(ig)=0.
            ENDDO
         ENDIF
c        CALL inifrict(ptimestep)
         icount=0

         PRINT*,'FIRSTCALL B '
          print*,'INIIO avant iophys_ini '
          call iophys_ini('phys.nc    ',presnivs)

      ENDIF
      icount=icount+1

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

      DO l=1,nlayer
         DO ig=1,ngrid
            pdv(ig,l)=0.
            pdu(ig,l)=0.
            pdt(ig,l)=0.
         ENDDO
      ENDDO

      DO ig=1,ngrid
         pdpsrf(ig)=0.
         zflubid(ig)=0.
         zdtsrf(ig)=0.
      ENDDO

c-----------------------------------------------------------------------
c   calcul du geopotentiel aux niveaux intercouches
c   ponderation des altitudes au niveau des couches en dp/p

      DO l=1,nlayer
         DO ig=1,ngrid
            zzlay(ig,l)=pphi(ig,l)/g
         ENDDO
      ENDDO
      DO ig=1,ngrid
         zzlev(ig,1)=0.
      ENDDO
      DO l=2,nlayer
         DO ig=1,ngrid
            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

c-----------------------------------------------------------------------
c   Transformation de la temperature en temperature potentielle
      DO l=1,nlayer
         DO ig=1,ngrid
            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
         ENDDO
      ENDDO
      DO l=1,nlayer
         DO ig=1,ngrid
            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c    2. Calcul of the radiative tendencies :
c    ---------------------------------------

!        print*,'callrad0'
      IF(callrad) THEN
!     print*,'callrad'

c   WARNING !!! on calcule le ray a chaque appel
c        IF( MOD(icount,iradia).EQ.0) THEN

            CALL solarlong(zday,zls)
            CALL orbite(zls,dist_sol,declin)
c      declin=0.
!     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))
c           call dump2d(iim,jjm-1,lati(2),'LATI  0   ')
c           call dump2d(iim,jjm-1,long(2),'LONG  0   ')
c           call dump2d(iim,jjm-1,sinlon(2),'sinlon0   ')
c           call dump2d(iim,jjm-1,coslon(2),'coslon0   ')
c           call dump2d(iim,jjm-1,sinlat(2),'sinlat   ')
c           call dump2d(iim,jjm-1,coslat(2),'coslat   ')
               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)
           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,rad',declin,ngrid,rad

c           call dump2d(iim,jjm-1,lati(2),'LATI      ')
               CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
            ENDIF
c           call dump2d(iim,jjm-1,fract(2),'FRACT A   ')
c           call dump2d(iim,jjm-1,mu0(2),'MU0 A     ')


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,ngridmax,nlayer
            print*,'iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer'
c           call dump2d(iim,jjm-1,albedo(2),'ALBEDO    ')
c           call dump2d(iim,jjm-1,mu0(2),'MU0       ')
c           call dump2d(iim,jjm-1,fract(2),'FRACT     ')
c           call dump2d(iim,jjm-1,lati(2),'LATI      ')
            zps_av=1.e5
            if (firstcall) then
               print*,'WARNING ps_rad impose'
            endif
            CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,
     $              pplev,zps_av,
     $              mu0,fract,zinsol,
     $              zfluxsw,zdtsw,
     $              lverbose)
c           call dump2d(iim,jjm-1,zfluxsw(2),'SWS 1     ')
c           stop

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


c    2.4 total flux and tendencies:
c    ------------------------------

c    2.4.1 fluxes

            DO ig=1,ngrid
               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

c    2.4.2 temperature tendencies

            DO l=1,nlayer
               DO ig=1,ngrid
                  dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
               ENDDO
            ENDDO

c        ENDIF

c    2.5 Transformation of the radiative tendencies:
c    -----------------------------------------------

         DO l=1,nlayer
            DO ig=1,ngrid
               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
            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

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

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

         zdum1(:,:)=0.
         zdum2(:,:)=0.

         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 ig=1,ngrid
            zdtsrf(ig)=zdtsrf(ig)+
     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
         ENDDO
      ENDIF
c
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

         zdufr(:,:)=0.
         zdvfr(:,:)=0.
         zdhfr(:,:)=0.

         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

c-----------------------------------------------------------------------
c   On incremente les tendances physiques de la temperature du sol:
c   ---------------------------------------------------------------

      DO ig=1,ngrid
         tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 
      ENDDO

      WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)

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   --------

c           call dump2d(iim,jjm-1,zfluxsw(2),'SWS 2     ')
      print*,'zday, zday_last ',zday,zday_last,icount
      if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
      print*,'zday, zday_last SORTIE ',zday,zday_last
       zday_last=zday
c  Ecriture/extension de la coordonnee temps

         do ig=1,ngridmax
            zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*285.))
         enddo

       call iophys_ecrit('u',llm,'Vent zonal moy','m/s',pu)
       call iophys_ecrit('v',llm,'Vent meridien moy','m/s',pv)
       call iophys_ecrit('temp',llm,'Temperature','K',pt)
       call iophys_ecrit('geop',llm,'Geopotential','m2/s2',pphi)
       call iophys_ecrit('plev',llm,'plev','Pa',pplev(:,1:nlayer))

       call iophys_ecrit('du',llm,'du',' ',pdu)
       call iophys_ecrit('dv',llm,'du',' ',pdv)
       call iophys_ecrit('dt',llm,'du',' ',pdt)
       call iophys_ecrit('dtsw',llm,'dtsw',' ',zdtsw)
       call iophys_ecrit('dtlw',llm,'dtlw',' ',zdtlw)

       do iq=1,nq
          nomq="tr."
          write(nomq(2:3),'(i1.1)') iq
          call iophys_ecrit(nomq,llm,nomq,' ',pq(:,:,iq))
       enddo

       call iophys_ecrit('ts',1,'Surface temper','K',tsurf)
       call iophys_ecrit('coslon',1,'coslon',' ',coslon)
       call iophys_ecrit('sinlon',1,'sinlon',' ',sinlon)
       call iophys_ecrit('coslat',1,'coslat',' ',coslat)
       call iophys_ecrit('sinlat',1,'sinlat',' ',sinlat)
       call iophys_ecrit('mu0',1,'mu0',' ',mu0)
       call iophys_ecrit('alb',1,'alb',' ',albedo)
       call iophys_ecrit('fract',1,'fract',' ',fract)
       call iophys_ecrit('ps',1,'Surface pressure','Pa',pplev)
       call iophys_ecrit('slp',1,'Sea level pressure','Pa',zpmer)
       call iophys_ecrit('swsurf',1,'SW surf','Pa',zfluxsw)
       call iophys_ecrit('lwsurf',1,'LW surf','Pa',zfluxlw)

      endif

c-----------------------------------------------------------------------
      IF(lastcall) THEN
         call iotd_fin
         PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
!        if (ierr.ne.0) then
!          write(6,*)' Pb d''ouverture du fichier restart'
!          write(6,*)' ierr = ', ierr
!          call exit(1)
!        endif
         write(75)  tsurf,tsoil
c    s             (tsurf(1),ig=1,iim+1),
c    s             ( (tsurf(ig),ig=(j-2)*iim+2,(j-1)*iim+1),
c    s              tsurf((j-2)*iim+2) ,j=2,jjm),
c    s              (tsurf(ngridmax),ig=1,iim+1),
c    s         (   (tsoil(1,l),ig=1,iim+1),
c    s             ( (tsoil(ig,l),ig=(j-2)*iim+2,(j-1)*iim+1),
c    s              tsoil((j-2)*iim+2,l) ,ig=2,jjm),
c    s              (tsoil(ngridmax,l),ig=1,iim+1)
c    s          ,l=1,nsoilmx)
      ENDIF


      RETURN
      END