!
MODULE surf_land_mod

  IMPLICIT NONE

CONTAINS
!
!****************************************************************************************
!  
  SUBROUTINE surf_land(itime, dtime, date0, jour, knon, knindex, &
       rlon, rlat, yrmu0, &
       debut, lafin, zlev, ccanopy, swnet, lwnet, albedo, &
       tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
       AcoefH, AcoefQ, BcoefH, BcoefQ, &
       AcoefU, AcoefV, BcoefU, BcoefV, & 
       pref, u1, v1, gustiness, rugoro, pctsrf, &
       lwdown_m, q2m, t2m, &
       snow, qsol, agesno, tsoil, &
       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, fluxbs, &   
       qsurf, tsurf_new, dflux_s, dflux_l, &
       flux_u1, flux_v1 , & 
       veget,lai,height &
#ifdef ISO
       ,xtprecip_rain, xtprecip_snow,xtspechum, &
       xtsnow, xtsol,xtevap,h1, &
       runoff_diag,xtrunoff_diag,Rland_ice &
#endif               
               )

    USE dimphy
    USE surface_data, ONLY    : ok_veget
    USE carbon_cycle_mod


    ! See comments in each module surf_land_orchidee_xxx for compatiblity with ORCHIDEE
#ifdef ORCHIDEE_NOOPENMP
    ! Compilation with cpp key ORCHIDEE NOOPENMP
    USE surf_land_orchidee_noopenmp_mod
#else
#if ORCHIDEE_NOZ0H
    ! Compilation with cpp key ORCHIDEE NOZ0H
    USE surf_land_orchidee_noz0h_mod
#else
#if ORCHIDEE_NOFREIN
    ! Compilation with cpp key ORCHIDEE_NOFREIN
    USE surf_land_orchidee_nofrein_mod
#else
#if ORCHIDEE_NOUNSTRUCT
    ! Compilation with cpp key ORCHIDEE_NOUNSTRUCT
    USE surf_land_orchidee_nounstruct_mod
#else
#if ORCHIDEE_NOLIC
    ! Compilation with cpp key ORCHIDEE_NOLIC
    USE surf_land_orchidee_nolic_mod
#else
    ! Default version
    USE surf_land_orchidee_mod
#endif
#endif
#endif
#endif
#endif
    
    USE surf_land_bucket_mod
    USE calcul_fluxs_mod
    USE indice_sol_mod
#ifdef ISO
    use infotrac_phy, ONLY: ntiso,niso
    use isotopes_mod, ONLY: nudge_qsol, iso_eau
#ifdef ISOVERIF
    use isotopes_verif_mod
#endif
#endif

    USE clesphys_mod_h
    USE yomcst_mod_h
USE print_control_mod, ONLY: lunout
    USE dimsoil_mod_h, ONLY: nsoilmx

    INCLUDE "dimpft.h"

! Input variables  
!****************************************************************************************
    INTEGER, INTENT(IN)                     :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
    REAL, INTENT(IN)                        :: date0
    REAL, DIMENSION(klon), INTENT(IN)       :: rlon, rlat
    REAL, DIMENSION(klon), INTENT(IN)       :: yrmu0  ! cosine of solar zenith angle
    LOGICAL, INTENT(IN)                     :: debut, lafin
    REAL, INTENT(IN)                        :: dtime
    REAL, DIMENSION(klon), INTENT(IN)       :: zlev, ccanopy
    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
    REAL, DIMENSION(klon), INTENT(IN)       :: albedo  ! albedo for whole short-wave interval
    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)       :: cdragh, cdragm
    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow, precip_bs
    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)       :: AcoefH, AcoefQ, BcoefH, BcoefQ
    REAL, DIMENSION(klon), INTENT(IN)       :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)       :: pref   ! pressure reference
    REAL, DIMENSION(klon), INTENT(IN)       :: u1, v1, gustiness
    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
    REAL, DIMENSION(klon), INTENT(IN)       :: lwdown_m  ! downwelling longwave radiation at mean surface
                                                         ! corresponds to previous sollwdown
    REAL, DIMENSION(klon), INTENT(IN)       :: q2m, t2m
