|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE INTERPOL_SPLINES(KLUOUT,KCODE,PX,PY,PFIELD,KREP) 00003 ! ######################################################### 00004 ! 00005 !!**** *INTERPOL_SPLINES* interpolates with spline f77 programs a 2D field 00006 !! from all grid points valid values 00007 !! 00008 !! PURPOSE 00009 !! ------- 00010 !! 00011 !! The points are all on only one grid (defined with the coordinates 00012 !! of all the points). The code to apply for each point is: 00013 !! 00014 !! KCODE>0 : data point (with field valid for interpolation) 00015 !! KCODE=-1: point to ignore 00016 !! KCODE=0 : point to interpolate 00017 !! 00018 !! 00019 !! 00020 !! METHOD 00021 !! ------ 00022 !! 00023 !! EXTERNAL 00024 !! -------- 00025 !! 00026 !! IMPLICIT ARGUMENTS 00027 !! ------------------ 00028 !! 00029 !! 00030 !! 00031 !! REFERENCE 00032 !! --------- 00033 !! 00034 !! AUTHOR 00035 !! ------ 00036 !! 00037 !! V. Masson Meteo-France 00038 !! 00039 !! MODIFICATION 00040 !! ------------ 00041 !! 00042 !! Original 19/03/95 00043 !! Modification 00044 !! 25/05/96 (V Masson) W and IW defined in MODD_SPLINESWORK 00045 !! 04/07/96 (V Masson) call with 1d dummy arguments 00046 !! 30/10/96 (V Masson) add deallocations 00047 !! 05/08/97 (V Masson) output listing as dummy argument 00048 !! 17/09/97 (V Masson) routine performs ONLY the splines interpolation 00049 !! 19/12/97 (V Masson) possibility to use the same splines for collocated 00050 !! fields 00051 !---------------------------------------------------------------------------- 00052 ! 00053 !* 0. DECLARATION 00054 ! ----------- 00055 ! 00056 USE MODD_SURF_PAR, ONLY : XUNDEF 00057 ! 00058 USE MODE_SPLINES 00059 ! 00060 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00061 USE PARKIND1 ,ONLY : JPRB 00062 ! 00063 USE MODI_ABOR1_SFX 00064 ! 00065 IMPLICIT NONE 00066 ! 00067 !* 0.1 Declaration of arguments 00068 ! ------------------------ 00069 ! 00070 INTEGER, INTENT(IN) :: KLUOUT ! output listing 00071 INTEGER,DIMENSION(:), INTENT(INOUT) :: KCODE ! code for each point 00072 ! >0 point used for interpolation 00073 ! 0 point to interpolate 00074 ! -1 point not used 00075 ! -2 point not used 00076 ! ! -3 if spline is no computed 00077 ! ! for this point 00078 REAL, DIMENSION(:), INTENT(IN) :: PX ! x of each grid mesh. 00079 REAL, DIMENSION(:), INTENT(IN) :: PY ! y of each grid mesh. 00080 REAL, DIMENSION(:,:),INTENT(INOUT) :: PFIELD ! pgd field on grid mesh. 00081 INTEGER, INTENT(OUT) :: KREP ! error flag 00082 ! 00083 !* 0.2 Declaration of local variables 00084 ! ------------------------------ 00085 ! 00086 INTEGER, PARAMETER :: ND = 2 ! 00087 INTEGER :: IOPT ! 0: splines recomputed 00088 ! ! 1: splines deduced from W and IW 00089 INTEGER :: JK ! loop control 00090 INTEGER :: JIN ! loop control 00091 INTEGER :: JOUT ! loop control 00092 INTEGER :: IVERB = 10 ! verbosity level 00093 INTEGER :: IFIELD ! number of colocated fields 00094 INTEGER, PARAMETER :: IDATAMAX = 2000 ! maximum number of data per blok 00095 ! ! WARNING, if IDATAMAX is modified 00096 ! ! MODIFY nmax in all splines.f routines 00097 INTEGER, PARAMETER :: ISDMAX = 20 ! maximum number of subdomain in one direction 00098 ! ! WARNING, if ISDMAX is modified 00099 ! ! MODIFY ndmax in all splines.f routines 00100 INTEGER, PARAMETER :: IMMAX = 3 ! maximum number of monomes (order 2) 00101 INTEGER :: IDATA ! total number of data 00102 REAL, DIMENSION(2,ND) :: ZXDOM ! coordinates of the limits of the domain 00103 REAL, DIMENSION(:,:), ALLOCATABLE :: ZXDATA ! coordinates of the available data 00104 REAL, DIMENSION(:), ALLOCATABLE :: ZDATA ! available data 00105 INTEGER,DIMENSION(ISDMAX,ISDMAX) :: NDOMDATA ! number of data in a subdomain 00106 00107 INTEGER :: IOUT ! total number of points to interpolate 00108 REAL, DIMENSION(:,:), ALLOCATABLE :: ZXOUT ! coordinates of these points 00109 REAL, DIMENSION(:), ALLOCATABLE :: ZVALOUT ! values of these points 00110 REAL :: ZMINVAL ! minimum value of the field 00111 REAL :: ZMAXVAL ! maximum value of the field 00112 ! 00113 !* 0.3 Declaration for f77 routines 00114 ! ---------------------------- 00115 ! 00116 INTEGER, PARAMETER :: IORDER = 2 ! order of the spline 00117 INTEGER, PARAMETER :: IM = (IORDER*(IORDER+1))/2 ! number of monomes 00118 INTEGER :: ISDI,ISDJ ! number of subdomains in each 00119 ! ! direction for spline computation 00120 INTEGER :: IIORDER 00121 INTEGER :: IIM 00122 INTEGER :: JSDI,JSDJ ! loop controls on subdomains 00123 INTEGER :: IINTER ! inverse ratio of the intersection of 00124 ! ! 2 subdomains over one of the subdomain 00125 REAL, DIMENSION(ND,ND) :: ZG ! matrix of correspondance between 00126 ! ! x and y units 00127 REAL :: ZP, ZS2 00128 ! 00129 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZC ! coefficients of the spline 00130 REAL, DIMENSION(IDATAMAX+3*IMMAX) :: ZWE ! work array 00131 ! 00132 !INTEGER, DIMENSION(2,IM) :: IW 00133 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00134 !------------------------------------------------------------------------------- 00135 ! 00136 !* 1. Miscellaneous Initializations 00137 ! ----------------------------- 00138 ! 00139 IF (LHOOK) CALL DR_HOOK('INTERPOL_SPLINES',0,ZHOOK_HANDLE) 00140 IFIELD=SIZE(PFIELD,2) 00141 ! 00142 !IW = 0 00143 IIORDER=IORDER 00144 IIM=IM 00145 !------------------------------------------------------------------------------- 00146 ! 00147 !* 2. Splines initializations 00148 ! ----------------------- 00149 ! 00150 ! 00151 !* 2.1 Coordinates of input data points 00152 ! -------------------------------- 00153 ! 00154 IDATA=COUNT(KCODE>0) 00155 ALLOCATE(ZXDATA(ND,IDATA)) 00156 ALLOCATE(ZDATA(IDATA)) 00157 ! 00158 ZXDATA(1,:)=PACK(PX(:),MASK=(KCODE(:)>0)) 00159 ZXDATA(2,:)=PACK(PY(:),MASK=(KCODE(:)>0)) 00160 ! 00161 !* 2.2 Coordinates of domain limits 00162 ! ---------------------------- 00163 ! 00164 ZXDOM(1,1)=MINVAL(PX(:),MASK=KCODE(:)>=0) 00165 ZXDOM(2,1)=MAXVAL(PX(:),MASK=KCODE(:)>=0) 00166 ZXDOM(1,2)=MINVAL(PY(:),MASK=KCODE(:)>=0) 00167 ZXDOM(2,2)=MAXVAL(PY(:),MASK=KCODE(:)>=0) 00168 ! 00169 !* 2.3 Coordinates of output interpolated points 00170 ! ----------------------------------------- 00171 ! 00172 IOUT=COUNT(KCODE==0) 00173 ALLOCATE(ZXOUT(ND,IOUT)) 00174 ALLOCATE(ZVALOUT(IOUT)) 00175 ! 00176 ZXOUT(1,:)=PACK(PX(:),MASK=(KCODE(:)==0)) 00177 ZXOUT(2,:)=PACK(PY(:),MASK=(KCODE(:)==0)) 00178 ! 00179 !* 2.4 Matrix of unit correspondance 00180 ! ----------------------------- 00181 ! 00182 ZG(1,1)=1. 00183 ZG(2,2)=1. 00184 ZG(1,2)=0. 00185 ZG(2,1)=0. 00186 ! 00187 !------------------------------------------------------------------------------- 00188 ! 00189 !* 3. Subdomains generation 00190 ! --------------------- 00191 ! 00192 IINTER=4 00193 ! 00194 ISDI=MAX(1,MIN(INT( SQRT ( IDATA/100+1 -1.E-6 )),ISDMAX)) 00195 ISDJ=ISDI 00196 ! 00197 ! 00198 subdomains: DO 00199 DO JSDI=1,ISDI 00200 DO JSDJ=1,ISDJ 00201 NDOMDATA(JSDI,JSDJ)= COUNT( (KCODE(:)>0) .AND. & 00202 ( ( PX(:) >= ZXDOM(1,1) + (JSDI-1-0.5/(IINTER-1))*(ZXDOM(2,1)-ZXDOM(1,1))/ISDI) & 00203 .AND. ( PX(:) <= ZXDOM(1,1) + (JSDI +0.5/(IINTER-1))*(ZXDOM(2,1)-ZXDOM(1,1))/ISDI) & 00204 .AND. ( PY(:) >= ZXDOM(1,2) + (JSDJ-1-0.5/(IINTER-1))*(ZXDOM(2,2)-ZXDOM(1,2))/ISDJ) & 00205 .AND. ( PY(:) <= ZXDOM(1,2) + (JSDJ +0.5/(IINTER-1))*(ZXDOM(2,2)-ZXDOM(1,2))/ISDJ) ) & 00206 ) 00207 END DO 00208 END DO 00209 IF ( (ISDI<ISDMAX) .AND. (ISDJ<ISDMAX) .AND. ANY(NDOMDATA(1:ISDI,1:ISDJ) > IDATAMAX/2) ) THEN 00210 ISDI=ISDI+1 00211 ISDJ=ISDJ+1 00212 CYCLE subdomains 00213 END IF 00214 EXIT subdomains 00215 END DO subdomains 00216 ! 00217 IF ( ANY(NDOMDATA(1:ISDI,1:ISDJ) > IDATAMAX) ) THEN 00218 WRITE(KLUOUT,*) '**********************************************************' 00219 WRITE(KLUOUT,*) '* ERROR during data interpolation *' 00220 WRITE(KLUOUT,*) '* interpolation by splines routines is impossible *' 00221 WRITE(KLUOUT,*) '* please use an input data file with a better resolution *' 00222 WRITE(KLUOUT,*) '**********************************************************' 00223 CALL ABOR1_SFX('INTERPOL_SPLINES: ERROR DURING DATA INTERPOLATION BY SPLINES') 00224 END IF 00225 ! 00226 WRITE(KLUOUT,*) '----------------------------------' 00227 ! 00228 !------------------------------------------------------------------------------- 00229 ! 00230 !* 4. Computation of the spline in the subdomains (f77 routine) 00231 ! ------------------------------------------- 00232 ! 00233 !* 4.1 Output printing 00234 ! --------------- 00235 ! 00236 IF (IVERB>=7) THEN 00237 DO JSDI=1,ISDI 00238 DO JSDJ=1,ISDJ 00239 WRITE(KLUOUT,*) ' domain sdi=',JSDI,' sdj=',JSDJ,' : ', & 00240 NDOMDATA(JSDI,JSDJ),' data points available' 00241 END DO 00242 END DO 00243 END IF 00244 ! 00245 !* 4.2 Choice between splines computation or deduction 00246 ! ----------------------------------------------- 00247 ! 00248 IF (IVERB>=7) WRITE(KLUOUT,*) ' splines computations:' 00249 ! 00250 !* 4.3 Allocations 00251 ! ----------- 00252 ! 00253 ALLOCATE(ZC( IDATAMAX+IMMAX , ISDI , ISDJ )) 00254 IOPT=0 00255 ! 00256 !* 4.4 Loop on fields 00257 ! -------------- 00258 ! 00259 DO JK=1,IFIELD 00260 ! 00261 ZDATA(:)=PACK(PFIELD(:,JK),MASK=(KCODE(:)>0)) 00262 ! 00263 !* 4.5 special case of uniform input data 00264 ! ----------------------------- 00265 ! 00266 ZMINVAL = MINVAL(ZDATA) 00267 ZMAXVAL = MAXVAL(ZDATA) 00268 !------------------------------------------------------------------------------- 00269 IF ( ABS(ZMAXVAL - ZMINVAL) <= 1.E-10 * MAX(ABS(ZMAXVAL),ABS(ZMINVAL)) ) THEN 00270 !------------------------------------------------------------------------------- 00271 ! 00272 ZVALOUT(:)=MAXVAL(ZDATA) 00273 ! 00274 !------------------------------------------------------------------------------- 00275 ELSE 00276 !------------------------------------------------------------------------------- 00277 ! 00278 !* 4.6 Call for splines coefficients 00279 ! ----------------------------- 00280 ! 00281 ! 00282 ZP = 0. 00283 ZS2 = 0. 00284 ! 00285 CALL SPLB2C(IIORDER,IIM,ZXDATA,ZG,ZDATA,ZS2,ZP,0,IOPT,ISDI,ISDJ,IINTER,ZXDOM,ZC,KREP) 00286 ! CALL SPLB2C(IDATA,ZXDATA,ZG,ZDATA,0,0.,ZP,IORDER,IM,IOPT,ISDI,ISDJ,IINTER,ZXDOM, & 00287 ! ZC,KREP,IW,W) 00288 ! 00289 IOPT=1 00290 ! 00291 IF (KREP/=0) THEN 00292 WRITE(KLUOUT,*) 'Routine INTERPOL_SPLINES: error in SPLB2C, KREP=',KREP 00293 WRITE(KLUOUT,*) 'The field is interpolated from 3 nearest points' 00294 IF (ALLOCATED(ZXDATA)) DEALLOCATE(ZXDATA) 00295 IF (ALLOCATED(ZDATA)) DEALLOCATE(ZDATA) 00296 IF (ALLOCATED(ZC)) DEALLOCATE(ZC) 00297 IF (ALLOCATED(ZVALOUT)) DEALLOCATE(ZVALOUT) 00298 IF (ALLOCATED(ZXOUT)) DEALLOCATE(ZXOUT) 00299 IF (LHOOK) CALL DR_HOOK('INTERPOL_SPLINES',1,ZHOOK_HANDLE) 00300 RETURN 00301 END IF 00302 ! 00303 !------------------------------------------------------------------------------- 00304 ! 00305 !* 5. Evaluation of the spline (f77 routine) 00306 ! ------------------------ 00307 ! 00308 ZVALOUT(:)=XUNDEF 00309 ! 00310 IF (IVERB>=7) WRITE(KLUOUT,*) ' splines evaluations' 00311 !CALL SPLB2E(IORDER,IM,ZG,ISDI,ISDJ,ZC,IOUT,ZXOUT,ZVALOUT,ZWE,IW,W) 00312 CALL SPLB2E(IIORDER,IIM,ISDI,ISDJ,ZG,ZC,ZXOUT,ZVALOUT) 00313 ! 00314 !------------------------------------------------------------------------------- 00315 END IF 00316 !------------------------------------------------------------------------------- 00317 ! 00318 !* 6. filling grid points 00319 ! ------------------- 00320 ! 00321 JOUT=0 00322 DO JIN=1,SIZE(PFIELD,1) 00323 IF (KCODE(JIN)/=0) CYCLE 00324 JOUT=JOUT+1 00325 PFIELD(JIN,JK) = ZVALOUT(JOUT) 00326 PFIELD(JIN,JK) = MAX (ZMINVAL, MIN(ZMAXVAL,PFIELD(JIN,JK))) ! control of extrema 00327 END DO 00328 ! 00329 !------------------------------------------------------------------------------- 00330 END DO 00331 !------------------------------------------------------------------------------- 00332 !------------------------------------------------------------------------------- 00333 ! 00334 !* 8. Deallocations 00335 ! ------------- 00336 ! 00337 IF (ALLOCATED(ZXDATA)) DEALLOCATE(ZXDATA) 00338 IF (ALLOCATED(ZDATA)) DEALLOCATE(ZDATA) 00339 IF (ALLOCATED(ZC)) DEALLOCATE(ZC) 00340 IF (ALLOCATED(ZVALOUT)) DEALLOCATE(ZVALOUT) 00341 IF (ALLOCATED(ZXOUT)) DEALLOCATE(ZXOUT) 00342 IF (LHOOK) CALL DR_HOOK('INTERPOL_SPLINES',1,ZHOOK_HANDLE) 00343 ! 00344 !------------------------------------------------------------------------------- 00345 ! 00346 END SUBROUTINE INTERPOL_SPLINES
1.8.0
