snowcro_diag.F90

Go to the documentation of this file.

!     #########

SUBROUTINE snowcro_diag(HSNOWMETAMO, &
                        PSNOWDZ, PSNOWSWE, PSNOWRHO, PSNOWGRAN1, PSNOWGRAN2, PSNOWAGE, &
                        PSNOWHIST, PSNOWTEMP, PSNOWLIQ, PDIRCOSZW, PSNOWDEND, PSNOWSPHER, &
                        PSNOWSIZE, PSNOWSSA, PSNOWTYPEMEPRA, PSNOWRAM, PSNOWSHEAR, &
                        PSNOWDEPTH_1DAYS, PSNOWDEPTH_3DAYS, PSNOWDEPTH_5DAYS, &
                        PSNOWDEPTH_7DAYS, PSNOWSWE_1DAYS, PSNOWSWE_3DAYS, PSNOWSWE_5DAYS,&
                        PSNOWSWE_7DAYS, PSNOWRAM_SONDE, PSNOW_WETTHICKNESS, PSNOW_REFROZENTHICKNESS)

! Diagnostics of Crocus snowpack model
! Author: M. Lafaysse, Meteo-France, October 2015

USE modd_surf_par,      ONLY : xundef

USE modd_csts,ONLY : xrholi, xrholw

USE modd_snow_par,ONLY : icris_dend1d, icris_nondend1d, &
                         ifr, ifr_lb, ilb, ilb_fin, ilb_ang, iroul, ifin, ifin_ar, ifin_ang, &
                         ipl, ipl_gob, igob, igel, igob_fon, iron_ang, xx, xd1, xd2, xd3

IMPLICIT NONE

 CHARACTER(3), INTENT(IN)         :: HSNOWMETAMO ! metamorphism option
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWSWE
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWGRAN1
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWGRAN2
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWAGE
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWHIST
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP
REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWLIQ
REAL, DIMENSION(:),   INTENT(IN) :: PDIRCOSZW !cosine of slope
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWDEND
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWSPHER
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWSIZE
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWSSA
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWTYPEMEPRA
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWRAM
REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWSHEAR
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWDEPTH_1DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWDEPTH_3DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWDEPTH_5DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWDEPTH_7DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWSWE_1DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWSWE_3DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWSWE_5DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWSWE_7DAYS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOWRAM_SONDE
REAL, DIMENSION(:), INTENT(OUT) :: PSNOW_WETTHICKNESS
REAL, DIMENSION(:), INTENT(OUT) :: PSNOW_REFROZENTHICKNESS

REAL :: ZG1, ZG2, ZRFIN, ZRDEN, ZRFGF

REAL :: ZDIAM

LOGICAL,DIMENSION(SIZE(PSNOWSWE,1)) :: GRAM, GWET, GREFROZEN

INTEGER :: ICLASS_DEND, ICLASS_SPHER, ICLASS_SIZE, ICLASS_HIST
INTEGER :: ICLASS
LOGICAL :: LTHERM

INTEGER :: JJ,JST

! PRINT*,ASSOCIATED(PSNOWDEPTH_1DAYS),ASSOCIATED(PSNOWDEPTH_3DAYS)
! PRINT*,ASSOCIATED(PSNOWDEPTH_5DAYS),ASSOCIATED(PSNOWDEPTH_7DAYS)
! PRINT*,"in snowcrodiag"
! PRINT*,ALLOCATED(PSNOWDEPTH_1DAYS),ALLOCATED(PSNOWDEPTH_3DAYS)
! PRINT*,ALLOCATED(PSNOWDEPTH_5DAYS),ALLOCATED(PSNOWDEPTH_7DAYS)
! PRINT*,SIZE(PSNOWDEPTH_1DAYS),SIZE(PSNOWDEPTH_3DAYS)
! PRINT*,SIZE(PSNOWDEPTH_5DAYS),SIZE(PSNOWDEPTH_7DAYS)

! Initializations

psnowdend      = xundef
psnowspher     = xundef
psnowsize      = xundef
psnowssa       = xundef
psnowtypemepra = xundef
psnowram       = xundef
psnowshear     = xundef

psnowdepth_1days = 0.
psnowdepth_3days = 0.
psnowdepth_5days = 0.
psnowdepth_7days = 0.
psnowswe_1days = 0.
psnowswe_3days = 0.
psnowswe_5days = 0.
psnowswe_7days = 0.
psnowram_sonde = 0.
psnow_wetthickness = 0.
psnow_refrozenthickness = 0.
gram = .true.
gwet = .true.
grefrozen = .true.