#ifdef ISO
    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtspechum
#endif 
! In/Output variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
#ifdef ISO
    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsnow, xtsol
#endif

! Output variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
!albedo SB >>>
!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new ! albdeo for shortwave interval 1(visible)
!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new ! albedo for shortwave interval 2(near infrared)
    REAL, DIMENSION(6), INTENT(IN) :: SFRWL
    REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
!albedo SB <<<
    REAL, DIMENSION(klon), INTENT(OUT)       :: evap
    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat, fluxbs
    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1  ! flux for U and V at first model level
    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: veget,lai
    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: height
#ifdef ISO
    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
    REAL, DIMENSION(klon), INTENT(OUT)      :: h1
    REAL, DIMENSION(niso,klon), INTENT(OUT)      :: xtrunoff_diag
    REAL, DIMENSION(klon), INTENT(OUT)      :: runoff_diag
    REAL, DIMENSION(niso,klon), INTENT(IN)        :: Rland_ice
#endif

! Local variables
!****************************************************************************************
    REAL, DIMENSION(klon) :: p1lay_tmp
    REAL, DIMENSION(klon) :: pref_tmp
    REAL, DIMENSION(klon) :: swdown     ! downwelling shortwave radiation at land surface
    REAL, DIMENSION(klon) :: epot_air           ! potential air temperature
    REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
    REAL, DIMENSION(klon) :: u0, v0     ! surface speed
    REAL, DIMENSION(klon) :: precip_totsnow     ! total solid precip
    INTEGER               :: i

!albedo SB >>>
    REAL, DIMENSION(klon)      :: alb1_new,alb2_new
!albedo SB <<<

#ifdef ISO       
      real, parameter :: t_coup = 273.15
      real, dimension(klon) :: fqfonte_diag
      real, dimension(klon) :: snow_evap_diag 
      real, dimension(klon) :: fqcalving_diag 
      integer :: ixt
#endif 
!****************************************************************************************
!Total solid precip

IF (ok_bs) THEN
precip_totsnow(:)=precip_snow(:)+precip_bs(:)
ELSE
precip_totsnow(:)=precip_snow(:)
ENDIF
!****************************************************************************************
#ifdef ISO
#ifdef ISOVERIF
!        write(*,*) 'surf_land_mod 162'
        do i=1,knon
          if (iso_eau.gt.0) then
            call iso_verif_egalite_choix(precip_snow(i), &
     &          xtprecip_snow(iso_eau,i),'surf_land_mod 129', &
     &          errmax,errmaxrel)
            call iso_verif_egalite_choix(qsol(i), &
     &          xtsol(iso_eau,i),'surf_land_mod 139', &
     &          errmax,errmaxrel)
          endif  
        enddo
#endif 
#ifdef ISOVERIF
!       write(*,*) 'surf_land 169: ok_veget=',ok_veget
        do i=1,knon
         do ixt=1,ntiso
           call iso_verif_noNaN(xtprecip_snow(ixt,i),'surf_land 146')
         enddo
        enddo
#endif
#endif


!**************************************************************************************** 
! Choice between call to vegetation model (ok_veget=true) or simple calculation below
!
!****************************************************************************************
   IF (ok_veget) THEN
!****************************************************************************************
!  Call model sechiba in model ORCHIDEE
!
!****************************************************************************************
       p1lay_tmp(:)      = 0.0
       pref_tmp(:)       = 0.0
       p1lay_tmp(1:knon) = p1lay(1:knon)/100.
       pref_tmp(1:knon)  = pref(1:knon)/100.
! 
!* Calculate incoming flux for SW and LW interval: swdown
!
       swdown(:) = 0.0
       DO i = 1, knon
          swdown(i) = swnet(i)/(1-albedo(i))
       END DO
!
!* Calculate potential air temperature
!
       epot_air(:) = 0.0
       DO i = 1, knon
          epot_air(i) = RCPD*temp_air(i)*(pref(i)/p1lay(i))**RKAPPA
       END DO

#ifdef ISO
      CALL abort_physic('surf_land_mod 220','isos pas prevus dans orchidee',1)
