      SUBROUTINE SISVAT_TS2
c #ES.                     (ETSo_0,ETSo_1,ETSo_d)                               
                                                                                
C +------------------------------------------------------------------------+    
C | MAR          SISVAT_TS2                            Mon 16-08-2009  MAR |    
C |   SubRoutine SISVAT_TS2 computes the Soil/Snow temperature and fluxes  |
C |   using the same method as in LMDZ for consistency.                    |
C |   The corresponding LMDZ routines are soil (soil.F90) and calcul_fluxs |
C |   (calcul_fluxs_mod.F90).                                              |   
C +------------------------------------------------------------------------+    
C |                                                                        |    
C |                                                                        |    
C |   PARAMETERS:  klonv: Total Number of columns =                        |    
C |   ^^^^^^^^^^        = Total Number of grid boxes of surface type       |    
C |                       (land ice for now)                               |    
C |                                                                        |    
C |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |    
C |   ^^^^^    sol_SV   : Downward Solar Radiation                  [W/m2] |    
C |            IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |       
C |            VV__SV   : SBL Top    Wind Speed                      [m/s] |    
C |            TaT_SV   : SBL Top    Temperature                       [K] |     
C |            QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |      
C |            dzsnSV   : Snow Layer Thickness                         [m] |        
C |            dt__SV   : Time Step                                    [s] |    
C |                                                                        |    
C |            SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |    
C |            Eso_sv   : Soil+Snow       Emissivity                   [-] |    
C |   ?        rah_sv   : Aerodynamic Resistance for Heat            [s/m] |    
C |                                                                        | 
C |            dz1_SV    : "inverse" layer thickness (centered)      [1/m] |
C |            dz2_SV    : layer thickness (layer above (?))           [m] |
C |            AcoHSV    : coefficient for Enthalpy evolution, from atm.   |
C |            AcoHSV    : coefficient for Enthalpy evolution, from atm.   |
C |            AcoQSV    : coefficient for Humidity evolution, from atm.   | 
C |            BcoQSV    : coefficient for Humidity evolution, from atm.   |
C |            ps__SV    : surface pressure                           [Pa] |
C |            p1l_SV    : 1st atmospheric layer pressure             [Pa] |     
C |            cdH_SV    : drag coeff Energy (?)                           |
C |            rsolSV    : Radiation balance surface                [W/m2] |
C |            lambSV    : Coefficient for soil layer geometry         [-] |
C |                                                                        |    
C |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|    
C |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |    
C |   ^^^^^^   rsolSV   : Radiation balance surface                 [W/m2] |     
C |                                                                        |    
C |   OUTPUT:  IRs_SV   : Soil      IR Radiation                    [W/m2] |    
C |   ^^^^^^   HSs_sv   : Sensible  Heat Flux                       [W/m2] |    
C |            HLs_sv   : Latent    Heat Flux                       [W/m2] |  
C |            TsfnSV   : new surface temperature                      [K] |
C |            Evp_sv   : Evaporation                              [kg/m2] |
C |            dSdTSV   : Sensible Heat Flux temp. derivative     [W/m2/K] |
C |            dLdTSV   : Latent Heat Flux temp. derivative       [W/m2/K] |
C |                                                                        | 
C |   ?        ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |    
C |   ?        ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |    
C |   ?        ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |    
C |                                                                        | 
C |________________________________________________________________________|

      USE VAR_SV
      USE VARdSV
 
      USE VARySV                      
      USE VARtSV                      
      USE VARxSV 
      USE VARphy
      USE indice_sol_mod


      IMPLICIT NONE                                                                    
                                                                                
                                                                                
C +--Global Variables                                                           
C +  ================                                                          

      INCLUDE "YOMCST.h"
      INCLUDE "YOETHF.h"
      INCLUDE "FCTTRE.h"
!      INCLUDE "indicesol.h"
      INCLUDE "comsoil.h"
!      include  "LMDZphy.inc"   
                                                                        
