Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / DIAPRO / interp_fordiachro.f90
1 !     ######spl
2       MODULE MODI_INTERP_FORDIACHRO
3 !     #############################
4 !
5 INTERFACE
6 !
7 SUBROUTINE INTERP_FORDIACHRO(KLREF,KD,KF,PTAB,PTABREF)
8 REAL,DIMENSION(:,:,:), INTENT(IN)         :: PTAB    
9 REAL,DIMENSION(SIZE(PTAB,1),SIZE(PTAB,2)) :: PTABREF
10 INTEGER          :: KLREF
11 INTEGER          :: KD, KF
12 END SUBROUTINE INTERP_FORDIACHRO
13 !
14 END INTERFACE
15 !
16 END MODULE MODI_INTERP_FORDIACHRO
17 !     ######spl
18       SUBROUTINE INTERP_FORDIACHRO(KLREF,KD,KF,PTAB,PTABREF)
19 !     ######################################################
20 !
21 !!****  *INTERP_FORDIACHRO* - Horizontal cross-section interpolation
22 !!
23 !!    PURPOSE
24 !!    -------
25 !       Interpolates 2D horizontal cross-sections within the Meso-NH 3D
26 !     arrays. These horizontal sections can be:
27 !     -> constant model-level sections (no interpolation, only sampling
28 !        of a particular level);
29 !     -> constant Z (sea-level  altitude) sections;
30 !     -> constant P (hydrostatic pressure) sections 
31 !     -> isentropic (constant potential temperature) 
32 !                                           sections
33 !
34 !!**  METHOD
35 !!    ------
36 !!      
37 !!      Linear interpolation of the model field with
38 !!    respect to "height"  when required
39 !!
40 !!    EXTERNAL
41 !!    --------
42 !!      None
43 !!
44 !!    IMPLICIT ARGUMENTS
45 !!    ------------------
46 !!      Module MODN_NCAR : defines NAM_DIRTRA_POS namelist (former NCAR common)
47 !!         CTYPHOR :  Horizontal cross-section type
48 !!                    (='K' --> model level section;
49 !!                     ='Z' --> constant-altitude section;
50 !!                     ='P' --> isobar section (planned)
51 !!                     ='T' --> isentrope section (planned))
52 !!         XSPVAL  : Special value
53 !!
54 !!      Module MODN_PARA  : Defines NAM_DOMAIN_POS namelist (former PARA common)
55 !!         Module MODD_DIM1 : contains dimensions of data arrays
56 !!                               NKMAX       : z array dimension
57 !!                               NIINF, NISUP: lower and upper bounds of arrays
58 !!                                             to be plotted in x direction
59 !!                               NJINF, NJSUP: lower and upper bounds of arrays
60 !!                                             to be plotted in y direction
61 !! 
62 !!      Module MODD_PARAMETERS : Contains array border depths
63 !!          JPHEXT : Horizontal external points number
64 !!          JPVEXT : Vertical external points number
65 !!         
66 !!      Module MODD_GRID1      : declares grid variables (Model module)
67 !!          XZZ    : true gridpoint z altitude
68 !!
69 !!
70 !!    REFERENCE
71 !!    ---------
72 !!
73 !!      MESO-NH User's Manual, TRACE Post Processing sections, Version 1.0:
74 !!       + Book1: Concepts and Fundamentals, to appear in 1994;
75 !!       + Book2: Technical Reference and Flowcharts, to appear in 1994;
76 !!       + Book3: Tutorial, November 1994.
77 !!
78 !!    AUTHOR
79 !!    ------
80 !!      
81 !!      J. Duron    * Laboratoire d'Aerologie *
82 !!
83 !!    MODIFICATIONS
84 !!    -------------
85 !!      Original       06/06/94
86 !!      Updated   PM   02/12/94
87 !-------------------------------------------------------------------------------
88 !
89 !*       0.    DECLARATIONS
90 !              ------------
91 !
92 USE MODN_NCAR
93 USE MODN_PARA           !NOTICE: MODN_PARA includes MODD_DIM1
94 USE MODD_PARAMETERS
95 USE MODD_MASK3D
96 USE MODD_GRID1
97 USE MODD_TYPE_AND_LH
98 !  07/08/96  !
99 USE MODD_NMGRID
100 USE MODD_PT_FOR_CH_FORDIACHRO
101 !  07/08/96  !
102 USE MODD_RESOLVCAR
103
104 IMPLICIT NONE
105 INTERFACE
106   SUBROUTINE COMPLAT(PLAT)
107   REAL,DIMENSION(:,:) :: PLAT
108   END SUBROUTINE
109 END INTERFACE
110 !
111 !*       0.1   Declaration of arguments and results
112 !
113 REAL,DIMENSION(:,:,:), INTENT(IN)         :: PTAB    !Input arrays where the 
114                                                      !horizontal section is cut 
115 REAL,DIMENSION(SIZE(PTAB,1),SIZE(PTAB,2)) :: PTABREF !Output array containing                                                        !the sampled plane
116 INTEGER :: KLREF  !Sampled level location:
117                   !If CTYPHOR='K'-> model level index given,
118                   !If CTYPHOR='Z'-> sea-level altitude given in meters,
119                   !If CTYPHOR='P'-> pressure level given in hPa,
120                   !If CTYPHOR='T'-> potential temperature level given in K. 
121 INTEGER :: KD, KF ! K Bounds
122 !
123 !*       0.2   Local Variables 
124 !
125 INTEGER :: IID,IJD
126 INTEGER :: II, IJ, IK
127 INTEGER :: JILOOP, JJLOOP, JKLOOP, IKB, IKE
128 REAL    :: ZREF, ZDIXEPS, ZXM, ZXP
129 !  07/08/96  !
130 INTEGER :: IIUP,IJUP,IKU
131 INTEGER :: IND1, IND2
132 REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: ZPTH, ZPTHPROV
133 REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: ZLAT
134 !  07/08/96  !
135 !
136 !-------------------------------------------------------------------------------
137 !
138 !*      1.    Preliminary calculations
139 !             ------------------------
140 !
141 IKB=1+JPVEXT
142 ZDIXEPS=10.*EPSILON(1.)
143 !  07/08/96  !
144 IKU=NKMAX+2*JPVEXT
145 !IKU=SIZE(PTAB,3)
146 IKE=IKU-JPVEXT
147 IIUP=NIMAX+2*JPHEXT
148 !IIUP=SIZE(PTAB,1)
149 !IJUP=SIZE(PTAB,2)
150 IJUP=NJMAX+2*JPHEXT
151 !print *,' IIUP,IJUP,IKU ',IIUP,IJUP,IKU
152 !print *,' INTERP_FORDIACHRO SIZE(XPRES,1),SIZE(XPRES,2),SIZE(XPRES,3) ',SIZE(XPRES,1),SIZE(XPRES,2),SIZE(XPRES,3)
153 !  07/08/96  !
154 IF(LPR)THEN
155   IF(ALLOCATED(ZPTH)) DEALLOCATE(ZPTH)
156   ALLOCATE(ZPTH(SIZE(XPRES,1),SIZE(XPRES,2),SIZE(XPRES,3)))
157 ENDIF
158 IF(LTK .OR. LEV .OR. LSV3)THEN
159   IF(ALLOCATED(ZPTH)) DEALLOCATE(ZPTH)
160   ALLOCATE(ZPTH(SIZE(XTH,1),SIZE(XTH,2),SIZE(XTH,3)))
161 ENDIF
162 if(nverbia > 0)then
163 print *,' INTERP_FORDIACHRO LPR,LTK,LEV,LSV3 ',LPR,LTK,LEV,LSV3
164 endif
165 !
166 ! If not a model level request, convert KLREF to 
167 ! the appropriate variable for interpolation
168 !
169 print *,' *** Interp KLREF, XLOOPZ ',KLREF,XLOOPZ
170 IF(CTYPHOR.EQ.'P')THEN                        ! 'P' requested
171 !>>>>>>>>>>>>>YET TO BE COMPLETED
172 !Mars 2000
173   IF(LCHREEL)THEN
174     ZREF=ALOG10(XLOOPZ*100.)
175   ELSE
176 !Mars 2000
177 !  07/08/96  !
178     ZREF=ALOG10(FLOAT(KLREF)*100.)
179 !  07/08/96  !
180 !Mars 2000
181   ENDIF
182 !Mars 2000
183 ELSE                                          ! 'Z' requested
184 !Mars 2000
185   IF(LCHREEL)THEN
186     ZREF=XLOOPZ
187   ELSE
188 !Mars 2000
189     ZREF=FLOAT(KLREF)
190 !Mars 2000
191   ENDIF
192 !Mars 2000
193 END IF
194 !
195 !-------------------------------------------------------------------------------
196 !
197 !*      2.    Sampling of the requested horizontal section
198 !             --------------------------------------------
199 !
200 !*      2.1   Sampling of a model level: no interpolation necessary
201 !
202 CALL CPSETC('CFT','CONSTANT FIELD - VALUE IS $ZDV$')
203 IF(CTYPHOR.EQ.'K')THEN
204   if(nverbia >0)then
205    print *, ' ** INTERP CTYPHOR.EQ. K, KLREF,KD,KF ', KLREF,KD,KF
206   endif
207   IF(LMSKTOP)THEN
208   if(nverbia >0)then
209   print *,' INTERP MSKTOP NLOOPT ',NLOOPT
210   endif
211     DO JILOOP=NIINF,NISUP
212       DO JJLOOP=NJINF,NJSUP
213         DO JKLOOP=KF,KD,-1
214           IID=JILOOP-NIINF+1
215           IJD=JJLOOP-NJINF+1
216             PTABREF(IID,IJD)=XSPVAL
217             IF(LMASK3(JILOOP,JJLOOP,JKLOOP,NLOOPT))THEN
218               PTABREF(IID,IJD)=PTAB(IID,IJD,JKLOOP-KD+1)
219               EXIT
220             ENDIF
221         ENDDO
222       ENDDO
223     ENDDO
224   if(nverbia >0)then
225      print *,' ** interp CTYPHOR=K AV RETURN DANS LMSKTOP, LPR=',LPR
226   endif
227     RETURN
228   ELSE
229   if(nverbia >0)then
230    print *, ' ** INTERP AV IF(KLREF < KD)THEN KLREF,KD,KF ', KLREF,KD,KF
231   endif
232   IF(KLREF < KD)THEN
233 ! Ajout LKCP Avril 2001 -> prise en compte bilans compresses en K
234   IF(LKCP)THEN
235     PTABREF(:,:)=PTAB(:,:,1)
236   ELSE
237   CALL CPSETC('CFT','UNDER IKB or 1st recorded LEVEL')
238   PTABREF(:,:)=XSPVAL
239   ENDIF
240   RETURN
241   ELSE IF(KLREF > KF)THEN
242   CALL CPSETC('CFT','OVER IKE or last recorded LEVEL')
243   PTABREF(:,:)=XSPVAL
244   RETURN
245   ELSE
246   PTABREF(:,:)=PTAB(:,:,KLREF-KD+1)
247   IND1=0
248   DO JILOOP=1,SIZE(PTABREF,1)
249   DO JJLOOP=1,SIZE(PTABREF,2)
250     IF(PTABREF(JILOOP,JJLOOP) /= XSPVAL)THEN
251       IND1=1
252       EXIT
253     ENDIF
254   ENDDO
255   ENDDO
256   IF(IND1 == 0)THEN
257     CALL CPSETC('CFT','No Value')
258     IND1=0
259   ENDIF
260   RETURN
261   ENDIF
262   ENDIF
263 END IF
264 !  07/08/96  !
265 IF(CTYPHOR.EQ.'P')THEN
266     ZPTH=XPRES(:,:,:,NLOOPT,1,1)
267 ENDIF
268 IF(CTYPHOR.EQ.'T' .OR. CTYPHOR.EQ.'E' .OR. CTYPHOR.EQ.'V')THEN
269     IF(CTYPHOR.EQ.'E')THEN
270       II=SIZE(XTH,1)
271       IJ=SIZE(XTH,2)
272       ALLOCATE(ZLAT(II,IJ))
273       CALL COMPLAT(ZLAT)
274       IK=SIZE(XTH,3)
275 ! 7 Mars 2000
276       print *,' interpol *** NLOOPT ',NLOOPT
277       ZPTH=XTH(:,:,:,NLOOPT,1,1)
278       DO JKLOOP=1,IK
279         WHERE(ZPTH(:,:,JKLOOP) /= XSPVAL)
280           ZPTH(:,:,JKLOOP)=ZPTH(:,:,JKLOOP)* &
281           SIGN(1.,ZLAT(:,:))
282         ENDWHERE
283 !       WHERE(XTH(:,:,JKLOOP,NLOOPT,1,1) /= XSPVAL)
284 !         XTH(:,:,JKLOOP,NLOOPT,1,1)=XTH(:,:,JKLOOP,NLOOPT,1,1)* &
285 !         SIGN(1.,ZLAT(:,:))
286 !       ENDWHERE
287       ENDDO
288       DEALLOCATE(ZLAT)
289     ELSE
290       ZPTH=XTH(:,:,:,NLOOPT,1,1)
291     ENDIF
292 !   ZPTH=XTH(:,:,:,NLOOPT,1,1)
293 ENDIF
294 IF(CTYPHOR == 'T' .OR. CTYPHOR == 'P' .OR. CTYPHOR == 'E' .OR. CTYPHOR.EQ.'V')THEN
295 !Mars 2000
296 IF(LSV3 .OR. LXYZ)THEN
297   IF(ALLOCATED(ZPTHPROV))THEN
298     DEALLOCATE(ZPTHPROV)
299   ENDIF
300   ALLOCATE(ZPTHPROV(SIZE(ZPTH,1),SIZE(ZPTH,2),SIZE(ZPTH,3)))
301   ZPTHPROV=XSPVAL
302 SELECT CASE (NMGRID)
303   CASE(1)
304   CASE(2)
305     WHERE(ZPTH(2:IIUP,:,:) /= XSPVAL .AND. ZPTH(1:IIUP-1,:,:) /= XSPVAL)
306       ZPTHPROV(2:IIUP,:,:)=.5*(ZPTH(2:IIUP,:,:) + ZPTH(1:IIUP-1,:,:))
307     ENDWHERE
308     WHERE(ZPTHPROV(2,:,:) /= XSPVAL .AND. ZPTHPROV(3,:,:) /= XSPVAL)
309       ZPTHPROV(1,:,:)=2.*ZPTHPROV(2,:,:)-ZPTHPROV(3,:,:)
310     ENDWHERE
311     ZPTH=ZPTHPROV
312   CASE(3)
313     WHERE(ZPTH(:,2:IJUP,:) /= XSPVAL .AND. ZPTH(:,1:IJUP-1,:) /= XSPVAL)
314       ZPTHPROV(:,2:IJUP,:)=.5*(ZPTH(:,2:IJUP,:) + ZPTH(:,1:IJUP-1,:))
315     ENDWHERE
316     WHERE(ZPTHPROV(:,2,:) /= XSPVAL .AND. ZPTHPROV(:,3,:) /= XSPVAL)
317       ZPTHPROV(:,1,:)=2.*ZPTHPROV(:,2,:)-ZPTHPROV(:,3,:)
318     ENDWHERE
319     ZPTH=ZPTHPROV
320   CASE(4)
321     WHERE(ZPTH(:,:,2:IKU) /= XSPVAL .AND. ZPTH(:,:,1:IKU-1) /= XSPVAL)
322       ZPTHPROV(:,:,2:IKU)=.5*(ZPTH(:,:,2:IKU) + ZPTH(:,:,1:IKU-1))
323     ENDWHERE
324     WHERE(ZPTHPROV(:,:,2) /= XSPVAL .AND. ZPTHPROV(:,:,3) /= XSPVAL)
325       ZPTHPROV(:,:,1)=2.*ZPTHPROV(:,:,2)-ZPTHPROV(:,:,3)
326     ENDWHERE
327     ZPTH=ZPTHPROV
328   CASE(5)
329     WHERE(ZPTH(2:IIUP,:,:) /= XSPVAL .AND. ZPTH(1:IIUP-1,:,:) /= XSPVAL)
330       ZPTHPROV(2:IIUP,:,:)=.5*(ZPTH(2:IIUP,:,:) + ZPTH(1:IIUP-1,:,:))
331     ENDWHERE
332     WHERE(ZPTHPROV(2,:,:) /= XSPVAL .AND. ZPTHPROV(3,:,:) /= XSPVAL)
333       ZPTHPROV(1,:,:)=2.*ZPTHPROV(2,:,:)-ZPTHPROV(3,:,:)
334     ENDWHERE
335     ZPTH=ZPTHPROV
336     ZPTHPROV=XSPVAL
337     WHERE(ZPTH(:,2:IJUP,:) /= XSPVAL .AND. ZPTH(:,1:IJUP-1,:) /= XSPVAL)
338       ZPTHPROV(:,2:IJUP,:)=.5*(ZPTH(:,2:IJUP,:) + ZPTH(:,1:IJUP-1,:))
339     ENDWHERE
340     WHERE(ZPTHPROV(:,2,:) /= XSPVAL .AND. ZPTHPROV(:,3,:) /= XSPVAL)
341       ZPTHPROV(:,1,:)=2.*ZPTHPROV(:,2,:)-ZPTHPROV(:,3,:)
342     ENDWHERE
343     ZPTH=ZPTHPROV
344   CASE(6)
345     WHERE(ZPTH(:,:,2:IKU) /= XSPVAL .AND. ZPTH(:,:,1:IKU-1) /= XSPVAL)
346       ZPTHPROV(:,:,2:IKU)=.5*(ZPTH(:,:,2:IKU) + ZPTH(:,:,1:IKU-1))
347     ENDWHERE
348     WHERE(ZPTHPROV(:,:,2) /= XSPVAL .AND. ZPTHPROV(:,:,3) /= XSPVAL)
349       ZPTHPROV(:,:,1)=2.*ZPTHPROV(:,:,2)-ZPTHPROV(:,:,3)
350     ENDWHERE
351     ZPTH=ZPTHPROV
352     ZPTHPROV=XSPVAL
353     WHERE(ZPTH(2:IIUP,:,:) /= XSPVAL .AND. ZPTH(1:IIUP-1,:,:) /= XSPVAL)
354       ZPTHPROV(2:IIUP,:,:)=.5*(ZPTH(2:IIUP,:,:) + ZPTH(1:IIUP-1,:,:))
355     ENDWHERE
356     WHERE(ZPTHPROV(2,:,:) /= XSPVAL .AND. ZPTHPROV(3,:,:) /= XSPVAL)
357       ZPTHPROV(1,:,:)=2.*ZPTHPROV(2,:,:)-ZPTHPROV(3,:,:)
358     ENDWHERE
359     ZPTH=ZPTHPROV
360   CASE(7)
361     WHERE(ZPTH(:,:,2:IKU) /= XSPVAL .AND. ZPTH(:,:,1:IKU-1) /= XSPVAL)
362       ZPTHPROV(:,:,2:IKU)=.5*(ZPTH(:,:,2:IKU) + ZPTH(:,:,1:IKU-1))
363     ENDWHERE
364     WHERE(ZPTHPROV(:,:,2) /= XSPVAL .AND. ZPTHPROV(:,:,3) /= XSPVAL)
365       ZPTHPROV(:,:,1)=2.*ZPTHPROV(:,:,2)-ZPTHPROV(:,:,3)
366     ENDWHERE
367     ZPTH=ZPTHPROV
368     ZPTHPROV=XSPVAL
369     WHERE(ZPTH(:,2:IJUP,:) /= XSPVAL .AND. ZPTH(:,1:IJUP-1,:) /= XSPVAL)
370       ZPTHPROV(:,2:IJUP,:)=.5*(ZPTH(:,2:IJUP,:) + ZPTH(:,1:IJUP-1,:))
371     ENDWHERE
372     WHERE(ZPTHPROV(:,2,:) /= XSPVAL .AND. ZPTHPROV(:,3,:) /= XSPVAL)
373       ZPTHPROV(:,1,:)=2.*ZPTHPROV(:,2,:)-ZPTHPROV(:,3,:)
374     ENDWHERE
375     ZPTH=ZPTHPROV
376 END SELECT
377 DEALLOCATE(ZPTHPROV)
378
379 ELSE
380 !Mars 2000
381
382 SELECT CASE (NMGRID)
383   CASE(1)
384   CASE(2)
385     ZPTH(2:IIUP,:,:)=.5*(ZPTH(2:IIUP,:,:) + ZPTH(1:IIUP-1,:,:))
386     ZPTH(1,:,:)=2.*ZPTH(2,:,:)-ZPTH(3,:,:)
387   CASE(3)
388     ZPTH(:,2:IJUP,:)=.5*(ZPTH(:,2:IJUP,:) + ZPTH(:,1:IJUP-1,:))
389     ZPTH(:,1,:)=2.*ZPTH(:,2,:)-ZPTH(:,3,:)
390   CASE(4)
391     ZPTH(:,:,2:IKU)=.5*(ZPTH(:,:,2:IKU) + ZPTH(:,:,1:IKU-1))
392     ZPTH(:,:,1)=2.*ZPTH(:,:,2)-ZPTH(:,:,3)
393   CASE(5)
394     ZPTH(2:IIUP,:,:)=.5*(ZPTH(2:IIUP,:,:) + ZPTH(1:IIUP-1,:,:))
395     ZPTH(1,:,:)=2.*ZPTH(2,:,:)-ZPTH(3,:,:)
396     ZPTH(:,2:IJUP,:)=.5*(ZPTH(:,2:IJUP,:) + ZPTH(:,1:IJUP-1,:))
397     ZPTH(:,1,:)=2.*ZPTH(:,2,:)-ZPTH(:,3,:)
398   CASE(6)
399     ZPTH(:,:,2:IKU)=.5*(ZPTH(:,:,2:IKU) + ZPTH(:,:,1:IKU-1))
400     ZPTH(:,:,1)=2.*ZPTH(:,:,2)-ZPTH(:,:,3)
401     ZPTH(2:IIUP,:,:)=.5*(ZPTH(2:IIUP,:,:) + ZPTH(1:IIUP-1,:,:))
402     ZPTH(1,:,:)=2.*ZPTH(2,:,:)-ZPTH(3,:,:)
403   CASE(7)
404     ZPTH(:,:,2:IKU)=.5*(ZPTH(:,:,2:IKU) + ZPTH(:,:,1:IKU-1))
405     ZPTH(:,:,1)=2.*ZPTH(:,:,2)-ZPTH(:,:,3)
406     ZPTH(:,2:IJUP,:)=.5*(ZPTH(:,2:IJUP,:) + ZPTH(:,1:IJUP-1,:))
407     ZPTH(:,1,:)=2.*ZPTH(:,2,:)-ZPTH(:,3,:)
408 END SELECT
409
410 !Mars 2000
411 ENDIF
412 !Mars 2000
413 !IF(CTYPHOR == 'P')print *,' ZPTH AP MISE SUR GRILLE '
414 ENDIF
415 !  07/08/96  !
416 !
417 !*      2.2   Not a model level request: interpolation necessary
418 !
419 DO JILOOP=NIINF,NISUP
420   DO JJLOOP=NJINF,NJSUP
421
422     IF((CTYPHOR.EQ.'E' .OR. CTYPHOR.EQ.'V') .AND. LINTERPTOP)THEN
423     DO JKLOOP=KF,KD,-1
424 !
425       IID=JILOOP-NIINF+1
426       IJD=JJLOOP-NJINF+1
427 !
428 !*     2.2.3  Potential vorticity request: prepares PV interpolation
429 !
430 !  07/08/96  !
431         ZXM=ZPTH(JILOOP,JJLOOP,JKLOOP)
432         ZXP=ZPTH(JILOOP,JJLOOP,MIN(KF,JKLOOP+1))
433 !  07/08/96  !
434 !
435 !*     2.3    Selects points within the TRACE display window
436 !
437 !
438 !  07/08/96  !
439       PTABREF(IID,IJD)=XSPVAL
440 !  18/02/2000 Essai pour prise en compte des valeurs speciales
441       IF(LSV3 .AND. LXYZ00)THEN
442         IF(ZXP == XSPVAL .OR. ZXM  == XSPVAL .OR. (ZXP == XSPVAL .AND. &
443           ZXM  == XSPVAL))THEN
444           if(nverbia == 20)then
445           print *,' ***interp JILOOP JJLOOP JKLOOP ZXP ZXM ',JILOOP,&
446           JJLOOP,JKLOOP,ZXP,ZXM
447           endif
448           CYCLE
449         ENDIF
450       ENDIF
451 !  18/02/2000 Essai pour prise en compte des valeurs speciales
452       IF((ZXP-ZREF)*(ZREF-ZXM).GE.0.)THEN
453         IF(JKLOOP+1 <= IKB .OR. JKLOOP+1 > IKE)THEN
454           CYCLE
455         ELSE
456           GO TO 4
457         ENDIF
458       ELSE IF(ZXP.GE.ZXM-ZDIXEPS.AND.ZXP.LE.ZXM+ZDIXEPS.AND.  &
459       ZREF.GE.ZXM-ZDIXEPS.AND.ZREF.LE.ZXM+ZDIXEPS)THEN
460         IF(JKLOOP+1 <= IKB .OR. JKLOOP+1 > IKE)THEN
461           CYCLE
462         ELSE
463           GO TO 4
464         ENDIF
465       ENDIF
466 !  07/08/96  !
467 !
468     ENDDO
469
470     ELSE
471
472     DO JKLOOP=KD,KF
473 !
474       IID=JILOOP-NIINF+1
475       IJD=JJLOOP-NJINF+1
476 !
477 !*     2.2.1  Pressure level request: prepares Log(P) interpolation
478 !
479       IF(CTYPHOR.EQ.'P')THEN
480 !>>>>>>>>>>>>YET TO BE DEVELOPED
481         ZXM=ALOG10(ZPTH(JILOOP,JJLOOP,JKLOOP))
482         ZXP=ALOG10(ZPTH(JILOOP,JJLOOP,MIN(KF,JKLOOP+1)))
483 !
484 !*     2.2.2  Altitude level request: prepares Z interpolation
485 !
486       ELSE IF (CTYPHOR.EQ.'Z')THEN
487         ZXM=XZZ(JILOOP,JJLOOP,JKLOOP)
488         ZXP=XZZ(JILOOP,JJLOOP,MIN(KF,JKLOOP+1))
489 !
490 !*     2.2.3  Potential temperature request: prepares Theta interpolation
491 !
492       ELSE IF(CTYPHOR.EQ.'T')THEN
493 !>>>>>>>>>>>>YET TO BE DEVELOPED
494 !  07/08/96  !
495         ZXM=ZPTH(JILOOP,JJLOOP,JKLOOP)
496         ZXP=ZPTH(JILOOP,JJLOOP,MIN(KF,JKLOOP+1))
497 !  07/08/96  !
498 ! Mars 2000 Ajout possibilite de faire interpolation a partir du bas
499 ! pour la vorticite potentielle et SV3
500       ELSE IF(CTYPHOR.EQ.'E' .AND. .NOT.LINTERPTOP)THEN
501         ZXM=ZPTH(JILOOP,JJLOOP,JKLOOP)
502         ZXP=ZPTH(JILOOP,JJLOOP,MIN(KF,JKLOOP+1))
503       ELSE IF(CTYPHOR.EQ.'V' .AND. .NOT.LINTERPTOP)THEN
504         ZXM=ZPTH(JILOOP,JJLOOP,JKLOOP)
505         ZXP=ZPTH(JILOOP,JJLOOP,MIN(KF,JKLOOP+1))
506 ! Mars 2000 
507       END IF
508 !
509 !*     2.3    Selects points within the TRACE display window
510 !
511 !
512 !  07/08/96  !
513       PTABREF(IID,IJD)=XSPVAL
514 !  23/03/2000 Essai pour prise en compte des valeurs speciales
515       IF(LSV3 .AND. LXYZ00)THEN
516         IF(ZXP == XSPVAL .OR. ZXM  == XSPVAL .OR. (ZXP == XSPVAL .AND. &
517           ZXM  == XSPVAL))THEN
518           if(nverbia == 20)then
519           print *,' ***interp JILOOP JJLOOP JKLOOP ZXP ZXM ',JILOOP,&
520           JJLOOP,JKLOOP,ZXP,ZXM
521           endif
522           CYCLE
523         ENDIF
524       ENDIF
525 !  23/03/2000 Essai pour prise en compte des valeurs speciales
526       IF((ZXP-ZREF)*(ZREF-ZXM).GE.0.)THEN
527         IF(JKLOOP+1 <= IKB .OR. JKLOOP+1 > IKE)THEN
528           CYCLE
529         ELSE
530           GO TO 4
531         ENDIF
532       ELSE IF(ZXP.GE.ZXM-ZDIXEPS.AND.ZXP.LE.ZXM+ZDIXEPS.AND.  &
533       ZREF.GE.ZXM-ZDIXEPS.AND.ZREF.LE.ZXM+ZDIXEPS)THEN
534         IF(JKLOOP+1 <= IKB .OR. JKLOOP+1 > IKE)THEN
535           CYCLE
536         ELSE
537           GO TO 4
538         ENDIF
539       ENDIF
540 !  07/08/96  !
541 !
542     ENDDO
543     ENDIF
544 !
545 !*    2.4    Out of display window: inserts a NCAR special value 
546 !            to suppress display
547 !
548     PTABREF(IID,IJD)=XSPVAL
549     GO TO 5
550 !
551 4   CONTINUE
552 !
553 !*    2.5   Requested level colocated with a model level: no interpolation
554
555     IF(ZXP==ZXM)THEN
556       PTABREF(IID,IJD)=PTAB(IID,IJD,JKLOOP-KD+1)
557 !       print *,' INTERP_FORDIACHRO ZXM ZXP ',ZXM,ZXP
558 !     IF(CTYPHOR == 'P')THEN
559 !       print *,' CAS ZXM=ZXP '
560 !     ENDIF 
561 !
562 !*    2.6   Requested level located between model levels: linear interpolation
563 !
564     ELSE
565       SELECT CASE(CTYPHOR)
566         CASE('Z')
567 !         print *,' ZXP - ZXM ',ZXP-ZXM
568           PTABREF(IID,IJD)=(PTAB(IID,IJD,JKLOOP-KD+1)*(ZXP-ZREF)+  &
569           PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1))*  &
570           (ZREF-ZXM))/MAX(1.E-8,(ZXP-ZXM))
571         CASE('T','E','V')
572 !         print *,' ZXP - ZXM ',ZXP-ZXM
573           LTHSTAB=.TRUE.
574           IF(JKLOOP+1 > IKB)THEN
575             IF(ZXP-ZXM >= 0.)THEN
576               LTHSTAB=.TRUE.
577             ELSE
578               LTHSTAB=.FALSE.
579 !             print *,' JKLOOP, ZXP, ZXM ',JKLOOP,ZXP,ZXM
580           ENDIF
581           ENDIF
582           PTABREF(IID,IJD)=(PTAB(IID,IJD,JKLOOP-KD+1)*(ZXP-ZREF)+  &
583           PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1))*  &
584           (ZREF-ZXM))/(ZXP-ZXM)
585         CASE('P')
586           PTABREF(IID,IJD)=(PTAB(IID,IJD,JKLOOP-KD+1)*(ZXP-ZREF)+  &
587           PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1))*  &
588           (ZREF-ZXM))/MIN(-1.E-8,(ZXP-ZXM))
589       END SELECT
590 !     IF(CTYPHOR == 'P' .AND. IID == 4 .AND. IJD == 8)THEN
591 !       print *,' IID,IJD,JKLOOP-KD+1,PTAB,ZXP-ZREF ',IID,IJD,JKLOOP-KD+1,PTAB(IID,IJD,JKLOOP-KD+1),ZXP-ZREF
592 !       print *,' IID,IJD,JKLOOP-KD+1,PTAB,ZXP-ZREF ZXP-ZXM',IID,IJD,JKLOOP-KD+1,PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP-KD+1+1)),ZREF-ZXM,ZXP-ZXM
593 !     ENDIF
594     END IF
595 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
596 ! Sept 2000 test suivant supprime
597 !   IF(LXYZ .OR. LSV3)THEN
598     IF(PTAB(IID,IJD,JKLOOP-KD+1) == XSPVAL .AND. &
599       PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1)) == XSPVAL)THEN
600       PTABREF(IID,IJD)=XSPVAL
601     ELSE IF(PTAB(IID,IJD,JKLOOP-KD+1) /= XSPVAL .AND. &
602       PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1)) == XSPVAL)THEN
603       PTABREF(IID,IJD)=XSPVAL
604 !     PTABREF(IID,IJD)=PTAB(IID,IJD,JKLOOP-KD+1)
605     ELSE IF(PTAB(IID,IJD,JKLOOP-KD+1) == XSPVAL .AND. &
606       PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1)) /= XSPVAL)THEN
607 !     PTABREF(IID,IJD)=PTAB(IID,IJD,MIN(KF-KD+1,JKLOOP+1-KD+1))
608       PTABREF(IID,IJD)=XSPVAL
609     ENDIF
610 !   ENDIF
611 ! Sept 2000
612 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
613 !
614 5 CONTINUE
615 !
616   ENDDO
617 ENDDO
618 !
619 IND1=0
620 IND2=0
621 DO JILOOP=1,SIZE(PTABREF,1)
622   DO JJLOOP=1,SIZE(PTABREF,2)
623     IF(PTABREF(JILOOP,JJLOOP) /= XSPVAL)THEN
624       IND1=1
625       EXIT
626     ELSE
627       IND2=1
628     ENDIF
629   ENDDO
630 ENDDO
631 !print *,' PTABREF 1-8 '
632 !DO JJLOOP=1,SIZE(PTABREF,2)
633 !  print *,(PTABREF(JILOOP,JJLOOP),JILOOP=1,8)
634 !ENDDO
635 !print *,' PTABREF 9-16 '
636 !DO JJLOOP=1,SIZE(PTABREF,2)
637 !  print *,(PTABREF(JILOOP,JJLOOP),JILOOP=9,16)
638 !ENDDO
639 !print *,' PTABREF 17-24 '
640 !DO JJLOOP=1,SIZE(PTABREF,2)
641 !  print *,(PTABREF(JILOOP,JJLOOP),JILOOP=17,24)
642 !ENDDO
643 IF(IND1 == 0 .AND. IND2 /= 0)THEN
644   CALL CPSETC('CFT','<IKB or 1st recorded LEVEL or >IKE LEVEL')
645   IF(LSV3 .OR. LXYZ)THEN
646     CALL CPSETC('CFT','No value')
647   ENDIF
648 ELSE 
649   CALL CPSETC('CFT','CONSTANT FIELD - VALUE IS $ZDV$')
650 ENDIF
651 if(nverbia > 0)then
652   print *,' INTERP_FORDIACHRO end: LPR,LTK,LEV,LSV3 ',LPR,LTK,LEV,LSV3
653 endif
654     
655 !
656 !----------------------------------------------------------------------------
657 !
658 !*     3.    EXIT
659 !            ----
660 !
661 RETURN
662 END SUBROUTINE INTERP_FORDIACHRO