!
! $Header$
!
MODULE surf_land_bucket_mod
!
! Surface land bucket module
!
! This module is used when no external land model is choosen.
!
  USE fonte_neige_mod
  USE calcul_fluxs_mod
  USE dimphy
  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_para
 
  IMPLICIT NONE

CONTAINS

  SUBROUTINE surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
       spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, ps, &
       u1_lay, v1_lay, rugoro, &
       radsol, snow, qsol, agesno, tsoil, &
       qsurf, z0_new, alblw, evap, fluxsens, fluxlat, &
       tsurf_new, alb_new, dflux_s, dflux_l)

!****************************************************************************************
! Bucket calculations for surface. 
!
    INCLUDE "clesphys.h"
    INCLUDE "indicesol.h"
    INCLUDE "dimsoil.h"
    INCLUDE "YOMCST.h"

! Input variables  
!****************************************************************************************
    INTEGER, INTENT(IN)                     :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
    LOGICAL, INTENT(IN)                     :: debut
    REAL, INTENT(IN)                        :: dtime
    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)       :: tq_cdrag
    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)       :: petAcoef, peqAcoef
    REAL, DIMENSION(klon), INTENT(IN)       :: petBcoef, peqBcoef
    REAL, DIMENSION(klon), INTENT(IN)       :: ps
    REAL, DIMENSION(klon), INTENT(IN)       :: u1_lay, v1_lay
    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro

! In/Output variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil

! Output variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      

! Local variables
!****************************************************************************************
    REAL, DIMENSION(klon) :: soilcap, soilflux
    REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
    REAL, DIMENSION(klon) :: alb_neig
    REAL, DIMENSION(klon) :: zfra
    INTEGER               :: i
!
!****************************************************************************************


!
!* Read from limit.nc : albedo(alb_new), length of rugosity(z0_new)
!
    CALL interfsur_lim(itime, dtime, jour, &
         knon, knindex, debut,  &
         alb_new, z0_new)
    
!
!* Calcultaion of fluxes 
!
! calculate constants
    CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
       
