mode_tartes.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
MODULE mode_tartes

!##########################
!
!! *MODE_TARTES*
!!
!! Radiative transfer in snowpack

!!
!!**  IMPLICIT ARGUMENTS
!!    ------------------
!!       NONE          
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!    M. Lafaysse       * Meteo France *
!!    translated from python codes of G. Picard, Q. Libois, LGGE.
!!
!! Main differences relatively to the python code :
!! ------------------------------------------------
!!
!! For optimization on large domains :
!!     * Loops are inside the subroutines 
!!     * All variables are multi-points (first dimension)
!!     * Loop on points is the last loop
!!     * Number of active or effective layers : argument of subroutines
!!
!! New routine to interpolate ice refractive index on the given wavelengths
!!
!! Wavelengths are assumed to be a parameter of the model and not an argument (if necessary, it will be in surfex namelists)
!!
!! Impurities : the python object is replaced by two variables (density and content). Last dimension represents impurity type.
!! For now : it is just implemented for soot (indice 1).
!!
!! To solve the linear system MX=Y :
!!     * Indices of upper diagonal (PDP) shifted one step right
!!     * Indices of lower diagonal (PDM) shifted one step left
!!
!! In routine INFINITE_MEDIUM_OPTICAL_PARAMETERS add a control IF ABS(PSNOWALB)<UEPSI
!!
!! Add a control of absorbed energy output to avoid occasional numerical problem in infra-red (not used)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    24/07/2013
!!      Matthieu Lafaysse interface with SURFEX 23/08/2013
!--------------------------------------------------------------------------------

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb


 CONTAINS

SUBROUTINE tartes(PSNOWSSA,PSNOWRHO,PSNOWDZ,PSNOWG0,PSNOWY0,PSNOWW0,PSNOWB0,PSNOWIMP_DENSITY,&
                  psnowimp_content,palb,psw_rad_dif,psw_rad_dir,pcoszen,knlvls_use,psnowalb, &
                  psnowenergy,psoilenergy)
!
USE modd_const_tartes, ONLY: npnbands,xpwavelengths,xrefice_r,xrefice_i,xrefimp_i,xp_mudiff
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowssa !snow specific surface area (m^2/kg) (npoints,nlayer) 
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowrho !snow density (kg/m^3) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowg0 ! asymmetry parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowy0 ! Value of y of snow grains at nr=1.3 (no unit
REAL, DIMENSION(:,:), INTENT(IN)   :: psnoww0 ! Value of W of snow grains at nr=1.3 (no unit)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowb0 ! absorption enhancement parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowdz !snow layers thickness (m) (npoints,nlayer)
REAL, DIMENSION(:,:,:), INTENT(IN) :: psnowimp_density !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(:,:,:), INTENT(IN) :: psnowimp_content !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
!
REAL, DIMENSION(:,:), INTENT(IN)   :: palb ! soil/vegetation albedo (npoints,nbands)
!
REAL, DIMENSION(:,:), INTENT(IN)   :: psw_rad_dif ! spectral diffuse incident light (W/m^2) (npoints,nbands)
REAL, DIMENSION(:,:), INTENT(IN)   :: psw_rad_dir ! spectral direct incident light (W/m^2) (npoints,nbands)
REAL, DIMENSION(:), INTENT(IN)     :: pcoszen ! cosine of zenithal solar angle (npoints)
!
INTEGER, DIMENSION(:), INTENT(IN)  :: knlvls_use ! number of effective snow layers (npoints)
!
REAL, DIMENSION(:,:), INTENT(OUT)   :: psnowalb !(npoints,nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowenergy !(npoints,nlayer,nbands)
REAL, DIMENSION(:,:), INTENT(OUT)   :: psoilenergy !(npoints,nbands)
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2)*2,NPNBANDS) :: zdm,zd,zdp !3 diagonals of the matrix
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2)*2,NPNBANDS) :: zvector_dir,zvector_dif
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zsnowssalb !total single scattering albedo (npoints,nlayer,nbands) 
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zsnowg !asymmetry factor (npoints,nlayer,nbands)
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zsnowalbedo! Albedo (npoints,nlayer,nbands)
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands) 
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zg_star,zssalb_star,zgamma1,zgamma2
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zdtaustar !Optical depth of each layer
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: ztaustar !Cumulated optical depth
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zgp_dir,zgm_dir ! Gp Gm vectors for direct radiation
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zgp_dif,zgm_dif ! Gp Gm vectors for diffuse radiation
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zxa_dir,zxa_dif,zxb_dir,zxb_dif,zxc_dir,zxc_dif,zxd_dir,zxd_dif
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zeprofile_dir,zeprofile_dif
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),NPNBANDS) :: zsoilabs_dir,zsoilabs_dif
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),NPNBANDS) :: zalb ! same as PALB but possibly modified if snowpack is truncated
!
REAL, DIMENSION(SIZE(PSNOWSSA,1))::zmudiff
!
INTEGER, DIMENSION(SIZE(PSNOWSSA,1),NPNBANDS) :: inlvls_eff ! number of effective snow layers
!
INTEGER, DIMENSION(NPNBANDS) :: imax_eff ! maximum number of effective layers over each band
!
INTEGER :: jb,ji,jl !loop counters
INTEGER :: imax_use ! maximum number of layers over the domain
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('TARTES',0,zhook_handle)
!
! Initialization
psnowalb    = 0.
psnowenergy = 0.
psoilenergy = 0.
zalb        = palb
!
zmudiff = xp_mudiff !the diffuse incident flux is treated as direct flux at zenithal angle 53deg
!
imax_use = maxval(knlvls_use)
!
!1 compute optical properties for each wavelength
 CALL single_scattering_optical_parameters(psnowssa,psnowrho,psnowg0,psnowy0,psnoww0,psnowb0,     &
                                           psnowimp_density,psnowimp_content,knlvls_use,imax_use, &
                                           zsnowssalb,zsnowg)

 CALL infinite_medium_optical_parameters(zsnowssalb,zsnowg,knlvls_use,imax_use,zsnowalbedo,zkestar, &
                                         zg_star,zssalb_star,zgamma1,zgamma2)

! 2 computation on every wavelength and layer of the optical depth
 CALL taustar_vector(psnowssa,psnowrho,psnowdz,zsnowssalb,zsnowg,zkestar,knlvls_use,imax_use, &
                     zdtaustar,ztaustar)

! estimate the effective layers for absorption
 CALL estimate_effective_layer_number(zkestar,zdtaustar,knlvls_use,imax_use,inlvls_eff,imax_eff)

!3 solve the radiative transfer for each wavelength
!3.1 Compte G+ and G- vectors for direct and diffuse radiations
 CALL gp_gm_vectors(zsnowssalb,zkestar,zg_star,zssalb_star,zgamma1,zgamma2,pcoszen,psw_rad_dir, &
                    inlvls_eff,imax_eff,zgp_dir,zgm_dir)
 CALL gp_gm_vectors(zsnowssalb,zkestar,zg_star,zssalb_star,zgamma1,zgamma2,zmudiff,psw_rad_dif, &
                    inlvls_eff,imax_eff,zgp_dif,zgm_dif)
!
!3.2 If the snowpack has been truncated, add a thick layer (optical) and force soil albedo to 1
DO jb = 1,npnbands
  DO ji = 1,SIZE(knlvls_use)
    IF ( inlvls_eff(ji,jb)<knlvls_use(ji) ) THEN
      zdtaustar(ji,inlvls_eff(ji,jb)+1,jb) = 30. / zkestar(ji,inlvls_eff(ji,jb)+1,jb)
      zalb(ji,jb) = 1.
    END IF
  END DO
END DO
!
!3.3 Compute the matrix and vectors
 CALL two_stream_matrix(zsnowalbedo,zalb,zkestar,zdtaustar,inlvls_eff,imax_eff,zdm,zd,zdp)
 CALL two_stream_vector(zsnowalbedo,zalb,zdtaustar,ztaustar,zgm_dir,zgp_dir,pcoszen,inlvls_eff,imax_eff,zvector_dir)
 CALL two_stream_vector(zsnowalbedo,zalb,zdtaustar,ztaustar,zgm_dif,zgp_dif,zmudiff,inlvls_eff,imax_eff,zvector_dif)
!
! DO JB=1,NPNBANDS,30
!     PRINT*,"band ",JB
!     PRINT*,ZDM(:,:,JB)
!     PRINT*,ZD(:,:,JB)
!     PRINT*,ZDP(:,:,JB)
!     PRINT*,ZVECTOR_DIR(:,:,JB)    
!     PRINT*,ZVECTOR_DIF(:,:,JB)    
! END DO
!
!3.4 solve the system
 CALL solves_two_stream2(zdm,zd,zdp,zvector_dir,zvector_dif,zsnowalbedo,psw_rad_dir,psw_rad_dif,inlvls_eff, &
                         imax_eff,zxa_dir,zxa_dif,zxb_dir,zxb_dif,zxc_dir,zxc_dif,zxd_dir,zxd_dif)
!
! DO JB=1,NPNBANDS,30
!     PRINT*,"solution band ",JB
!     PRINT*,ZXA_DIR(:,:,JB)
!     PRINT*,ZXA_DIF(:,:,JB)  
! END DO
!
!4 Diagnostics
!4.1 Albedo
 CALL snowpack_albedo(zxc_dir(:,1,:),zxc_dif(:,1,:),zxd_dir(:,1,:),zxd_dif(:,1,:),  &
                     zgp_dir(:,1,:),zgp_dif(:,1,:),pcoszen,zmudiff,psw_rad_dir,    &
                     psw_rad_dif,psnowalb)
