pgd_topo_index.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 pgd_topo_index (DTCO, UG, U, USS, S, OCTI, &
                                 HPROGRAM,KLU,HCTI,HCTIFILETYPE,OIMP_CTI)
!     ##################################################################
!
!!**** *PGD_TOPO_INDEX* monitor for computing topographic index statistics used by TOIPMODEL
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    B. Decharme        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    06/2009
!!    B. Decharme 05/2013 New topographic index linear regression for Topmodel
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_sso_n, ONLY : sso_t
USE modd_isba_n, ONLY : isba_s_t
!
USE modd_surfex_mpi, ONLY : nrank, npio
!
USE modd_pgd_grid,       ONLY : nl
!
!
USE modd_pgdwork,        ONLY : xall, xext_all, nsize_all, &
                                xmin_work, xmax_work,     &
                                xmean_work, xstd_work,    &
                                xskew_work, nsize, xsumval 
!
USE modd_surf_par,       ONLY : xundef
!
!
USE modd_sgh_par,        ONLY : xregp, xrega
!
USE modi_get_luout
USE modi_get_grid_coord
USE modi_read_surf
USE modi_treat_field
USE modi_pack_same_rank
USE modi_interpol_field
!
USE modi_init_io_surf_n
USE modi_end_io_surf_n
#ifdef SFX_ASC
USE modd_io_surf_asc, ONLY : cfilein
#endif
#ifdef SFX_FA
USE modd_io_surf_fa,  ONLY : cfilein_fa
#endif
#ifdef SFX_LFI
USE modd_io_surf_lfi, ONLY : cfilein_lfi
#endif
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
USE modi_get_surf_mask_n
!
USE modi_get_type_dim_n
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
TYPE(isba_s_t), INTENT(INOUT) :: S
!
LOGICAL, INTENT(INOUT) :: OCTI
!
 CHARACTER(LEN=6),     INTENT(IN)  :: HPROGRAM     ! program calling
INTEGER,              INTENT(IN)  :: KLU          ! number of nature points
 CHARACTER(LEN=28),    INTENT(IN)  :: HCTI         ! topographic index file name
 CHARACTER(LEN=6),     INTENT(IN)  :: HCTIFILETYPE ! topographic index file type
LOGICAL,              INTENT(IN)  :: OIMP_CTI     ! .true. if topographic index statistics is imposed
!
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
REAL, DIMENSION(:), ALLOCATABLE :: ZLAT, ZDELTA, ZMEAN_INI, &
                                   ZTI_MEAN, ZTI_STD, ZTI_SKEW
!
INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK
!
LOGICAL :: LREG,LREG10,LREG2
!
INTEGER :: IFULL     ! total number of points
INTEGER :: I_DIM
INTEGER :: IRET      ! error code
INTEGER :: ILUOUT    ! output listing logical unit
!
 CHARACTER(LEN=6  ) :: YFILETYPE, YSCHEME, YSUBROUTINE
 CHARACTER(LEN=20)  :: YFIELD        ! Name of the field.
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!*    1.      Initializations
!             ---------------
!
IF (lhook) CALL dr_hook('PGD_TOPO_INDEX',0,zhook_handle)
IF(len_trim(hcti)==0)THEN
!
  ALLOCATE(s%XTI_MIN (0))
  ALLOCATE(s%XTI_MAX (0))
  ALLOCATE(s%XTI_MEAN(0))
  ALLOCATE(s%XTI_STD (0))
  ALLOCATE(s%XTI_SKEW(0))
!        
!-------------------------------------------------------------------------------
ELSE
!-------------------------------------------------------------------------------
!
  octi = .true.
!
!*    2.      Find LUOUT
!             ----------
!
  CALL get_luout(hprogram,iluout)
!
  WRITE(iluout,*) '*****************************************'
  WRITE(iluout,*) 'Comput Topographic indexes for TOPMODEL  '
  WRITE(iluout,*) '*****************************************'
!
!*    3.      Allocations of statistics arrays
!             --------------------------------
!
  ALLOCATE(s%XTI_MIN (klu))
  ALLOCATE(s%XTI_MAX (klu))
  ALLOCATE(s%XTI_MEAN(klu))
  ALLOCATE(s%XTI_STD (klu))
  ALLOCATE(s%XTI_SKEW(klu))
!
  s%XTI_MIN (:) = xundef
  s%XTI_MAX (:) = xundef
  s%XTI_MEAN(:) = xundef
  s%XTI_STD (:) = xundef
  s%XTI_SKEW(:) = xundef