#endif
       ! temporary for keeping same results using lwdown_m instead of lwdown
       CALL surf_land_orchidee(itime, dtime, date0, knon, &
            knindex, rlon, rlat, yrmu0, pctsrf, &
            debut, lafin, &
            zlev,  u1, v1, gustiness, temp_air, spechum, epot_air, ccanopy, & 
            cdragh, AcoefH, AcoefQ, BcoefH, BcoefQ, &
            precip_rain, precip_totsnow, lwdown_m, swnet, swdown, &
            pref_tmp, q2m, t2m, &
            evap, fluxsens, fluxlat,  &              
            tsol_rad, tsurf_new, alb1_new, alb2_new, &
            emis_new, z0m, z0h, qsurf, &
            veget, lai, height &
!#ifdef ISO
!            , xtprecip_rain, xtprecip_snow, xtspechum, xtevap &
!#endif
            )                  

#ifdef ISO
#ifdef ISOVERIF
     write(*,*) 'surf_land 193: apres surf_land_orchidee'   
     do i=1,knon
        if (iso_eau.gt.0) then
             call iso_verif_egalite_choix(xtevap(iso_eau,i),evap(i), &
    &            'surf_land 197',errmax,errmaxrel)
        endif !if (iso_eau.gt.0) then      
      enddo !do i=1,knon  
#endif
#endif
!  
!* Add contribution of relief to surface roughness
!  
       DO i=1,knon
          z0m(i) = MAX(1.5e-05,SQRT(z0m(i)**2 + rugoro(i)**2))
       ENDDO

    ELSE  ! not ok_veget
!****************************************************************************************
! No extern vegetation model choosen, call simple bucket calculations instead.
!
!****************************************************************************************
#ifdef ISO
#ifdef ISOVERIF
!       write(*,*) 'surf_land 247'
        call iso_verif_egalite_vect1D( &
     &           xtsnow,snow,'surf_land_mod 207',niso,klon)
#endif
#endif

#ifdef ISO
        if (nudge_qsol.eq.1) then
          call surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)
        endif
        !write(*,*) 'surf_land 258'
#endif
       CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
            tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
            spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, pref, &
            u1, v1, gustiness, rugoro, swnet, lwnet, &
            snow, qsol, agesno, tsoil, &
            qsurf, z0m, alb1_new, alb2_new, evap, &
            fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l &
#ifdef ISO
            ,xtprecip_rain, xtprecip_snow,xtspechum, &
            xtsnow, xtsol,xtevap,h1, &
     &      runoff_diag, xtrunoff_diag,Rland_ice &
#endif           
     &       )
        z0h(1:knon)=z0m(1:knon) ! En attendant mieux


    ENDIF ! ok_veget

        ! blowing snow not treated yet over land
        fluxbs(:)=0.


!****************************************************************************************
! Calculation for all land models
! - Flux calculation at first modele level for U and V
!****************************************************************************************
! Suppose zero surface speed
    u0(:)=0.0
    v0(:)=0.0
    CALL calcul_flux_wind(knon, dtime, &
         u0, v0, u1, v1, gustiness, cdragm, &
         AcoefU, AcoefV, BcoefU, BcoefV, &
         p1lay, temp_air, &
         flux_u1, flux_v1)

#ifdef ISO
#ifdef ISOVERIF
!     write(*,*) 'surf_land 237: sortie'   
      DO i=1,knon
        IF (iso_eau >= 0) THEN
             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
    &            'surf_land 241',errmax,errmaxrel)
        ENDIF !if (iso_eau.gt.0) then      
      ENDDO !do i=1,knon  
#endif
#endif

!albedo SB >>>
     SELECT CASE(NSW)
     CASE(2)
       alb_dir_new(1:knon,1)=alb1_new(1:knon)
       alb_dir_new(1:knon,2)=alb2_new(1:knon)
     CASE(4)
       alb_dir_new(1:knon,1)=alb1_new(1:knon)
       alb_dir_new(1:knon,2)=alb2_new(1:knon)
       alb_dir_new(1:knon,3)=alb2_new(1:knon)
       alb_dir_new(1:knon,4)=alb2_new(1:knon)
     CASE(6)
       alb_dir_new(1:knon,1)=alb1_new(1:knon)
       alb_dir_new(1:knon,2)=alb1_new(1:knon)
       alb_dir_new(1:knon,3)=alb1_new(1:knon)
       alb_dir_new(1:knon,4)=alb2_new(1:knon)
       alb_dir_new(1:knon,5)=alb2_new(1:knon)
       alb_dir_new(1:knon,6)=alb2_new(1:knon)
     END SELECT

     alb_dif_new=alb_dir_new