C +--OUTPUT for Stand Alone NetCDF File                                         
C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                                         
c #NC      real*8        SOsoKL(klonv)             ! Absorbed Solar Radiation        
c #NC      real*8        IRsoKL(klonv)             ! Absorbed IR    Radiation        
c #NC      real*8        HSsoKL(klonv)             ! Absorbed Sensible Heat Flux     
c #NC      real*8        HLsoKL(klonv)             ! Absorbed Latent   Heat Flux     
c #NC      real*8        HLs_KL(klonv)             ! Evaporation                     
c #NC      real*8        HLv_KL(klonv)             ! Transpiration                   
c #NC      common/DumpNC/SOsoKL,IRsoKL                                               
c #NC     .             ,HSsoKL,HLsoKL                                               
c #NC     .             ,HLs_KL,HLv_KL          

C +--Internal Variables                                                         
C +  ==================    

      integer ig,jk,isl
      real mu     
      real Tsrf(klonv)               ! surface temperature as extrapolated from soil
      real mug(klonv)                 !hj coef top layers 
      real ztherm_i(klonv),zdz2(klonv,-nsol:nsno),z1s
      real pfluxgrd(klonv), pcapcal(klonv), cal(klonv)
      real beta(klonv), dif_grnd(klonv)
      real C_coef(klonv,-nsol:nsno),D_coef(klonv,-nsol:nsno)

      REAL, DIMENSION(klonv)   :: zx_mh, zx_nh, zx_oh
      REAL, DIMENSION(klonv)   :: zx_mq, zx_nq, zx_oq
      REAL, DIMENSION(klonv)   :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
      REAL, DIMENSION(klonv)   :: zx_sl, zx_k1
      REAL, DIMENSION(klonv)   :: d_ts
      REAL                     :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
      REAL                     :: qsat_new, q1_new
C      REAL, PARAMETER          :: t_grnd = 271.35, t_coup = 273.15
C      REAL, PARAMETER          :: max_eau_sol = 150.0
      REAL, DIMENSION(klonv)   :: IRs__D, dIRsdT


      REAL t_grnd                      ! not used
      parameter(t_grnd = 271.35)       !
      REAL t_coup                      ! distinguish evap/sublimation
      parameter(t_coup = 273.15)       !
      REAL max_eau_sol
      parameter(max_eau_sol = 150.0)


!        write(*,*)'T check'

!        DO  ig = 1,knonv  
!            DO  jk = 1,isnoSV(ig) !nsno       
!              IF (TsisSV(ig,jk) <= 1.) THEN          !hj check
!                TsisSV(ig,jk) = TsisSV(ig,isnoSV(ig))
!              ENDIF  

!              IF (TsisSV(ig,jk) <= 1.) THEN          !hj check
!                TsisSV(ig,jk) = 273.15
!              ENDIF   
!            END DO
!        END DO     

C!=======================================================================
C! I. First part: corresponds to soil.F90 in LMDZ
C!=======================================================================

      DO ig = 1,knonv      
        DO jk =1,isnoSV(ig)
          dz2_SV(ig,jk)=dzsnSV(ig,jk)
C! use arithmetic center between layers to derive dz1 for snow layers for simplicity: 
          dz1_SV(ig,jk)=2./(dzsnSV(ig,jk)+dzsnSV(ig,jk-1)) 
        ENDDO
      ENDDO

      DO ig = 1,knonv
        ztherm_i(ig)   = inertie_lic
        IF (isnoSV(ig) > 0) ztherm_i(ig)   = inertie_sno
      ENDDO