!
!-------------------------------------------------------------------------------
!
!*    4.      Allocations of work arrays
!             --------------------------
!
  CALL get_type_dim_n(dtco, u, 'NATURE',i_dim)
  IF (i_dim/=klu) THEN
     WRITE(iluout,*)'PGD_TOPO_INDEX: Wrong dimension of MASK: ',i_dim,klu
     CALL abor1_sfx('PGD_TOPO_INDEX: WRONG DIMENSION OF MASK')
  ENDIF
!
  ALLOCATE(imask(klu))
  ifull=0
  CALL get_surf_mask_n(dtco, u, &
                       'NATURE',klu,imask,ifull,iluout)
  IF (ifull/=nl) THEN
     WRITE(iluout,*)'PGD_TOPO_INDEX: Wrong dimension of FULL: ',ifull,nl
     CALL abor1_sfx('PGD_TOPO_INDEX: WRONG DIMENSION OF FULL')
  ENDIF
!
  ALLOCATE(xmin_work(ifull))
  ALLOCATE(xmax_work(ifull))
  ALLOCATE(xmean_work(ifull))
  ALLOCATE(xstd_work(ifull))
  ALLOCATE(xskew_work(ifull))
!
  xmin_work(:)=xundef
  xmax_work(:)=xundef
  xmean_work(:)=xundef
  xstd_work(:)=xundef
  xskew_work(:)=xundef
!
!-------------------------------------------------------------------------------
!
!*    5.      Use of imposed field
!             --------------------
!
  IF (oimp_cti) THEN
!
     yfiletype=hctifiletype
     IF(hctifiletype=='NETCDF')THEN
        CALL abor1_sfx('Use another format than netcdf for cti input file with LIMP_CTI')
     ELSE