!albedo SB <<<
    
  END SUBROUTINE surf_land


#ifdef ISO
  SUBROUTINE surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)

    USE dimphy   
    USE infotrac_phy, ONLY: niso
    USE isotopes_mod, ONLY: region_nudge_qsol    
    INTEGER, INTENT(IN)                       :: knon         
    REAL, DIMENSION(klon), INTENT(IN)         :: rlon, rlat
    REAL, DIMENSION(klon), INTENT(INOUT)      :: qsol
    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex    
    REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsol 
    REAL :: lat_min_nudge_qsol,lat_max_nudge_qsol
    REAL :: lon_min_nudge_qsol,lon_max_nudge_qsol
    INTEGER :: i,ixt
    REAL :: qsol_new

    IF (region_nudge_qsol == 1) THEN
        ! Aamzonie du Sud
        lat_min_nudge_qsol=-15.0
        lat_max_nudge_qsol=-5.0
        lon_min_nudge_qsol=-70.0
        lon_max_nudge_qsol=-50.0
    ELSE IF (region_nudge_qsol == 2) THEN
        ! Aamzonie du Nord
        lat_min_nudge_qsol=-5.0
        lat_max_nudge_qsol=5.0
        lon_min_nudge_qsol=-70.0
        lon_max_nudge_qsol=-50.0
    ELSE
        WRITE(*,*) 'surf_land 298: cas pas prevu'
        WRITE(*,*) 'region_nudge_qsol=',region_nudge_qsol
        stop
    ENDIF

!    write(*,*) 'surf_land 314: knon=',knon
!    write(*,*) 'rlat=',rlat
!    write(*,*) 'rlon=',rlon
!    write(*,*) 'region_nudge_qsol=',region_nudge_qsol

    DO i=1,knon
      IF ((rlat(knindex(i)) >= lat_min_nudge_qsol).and. &
  &       (rlat(knindex(i)) <= lat_max_nudge_qsol).and. &
  &       (rlon(knindex(i)) >= lon_min_nudge_qsol).and. &
  &       (rlon(knindex(i)) <= lon_max_nudge_qsol)) THEN
!        write(*,*) 'surf_land 324: bon domaine: rlat,rlon,qsol=', &
!  &             rlat(knindex(i)),rlon(knindex(i)),qsol(knindex(i))
        qsol_new=qsol(i)
        IF (region_nudge_qsol == 1) THEN   
           qsol_new=max(qsol(i),50.0)   
        ELSE IF (region_nudge_qsol == 2) THEN      
           qsol_new=max(qsol(i),120.0)
        ELSE !if (region_nudge_qsol.eq.1) then
           WRITE(*,*) 'surf_land 317: cas pas prevu'
           WRITE(*,*) 'region_nudge_qsol=',region_nudge_qsol
           STOP
        ENDIF !if (region_nudge_qsol.eq.1) then
        IF (qsol(i) > 0.0) THEN
           DO ixt=1,niso
              xtsol(ixt,i)=xtsol(ixt,i)*qsol_new/qsol(i)
           ENDDO
        ELSE !IF (qsol(i) > 0.0) THEN
           DO ixt=1,niso
             xtsol(ixt,i)=0.0
           ENDDO
        ENDIF !IF (qsol(i) > 0.0) THEN
        qsol(i)=qsol_new
        WRITE(*,*) 'surf_land 346: qsol_new=',qsol(i)      
     ENDIF ! if ((rlat(i).ge.lat_min_nudge_qsol).and.
  ENDDO !DO i=1,knon

  END SUBROUTINE surf_land_nudge_qsol
#endif

!
!****************************************************************************************
!  
END MODULE surf_land_mod
!
!****************************************************************************************
!