DevelopmentActivities/MergeHydro/Martial_notes_on_merge: correction_neige

 

SUBROUTINE thermosoil_humlev(kjpindex, shumdiag, snow)

    ! input scalar 
    INTEGER(i_std), INTENT(in)                            :: kjpindex      !! Domain size
    ! input fields
!Isa
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: snow
    REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in)    :: shumdiag      !! Diagnostoc profile
    !
    ! modified fields
    !
    ! output fields

    !
    ! local variable
    !
    INTEGER(i_std)  :: ji, jd, jg
    REAL(r_std)  :: lev_diag, prev_diag, lev_prog, prev_prog
    REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfactw
!Isa
    REAL(r_std), DIMENSION(kjpindex)                    :: snow_h
    !
    !
    IF ( .NOT. ALLOCATED(intfactw)) THEN
        ALLOCATE(intfactw(ngrnd, nbdl))
   ENDIF
!Isa
snow_h(:)=snow(:)/sn_dens

do ji=1,kjpindex
        prev_diag = 0.0
        DO jd = 1, ngrnd
           lev_diag = prev_diag + dz2(jd)
           prev_prog = 0.0+snow_h(ji)
           DO jg = 1, nbdl
                IF ( jg == nbdl .AND. diaglev(jg)+snow_h(ji) < lev_diag ) THEN
                   !! Just make sure we cover the deepest layers
                   lev_prog = lev_diag+snow_h(ji)
                ELSE
                   lev_prog = diaglev(jg)+snow_h(ji)
                ENDIF
                intfactw(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), 0.0)/(lev_diag-prev_diag)
                prev_prog = lev_prog
          ENDDO
          prev_diag = lev_diag
        ENDDO

    DO jd = 1, nbdl
      DO jg = 1, ngrnd
            wetdiag(ji,jg) = wetdiag(ji,jg) + shumdiag(jd)*intfactw(jg,jd)
        ENDDO
    ENDDO
enddo
END SUBROUTINE thermosoil_humlev