Gaelle 2/5/2016 : correction instant M et T
[MNH-git_open_source-lfs.git] / src / MNH / compute_r00.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !--------------- special set of characters for RCS information
7 !-----------------------------------------------------------------
8 ! $Source$ $Revision$
9 ! MASDEV4_7 diag 2006/05/18 13:07:25
10 !-----------------------------------------------------------------
11 !     ###############################
12       MODULE MODI_COMPUTE_R00
13 !     ###############################
14 !
15 INTERFACE
16 SUBROUTINE COMPUTE_R00(HFMFILE)       
17 !
18 CHARACTER (LEN=28), INTENT(IN) :: HFMFILE   ! name of the OUTPUT FM-file
19 !
20 END SUBROUTINE COMPUTE_R00
21 END INTERFACE
22 END MODULE MODI_COMPUTE_R00
23 !
24 !     ###############################
25       SUBROUTINE COMPUTE_R00(HFMFILE)       
26 !     ###############################
27 !
28 !!**** 
29 !!
30 !!    PURPOSE
31 !!    -------
32 !     
33 !!**  METHOD 
34 !!    ------
35 !!    
36 !!    EXTERNAL
37 !!    --------
38 !!
39 !!    IMPLICIT ARGUMENTS
40 !!    ------------------ 
41 !!     MODD_STO_FILE : CFILES
42 !!     MODD_GRID1 : XZZ, XXHAT,XYHAT
43 !!     MODD_LUNIT1: CINIFILE, CLUOUT
44 !!     MODD_FIELD1: XSVM
45 !!     MODD_CONF : NVERB
46 !!     MODD_PARAMETERS : NUNDEF
47 !!     
48 !!    REFERENCE
49 !!    ---------
50 !!      
51 !!    AUTHOR
52 !!    ------
53 !!       F. Gheusi and J. Stein  * Meteo France *
54 !!
55 !!    MODIFICATIONS
56 !!    -------------
57 !!     J. Stein  Jan. 2001  add supplementary starts and spare some memory
58 !!     October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
59 !!                              change of YCOMMENT
60 !!     Mai 2016 (G.Delautier) replace LG?M by LG?T
61 !-------------------------------------------------------------------------------
62 !
63 !*       0.     DECLARATIONS
64 !               ------------
65 !    
66 USE MODD_FIELD_n
67 USE MODD_GRID_n
68 USE MODD_LUNIT_n
69 USE MODD_GRID_n
70 USE MODD_STO_FILE
71 USE MODD_CONF
72 USE MODD_PARAMETERS
73 USE MODD_NSV,            ONLY : NSV_LGBEG,NSV_LGEND
74 !
75 USE MODI_SHUMAN
76 !
77 USE MODD_VAR_ll
78 !
79 USE MODE_FM
80 USE MODE_FMWRIT
81 USE MODE_FMREAD
82 USE MODE_IO_ll
83 USE MODE_ll
84 USE MODD_TYPE_DATE
85 !
86 IMPLICIT NONE
87 !
88 !*       0.1   declarations of arguments
89 !
90 CHARACTER (LEN=28), INTENT(IN) :: HFMFILE   ! name of the OUTPUT FM-file
91 !
92 !*       0.2   declarations of local variables
93 !
94 INTEGER  :: IRESP                ! return code in FM routines
95 INTEGER  :: INPRAR               ! number of articles predicted  in
96                                  !  the LFIFM file
97 INTEGER  :: ININAR               ! number of articles  present in
98                                  !  the LFIFM file
99 INTEGER  :: ITYPE                ! type of file (conv2dia and transfer)
100 !
101 CHARACTER (LEN=100)                :: YCOMMENT
102 CHARACTER (LEN=16)                 :: YRECFM
103 INTEGER                            :: IFILECUR,JFILECUR,NIU,NJU,NKU,IGRID,ILENCH
104 INTEGER                            :: NFILES,JLOOP
105 REAL                               :: ZXOR,ZYOR,ZDX,ZDY
106 REAL                               :: ZSPVAL
107 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZX0, ZY0, ZZ0        ! origin of the 
108        ! particules colocated with the mesh-grid points read in the file
109 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZX00, ZY00, ZZ00, ZZL ! cumulative
110        ! origin for more than one restart of the tracers 
111 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZTH0          ! same fields 
112        ! for Theta as for the coordinates of the origin
113 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZRV0          ! same fields 
114        ! for Rv as for the coordinates of the origin
115 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZWORK1,ZWORK2,ZWORK3       
116 TYPE(DATE_TIME)                    :: TDTCUR_START
117 CHARACTER(LEN=24)                  :: YDATE 
118 INTEGER                            :: IHOUR, IMINUTE
119 REAL                               :: ZSECOND, ZREMAIN
120 LOGICAL                            :: GSTART
121 INTEGER                            :: INBR_START
122 REAL                               :: ZXMAX,ZYMAX,ZZMAX  ! domain extrema
123 INTEGER, DIMENSION(100)            :: NBRFILES 
124 INTEGER                            :: IKU
125 !
126 !-------------------------------------------------------------------------------
127 !
128 !*       1.0    INITIALIZATION
129 !               --------------
130 !
131 ITYPE=2
132 ZSPVAL=-1.E+11
133 IKU=SIZE(XZHAT)
134 !
135 !-------------------------------------------------------------------------------
136 !
137 !*       2.0    FIND THE FILE TO BE TREATED AND THE INIT-SV FILES
138 !               -------------------------------------------------
139 !
140 ! Search the number of the file to be treated
141 IFILECUR=0
142 DO JFILECUR=1,100
143   IF (CINIFILE==CFILES(JFILECUR)) THEN
144     IFILECUR=JFILECUR
145     EXIT
146   END IF
147 END DO
148 !
149 IF (IFILECUR==0) THEN
150   PRINT*,'PROBLEM WITH THE FOLLOWING FILE: ',CINIFILE
151   PRINT*,CFILES
152 !callabortstop
153   CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
154   CALL ABORT
155   STOP
156 ENDIF
157 !
158 ! Search the number of the files(NFILES), where the Lagrangian tracers 
159 !have been reinitialized 
160 NFILES=0
161 DO JFILECUR=IFILECUR+1,100
162   IF (LEN_TRIM(CFILES(JFILECUR)) /= 0 .AND.        &
163       CFILES_STA(JFILECUR) == 'INIT_SV') THEN
164     NFILES= NFILES +1 
165     NBRFILES(NFILES)=JFILECUR       ! contains the number of the files where
166                                     ! the Lag. tracers have been restarted
167   ENDIF
168 END DO
169 !
170 ! compute the number of supplementary cumulative starts
171 INBR_START=1
172 DO JLOOP=1,NFILES-1
173   IF (NSTART_SUPP(JLOOP)/=NUNDEF .AND. NSTART_SUPP(JLOOP)> IFILECUR ) THEN
174     INBR_START=INBR_START+1
175   END IF
176 END DO
177 !
178 !-------------------------------------------------------------------------------
179 !
180 !*       3.0    ALLOCATIONS OF THE ARRAYS AND CONVERSIONS
181 !               -----------------------------------------
182 !
183 !
184 NIU=SIZE(XZZ,1)
185 NJU=SIZE(XZZ,2)
186 NKU=SIZE(XZZ,3)
187 !
188 ALLOCATE(ZX0(NIU,NJU,NKU))
189 ALLOCATE(ZY0(NIU,NJU,NKU))
190 ALLOCATE(ZZ0(NIU,NJU,NKU))
191 ALLOCATE(ZWORK1(NIU,NJU,NKU))
192 ALLOCATE(ZWORK2(NIU,NJU,NKU))
193 ALLOCATE(ZWORK3(NIU,NJU,NKU))
194 ALLOCATE(ZX00(NIU,NJU,NKU))
195 ALLOCATE(ZY00(NIU,NJU,NKU))
196 ALLOCATE(ZZ00(NIU,NJU,NKU))
197 ALLOCATE(ZZL(NIU,NJU,NKU))
198 ALLOCATE(ZTH0(NIU,NJU,NKU))
199 ALLOCATE(ZRV0(NIU,NJU,NKU))
200 ! initial values
201 ZXOR=0.5 * (XXHAT(2)+XXHAT(3)) 
202 ZYOR=0.5 * (XYHAT(2)+XYHAT(3))
203 ZDX= XXHAT(3)-XXHAT(2)
204 ZDY= XYHAT(3)-XYHAT(2)
205 ZZL=MZF(1,IKU,1,XZZ)
206 ZZL(:,:,NKU)=2*XZZ(:,:,NKU)-ZZL(:,:,NKU-1)
207 ZXMAX=ZXOR+(NIU-3)*ZDX
208 ZYMAX=ZYOR+(NJU-3)*ZDY
209 ZZMAX=ZZL(2,2,NKU-1)
210 !  conversion from km to meters
211 ZXOR=ZXOR*1.E-3
212 ZYOR=ZYOR*1.E-3
213 ZDX=ZDX*1.E-3
214 ZDY=ZDY*1.E-3
215 ZZL(:,:,:)=ZZL(:,:,:)*1.E-3
216 ZXMAX=ZXMAX*1.E-3
217 ZYMAX=ZYMAX*1.E-3
218 ZZMAX=ZZMAX*1.E-3
219 !
220 ZX00(:,:,:)=XSVT(:,:,:,NSV_LGBEG)*1.E-3   ! ZX0 in km
221 ZY00(:,:,:)=XSVT(:,:,:,NSV_LGBEG+1)*1.E-3 ! ZY0 in km
222 ZZ00(:,:,:)=XSVT(:,:,:,NSV_LGEND)*1.E-3   ! ZZ0 in km
223 !
224 IF (L2D) THEN
225   WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
226           ZZ00>ZZMAX)
227     ZX00=ZSPVAL
228     ZZ00=ZSPVAL
229   END WHERE
230 ELSE
231   WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
232           ZY00<ZYOR .OR. ZY00>ZYMAX .OR. &
233               ZZ00>ZZMAX)
234     ZX00=ZSPVAL
235     ZY00=ZSPVAL
236     ZZ00=ZSPVAL
237   END WHERE
238 END IF
239 !
240 !-------------------------------------------------------------------------------
241 !
242 !*       4.0    COMPUTE THE ORIGIN STEP BY STEP
243 !               -------------------------------
244 !
245 !
246 ! General loop for the files where a reinitialisation of the tracers 
247 ! is performed
248 DO JFILECUR=1,NFILES
249   !
250   CALL FMOPEN_ll(CFILES(NBRFILES(JFILECUR)),'READ',CLUOUT,   &
251                  INPRAR,ITYPE,NVERB,ININAR,IRESP)
252 !
253 !*       4.1  check if this file is a start instant
254 !
255   GSTART=.FALSE.
256   DO JLOOP=1,NFILES
257     IF (NBRFILES(JFILECUR)==NSTART_SUPP(JLOOP) .OR. JFILECUR==NFILES) THEN
258       INBR_START=INBR_START-1
259       GSTART=.TRUE.
260       EXIT
261     END IF
262   ENDDO
263 !
264 !*       4.2 read the potential temp or the water vapor at the start instant      
265 !
266   IF (GSTART) THEN
267     !
268     YRECFM='DTCUR'
269     CALL FMREAD(CFILES(NBRFILES(JFILECUR)),YRECFM,CLUOUT,'--',TDTCUR_START, &
270                 IGRID,ILENCH,YCOMMENT,IRESP)
271     IHOUR   = INT(TDTCUR_START%TIME/3600.)
272     ZREMAIN = MOD(TDTCUR_START%TIME,3600.)
273     IMINUTE = INT(ZREMAIN/60.)
274     ZSECOND = MOD(ZREMAIN,60.)
275     WRITE(YDATE,FMT='(1X,I4.4,I2.2,I2.2,2X,I2.2,"H",I2.2,"M", &
276          & F5.2,"S")') TDTCUR_START%TDATE, IHOUR,IMINUTE,ZSECOND  
277     !
278     YRECFM='THT'
279     CALL FMREAD(CFILES(NBRFILES(JFILECUR)),YRECFM,CLUOUT,'XY', &
280                 ZTH0(:,:,:),IGRID,ILENCH,YCOMMENT,IRESP)
281     !
282     YRECFM='RVT'
283     CALL FMREAD(CFILES(NBRFILES(JFILECUR)),YRECFM,CLUOUT,'XY', &
284                 ZRV0(:,:,:),IGRID,ILENCH,YCOMMENT,IRESP)
285     !
286     ZRV0(:,:,:)=ZRV0(:,:,:)*1.E+3  ! ZRV0 in g/kg
287     !
288   END IF
289 !
290 !*       4.3  store the X0,Y0,Z0 field for the current start before 
291 !             computing the new origin
292 !
293   IF (GSTART) THEN
294     IGRID=1
295     PRINT *,'INBR_START',INBR_START,' NBRFILES(JFILECUR)',NBRFILES(JFILECUR)
296     WRITE(YRECFM,'(A2,I2.2)')'X0',INBR_START
297     WRITE(YCOMMENT,'(A8,I2.2)')'X_Y_Z_X0',INBR_START
298     YCOMMENT=YCOMMENT(1:10)//YDATE//' (KM)'
299     PRINT *,'YCOMMENT = ',YCOMMENT
300     ILENCH=LEN(YCOMMENT)
301     CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',ZX00(:,:,:),IGRID,ILENCH,   &
302                 YCOMMENT,IRESP)
303     !
304     WRITE(YRECFM,'(A2,I2.2)')'Y0',INBR_START
305     WRITE(YCOMMENT,'(A8,I2.2)')'X_Y_Z_Y0',INBR_START
306     YCOMMENT=YCOMMENT(1:10)//YDATE//' (KM)'
307     PRINT *,'YCOMMENT = ',YCOMMENT
308     ILENCH=LEN(YCOMMENT)
309     CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',ZY00(:,:,:),IGRID,ILENCH,   &
310                 YCOMMENT,IRESP)
311     !
312     WRITE(YRECFM,'(A2,I2.2)')'Z0',INBR_START
313     WRITE(YCOMMENT,'(A8,I2.2)')'X_Y_Z_Z0',INBR_START
314     YCOMMENT=YCOMMENT(1:10)//YDATE//' (KM)'
315     PRINT *,'YCOMMENT = ',YCOMMENT
316     ILENCH=LEN(YCOMMENT)
317     CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',ZZ00(:,:,:),IGRID,ILENCH,   &
318                 YCOMMENT,IRESP)
319   END IF
320 !
321 !
322 !*       4.6   compute and store potential temp and water vapor at the origin
323 !
324   IF (GSTART) THEN
325     !
326     CALL INTERPXYZ(ZX00,ZY00,ZZ00,     &
327                    ZTH0,ZWORK1         )
328     !
329     CALL INTERPXYZ(ZX00,ZY00,ZZ00,     &
330                    ZRV0,ZWORK2         )
331     !
332   END IF
333 !
334   IF (GSTART) THEN
335     !
336     WRITE(YRECFM,'(A3,I2.2)')'TH0',INBR_START
337     WRITE(YCOMMENT,'(A9,I2.2)')'X_Y_Z_TH0',INBR_START
338     YCOMMENT=YCOMMENT(1:11)//YDATE//' (K)'
339     PRINT *,'YCOMMENT = ',YCOMMENT
340     ILENCH=LEN(YCOMMENT)
341     CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',   &
342                 ZWORK1(:,:,:),IGRID,ILENCH,  &
343                 YCOMMENT,IRESP)
344     !
345     WRITE(YRECFM,'(A3,I2.2)')'RV0',INBR_START
346     WRITE(YCOMMENT,'(A9,I2.2)')'X_Y_Z_RV0',INBR_START
347     YCOMMENT=YCOMMENT(1:11)//YDATE//' (G/KG)'
348     PRINT *,'YCOMMENT = ',YCOMMENT
349     ILENCH=LEN(YCOMMENT)
350     CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',   &
351                 ZWORK2(:,:,:),IGRID,ILENCH,  &
352                 YCOMMENT,IRESP)
353   ENDIF
354 !*       4.4   compute the origin of the particules using one more segment
355 !
356   IF (JFILECUR /= NFILES) THEN
357     YRECFM='LGXT'
358     CALL FMREAD(CFILES(NBRFILES(JFILECUR)),YRECFM,CLUOUT,'XY', &
359               ZX0(:,:,:),IGRID,ILENCH,YCOMMENT,IRESP)
360     ZX0(:,:,:)=ZX0(:,:,:)*1.E-3   ! ZX0 in km
361     !
362     YRECFM='LGYT'
363     CALL FMREAD(CFILES(NBRFILES(JFILECUR)),YRECFM,CLUOUT,'XY', &
364               ZY0(:,:,:),IGRID,ILENCH,YCOMMENT,IRESP)
365     ZY0(:,:,:)=ZY0(:,:,:)*1.E-3   ! ZY0 in km
366     !
367     YRECFM='LGZT'
368     CALL FMREAD(CFILES(NBRFILES(JFILECUR)),YRECFM,CLUOUT,'XY', &
369               ZZ0(:,:,:),IGRID,ILENCH,YCOMMENT,IRESP)
370     ZZ0(:,:,:)=ZZ0(:,:,:)*1.E-3   ! ZZ0 in km
371     !
372     ! old position of the set of particles
373     ZWORK1=ZX00
374     ZWORK2=ZY00
375     ZWORK3=ZZ00
376     !
377     IF (L2D) THEN
378       CALL INTERPXYZ(ZWORK1,ZWORK2,ZWORK3,         &
379                      ZX0,ZX00,ZZ0,ZZ00             )
380     ELSE
381       CALL INTERPXYZ(ZWORK1,ZWORK2,ZWORK3,         &
382                      ZX0,ZX00,ZY0,ZY00,ZZ0,ZZ00    )
383     END IF
384     !
385     IF (L2D) THEN
386       WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
387               ZZ00>ZZMAX)
388         ZX00=ZSPVAL
389         ZZ00=ZSPVAL
390       END WHERE
391     ELSE
392       WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
393               ZY00<ZYOR .OR. ZY00>ZYMAX .OR. &
394               ZZ00>ZZMAX)
395         ZX00=ZSPVAL
396         ZY00=ZSPVAL
397         ZZ00=ZSPVAL
398       END WHERE
399     END IF
400     !
401   END IF
402 !
403 !*       4.5   close the input file
404 !
405   CALL FMCLOS_ll(CFILES(NBRFILES(JFILECUR)),'KEEP',CLUOUT,IRESP)
406 !
407 END DO
408 !
409 PRINT*, ' '
410 PRINT*, 'DIAG AFTER ORIGIN COMPUTATIONS AND STORAGE'
411 !
412 !-------------------------------------------------------------------------------
413 !
414 !
415 CONTAINS
416 !
417 !
418 !-------------------------------------------------------------------------------
419 !
420 !
421 SUBROUTINE INTERPXYZ(PX,PY,PZ,PIN1,POUT1,PIN2,POUT2,PIN3,POUT3)
422 !
423 !
424 !*      0. DECLARATIONS
425 !          ------------
426 !
427 !*       0.1  declaration of arguments
428 !
429 REAL, INTENT(IN),  DIMENSION(:,:,:)           :: PX,PY,PZ
430 REAL, INTENT(IN),  DIMENSION(:,:,:)           :: PIN1
431 REAL, INTENT(OUT), DIMENSION(:,:,:)           :: POUT1
432 REAL, INTENT(IN),  DIMENSION(:,:,:), OPTIONAL :: PIN2,PIN3
433 REAL, INTENT(OUT), DIMENSION(:,:,:), OPTIONAL :: POUT2,POUT3   
434 !
435 !*       0.2  declaration of local variables
436 !
437 INTEGER  :: JI,JJ,JK,JKK    ! loop index
438 INTEGER  :: II,IJ,IK        ! grid index for the interpolation
439 REAL     :: ZXREL,ZYREL     ! fractional grid index for the interpolation
440 REAL, DIMENSION(SIZE(PIN1,3)) :: ZZLXY ! vertical grid at the interpolated point
441 REAL     :: ZEPS1,ZEPS2,ZEPS3          ! coeff. for the interpolation
442 REAL     :: ZX,ZY,ZZ
443 LOGICAL  :: GEXT
444 !
445 !-------------------------------------------------------------------------------
446 !
447 DO JK=1,NKU
448   DO JJ=1,NJU
449     DO JI=1,NIU
450       !
451       ZX=PX(JI,JJ,JK) 
452       ZY=PY(JI,JJ,JK)
453       ZZ=PZ(JI,JJ,JK)
454       !
455       ! remove external points
456       IF (L2D) THEN
457         GEXT=(ZX==ZSPVAL).OR.(ZZ==ZSPVAL)
458       ELSE
459         GEXT=(ZX==ZSPVAL).OR.(ZY==ZSPVAL).OR.(ZZ==ZSPVAL)
460       END IF
461       IF (GEXT) THEN
462         POUT1(JI,JJ,JK) = ZSPVAL
463         IF (PRESENT(PIN2)) THEN
464           POUT2(JI,JJ,JK) = ZSPVAL
465         END IF
466         IF (PRESENT(PIN3)) THEN
467           POUT3(JI,JJ,JK) = ZSPVAL
468         ENDIF
469         !
470         CYCLE
471         !
472       END IF
473       !
474       ZXREL=(ZX-ZXOR)/ZDX+2
475       ZYREL=(ZY-ZYOR)/ZDY+2
476       !
477       II=FLOOR(ZXREL)
478       IJ=FLOOR(ZYREL)
479       !
480       ZEPS1=ZXREL-REAL(II)
481       ZEPS2=ZYREL-REAL(IJ)
482       IF (L2D) ZEPS2=0.
483       !
484       DO JKK=1,NKU
485         ZZLXY(JKK)=ZEPS2*(ZEPS1*(ZZL(II+1,IJ+1,JKK))+(1-ZEPS1)*(ZZL(II,IJ+1,JKK)))     &
486              + (1-ZEPS2)*(ZEPS1*(ZZL(II+1,IJ,JKK))+(1-ZEPS1)*(ZZL(II,IJ,JKK)))
487       ENDDO
488       !
489       IK=999
490       DO JKK=2,NKU
491         IF (ZZLXY(JKK).GE.ZZ) THEN
492           IK=JKK-1
493           EXIT 
494         ENDIF
495       ENDDO
496       !
497       IF (IK==999) THEN
498         PRINT*,'PROBLEM AT POINT',II,IJ
499         PRINT*,'XREL, YREL, Z =',ZXREL,ZYREL,ZZ
500         PRINT*,'ZZLXY(NKU)',ZZLXY(NKU)
501 !callabortstop
502         CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
503         CALL ABORT
504         STOP
505       END IF 
506       !
507       ZEPS3=(ZZ-ZZLXY(IK))/(ZZLXY(IK+1)-ZZLXY(IK))
508       !
509       POUT1(JI,JJ,JK) =                                                       & 
510         ZEPS3 *                                                               &
511       (  ZEPS2*(ZEPS1*(PIN1(II+1,IJ+1,IK+1))+(1-ZEPS1)*(PIN1(II,IJ+1,IK+1)))  &
512        + (1-ZEPS2)*(ZEPS1*(PIN1(II+1,IJ,IK+1))+(1-ZEPS1)*(PIN1(II,IJ,IK+1)))  &
513       )                                                                       & 
514       + (1-ZEPS3) *                                                           &
515       (  ZEPS2*(ZEPS1*(PIN1(II+1,IJ+1,IK))+(1-ZEPS1)*(PIN1(II,IJ+1,IK)))      &
516        + (1-ZEPS2)*(ZEPS1*(PIN1(II+1,IJ,IK))+(1-ZEPS1)*(PIN1(II,IJ,IK)))      &
517       )
518       IF (PRESENT(POUT2)) THEN
519         POUT2(JI,JJ,JK) =                                                     & 
520           ZEPS3 *                                                             &
521         (  ZEPS2*(ZEPS1*(PIN2(II+1,IJ+1,IK+1))+(1-ZEPS1)*(PIN2(II,IJ+1,IK+1)))&
522          + (1-ZEPS2)*(ZEPS1*(PIN2(II+1,IJ,IK+1))+(1-ZEPS1)*(PIN2(II,IJ,IK+1)))&
523         )                                                                     & 
524         + (1-ZEPS3) *                                                         &
525         (  ZEPS2*(ZEPS1*(PIN2(II+1,IJ+1,IK))+(1-ZEPS1)*(PIN2(II,IJ+1,IK)))    &
526          + (1-ZEPS2)*(ZEPS1*(PIN2(II+1,IJ,IK))+(1-ZEPS1)*(PIN2(II,IJ,IK)))    &
527         )
528       ENDIF
529         !
530       IF (PRESENT(POUT3)) THEN
531         POUT3(JI,JJ,JK) =                                                     & 
532           ZEPS3 *                                                             &
533         (  ZEPS2*(ZEPS1*(PIN3(II+1,IJ+1,IK+1))+(1-ZEPS1)*(PIN3(II,IJ+1,IK+1)))&
534          + (1-ZEPS2)*(ZEPS1*(PIN3(II+1,IJ,IK+1))+(1-ZEPS1)*(PIN3(II,IJ,IK+1)))&
535         )                                                                     &
536         + (1-ZEPS3) *                                                         &
537         (  ZEPS2*(ZEPS1*(PIN3(II+1,IJ+1,IK))+(1-ZEPS1)*(PIN3(II,IJ+1,IK)))    &
538          + (1-ZEPS2)*(ZEPS1*(PIN3(II+1,IJ,IK))+(1-ZEPS1)*(PIN3(II,IJ,IK)))    &
539         )
540       ENDIF
541       !
542     END DO
543   END DO
544 END DO
545 !
546 !-------------------------------------------------------------------------------
547 !
548 !
549 END SUBROUTINE INTERPXYZ
550 !
551 !-------------------------------------------------------------------------------
552 !
553 END SUBROUTINE COMPUTE_R00