C!-----------------------------------------------------------------------
C! 1)
C! Calculation of Cgrf and Dgrd coefficients using soil temperature from 
C! previous time step.
C!
C! These variables are recalculated on the local compressed grid instead 
C! of saved in restart file.
C!-----------------------------------------------------------------------
      DO ig=1,knonv
        DO jk=-nsol,nsno
          zdz2(ig,jk)=dz2_SV(ig,jk)/dt__SV                                !ptimestep
        ENDDO
      ENDDO
 
      DO ig=1,knonv
        z1s = zdz2(ig,-nsol)+dz1_SV(ig,-nsol+1)
        C_coef(ig,-nsol+1)=zdz2(ig,-nsol)*TsisSV(ig,-nsol)/z1s
        D_coef(ig,-nsol+1)=dz1_SV(ig,-nsol+1)/z1s
      ENDDO

      DO ig=1,knonv
        DO jk=-nsol+1,isnoSV(ig)-1,1
          z1s = 1./(zdz2(ig,jk)+dz1_SV(ig,jk+1)+dz1_SV(ig,jk)           &
     &          *(1.-D_coef(ig,jk)))
          C_coef(ig,jk+1)=                                              &
     &          (TsisSV(ig,jk)*zdz2(ig,jk)                              &
     &          +dz1_SV(ig,jk)*C_coef(ig,jk)) * z1s
          D_coef(ig,jk+1)=dz1_SV(ig,jk+1)*z1s
        ENDDO
      ENDDO

C!-----------------------------------------------------------------------
C! 2)
C! Computation of the soil temperatures using the Cgrd and Dgrd
C! coefficient computed above
C!
C!-----------------------------------------------------------------------
C! Extrapolate surface Temperature                   !hj check
      mu=1./((2.**1.5-1.)/(2.**(0.5)-1.)-1.)

!     IF (knonv>0) THEN 
!      DO ig=1,8 
!        write(*,*)ig,'sisvat: Tsis ',TsisSV(ig,isnoSV(ig))
!        write(*,*)'max-1            ',TsisSV(ig,isnoSV(ig)-1)
!        write(*,*)'max-2            ',TsisSV(ig,isnoSV(ig)-2)
!        write(*,*)'0                ',TsisSV(ig,0)
!!        write(*,*)min(max(isnoSV(ig),0),1),max(1-isnoSV(ig),0)
!      ENDDO
!     END IF

      DO ig=1,knonv    
        IF (isnoSV(ig)>0) THEN
          IF (isnoSV(ig)>1) THEN
           mug(ig)=1./(1.+dzsnSV(ig,isnoSV(ig)-1)/dzsnSV(ig,isnoSV(ig))) !mu
          ELSE
           mug(ig) = 1./(1.+dzsnSV(ig,isnoSV(ig)-1)/dz_dSV(0)) !mu
          ENDIF
        ELSE
          mug(ig) = lambSV
        ENDIF

        IF (mug(ig)  <= 0.05) THEN
         write(*,*)'Attention mu low', mug(ig)
        ENDIF
        IF (mug(ig)  >= 0.98) THEN
         write(*,*)'Attention mu high', mug(ig)
        ENDIF

        Tsrf(ig)=(1.5*TsisSV(ig,isnoSV(ig))-0.5*TsisSV(ig,isnoSV(ig)-1))&
     &           *min(max(isnoSV(ig),0),1)+                             &
     &           ((mug(ig)+1)*TsisSV(ig,0)-mug(ig)*TsisSV(ig,-1))       &
     &           *max(1-isnoSV(ig),0)  
      ENDDO

  

C! Surface temperature
      DO ig=1,knonv
      TsisSV(ig,isnoSV(ig))=(mug(ig)*C_coef(ig,isnoSV(ig))+Tsf_SV(ig))/ &
     &        (mug(ig)*(1.-D_coef(ig,isnoSV(ig)))+1.)
      ENDDO
  
C! Other temperatures
      DO ig=1,knonv
        DO jk=isnoSV(ig),-nsol+1,-1
          TsisSV(ig,jk-1)=C_coef(ig,jk)+D_coef(ig,jk)                   &
     &          *TsisSV(ig,jk)
        ENDDO
      ENDDO
C      write(*,*)ig,'Tsis',TsisSV(ig,0)

C      IF (indice == is_sic) THEN
C        DO ig = 1,knonv
C          TsisSV(ig,-nsol) = RTT - 1.8
C        END DO
C      ENDIF

CC      !hj new 11 03 2010
        DO ig=1,knonv                                                          
          isl         = isnoSV(ig)                                             
