!
MODULE surf_land_bucket_hetero_mod
!
! 2025/04 A. Maison (adapted from surf_land_bucket_mod)
! Surface land bucket module with heterogeneous continental sub-surfaces
! This module is used when no external land model is choosen and iflag_hetero_surf = 1 or 2.
!
  IMPLICIT NONE

CONTAINS

  SUBROUTINE surf_land_bucket_hetero(itime, jour, knon, knindex, debut, dtime,&
       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
       spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, pref, plev, &
       u1, v1, gustiness, rugoro, swnet, lwnet, &
       snow, qsol, agesno, tsoil, &
       qsurf, z0m, z0h, alb1_new, alb2_new, evap, &
       fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l, &
       tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
       cdragm_tersrf, cdragh_tersrf, &
       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)

    USE clesphys_mod_h
    USE dimsoil_mod_h, ONLY: nsoilmx
    USE yomcst_mod_h, ONLY: RD, RG, RCPD, RSIGMA
    USE compbl_mod_h
    USE dimpft_mod_h
    USE limit_read_mod
    USE surface_data
    USE fonte_neige_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
    USE phys_state_var_mod, ONLY: frac_tersrf, ratio_z0m_z0h_tersrf, z0m_tersrf, &
                                  albedo_tersrf, beta_tersrf, inertie_tersrf, &
                                  hcond_tersrf
    USE surf_param_mod, ONLY: eff_surf_param, average_surf_var
    USE cdrag_mod

#ifdef ISO
    use infotrac_phy, ONLY: niso
#endif

!****************************************************************************************
! Bucket calculations for surface. 
!****************************************************************************************
!
! 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)       :: pref
    REAL, DIMENSION(klon), INTENT(IN)       :: u1, v1, gustiness
    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
    REAL, DIMENSION(klon, nbtersrf), INTENT(IN) :: tsurf_tersrf


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


! Output variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new, alb2_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
    !
    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: tsurf_new_tersrf
    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: qsurf_tersrf
    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragm_tersrf, cdragh_tersrf
    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: swnet_tersrf, lwnet_tersrf
    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxsens_tersrf, fluxlat_tersrf

#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)  :: runoff_diag
!    REAL, DIMENSION(niso,klon),  INTENT(IN)  :: Rland_ice
#endif
! Local variables
!****************************************************************************************
    REAL, DIMENSION(klon) :: soilcap, soilflux
    REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
    REAL, DIMENSION(klon) :: alb_neig, alb_lim, icesub
    REAL, DIMENSION(klon) :: zfra
    REAL, DIMENSION(klon) :: radsol
    REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay
    REAL, DIMENSION(klon) :: dummy_riverflow,dummy_coastalflow 
    INTEGER               :: i, j, k
#ifdef ISO
    INTEGER               :: ixt
    REAL, PARAMETER       :: t_coup = 273.15
    REAL, DIMENSION(klon) :: fq_fonte_diag
    REAL, DIMENSION(klon) :: fqfonte_diag
    REAL, DIMENSION(klon) :: snow_evap_diag 
    REAL, DIMENSION(klon) :: fqcalving_diag 
    REAL                  :: max_eau_sol_diag  
    REAL, DIMENSION(klon) :: run_off_lic_diag 
    REAL :: coeff_rel_diag
