Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / EXTRACTDIA / compute_r00_pc.f90
1 PROGRAM COMPUTE_R00       
2 !     ###############################
3 !
4 ! ce programme est la version PC du programme compute_r00.f90 de mesoNH utilisee
5 ! dans DIAG pouvant tourner sur PC seul afin de pouvoir se passer du
6 ! super-calculateur pour reconstituer a loisir des lachers de particules
7 ! arbitraires.
8 !
9 ! on garde la structure Fortran 90 et les routines d'interpolation mais on 
10 ! saisit les noms des fichiers a travers le fichier compute_r00.nam
11 !-------------------------------------------------------------------------------
12 !
13 !*       0.     DECLARATIONS
14 !               ------------
15 !
16 !  modules commun pour la lecture et l'ecriture
17 !
18 !                    NIMAX,NJMAX,NKMAX, NIINF, NISUP
19 USE MODD_DIM1
20 !                    grille : XXDXHAT(:,1:7) et XXX(:,1:7), XXZS(:,:,1:7)
21 USE MODD_COORD
22 !                    ref grille: XLON0,XLAT0,XBETA,XRPK
23 USE MODD_GRID
24 !                    descriptif grille: XXHAT(:) ,XLAT(:,:),XDXHAT(:),XMAP(:,:)
25 !                    ,XZS(:,:),XZZ(:,:,:) ,XCOSSLOPE(:,:),XDIRCOSXW(:,:)
26 USE MODD_GRID1
27 !                    XVAR(i,j,k,,,), XMASK,XTRAJ ,XDATIME(16,t)
28 USE MODD_ALLOC_FORDIACHRO
29 !                    NBFILES + nom des fichiers CFILEDIAS, CLUOUTDIAS
30 USE MODD_DIACHRO, ONLY:CFILEDIA,CLUOUTDIA, &
31                    NLUOUTDIA,NRESPDIA,NNPRARDIA,NFTYPEDIA,NVERBDIA, NNINARDIA
32 !
33 USE MODI_WRITEVAR
34 !
35 IMPLICIT NONE
36 !
37 TYPE DATE
38 INTEGER :: YEAR
39 INTEGER :: MONTH
40 INTEGER :: DAY
41 END TYPE DATE
42 !
43 TYPE DATE_TIME
44 TYPE (DATE) :: TDATE
45 REAL :: TIME
46 END TYPE DATE_TIME 
47 !
48 !
49 CHARACTER (LEN=28) :: HFMFILE   ! name of the OUTPUT FM-file
50 CHARACTER (LEN=31) :: HFMFILE_sto 
51 !
52 !*       0.2   declarations of local variables
53 !
54 INTEGER  :: IRESP                ! return code in FM routines
55 INTEGER  :: INPRAR               ! number of articles predicted  in
56                                  !  the LFIFM file
57 INTEGER  :: ININAR               ! number of articles  present in
58                                  !  the LFIFM file
59 INTEGER  :: ITYPE                ! type of file (conv2dia and transfer)
60 !
61 ! **** la longueur du nom ne doit pas depasser 13 car. si le fichier
62 ! contient des groupes a un seul PROCessus, ou 9 si plusieurs PROCessus ****
63 CHARACTER (LEN=13)                 :: YRECFM
64 CHARACTER (LEN=100)                :: YCOMMENT
65 !
66 INTEGER                            :: IFILECUR,JFILECUR,NIU,NJU,NKU,IGRID,ILENCH
67 INTEGER                            :: NFILES,JLOOP
68 REAL                               :: ZXOR,ZYOR,ZDX,ZDY
69 REAL                               :: ZSPVAL
70 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZX0, ZY0, ZZ0        ! origin of the 
71        ! particules colocated with the mesh-grid points read in the file
72 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZX00, ZY00, ZZ00, ZZL ! cumulative
73        ! origin for more than one restart of the tracers 
74 REAL, ALLOCATABLE, DIMENSION(:,:,:,:):: ZWORK 
75 TYPE(DATE_TIME)                    :: TDTCUR_START
76 CHARACTER(LEN=24)                  :: YDATE 
77 INTEGER                            :: IHOUR, IMINUTE
78 REAL                               :: ZSECOND, ZREMAIN
79 LOGICAL                            :: GSTART
80 INTEGER                            :: INBR_START
81 REAL                               :: ZXMAX,ZYMAX,ZZMAX  ! domain extrema
82 INTEGER, DIMENSION(100)            :: NBRFILES
83 !  declarations supplementaires
84 INTEGER                            :: iret     ! code de retour de lecture  
85 CHARACTER (LEN=3), SAVE :: CNAME_SUP
86 !-----------------------------------------------------------------------
87 !  definitions des noms de fichiers venant de modd_sto_file de Meso-NH
88 !  et definition de la namelist prise dans diag.f90
89 CHARACTER (LEN=28), SAVE :: CFILES(100)       ! names of the files to be treated
90 CHARACTER (LEN=28), SAVE :: CFILES_STA(100)   ! status of these files 'INIT_SV'
91                                               ! if a restart of the lagrangian
92                                               ! tracers has been performed
93 INTEGER           , SAVE :: NSTART_SUPP(100)  ! supplementary starts 
94                                               ! for the lagrangian trajectories
95 !-----------------------------------------------------------------------
96 !                                              
97 !  article supplementaire
98 CHARACTER (LEN=28), SAVE :: CFIELD_LAG(100)   ! tableau de noms de record devant
99 CHARACTER (LEN=100),DIMENSION(:),ALLOCATABLE  :: YUNITE
100 ! etre etudies lagrangiennement THM RVM RRM...
101 INTEGER                  :: NUNDEF, inbr_field, NFILES_tot, k, ifield
102 LOGICAL                  :: L2D
103 CHARACTER (LEN=3), SAVE :: CFLAGFILE
104 !
105 !-----------------------------------------------------------------------
106 !
107 NAMELIST/NAM_STO_FILE/ CFILES, NSTART_SUPP
108 !-----------------------------------------------------------------------
109 !
110 NAMELIST/NAM_FIELD /CFIELD_LAG
111 !
112 !-------------------------------------------------------------------------------
113 !
114 !*       1.0    Lecture des noms des fichiers et initialisation
115 !               -----------------------------------------------
116 !
117 ! ouverture du fichier contenant les noms des fichiers diachroniques a
118 ! traiter
119 !
120 ! ecrire les fichiers dans le meme ordre que pour DIAG1.nam (cf doc Gheusi +
121 !   Stein) dans NAM_STO_FILES i.e. ordre inverse chrono
122 !
123 !
124 open (unit=104,FILE='compute_r00.nam',FORM='FORMATTED')
125 !
126 !
127 nverbdia=1
128 ITYPE=2
129 ZSPVAL=-1.E+11
130 NUNDEF=-9999
131 CFILES(:) = '                         '
132 NSTART_SUPP(:) = NUNDEF
133 CFILES_STA(:) = 'INIT_SV'
134 CFIELD_LAG(:) = '                         '
135 CNAME_SUP='SAM'
136 !
137 READ(104,NML=NAM_STO_FILE)
138 !
139 READ(104,NML=NAM_FIELD)
140 !
141 close(104)
142 !
143 !-------------------------------------------------------------------------------
144 !
145 !*       2.0    FIND THE FILE TO BE TREATED AND THE INIT-SV FILES
146 !               -------------------------------------------------
147 !
148 !
149 ! determination du nombre de champs de var a traiter lagrangiennement
150 inbr_field=0
151 DO JLOOP=1,100
152   IF (LEN_TRIM(CFIELD_LAG(JLOOP))/= 0) THEN
153     inbr_field=inbr_field+1
154   END IF
155 END DO
156 !
157 !
158 ! recherche du nombre total de fichier a traiter
159 NFILES_tot=0
160 DO JFILECUR=1,100
161   IF (LEN_TRIM(CFILES(JFILECUR)) /= 0 ) THEN
162     NFILES_tot= NFILES_tot +1 
163   ENDIF
164 END DO
165 !
166 ! ouverture des fichiers
167 do jfilecur=1,NFILES_tot
168   CFLAGFILE='OPE'
169   CALL READVAR('ZSBIS',CFILES(jfilecur),CFLAGFILE,nverbdia,iret)
170 end do
171 !
172 if (nverbdia>0) then
173   print *,'nbre de fichiers a traiter',NFILES_tot
174   print *,'nombre de champs de var a traiter lagrangiennement',inbr_field
175 end if
176 !
177 !**************************************************
178 !**************************************************
179 ! pour coller a la version MESONH je prends les memes noms
180 ! cette boucle correspond au traitement de diag pour chacun des fichiers 
181 ! a traiter
182 do ifilecur=1,NFILES_tot
183 HFMFILE=CFILES(ifilecur)
184 print *,'fichier traite HFMFILE = ',HFMFILE
185 !**************************************************
186 !**************************************************
187 !   rem on n'indente pas la boucle ifilecur
188 !pour garder le code commun avec compute_r00.f90 sur VPP
189 !
190 !
191 ! Search the number of the files(NFILES), where the Lagrangian tracers 
192 !have been reinitialized 
193 NFILES=0
194 DO JFILECUR=IFILECUR+1,100
195   IF (LEN_TRIM(CFILES(JFILECUR)) /= 0 .AND.        &
196       CFILES_STA(JFILECUR) == 'INIT_SV') THEN
197     NFILES= NFILES +1 
198     NBRFILES(NFILES)=JFILECUR       ! contains the number of the files where
199                                     ! the Lag. tracers have been restarted
200   ENDIF
201 END DO
202 !
203 ! compute the number of supplementary cumulative starts
204 INBR_START=1
205 DO JLOOP=1,NFILES-1
206   IF (NSTART_SUPP(JLOOP)/=NUNDEF .AND. NSTART_SUPP(JLOOP)> IFILECUR ) THEN
207     INBR_START=INBR_START+1
208   END IF
209 END DO
210 !
211 if (nverbdia >0) then
212   print *,'INBR_START = ',INBR_START,' pour le fichier ',IFILECUR
213 end if
214 !-------------------------------------------------------------------------------
215 !
216 !*       3.0    ALLOCATIONS OF THE ARRAYS AND CONVERSIONS
217 !               -----------------------------------------
218 !
219 NIU=SIZE(XZZ,1)
220 NJU=SIZE(XZZ,2)
221 NKU=SIZE(XZZ,3)
222 if (nju==3) then
223   L2D=.TRUE.
224 else
225   L2D=.FALSE.
226 end if
227 if (nverbdia >0) print *,'L2D = ',L2D
228 !
229 if (.NOT. allocated(ZX0)) then  ! pas d'indentation pour garder la possibilite 
230                                ! de faire un diff des compute_r00
231 ALLOCATE(ZX0(NIU,NJU,NKU))
232 ALLOCATE(ZY0(NIU,NJU,NKU))
233 ALLOCATE(ZZ0(NIU,NJU,NKU))
234 ALLOCATE(ZWORK(NIU,NJU,NKU,inbr_field+3))
235 ALLOCATE(YUNITE(inbr_field))
236 ALLOCATE(ZX00(NIU,NJU,NKU))
237 ALLOCATE(ZY00(NIU,NJU,NKU))
238 ALLOCATE(ZZ00(NIU,NJU,NKU))
239 ALLOCATE(ZZL(NIU,NJU,NKU))
240 !
241 end if
242 ! initial values
243 ZXOR=0.5 * (XXHAT(2)+XXHAT(3)) 
244 ZYOR=0.5 * (XYHAT(2)+XYHAT(3))
245 ZDX= XXHAT(3)-XXHAT(2)
246 ZDY= XYHAT(3)-XYHAT(2)
247 !ZZL=MZF(XZZ)
248 do k=1,nku-1
249   zzl(:,:,k)=(XZZ(:,:,k)+XZZ(:,:,k+1))*0.5
250 end do
251 ZZL(:,:,NKU)=2*XZZ(:,:,NKU)-ZZL(:,:,NKU-1)
252 ZXMAX=ZXOR+(NIU-3)*ZDX
253 ZYMAX=ZYOR+(NJU-3)*ZDY
254 ZZMAX=ZZL(2,2,NKU-1)
255 !  conversion from m to km
256 ZXOR=ZXOR*1.E-3
257 ZYOR=ZYOR*1.E-3
258 ZDX=ZDX*1.E-3
259 ZDY=ZDY*1.E-3
260 ZZL(:,:,:)=ZZL(:,:,:)*1.E-3
261 ZXMAX=ZXMAX*1.E-3
262 ZYMAX=ZYMAX*1.E-3
263 ZZMAX=ZZMAX*1.E-3
264 !
265 CALL READVAR('LGXM',CFILES(ifilecur),CFLAGFILE,nverbdia,iret)
266 ZX00(:,:,:)=XVAR(:,:,:,1,1,1)   
267 CALL READVAR('LGYM',CFILES(ifilecur),CFLAGFILE,nverbdia,iret)
268 ZY00(:,:,:)=XVAR(:,:,:,1,1,1)  
269 CALL READVAR('LGZM',CFILES(ifilecur),CFLAGFILE,nverbdia,iret)
270 ZZ00(:,:,:)=XVAR(:,:,:,1,1,1) 
271 ! what is the unit of Lag. var. (km after DIAG, m after MODEL) ?
272 IF (INDEX(CCOMMENT(1),'KM')/=0 .OR. &
273     MAXVAL(ZZ00(:,:,:))<100.         ) THEN
274   print*,'unit of Lagrangian variables in ',TRIM(CFILES(ifilecur)),' is KM' 
275 ELSE
276   print*,'unit of Lagrangian variables in ',TRIM(CFILES(ifilecur)),' is M' 
277   ZX00(:,:,:)=ZX00(:,:,:)*1.E-3  !  conversion from m to km
278   ZY00(:,:,:)=ZY00(:,:,:)*1.E-3
279   ZZ00(:,:,:)=ZZ00(:,:,:)*1.E-3
280 ENDIF
281 !
282 !
283 IF (L2D) THEN
284   WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
285           ZZ00>ZZMAX)
286     ZX00=ZSPVAL
287     ZZ00=ZSPVAL
288   END WHERE
289 ELSE
290   WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
291           ZY00<ZYOR .OR. ZY00>ZYMAX .OR. &
292               ZZ00>ZZMAX)
293     ZX00=ZSPVAL
294     ZY00=ZSPVAL
295     ZZ00=ZSPVAL
296   END WHERE
297 END IF
298 !
299 !-------------------------------------------------------------------------------
300 !
301 !*       4.0    COMPUTE THE ORIGIN STEP BY STEP
302 !               -------------------------------
303 !
304 !
305 ! General loop for the files where a reinitialisation of the tracers 
306 ! is performed
307 DO JFILECUR=1,NFILES
308   !
309   !CALL FMOPEN_ll(CFILES(NBRFILES(JFILECUR)),'READ',CLUOUT,   &
310   !               INPRAR,ITYPE,NVERB,ININAR,IRESP)
311 !
312 !*       4.1  check if this file is a start instant
313 !
314   GSTART=.FALSE.
315   DO JLOOP=1,NFILES
316     IF (NBRFILES(JFILECUR)==NSTART_SUPP(JLOOP) .OR. JFILECUR==NFILES) THEN
317       INBR_START=INBR_START-1
318       GSTART=.TRUE.
319       EXIT
320     END IF
321   ENDDO
322   !
323   if (nverbdia>0) then
324    print *, 'fichier pour la reconstitution ',JFILECUR,' GSTART =',GSTART
325   end if
326 !
327 !*       4.2 read the potential temp or the water vapor at the start instant      
328 !
329   IF (GSTART) THEN
330     !
331     if(inbr_field>0) then
332      do ifield=1,inbr_field
333       YRECFM=CFIELD_LAG(ifield)
334       CALL READVAR(YRECFM,CFILES(NBRFILES(JFILECUR)),CFLAGFILE &
335                     ,nverbdia,iret)
336       ZWORK(:,:,:,ifield)=XVAR(:,:,:,1,1,1)
337       YUNITE(ifield)=CUNITE(1)
338      end do
339     else
340       CALL READVAR('PABSM',CFILES(NBRFILES(JFILECUR)),CFLAGFILE &
341                     ,nverbdia,iret)
342     endif
343     TDTCUR_START%TDATE%YEAR=XDATIME(5,1)
344     TDTCUR_START%TDATE%MONTH=XDATIME(6,1)
345     TDTCUR_START%TDATE%DAY=XDATIME(7,1)
346     TDTCUR_START%TIME=XDATIME(8,1)
347     IHOUR   = INT(TDTCUR_START%TIME/3600.)
348     ZREMAIN = MOD(TDTCUR_START%TIME,3600.)
349     IMINUTE = INT(ZREMAIN/60.)
350     ZSECOND = MOD(ZREMAIN,60.)
351     WRITE(YDATE,FMT='(1X,I4.4,I2.2,I2.2,2X,I2.2,"H",I2.2,"M", &
352          & F5.2,"S")') TDTCUR_START%TDATE, IHOUR,IMINUTE,ZSECOND 
353   END IF
354 !
355 !*       4.3  store the X0,Y0,Z0 field for the current start before 
356 !             computing the new origin
357 !
358   IF (GSTART) THEN
359     IGRID=1
360     PRINT *,'INBR_START',INBR_START,' NBRFILES(JFILECUR)',NBRFILES(JFILECUR)
361     WRITE(YRECFM,'(A2,I2.2)')'X0',INBR_START
362     WRITE(YCOMMENT,'(A8,I2.2)')'X_Y_Z_X0',INBR_START
363     CTITRE(1)=YRECFM
364     CUNITE(1)='(KM)' 
365     CCOMMENT(1)=YCOMMENT(1:10)//YDATE//' (KM)'
366     PRINT *,'COMMENT = ',CCOMMENT(1)
367     XVAR(:,:,:,1,1,1)=ZX00(:,:,:)
368     CALL WRITEVAR(1,NIU,1,NJU,1,NKU,1,1,1,1,1,1, &
369          YRECFM,HFMFILE,'OLD',CNAME_SUP,nverbdia,iret)
370     !
371     WRITE(YRECFM,'(A2,I2.2)')'Y0',INBR_START
372     WRITE(YCOMMENT,'(A8,I2.2)')'X_Y_Z_Y0',INBR_START
373     CTITRE(1)=YRECFM
374     CCOMMENT(1)=YCOMMENT(1:10)//YDATE//' (KM)'
375     CUNITE(1)='(KM)' 
376     PRINT *,'COMMENT = ',CCOMMENT(1)
377     XVAR(:,:,:,1,1,1)=ZY00(:,:,:)
378     CALL WRITEVAR(1,NIU,1,NJU,1,NKU,1,1,1,1,1,1, &
379          YRECFM,HFMFILE,'OLD',CNAME_SUP,nverbdia,iret)
380     !
381     WRITE(YRECFM,'(A2,I2.2)')'Z0',INBR_START
382     WRITE(YCOMMENT,'(A8,I2.2)')'X_Y_Z_Z0',INBR_START
383     CTITRE(1)=YRECFM
384     CCOMMENT(1)=YCOMMENT(1:10)//YDATE//' (KM)'
385     CUNITE(1)='(KM)' 
386     PRINT *,'COMMENT = ',CCOMMENT(1)
387     XVAR(:,:,:,1,1,1)=ZZ00(:,:,:)
388     CALL WRITEVAR(1,NIU,1,NJU,1,NKU,1,1,1,1,1,1, &
389          YRECFM,HFMFILE,'OLD',CNAME_SUP,nverbdia,iret)
390   END IF
391 !
392 !*       4.4   compute the origin of the particules using one more segment
393 !
394   IF (JFILECUR /= NFILES) THEN
395     CALL READVAR('LGXM',CFILES(NBRFILES(JFILECUR)),  &
396                   CFLAGFILE,nverbdia,iret)
397     ZX0(:,:,:)=XVAR(:,:,:,1,1,1) 
398     CALL READVAR('LGYM',CFILES(NBRFILES(JFILECUR)),  &
399                   CFLAGFILE,nverbdia,iret)
400     ZY0(:,:,:)=XVAR(:,:,:,1,1,1) 
401     CALL READVAR('LGZM',CFILES(NBRFILES(JFILECUR)),  &
402                   CFLAGFILE,nverbdia,iret)
403     ZZ0(:,:,:)=XVAR(:,:,:,1,1,1)  
404     ! what is the unit of Lag. var. (km after DIAG, m after MODEL) ?
405     IF (INDEX(CCOMMENT(1),'KM')/=0 .OR. &
406        MAXVAL(ZZ00(:,:,:))<100.         ) THEN
407       print*,'unit of Lagrangian variables in ', &
408              TRIM(CFILES(NBRFILES(jfilecur))),' is KM' 
409     ELSE
410       print*,'unit of Lagrangian variables in ', &
411              TRIM(CFILES(NBRFILES(jfilecur))),' is M' 
412       ZX00(:,:,:)=ZX00(:,:,:)*1.E-3  !  conversion from m to km
413       ZY00(:,:,:)=ZY00(:,:,:)*1.E-3
414       ZZ00(:,:,:)=ZZ00(:,:,:)*1.E-3
415     ENDIF
416     !
417     ! old position of the set of particles
418     ZWORK(:,:,:,inbr_field+1)=ZX00
419     ZWORK(:,:,:,inbr_field+2)=ZY00
420     ZWORK(:,:,:,inbr_field+3)=ZZ00
421     !
422     IF (L2D) THEN
423       CALL INTERPXYZ(ZWORK(:,:,:,inbr_field+1),ZWORK(:,:,:,inbr_field+2),&
424                      ZWORK(:,:,:,inbr_field+3),ZX0,ZX00,ZZ0,ZZ00             )
425     ELSE
426       CALL INTERPXYZ(ZWORK(:,:,:,inbr_field+1),ZWORK(:,:,:,inbr_field+2),&
427                      ZWORK(:,:,:,inbr_field+3),ZX0,ZX00,ZY0,ZY00,ZZ0,ZZ00    )
428     END IF
429     !
430     IF (L2D) THEN
431       WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
432               ZZ00>ZZMAX)
433         ZX00=ZSPVAL
434         ZZ00=ZSPVAL
435       END WHERE
436     ELSE
437       WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
438               ZY00<ZYOR .OR. ZY00>ZYMAX .OR. &
439               ZZ00>ZZMAX)
440         ZX00=ZSPVAL
441         ZY00=ZSPVAL
442         ZZ00=ZSPVAL
443       END WHERE
444     END IF
445     !
446     !
447   END IF
448 !
449 !*       4.5   close the input file
450 !
451   !!CALL FMCLOS_ll(CFILES(NBRFILES(JFILECUR)),'KEEP',CLUOUT,IRESP)
452 !
453 !
454 !*       4.6   compute and store potential temp and water vapor at the origin
455 !
456   IF (GSTART) THEN
457     !
458     do ifield=1,inbr_field
459     !
460       CALL INTERPXYZ(ZX00,ZY00,ZZ00,     &
461       ZWORK(:,:,:,ifield),ZWORK(:,:,:,inbr_field+1)         )
462     !
463     WRITE(YRECFM,'(A3,I2.2)')CFIELD_LAG(ifield),INBR_START
464     CTITRE(1)=YRECFM
465     print*,'CFIELD_LAG ',ifield,' TITRE= ',TRIM(CTITRE(1))
466     WRITE(YCOMMENT,'(A6,A3,I2.2)')'X_Y_Z_',CFIELD_LAG(ifield),INBR_START
467     CCOMMENT(1)=YCOMMENT(1:10)//YDATE//' (USI)'
468     PRINT *,'COMMENT = ',TRIM(CCOMMENT(1))
469     CUNITE(1)=YUNITE(ifield)
470     PRINT *,'CUNIT = ',TRIM(CUNITE(1))
471     XVAR(:,:,:,1,1,1)=ZWORK(:,:,:,ifield)
472     CALL WRITEVAR(1,NIU,1,NJU,1,NKU,1,1,1,1,1,1, &
473          YRECFM,HFMFILE,'OLD',CNAME_SUP,nverbdia,iret)
474     !
475     !
476     end do
477     !
478   END IF
479 !
480 !
481 END DO
482 !
483 ! fermeture du fichier diachronique
484 IF (GSTART) call WRITEVAR(1,NIU,1,NJU,1,NKU,1,1,1,1,1,1, &
485                  YRECFM,HFMFILE,'CLO',CNAME_SUP,nverbdia,iret)
486 end do
487 !***********************************************
488 !***********************************************
489 !
490 PRINT*, ' '
491 PRINT*, 'COMPUTE_R00 AFTER ORIGIN COMPUTATIONS AND STORAGE'
492 !
493 !-------------------------------------------------------------------------------
494 !!
495 CONTAINS
496 !
497 !
498 !-------------------------------------------------------------------------------
499 !
500 !
501 SUBROUTINE INTERPXYZ(PX,PY,PZ,PIN1,POUT1,PIN2,POUT2,PIN3,POUT3)
502 !
503 !
504 !*      0. DECLARATIONS
505 !          ------------
506 !
507 !*       0.1  declaration of arguments
508 !
509 REAL, INTENT(IN),  DIMENSION(:,:,:)           :: PX,PY,PZ
510 REAL, INTENT(IN),  DIMENSION(:,:,:)           :: PIN1
511 REAL, INTENT(OUT), DIMENSION(:,:,:)           :: POUT1
512 REAL, INTENT(IN),  DIMENSION(:,:,:), OPTIONAL :: PIN2,PIN3
513 REAL, INTENT(OUT), DIMENSION(:,:,:), OPTIONAL :: POUT2,POUT3   
514 !
515 !*       0.2  declaration of local variables
516 !
517 INTEGER  :: JI,JJ,JK,JKK    ! loop index
518 INTEGER  :: II,IJ,IK        ! grid index for the interpolation
519 REAL     :: ZXREL,ZYREL     ! fractional grid index for the interpolation
520 REAL, DIMENSION(SIZE(PIN1,3)) :: ZZLXY ! vertical grid at the interpolated point
521 REAL     :: ZEPS1,ZEPS2,ZEPS3          ! coeff. for the interpolation
522 REAL     :: ZX,ZY,ZZ
523 LOGICAL  :: GEXT
524 !
525 !-------------------------------------------------------------------------------
526 !
527 DO JK=1,NKU
528   DO JJ=1,NJU
529     DO JI=1,NIU
530       !
531       ZX=PX(JI,JJ,JK) 
532       ZY=PY(JI,JJ,JK)
533       ZZ=PZ(JI,JJ,JK)
534       !
535       ! remove external points
536       IF (L2D) THEN
537         GEXT=(ZX==ZSPVAL).OR.(ZZ==ZSPVAL)
538       ELSE
539         GEXT=(ZX==ZSPVAL).OR.(ZY==ZSPVAL).OR.(ZZ==ZSPVAL)
540       END IF
541       IF (GEXT) THEN
542         POUT1(JI,JJ,JK) = ZSPVAL
543         IF (PRESENT(PIN2)) THEN
544           POUT2(JI,JJ,JK) = ZSPVAL
545         END IF
546         IF (PRESENT(PIN3)) THEN
547           POUT3(JI,JJ,JK) = ZSPVAL
548         ENDIF
549         !
550         CYCLE
551         !
552       END IF
553       !
554       ZXREL=(ZX-ZXOR)/ZDX+2
555       ZYREL=(ZY-ZYOR)/ZDY+2
556       !
557       II=FLOOR(ZXREL)
558       IJ=FLOOR(ZYREL)
559       !
560       ZEPS1=ZXREL-REAL(II)
561       ZEPS2=ZYREL-REAL(IJ)
562       IF (L2D) ZEPS2=0.
563       !
564       DO JKK=1,NKU
565         ZZLXY(JKK)=ZEPS2*(ZEPS1*(ZZL(II+1,IJ+1,JKK))+(1-ZEPS1)*(ZZL(II,IJ+1,JKK)))     &
566              + (1-ZEPS2)*(ZEPS1*(ZZL(II+1,IJ,JKK))+(1-ZEPS1)*(ZZL(II,IJ,JKK)))
567       ENDDO
568       !
569       IK=999
570       DO JKK=2,NKU
571         IF (ZZLXY(JKK).GE.ZZ) THEN
572           IK=JKK-1
573           EXIT 
574         ENDIF
575       ENDDO
576       !
577       IF (IK==999) THEN
578         PRINT*,'PROBLEM AT POINT',II,IJ
579         PRINT*,'XREL, YREL, Z =',ZXREL,ZYREL,ZZ
580         PRINT*,'ZZLXY(NKU)',ZZLXY(NKU)
581         STOP
582       END IF 
583       !
584       ZEPS3=(ZZ-ZZLXY(IK))/(ZZLXY(IK+1)-ZZLXY(IK))
585       !
586       POUT1(JI,JJ,JK) =                                                       & 
587         ZEPS3 *                                                               &
588       (  ZEPS2*(ZEPS1*(PIN1(II+1,IJ+1,IK+1))+(1-ZEPS1)*(PIN1(II,IJ+1,IK+1)))  &
589        + (1-ZEPS2)*(ZEPS1*(PIN1(II+1,IJ,IK+1))+(1-ZEPS1)*(PIN1(II,IJ,IK+1)))  &
590       )                                                                       & 
591       + (1-ZEPS3) *                                                           &
592       (  ZEPS2*(ZEPS1*(PIN1(II+1,IJ+1,IK))+(1-ZEPS1)*(PIN1(II,IJ+1,IK)))      &
593        + (1-ZEPS2)*(ZEPS1*(PIN1(II+1,IJ,IK))+(1-ZEPS1)*(PIN1(II,IJ,IK)))      &
594       )
595       IF (PRESENT(POUT2)) THEN
596         POUT2(JI,JJ,JK) =                                                     & 
597           ZEPS3 *                                                             &
598         (  ZEPS2*(ZEPS1*(PIN2(II+1,IJ+1,IK+1))+(1-ZEPS1)*(PIN2(II,IJ+1,IK+1)))&
599          + (1-ZEPS2)*(ZEPS1*(PIN2(II+1,IJ,IK+1))+(1-ZEPS1)*(PIN2(II,IJ,IK+1)))&
600         )                                                                     & 
601         + (1-ZEPS3) *                                                         &
602         (  ZEPS2*(ZEPS1*(PIN2(II+1,IJ+1,IK))+(1-ZEPS1)*(PIN2(II,IJ+1,IK)))    &
603          + (1-ZEPS2)*(ZEPS1*(PIN2(II+1,IJ,IK))+(1-ZEPS1)*(PIN2(II,IJ,IK)))    &
604         )
605       ENDIF
606         !
607       IF (PRESENT(POUT3)) THEN
608         POUT3(JI,JJ,JK) =                                                     & 
609           ZEPS3 *                                                             &
610         (  ZEPS2*(ZEPS1*(PIN3(II+1,IJ+1,IK+1))+(1-ZEPS1)*(PIN3(II,IJ+1,IK+1)))&
611          + (1-ZEPS2)*(ZEPS1*(PIN3(II+1,IJ,IK+1))+(1-ZEPS1)*(PIN3(II,IJ,IK+1)))&
612         )                                                                     &
613         + (1-ZEPS3) *                                                         &
614         (  ZEPS2*(ZEPS1*(PIN3(II+1,IJ+1,IK))+(1-ZEPS1)*(PIN3(II,IJ+1,IK)))    &
615          + (1-ZEPS2)*(ZEPS1*(PIN3(II+1,IJ,IK))+(1-ZEPS1)*(PIN3(II,IJ,IK)))    &
616         )
617       ENDIF
618       !
619     END DO
620   END DO
621 END DO
622 !
623 !-------------------------------------------------------------------------------
624 !
625 !
626 END SUBROUTINE INTERPXYZ
627 !
628 !-------------------------------------------------------------------------------
629 !
630 END program