!
!4.2 Energy profile
 CALL energy_profile(zxa_dir,zxb_dir,zxc_dir,zxd_dir,zkestar,zdtaustar,ztaustar,zgm_dir,zgp_dir,pcoszen, &
                     inlvls_eff,imax_eff,zeprofile_dir)
 CALL energy_profile(zxa_dif,zxb_dif,zxc_dif,zxd_dif,zkestar,zdtaustar,ztaustar,zgm_dif,zgp_dif,zmudiff, &
                     inlvls_eff,imax_eff,zeprofile_dif)
!
DO jb = 1,npnbands
  DO jl = 1,SIZE(psnowssa,2)
    WHERE ( jl<=inlvls_eff(:,jb) )
      psnowenergy(:,jl,jb) = psw_rad_dir(:,jb) * zeprofile_dir(:,jl,jb) + &
                             psw_rad_dif(:,jb) * zeprofile_dif(:,jl,jb)
        ENDWHERE
    END DO
END DO
!
!4.3 Soil absorption
 CALL soil_absorption(zxa_dir,zxb_dir,zkestar,zdtaustar,ztaustar,zgm_dir,pcoszen,palb,inlvls_eff,zsoilabs_dir)
 CALL soil_absorption(zxa_dif,zxb_dif,zkestar,zdtaustar,ztaustar,zgm_dif,zmudiff,palb,inlvls_eff,zsoilabs_dif)
!
psoilenergy = psw_rad_dir * zsoilabs_dir + psw_rad_dif * zsoilabs_dif
!
IF (lhook) CALL dr_hook('TARTES',1,zhook_handle)
!
END SUBROUTINE tartes
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE init_tartes()
!
!In surfex this routine has to been called only once (Crocus init).
!
IMPLICIT NONE
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('INIT_TARTES',0,zhook_handle)
!
 CALL refice() ! Interpolate refractive index for pure ice on the prescribed wavelengths
 CALL refsoot_imag() ! Compute refractive index of soot according to wavelengths from Chang (1990)
!
IF (lhook) CALL dr_hook('INIT_TARTES',1,zhook_handle)
!
END SUBROUTINE init_tartes
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE refice()
!
! Interpolate refractive index for pure ice on the prescribed wavelengths
!
USE modd_const_tartes, ONLY: npnbands,xpwavelengths,xpwavelengths_m,xrefice_r,xrefice_i, &
                             npnbands_ref,xpwavelengths_ref,xprefice_r,xprefice_i,       & 
                             xrefice_norm,xginf,xconst_c
USE modd_csts, ONLY: xpi, xrholi                           
!
USE modi_abor1_sfx
!
IMPLICIT NONE                             
!
! Log of PPWAVELENGTHS PPWAVELENGTHS_REF PPREFICE_I for interpolation
REAL, DIMENSION(NPNBANDS)     :: zlog_wl
REAL, DIMENSION(NPNBANDS_REF) :: zlog_wl_ref
REAL, DIMENSION(NPNBANDS_REF) :: zlog_refice_i
!
INTEGER :: jb,jbref !loop counters (wl band)
!
LOGICAL :: ginf !logical for interpolation
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('REFICE',0,zhook_handle)
!
zlog_wl       = log(xpwavelengths)
zlog_wl_ref   = log(xpwavelengths_ref)
zlog_refice_i = log(xprefice_i)

! Interpolate retractive index
DO jb = 1,npnbands
  !
  ginf = .true.
  !
  DO jbref = 1,npnbands_ref
    !
    IF ( xpwavelengths_ref(jbref)>xpwavelengths(jb) ) THEN
      !
      IF ( jbref<2 ) CALL abor1_sfx("FATAL ERROR INIT_TARTES (interpolation of refractive indexs)")
      !
      ginf = .false.
      !
      xrefice_r(jb) = ( (xpwavelengths(jb)    - xpwavelengths_ref(jbref-1)) * xprefice_r(jbref)   +   &
                        (xpwavelengths_ref(jbref) - xpwavelengths(jb)     ) * xprefice_r(jbref-1) ) / &
                       ( xpwavelengths_ref(jbref) - xpwavelengths_ref(jbref-1) )
      xrefice_i(jb) = exp( ( (zlog_wl(jb)    - zlog_wl_ref(jbref-1)) * zlog_refice_i(jbref)   +   &
                             (zlog_wl_ref(jbref) - zlog_wl(jb)     ) * zlog_refice_i(jbref-1) ) / &
                            ( zlog_wl_ref(jbref) - zlog_wl_ref(jbref-1) ) )
      !
      xrefice_norm(jb) = xrefice_r(jb) - 1.3 !factor in eq 72-73-74
      xginf(jb) = 0.9751 - 0.105 * xrefice_norm(jb) !doc equation 72
      xconst_c(jb) = 24. * xpi * xrefice_i(jb) / ( xrholi * xpwavelengths_m(jb) ) !constant c*SSA in equation 71
      !
      EXIT
      !
    END IF
  !
  END DO
  !
  IF ( ginf ) CALL abor1_sfx("FATAL ERROR INIT_TARTES (interpolation of refractive indexs)")
  !
END DO
!
IF (lhook) CALL dr_hook('REFICE',1,zhook_handle)
!
END SUBROUTINE refice
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE refsoot_imag()

! Compute refractive index of soot according to wavelengths from Chang (1990)

USE modd_const_tartes, ONLY : npnbands,xpwavelengths,xrefimp_i
!
!PPWAVELENGTHS nanometers
!
IMPLICIT NONE  
!
REAL, DIMENSION    (NPNBANDS) :: zwl_um ! Wavelengths in micrometers (Chang, 1990 formulas)
REAL, DIMENSION    (NPNBANDS) :: zindex_soot_real,zindex_soot_imag ! real and imaginary components of refractive index
 COMPLEX, DIMENSION(NPNBANDS) :: zindex_soot !complex refractive index
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('REFSOOT_IMAG',0,zhook_handle)
!
zwl_um = xpwavelengths / 1000.
!
zindex_soot_real = 1.811  + 0.1263*log(zwl_um) + 0.027 *log(zwl_um)**2 + 0.0417*log(zwl_um)**3
zindex_soot_imag = 0.5821 + 0.1213*log(zwl_um) + 0.2309*log(zwl_um)**2 - 0.01  *log(zwl_um)**3
!
zindex_soot = zindex_soot_real - cmplx(0,1) * zindex_soot_imag
!
! absorption cross section of small particles (Bohren and Huffman, 1983)
xrefimp_i(:,1) = aimag( (zindex_soot**2-1.) / (zindex_soot**2 + 2.) )
!
IF (lhook) CALL dr_hook('REFSOOT_IMAG',1,zhook_handle)
!
END SUBROUTINE refsoot_imag
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE shape_parameter_variations(PSNOWG0,PSNOWY0,PSNOWW0,PSNOWB0,PSNOWG00,PSNOWY,PSNOWW,PSNOWB)

!compute shape parameter variations as a function of the the refraction index with respect to the value in the visible range.
!These variation equations were obtained for sphere (Light Scattering Media Optics, Kokhanovsky, A., p.61) but should also apply to other shapes in a first approximation.
!see doc Section 2
!
USE modd_const_tartes, ONLY: npnbands,xrefice_norm !number of spectral bands
!
IMPLICIT NONE  
!
REAL, DIMENSION(:,:), INTENT(IN) :: psnowg0 !asymmetry parameter of snow grains at refractive index=1.3 and at non absorbing wavelengths (no unit)   (npoints*nlayers)
REAL, DIMENSION(:,:), INTENT(IN) :: psnowy0 !Value of y of snow grains at  refractive index=1.3 (no unit) (npoints*nlayers)
REAL, DIMENSION(:,:), INTENT(IN) :: psnoww0 !Value of W of snow grains at  refractive index=1.3 (no unit) (npoints*nlayers)
REAL, DIMENSION(:,:), INTENT(IN) :: psnowb0 !absorption enhancement parameter of snow grains at refractive index=1.3 and at non absorbing wavelengths (no unit) (npoints*nlayers)
!
! Spectral parameters necessary to compute the asymmetry parameter and single scattering albedo of snow. For now, those parameters do not evolve with time. They depend on shape only
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowg00 !(npoints*nlayers*nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowy !(npoints*nlayers*nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnoww !(npoints*nlayers*nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowb !(npoints*nlayers*nbands)
!
INTEGER :: jb !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SHAPE_PARAMETER_VARIATIONS',0,zhook_handle)
!
DO jb = 1,npnbands
  psnowg00(:,:,jb) = psnowg0(:,:) - 0.38  * xrefice_norm(jb) !doc equation 73
  psnowb(:,:,jb) = psnowb0(:,:) + 0.4   * xrefice_norm(jb) !doc equation 77
  psnoww(:,:,jb) = psnoww0(:,:) + 0.17  * xrefice_norm(jb) 
  psnowy(:,:,jb) = psnowy0(:,:) + 0.752 * xrefice_norm(jb) !doc equation 74
END DO
!
IF (lhook) CALL dr_hook('SHAPE_PARAMETER_VARIATIONS',1,zhook_handle)
!
END SUBROUTINE shape_parameter_variations
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE impurities_co_single_scattering_albedo(PSNOWSSA,PSNOWIMP_DENSITY,PSNOWIMP_CONTENT, &
                                                  knlvls_use,kmax_use,pcossalb)
!
USE modd_const_tartes, ONLY: npnbands,npnimp,xpwavelengths_m,xrefimp_i
USE modd_csts, ONLY: xpi
!
IMPLICIT NONE
!
!return the spectral absorption of snow due to the impurities
!see doc Section 2.6
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowssa !snow specific surface area (m^2/kg) (npoints,nlayer) 
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowimp_density !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowimp_content !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use !number of active layers
INTEGER, INTENT(IN)                 :: kmax_use !maximum number of active layers over the domain
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pcossalb !co single scattering albedo of impurities
!
REAL,DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2)) :: zabs_imp
!
INTEGER :: jimp !loop counter
INTEGER :: jb, jl,jj !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('IMPURITIES_CO_SINGLE_SCATTERING_ALBEDO',0,zhook_handle)
!
pcossalb = 0.
!
DO jb = 1,npnbands
  DO jimp = 1,npnimp
    DO jl = 1,kmax_use
      DO jj = 1,SIZE(knlvls_use)
        !
        IF ( knlvls_use(jj)>=jl ) THEN
          !
          zabs_imp(jj,jl)       = -xrefimp_i(jb,jimp)
          pcossalb(jj,jl,jb) = pcossalb(jj,jl,jb) + &
                                      12. * xpi / ( xpwavelengths_m(jb)*psnowssa(jj,jl) ) * &
                                      psnowimp_content(jj,jl,jimp) / psnowimp_density(jj,jl,jimp) * &
                                      zabs_imp(jj,jl) !doc equation 79
          !
        ENDIF
        !
      ENDDO
    ENDDO
  ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('IMPURITIES_CO_SINGLE_SCATTERING_ALBEDO',1,zhook_handle)
!
END SUBROUTINE impurities_co_single_scattering_albedo
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE single_scattering_optical_parameters(PSNOWSSA,PSNOWRHO,PSNOWG0,PSNOWY0,PSNOWW0,PSNOWB0,     &
                                                psnowimp_density,psnowimp_content,knlvls_use,kmax_use, &
                                                psnowssalb,psnowg)
!
!see doc Section 2.3, 2.5, 2.6
USE modd_const_tartes, ONLY: npnbands,xginf,xconst_c
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowssa !snow specific surface area (m^2/kg) (npoints,nlayer) 
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowrho !snow density (kg/m^3) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowg0 ! asymmetry parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowy0 ! Value of y of snow grains at nr=1.3 (no unit
REAL, DIMENSION(:,:), INTENT(IN)    :: psnoww0 ! Value of W of snow grains at nr=1.3 (no unit)
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowb0 ! absorption enhancement parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowimp_density !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowimp_content !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use !number of active layers
INTEGER, INTENT(IN)                 :: kmax_use !maximum number of active layers over the domain
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowssalb !total single scattering albedo (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowg !asymmetry factor (npoints,nlayer,nbands) 
!
!Local variables
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zsnowg00,zsnowy,zsnoww,zsnowb
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zsnowcossalb ! co- single scattering albedo of pure snow
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2),NPNBANDS) :: zimpcossalb ! co- single scattering albedo of impurities
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2)) :: zc,zphi
!
INTEGER :: jb,jl,jj !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SINGLE_SCATTERING_OPTICAL_PARAMETERS',0,zhook_handle)
!
 CALL shape_parameter_variations(psnowg0,psnowy0,psnoww0,psnowb0,zsnowg00,zsnowy,zsnoww,zsnowb)
