SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/SURFEX/pgd_topo

Go to the documentation of this file.
00001 !     #########
00002       SUBROUTINE PGD_TOPO_INDEX(HPROGRAM,KLU,HCTI,HCTIFILETYPE,OIMP_CTI)
00003 !     ##################################################################
00004 !
00005 !!**** *PGD_TOPO_INDEX* monitor for computing topographic index statistics used by TOIPMODEL
00006 !!
00007 !!    PURPOSE
00008 !!    -------
00009 !!
00010 !!    METHOD
00011 !!    ------
00012 !!   
00013 !
00014 !!    EXTERNAL
00015 !!    --------
00016 !!
00017 !!    IMPLICIT ARGUMENTS
00018 !!    ------------------
00019 !!
00020 !!    REFERENCE
00021 !!    ---------
00022 !!
00023 !!    AUTHOR
00024 !!    ------
00025 !!
00026 !!    B. Decharme        Meteo-France
00027 !!
00028 !!    MODIFICATION
00029 !!    ------------
00030 !!
00031 !!    Original    06/2009
00032 !!
00033 !----------------------------------------------------------------------------
00034 !
00035 !*    0.     DECLARATION
00036 !            -----------
00037 !
00038 USE MODD_PGD_GRID,       ONLY : NL
00039 !
00040 USE MODD_SURF_ATM_n,     ONLY : XNATURE
00041 !
00042 USE MODD_PGDWORK,        ONLY : XSUMVAL, XSUMVAL2, NSIZE, &
00043                                   XMIN_WORK, XMAX_WORK,     &
00044                                   XMEAN_WORK, XSTD_WORK,    &
00045                                   XSKEW_WORK, XSUMVAL3  
00046 !
00047 USE MODD_SURF_PAR,       ONLY : XUNDEF
00048 !
00049 USE MODD_ISBA_n,         ONLY : LCTI,XTI_MIN,XTI_MAX,XTI_MEAN,XTI_STD,XTI_SKEW
00050 !
00051 USE MODI_GET_LUOUT
00052 USE MODI_GET_GRID_COORD
00053 USE MODI_READ_SURF
00054 USE MODI_TREAT_FIELD
00055 USE MODI_PACK_SAME_RANK
00056 USE MODI_INTERPOL_FIELD
00057 !
00058 USE MODI_INIT_IO_SURF_n
00059 USE MODI_END_IO_SURF_n
00060 #ifdef ASC
00061 USE MODD_IO_SURF_ASC, ONLY : CFILEIN
00062 #endif
00063 #ifdef FA
00064 USE MODD_IO_SURF_FA,  ONLY : CFILEIN_FA
00065 #endif
00066 #ifdef LFI
00067 USE MODD_IO_SURF_LFI, ONLY : CFILEIN_LFI
00068 #endif
00069 !
00070 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00071 USE PARKIND1  ,ONLY : JPRB
00072 !
00073 USE MODI_ABOR1_SFX
00074 !
00075 USE MODI_GET_SURF_MASK_n
00076 !
00077 USE MODI_GET_TYPE_DIM_n
00078 !
00079 IMPLICIT NONE
00080 !
00081 !*    0.1    Declaration of arguments
00082 !            ------------------------
00083 !
00084  CHARACTER(LEN=6),     INTENT(IN)  :: HPROGRAM     ! program calling
00085 INTEGER,              INTENT(IN)  :: KLU          ! number of nature points
00086  CHARACTER(LEN=28),    INTENT(IN)  :: HCTI         ! topographic index file name
00087  CHARACTER(LEN=6),     INTENT(IN)  :: HCTIFILETYPE ! topographic index file type
00088 LOGICAL,              INTENT(IN)  :: OIMP_CTI     ! .true. if topographic index statistics is imposed
00089 !
00090 !
00091 !*    0.2    Declaration of local variables
00092 !            ------------------------------
00093 !
00094 REAL, DIMENSION(:), ALLOCATABLE :: ZLAT
00095 !
00096 INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK
00097 !
00098 INTEGER :: IFULL     ! total number of points
00099 INTEGER :: I_DIM
00100 INTEGER :: IRET      ! error code
00101 INTEGER :: ILUOUT    ! output listing logical unit
00102 !
00103  CHARACTER(LEN=6  ) :: YFILETYPE, YSCHEME, YSUBROUTINE
00104  CHARACTER(LEN=20)  :: YFIELD        ! Name of the field.
00105 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00106 !
00107 !-------------------------------------------------------------------------------
00108 !
00109 !*    1.      Initializations
00110 !             ---------------
00111 !
00112 IF (LHOOK) CALL DR_HOOK('PGD_TOPO_INDEX',0,ZHOOK_HANDLE)
00113 IF(LEN_TRIM(HCTI)==0)THEN
00114 !
00115   ALLOCATE(XTI_MIN (0))
00116   ALLOCATE(XTI_MAX (0))
00117   ALLOCATE(XTI_MEAN(0))
00118   ALLOCATE(XTI_STD (0))
00119   ALLOCATE(XTI_SKEW(0))
00120 !        
00121 !-------------------------------------------------------------------------------
00122 ELSE
00123 !-------------------------------------------------------------------------------
00124 !
00125   LCTI = .TRUE.
00126 !
00127 !*    2.      Find LUOUT
00128 !             ----------
00129 !
00130   CALL GET_LUOUT(HPROGRAM,ILUOUT)
00131 !
00132   WRITE(ILUOUT,*) '*****************************************'
00133   WRITE(ILUOUT,*) 'Comput Topographic indexes for TOPMODEL  '
00134   WRITE(ILUOUT,*) '*****************************************'
00135 !
00136 !*    3.      Allocations of statistics arrays
00137 !             --------------------------------
00138 !
00139   ALLOCATE(XTI_MIN (KLU))
00140   ALLOCATE(XTI_MAX (KLU))
00141   ALLOCATE(XTI_MEAN(KLU))
00142   ALLOCATE(XTI_STD (KLU))
00143   ALLOCATE(XTI_SKEW(KLU))
00144 !
00145   XTI_MIN (:) = XUNDEF
00146   XTI_MAX (:) = XUNDEF
00147   XTI_MEAN (:) = XUNDEF
00148   XTI_STD (:) = XUNDEF
00149   XTI_SKEW(:) = XUNDEF
00150 !
00151 !-------------------------------------------------------------------------------
00152 !
00153 !*    4.      Allocations of work arrays
00154 !             --------------------------
00155 !
00156   CALL GET_TYPE_DIM_n('NATURE',I_DIM)
00157   IF (I_DIM/=KLU) THEN
00158      WRITE(ILUOUT,*)'PGD_TOPO_INDEX: Wrong dimension of MASK: ',I_DIM,KLU
00159      CALL ABOR1_SFX('PGD_TOPO_INDEX: WRONG DIMENSION OF MASK')
00160   ENDIF
00161 !
00162   ALLOCATE(IMASK(KLU))
00163   IFULL=0
00164   CALL GET_SURF_MASK_n('NATURE',KLU,IMASK,IFULL,ILUOUT)
00165   IF (IFULL/=NL) THEN
00166      WRITE(ILUOUT,*)'PGD_TOPO_INDEX: Wrong dimension of FULL: ',IFULL,NL
00167      CALL ABOR1_SFX('PGD_TOPO_INDEX: WRONG DIMENSION OF FULL')
00168   ENDIF
00169 !
00170   ALLOCATE(XMIN_WORK  (IFULL))
00171   ALLOCATE(XMAX_WORK  (IFULL))
00172   ALLOCATE(XMEAN_WORK (IFULL))
00173   ALLOCATE(XSTD_WORK  (IFULL))
00174   ALLOCATE(XSKEW_WORK (IFULL))
00175 !
00176   XMIN_WORK (:)=XUNDEF
00177   XMAX_WORK (:)=XUNDEF
00178   XMEAN_WORK(:)=XUNDEF
00179   XSTD_WORK (:)=XUNDEF
00180   XSKEW_WORK(:)=XUNDEF
00181 !
00182 !-------------------------------------------------------------------------------
00183 !
00184 !*    5.      Use of imposed field
00185 !             --------------------
00186 !
00187   IF (OIMP_CTI) THEN
00188 !
00189      YFILETYPE=HCTIFILETYPE
00190      IF(HCTIFILETYPE=='NETCDF')THEN
00191         CALL ABOR1_SFX('Use another format than netcdf for cti input file with LIMP_CTI')
00192      ELSE
00193 #ifdef ASC
00194        CFILEIN     = ADJUSTL(ADJUSTR(HCTI)//'.txt')
00195 #endif
00196 #ifdef FA
00197        CFILEIN_FA  = ADJUSTL(ADJUSTR(HCTI)//'.fa')
00198 #endif
00199 #ifdef LFI
00200        CFILEIN_LFI = ADJUSTL(HCTI)
00201 #endif
00202        CALL INIT_IO_SURF_n(YFILETYPE,'FULL  ','SURF  ','READ ')
00203      ENDIF     
00204 !   
00205      CALL READ_SURF(YFILETYPE,'TI_MIN' ,XMIN_WORK ,IRET) 
00206      CALL READ_SURF(YFILETYPE,'TI_MAX' ,XMAX_WORK ,IRET)
00207      CALL READ_SURF(YFILETYPE,'TI_MEAN',XMEAN_WORK,IRET)
00208      CALL READ_SURF(YFILETYPE,'TI_STD' ,XSTD_WORK ,IRET) 
00209      CALL READ_SURF(YFILETYPE,'TI_SKEW',XSKEW_WORK,IRET) 
00210 !
00211      CALL END_IO_SURF_n(YFILETYPE)
00212 !
00213   ELSE
00214 !
00215 !-------------------------------------------------------------------------------
00216 !
00217 !*    6.      Use of cti file
00218 !             ---------------
00219 !
00220      ALLOCATE(NSIZE   (IFULL))
00221      ALLOCATE(XSUMVAL (IFULL))
00222      ALLOCATE(XSUMVAL2(IFULL))
00223      ALLOCATE(XSUMVAL3(IFULL))
00224 !
00225      NSIZE    (:) = 0.
00226      XSUMVAL  (:) = 0.
00227      XSUMVAL2 (:) = 0.
00228      XSUMVAL3 (:) = 0.
00229 !
00230      XMAX_WORK(:) =-99999.
00231 !
00232      YFIELD      = 'CTI'
00233      YSCHEME     = 'SURF  '
00234      YSUBROUTINE = 'A_CTI '
00235      CALL TREAT_FIELD(HPROGRAM,YSCHEME,HCTIFILETYPE,YSUBROUTINE,HCTI,YFIELD)
00236 !
00237 !-------------------------------------------------------------------------------
00238 !
00239 !*    7.      Coherence
00240 !             ---------
00241 !
00242      WHERE(NSIZE(:)<36.OR.XSTD_WORK(:)==0.0)
00243           XMIN_WORK (:) = XUNDEF
00244           XMAX_WORK (:) = XUNDEF
00245           XMEAN_WORK(:) = XUNDEF
00246           XSTD_WORK (:) = XUNDEF
00247           XSKEW_WORK(:) = XUNDEF
00248           NSIZE     (:) = 0
00249      ENDWHERE 
00250 !
00251      WHERE(XNATURE(:)>0.0.AND.XSKEW_WORK(:)<=-8.0)
00252           XMIN_WORK (:) = XUNDEF
00253           XMAX_WORK (:) = XUNDEF
00254           XMEAN_WORK(:) = XUNDEF
00255           XSTD_WORK (:) = XUNDEF
00256           XSKEW_WORK(:) = XUNDEF
00257           NSIZE     (:) = 0
00258      ENDWHERE             
00259 !
00260      WHERE(XNATURE(:)==0.)
00261           XMIN_WORK (:) = XUNDEF
00262           XMAX_WORK (:) = XUNDEF
00263           XMEAN_WORK(:) = XUNDEF
00264           XSTD_WORK (:) = XUNDEF
00265           XSKEW_WORK(:) = XUNDEF
00266           NSIZE     (:) = 0
00267      ENDWHERE   
00268 !
00269 !-------------------------------------------------------------------------------
00270 !
00271 !*    5.      Interpolation if some points are not initialized (no data for these points)
00272 !             ------------------------------------------------
00273 !
00274     WRITE(ILUOUT,*) '*********************************************'
00275     WRITE(ILUOUT,*) 'Interpolation if some index not initialized  '
00276     WRITE(ILUOUT,*) '*********************************************'
00277 !
00278     ALLOCATE(ZLAT(NL))
00279     CALL GET_GRID_COORD(ILUOUT,PY=ZLAT)
00280 !
00281     WHERE (XNATURE(:)==0..AND.NSIZE(:)==0) NSIZE(:) = -1
00282 !
00283 !   No Antarctic
00284     WHERE(XNATURE(:)>0..AND.ZLAT(:)<-60.)
00285           XMIN_WORK (:) = XUNDEF
00286           XMAX_WORK (:) = XUNDEF
00287           XMEAN_WORK(:) = XUNDEF
00288           XSTD_WORK (:) = XUNDEF
00289           XSKEW_WORK(:) = XUNDEF
00290           NSIZE     (:) = -1
00291     ENDWHERE   
00292 !
00293     CALL INTERPOL_FIELD(HPROGRAM,ILUOUT,NSIZE,XMIN_WORK (:),'TI_MIN ',PDEF=XUNDEF,KNPTS=1)
00294     CALL INTERPOL_FIELD(HPROGRAM,ILUOUT,NSIZE,XMAX_WORK (:),'TI_MAX ',PDEF=XUNDEF,KNPTS=1)
00295     CALL INTERPOL_FIELD(HPROGRAM,ILUOUT,NSIZE,XMEAN_WORK(:),'TI_MEAN',PDEF=XUNDEF,KNPTS=1)
00296     CALL INTERPOL_FIELD(HPROGRAM,ILUOUT,NSIZE,XSTD_WORK (:),'TI_STD ',PDEF=XUNDEF,KNPTS=1)
00297     CALL INTERPOL_FIELD(HPROGRAM,ILUOUT,NSIZE,XSKEW_WORK(:),'TI_SKEW ',PDEF=XUNDEF,KNPTS=1)
00298 !
00299     DEALLOCATE(NSIZE     )
00300     DEALLOCATE(XSUMVAL   )
00301     DEALLOCATE(XSUMVAL2  )
00302     DEALLOCATE(XSUMVAL3  )
00303     DEALLOCATE(ZLAT      )
00304 !
00305   ENDIF
00306 !-------------------------------------------------------------------------------
00307 !
00308 !*    8.      Asign parameters
00309 !             ----------------
00310 !
00311   CALL PACK_SAME_RANK(IMASK,XMIN_WORK ,XTI_MIN)
00312   CALL PACK_SAME_RANK(IMASK,XMAX_WORK ,XTI_MAX)
00313   CALL PACK_SAME_RANK(IMASK,XMEAN_WORK,XTI_MEAN)
00314   CALL PACK_SAME_RANK(IMASK,XSTD_WORK ,XTI_STD)
00315   CALL PACK_SAME_RANK(IMASK,XSKEW_WORK,XTI_SKEW)
00316 !
00317   WRITE(ILUOUT,*) '******************************'
00318   WRITE(ILUOUT,*) 'End Comput Topographic indexes'
00319   WRITE(ILUOUT,*) '******************************'
00320 !
00321 !-------------------------------------------------------------------------------
00322 !
00323   DEALLOCATE(XMIN_WORK )
00324   DEALLOCATE(XMAX_WORK )
00325   DEALLOCATE(XMEAN_WORK)
00326   DEALLOCATE(XSTD_WORK )
00327   DEALLOCATE(XSKEW_WORK)
00328 !
00329 !-------------------------------------------------------------------------------
00330 ENDIF
00331 IF (LHOOK) CALL DR_HOOK('PGD_TOPO_INDEX',1,ZHOOK_HANDLE)
00332 !-------------------------------------------------------------------------------
00333 !
00334 END SUBROUTINE PGD_TOPO_INDEX