#endif
    !
    REAL, DIMENSION(klon)           :: zlev, geop1, speed, pblh, ri_in, sst
    REAL, DIMENSION(klon)           :: beta_eff, inertie_eff, conv_ratio_eff
    REAL, DIMENSION(klon)           :: meansqT
    REAL, DIMENSION(klon, nbtersrf) :: z0h_tersrf, emis_tersrf, conv_ratio_tersrf
    REAL, DIMENSION(klon, nbtersrf) :: evap_tersrf
    REAL, DIMENSION(klon, nbtersrf) :: dflux_s_tersrf, dflux_l_tersrf
    REAL, DIMENSION(klon, nbtersrf) :: radsol_tersrf
    REAL, DIMENSION(klon, nbtersrf) :: zri_tersrf
    REAL, PARAMETER                 :: klon_1D = 1

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

    ! *** Calculations common to the two flag values ***

    ! average albedo
    alb_lim = eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI')

    ! suppose zero surface speed
    u0(:)=0.0
    v0(:)=0.0
    u1_lay(:) = u1(:) - u0(:)
    v1_lay(:) = v1(:) - v0(:)
    speed(:) = (u1_lay(:)**2 + v1_lay(:)**2)**0.5
    !
    geop1(1:knon) = RD * temp_air(1:knon) / (0.5*(pref(1:knon)+p1lay(1:knon))) &
            * (pref(1:knon)-p1lay(1:knon))
    zlev(1:knon) = (plev(1:knon)*RD*temp_air(1:knon)/(pref(1:knon)*RG))/2.
    !
    ! compute roughness lengths
    DO i=1, knon
      DO j=1, nbtersrf
        IF (ratio_z0m_z0h_tersrf(i,j) .NE. 0.) THEN
          z0h_tersrf(i,j) = z0m_tersrf(i,j) / ratio_z0m_z0h_tersrf(i,j)
        ELSE
          z0h_tersrf(i,j) = 1.E-12
        ENDIF
      ENDDO
    ENDDO

    z0m = eff_surf_param(klon, nbtersrf, z0m_tersrf, frac_tersrf, 'CDN', zlev)
    z0h = eff_surf_param(klon, nbtersrf, z0h_tersrf, frac_tersrf, 'CDN', zlev)

    DO i=1, knon
       z0m(i) = MAX(1.5e-05,SQRT(z0m(i)**2 + rugoro(i)**2))
    END DO

    ! compute the ratio to convert and print soil depths in meters (conv_ratio = (cond/cap)^0.5 and cap = I^2/cond)
    DO j=1, nbtersrf
       conv_ratio_tersrf(:,j) = hcond_tersrf(:,j)/inertie_tersrf(:,j)
    ENDDO

    !
    ! *** Surface parameter aggregation ***
    !
    IF (iflag_hetero_surf == 1) THEN
      !* Calcultaion of fluxes 
      
      ! calculate total absorbed radiance at surface
      radsol(:) = 0.0
      radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)

      ! calculate constants (needeed for capsol and dif_grnd)
      CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
      IF (type_veget=='betaclim') THEN
        CALL calbeta_clim(knon,jour,latitude(knindex(1:knon)),beta)
      ENDIF

      ! mean evapotranspiration coefficient
      beta_eff = eff_surf_param(klon, nbtersrf, beta_tersrf, frac_tersrf, 'ARI')
      beta = beta_eff

      ! calculate temperature, heat capacity and conduction flux in soil
      IF (soil_model) THEN
         inertie_eff = eff_surf_param(klon, nbtersrf, inertie_tersrf, frac_tersrf, 'ARI')
         conv_ratio_eff = eff_surf_param(klon, nbtersrf, conv_ratio_tersrf, frac_tersrf, 'ARI')
         !
         CALL soil_hetero(dtime, is_ter, knon, snow, tsurf, qsol,  & 
           & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux, &
           & inertie_eff, conv_ratio_eff)
         !
         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

      ! calculate fluxes
      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)
    
      
      !* Calculate snow height, run_off, age of snow
#ifdef ISO
    DO i=1,knon
      ! initialisation:
      fqfonte_diag(i)  =0.0
      fq_fonte_diag(i) =0.0
      snow_evap_diag(i)=0.0
    ENDDO !DO i=1,knon 
