!
MODULE surf_land_bucket_mod
!
! Surface land bucket module
!
! This module is used when no external land model is choosen.
!
  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, pref, &
       u1, v1, gustiness, rugoro, swnet, lwnet, &
       snow, qsol, agesno, tsoil, &
       qsurf, z0_new, 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           
            )
!$gpum horizontal knon klon
    USE limit_read_mod
    USE surface_data
    USE simplehydrol_mod
    USE calcul_fluxs_mod
    USE cpl_mod
    USE dimphy
    USE geometry_mod, ONLY: longitude,latitude 
    USE mod_grid_phy_lmdz
    USE mod_phys_lmdz_para
    USE indice_sol_mod
#ifdef ISO
    use infotrac_phy, ONLY: ntiso,niso
    USE isotopes_mod, ONLY: iso_eau, iso_HDO, iso_O18, iso_O17, &
        ridicule_qsol
    USE isotopes_routines_mod, ONLY: calcul_iso_surf_ter_vectall
#ifdef ISOVERIF
    USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_noNaN, &
        iso_verif_aberrant_o17,iso_verif_egalite_choix,iso_verif_egalite
#endif
#endif
    USE clesphys_mod_h
    USE yomcst_mod_h
    USE dimsoil_mod_h, ONLY: nsoilmx
    USE albsno_mod, ONLY :  albsno
    USE soil_mod, ONLY :  soil
    USE calbeta_mod, ONLY :  calbeta
    USE calbeta_clim_mod, ONLY :  calbeta_clim


!****************************************************************************************
! Bucket calculations for surface.
!

! Input variables  
!****************************************************************************************
    INTEGER, INTENT(IN)                     :: itime, jour, knon
    INTEGER, DIMENSION(knon), INTENT(IN)    :: knindex
    LOGICAL, INTENT(IN)                     :: debut
    REAL, INTENT(IN)                        :: dtime
    REAL, DIMENSION(knon), INTENT(IN)       :: tsurf
    REAL, DIMENSION(knon), INTENT(IN)       :: p1lay
    REAL, DIMENSION(knon), INTENT(IN)       :: tq_cdrag
    REAL, DIMENSION(knon), INTENT(IN)       :: precip_rain, precip_snow
    REAL, DIMENSION(knon), INTENT(IN)       :: temp_air, spechum
    REAL, DIMENSION(knon), INTENT(IN)       :: petAcoef, peqAcoef
    REAL, DIMENSION(knon), INTENT(IN)       :: petBcoef, peqBcoef
    REAL, DIMENSION(knon), INTENT(IN)       :: pref
    REAL, DIMENSION(knon), INTENT(IN)       :: u1, v1, gustiness
    REAL, DIMENSION(knon), INTENT(IN)       :: rugoro
    REAL, DIMENSION(knon), INTENT(IN)       :: swnet, lwnet