!
DO jl = 1,kmax_use
  !
  DO jj =1,SIZE(psnowssa,1)
    !
    IF ( knlvls_use(jj)>=jl ) THEN
      !
      DO jb = 1,npnbands
        !
        ! calculation of the spectral asymmetry parameter of snow
        zc(jj,jl) = xconst_c(jb) / psnowssa(jj,jl)
        !    
        psnowg(jj,jl,jb) = xginf(jb) - ( xginf(jb)-zsnowg00(jj,jl,jb) ) * exp( -zsnowy(jj,jl,jb)*zc(jj,jl) )
       !
        ! co- single scattering albedo of pure snow
        zphi(jj,jl)       = 2./3. * zsnowb(jj,jl,jb) / ( 1.-zsnoww(jj,jl,jb) )
        zsnowcossalb(jj,jl,jb) = 0.5 * ( 1.-zsnoww(jj,jl,jb) ) * ( 1.-exp( -zphi(jj,jl)*zc(jj,jl) ) ) !doc equation 76
        !    
      ENDDO
      !
    ENDIF
    !
  ENDDO
  !
ENDDO
!
!adding co- single scattering albedo for impureties
 CALL impurities_co_single_scattering_albedo(psnowssa,psnowimp_density,psnowimp_content,&
                                            knlvls_use,kmax_use,zimpcossalb)
!
zsnowcossalb = zsnowcossalb + zimpcossalb
!
!total single scattering albedo
psnowssalb = 1.-zsnowcossalb
!
IF (lhook) CALL dr_hook('SINGLE_SCATTERING_OPTICAL_PARAMETERS',1,zhook_handle)
!
END SUBROUTINE single_scattering_optical_parameters
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE infinite_medium_optical_parameters(PSNOWSSALB,PSNOWG,KNLVLS_USE,KMAX_USE,PSNOWALBEDO,PKESTAR,&
                                              pg_star,pssalb_star,pgamma1,pgamma2)