C          dIRsdT(ig) = Eso_sv(ig)* SteBo   * 4.                        & ! - d(IR)/d(T)   
C     &                             * Tsf_SV(ig)                         & !T TsisSV(ig,isl)           !                
C     &                             * Tsf_SV(ig)                         & !TsisSV(ig,isl)           !                
C     &                             * Tsf_SV(ig) !TsisSV(ig,isl)           !                
C          IRs__D(ig) = dIRsdT(ig)* Tsf_SV(ig) * 0.75 !TsisSV(ig,isl) * 0.75    !: 
          dIRsdT(ig) = Eso_sv(ig)* StefBo   * 4.                        & ! - d(IR)/d(T)   
     &                             * TsisSV(ig,isl)                     & !
     &                             * TsisSV(ig,isl)                     & !
     &                             * TsisSV(ig,isl)                     & !
          IRs__D(ig) = dIRsdT(ig)* TsisSV(ig,isl) * 0.75                  !: 
         END DO
       !hj
C!-----------------------------------------------------------------------
C! 3)
C! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil 
C! temperature
C!-----------------------------------------------------------------------
      DO ig=1,knonv
        z1s = zdz2(ig,-nsol)+dz1_SV(ig,-nsol+1)
        C_coef(ig,-nsol+1) = zdz2(ig,-nsol)*TsisSV(ig,-nsol)/z1s
        D_coef(ig,-nsol+1) = dz1_SV(ig,-nsol+1)/z1s
      ENDDO

      DO ig=1,knonv  
        DO jk=-nsol+1,isnoSV(ig)-1,1
          z1s = 1./(zdz2(ig,jk)+dz1_SV(ig,jk+1)+dz1_SV(ig,jk)           &
     &          *(1.-D_coef(ig,jk)))
          C_coef(ig,jk+1) = (TsisSV(ig,jk)*zdz2(ig,jk)+                 &
     &          dz1_SV(ig,jk)*C_coef(ig,jk)) * z1s
          D_coef(ig,jk+1) = dz1_SV(ig,jk+1)*z1s
        ENDDO
      ENDDO

C!-----------------------------------------------------------------------
C! 4)
C! Computation of the surface diffusive flux from ground and
C! calorific capacity of the ground
C!-----------------------------------------------------------------------
      DO ig=1,knonv
C! (pfluxgrd)
        pfluxgrd(ig) = ztherm_i(ig)*dz1_SV(ig,isnoSV(ig))*              &
     &      (C_coef(ig,isnoSV(ig))+(D_coef(ig,isnoSV(ig))-1.)           &
     &      *TsisSV(ig,isnoSV(ig)))
C! (pcapcal)
        pcapcal(ig)  = ztherm_i(ig)*                                    &
     &      (dz2_SV(ig,isnoSV(ig))+dt__SV*(1.-D_coef(ig,isnoSV(ig)))    &
     &      *dz1_SV(ig,isnoSV(ig)))
        z1s = mug(ig)*(1.-D_coef(ig,isnoSV(ig)))+1.
        pcapcal(ig)  = pcapcal(ig)/z1s
        pfluxgrd(ig) = ( pfluxgrd(ig)                                   &
     &       + pcapcal(ig) * (TsisSV(ig,isnoSV(ig)) * z1s               &
     &       - mug(ig)* C_coef(ig,isnoSV(ig))                           &
     &       - Tsf_SV(ig))       /dt__SV )
      ENDDO
        
      
      cal(1:knonv) = RCPD / pcapcal(1:knonv)
      rsolSV(1:knonv)  = rsolSV(1:knonv) + pfluxgrd(1:knonv)
C!=======================================================================
C! II. Second part: corresponds to calcul_fluxs_mod.F90 in LMDZ
C!=======================================================================

      Evp_sv = 0.
c #NC HSsoKL=0.
c #NC HLsoKL=0.
      dSdTSV = 0.
      dLdTSV = 0.
	
      beta(:) = 1.0
      dif_grnd(:) = 0.0

C! zx_qs = qsat en kg/kg
C!**********************************************************************x***************

      DO ig = 1,knonv
        IF (ps__SV(ig)<1.) THEN