#endif    

      CALL fonte_neige( knon, is_ter, knindex, dtime, &
           tsurf, precip_rain, precip_snow, &
           snow, qsol, tsurf_new, evap, icesub &
#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
     &  )
    
      ! 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_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)
       
      !* 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
         CALL cpl_send_land_fields(itime, knon, knindex, &
              dummy_riverflow, dummy_coastalflow)
      ENDIF

    !
    ! ***  Flux aggregation ***
    !
    ELSE IF (iflag_hetero_surf == 2) THEN
      ! initialize output tables
      evap_tersrf(:,:)      = 0.
      fluxsens_tersrf(:,:)  = 0.
      fluxlat_tersrf(:,:)   = 0.
      tsurf_new_tersrf(:,:) = 0.
      dflux_s_tersrf(:,:)   = 0.
      dflux_l_tersrf(:,:)   = 0.
      radsol_tersrf(:,:)    = 0.
      swnet_tersrf(:,:)     = 0.
      lwnet_tersrf(:,:)     = 0.
      ! hyp: surface emissivity = 1
      emis_tersrf(:,:)      = 1.
 
      ! * calculate total absorbed radiance at surface
      DO j=1, nbtersrf
        ! SW
        swnet_tersrf(klon_1D,j) = (1. - albedo_tersrf(klon_1D,j)) / (1. - alb_lim(klon_1D)) * swnet(klon_1D)
        ! LW
        ! first order
        lwnet_tersrf(klon_1D,j) = lwnet(klon_1D) + 4. * emis_tersrf(klon_1D,j) * RSIGMA * tsurf(klon_1D)**3 * &
                                  (tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))
      ENDDO
      ! LW second order corrections
      !- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
      IF (iflag_order2_sollw == 1) THEN
        meansqT(:) = 0. ! as working buffer
        !
        DO j=1, nbtersrf
          meansqT(klon_1D) = meansqT(klon_1D) + (tsurf_tersrf(klon_1D,j) - tsurf(klon_1D))**2 * frac_tersrf(klon_1D,j)
        ENDDO
        DO j=1, nbtersrf
          lwnet_tersrf(klon_1D,j) = lwnet_tersrf(klon_1D,j) + 6. * RSIGMA * tsurf(klon_1D)**2 * (meansqT(klon_1D) - &
                              (tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))**2)
        ENDDO
      ENDIF
      ! net radiation
      radsol_tersrf(:,:) = swnet_tersrf(:,:) + lwnet_tersrf(:,:)
 
      ! * compute evapotranspiration coefficient
      capsol(:) = 1.0/(2.5578E+06*0.15)
      dif_grnd(:) = 0.

      ! unused variables in cdrag routine
      pblh(:) = 0.
      ri_in(:) = 0.
      sst(:) = 0.

      ! Loop on sub-surfaces
      DO j=1, nbtersrf
        ! * drag coefficients
        CALL cdrag(knon, is_ter, speed, temp_air, spechum, geop1, pref, pblh, &
          tsurf_tersrf(:,j), qsurf_tersrf(:,j), z0m_tersrf(:,j), z0h_tersrf(:,j), ri_in, 0, &
          cdragm_tersrf(:,j), cdragh_tersrf(:,j), zri_tersrf(:,j), plev, precip_rain, sst, p1lay)

        ! * calculate temperature, heat capacity and conduction flux in soil
        IF (soil_model) THEN
          !
          CALL soil_hetero(dtime, is_ter, knon, snow, tsurf_tersrf(:,j), qsol,  &
            longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil_tersrf(:,:,j), soilcap, soilflux, &
            inertie_tersrf(:,j), conv_ratio_tersrf(:,j))
          !
          cal(:) = RCPD / soilcap(:)
          radsol_tersrf(:,j) = radsol_tersrf(:,j)  + soilflux(:)
          !
        ELSE
          cal = RCPD * capsol
          IF (klon_glo .EQ. 1) THEN
            cal = 0.
          ENDIF
        ENDIF

        ! * calcultaion of fluxes 
        CALL calcul_fluxs(knon, is_ter, dtime, &
             tsurf_tersrf(klon_1D,j), p1lay, cal, beta_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), pref, &
             precip_rain, precip_snow, snow, qsurf_tersrf(klon_1D,j),  &
             radsol_tersrf(klon_1D,j), dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
             1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
             tsurf_new_tersrf(klon_1D,j), evap_tersrf(klon_1D,j), fluxlat_tersrf(klon_1D,j), fluxsens_tersrf(klon_1D,j), &
             dflux_s_tersrf(klon_1D,j), dflux_l_tersrf(klon_1D,j))

        ! if snow > 0
        ! calculate snow height, run_off, age of snow
        CALL fonte_neige( knon, is_ter, knindex, dtime, &
             tsurf_tersrf(:,j), precip_rain, precip_snow, &
             snow, qsol, tsurf_new_tersrf(:,j), evap_tersrf(:,j), icesub &
#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
     &  )

      ENDDO ! loop on sub-surfaces

      ! 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_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)
       
      ! 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
         CALL cpl_send_land_fields(itime, knon, knindex, &
              dummy_riverflow, dummy_coastalflow)
      ENDIF

      ! * average of fluxes and surface variables
      qsurf = average_surf_var(klon, nbtersrf, qsurf_tersrf, frac_tersrf, 'ARI')
      tsurf_new = average_surf_var(klon, nbtersrf, tsurf_new_tersrf, frac_tersrf, 'ARI')
      evap = average_surf_var(klon, nbtersrf, evap_tersrf, frac_tersrf, 'ARI')
      fluxlat = average_surf_var(klon, nbtersrf, fluxlat_tersrf, frac_tersrf, 'ARI')
      fluxsens = average_surf_var(klon, nbtersrf, fluxsens_tersrf, frac_tersrf, 'ARI')
      dflux_l = average_surf_var(klon, nbtersrf, dflux_l_tersrf, frac_tersrf, 'ARI')
      dflux_s = average_surf_var(klon, nbtersrf, dflux_s_tersrf, frac_tersrf, 'ARI')
      DO k=1, nsoilmx
        tsoil(:,k) = average_surf_var(klon, nbtersrf, tsoil_tersrf(:,k,:), frac_tersrf, 'ARI')
      ENDDO

      ! order 2 correction to tsurf_new, for radiation computations (main atm effect of Ts)
      IF (iflag_order2_sollw == 1) THEN
        meansqT(:) = 0. ! as working buffer
        DO j=1, nbtersrf
          meansqT(klon_1D) = meansqT(klon_1D)+(tsurf_tersrf(klon_1D,j)-tsurf_new(klon_1D))**2 *frac_tersrf(klon_1D,j)
        ENDDO
        tsurf_new(:) = tsurf_new(:)+1.5*meansqT(:)/tsurf_new(:)
      ENDIF   ! iflag_order2_sollw == 1

    ENDIF ! iflag_hetero_surf
!
!* End
!
  END SUBROUTINE surf_land_bucket_hetero
!
!****************************************************************************************
!
END MODULE surf_land_bucket_hetero_mod
