Philippe 23/02/2016: lfi2cdf: modif for JPHEXT/=1
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / mesonh / hor_interp_4pts.f90
1 !-----------------------------------------------------------------
2 !--------------- special set of characters for RCS information
3 !-----------------------------------------------------------------
4 ! $Source$ $Revision$ $Date$
5 !-----------------------------------------------------------------
6 !-----------------------------------------------------------------
7 !-----------------------------------------------------------------
8 !     ###########################
9       MODULE MODI_HOR_INTERP_4PTS
10 !     ###########################
11 INTERFACE HOR_INTERP_4PTS
12       SUBROUTINE HOR_INTERP_4PTS_2D(PX1,PY1,PFIELD1,PX2,PY2,PFIELD2)
13       
14 !
15 REAL,   DIMENSION(:),INTENT(IN)       :: PX1       ! x of each grid mesh.
16 REAL,   DIMENSION(:),INTENT(IN)       :: PY1       ! y of each grid mesh.
17 REAL,   DIMENSION(:,:),INTENT(IN)     :: PFIELD1   ! field on grid mesh
18 !
19 REAL,   DIMENSION(:,:),INTENT(IN)     :: PX2       ! x of each new grid mesh.
20 REAL,   DIMENSION(:,:),INTENT(IN)     :: PY2       ! y of each new grid mesh.
21 REAL,   DIMENSION(:,:),INTENT(OUT)    :: PFIELD2   ! field on new grid mesh
22 !
23 END SUBROUTINE HOR_INTERP_4PTS_2D
24 !
25       SUBROUTINE HOR_INTERP_4PTS_3D(PX1,PY1,PFIELD1,PX2,PY2,PFIELD2)
26       
27 !
28 REAL,   DIMENSION(:),INTENT(IN)       :: PX1       ! x of each grid mesh.
29 REAL,   DIMENSION(:),INTENT(IN)       :: PY1       ! y of each grid mesh.
30 REAL,   DIMENSION(:,:,:),INTENT(IN)   :: PFIELD1   ! field on grid mesh
31 !
32 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PX2       ! x of each new grid mesh.
33 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PY2       ! y of each new grid mesh.
34 REAL,   DIMENSION(:,:,:),INTENT(OUT)  :: PFIELD2   ! field on new grid mesh
35 !
36 END SUBROUTINE HOR_INTERP_4PTS_3D
37 END INTERFACE
38 END MODULE MODI_HOR_INTERP_4PTS
39 !
40 !
41 !     ##############################
42       MODULE MODI_HOR_INTERP_4PTS_3D
43 !     ##############################
44 INTERFACE HOR_INTERP_4PTS_3D
45       SUBROUTINE HOR_INTERP_4PTS_3D(PX1,PY1,PFIELD1,PX2,PY2,PFIELD2)
46       
47 !
48 REAL,   DIMENSION(:),INTENT(IN)       :: PX1       ! x of each grid mesh.
49 REAL,   DIMENSION(:),INTENT(IN)       :: PY1       ! y of each grid mesh.
50 REAL,   DIMENSION(:,:,:),INTENT(IN)   :: PFIELD1   ! field on grid mesh
51 !
52 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PX2       ! x of each new grid mesh.
53 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PY2       ! y of each new grid mesh.
54 REAL,   DIMENSION(:,:,:),INTENT(OUT)  :: PFIELD2   ! field on new grid mesh
55 !
56 END SUBROUTINE HOR_INTERP_4PTS_3D
57 END INTERFACE
58 END MODULE MODI_HOR_INTERP_4PTS_3D
59 !
60 !     ##############################################################
61       SUBROUTINE HOR_INTERP_4PTS_3D(PX1,PY1,PFIELD1,PX2,PY2,PFIELD2)
62 !     ##############################################################
63 !
64 !!**** *HOR_INTERP_4PTS* interpolates horizontally a 3D field from a
65 !!                       REGULAR horizontal grid to any other grid
66 !!
67 !!    PURPOSE
68 !!    -------
69 !!
70 !!
71 !!    METHOD
72 !!    ------
73 !!   
74 !!    Bogus value of input field is XUNDEF
75 !!
76 !!    The routine uses only the points with physical values for interpolation:
77 !!       4pts available: interpolations linear in the 2 directions
78 !!       3pts available: plane interpolation
79 !!       2pts available: linear interpolation
80 !!       1pt  available: copy
81 !!
82 !!    Bogus value returned where field could not be interpolated is XUNDEF
83 !!
84 !!    EXTERNAL
85 !!    --------
86 !!
87 !!    IMPLICIT ARGUMENTS
88 !!    ------------------
89 !!
90 !!
91 !!    REFERENCE
92 !!    ---------
93 !!
94 !!    AUTHOR
95 !!    ------
96 !!
97 !!    V. Masson          Meteo-France
98 !!
99 !!    MODIFICATION
100 !!    ------------
101 !!
102 !!    Original    19/03/95
103 !----------------------------------------------------------------------------
104 !
105 !*    0.     DECLARATION
106 !            -----------
107 !
108 !USE MODE_FM
109 !USE MODD_LUNIT
110 USE MODD_PARAMETERS, ONLY: XUNDEF
111 !
112 IMPLICIT NONE
113 !
114 !*    0.1    Declaration of arguments
115 !            ------------------------
116 !
117 REAL,   DIMENSION(:),INTENT(IN)       :: PX1       ! x of each grid mesh.
118 REAL,   DIMENSION(:),INTENT(IN)       :: PY1       ! y of each grid mesh.
119 REAL,   DIMENSION(:,:,:),INTENT(IN)   :: PFIELD1   ! field on grid mesh
120 !
121 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PX2       ! x of each new grid mesh.
122 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PY2       ! y of each new grid mesh.
123 REAL,   DIMENSION(:,:,:),INTENT(OUT)  :: PFIELD2   ! field on new grid mesh
124 !
125 !*    0.2    Declaration of local variables
126 !            ------------------------------
127 !
128 INTEGER                 :: ILUOUT0              ! logical unit
129 INTEGER                 :: IRESP                ! return codes
130 INTEGER                 :: JK
131 INTEGER                 :: IIU,IJU,IIOUT,IJOUT,II,IJ
132 INTEGER                 :: JI,JJ,JIOUT,JJOUT
133 REAL                    :: ZEPS
134 REAL :: ZXA,ZXB,ZXC,ZXD,ZYA,ZYB,ZYC,ZYD,ZA,ZB,ZC,ZD
135 REAL, DIMENSION(3) :: ZX,ZY,ZF
136 INTEGER :: JLOOP
137 REAL :: ZDET,ZALPHA,ZBETA,ZGAMMA
138 !
139 !-------------------------------------------------------------------------------
140 !
141 !*    1.     Initializations
142 !            ---------------
143 !
144 print *,'HOR_INTERP_4PTS: old grid',SIZE(PX1),SIZE(PY1), &
145                          'new grid ',SIZE(PX2,1),SIZE(PY2,2)
146 !CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRESP)
147 IIOUT=SIZE(PX2,1)
148 IJOUT=SIZE(PY2,2)
149 IIU=SIZE(PX1)
150 IJU=SIZE(PY1)
151 ZEPS=1.E-10
152 !
153 !-------------------------------------------------------------------------------
154 !
155 !print * ,' avant boucle JK= k i j iold jold',SIZE(PFIELD1,3) ,IIOUT,IJOUT,SIZE(PX1),SIZE(PY1)
156 !print *, 'PX1(fin), PY1(fin)', PX1(SIZE(PX1)), PY1(SIZE(PY1))
157 !print *, 'PX2(IIOUT,IJOUT), PY2(IIOUT,IJOUT)', PX2(IIOUT,IJOUT), PY2(IIOUT,IJOUT)
158 DO JK=1,SIZE(PFIELD1,3)
159   DO JIOUT=1,IIOUT
160     DO JJOUT=1,IJOUT
161       II=COUNT(PX1(:)<PX2(JIOUT,JJOUT))
162       IJ=COUNT(PY1(:)<PY2(JIOUT,JJOUT))
163       IF ( II<1 .OR. II>=IIU .OR. IJ<1 .OR. IJ>=IJU) THEN
164         PFIELD2(JIOUT,JJOUT,:)=XUNDEF
165         !print *,'pt nouvelle grille hors ancienne grille i j nbi nbj:',JIOUT,JJOUT,II,IJ
166         !print *,'PX2(JIOUT,JJOUT),PY2(JIOUT,JJOUT)' ,PX2(JIOUT,JJOUT),PY2(JIOUT,JJOUT) 
167         CYCLE
168       END IF
169 !
170       !print *,' valeur non indef i j nbi nbj:',JIOUT,JJOUT, II,IJ
171       ZXA=PX1(II)
172       ZXB=PX1(II)
173       ZXC=PX1(II+1)
174       ZXD=PX1(II+1)
175 !
176       ZYA=PY1(IJ)
177       ZYB=PY1(IJ+1)
178       ZYC=PY1(IJ)
179       ZYD=PY1(IJ+1)
180 !
181       ZA=PFIELD1(II,IJ,JK)
182       ZB=PFIELD1(II,IJ+1,JK)
183       ZC=PFIELD1(II+1,IJ,JK)
184       ZD=PFIELD1(II+1,IJ+1,JK)
185 !
186       IF (ALL(ABS(PFIELD1(II:II+1,IJ:IJ+1,JK)-XUNDEF)<ZEPS) ) THEN
187         !print * ,' 4 points a indef  :', PFIELD1(II:II+1,IJ:IJ+1,JK)
188         PFIELD2(JIOUT,JJOUT,JK)=XUNDEF
189         CYCLE
190       ELSE IF (ALL(ABS(PFIELD1(II:II+1,IJ:IJ+1,JK)-XUNDEF)>=ZEPS) ) THEN
191         ZALPHA=ZA+(ZB-ZA)*(PY2(JIOUT,JJOUT)-ZYA)/(ZYB-ZYA)
192         ZBETA =ZC+(ZD-ZC)*(PY2(JIOUT,JJOUT)-ZYC)/(ZYD-ZYC)
193         PFIELD2(JIOUT,JJOUT,JK)=ZALPHA+(ZBETA-ZALPHA)*(PX2(JIOUT,JJOUT)-ZXA)/(ZXC-ZXA)
194       ELSE
195         JLOOP=0
196         DO JI=II,II+1
197           DO JJ=IJ,IJ+1
198             IF (ABS(PFIELD1(JI,JJ,JK)-XUNDEF)>ZEPS) THEN
199               JLOOP=JLOOP+1
200               ZX(JLOOP)=PX1(JI)
201               ZY(JLOOP)=PY1(JJ)
202               ZF(JLOOP)=PFIELD1(JI,JJ,JK)
203             END IF
204           END DO
205         END DO
206         IF (JLOOP==1) THEN
207           PFIELD2(JIOUT,JJOUT,JK)=ZF(1)
208         ELSE IF (JLOOP==2) THEN
209           IF (ABS(ZX(1)-ZX(2))>ZEPS) THEN
210             PFIELD2(JIOUT,JJOUT,JK)=ZF(1)+(ZF(2)-ZF(1))*(PX2(JIOUT,JJOUT)-ZX(1))/(ZX(2)-ZX(1))
211           ELSE
212             PFIELD2(JIOUT,JJOUT,JK)=ZF(1)+(ZF(2)-ZF(1))*(PY2(JIOUT,JJOUT)-ZY(1))/(ZY(2)-ZY(1))
213           END IF
214         ELSE IF (JLOOP==3) THEN
215           ZDET=(ZX(1)-ZX(3))*(ZY(2)-ZY(3))-(ZX(2)-ZX(3))*(ZY(1)-ZY(3))
216           ZALPHA=( (ZF(1)-ZF(3))*(ZY(2)-ZY(3))-(ZF(2)-ZF(3))*(ZY(1)-ZY(3)) )/ZDET
217           ZBETA=-( (ZF(1)-ZF(3))*(ZX(2)-ZX(3))-(ZF(2)-ZF(3))*(ZX(1)-ZX(3)) )/ZDET
218           ZGAMMA=ZF(1)-ZALPHA*ZX(1)-ZBETA*ZY(1)
219           PFIELD2(JIOUT,JJOUT,JK)=ZALPHA*PX2(JIOUT,JJOUT) &
220                                    +ZBETA *PY2(JIOUT,JJOUT) &
221                                    +ZGAMMA
222         END IF
223       END IF
224     END DO
225   END DO
226 END DO
227 print *, 'fin routine HOR_INTERP_4PTS_3D'
228 !-------------------------------------------------------------------------------
229 !
230 !WRITE(ILUOUT0,*) ' Routine HOR_INTERP_4PTS completed'
231 !
232 END SUBROUTINE HOR_INTERP_4PTS_3D
233 !
234 !     ##############################################################
235       SUBROUTINE HOR_INTERP_4PTS_2D(PX1,PY1,PFIELD1,PX2,PY2,PFIELD2)
236 !     ##############################################################
237 !
238 !!**** *HOR_INTERP_4PTS* interpolates horizontally a 2D field from a
239 !!                       REGULAR horizontal grid to any other grid
240 !!
241 !!    PURPOSE
242 !!    -------
243 !!
244 !!
245 !!    METHOD
246 !!    ------
247 !!   
248 !!    Bogus value of input field is XUNDEF
249 !!
250 !!    The routine uses only the points with physical values for interpolation:
251 !!       4pts available: interpolations linear in the 2 directions
252 !!       3pts available: plane interpolation
253 !!       2pts available: linear interpolation
254 !!       1pt  available: copy
255 !!
256 !!    Bogus value returned where field could not be interpolated is XUNDEF
257 !!
258 !!    EXTERNAL
259 !!    --------
260 !!
261 !!    IMPLICIT ARGUMENTS
262 !!    ------------------
263 !!
264 !!
265 !!    REFERENCE
266 !!    ---------
267 !!
268 !!    AUTHOR
269 !!    ------
270 !!
271 !!    V. Masson          Meteo-France
272 !!
273 !!    MODIFICATION
274 !!    ------------
275 !!
276 !!    Original    04/07/96
277 !----------------------------------------------------------------------------
278 !
279 !*    0.     DECLARATION
280 !            -----------
281 !
282 USE MODI_HOR_INTERP_4PTS_3D
283 !
284 IMPLICIT NONE
285 !
286 !*    0.1    Declaration of arguments
287 !            ------------------------
288 !
289 REAL,   DIMENSION(:),INTENT(IN)       :: PX1       ! x of each grid mesh.
290 REAL,   DIMENSION(:),INTENT(IN)       :: PY1       ! y of each grid mesh.
291 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PFIELD1   ! field on grid mesh
292 !
293 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PX2       ! x of each new grid mesh.
294 REAL,   DIMENSION(:,:),  INTENT(IN)   :: PY2       ! y of each new grid mesh.
295 REAL,   DIMENSION(:,:),  INTENT(OUT)  :: PFIELD2   ! field on new grid mesh
296 !
297 !*    0.2    Declaration of local variables
298 !            ------------------------------
299 !
300 REAL, DIMENSION(SIZE(PFIELD1,1),SIZE(PFIELD1,2),1) :: ZFIELD1
301 REAL, DIMENSION(SIZE(PFIELD2,1),SIZE(PFIELD2,2),1) :: ZFIELD2
302 !
303 !-------------------------------------------------------------------------------
304 !
305 ZFIELD1(:,:,1)=PFIELD1(:,:)
306 CALL HOR_INTERP_4PTS_3D(PX1(:),PY1(:),ZFIELD1(:,:,:), &
307                         PX2(:,:),PY2(:,:),ZFIELD2(:,:,:))
308 PFIELD2(:,:)=ZFIELD2(:,:,1)
309 !-------------------------------------------------------------------------------
310 !
311 END SUBROUTINE HOR_INTERP_4PTS_2D