!          write(*,*)'ig',ig,'ps',ps__SV(ig)
          ps__SV(ig)=max(ps__SV(ig),1.e-8)
        ENDIF
        IF (p1l_SV(ig)<1.) THEN
!          write(*,*)'ig',ig,'p1l',p1l_SV(ig)
          p1l_SV(ig)=max(p1l_SV(ig),1.e-8)
        ENDIF
        IF (TaT_SV(ig)<180.) THEN
!          write(*,*)'ig',ig,'TaT',TaT_SV(ig)
          TaT_SV(ig)=max(TaT_SV(ig),180.)
        ENDIF
        IF (QaT_SV(ig)<1.e-8) THEN
!          write(*,*)'ig',ig,'QaT',QaT_SV(ig)
          QaT_SV(ig)=max(QaT_SV(ig),1.e-8)
        ENDIF
        IF (Tsf_SV(ig)<100.) THEN
!          write(*,*)'ig',ig,'Tsf',Tsf_SV(ig)
          Tsf_SV(ig)=max(Tsf_SV(ig),180.)
        ENDIF
        IF (Tsf_SV(ig)>500.) THEN
!          write(*,*)'ig',ig,'Tsf',Tsf_SV(ig)
          Tsf_SV(ig)=min(Tsf_SV(ig),400.)
        ENDIF
!         IF (Tsrf(ig).LT.1.) THEN
!!          write(*,*)'ig',ig,'Tsrf',Tsrf(ig)  
!           Tsrf(ig)=max(Tsrf(ig),TaT_SV(ig)-20.)
!         ENDIF
         IF (cdH_SV(ig)<1.e-10) THEN
!           IF (ig.le.3)   write(*,*)'ig',ig,'cdH',cdH_SV(ig)
           cdH_SV(ig)=.5
         ENDIF
      ENDDO  


      DO ig = 1,knonv
        zx_pkh(ig) = 1. ! (ps__SV(ig)/ps__SV(ig))**RKAPPA
        IF (thermcep) THEN
          zdelta=MAX(0.,SIGN(1.,rtt-Tsf_SV(ig)))
          zcvm5 = R5LES*LhvH2O*(1.-zdelta) + R5IES*LhsH2O*zdelta
          zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*QaT_SV(ig))
          zx_qs= r2es * FOEEW(Tsf_SV(ig),zdelta)/ps__SV(ig)
          zx_qs=MIN(0.5,zx_qs)
          !write(*,*)'zcor',retv*zx_qs
          zcor=1./(1.-retv*zx_qs)
          zx_qs=zx_qs*zcor
          zx_dq_s_dh = FOEDE(Tsf_SV(ig),zdelta,zcvm5,zx_qs,zcor)        &
     &          /LhvH2O / zx_pkh(ig)
        ELSE
          IF (Tsf_SV(ig)<t_coup) THEN
             zx_qs = qsats(Tsf_SV(ig)) / ps__SV(ig)
             zx_dq_s_dh = dqsats(Tsf_SV(ig),zx_qs)/LhvH2O               &
     &             / zx_pkh(ig)
          ELSE
             zx_qs = qsatl(Tsf_SV(ig)) / ps__SV(ig)
             zx_dq_s_dh = dqsatl(Tsf_SV(ig),zx_qs)/LhvH2O               &
     &             / zx_pkh(ig)
          ENDIF
        ENDIF
        zx_dq_s_dt(ig) = RCPD * zx_pkh(ig) * zx_dq_s_dh
        zx_qsat(ig) = zx_qs
C        zx_coef(ig) = cdH_SV(ig) *                                     &
C     &       (1.0+SQRT(u1lay(ig)**2+v1lay(ig)**2)) *                   &
C     &       p1l_SV(ig)/(RD*t1lay(ig))
        zx_coef(ig) = cdH_SV(ig) *                                      &
     &       (1.0+VV__SV(ig)) *                                         &
     &       p1l_SV(ig)/(RD*TaT_SV(ig))
       
      ENDDO


