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