read_gr_snow.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.
!     #########
      SUBROUTINE read_gr_snow(HPROGRAM, HSURFTYPE, HPREFIX,     &
                              KLU, KSIZE_P, KMASK_P, KPATCH, TPSNOW, HDIR, KVERSION, KBUGFIX, KNPATCH)  
!     ##########################################################
!
!!****  *READ_GR_SNOW* - routine to read snow surface fields
!!
!!    PURPOSE
!!    -------
!       Initialize snow surface fields.
!
!!**  METHOD
!!    ------
!!    
!!    
!!
!!    EXTERNAL
!!    --------
!!      
!!       
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!    REFERENCE
!!    ---------
!!      
!!      
!!
!!    AUTHOR
!!    ------
!!      V. Masson       * Meteo France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       20/01/99
!       F.solmon       06/00 adaptation for patch
!       V.Masson       01/03 new version of ISBA
!       B. Decharme    2008  If no WSNOW, WSNOW = XUNDEF
!-----------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
USE modd_type_snow
!
USE modi_read_surf
!
USE modi_allocate_gr_snow
USE modi_pack_same_rank
!
USE modd_surf_par, ONLY : xundef
USE modd_prep_snow, ONLY : lsnow_frac_tot
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   declarations of arguments
!
 CHARACTER(LEN=6),   INTENT(IN)           :: HPROGRAM  ! calling program
 CHARACTER (LEN=*),  INTENT(IN)           :: HSURFTYPE ! generic name used for
                                                      ! snow characteristics
                                                      ! storage in file
 CHARACTER (LEN=3),  INTENT(IN)           :: HPREFIX   ! generic name for patch
!                                                     ! identification                      
INTEGER,            INTENT(IN)           :: KLU       ! horizontal size of snow var.
INTEGER,            INTENT(IN)           :: KSIZE_P
INTEGER, DIMENSION(:), INTENT(IN) :: KMASK_P
INTEGER,            INTENT(IN)           :: KPATCH    ! number of tiles
TYPE(surf_snow), INTENT(INOUT)           :: TPSNOW    ! snow characteristics
 CHARACTER (LEN=1),  INTENT(IN), OPTIONAL :: HDIR      ! type of reading
!                                                     ! HDIR = 'A' : entire field on All processors
!                                                     ! HDIR = 'H' : distribution on each processor
!
INTEGER,            INTENT(IN), OPTIONAL :: KVERSION
INTEGER,            INTENT(IN), OPTIONAL :: KBUGFIX
INTEGER, INTENT(IN), OPTIONAL :: KNPATCH
!
!*       0.2   declarations of local variables
!
 CHARACTER (LEN=7) :: YFMT0               ! format for writing
 CHARACTER (LEN=100) :: YFMT                ! format for writing
 CHARACTER(LEN=16)   :: YRECFM2 
 CHARACTER(LEN=12)   :: YRECFM              ! Name of the article to be read
 CHARACTER(LEN=4)    :: YNLAYER     !Format depending on the number of layers
 CHARACTER(LEN=1)    :: YDIR                ! type of reading
 CHARACTER(LEN=3) :: YPAT
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK
!
INTEGER             :: IRESP, JI, JP              ! Error code after redding
INTEGER             :: ISURFTYPE_LEN, IPAT_LEN     ! 
INTEGER             :: JL              ! loop counter
INTEGER             :: IVERSION, IBUGFIX
INTEGER :: INPATCH
!
LOGICAL :: GVERSION, GDIM, GDIM2
LOGICAL             :: GSNOW               ! snow written in the file
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('READ_GR_SNOW_1',0,zhook_handle)
!
ydir = 'H'
IF (PRESENT(hdir)) ydir = hdir
!
inpatch = 1
IF (PRESENT(knpatch)) inpatch = knpatch
!
!-------------------------------------------------------------------------------
IF(PRESENT(kversion))THEN
  iversion=kversion
ELSE
  CALL read_surf(hprogram,'VERSION',iversion,iresp)
