|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE READ_DIRECT(HPROGRAM,HSCHEME,HSUBROUTINE,HFILENAME,HFIELD) 00003 ! ######################################################### 00004 ! 00005 !!**** *READ_DIRECT1* reads a latlon file and call treatment subroutine 00006 !! 00007 !! PURPOSE 00008 !! ------- 00009 !! 00010 !! METHOD 00011 !! ------ 00012 !! 00013 !! EXTERNAL 00014 !! -------- 00015 !! 00016 !! IMPLICIT ARGUMENTS 00017 !! ------------------ 00018 !! 00019 !! REFERENCE 00020 !! --------- 00021 !! 00022 !! AUTHOR 00023 !! ------ 00024 !! 00025 !! V. Masson Meteo-France 00026 !! 00027 !! MODIFICATION 00028 !! ------------ 00029 !! 00030 !! Original 11/09/95 00031 !! 00032 !! V. Masson, March 2010 Optimization of some lat/lon boundaries computations 00033 !---------------------------------------------------------------------------- 00034 ! 00035 !* 0. DECLARATION 00036 ! ----------- 00037 ! 00038 USE MODD_PGD_GRID, ONLY : LLATLONMASK, XMESHLENGTH 00039 ! 00040 USE MODD_ARCH, ONLY : LITTLE_ENDIAN_ARCH 00041 ! 00042 USE MODI_GET_LUOUT 00043 USE MODI_OPEN_NAMELIST 00044 USE MODI_CLOSE_NAMELIST 00045 USE MODI_OPEN_FILE 00046 USE MODI_CLOSE_FILE 00047 USE MODI_READHEAD 00048 USE MODI_INI_SSOWORK 00049 USE MODI_PT_BY_PT_TREATMENT 00050 USE MODE_CHAR2REAL 00051 ! 00052 ! 00053 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00054 USE PARKIND1 ,ONLY : JPRB 00055 ! 00056 USE MODI_ABOR1_SFX 00057 ! 00058 USE MODI_REFRESH_PGDWORK 00059 ! 00060 IMPLICIT NONE 00061 ! 00062 !* 0.1 Declaration of arguments 00063 ! ------------------------ 00064 ! 00065 CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! Type of program 00066 CHARACTER(LEN=6), INTENT(IN) :: HSCHEME ! Scheme treated 00067 CHARACTER(LEN=6), INTENT(IN) :: HSUBROUTINE ! Name of the subroutine to call 00068 CHARACTER(LEN=28), INTENT(IN) :: HFILENAME ! Name of the field file. 00069 CHARACTER(LEN=20), INTENT(IN) :: HFIELD ! Name of the field. 00070 ! 00071 !* 0.2 Declaration of local variables read in the data file head 00072 ! --------------------------------------------------------- 00073 ! 00074 REAL :: ZGLBLATMIN ! minimum latitude of data box in the file 00075 REAL :: ZGLBLONMIN ! minimum longitude of data box in the file 00076 REAL :: ZGLBLATMAX ! maximum latitude of data box in the file 00077 REAL :: ZGLBLONMAX ! maximum longitude of data box in the file 00078 INTEGER :: INBLINE ! number of latitude rows (number of lines 00079 INTEGER :: INBCOL ! number of longitude rows (number of columns) 00080 REAL :: ZNODATA ! value below which data are not considered 00081 ! 00082 !* 0.3 Declaration of local variables 00083 ! ------------------------------ 00084 ! 00085 INTEGER :: IGLB, IGLBHDR ! logical units 00086 INTEGER :: ILUOUT ! output listing logical unit 00087 INTEGER :: IERR ! return codes 00088 ! 00089 CHARACTER(LEN=28) :: YFILENAME ! Name of the field file without header 00090 CHARACTER(LEN=28) :: YFILEHDR ! Name of the field file header 00091 ! 00092 CHARACTER(LEN=7) :: YTYPE ! type of numerical field stored in the 00093 ! ! direct access file ('INTEGER','REAL ') 00094 INTEGER :: IBITS ! number of bits of a record in the 00095 ! ! direct access file (16,32,64) 00096 INTEGER :: JLOOP ! loop index 00097 CHARACTER(LEN=100):: YSTRING ! string 00098 CHARACTER(LEN=88 ):: YSTRING1 ! part of string STRING 00099 INTEGER :: IINDEX ! index of a character in string STRING1 00100 ! 00101 REAL :: ZDLAT ! latitude mesh in the data file 00102 REAL :: ZDLON ! longitude mesh in the data file 00103 INTEGER :: JLINE ! index of line 00104 INTEGER :: JCOL ! index of column 00105 INTEGER, DIMENSION(2) :: ICOL1, ICOL2 ! limits of index of columns 00106 INTEGER :: ILINE1,ILINE2! limits of index of lines 00107 INTEGER :: ICOL ! number of columns in mask domain 00108 INTEGER :: ICOLINDEX ! column index in record 00109 REAL, DIMENSION(:), POINTER :: ZLAT ! latitude of data points 00110 REAL, DIMENSION(:), POINTER :: ZLON ! longitude of data points 00111 ! 00112 INTEGER :: JLON, JLAT ! loop counters 00113 REAL :: ZLONMIN ! minimum longitude of mask mesh 00114 REAL :: ZLONMAX ! maximum longitude of mask mesh 00115 REAL :: ZLATMIN ! minimum latitude of mask mesh 00116 REAL :: ZLATMAX ! maximum latitude of mask mesh 00117 REAL :: ZSHIFT ! shift on longitudes 00118 REAL :: ZFACT ! Factor integer to real 00119 ! 00120 CHARACTER, DIMENSION(:), ALLOCATABLE :: IVALUE8 ! value of a data point 00121 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: IVALUE16 ! value of a data point 00122 CHARACTER(LEN=4), DIMENSION(:), ALLOCATABLE :: IVALUE32R ! value of a data point 00123 INTEGER (KIND=4), DIMENSION(:), ALLOCATABLE :: IVALUE32 ! value of a data point 00124 INTEGER (KIND=8), DIMENSION(:), ALLOCATABLE :: IVALUE64 ! value of a data point 00125 REAL (KIND=4), DIMENSION(:), ALLOCATABLE :: ZVALUE32 ! value of a data point 00126 CHARACTER(LEN=8), DIMENSION(:), ALLOCATABLE :: ZVALUE64 ! value of a data point 00127 ! 00128 REAL, DIMENSION(:), ALLOCATABLE :: ZVALUE ! value of a record of data points 00129 ! 00130 REAL, DIMENSION(:), ALLOCATABLE :: ZVALUE_WORK ! value of a valid data points 00131 REAL, DIMENSION(:), ALLOCATABLE :: ZLAT_WORK ! latitude of a valid data points 00132 REAL, DIMENSION(:), ALLOCATABLE :: ZLON_WORK ! longitude of a valid data points 00133 INTEGER :: IWORK ! index of these data 00134 LOGICAL :: GSWAP ! T: swap has been done 00135 ! 00136 ! 00137 INTEGER :: IRECLENGTH ! record length 00138 INTEGER :: IREC ! record number 00139 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00140 !---------------------------------------------------------------------------- 00141 ! 00142 IF (LHOOK) CALL DR_HOOK('READ_DIRECT',0,ZHOOK_HANDLE) 00143 CALL GET_LUOUT(HPROGRAM,ILUOUT) 00144 ! 00145 !* 1. Openning of global field 00146 ! ------------------------ 00147 ! 00148 !* 1.1 Logical unit attributions 00149 ! ------------------------- 00150 ! 00151 YFILENAME=ADJUSTL(ADJUSTR(HFILENAME)//'.dir') 00152 YFILEHDR =ADJUSTL(ADJUSTR(HFILENAME)//'.hdr') 00153 ! 00154 !* 1.2 Openning of header 00155 ! ------------------ 00156 ! 00157 CALL OPEN_NAMELIST(HPROGRAM,IGLBHDR,YFILEHDR) 00158 ! 00159 !* 1.3 Reading in header of direct access characteristics 00160 ! -------------------------------------------------- 00161 ! 00162 DO JLOOP=1,9 00163 READ(IGLBHDR,'(A100)') YSTRING 00164 IF (YSTRING(1:10)=='recordtype') EXIT 00165 END DO 00166 ! 00167 REWIND(IGLBHDR) 00168 ! 00169 YSTRING1=YSTRING(12:100) 00170 ! 00171 !* string analysis 00172 ! 00173 IINDEX=INDEX(YSTRING1,'n') ! n for integer 00174 IF (IINDEX/=0) THEN 00175 YTYPE='INTEGER' 00176 ELSE 00177 YTYPE='REAL ' 00178 END IF 00179 IINDEX=INDEX(YSTRING1,'8') 00180 IF (IINDEX/=0) IBITS=8 00181 IINDEX=INDEX(YSTRING1,'1') 00182 IF (IINDEX/=0) IBITS=16 00183 IINDEX=INDEX(YSTRING1,'3') 00184 IF (IINDEX/=0) IBITS=32 00185 IINDEX=INDEX(YSTRING1,'4') 00186 IF (IINDEX/=0) IBITS=64 00187 ! 00188 IF(YTYPE=='INTEGER')THEN 00189 IF(HFIELD=='CTI'.OR.HFIELD=='sand fraction'.OR. & 00190 HFIELD=='clay fraction'.OR.HFIELD=='organic carbon')THEN 00191 ZFACT=100.0 00192 ELSEIF (HFIELD=='water depth') THEN 00193 ZFACT=10.0 00194 ELSE 00195 ZFACT=1.0 00196 ENDIF 00197 ELSE 00198 ZFACT=1.0 00199 ENDIF 00200 ! 00201 !---------------------------------------------------------------------------- 00202 ! 00203 !* 2. Reading of the global field 00204 ! --------------------------- 00205 ! 00206 !* 2.1 Head of data file 00207 ! ----------------- 00208 ! 00209 CALL READHEAD(IGLBHDR,ZGLBLATMIN,ZGLBLATMAX,ZGLBLONMIN,ZGLBLONMAX, & 00210 INBLINE,INBCOL,ZNODATA,ZDLAT,ZDLON,ZLAT,ZLON,IERR) 00211 IF (IERR/=0) CALL ABOR1_SFX('READ_DIRECT: PB IN FILE HEADER') 00212 ! 00213 !* 2.2 Closing of header 00214 ! ----------------- 00215 ! 00216 CALL CLOSE_NAMELIST(HPROGRAM,IGLBHDR) 00217 ! 00218 !* 2.3 Dimension of work arrays 00219 ! ------------------------ 00220 ! 00221 ALLOCATE(ZLAT_WORK (INBCOL)) 00222 ALLOCATE(ZLON_WORK (INBCOL)) 00223 ALLOCATE(ZVALUE_WORK(INBCOL)) 00224 ! 00225 !---------------------------------------------------------------------------- 00226 ! 00227 !* 3. Adapt subgrid mesh to input file resolution 00228 ! ------------------------------------------- 00229 ! 00230 IF (HSUBROUTINE=='A_OROG') CALL INI_SSOWORK(XMESHLENGTH,ZDLAT,ZDLON) 00231 ! 00232 !---------------------------------------------------------------------------- 00233 ! 00234 !* 7. Openning of direct access file 00235 ! ------------------------------ 00236 ! 00237 !* 7.1 Record length 00238 ! ------------- 00239 ! 00240 IRECLENGTH = IBITS/8 * INBCOL 00241 ALLOCATE (ZVALUE(INBCOL)) 00242 ALLOCATE (IVALUE8 (INBCOL)) 00243 ALLOCATE (IVALUE16(INBCOL)) 00244 ALLOCATE (IVALUE32(INBCOL)) 00245 ALLOCATE (IVALUE32R(INBCOL)) 00246 ALLOCATE (IVALUE64(INBCOL)) 00247 ALLOCATE (ZVALUE32(INBCOL)) 00248 ALLOCATE (ZVALUE64(INBCOL)) 00249 ! 00250 !* 7.2 Openning of direct access file 00251 ! ------------------------------ 00252 ! 00253 CALL OPEN_FILE(HPROGRAM,IGLB,YFILENAME,'UNFORMATTED', & 00254 HACTION='READ',HACCESS='DIRECT',KRECL=IRECLENGTH ) 00255 ! 00256 !---------------------------------------------------------------------------- 00257 ! 00258 !* 4. loop on mask meshes (lat) 00259 ! ------------------- 00260 ! 00261 GSWAP = .FALSE. 00262 ! 00263 JLAT = 0 00264 ! 00265 DO 00266 ! 00267 JLAT = JLAT + 1 00268 IF (JLAT==361) EXIT 00269 ! 00270 IF ( .NOT. ANY(LLATLONMASK(:,JLAT)) ) CYCLE 00271 ! 00272 ZLATMIN = (JLAT-180)/2. - 0.5 00273 ZLATMAX = (JLAT-180)/2. 00274 ! 00275 !---------------------------------------------------------------------------- 00276 ! 00277 !* 5. index limits on latitude 00278 ! ------------------------ 00279 ! 00280 ILINE1=MAX(MIN(INT((ZGLBLATMAX-ZDLAT/2.-ZLATMAX)/ZDLAT+1.),INBLINE),0)+1 00281 ILINE2=MAX(MIN(INT((ZGLBLATMAX-ZDLAT/2.-ZLATMIN)/ZDLAT+1.),INBLINE),0) 00282 IF ( .NOT. ANY(ZLAT(:)<ZLATMAX .AND. ZLAT(:)>=ZLATMIN) ) CYCLE 00283 ! 00284 !---------------------------------------------------------------------------- 00285 ! 00286 !* 8. Loop on lines 00287 ! ------------- 00288 ! 00289 DO JLINE = ILINE1,ILINE2 00290 ! 00291 !---------------------------------------------------------------------------- 00292 ! 00293 !* 10. Reading in the direct access file 00294 ! --------------------------------- 00295 ! 00296 !* 10.1 Record number 00297 ! ------------- 00298 ! 00299 IREC=JLINE 00300 ! 00301 !* 10.2 Reading the correct data type and conversion into real 00302 ! ------------------------------------------------------ 00303 ! 00304 IF (YTYPE=='INTEGER' .AND. IBITS== 8) THEN 00305 READ(IGLB,REC=IREC) IVALUE8(:) 00306 ZVALUE(:)=IVALUE8(:) 00307 ! negative values are shifted to positive values according to binary coding 00308 WHERE (ZVALUE(:)<0.) ZVALUE(:) = NINT(256.+ZVALUE(:)) 00309 ! 00310 ELSE IF (YTYPE=='INTEGER' .AND. IBITS==16) THEN 00311 READ(IGLB,REC=IREC) IVALUE16(:) 00312 ZVALUE(:)=IVALUE16(:) 00313 IF ( ANY(ABS(ZVALUE)>15000) ) THEN 00314 IF (GSWAP) CALL ABOR1_SFX('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE') 00315 LITTLE_ENDIAN_ARCH = .NOT. LITTLE_ENDIAN_ARCH 00316 GSWAP = .TRUE. 00317 WRITE(ILUOUT,*) '*******************************************************************' 00318 WRITE(ILUOUT,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', & 00319 LITTLE_ENDIAN_ARCH 00320 WRITE(ILUOUT,*) '*******************************************************************' 00321 JLAT = 0 00322 CALL REFRESH_PGDWORK 00323 EXIT ! rereads the file 00324 END IF 00325 00326 ELSE IF (YTYPE=='INTEGER' .AND. IBITS==32) THEN 00327 READ(IGLB,REC=IREC) IVALUE32(:) 00328 ZVALUE(:)=IVALUE32(:) 00329 ELSE IF (YTYPE=='INTEGER' .AND. IBITS==64) THEN 00330 READ(IGLB,REC=IREC) IVALUE64(:) 00331 ZVALUE(:)=IVALUE64(:) 00332 ELSE IF (YTYPE=='REAL ' .AND. IBITS==32) THEN 00333 READ(IGLB,REC=IREC) IVALUE32R(:) 00334 ZVALUE(:)=IVALUE32R(:) 00335 IF ( ANY(ABS(ZVALUE)>0. .AND. ABS(ZVALUE)<1.E-50) & 00336 .OR. ANY(ABS(ZVALUE)>1.E20) ) THEN 00337 IF (GSWAP) CALL ABOR1_SFX('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE') 00338 LITTLE_ENDIAN_ARCH = .NOT. LITTLE_ENDIAN_ARCH 00339 GSWAP = .TRUE. 00340 WRITE(ILUOUT,*) '*******************************************************************' 00341 WRITE(ILUOUT,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', & 00342 LITTLE_ENDIAN_ARCH 00343 WRITE(ILUOUT,*) '*******************************************************************' 00344 JLAT = 0 00345 CALL REFRESH_PGDWORK 00346 EXIT 00347 END IF 00348 ELSE IF (YTYPE=='REAL ' .AND. IBITS==64) THEN 00349 READ(IGLB,REC=IREC) ZVALUE64(:) 00350 ZVALUE(:)=ZVALUE64(:) 00351 IF ( ANY(ABS(ZVALUE)>0. .AND. ABS(ZVALUE)<1.E-50) & 00352 .OR. ANY(ABS(ZVALUE)>1.E20) ) THEN 00353 IF (GSWAP) CALL ABOR1_SFX('READ_DIRECT: SWAP ALREADY DONE, CANNOT BE REDONE') 00354 LITTLE_ENDIAN_ARCH = .NOT. LITTLE_ENDIAN_ARCH 00355 GSWAP = .TRUE. 00356 WRITE(ILUOUT,*) '*******************************************************************' 00357 WRITE(ILUOUT,*) 'Architecture of the machine needs to swap LITTLE_ENDIAN_ARCH to ', & 00358 LITTLE_ENDIAN_ARCH 00359 WRITE(ILUOUT,*) '*******************************************************************' 00360 JLAT = 0 00361 CALL REFRESH_PGDWORK 00362 EXIT 00363 END IF 00364 ELSE 00365 CALL ABOR1_SFX('READ_DIRECT1: DATA TYPE NOT SUPPORTED') 00366 END IF 00367 ! 00368 IF(HFIELD=='CTI')THEN 00369 WHERE(ZVALUE(:)<0.0)ZVALUE(:)=ZNODATA 00370 ENDIF 00371 ! 00372 WHERE(ZVALUE(:)/=ZNODATA)ZVALUE(:)=ZVALUE(:)/ZFACT 00373 ! 00374 !---------------------------------------------------------------------------- 00375 ! 00376 !* 4. loop on mask meshes (lon) 00377 ! ------------------- 00378 ! 00379 DO JLON=1,720 00380 ! 00381 IF (.NOT. LLATLONMASK(JLON,JLAT)) CYCLE 00382 ! 00383 ZLONMIN = JLON /2. - 0.5 00384 ZLONMAX = JLON /2. 00385 ! 00386 !---------------------------------------------------------------------------- 00387 ! 00388 !* 5. limits on longitude 00389 ! ------------------- 00390 ! 00391 !* 5.1 left domain border is set just higher than the global field min. longitude 00392 ! ----------------------------------------------------------------- 00393 ! 00394 ZSHIFT = 360. * NINT((ZLONMIN-ZGLBLONMIN-180.+1.E-10)/360.) 00395 00396 ! 00397 ZGLBLONMIN = ZGLBLONMIN + ZSHIFT 00398 ZGLBLONMAX = ZGLBLONMAX + ZSHIFT 00399 00400 ! 00401 ! 00402 !* 5.2 index limits on longitude 00403 ! ------------------------- 00404 ! 00405 ICOL1(1)=MAX(MIN(INT((ZLONMIN-ZGLBLONMIN-ZDLON/2.)/ZDLON+1.),INBCOL),0)+1 00406 ICOL2(1)=MAX(MIN(INT((ZLONMAX-ZGLBLONMIN-ZDLON/2.)/ZDLON+1.),INBCOL),0) 00407 ! 00408 !* Does right domain border goes outside the global field longitude range? 00409 !* Does it then goes into the global field domain by the other side? 00410 !* Then a second part of the global field domain must be considered 00411 ! 00412 ICOL1(2)=1 00413 ICOL2(2)=MAX(MIN(INT((ZLONMAX-ZGLBLONMIN-ZDLON/2.-360.)/ZDLON+1.),INBCOL),0) 00414 ! 00415 !---------------------------------------------------------------------------- 00416 ! 00417 !* 6. Loop on longitude limits 00418 ! ------------------------ 00419 ! 00420 DO JLOOP=1,2 00421 ! 00422 ICOL = ICOL2(JLOOP) - ICOL1(JLOOP) + 1 00423 ! 00424 IF (ICOL<1) CYCLE 00425 ! 00426 !---------------------------------------------------------------------------- 00427 ! 00428 !* 11. Loop on columns 00429 ! --------------- 00430 ! 00431 IWORK=0 00432 ! 00433 DO JCOL=1,ICOL 00434 ! 00435 !* 11.1 Recovers point value 00436 ! -------------------- 00437 ! 00438 ICOLINDEX = JCOL+ICOL1(JLOOP)-1 00439 00440 ! 00441 !* 11.2 Test with respect to the 'no data' value 00442 ! ---------------------------------------- 00443 ! 00444 IF (ABS(ZVALUE(ICOLINDEX)-ZNODATA)<=1.E-10) CYCLE 00445 ! 00446 ! 00447 !* 11.3 copy of the correct values in a work array 00448 ! ------------------------------------------ 00449 ! 00450 IWORK = IWORK + 1 00451 ZLAT_WORK (IWORK) = ZLAT (JLINE) 00452 ZLON_WORK (IWORK) = ZLON (ICOLINDEX) 00453 ZVALUE_WORK(IWORK) = ZVALUE(ICOLINDEX) 00454 ! 00455 END DO 00456 ! 00457 !------------------------------------------------------------------------------- 00458 ! 00459 !* 12. Call to the adequate subroutine (point by point treatment) 00460 ! ---------------------------------------------------------- 00461 ! 00462 IF (IWORK>0) & 00463 CALL PT_BY_PT_TREATMENT(ILUOUT, ZLAT_WORK(1:IWORK),ZLON_WORK(1:IWORK), & 00464 ZVALUE_WORK(1:IWORK), & 00465 HSUBROUTINE ) 00466 ! 00467 !------------------------------------------------------------------------------- 00468 END DO 00469 !------------------------------------------------------------------------------- 00470 END DO 00471 !------------------------------------------------------------------------------- 00472 END DO 00473 !------------------------------------------------------------------------------- 00474 END DO 00475 ! 00476 !------------------------------------------------------------------------------- 00477 ! 00478 00479 DEALLOCATE(ZLAT) 00480 DEALLOCATE(ZLON) 00481 ! 00482 DEALLOCATE(ZLAT_WORK ) 00483 DEALLOCATE(ZLON_WORK ) 00484 DEALLOCATE(ZVALUE_WORK) 00485 ! 00486 DEALLOCATE (ZVALUE) 00487 DEALLOCATE (IVALUE8 ) 00488 DEALLOCATE (IVALUE16) 00489 DEALLOCATE (IVALUE32) 00490 DEALLOCATE (IVALUE32R) 00491 DEALLOCATE (IVALUE64) 00492 DEALLOCATE (ZVALUE32) 00493 DEALLOCATE (ZVALUE64) 00494 ! 00495 CALL CLOSE_FILE(HPROGRAM,IGLB) 00496 IF (LHOOK) CALL DR_HOOK('READ_DIRECT',1,ZHOOK_HANDLE) 00497 ! 00498 !------------------------------------------------------------------------------- 00499 ! 00500 END SUBROUTINE READ_DIRECT
1.8.0