! calculate temperature, heat capacity and conduction flux in soil
    IF (soil_model) THEN 
       CALL soil(dtime, is_ter, knon, snow, tsurf, tsoil, soilcap, soilflux)
       DO i=1, knon
          cal(i) = RCPD / soilcap(i)
          radsol(i) = radsol(i)  + soilflux(i)
       END DO
    ELSE 
       cal(:) = RCPD * capsol(:)
    ENDIF
    
    CALL calcul_fluxs(knon, is_ter, dtime, &
         tsurf, p1lay, cal, beta, tq_cdrag, ps, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
         petAcoef, peqAcoef, petBcoef, peqBcoef, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    
!
!* Calculate snow height, run_off, age of snow
!      
    CALL fonte_neige( knon, is_ter, knindex, dtime, &
         tsurf, precip_rain, precip_snow, &
         snow, qsol, tsurf_new, evap)
!
!* Calculate the age of snow
!
    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
    
    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
    
    DO i=1, knon
       zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
       alb_new(i)  = alb_neig(i) *zfra(i) + alb_new(i)*(1.0-zfra(i))
    END DO

    alblw(:) = 0.0
    alblw(1:knon) = alb_new(1:knon)

!
!* Calculate the rugosity
!
    z0_new = SQRT(z0_new**2+rugoro**2)
       
!
!* End
!
  END SUBROUTINE surf_land_bucket
!
!****************************************************************************************
!
  SUBROUTINE interfsur_lim(itime, dtime, jour, &
       knon, knindex, debut, &
       lmt_alb_p, lmt_rug_p)
    
! Cette routine sert d'interface entre le modele atmospherique et un fichier
! de conditions aux limites
!
! L. Fairhead 02/2000
!
! input:
!   itime        numero du pas de temps courant
!   dtime        pas de temps de la physique (en s)
!   jour         jour a lire dans l'annee
!   knon         nombre de points dans le domaine a traiter
!   knindex      index des points de la surface a traiter
!   debut        logical: 1er appel a la physique (initialisation)
!
! output:
!   lmt_alb_p      Albedo lu 
!   lmt_rug_p      longueur de rugosite lue

    INCLUDE "netcdf.inc"

! Input variables
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: itime
    REAL   , INTENT(IN)                      :: dtime
    INTEGER, INTENT(IN)                      :: jour
    INTEGER, INTENT(IN)                      :: knon
    INTEGER, DIMENSION(klon_loc), INTENT(IN) :: knindex
    LOGICAL, INTENT(IN)                      :: debut

! Output variables
!****************************************************************************************
    REAL, INTENT(out), DIMENSION(klon_loc)   :: lmt_alb_p
    REAL, INTENT(out), DIMENSION(klon_loc)   :: lmt_rug_p

! Local variables with attribute SAVE
!****************************************************************************************
    INTEGER,SAVE   :: lmt_pas     ! frequence de lecture des conditions limites 
                                  ! (en pas de physique)
    !$OMP THREADPRIVATE(lmt_pas) 
    LOGICAL,SAVE   :: deja_lu_sur ! pour indiquer que le jour a lire a deja
                                  ! lu pour une surface precedente
    !$OMP THREADPRIVATE(deja_lu_sur) 
    INTEGER,SAVE                           :: jour_lu_sur 
    !$OMP THREADPRIVATE(jour_lu_sur) 
    CHARACTER (len = 20),SAVE              :: fich ='limit.nc'
    !$OMP THREADPRIVATE(fich) 
    LOGICAL,SAVE                           :: check = .FALSE.
    !$OMP THREADPRIVATE(check) 
! Champs lus dans le fichier de CL
    REAL, ALLOCATABLE , SAVE, DIMENSION(:) :: alb_lu_p, rug_lu_p
    !$OMP THREADPRIVATE(alb_lu_p, rug_lu_p) 

! quelques variables pour netcdf
    INTEGER ,SAVE                          :: nid, nvarid
    !$OMP THREADPRIVATE(nid, nvarid) 
    INTEGER, DIMENSION(2),SAVE             :: start, epais
    !$OMP THREADPRIVATE(start, epais) 

! Other local variables
!****************************************************************************************
    INTEGER                                :: ii, ierr
    CHARACTER (len = 20)                   :: modname = 'interfsur_lim'
    CHARACTER (len = 80)                   :: abort_message
    REAL, DIMENSION(klon_glo)              :: alb_lu
    REAL, DIMENSION(klon_glo)              :: rug_lu

!
! End delaration
!****************************************************************************************

    IF (debut) THEN
       lmt_pas = NINT(86400./dtime * 1.0) ! pour une lecture une fois par jour
       jour_lu_sur = jour - 1
       ALLOCATE(alb_lu_p(klon_loc))
       ALLOCATE(rug_lu_p(klon_loc))
    ENDIF
    
    IF ((jour - jour_lu_sur) /= 0) deja_lu_sur = .FALSE.
  
    IF (check) WRITE(*,*) modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur 
    IF (check) WRITE(*,*) modname,':: itime, lmt_pas', itime, lmt_pas
    IF (check) CALL flush(6)

!    
! Tester d'abord si c'est le moment de lire le fichier
!
    IF (MOD(itime-1, lmt_pas) == 0 .AND. .NOT. deja_lu_sur) THEN

!
! Ouverture et lecture du fichier
!
!$OMP MASTER
       IF (is_mpi_root) THEN
          fich = TRIM(fich)
          IF (check) WRITE(*,*)modname,' ouverture fichier ',fich
          IF (check) CALL flush(6)
          ierr = NF_OPEN (fich, NF_NOWRITE,nid)
          IF (ierr.NE.NF_NOERR) THEN
             abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF
! 
! La tranche de donnees a lire:
          start(1) = 1
          start(2) = jour
          epais(1) = klon_glo
          epais(2) = 1
!
! Lecture albedo
          ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
          IF (ierr /= NF_NOERR) THEN
             abort_message = 'Le champ <ALB> est absent'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF
#ifdef NC_DOUBLE
          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
#else
          ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
#endif
          IF (ierr /= NF_NOERR) THEN
             abort_message = 'Lecture echouee pour <ALB>'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF
!
! Lecture rugosite!
          ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
          IF (ierr /= NF_NOERR) THEN
             abort_message = 'Le champ <RUG> est absent'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF
#ifdef NC_DOUBLE
          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
#else
          ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
#endif
          IF (ierr /= NF_NOERR) THEN
             abort_message = 'Lecture echouee pour <RUG>'
             CALL abort_gcm(modname,abort_message,1)
          ENDIF

!
! Fin de lecture
          ierr = NF_CLOSE(nid)

       ENDIF ! is_mpi_root
!$OMP END MASTER

       CALL Scatter(alb_lu,alb_lu_p)
       CALL Scatter(rug_lu,rug_lu_p)

       deja_lu_sur = .TRUE.
       jour_lu_sur = jour       

    ENDIF
  
!
! Recopie des variables dans les champs de sortie
!
    lmt_alb_p(:) = 999999.
    lmt_rug_p(:) = 999999.
    DO ii = 1, knon
       lmt_alb_p(ii) = alb_lu_p(knindex(ii))
       lmt_rug_p(ii) = rug_lu_p(knindex(ii))
    ENDDO
    

  END SUBROUTINE interfsur_lim
!
!****************************************************************************************
!
END MODULE surf_land_bucket_mod