DO jst=1,SIZE(psnowswe,2)
  DO jj=1,SIZE(psnowswe,1)
    
    IF (psnowswe(jj,jst)>0) THEN
    
      ! In this routine Crocus diagnostics are perpendicular to the slope.
      ! The projection is done in diag_misc_isban             
    
      zg1 = psnowgran1(jj,jst)/99.
      zrfin = 0.17*psnowrho(jj,jst)-31
    
      IF (psnowgran1(jj,jst)>=0) THEN
        !Non dendritic case

        !Dendricity,sphericty and grain size
        psnowsize(jj,jst)  = psnowgran2(jj,jst)
        psnowdend(jj,jst)  = 0
        psnowspher(jj,jst) = psnowgran1(jj,jst) / xx

        !Optical diameter for SSA diagnostic
        zdiam = psnowsize(jj,jst) * psnowspher(jj,jst) + &
        max( 0.0004, 0.5*psnowsize(jj,jst) ) * ( 1.-psnowspher(jj,jst) )

        !10 classes of sphericity 0:[0,0.05[, 1:[0.05,0.15[, ..., 9:[0.85,1.0]
        !###########Strange way of defining sphericity classes -> Check with very old versions
        !###########ICLASS_SPHER = MIN(NINT(10 * PSNOWSPHER(JJ,JST)),9) Not exactly the same
        iclass_spher = min(int(10 * psnowspher(jj,jst) + 0.05),9)


        !6 classes of historical variable {0,1,...,5}
        !##########Why is PSNOWHIST stored as float64. It always takes only 6 integer values
        iclass_hist = nint(psnowhist(jj,jst))

        !3 classes of grain size in mm 0:[0,0.55[, 1:[0.55,1.05[, 2:[1.05, +inf[
        !#########Strange +0.05
        IF (psnowsize(jj,jst) < 0.00055) THEN
          iclass_size = 0
        ELSEIF (psnowsize(jj,jst) < 0.00105) THEN
          iclass_size = 1
        ELSE
          iclass_size = 2
        ENDIF

        !Overall 10x3x6 classes from 1 to 180 (included)
        iclass = 1 + iclass_spher + iclass_size*10 + iclass_hist*30

        !Snow type obtained in table ICRIS_NONDEND1D
        psnowtypemepra(jj,jst) = icris_nondend1d(iclass)


        ! Ram resistance (non dendritic case)
        ! PSNOWLIQ tel massique en m
        ! PSNOWLIQ*XHROLW/PSNOWDZ tel volumique en kg m-3
        ! PSNOWLIQ/PSNOWDZ : rapport sans unité
        ! Seuil à 0.5% soit 0.005
        ltherm=((psnowtemp(jj,jst)<272.96).OR.(psnowliq(jj,jst)/psnowdz(jj,jst)<=0.005))

        SELECT CASE (nint(psnowtypemepra(jj,jst)))
        CASE (ifin)
          psnowram(jj,jst)=max(3.,zrfin)
        CASE (ifin_ang)
          IF (psnowrho(jj,jst)<200) THEN
            psnowram(jj,jst)=zrfin*psnowspher(jj,jst)+&
            (1- psnowspher(jj,jst))*(zrfin*(0.8-psnowsize(jj,jst))+2*psnowsize(jj,jst))
          ELSE
            psnowram(jj,jst)=2
          ENDIF
        CASE (ifin_ar,igel,igob_fon,iron_ang)
      
          IF (ltherm) THEN
            psnowram(jj,jst)=max(10.,0.103*psnowrho(jj,jst)-19.666)
          ELSE
            IF (psnowrho(jj,jst)<250) THEN
              psnowram(jj,jst)=1
            ELSE
              psnowram(jj,jst)=max(2.,0.16*psnowrho(jj,jst)-54)
            ENDIF
          END IF
      
        CASE (ipl,ipl_gob)
          IF (psnowsize(jj,jst)>0.8) THEN
            psnowram(jj,jst)=max(3.,zrfin)*(0.8-psnowsize(jj,jst))+2*psnowsize(jj,jst)
          ELSE
            psnowram(jj,jst)=2
          ENDIF
        CASE DEFAULT
        END SELECT

      ELSE
        !Dendritic case

        !Dendricity,sphericty and grain size
        psnowsize(jj,jst)  =  xundef  !Grain size not defined for dendritic snow
        psnowdend(jj,jst)  = -psnowgran1(jj,jst) / xx
        psnowspher(jj,jst) =  psnowgran2(jj,jst) / xx

        !Optical diameter for SSA diagnostic
        zdiam = psnowdend(jj,jst) * xd1 + (1 - psnowdend(jj,jst)) * &
        (psnowspher(jj,jst) * xd2 + (1 - psnowspher(jj,jst)) * xd3)
        !ZDIAM =  -PSNOWGRAN1(JJ,JST)*XD1/XX + (1.+PSNOWGRAN1(JJ,JST)/XX) * &
        !      ( PSNOWGRAN2(JJ,JST)*XD2/XX + (1.-PSNOWGRAN2(JJ,JST)/XX) * XD3 )
        zdiam = zdiam/10000.

        !10 classes of dendricity 0:[0,0.1[, ..., 9:[0.9,1.0[ (value 1.0 does not exist)
        iclass_dend = int(10 * psnowdend(jj,jst))

        !10 classes of sphericity 0:[0,0.05[, 1:[0.05,0.15[, ..., 9:[0.85,1.0]
        !###########Strange way of defining sphericity classes -> Check with very old versions
        !###########ICLASS_SPHER = MIN(NINT(10 * PSNOWSPHER(JJ,JST)),9) Not exactly the same
        iclass_spher = min(int(10 * psnowspher(jj,jst) + 0.05),9)

        !Overall 10x10 classes from 1 to 100 (included)
        iclass = 1 + iclass_dend + iclass_spher*10

        !Snow type obtained in table ICRIS_DEND1D
        psnowtypemepra(jj,jst) = icris_dend1d(iclass)

      ENDIF

      ! All cases
      ! Compute depth and SWE of recent snow
      IF(psnowage(jj,jst)<=1)THEN
        psnowdepth_1days(jj) = psnowdepth_1days(jj) + psnowdz(jj,jst)
        psnowswe_1days(jj) = psnowswe_1days(jj) + psnowswe(jj,jst)
      ENDIF

      IF(psnowage(jj,jst)<=3)THEN
        psnowdepth_3days(jj) = psnowdepth_3days(jj) + psnowdz(jj,jst)    
        psnowswe_3days(jj) = psnowswe_3days(jj) + psnowswe(jj,jst)
      ENDIF
    
      IF(psnowage(jj,jst)<=5)THEN
        psnowdepth_5days(jj) = psnowdepth_5days(jj) + psnowdz(jj,jst)    
        psnowswe_5days(jj) = psnowswe_5days(jj) + psnowswe(jj,jst)
      ENDIF
  
      IF(psnowage(jj,jst)<=7)THEN
        psnowdepth_7days(jj) = psnowdepth_7days(jj) + psnowdz(jj,jst)    
        psnowswe_7days(jj) = psnowswe_7days(jj) + psnowswe(jj,jst)
      END IF
    
      ! Ram sonde penetration
      IF ((gram(jj)).AND.(psnowram(jj,jst)<=2.)) THEN
        psnowram_sonde(jj)=psnowram_sonde(jj)+psnowdz(jj,jst)
      ELSE
        gram(jj)=.false.
      ENDIF

      ! Depth of wet snow
      IF ((gwet(jj)).AND.(psnowliq(jj,jst)>0.)) THEN
        psnow_wetthickness(jj)=psnow_wetthickness(jj)+psnowdz(jj,jst)
      ELSE
        gwet(jj)=.false.
      ENDIF
      ! Depth of refrozen snow
      IF ((grefrozen(jj)).AND.(psnowhist(jj,jst)>=2).AND.(psnowtemp(jj,jst)<273.15)) THEN
        psnow_refrozenthickness(jj)=psnow_refrozenthickness(jj)+psnowdz(jj,jst)
      ELSE
        grefrozen(jj)=.false.
      ENDIF    
    
      ! Specific surface area
      IF ( hsnowmetamo=='B92' ) THEN
        psnowssa(jj,jst) = 6. / (xrholi*zdiam)
      ELSE
        psnowssa(jj,jst) = 6. / (xrholi*psnowgran1(jj,jst))
      END IF

    END IF
  END DO
END DO
  
END SUBROUTINE snowcro_diag