C! === Calcul de la temperature de surface ===
C! zx_sl = chaleur latente d'evaporation ou de sublimation
C!****************************************************************************************

      DO ig = 1,knonv
        zx_sl(ig) = LhvH2O
        IF (Tsf_SV(ig) < RTT) zx_sl(ig) = LhsH2O
        zx_k1(ig) = zx_coef(ig)
      ENDDO
    

      DO ig = 1,knonv
C! Q
        zx_oq(ig) = 1. - (beta(ig) * zx_k1(ig) * BcoQSV(ig) * dt__SV)
        zx_mq(ig) = beta(ig) * zx_k1(ig) *                              &
     &       (AcoQSV(ig) - zx_qsat(ig) +                                &
     &       zx_dq_s_dt(ig) * Tsf_SV(ig))                               &
     &       / zx_oq(ig)
        zx_nq(ig) = beta(ig) * zx_k1(ig) * (-1. * zx_dq_s_dt(ig))       &
     &       / zx_oq(ig)
       
C! H
        zx_oh(ig) = 1. - (zx_k1(ig) * BcoHSV(ig) * dt__SV)
        zx_mh(ig) = zx_k1(ig) * AcoHSV(ig) / zx_oh(ig)
        zx_nh(ig) = - (zx_k1(ig) * RCPD * zx_pkh(ig))/ zx_oh(ig)
     
C! surface temperature
        TsfnSV(ig) = (Tsf_SV(ig) + cal(ig)/RCPD * zx_pkh(ig) * dt__SV * &
     &       (rsolSV(ig) + zx_mh(ig) + zx_sl(ig) * zx_mq(ig))           & 
     &       + dif_grnd(ig) * t_grnd * dt__SV)/                         &
     &       ( 1. - dt__SV * cal(ig)/(RCPD * zx_pkh(ig)) *              &
     &       (zx_nh(ig) + zx_sl(ig) * zx_nq(ig))                        &  
     &       + dt__SV * dif_grnd(ig))

!hj rajoute 22 11 2010 tuning...
        TsfnSV(ig) = min(RTT+0.02,TsfnSV(ig))
       
        d_ts(ig) = TsfnSV(ig) - Tsf_SV(ig)


C!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
C!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
        Evp_sv(ig) = - zx_mq(ig) - zx_nq(ig) * TsfnSV(ig) 
        HLs_sv(ig) = - Evp_sv(ig) * zx_sl(ig)
        HSs_sv(ig) = zx_mh(ig) + zx_nh(ig) * TsfnSV(ig)
       
C! Derives des flux dF/dTs (W m-2 K-1):
        dSdTSV(ig) = zx_nh(ig)
        dLdTSV(ig) = zx_sl(ig) * zx_nq(ig)


!hj  new 11 03 2010                                      
        isl         = isnoSV(ig)
!        TsisSV(ig,isl) = TsfnSV(ig)                                         
        IRs_SV(ig) = IRs__D(ig)                                         &!     
     &                - dIRsdT(ig) * TsfnSV(ig) !TsisSV(ig,isl)?      !  

! hj
c #NC   SOsoKL(ig) = sol_SV(ig) * SoSosv(ig)              ! Absorbed Sol.    
c #NC   IRsoKL(ig) =               IRs_SV(ig)                           & !Up Surf. IR
c #NC&        +     tau_sv(ig)      *IRd_SV(ig)*Eso_sv(ig)              & !Down Atm IR
c #NC&        -(1.0-tau_sv(ig)) *0.5*IRv_sv(ig)            ! Down Veg IR  
c #NC   HLsoKL(ig) = HLs_sv(ig) 
c #NC   HSsoKL(ig) = HSs_sv(ig) 
c #NC   HLs_KL(ig) = Evp_sv(ig)

C! Nouvelle valeure de l'humidite au dessus du sol
        qsat_new=zx_qsat(ig) + zx_dq_s_dt(ig) * d_ts(ig)
        q1_new = AcoQSV(ig) - BcoQSV(ig)* Evp_sv(ig)*dt__SV
        QaT_SV(ig)=q1_new*(1.-beta(ig)) + beta(ig)*qsat_new 

      ENDDO

      end ! SUBROUTINE SISVAT_TS2