!return albedo and kestar using Delta-Eddington Approximation (The Delta-Eddington Approximation of Radiative Flux Transfer, Jospeh et al (1976)).  
! Fluxes in the snowpack depend on these 2 quantities
! see doc section 1.4
!
USE modd_const_tartes, ONLY: npnbands
USE modd_snow_metamo, ONLY: xuepsi
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowssalb !total single scattering albedo (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowg !asymmetry factor (npoints,nlayer,nbands) 
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use !number of active layers
INTEGER, INTENT(IN)                 :: kmax_use !maximum number of active layers over the domain
REAL, DIMENSION(:,:,:), INTENT(OUT) :: psnowalbedo ! Albedo (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pg_star,pssalb_star,pgamma1,pgamma2 !(npoints,nlayer,nbands) 
!
INTEGER :: jb,jl,jj !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('INFINITE_MEDIUM_OPTICAL_PARAMETERS',0,zhook_handle)
!
DO jb = 1,npnbands
  !
  DO jl = 1,kmax_use
    !
    DO jj =1,SIZE(psnowg,1)
      !
      IF ( knlvls_use(jj)>=jl ) THEN
        !
        pg_star(jj,jl,jb) = psnowg(jj,jl,jb) / ( 1. + psnowg(jj,jl,jb) ) !doc equation 12
        pssalb_star(jj,jl,jb) = psnowssalb(jj,jl,jb) * ( 1. - psnowg(jj,jl,jb)**2 ) / &
                                      ( 1. - psnowg(jj,jl,jb)**2 * psnowssalb(jj,jl,jb) ) !doc equation 16
        !
        ! Jimenez-Aquino, J. and Varela, J. R., (2005)
        pgamma1(jj,jl,jb) =  0.25 * ( 7. - pssalb_star(jj,jl,jb)*(4.+3.*pg_star(jj,jl,jb)) )      !doc equation 38
        pgamma2(jj,jl,jb) = -0.25 * ( 1. - pssalb_star(jj,jl,jb)*(4.-3.*pg_star(jj,jl,jb)) )      !doc equation 39
        !
        pkestar(jj,jl,jb) = sqrt( pgamma1(jj,jl,jb)**2 - pgamma2(jj,jl,jb)**2 )                     !doc equation 42
        psnowalbedo(jj,jl,jb) = ( pgamma1(jj,jl,jb)-pkestar(jj,jl,jb) ) / pgamma2(jj,jl,jb)  !doc equation 43
        !
        ! Modif M Lafaysse to avoid division by 0
        !NB JJ note that this variable can be negative in infra-red wavelengths (it represents the albedo only in smallest wavelengths
        IF ( abs(psnowalbedo(jj,jl,jb))<xuepsi ) THEN
          psnowalbedo(jj,jl,jb) = sign( xuepsi, psnowalbedo(jj,jl,jb) )
        ENDIF
        !
      ENDIF
      !
    ENDDO
    !
  ENDDO
  !
ENDDO
!
IF (lhook) CALL dr_hook('INFINITE_MEDIUM_OPTICAL_PARAMETERS',1,zhook_handle)
!
END SUBROUTINE infinite_medium_optical_parameters
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
!
SUBROUTINE taustar_vector(PSNOWSSA,PSNOWRHO,PSNOWDZ,PSNOWSSALB,PSNOWG,PKESTAR,KNLVLS_USE,KMAX_USE,&
                          pdtaustar,ptaustar)
!compute the taustar and dtaustar of the snowpack, the optical depth of each layer and cumulated optical depth
!see doc Section 1.2, 1.8, 2.4
!
USE modd_const_tartes, ONLY: npnbands,xpmax_opticaldepth
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowssa !snow specific surface area (m^2/kg) (npoints,nlayer) 
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowrho !snow density (kg/m^3) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowdz !snow depth (m) (npoints,nlayer)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowssalb  !total single scattering albedo (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowg ! asymmetry factor (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands) 
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use !number of active layers
INTEGER, INTENT(IN)                 :: kmax_use !maximum number of active layers over the domain
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pdtaustar !Optical depth of each layer
REAL, DIMENSION(:,:,:), INTENT(OUT) :: ptaustar !Cumulated optical depth
!
REAL, DIMENSION(SIZE(PSNOWSSA,1),SIZE(PSNOWSSA,2)) :: zsigext  ! Extinction coefficient (npoints,nlayer) 
!
INTEGER :: jb,jl !loop counters
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('TAUSTAR_VECTOR',0,zhook_handle)
!
!Compute extinction coefficient
zsigext = psnowrho*psnowssa/2. !doc equation 75
!
DO jb=1,npnbands
  !
  DO jl = 1,kmax_use
    !
    WHERE ( knlvls_use>=jl )
      !Optical depth of each layer with delta-eddington variable change, doc equation 15
      ! + optical depth threshold (doc section 1.8)
      pdtaustar(:,jl,jb) = min( zsigext(:,jl) * psnowdz(:,jl) * ( 1.- psnowssalb(:,jl,jb)*psnowg(:,jl,jb)**2 ), &
                                xpmax_opticaldepth / pkestar(:,jl,jb) )
    ENDWHERE
    !
  END DO
  !  
  !Cumulated optical depth
  !First layer
  ptaustar(:,1,jb) = pdtaustar(:,1,jb)
  !Other layers
  DO jl = 2,kmax_use
    WHERE ( knlvls_use>=jl )
      ptaustar(:,jl,jb) = ptaustar(:,jl-1,jb) + pdtaustar(:,jl,jb)
    ENDWHERE
  END DO
  !  
END DO

IF (lhook) CALL dr_hook('TAUSTAR_VECTOR',1,zhook_handle)

END SUBROUTINE taustar_vector
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE estimate_effective_layer_number(PKESTAR,PDTAUSTAR,KNLVLS_USE,KMAX_USE,KNLVLS_EFF,KMAX_EFF)
! estimate the number of layers to take into account at each wavelength
!doc section 1.8
!
USE modd_const_tartes, ONLY : npnbands,xpwavelengths,xptaumax
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)   :: pkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(IN)   :: pdtaustar !Optical depth of each layer
INTEGER, DIMENSION(:), INTENT(IN)    :: knlvls_use !number of active layers
INTEGER, INTENT(IN)                  :: kmax_use !maximum number of active layers over the domain
INTEGER, DIMENSION(:,:), INTENT(OUT) :: knlvls_eff !number of effective layers (npoints,nbands)
INTEGER, DIMENSION(:), INTENT(OUT)   :: kmax_eff !maximum number of effective layers over the domain (nbands)
!
REAL, DIMENSION(SIZE(PKESTAR,1),SIZE(PKESTAR,2),NPNBANDS) :: ztau
LOGICAL, DIMENSION(SIZE(PKESTAR,1)) :: geff
INTEGER :: jb,jl !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('ESTIMATE_EFFECTIVE_LAYER_NUMBER',0,zhook_handle)
!
DO jb = 1,npnbands
  !
  knlvls_eff(:,jb) = knlvls_use
  ztau(:,1,jb) = pkestar(:,1,jb) * pdtaustar(:,1,jb)
  geff = .true.
  !
  DO jl = 2,kmax_use
    !
    WHERE ( (knlvls_use>=jl) .AND. geff )
      ztau(:,jl,jb) = ztau(:,jl-1,jb) + pkestar(:,jl,jb) * pdtaustar(:,jl,jb)
    ELSEWHERE
      ztau(:,jl,jb) = 0.      
    ENDWHERE
    WHERE ( ztau(:,jl,jb)>xptaumax )
      knlvls_eff(:,jb) = max(1,jl-1)
      geff = .false.
    ENDWHERE
    !
  END DO
  !  
  kmax_eff(jb) = maxval(knlvls_eff(:,jb))
  !  
END DO

IF (lhook) CALL dr_hook('ESTIMATE_EFFECTIVE_LAYER_NUMBER',1,zhook_handle)

END SUBROUTINE estimate_effective_layer_number
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
!
SUBROUTINE gp_gm_vectors(PSNOWSSALB,PKESTAR,PG_STAR,PSSALB_STAR,PGAMMA1,PGAMMA2,PCOSZEN,PSW_RAD,KNLVLS_EFF,KMAX_EFF,PGP,PGM)
!
!return GP and GM vectors of equations 40/41 and 46/47
! (equations for the downward and upward fluxes in the snowpack for the 2 stream approximation)
!
USE modd_const_tartes, ONLY : npnbands
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowssalb !total single scattering albedo (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands) 
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pg_star ! asymmetry factor * (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pssalb_star
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pgamma1,pgamma2
REAL, DIMENSION(:,:), INTENT(IN)    :: psw_rad ! incident radiation (direct or diffuse) (npoints,nbands)
REAL, DIMENSION(:), INTENT(IN)      :: pcoszen ! cosine of zenithal solar angle (npoints)
INTEGER, DIMENSION(:,:), INTENT(IN) :: knlvls_eff !number of effective layers (npoints,nbands)
INTEGER, DIMENSION(:), INTENT(IN)   :: kmax_eff !maximum number of effective layers over the domain (nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pgp !GP vector (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pgm !GM vector (npoints,nlayer,nbands)
!
REAL :: zgamma3,zgamma4,zg ! intermediate terms
!
INTEGER :: jb,jl,jj !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('GP_GM_VECTORS',0,zhook_handle)
!
!Init
pgp = 0.
pgm = 0.
!
DO jb = 1,npnbands
  !
  DO jl = 1,kmax_eff(jb)
    !
    DO jj =1,SIZE(psw_rad,1)
      !
      IF ( psw_rad(jj,jb)>0. .AND. knlvls_eff(jj,jb)>=jl ) THEN
        !
        zgamma3 = 0.25 * ( 2. - 3.*pg_star(jj,jl,jb)*pcoszen(jj) ) !doc equation 28
        zgamma4 = 0.25 * ( 2. + 3.*pg_star(jj,jl,jb)*pcoszen(jj) ) !doc equation 27
        zg = pcoszen(jj)**2 * pssalb_star(jj,jl,jb) / ( (pkestar(jj,jl,jb)*pcoszen(jj))**2 - 1. ) !factor eq 44-45
        pgp(jj,jl,jb) = zg * ( (pgamma1(jj,jl,jb)-1./pcoszen(jj))*zgamma3 + pgamma2(jj,jl,jb)*zgamma4 ) !doc equation 45
        pgm(jj,jl,jb) = zg * ( (pgamma1(jj,jl,jb)+1./pcoszen(jj))*zgamma4 + pgamma2(jj,jl,jb)*zgamma3 ) !doc equation 44
        !
      ENDIF
      !
      ENDDO
    !
  END DO
  !
END DO
!
IF (lhook) CALL dr_hook('GP_GM_VECTORS',1,zhook_handle)
!
END SUBROUTINE gp_gm_vectors
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE two_stream_matrix(PSNOWALBEDO,PSOILALBEDO,PKESTAR,PDTAUSTAR,KNLVLS_EFF,KMAX_EFF,PDM,PD,PDP)
!compute the matrix in the system describing the continuity and boundary conditions at one point and one wavelength.
! see doc section 1.5
!
USE modd_const_tartes, ONLY : npnbands
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowalbedo ! snow albedo (npoints,nlayer,nbands)
REAL, DIMENSION(:,:), INTENT(IN)    :: psoilalbedo ! soil albedo (npoints,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pkestar     !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pdtaustar   !Optical depth of each layer
INTEGER, DIMENSION(:,:), INTENT(IN) :: knlvls_eff  !number of effective layers (npoints,nbands)
INTEGER, DIMENSION(:), INTENT(IN)   :: kmax_eff    !maximum number of effective layers over the domain (nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pdm,pd,pdp  ! three diagnonals of the matrix
!
REAL, DIMENSION(SIZE(PSNOWALBEDO,1)) :: zfdiag
!
REAL :: zfdiag2
!
INTEGER :: jb,jl,ji !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('TWO_STREAM_MATRIX',0,zhook_handle)
!
!Initialization
pdm = 0.
pd  = 0.
pdp = 0.
!
DO jb = 1,npnbands
  !
  DO jl = 1,kmax_eff(jb)-1
    !
!       IF(PSNOWALBEDO(1,JL,JB)<0.000001) THEN
!           PRINT*,"WARNING PSNOWALBEDO LAYER ",JL," BAND ",JB," : ",PSNOWALBEDO(1,JL,JB)
!       END IF


!       IF (EXP(-PKESTAR(1,JL,JB)*PDTAUSTAR(1,JL,JB))<0.000001) THEN
!           PRINT*,"WARNING ZFDIAG ",JL," BAND ",JB," : ",EXP(-PKESTAR(1,JL,JB)*PDTAUSTAR(1,JL,JB))     
!       END IF
!
    DO ji =1,SIZE(knlvls_eff,1)
      !
      IF ( jl<=knlvls_eff(ji,jb)-1 ) THEN
        !
        !See matrix documentation page 8 and formal expressions page 9
        zfdiag(ji) = exp( -pkestar(ji,jl,jb)*pdtaustar(ji,jl,jb) )
        !    
        !Décalage d'un indice vers la droite par rapport au code python
        pdm(ji,jl*2,jb)   = ( 1. - psnowalbedo(ji,jl,jb)*psnowalbedo(ji,jl+1,jb) ) * zfdiag(ji)
        pdm(ji,jl*2+1,jb) = ( 1./psnowalbedo(ji,jl,jb) - psnowalbedo(ji,jl,jb) )   * 1./zfdiag(ji)
        !  
        pd(ji,jl*2,jb)    = ( 1. - psnowalbedo(ji,jl+1,jb)/psnowalbedo(ji,jl,jb) ) * 1./zfdiag(ji)
        pd(ji,jl*2+1,jb)  = psnowalbedo(ji,jl,jb) - psnowalbedo(ji,jl+1,jb)
        !    
        !Décalage d'un indice vers la gauche par rapport au code python
        pdp(ji,jl*2,jb)   = psnowalbedo(ji,jl+1,jb) * psnowalbedo(ji,jl+1,jb) - 1.
        pdp(ji,jl*2+1,jb) = psnowalbedo(ji,jl,jb) - 1./psnowalbedo(ji,jl+1,jb)
        !
      ENDIF
      !
    ENDDO
    !
  ENDDO
  !  
  pdp(:,1,jb) = 1. !Décalage d'un indice vers la gauche par rapport au code python
  pd(:,1,jb) = 1.
  !  
  DO ji=1,SIZE(psnowalbedo,1)
    !
    zfdiag2 = exp( -pkestar(ji,knlvls_eff(ji,jb),jb) * pdtaustar(ji,knlvls_eff(ji,jb),jb) )
    !      
    !Décalage d'un indice vers la droite par rapport au code python
    pdm(ji,2*knlvls_eff(ji,jb),jb) = zfdiag2    * &
                                      ( psnowalbedo(ji,knlvls_eff(ji,jb),jb)    - psoilalbedo(ji,jb) )
    !
    pd(ji,2*knlvls_eff(ji,jb),jb) = 1./zfdiag2 * &
                                      ( 1./psnowalbedo(ji,knlvls_eff(ji,jb),jb) - psoilalbedo(ji,jb) )
  END DO
  !  
END DO
!
IF (lhook) CALL dr_hook('TWO_STREAM_MATRIX',1,zhook_handle)
!
END SUBROUTINE two_stream_matrix
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE two_stream_vector(PSNOWALBEDO,PSOILALBEDO,PDTAUSTAR,PTAUSTAR,PGM,PGP,PCOSZEN,KNLVLS_EFF,KMAX_EFF,PVECTOR)
!compute the V vector in the system describing the continuity and boundary conditions at one point and one wavelength.
! see doc section 1.5
!
USE modd_const_tartes, ONLY : npnbands
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowalbedo  ! snow albedo (npoints,nlayer,nbands)
REAL, DIMENSION(:,:), INTENT(IN)    :: psoilalbedo  ! soil albedo (npoints,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pdtaustar    !Optical depth of each layer
REAL, DIMENSION(:,:,:), INTENT(IN)  :: ptaustar     !Cumulated optical depth
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pgp          !GP vector (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pgm !GM vector (npoints,nlayer,nbands)
REAL, DIMENSION(:), INTENT(IN)      :: pcoszen ! cosine of zenithal solar angle (npoints)
INTEGER, DIMENSION(:,:), INTENT(IN) :: knlvls_eff !number of effective layers (npoints,nbands)
INTEGER, DIMENSION(:), INTENT(IN)   :: kmax_eff !maximum number of effective layers over the domain (nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pvector !output vector V
!
REAL :: zdgp,zdgm,zexp
!
INTEGER :: jb,jl,ji !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('TWO_STREAM_VECTOR',0,zhook_handle)
!
pvector(:,1,:) = -pgm(:,1,:)
!
DO jb = 1,npnbands
  !
  DO ji = 1,SIZE(psnowalbedo,1)
    !
    DO jl = 1,kmax_eff(jb)
      !
      IF ( jl<=knlvls_eff(ji,jb)-1 ) THEN
        !
        zdgp = pgp(ji,jl+1,jb) - pgp(ji,jl,jb) !doc equation 58
        zdgm = pgm(ji,jl+1,jb) - pgm(ji,jl,jb) !doc equation 58
        !
        zexp = exp( -ptaustar(ji,jl,jb)/pcoszen(ji) )
        !see expression doc page 9
        pvector(ji,2*jl,jb)   = ( zdgm - psnowalbedo(ji,jl+1,jb) * zdgp ) * zexp 
        pvector(ji,2*jl+1,jb) = ( zdgp - psnowalbedo(ji,jl,jb)   * zdgm ) * zexp 
        !
      END IF
      !
    END DO
    !
    pvector(ji,2*knlvls_eff(ji,jb),jb) = ( psoilalbedo(ji,jb) * &
                                            ( pgm(ji,knlvls_eff(ji,jb),jb) + pcoszen(ji) ) - &
                                           pgp(ji,knlvls_eff(ji,jb),jb) ) * &
                                            exp( -ptaustar(ji,knlvls_eff(ji,jb),jb) / pcoszen(ji) )
    !
  END DO
  !
END DO
!
IF (lhook) CALL dr_hook('TWO_STREAM_VECTOR',1,zhook_handle)
!
END SUBROUTINE two_stream_vector
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE solves_two_stream2(PDM,PD,PDP,PVECT_DIR,PVECT_DIF,PSNOWALBEDO,PSW_RAD_DIR,PSW_RAD_DIF,KNLVLS_EFF, &
                              kmax_eff,pxa_dir,pxa_dif,pxb_dir,pxb_dif,pxc_dir,pxc_dif,pxd_dir,pxd_dif)
!solve the two stream linear system for both direct and diffuse radiation
!
USE modd_const_tartes, ONLY : npnbands
!
USE modi_tridiag_ground_snowcro
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pdm,pd,pdp ! three diagnonals of the matrix
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pvect_dir,pvect_dif !two-stream vector V (2 vectors are used when ther is diffuse AND direct incident light)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: psnowalbedo! snow albedo (npoints,nlayer,nbands)
REAL, DIMENSION(:,:), INTENT(IN)    :: psw_rad_dir,psw_rad_dif !not used for now : useful if we want to pack the points excluding zero radiations
INTEGER, DIMENSION(:,:), INTENT(IN) :: knlvls_eff !number of effective layers (npoints,nbands)
INTEGER, DIMENSION(:), INTENT(IN)   :: kmax_eff !maximum number of effective layers over the domain (nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: pxa_dir,pxa_dif,pxb_dir,pxb_dif,pxc_dir,pxc_dif,pxd_dir,pxd_dif ! solutions (coeffs A and B in eq 46-47 or 48-49)
!
REAL, DIMENSION(SIZE(PDM,1),2*SIZE(PSNOWALBEDO,2),NPNBANDS) :: zx0_dir,zx0_dif
INTEGER :: jb,ji,jl !loop counters
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SOLVES_TWO_STREAM2',0,zhook_handle)
!
! for now we inverse the matrix twice :
! to be improved by adding a dimension in tridiag_ground_snowcro
 CALL tridiag_ground_snowcro(pdm(:,:,:),pd(:,:,:),pdp(:,:,:), &
                            pvect_dir(:,:,:),zx0_dir(:,:,:),  &
                            2*knlvls_eff(:,:),0)
!
 CALL tridiag_ground_snowcro(pdm(:,:,:),pd(:,:,:),pdp(:,:,:), &
                            pvect_dif(:,:,:),zx0_dif(:,:,:),  &
                            2*knlvls_eff(:,:),0)
!
!for now we always compute everything
DO jb = 1,npnbands
  !  
  DO jl=1,kmax_eff(jb)
    !
    DO ji=1,SIZE(pdm,1)
      !
      IF ( jl<=knlvls_eff(ji,jb) ) THEN
        pxa_dir(ji,jl,jb) = zx0_dir(ji,jl*2-1,jb)
        pxa_dif(ji,jl,jb) = zx0_dif(ji,jl*2-1,jb)
        pxb_dir(ji,jl,jb) = zx0_dir(ji,jl*2,jb)
        pxb_dif(ji,jl,jb) = zx0_dif(ji,jl*2,jb)
        !
        pxc_dir(ji,jl,jb) = pxa_dir(ji,jl,jb) * psnowalbedo(ji,jl,jb)
        pxc_dif(ji,jl,jb) = pxa_dif(ji,jl,jb) * psnowalbedo(ji,jl,jb)
        pxd_dir(ji,jl,jb) = pxb_dir(ji,jl,jb) / psnowalbedo(ji,jl,jb)
        pxd_dif(ji,jl,jb) = pxb_dif(ji,jl,jb) / psnowalbedo(ji,jl,jb)
      END IF
      !
    END DO
    !
  END DO
  !
END DO
!
IF (lhook) CALL dr_hook('SOLVES_TWO_STREAM2',1,zhook_handle)
!
END SUBROUTINE solves_two_stream2
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE snowpack_albedo(PXC_DIR,PXC_DIF,PXD_DIR,PXD_DIF,PGP_DIR,PGP_DIF,&
                           pcoszen_dir,pcoszen_dif,psw_rad_dir,psw_rad_dif,psnowalb)
! compute the albedo of the snowpack at one wavelength
!
USE modd_const_tartes, ONLY : npnbands
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pxc_dir,pxc_dif,pxd_dir,pxd_dif ! for first level (npoints*nbands)
REAL, DIMENSION(:,:), INTENT(IN)  :: pgp_dir,pgp_dif ! for first level (npoints*nbands)
REAL, DIMENSION(:), INTENT(IN)    :: pcoszen_dir,pcoszen_dif  ! cosine of zenithal solar angle (npoints)
REAL, DIMENSION(:,:), INTENT(IN)  :: psw_rad_dir,psw_rad_dif  ! incident radiation W/m^2 (npoints*nbands)
REAL, DIMENSION(:,:), INTENT(OUT) :: psnowalb ! albedo at one wavelength
!
REAL,DIMENSION(SIZE(PSW_RAD_DIR,1)) :: zref_dir,zref_dif ! reflected direct and diffuse radiations W/m^2 
                                                         ! for one band
REAL,DIMENSION(SIZE(PSW_RAD_DIR,1)) :: zinc ! incident radiation
!
INTEGER :: jb !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWPACK_ALBEDO',0,zhook_handle)
!
DO jb = 1,npnbands
  !
  ! Doc equation 66 (separated in direct and diffuse components)
  zref_dir = ( pxc_dir(:,jb)+pxd_dir(:,jb)+pgp_dir(:,jb) ) * psw_rad_dir(:,jb)
  zref_dif = ( pxc_dif(:,jb)+pxd_dif(:,jb)+pgp_dif(:,jb) ) * psw_rad_dif(:,jb) 
  zinc = psw_rad_dir(:,jb)*pcoszen_dir + psw_rad_dif(:,jb)*pcoszen_dif
  WHERE ( zinc>0. )
    psnowalb(:,jb) = (zref_dir+zref_dif) / zinc
  ELSEWHERE
    psnowalb(:,jb) = 0.
  ENDWHERE
  !  
END DO

IF (lhook) CALL dr_hook('SNOWPACK_ALBEDO',1,zhook_handle)

END SUBROUTINE snowpack_albedo
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE energy_profile(PXA,PXB,PXC,PXD,PKESTAR,PDTAUSTAR,PTAUSTAR,PGM,PGP,PCOSZEN,KNLVLS_EFF,KMAX_EFF,PEPROFILE)

    !compute energy absorption for each layer and wavelength
USE modd_const_tartes, ONLY : npnbands
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pxa,pxb,pxc,pxd ! solutions of linear system (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pdtaustar !Optical depth of each layer (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: ptaustar !Cumulated optical depth (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pgp !GP vector (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pgm !GM vector (npoints,nlayer,nbands)
REAL, DIMENSION(:), INTENT(IN)      :: pcoszen! cosine of zenithal solar angle (npoints)
INTEGER, DIMENSION(:,:), INTENT(IN) :: knlvls_eff !number of effective layers (npoints,nbands)
INTEGER, DIMENSION(:), INTENT(IN)   :: kmax_eff !maximum number of effective layers over the domain (nbands)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: peprofile ! energy absorbed by each layer (W/m^2) npoints,nlayer,nbands)
!
REAL :: zdexp, zfdu, zfdd, zstar
!
INTEGER::jb,jl,jj !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('ENERGY_PROFILE',0,zhook_handle)

DO jb = 1,npnbands
  !
  DO jj =1,SIZE(peprofile,1)
    !
    zstar = pkestar(jj,1,jb) * pdtaustar(jj,1,jb)
    !
    !surface layer doc equation 64
    peprofile(jj,1,jb) = ( pcoszen(jj) - ( pxc(jj,1,jb)+pxd(jj,1,jb)+pgp(jj,1,jb) ) ) + &
                            ( pxc(jj,1,jb) * exp(-zstar) + pxd(jj,1,jb) * exp(zstar) + &
                              pgp(jj,1,jb) * exp( -pdtaustar(jj,1,jb)/pcoszen(jj)) ) - &
                            ( pxa(jj,1,jb) * exp(-zstar) + pxb(jj,1,jb) * exp(zstar) + &
                              pgm(jj,1,jb) * exp( -pdtaustar(jj,1,jb)/pcoszen(jj)) + &
                               pcoszen(jj) * exp( -ptaustar(jj,1,jb)/pcoszen(jj)) ) 
    !
    !internal layers
    ! 
    DO jl = 2,kmax_eff(jb)
      !  
      zstar = pkestar(jj,jl,jb) * pdtaustar(jj,jl,jb)
      !
      IF ( jl<=knlvls_eff(jj,jb) ) THEN
        !
        !last factor in equations 62 and 63
        zdexp = exp( -ptaustar(jj,jl  ,jb)/pcoszen(jj) ) - exp( -ptaustar(jj,jl-1,jb)/pcoszen(jj) )
        !
        !doc equation 62
        zfdu = pxc(jj,jl,jb) * ( exp(-zstar) -1. ) + &
               pxd(jj,jl,jb) * ( exp( zstar) -1. ) + pgp(jj,jl,jb) * zdexp
        !
        !doc equation 63
        zfdd = pxa(jj,jl,jb) * ( exp(-zstar) -1. ) + &
               pxb(jj,jl,jb) * ( exp( zstar) -1. ) + ( pgm(jj,jl,jb) + pcoszen(jj) ) * zdexp
        !      
        peprofile(jj,jl,jb) = zfdu - zfdd !doc equation 61
        !
      ELSE
        !
        peprofile(jj,jl,jb) = 0.
        !
      ENDIF
      !
    ENDDO
    !
  ENDDO
  !
ENDDO
!
IF (lhook) CALL dr_hook('ENERGY_PROFILE',1,zhook_handle)
!
END SUBROUTINE energy_profile
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE soil_absorption(PXA,PXB,PKESTAR,PDTAUSTAR,PTAUSTAR,PGM,PCOSZEN,PALB,KNLVLS_EFF,PSOILENERGY)
!compute the energy absorbed by the soil at each wavelength
!
USE modd_const_tartes, ONLY : npnbands
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pxa,pxb ! solutions of linear system (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pkestar !Asymptotic Flux Extinction Coefficent (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pdtaustar !Optical depth of each layer (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: ptaustar !Cumulated optical depth (npoints,nlayer,nbands)
REAL, DIMENSION(:,:,:), INTENT(IN)  :: pgm !GM vector (npoints,nlayer,nbands)
REAL, DIMENSION(:), INTENT(IN)      :: pcoszen! cosine of zenithal solar angle (npoints)
REAL, DIMENSION(:,:), INTENT(IN)    :: palb! soil albedo (npoints,nbands)
INTEGER, DIMENSION(:,:), INTENT(IN) :: knlvls_eff !number of effective layers (npoints,nbands)
REAL, DIMENSION(:,:), INTENT(OUT)   :: psoilenergy
INTEGER :: jb,ji !loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SOIL_ABSORPTION',0,zhook_handle)
!
DO jb=1,npnbands
  !
  DO ji=1,SIZE(pxa,1)
    !
    !doc equation 65
    psoilenergy(ji,jb) = ( 1.-palb(ji,jb) ) * &
                           ( pxa(ji,knlvls_eff(ji,jb),jb) * &
                               exp( -pkestar(ji,knlvls_eff(ji,jb),jb) * pdtaustar(ji,knlvls_eff(ji,jb),jb) ) + &
                             pxb(ji,knlvls_eff(ji,jb),jb) * &
                               exp(  pkestar(ji,knlvls_eff(ji,jb),jb) * pdtaustar(ji,knlvls_eff(ji,jb),jb) ) + &
                           ( pgm(ji,knlvls_eff(ji,jb),jb)+pcoszen(ji) ) * &
                               exp( -ptaustar(ji,knlvls_eff(ji,jb),jb)/pcoszen(ji) ) )
    !
  END DO
  !
END DO

IF (lhook) CALL dr_hook('SOIL_ABSORPTION',1,zhook_handle)

END SUBROUTINE soil_absorption
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE spectral_repartition(PSW_RAD,PCOSZEN,PSW_RAD_DIF,PSW_RAD_DIR,PNIR_ABS)

USE modd_const_tartes, ONLY : npnbands,xpratio_dir,xpratio_dif,xpcoefnir_dir,xpcoefnir_dif,xp_mudiff

IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN)    :: psw_rad ! broadband global incident light (W/m^2) (npoints)
REAL, DIMENSION(:), INTENT(IN)    :: pcoszen ! cosine of zenithal solar angle (npoints)
REAL, DIMENSION(:,:), INTENT(OUT) :: psw_rad_dif ! spectral diffuse incident light (W/m^2) (npoints,nbands)
REAL, DIMENSION(:,:), INTENT(OUT) :: psw_rad_dir ! spectral direct incident light (W/m^2) (npoints,nbands)
REAL, DIMENSION(:), INTENT(OUT)   :: pnir_abs ! Near infrared radiation (2500-4000 nm) absorbed by snowpack (W/m^2) (npoints)
!
REAL, DIMENSION(SIZE(PSW_RAD)) :: zsw_rad_broaddir,zsw_rad_broaddif ! direct and diffuse broadband incident light (W/m^2) (npoints)
INTEGER :: jb !Loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SPECTRAL_REPARTITION',0,zhook_handle)
!
! Separate broadband global radiation in direct and diffuse (parametrization Marie Dumont)
! NB : thresold 1. to the factor when zenithal angle close to pi/2
zsw_rad_broaddif = min( exp( - 1.54991930344*pcoszen**3 + 3.73535795329*pcoszen**2 &
                             - 3.52421131883*pcoszen + 0.0299111951172 ), 1. ) * psw_rad
zsw_rad_broaddir = psw_rad - zsw_rad_broaddif
!
! Spectral decomposition
DO jb = 1,npnbands
  psw_rad_dif(:,jb) = xpratio_dif(jb) * zsw_rad_broaddif / xp_mudiff
  psw_rad_dir(:,jb) = xpratio_dir(jb) * zsw_rad_broaddir / pcoszen(:)
END DO

pnir_abs = zsw_rad_broaddif*xpcoefnir_dif + zsw_rad_broaddir*xpcoefnir_dir

IF (lhook) CALL dr_hook('SPECTRAL_REPARTITION',1,zhook_handle)

END SUBROUTINE spectral_repartition
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
SUBROUTINE snowcro_tartes(PSNOWGRAN1,PSNOWGRAN2,PSNOWRHO,PSNOWDZ,PSNOWG0,PSNOWY0,PSNOWW0,PSNOWB0, &
                          psnowimp_density,psnowimp_content,palb,psw_rad,pzenith,knlvls_use,      &
                          psnowalb,pradsink,pradxs,odebug,hsnowmetamo)
!
! Interface between Tartes and Crocus
! M. Lafaysse 26/08/2013
!
USE modd_const_tartes, ONLY : npnimp
USE modd_snow_metamo,  ONLY : xuepsi
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowgran1,psnowgran2  ! (npoints,nlayer) 
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowrho !snow density (kg/m^3) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowg0 ! asymmetry parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowy0 ! Value of y of snow grains at nr=1.3 (no unit
REAL, DIMENSION(:,:), INTENT(IN)   :: psnoww0 ! Value of W of snow grains at nr=1.3 (no unit)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowb0 ! absorption enhancement parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowdz !snow layers thickness (m) (npoints,nlayer)
REAL, DIMENSION(:,:,:), INTENT(IN) :: psnowimp_density !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(:,:,:), INTENT(IN) :: psnowimp_content !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
!
REAL, DIMENSION(:), INTENT(IN)     :: palb ! soil/vegetation albedo (npoints)
!
REAL, DIMENSION(:), INTENT(IN)     :: psw_rad ! global broadband incident light (W/m^2) (npoints)
REAL, DIMENSION(:), INTENT(IN)     :: pzenith ! zenithal solar angle (npoints)
!
INTEGER, DIMENSION(:), INTENT(IN)  :: knlvls_use ! number of effective snow layers (npoints)
!
!Same outputs as SNOWCRORAD and SNOWCROALB
REAL, DIMENSION(:,:), INTENT(OUT) :: pradsink !(npoints,nlayers)
REAL, DIMENSION(:), INTENT(OUT)   :: pradxs !(npoints,nlayers)
REAL, DIMENSION(:), INTENT(OUT)   :: psnowalb !(npoints,nlayers)
!
LOGICAL, INTENT(IN) :: odebug ! Print for debugging
 CHARACTER(3), INTENT(IN)          :: hsnowmetamo ! metamorphism scheme
!
!packed variables
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),NPNIMP) :: zsnowimp_density_p !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),NPNIMP) :: zsnowimp_content_p !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowgran1_p,zsnowgran2_p  ! (npoints,nlayer) 
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrho_p !snow density (kg/m^3) (npoints,nlayer)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowg0_p ! asymmetry parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit) (npoints,nlayer)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowy0_p ! Value of y of snow grains at nr=1.3 (no unit
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnoww0_p ! Value of W of snow grains at nr=1.3 (no unit)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowb0_p ! absorption enhancement parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowdz_p !snow layers thickness (m) (npoints,nlayer)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zradsink_p
!
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zalb_p ! soil/vegetation albedo (npoints)
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zsw_rad_p ! global broadband incident light (W/m^2) (npoints)
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zzenith_p ! zenithal solar angle (npoints)
!
!Same outputs as SNOWCRORAD and SNOWCROALB
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zradxs_p
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zsnowalb_p
!
INTEGER, DIMENSION(SIZE(PSNOWRHO,1)) :: inlvls_use_p ! number of effective snow layers (npoints)
!
INTEGER, DIMENSION(SIZE(PSNOWRHO,1)) :: idaymask ! mask for points where it's day
!
INTEGER :: imax_use ! maximum number of layers over the domain
INTEGER :: jl,jimp,jj,jj_p !Loop counter
INTEGER :: ipointday
INTEGER :: inpoints
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO_TARTES',0,zhook_handle)
!
!Default values (night)
pradsink = 0.
pradxs   = 0.
psnowalb = 1.
!
inpoints = SIZE(psnowrho,1)
!
!Mask
ipointday = 0
DO jj = 1,inpoints
  IF ( cos(pzenith(jj))>xuepsi .AND. psw_rad(jj)>xuepsi ) THEN
    !mask for day
    ipointday = ipointday + 1 
    idaymask(ipointday) = jj
  END IF
END DO
!
IF ( ipointday>=1 ) THEN
  !
  ! Pack 1D variables
  DO jj_p = 1,ipointday
    !
    jj = idaymask(jj_p)
    !
    zalb_p(jj_p) = palb(jj)
    zsw_rad_p(jj_p) = psw_rad(jj)
    zzenith_p(jj_p) = pzenith(jj)
    inlvls_use_p(jj_p) = knlvls_use(jj)
    !
  END DO
  !
  imax_use = maxval(knlvls_use)
  !
  ! Pack 2D variables
  DO jl = 1,imax_use
    !
    DO jj_p = 1,ipointday
      !
      jj = idaymask(jj_p)
      !
      zsnowgran1_p(jj_p,jl) = psnowgran1(jj,jl)
      zsnowgran2_p(jj_p,jl) = psnowgran2(jj,jl)
      zsnowrho_p(jj_p,jl) = psnowrho(jj,jl)
      zsnowg0_p(jj_p,jl) = psnowg0(jj,jl)
      zsnowy0_p(jj_p,jl) = psnowy0(jj,jl)
      zsnoww0_p(jj_p,jl) = psnoww0(jj,jl)
      zsnowb0_p(jj_p,jl) = psnowb0(jj,jl)      
      zsnowdz_p(jj_p,jl) = psnowdz(jj,jl)
      !
    END DO
    !
  END DO
  !
  ! Pack 3D variables
  DO jimp = 1,npnimp
    !
    DO jl = 1,imax_use
      !
      DO jj_p = 1,ipointday
        !
        jj = idaymask(jj_p)
        !
        zsnowimp_density_p(jj_p,jl,jimp) = psnowimp_density(jj,jl,jimp)
        zsnowimp_content_p(jj_p,jl,jimp) = psnowimp_content(jj,jl,jimp)      
        !
      END DO
      !
    END DO
    !
  END DO
  !
!RJ: fix fp-invalid(nan) trapping in ISBA_DIF8_SN3L_NIT_SNCRO8_C13_SNOWRAD_TAR, ISBA_DIF8_SN3L_NIT_SNCRO8_C13_SNOWRAD_TA2 tests
#ifdef RJ_OFIX
!RJ: temp fix to avoid accessing uninited values (NANS) in mode_tartes.F90 by explicit inited shape, no inpact for results, problem in array padding
  CALL snowcro_call_tartes(zsnowgran1_p(1:ipointday,1:imax_use),zsnowgran2_p(1:ipointday,1:imax_use), &
                           zsnowrho_p(1:ipointday,1:imax_use),zsnowdz_p(1:ipointday,1:imax_use),      &
                           zsnowg0_p(1:ipointday,1:imax_use),zsnowy0_p(1:ipointday,1:imax_use),       &
                           zsnoww0_p(1:ipointday,1:imax_use),zsnowb0_p(1:ipointday,1:imax_use),       &
                           zsnowimp_density_p(1:ipointday,1:imax_use,1:npnimp),                       &
                           zsnowimp_content_p(1:ipointday,1:imax_use,1:npnimp),                       &
                           zalb_p(1:ipointday),zsw_rad_p(1:ipointday),                                &
                           zzenith_p(1:ipointday),inlvls_use_p(1:ipointday),zsnowalb_p(1:ipointday),  & 
                           zradsink_p(1:ipointday,1:imax_use),zradxs_p(1:ipointday),odebug,hsnowmetamo)
#else
  CALL snowcro_call_tartes(zsnowgran1_p(1:ipointday,:),zsnowgran2_p(1:ipointday,:),zsnowrho_p(1:ipointday,:),     &
                           zsnowdz_p(1:ipointday,:),zsnowg0_p(1:ipointday,:),zsnowy0_p(1:ipointday,:),            &
                           zsnoww0_p(1:ipointday,:),zsnowb0_p(1:ipointday,:),zsnowimp_density_p(1:ipointday,:,:), &
                           zsnowimp_content_p(1:ipointday,:,:),zalb_p(1:ipointday),zsw_rad_p(1:ipointday),        &
                           zzenith_p(1:ipointday),inlvls_use_p(1:ipointday),zsnowalb_p(1:ipointday),              &
                           zradsink_p(1:ipointday,:),zradxs_p(1:ipointday),odebug,hsnowmetamo)
#endif
  !
  !Unpack 1d output variables
  !
  DO jj_p = 1,ipointday
    !
    jj = idaymask(jj_p)
    !
    pradxs(jj) = zradxs_p(jj_p)
    psnowalb(jj) = zsnowalb_p(jj_p)
    !
  END DO
  !
  !Unpack 2d output  variables
  DO jl = 1,imax_use
    !
    DO jj_p = 1,ipointday
      !
      jj = idaymask(jj_p)
      !
      pradsink(jj,jl) = zradsink_p(jj_p,jl)
      !
    END DO
    !
  END DO
  !
END IF
!
IF (lhook) CALL dr_hook('SNOWCRO_TARTES',1,zhook_handle)
!
END SUBROUTINE snowcro_tartes

!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------

SUBROUTINE snowcro_call_tartes(PSNOWGRAN1,PSNOWGRAN2,PSNOWRHO,PSNOWDZ,PSNOWG0,PSNOWY0,PSNOWW0,PSNOWB0, &
                               psnowimp_density,psnowimp_content,palb,psw_rad,pzenith,knlvls_use,      &
                               psnowalb,pradsink,pradxs,odebug,hsnowmetamo)
!
! Interface between Tartes and Crocus
! M. Lafaysse 26/08/2013
!
USE modd_const_tartes, ONLY : npnbands,xpwavelengths,xp_mudiff
USE modd_csts, ONLY : xrholi,xpi
!
USE mode_snow3l, ONLY : get_diam
!
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN)   :: psnowgran1,psnowgran2  ! (npoints,nlayer) 
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowrho !snow density (kg/m^3) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowg0 ! asymmetry parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit) (npoints,nlayer)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowy0 ! Value of y of snow grains at nr=1.3 (no unit
REAL, DIMENSION(:,:), INTENT(IN)   :: psnoww0 ! Value of W of snow grains at nr=1.3 (no unit)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowb0 ! absorption enhancement parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit)
REAL, DIMENSION(:,:), INTENT(IN)   :: psnowdz !snow layers thickness (m) (npoints,nlayer)
REAL, DIMENSION(:,:,:), INTENT(IN) :: psnowimp_density !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(:,:,:), INTENT(IN) :: psnowimp_content !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
!
REAL, DIMENSION(:), INTENT(IN)     :: palb ! soil/vegetation albedo (npoints)
!
REAL, DIMENSION(:), INTENT(IN)     :: psw_rad ! global broadband incident light (W/m^2) (npoints)
REAL, DIMENSION(:), INTENT(IN)     :: pzenith ! zenithal solar angle (npoints)
!
INTEGER, DIMENSION(:), INTENT(IN)  :: knlvls_use ! number of effective snow layers (npoints)
!
!Same outputs as SNOWCRORAD and SNOWCROALB
REAL, DIMENSION(:,:), INTENT(OUT) :: pradsink !(npoints,nlayers)
REAL, DIMENSION(:), INTENT(OUT)   :: pradxs !(npoints,nlayers)
REAL, DIMENSION(:), INTENT(OUT)   :: psnowalb !(npoints,nlayers)

LOGICAL,INTENT(IN) :: odebug ! Print for debugging
 CHARACTER(3), INTENT(IN)          :: hsnowmetamo ! metamorphism scheme

!Local variables
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),NPNBANDS) :: zsnowenergy !(npoints,nlayer,nbands)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowssa !snow specific surface area (m^2/kg) (npoints,nlayer) 
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowenergy_bb ! (W/m^2) (npoints,nlayers)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),NPNBANDS) :: zsw_rad_dif ! spectral diffuse incident light (W/m^2) (npoints,nbands)
REAL, DIMENSION(SIZE(PSNOWRHO,1),NPNBANDS) :: zsw_rad_dir ! spectral direct incident light (W/m^2) (npoints,nbands)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),NPNBANDS) :: zsnowalb  !(npoints,nbands)
REAL, DIMENSION(SIZE(PSNOWRHO,1),NPNBANDS) :: zalb ! soil/vegetation albedo (npoints,nbands)
REAL, DIMENSION(SIZE(PSNOWRHO,1),NPNBANDS) :: zsoilenergy !(npoints,nbands)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: znir_abs ! Near infrared radiation (2500-4000 nm) absorbed by snowpack (W/m^2) (npoints)
!
!broad band
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ztotsnowenergy ! (W/m^2) (npoints)
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zsoilenergy_bb ! (W/m^2) (npoints)
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zreflected_bb ! (W/m^2) (npoints)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zsnowenergy_cum,zsnowenergy_upper ! Cumulated absorbed energy for 1 wavelength W/m^2 (npoints)
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zmax ! maximum available energy for 1 wavelength  W/m^2 (npoints,nbands)
!
REAL :: zdiam !optical diameter
!
INTEGER :: jb,jl,jj !Loop counter
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO_CALL_TARTES',0,zhook_handle)
!
! Compute SSA from SNOWGRAN1 and SNOWGRAN2
DO jl = 1,SIZE(psnowrho,2)
  !
  DO jj = 1,SIZE(psnowrho,1)
    !
    IF ( jl<=knlvls_use(jj) ) THEN
      !
      CALL get_diam(psnowgran1(jj,jl),psnowgran2(jj,jl),zdiam,hsnowmetamo)
      zsnowssa(jj,jl) = 6. / (xrholi*zdiam)
      !
    ENDIF
    !
  ENDDO
  !
ENDDO
!
DO jb = 1,npnbands
  ! Soil-vegetation albedo homogeneous in wavelength
  zalb(:,jb) = palb(:)
END DO
!
!Spectral repartition of radiation
 CALL spectral_repartition(psw_rad,cos(pzenith),zsw_rad_dif,zsw_rad_dir,znir_abs)
!
IF ( odebug ) THEN
  WRITE(*,*) "ZSW_RAD_DIF=",zsw_rad_dif
  WRITE(*,*) "ZSW_RAD_DIR=",zsw_rad_dir
  WRITE(*,*) "PZENITH=",pzenith
END IF
!
!Call tartes model
! For test and debugging this routine can be called independently by a python interface   
!
 CALL tartes(zsnowssa,psnowrho,psnowdz,psnowg0,psnowy0,psnoww0,psnowb0,psnowimp_density,psnowimp_content,zalb,&
             zsw_rad_dif,zsw_rad_dir,cos(pzenith),knlvls_use,zsnowalb,zsnowenergy,zsoilenergy)
! 
    ! Modif ML : in some cases, Tartes is unstable in infra-red wavelengths : control of energy values and apply threshold if necessary
    ! --------------------------------------------------------------------------------------------------------------------
    
    
    ! This does not seem to be necessary at CDP when the effective number of layers is properly limited.
    ! However, we let the comment code because it might happen again.
    
!     
!     DO JB=1,NPNBANDS
!       ! maximum available energy at this wavelength
!       ZMAX=ZSW_RAD_DIF(:,JB)*PP_MUDIFF+ZSW_RAD_DIR(:,JB)*COS(PZENITH)
!       ZSNOWENERGY_CUM(:)=0.
!       DO JL=1,SIZE(PSNOWRHO,2)
!         DO JJ=1,SIZE(PSNOWRHO,1)        
!           IF (JL<=KNLVLS_USE(JJ)) THEN
!             ZSNOWENERGY_UPPER(JJ)=ZSNOWENERGY_CUM(JJ) !0 for surface layer
!             ! absorbed energy cumulated from the surface
!             ZSNOWENERGY_CUM(JJ)=ZSNOWENERGY_CUM(JJ)+ZSNOWENERGY(JJ,JL,JB)
!             
!             ! Case when the energy is negative but close to 0. It can occurs in short wavelengths. Force to 0.
!             IF ((ZSNOWENERGY(JJ,JL,JB)<0).AND. (ZSNOWENERGY(JJ,JL,JB)>-0.1)) THEN
!               ZSNOWENERGY(JJ,JL,JB)=0.
!             END IF
!             
!             !if the cumulated absorbed energy excess the available energy or if severe negative energy, numerical problem in Tartes : total absorption
!             IF ((ZSNOWENERGY_CUM(JJ)>ZMAX(JJ)).OR.(ZSNOWENERGY(JJ,JL,JB)<=-0.1)) THEN
! !               IF (PPWAVELENGTHS(JB)<=1000) THEN             
!               IF ((PPWAVELENGTHS(JB)<=1000) .AND.(ABS(ZSNOWENERGY_CUM(JJ)-ZMAX(JJ))>0.01)) THEN 
!                 ! Tolerance 0.01 W/m2 of excess energy in the visible
!                 ! Above, the problem should never happen in visible
!                 
!                 PRINT*,"JB=",JB,"JL=",JL
!                 IF (ZSNOWENERGY_CUM(JJ)>ZMAX(JJ)) PRINT*,"excess energy"
!                 IF (ZSNOWENERGY(JJ,JL,JB)<0) PRINT*, "negative energy"
! 
!                 PRINT*,ZSW_RAD_DIF(JJ,JB),ZSW_RAD_DIR(JJ,JB),ZMAX(JJ),ZSNOWENERGY_CUM(JJ)
!                 
!                 PRINT*,"profile :",ZSNOWENERGY(JJ,:,JB)
!                 
!                 STOP "FATAL ERROR TARTES !!"
!               END IF
!               ! The layer absorbes all the remaining energy at this wavelength
!               ZSNOWENERGY(JJ,JL,JB)=ZMAX(JJ)-ZSNOWENERGY_UPPER(JJ)
!               ! update cumulated energy
!               ZSNOWENERGY_CUM(JJ)=ZMAX(JJ)
!             END IF
!           END IF
!         END DO
!       END DO
!       
!       ! Threshold on soil absorbed energy
!       DO JJ=1,SIZE(PSNOWRHO,1)
!         ZSNOWENERGY_UPPER(JJ)=ZSNOWENERGY_CUM(JJ)
!         ZSNOWENERGY_CUM(JJ)=ZSNOWENERGY_CUM(JJ)+ZSOILENERGY(JJ,JB)
!         IF (ZSNOWENERGY_CUM(JJ)>ZMAX(JJ)) THEN
!           IF (PPWAVELENGTHS(JB)<=1200) THEN
!                  STOP "FATAL ERROR TARTES (soil excess energy in visible)!!"           
!           END IF
!           ! The layer absorbes all the remaining energy at this wavelength
!           ZSOILENERGY(JJ,JB)=ZMAX(JJ)-ZSNOWENERGY_UPPER(JJ)
!         END IF
!       END DO
!       
!     END DO

    ! End modif ML
    ! --------------------------------------------------------------------------------------------------------------------
!    
! Broadband absorbed energy by snowpack and soil
zsnowenergy_bb = 0.
zsoilenergy_bb = 0.
!
DO jb = 1,npnbands
  zsnowenergy_bb(:,:) = zsnowenergy_bb(:,:) + zsnowenergy(:,:,jb)
  zsoilenergy_bb(:)   = zsoilenergy_bb(:)   + zsoilenergy(:,jb)
END DO
!
!Add near infra-red wavelengths to first layer
zsnowenergy_bb(:,1) = zsnowenergy_bb(:,1) + znir_abs
!
! Total energy absorbed by snowpack
ztotsnowenergy(:)=0
DO jl = 1,SIZE(psnowrho,2)
  DO jj = 1,SIZE(psnowrho,1)
    IF ( jl<=knlvls_use(jj) ) THEN
      ztotsnowenergy(jj) = ztotsnowenergy(jj) + zsnowenergy_bb(jj,jl)
    END IF
  END DO
END DO
!    
! Reflected energy
zreflected_bb = psw_rad - ztotsnowenergy - zsoilenergy_bb
! 
! Broad band Albedo
! PSW_RAD is never 0 because this routine is not called during the night
psnowalb = zreflected_bb / psw_rad
!   
! Source term
pradsink(:,1) = -psw_rad(:) + zreflected_bb + zsnowenergy_bb(:,1)
!    
DO jl = 2,SIZE(psnowrho,2)
  pradsink(:,jl) = pradsink(:,jl-1) + zsnowenergy_bb(:,jl)
END DO
!
!Excess energy
pradxs = psw_rad - ztotsnowenergy - zreflected_bb
!
IF (lhook) CALL dr_hook('SNOWCRO_CALL_TARTES',1,zhook_handle)
!
END SUBROUTINE snowcro_call_tartes       

END MODULE mode_tartes