#ifdef SFX_ASC
       cfilein     = adjustl(adjustr(hcti)//'.txt')
#endif
#ifdef SFX_FA
       cfilein_fa  = adjustl(adjustr(hcti)//'.fa')
#endif
#ifdef SFX_LFI
       cfilein_lfi = adjustl(hcti)
#endif
 CALL init_io_surf_n(dtco, u, &
                           yfiletype,'FULL  ','SURF  ','READ ')
     ENDIF     
!   
     CALL read_surf(yfiletype,'TI_MIN' ,xmin_work ,iret) 
     CALL read_surf(yfiletype,'TI_MAX' ,xmax_work ,iret)
     CALL read_surf(yfiletype,'TI_MEAN',xmean_work,iret)
     CALL read_surf(yfiletype,'TI_STD' ,xstd_work ,iret) 
     CALL read_surf(yfiletype,'TI_SKEW',xskew_work,iret) 
!
     CALL end_io_surf_n(yfiletype)
!
  ELSE
!
!-------------------------------------------------------------------------------
!
!*    6.      Use of cti file
!             ---------------
!
     ALLOCATE(nsize_all(u%NDIM_FULL,1))
     ALLOCATE(xext_all(u%NDIM_FULL,2))
     ALLOCATE(xall(u%NDIM_FULL,3,1))     
!
     nsize_all(:,1) = 0.
     xext_all (:,1) = -99999.
     xext_all (:,2) = 99999.
     xall   (:,:,1) = 0.     
!
     xmax_work(:) =-99999.
!
     yfield      = 'CTI'
     yscheme     = 'SURF  '
     ysubroutine = 'A_CTI '
     CALL treat_field(ug, u, uss, &
                      hprogram,yscheme,hctifiletype,ysubroutine,hcti,yfield)
!
!-------------------------------------------------------------------------------
!
!*    7.      Coherence
!             ---------
!
     WHERE(nsize(:,1)<36.OR.xstd_work(:)==0.0)
          xmin_work(:) = xundef
          xmax_work(:) = xundef
          xmean_work(:) = xundef
          xstd_work(:) = xundef
          xskew_work(:) = xundef
          nsize(:,1) = 0
     ENDWHERE 
!
     WHERE(u%XNATURE(:)>0.0.AND.xskew_work(:)<=-8.0)
          xmin_work(:) = xundef
          xmax_work(:) = xundef
          xmean_work(:) = xundef
          xstd_work(:) = xundef
          xskew_work(:) = xundef
          nsize(:,1) = 0
     ENDWHERE             
!
     WHERE(u%XNATURE(:)==0.)
          xmin_work(:) = xundef
          xmax_work(:) = xundef
          xmean_work(:) = xundef
          xstd_work(:) = xundef
          xskew_work(:) = xundef
          nsize(:,1) = 0
     ENDWHERE   
!
!-------------------------------------------------------------------------------
!
!*    8.      Regression 1km to 100m or 2m
!             ----------------------------
!
     IF(hctifiletype=='DIRECT')THEN
!
!      Topographic index linear regression for Topmodel if topo ref at 1km
!      pan and king (2012) 1km to 2m
!    
       CALL ctireg(lreg,lreg10,lreg2)
!
       IF(all(xmean_work(:)==xundef))lreg=.false.
!
       IF(lreg)THEN
!          
         WRITE(iluout,*)'WITH DIF, TOPO INDEX USED REGRESSIONS OF ' 
!
         ALLOCATE(zdelta(ifull))
         ALLOCATE(zmean_ini(ifull))
         ALLOCATE(zti_mean(ifull))
         ALLOCATE(zti_std(ifull))
         ALLOCATE(zti_skew(ifull))
!
         zmean_ini(:)=xmean_work(:)
         zti_mean(:)=xundef
         zti_std(:)=xundef
         zti_skew(:)=xundef
         zdelta(:)= 0.0
!
!        1km to 10m
         IF(lreg10)THEN
           WRITE(iluout,*)' PAN AND KING (2012) 1km to 10m '
           WHERE(xmean_work(:)/=xundef.AND.(xmax_work(:)-xmin_work(:))>0.2)
               zti_mean(:)= 1.136+0.657*xmean_work(:)-0.640*xstd_work(:)+0.053*xskew_work(:)
               zti_std(:)=-0.128+0.187*xmean_work(:)+0.168*xstd_work(:)-0.261*xskew_work(:)
               zti_skew(:)= 3.768-0.246*xmean_work(:)+0.317*xstd_work(:)+0.222*xskew_work(:)
           ENDWHERE
             xmean_work(:)=zti_mean(:)
             xstd_work(:)=zti_std(:)
             xskew_work(:)=zti_skew(:)
         ELSE
             zti_mean(:)=xmean_work(:)
             zti_std(:)=xstd_work(:)
             zti_skew(:)=xskew_work(:)
         ENDIF
!
!          10m to 2m
         IF(lreg2)THEN
           WRITE(iluout,*)' PAN AND KING (2012) 10m to 2m '
           WHERE(xmean_work(:)/=xundef.AND.(xmax_work(:)-xmin_work(:))>0.2)
                   xmean_work(:)=-3.826+1.402*zti_mean(:)-0.434*zti_std(:)+0.328*zti_skew(:)
                   xstd_work(:)= 3.655-0.209*zti_mean(:)+0.440*zti_std(:)-0.091*zti_skew(:)
                   xskew_work(:)= 2.266-0.023*zti_mean(:)-0.245*zti_std(:)-0.240*zti_skew(:)
           ENDWHERE
         ENDIF
!
         WHERE(xmean_work(:)/=xundef.AND.(xmax_work(:)-xmin_work(:))>0.2)
               xstd_work(:)=max(0.2,xstd_work(:))
               xskew_work(:)=max(0.2,xskew_work(:))
         ENDWHERE         
!           
         WHERE(xmean_work(:)>0.0.AND.xmean_work(:)/=xundef)
               zdelta(:)= (xmean_work(:)-zmean_ini(:))
               xmin_work(:)= max(             0.0,xmin_work(:)+zdelta(:))
               xmax_work(:)= max(xmin_work(:)+0.2,xmax_work(:)+zdelta(:))
         ELSEWHERE
              xmin_work(:) = xundef
              xmax_work(:) = xundef
              xmean_work(:) = xundef
              xstd_work(:) = xundef
              xskew_work(:) = xundef
              nsize(:,1) = 0
         ENDWHERE
!
         DEALLOCATE(zdelta   )
         DEALLOCATE(zmean_ini)
         DEALLOCATE(zti_mean )
         DEALLOCATE(zti_std  )
         DEALLOCATE(zti_skew )
!     
       ENDIF
!  
     ENDIF
!
!-------------------------------------------------------------------------------
!
!*    10.     Interpolation if some points are not initialized (no data for these points)
!             ------------------------------------------------
!
    WRITE(iluout,*) '*********************************************'
    WRITE(iluout,*) 'Interpolation if some index not initialized  '
    WRITE(iluout,*) '*********************************************'
!
    ALLOCATE(zlat(nl))
    CALL get_grid_coord(ug%G%CGRID, ug%G%NGRID_PAR, ug%G%XGRID_PAR, u%NSIZE_FULL, &
                        iluout,py=zlat)
!
    WHERE (u%XNATURE(:)==0..AND.nsize(:,1)==0) nsize(:,1) = -1
!
!   No Antarctic
    WHERE(u%XNATURE(:)>0..AND.zlat(:)<-60.)
          xmin_work(:) = xundef
          xmax_work(:) = xundef
          xmean_work(:) = xundef
          xstd_work(:) = xundef
          xskew_work(:) = xundef
          nsize(:,1) = -1
    ENDWHERE   
!
    IF(all(nsize(:,1)==0.0))nsize(:,1)=-1
!
    CALL interpol_field(ug, u, hprogram,iluout,nsize(:,1),xmin_work(:),'TI_MIN ',pdef=xundef,knpts=1)
    CALL interpol_field(ug, u, hprogram,iluout,nsize(:,1),xmax_work(:),'TI_MAX ',pdef=xundef,knpts=1)
    CALL interpol_field(ug, u, hprogram,iluout,nsize(:,1),xmean_work(:),'TI_MEAN',pdef=xundef,knpts=1)
    CALL interpol_field(ug, u, hprogram,iluout,nsize(:,1),xstd_work(:),'TI_STD ',pdef=xundef,knpts=1)
    CALL interpol_field(ug, u, hprogram,iluout,nsize(:,1),xskew_work(:),'TI_SKEW',pdef=xundef,knpts=1)
!
    DEALLOCATE(nsize     )
    DEALLOCATE(xsumval   )
    DEALLOCATE(zlat      )
!
  ENDIF
!-------------------------------------------------------------------------------
!
!*    11.     Asign parameters
!             ----------------
!
  CALL pack_same_rank(imask,xmin_work ,s%XTI_MIN)
  CALL pack_same_rank(imask,xmax_work ,s%XTI_MAX)
  CALL pack_same_rank(imask,xmean_work,s%XTI_MEAN)
  CALL pack_same_rank(imask,xstd_work ,s%XTI_STD)
  CALL pack_same_rank(imask,xskew_work,s%XTI_SKEW)
!  
!-------------------------------------------------------------------------------
!
  DEALLOCATE(xmin_work )
  DEALLOCATE(xmax_work )
  DEALLOCATE(xmean_work)
  DEALLOCATE(xstd_work )
  DEALLOCATE(xskew_work)
!
  WRITE(iluout,*) '******************************'
  WRITE(iluout,*) 'End Comput Topographic indexes'
  WRITE(iluout,*) '******************************'
!
ENDIF
!
IF (lhook) CALL dr_hook('PGD_TOPO_INDEX',1,zhook_handle)
!-------------------------------------------------------------------------------
!
CONTAINS
!
SUBROUTINE ctireg(OREG,OREG10,OREG2)  
!      
!*    0.     DECLARATION
!            -----------
!
USE modi_open_namelist
USE modi_close_namelist
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
LOGICAL,           INTENT(OUT) :: OREG,OREG10,OREG2        ! regression key
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
REAL    :: ZDLAT                      ! latitude mesh in the data file
REAL    :: ZDLON                      ! longitude mesh in the data file
REAL    :: ZGLBLATMIN                 ! minimum latitude of data box in the file
REAL    :: ZGLBLONMIN                 ! minimum longitude of data box in the file
REAL    :: ZGLBLATMAX                 ! maximum latitude of data box in the file
REAL    :: ZGLBLONMAX                 ! maximum longitude of data box in the file
!
 CHARACTER(LEN=28)  :: YFILEHDR        ! Name of the field file header
!
INTEGER                    :: INBLAT
INTEGER                    :: INBLON
INTEGER                    :: IGLB
INTEGER                    :: JHEAD      ! loop control
INTEGER                    :: ININDEX    ! index of character 'N' in YSTRING1
INTEGER                    :: ISINDEX    ! index of character 'S' in YSTRING1
INTEGER                    :: IEINDEX    ! index of character 'E' in YSTRING1
INTEGER                    :: IWINDEX    ! index of character 'W' in YSTRING1
REAL, DIMENSION(7)         :: ZVAL       ! values of the head data
INTEGER                    :: IHEAD      ! index of the data in the array ZVAL
 CHARACTER(LEN=100)         :: YSTRING    ! total string in the head
 CHARACTER(LEN=100)         :: YSTRING1   ! string less the begining line descriptor
 CHARACTER(LEN=100)         :: YVAL       ! absolute value of the data of the line
INTEGER                    :: IPOINT     ! index of '.' in the string YVAL
INTEGER                    :: ILENGTH    ! length of the string YVAL
INTEGER                    :: IFRACLENGTH! length of the fractional part in string YVAL
 CHARACTER(LEN=2)           :: YLENGTH    ! length of the string YVAL
 CHARACTER(LEN=2)           :: YFRACLENGTH! length of the fractional part in string YVAL
 CHARACTER(LEN=10)          :: YINTERNALFORMAT ! format to read YVAL in real ZVAL
!
REAL    :: Z1000M, Z100M, Z10M
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('PGD_TOPO_INDEX:CTIREG',0,zhook_handle)
!
oreg  =.false.
oreg10=.false.
oreg2 =.false.
!
iglb=11
yfilehdr =adjustl(adjustr(hcti)//'.hdr')
 CALL open_namelist(hprogram,iglb,yfilehdr)
!
!*         1.    Line of comments
!                ----------------
!
READ (iglb,'(A100)',end=99) ystring
!
!-------------------------------------------------------------------------------
!
!*         2.    Other lines
!                -----------
!
DO jhead=1,7
  READ (iglb,'(A100)',end=99) ystring
  ystring=adjustl(ystring)
!
!*         2.1   Selection of the line
!                ---------------------
!
  SELECT CASE (ystring(1:5))
         CASE('cutva')
           ihead=1
           ystring1=ystring(10:100)
         CASE('nodat')
           ihead=1
           ystring1=ystring(8:100)
         CASE('north')
           ihead=2
           ystring1=ystring(7:100)           
         CASE('south')
           ihead=3
           ystring1=ystring(7:100) 
         CASE('east:')
           ihead=4
           ystring1=ystring(6:100)  
         CASE('west:')
           ihead=5
           ystring1=ystring(6:100)  
         CASE('rows:')
           ihead=6
           ystring1=ystring(6:100) 
         CASE('cols:')
           ihead=7
           ystring1=ystring(6:100) 
  END SELECT
!
!*         2.2   Test on presence of geographical descritor (N, E, S or W)
!                ---------------------------------------------------------
!
  inindex=index(ystring1,'N')
  isindex=index(ystring1,'S')
  ieindex=index(ystring1,'E')
  iwindex=index(ystring1,'W')
  yval=adjustl(ystring1)
  IF (inindex/=0) yval=adjustl(ystring1(1:inindex-1))
  IF (isindex/=0) yval='-'//adjustl(ystring1(1:isindex-1))
  IF (ieindex/=0) yval=adjustl(ystring1(1:ieindex-1))
  IF (iwindex/=0) yval='-'//adjustl(ystring1(1:iwindex-1))
!
!*         2.3   Transformation of the data in real
!                ----------------------------------
!
  ipoint=index(yval,'.')
  IF (ipoint==0) yval=adjustl(adjustr(yval)//'.')
!
!*         2.4   Definition of the format of the data
!                ------------------------------------
!
  ilength=len_trim(adjustl(adjustr(yval)))
  ifraclength=ilength-index(yval,'.')
  WRITE(ylength,'(I2)') ilength
  WRITE(yfraclength,'(I2)') ifraclength
  yinternalformat='(F'//ylength//'.'//yfraclength//')'
!
!*         2.5   Data definition
!                ---------------
!
  READ(yval,adjustl(yinternalformat)) zval(ihead)
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*         3.    Initialization of arguments, longitudes and latitudes
!                -----------------------------------------------------
!
zglblatmax=zval(2)
zglblatmin=zval(3)
zglblonmin=zval(5)
zglblonmax=zval(4)+nint((zval(5)-zval(4)+180.+1.e-10)/360.)*360.
inblat=nint(zval(6))
inblon=nint(zval(7))
!
zdlat=(zglblatmax-zglblatmin)/inblat
zdlon=(zglblonmax-zglblonmin)/inblon
!
z1000m = 30./3600. !(minute arc to m)
z100m  = z1000m/10.
z10m   = z1000m/100.
!
IF(sqrt(zdlat*zdlon)>z100m)THEN
   oreg  =.true.
   oreg10=.true.
   oreg2 =.true.
ENDIF
!
IF(sqrt(zdlat*zdlon)>=z10m.AND.sqrt(zdlat*zdlon)<=z100m)THEN
   oreg  =.true.
   oreg10=.false.
   oreg2 =.true.
ENDIF
!
 CALL close_namelist(hprogram,iglb)
!
IF (lhook) CALL dr_hook('PGD_TOPO_INDEX:CTIREG',1,zhook_handle)
RETURN
99 CONTINUE
 CALL abor1_sfx('CTIREG: PB READING TOPO INDEX FILE HEADER')
IF (lhook) CALL dr_hook('PGD_TOPO_INDEX:CTIREG',1,zhook_handle)
!
END SUBROUTINE ctireg
!
!-------------------------------------------------------------------------------
END SUBROUTINE pgd_topo_index