#ifdef ISO
    REAL, DIMENSION(ntiso,knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
    REAL, DIMENSION(ntiso,knon), INTENT(IN) :: xtspechum    
#endif

! In/Output variables
!****************************************************************************************
    REAL, DIMENSION(knon), INTENT(INOUT)          :: snow, qsol
    REAL, DIMENSION(knon), INTENT(INOUT)          :: agesno
    REAL, DIMENSION(knon, nsoilmx), INTENT(INOUT) :: tsoil
#ifdef ISO
    REAL, DIMENSION(niso,knon), INTENT(INOUT)       :: xtsnow,xtsol
#endif

! Output variables
!****************************************************************************************
    REAL, DIMENSION(knon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(knon), INTENT(OUT)       :: z0_new
    REAL, DIMENSION(knon), INTENT(OUT)       :: alb1_new, alb2_new
    REAL, DIMENSION(knon), INTENT(OUT)       :: evap, fluxsens, fluxlat
    REAL, DIMENSION(knon), INTENT(OUT)       :: tsurf_new
    REAL, DIMENSION(knon), INTENT(OUT)       :: dflux_s, dflux_l      
#ifdef ISO
    REAL, DIMENSION(ntiso,knon), INTENT(OUT) :: xtevap
    REAL, DIMENSION(knon),       INTENT(OUT) :: h1
    REAL, DIMENSION(niso,knon),  INTENT(OUT) :: xtrunoff_diag
    REAL, DIMENSION(knon),       INTENT(OUT) :: runoff_diag
    REAL, DIMENSION(niso,knon),  INTENT(IN)  :: Rland_ice
#endif

! Local variables
!****************************************************************************************
    REAL, DIMENSION(knon) :: soilcap, soilflux
    REAL, DIMENSION(knon) :: cal, beta, dif_grnd, capsol
    REAL, DIMENSION(knon) :: alb_neig, alb_lim, icesub, icemelt
    REAL, DIMENSION(knon) :: zfra
    REAL, DIMENSION(knon) :: radsol       ! total net radiance at surface
    REAL, DIMENSION(knon) :: u0, v0, u1_lay, v1_lay
    REAL, DIMENSION(knon) :: dummy_riverflow,dummy_coastalflow 
    INTEGER               :: i
#ifdef ISO
    INTEGER               :: ixt
    REAL, DIMENSION(niso,knon) :: xtsnow_prec,xtsol_prec
    REAL, DIMENSION(knon) :: snow_prec,qsol_prec
    REAL, PARAMETER       :: t_coup = 273.15
    REAL, DIMENSION(knon) :: fq_fonte_diag
    REAL, DIMENSION(knon) :: fqfonte_diag
    REAL, DIMENSION(knon) :: snow_evap_diag 
    REAL, DIMENSION(knon) :: fqcalving_diag 
    REAL                  :: max_eau_sol_diag  
    REAL, DIMENSION(knon) :: run_off_lic_diag 
    REAL :: coeff_rel_diag
#endif 
    REAL :: lat(knon), lon(knon)

!
!****************************************************************************************

#ifdef ISO
#ifdef ISOVERIF
        !write(*,*) 'surf_land_bucket 152'
        DO i=1,knon
          IF (iso_eau > 0) THEN
            CALL iso_verif_egalite_choix(precip_snow(i), &
     &                                   xtprecip_snow(iso_eau,i),'surf_land_bucket 131', &
     &                                   errmax,errmaxrel)
            CALL iso_verif_egalite_choix(qsol(i), &
     &                                   xtsol(iso_eau,i),'surf_land_bucket 134', &
     &                                   errmax,errmaxrel)
          ENDIF
        ENDDO
#endif 
#ifdef ISOVERIF
        DO i=1,knon
          DO ixt=1,niso
            CALL iso_verif_noNaN(xtsol(ixt,i),'surf_land_mod_bucket 142')
          ENDDO !do ixt=1,niso
        ENDDO !do i=1,knon
        !write(*,*) 'surf_land_bucket 152'
#endif
#endif

!
!* Read from limit.nc : albedo(alb_lim), length of rugosity(z0_new)
!
    CALL limit_read_rug_alb(itime, dtime, jour,&
         knon, knindex, &
         z0_new, alb_lim)
!
!* Calcultaion of fluxes 
!

! calculate total absorbed radiance at surface
       radsol(:) = 0.0
       radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)

! calculate constants
    CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
    if (type_veget=='betaclim') then
       lat(1:knon) = latitude(knindex(1:knon))
       CALL calbeta_clim(knon,jour,lat ,beta)
    endif
       
! calculate temperature, heat capacity and conduction flux in soil
    IF (soil_model) THEN
       lat(1:knon) = latitude(knindex(1:knon))
       lon(1:knon) = longitude(knindex(1:knon))
       CALL soil(dtime, is_ter, knon, snow, tsurf, qsol,  & 
        & lon, lat, tsoil, soilcap, soilflux)

       DO i=1, knon
          cal(i) = RCPD / soilcap(i)
          radsol(i) = radsol(i)  + soilflux(i)
       END DO
    ELSE 
       cal(:) = RCPD * capsol(:)
       IF (klon_glo .EQ. 1) THEN
         cal(:) = 0.
       ENDIF
    ENDIF
    
! Suppose zero surface speed
    u0(:)=0.0
    v0(:)=0.0
    u1_lay(:) = u1(:) - u0(:)
    v1_lay(:) = v1(:) - v0(:)

    CALL calcul_fluxs(knon, is_ter, dtime, &
         tsurf, p1lay, cal, beta, tq_cdrag, tq_cdrag, pref, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
         1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    
#ifdef ISO
   ! verif
#ifdef ISOVERIF
    !write(*,*) 'surf_land_bucket 211'
    DO i=1,knon
      IF (iso_eau > 0) THEN
        CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), &
     &           snow(i),'surf_land_bucket 522', &
     &           errmax,errmaxrel) 
      ENDIF !IF (iso_eau > 0) then
    ENDDO !DO i=1,knon 
#endif
   ! end verif
#endif         
#ifdef ISO
    DO i=1,knon
      snow_prec(i)=snow(i)
      qsol_prec(i)=qsol(i)
      DO ixt=1,niso
        xtsnow_prec(ixt,i)=xtsnow(ixt,i)
        xtsol_prec(ixt,i) =xtsol(ixt,i)
      ENDDO !DO ixt=1,niso
      ! initialisation:
      fqfonte_diag(i)  =0.0
      fq_fonte_diag(i) =0.0
      snow_evap_diag(i)=0.0
    ENDDO !DO i=1,knon 