ENDIF
IF(PRESENT(kbugfix))THEN
  ibugfix=kbugfix
ELSE
  CALL read_surf(hprogram,'BUG',ibugfix,iresp)
ENDIF
!
gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
gdim2 = gdim
IF (gdim) CALL read_surf(hprogram,'SPLIT_PATCH',gdim2,iresp)
!
!-------------------------------------------------------------------------------
!
gversion = (iversion>7 .OR. iversion==7 .AND. ibugfix>=3)
!
!*       1.    Type of snow scheme
!              -------------------
!
isurftype_len=len_trim(hsurftype)
!
IF (kpatch<=1) THEN

  IF (iversion <=2 .OR. (iversion==3 .AND. ibugfix<=4)) THEN
    WRITE(yfmt,'(A5,I1,A4)')     '(A5,A',isurftype_len,',A5)'
    WRITE(yrecfm2,yfmt) 'SNOW_',hsurftype,'_TYPE'
  ELSE
    IF (iversion<7 .OR. iversion==7 .AND. ibugfix<3) THEN
      WRITE(yfmt,'(A5,I1,A4)')     '(A3,A',isurftype_len,',A5)'
      WRITE(yrecfm2,yfmt) 'SN_',hsurftype,'_TYPE'
    ELSE
      WRITE(yfmt,'(A5,I1,A4)')     '(A3,A',isurftype_len,',A4)'
      WRITE(yrecfm2,yfmt) 'SN_',hsurftype,'_TYP'
      yrecfm2=adjustl(hprefix//yrecfm2)
    ENDIF
  END IF
  !
  CALL read_surf(hprogram,yrecfm2,tpsnow%SCHEME,iresp)
  !
  !*       2.    Snow levels
  !              -----------
  !
  !
  IF (iversion <=2 .OR. (iversion==3 .AND. ibugfix<=4)) THEN
    WRITE(yfmt,'(A5,I1,A4)')     '(A5,A',isurftype_len,',A6)'
    WRITE(yrecfm2,yfmt) 'SNOW_',hsurftype,'_LAYER'
  ELSE
    WRITE(yfmt,'(A5,I1,A4)')     '(A3,A',isurftype_len,',A2)'
    WRITE(yrecfm2,yfmt) 'SN_',hsurftype,'_N'
    IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm2=adjustl(hprefix//yrecfm2)
  END IF
  !
  CALL read_surf(hprogram,yrecfm2,tpsnow%NLAYER,iresp)
  !
ENDIF
!
!*       2.    Presence of snow fields in the file
!              -----------------------------------
!
IF (kpatch>0.AND.gdim2) THEN
  WRITE(ypat,'(I2)') kpatch
  ypat = "P"//adjustl(ypat)
  ipat_len = len_trim(adjustl(ypat))        
ELSE
  ypat = " "
  ipat_len=1
ENDIF
!
IF (iversion >6 .OR. (iversion==6 .AND. ibugfix>=1)) THEN
  WRITE(yfmt,'(A5,I1,A2,I1,A1)') '(A3,A',isurftype_len,',A1,A',ipat_len,')'
  WRITE(yrecfm,yfmt) 'SN_',adjustl(hsurftype(:len_trim(hsurftype))),adjustl(ypat(:len_trim(ypat)))
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm=adjustl(hprefix//yrecfm)
  CALL read_surf(hprogram,yrecfm,gsnow,iresp)
ELSE
  IF (tpsnow%NLAYER==0) THEN
    gsnow = .false.
    IF (tpsnow%SCHEME=='D95' .OR. tpsnow%SCHEME=='1-L' .OR. tpsnow%SCHEME=='EBA') tpsnow%NLAYER=1
    IF (tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO'                          ) tpsnow%NLAYER=3
  ELSE
    gsnow = .true.
  END IF
END IF
!
IF (.NOT.gdim2) THEN
  ypat = " "
  ipat_len=1
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       3.    Allocations
!              -----------
!
 CALL allocate_gr_snow(tpsnow,ksize_p)
!
IF (.NOT. gsnow) THEN
  IF (lhook) CALL dr_hook('READ_GR_SNOW_1',1,zhook_handle)
  RETURN
END IF
!-------------------------------------------------------------------------------
!
!*       4.    Additional key
!              ---------------
!
IF (iversion >= 7 .AND. hsurftype=='VEG'.AND.kpatch==1)  &
  CALL read_surf(hprogram,'LSNOW_FRAC_T',lsnow_frac_tot,iresp)
!
!-------------------------------------------------------------------------------
!
!
!*       5.    Snow reservoir
!              --------------
!
ALLOCATE(zwork(klu,inpatch))
!
IF (tpsnow%SCHEME=='1-L' .OR. tpsnow%SCHEME=='D95' .OR. tpsnow%SCHEME=='EBA' &
                         .OR. tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN 
  !
  WRITE(yfmt0,'(A5,I1,A1)') ',A1,A',isurftype_len
  !
  IF (gversion) THEN
    yfmt = '(A3'//yfmt0//')'
  ELSE
    yfmt = '(A5'//yfmt0//')'
  ENDIF
  CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"WSNOW",hsurftype,tpsnow%WSNOW)
  CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"RSNOW",hsurftype,tpsnow%RHO)
  !
  !*       7.    Snow temperature
  !              ----------------
  !
  IF (tpsnow%SCHEME=='1-L') THEN
    !
    CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"TSNOW",hsurftype,tpsnow%T)
    !
  ENDIF
  !
  !*       8.    Heat content
  !              ------------
  !
  IF (tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN
    !
    CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"HSNOW",hsurftype,tpsnow%HEAT)
    !
    IF (tpsnow%SCHEME=='CRO') THEN
      !
      CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"SHIST",hsurftype,tpsnow%HIST)
      !
      !*       9.    Snow Gran1
      !              ------------
      !
      IF (gversion) THEN
        yfmt = "(A2,A1"//yfmt0//')'       
      ELSE
        yfmt = "(A5"//yfmt0//')'
      ENDIF
      yfmt = yfmt//ynlayer//')'
      CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"SGRAN",hsurftype,tpsnow%GRAN1,hrec2="1")
      CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"SGRAN",hsurftype,tpsnow%GRAN2,hrec2="2")
      !
    ENDIF
    !
    IF ((tpsnow%SCHEME=='3-L'.AND.iversion>=8) .OR. tpsnow%SCHEME=='CRO') THEN
      !*       12.    Age parameter
      !              -------------------
      !
      IF (gversion) THEN
        yfmt = "(A3"//yfmt0//')'         
      ELSE
        yfmt = "(A4"//yfmt0//')'
      ENDIF
      yfmt = yfmt//ynlayer//')'
      CALL read_layers(gversion,tpsnow%NLAYER,ydir,hprefix,yfmt,"SAGE",hsurftype,tpsnow%AGE)
      !
    ELSE
      !
      DO jl = 1,tpsnow%NLAYER
        WHERE (tpsnow%WSNOW(:,1) >= 0.0) 
          tpsnow%AGE(:,jl) = 0.0
        ELSEWHERE
          tpsnow%AGE(:,jl) = xundef
        ENDWHERE
      ENDDO
      !
    END IF    
    !
  ENDIF
  !
  WRITE(yfmt,'(A5,I1,A2,I1,A1)') '(A4,A',isurftype_len,',A',ipat_len,')'
  WRITE(yrecfm,yfmt) 'ASN_',adjustl(hsurftype(:len_trim(hsurftype))),adjustl(ypat)
  IF (gversion) yrecfm=adjustl(hprefix//yrecfm)
  IF (gdim2) THEN
    CALL read_surf(hprogram,yrecfm,zwork(:,1),iresp,hdir=ydir)
    CALL pack_same_rank(kmask_p,zwork(:,1),tpsnow%ALB(:))
  ELSE
    CALL read_surf(hprogram,yrecfm,zwork,iresp,hdir=ydir)
    CALL pack_same_rank(kmask_p,zwork(:,max(1,kpatch)),tpsnow%ALB(:))
  ENDIF
  !
ENDIF
!
DEALLOCATE(zwork)
!
IF (lhook) CALL dr_hook('READ_GR_SNOW_1',1,zhook_handle)
!
IF (tpsnow%SCHEME=='1-L' .OR. tpsnow%SCHEME=='D95' .OR. tpsnow%SCHEME=='EBA' &
                         .OR. tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN 
  !
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('READ_GR_SNOW_2',0,zhook_handle_omp)
!$OMP DO PRIVATE(JI,JL)
  DO ji = 1,SIZE(tpsnow%WSNOW,1)
    !
    IF (tpsnow%WSNOW(ji,1) == 0.0 ) THEN
      !
      tpsnow%ALB(ji) = xundef
      !
      DO jl = 1,tpsnow%NLAYER
        !
        tpsnow%RHO(ji,jl)=xundef
        IF (tpsnow%SCHEME=='1-L') THEN
          tpsnow%T(ji,jl) = xundef
        ELSEIF (tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN
          tpsnow%HEAT(ji,jl) = xundef
          IF (tpsnow%SCHEME=='CRO') THEN
            tpsnow%HIST (ji,jl) = xundef
            tpsnow%GRAN1(ji,jl) = xundef
            tpsnow%GRAN2(ji,jl) = xundef
            tpsnow%AGE  (ji,jl) = xundef
          ENDIF
        ENDIF
        !
      ENDDO
    ENDIF
  ENDDO
!$OMP ENDDO
IF (lhook) CALL dr_hook('READ_GR_SNOW_2',1,zhook_handle_omp)
!$OMP END PARALLEL
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
CONTAINS
!
SUBROUTINE read_layers(OVERSION,KNL,HDIRIN,HPREF,HFMT,HREC,HSURF,PTAB,HREC2)
!
USE mode_read_surf_layers
!
IMPLICIT NONE
!
LOGICAL, INTENT(IN) :: OVERSION
INTEGER, INTENT(IN) :: KNL
 CHARACTER(LEN=*), INTENT(IN) :: HDIRIN
 CHARACTER(LEN=*), INTENT(IN) :: HPREF
 CHARACTER(LEN=*), INTENT(IN) :: HFMT
 CHARACTER(LEN=*), INTENT(IN) :: HREC
 CHARACTER(LEN=*), INTENT(IN) :: HSURF
REAL, DIMENSION(:,:), INTENT(OUT) :: PTAB
 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: HREC2
!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWORK3D
 CHARACTER(LEN=1) :: YREC2
 CHARACTER(LEN=12)   :: YRECFM   ! Name of the article to be read
INTEGER :: JL, IRESP
!
IF (PRESENT(hrec2)) THEN
  yrec2=trim(hrec2)
ELSE
  yrec2=""
ENDIF
!
IF (yrec2/="") THEN
  WRITE(yrecfm,hfmt) trim(hrec),trim(yrec2),'_',trim(hsurf)
ELSE
  WRITE(yrecfm,hfmt) trim(hrec),'_',trim(hsurf)
ENDIF
IF (oversion) yrecfm=adjustl(trim(hpref)//yrecfm)
!
IF (gdim2) THEN
  ALLOCATE(zwork3d(klu,SIZE(ptab,2),1))
ELSE
  ALLOCATE(zwork3d(klu,SIZE(ptab,2),inpatch))
ENDIF
!
 CALL read_surf_layers(hprogram,yrecfm,gdim2,zwork3d,iresp,kpatch=kpatch,hdir=ydir)
!
DO jl = 1,knl
  IF (gdim2) THEN
    CALL pack_same_rank(kmask_p,zwork3d(:,jl,1),ptab(:,jl))
  ELSE
    CALL pack_same_rank(kmask_p,zwork3d(:,jl,max(1,kpatch)),ptab(:,jl))
  ENDIF
ENDDO
!
DEALLOCATE(zwork3d)
!
END SUBROUTINE read_layers
!
END SUBROUTINE read_gr_snow