#ifdef ISOVERIF
    ! write(*,*) 'surf_land_bucket 235'
    DO i=1,knon  
      IF (iso_eau > 0) THEN
        CALL iso_verif_egalite(qsol_prec(i),xtsol_prec(iso_eau,i), &
    &                              'surf_land_bucket 141')
      ENDIF
    ENDDO !DO i=1,knon
#endif    
#endif    
!
!* Calculate snow height, run_off, age of snow
!      
    CALL simplehydrol( knon, is_ter, knindex, dtime, &
         tsurf, precip_rain, precip_snow, &
         snow, qsol, tsurf_new, evap, icesub, icemelt &
#ifdef ISO    
     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
     & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
#endif
     &   )

#ifdef ISO
#ifdef ISOVERIF
        DO i=1,knon
          DO ixt=1,niso
            CALL iso_verif_noNaN(xtsol_prec(ixt,i),'surf_land_burcket 237')
          ENDDO
        ENDDO
#endif
#ifdef ISOVERIF
        !write(*,*) 'surf_land_bucket 235'
        DO i=1,knon
          IF (iso_eau > 0) THEN
            CALL iso_verif_egalite_choix(qsol_prec(i), &
     &                                   xtsol_prec(iso_eau,i),'surf_land_bucket 628', &
     &                                   errmax,errmaxrel)
            CALL iso_verif_egalite_choix(precip_snow(i), &
     &                                   xtprecip_snow(iso_eau,i),'surf_land_bucket 227', &
     &                                   errmax,errmaxrel)
             ! attention, dans simplehydrol, on modifie snow sans modifier
             ! xtsnow
             ! c'est fait plus tard dans gestion_neige
!            write(*,*) 'surf_land_bucket 287: i=',i
!            write(*,*) 'snow(i)=',snow(i)
            CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), &
     &                                   snow_prec(i),'surf_land_bucket 245', &
     &                                   errmax,errmaxrel) 
          ENDIF  
          IF ((iso_O17 > 0).AND.(iso_O18 > 0)) THEN
              IF (qsol_prec(i) > ridicule_qsol) THEN
                CALL iso_verif_aberrant_o17(xtsol_prec(iso_O17,i)/qsol_prec(i) &
     &                                     ,xtsol_prec(iso_O18,i)/qsol_prec(i) &
     &                                     ,'surf_land_bucket 642')
              ENDIF !IF ((qsol_prec(i) > ridicule_qsol) &
          ENDIF !IF ((iso_O17 > 0).AND.(iso_O18 > 0)) THEN
        ENDDO  !DO i=1,knon 
        !write(*,*) 'surf_land_mod 291'
        !write(*,*) 'snow_evap_diag(1)=',snow_evap_diag(1)
#endif          
        CALL calcul_iso_surf_ter_vectall(knon,knon, &
     &           evap,snow_evap_diag,snow, &
     &           fq_fonte_diag,fqfonte_diag,dtime,precip_rain,xtprecip_rain, &
     &           precip_snow,xtprecip_snow, snow_prec,xtsnow_prec, &
     &           tsurf_new,xtspechum,pref,spechum,t_coup,u1_lay,v1_lay,p1lay, &
     &           qsol,xtsol,qsol_prec,xtsol_prec, &
     &           max_eau_sol_diag, &
     &           xtevap,xtsnow,h1,runoff_diag,xtrunoff_diag,fqcalving_diag, &
     &           knindex,is_ter,run_off_lic_diag,coeff_rel_diag,Rland_ice &
     &   )
!#ifdef ISOVERIF
!        write(*,*) 'surf_land_bucket 303'
!#endif
#endif

!
!* Calculate the age of snow
!
    CALL albsno(knon,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_lim(i)  = alb_neig(i) *zfra(i) + alb_lim(i)*(1.0-zfra(i))
    END DO

!
!* Return albedo : 
!    alb1_new and alb2_new are here given the same values
!
    alb1_new(:) = 0.0
    alb2_new(:) = 0.0
    alb1_new(1:knon) = alb_lim(1:knon)
    alb2_new(1:knon) = alb_lim(1:knon)
       
!
!* Calculate the rugosity
!
    DO i = 1, knon
       z0_new(i) = MAX(1.5e-05,SQRT(z0_new(i)**2 + rugoro(i)**2))
    END DO

!* Send to coupler
!  The run-off from river and coast are not calculated in the bucket modele.
!  For testing purpose of the coupled modele we put the run-off to zero.
    IF (type_ocean=='couple') THEN
       dummy_riverflow(:)   = 0.0
       dummy_coastalflow(:) = 0.0
!$gpum nocall
       CALL cpl_send_land_fields(itime, knon, knindex, &
            dummy_riverflow, dummy_coastalflow)
    ENDIF

!
!* End
!
  END SUBROUTINE surf_land_bucket
!
!****************************************************************************************
!
END MODULE surf_land_bucket_mod
