Newer
Older
!SFX_LIC Copyright 2004-2018 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
SUBROUTINE INTERPOL_NPTS (UG, U, HPROGRAM,KLUOUT,KNPTS,KCODE,PX,PY,PFIELD,KNEAR_NBR)
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
! #########################################################
!
!!**** *INTERPOL_NPTS* interpolates with ###ine f77 programs a 2D field
!! from all grid points valid values
!!
!! PURPOSE
!! -------
!!
!! The points are all on only one grid (defined with the coordinates
!! of all the points). The code to apply for each point is:
!!
!! KCODE>0 : data point (with field valid for interpolation)
!! KCODE=-1: point to ignore
!! KCODE=0 : point to interpolate
!!
!!
!!
!! METHOD
!! ------
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!!
!!
!! REFERENCE
!! ---------
!!
!! AUTHOR
!! ------
!!
!! V. Masson Meteo-France
!!
!! MODIFICATION
!! ------------
!!
!! Original 03/2004
!! Modification
!! B. Decharme 2014 scan all point case if gaussien grid or NHALO = 0
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!----------------------------------------------------------------------------
!
!* 0. DECLARATION
! -----------
!
!
USE MODD_SURF_ATM_GRID_n, ONLY : SURF_ATM_GRID_t
USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
!
USE MODD_SURFEX_OMP, ONLY : NBLOCK
USE MODD_SURFEX_MPI, ONLY : NRANK, NPIO, NCOMM, NPROC, IDX_I, NINDEX,NNUM
USE MODD_SURF_PAR, ONLY : XUNDEF
!
USE MODI_GATHER_AND_WRITE_MPI
USE MODI_READ_AND_SEND_MPI
!
USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
USE PARKIND1 ,ONLY : JPRB
!
USE MODD_IO_ll, ONLY : ISIOP, ISP, ISNPROC
USE MODD_VAR_ll, ONLY : NMNH_COMM_WORLD
USE MODE_GATHER_ll
USE MODE_TOOLS_ll, ONLY : GET_GLOBALDIMS_ll
USE MODD_PARAMETERS_ll,ONLY : JPHEXT
!
USE MODD_IO_SURF_MNH, ONLY : NIU, NJU
#endif
!
#if defined(SFX_MPI) || defined(SFX_MNH)
INCLUDE "mpif.h"
#endif
!
!* 0.1 Declaration of arguments
! ------------------------
!
!
TYPE(SURF_ATM_GRID_t), INTENT(INOUT) :: UG
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! host program
INTEGER, INTENT(IN) :: KLUOUT ! output listing
INTEGER, INTENT(IN) :: KNPTS ! number of points to interpolate with
INTEGER,DIMENSION(:), INTENT(INOUT) :: KCODE ! code for each point
! >0 point used for interpolation
! 0 point to interpolate
! -1 point not used
! -2 point not used
! ! -3 if spline is no computed
! ! for this point
REAL, DIMENSION(:), INTENT(IN) :: PX ! x of each grid mesh.
REAL, DIMENSION(:), INTENT(IN) :: PY ! y of each grid mesh.
REAL, DIMENSION(:,:),INTENT(INOUT) :: PFIELD ! pgd field on grid mesh.
!
!* 0.2 Declaration of local variables
! ------------------------------
!
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZFIELD, ZFIELD2
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFIELD3
REAL, DIMENSION(:,:), ALLOCATABLE :: ZNDIST ! 3 nearest square distances
REAL, DIMENSION(:,:), ALLOCATABLE :: ZNVAL ! 3 corresponding field values
REAL, DIMENSION(:), ALLOCATABLE :: ZX, ZY
REAL, DIMENSION(:), ALLOCATABLE :: ZDIST ! square distance between two interpolating and interpolated points
!
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: ININD0, ININD_ALL ! 3 corresponding field values
INTEGER, DIMENSION(:,:), ALLOCATABLE :: ININD
INTEGER, DIMENSION(:), ALLOCATABLE :: IINDEX ! list of index to scan
INTEGER, DIMENSION(:), ALLOCATABLE :: ISIZE, ISIZE_TOT
!
REAL :: ZSUM
INTEGER :: IP, ICPT, IL1, JL, JP, JK, JKK
!
INTEGER :: ICOUNT, IL2 ! counter
INTEGER :: INPTS
INTEGER :: INFOMPI, IDIM_FULL, INEAR_NBR, IOLD
!
#if defined(SFX_MPI) || defined(SFX_MNH)
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
#endif
!
#ifdef SFX_MNH
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOORD_2D, ZCOORD_2D_ALL
INTEGER, DIMENSION(:,:), ALLOCATABLE :: ISIZE_2D, ISIZE_2D_ALL, INUM_2D, INUM_2D_ALL, IINDEX_2D, IINDEX_2D_ALL
INTEGER, DIMENSION(:), ALLOCATABLE :: INUM_1D, IINDEX_1D, INUM_TOT, IINDEX_TOT
INTEGER :: IIMAX, IJMAX, IIU, IJU, JI
INTEGER :: IRANK_SAVE, IPROC_SAVE, IPIO_SAVE, ICOMM_SAVE
#endif
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_1',0,ZHOOK_HANDLE)
IOLD = 0
!NHALO /= 0
IF ( KNEAR_NBR/=U%NDIM_FULL ) THEN
IOLD = 1
ELSE
! case mesonh and NHALO = 0
#ifdef SFX_MNH
IOLD = 2
! case NHALO = 0
#else
IOLD = 1
#endif
ENDIF
IF (IOLD==2) THEN
!
#ifdef SFX_MNH
!physical dimensions of the task
CALL GET_DIM_PHYS_ll('B',IIU,IJU)
! total dimensions
CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
IDIM_FULL = (IIMAX) * (IJMAX)
INEAR_NBR = IDIM_FULL
!
! on sauve les infos de bases
IRANK_SAVE = NRANK
IPROC_SAVE = NPROC
IPIO_SAVE = NPIO
ICOMM_SAVE = NCOMM
!
! on met les infos de mésonh
NRANK = ISP-1
NPROC = ISNPROC
NPIO = ISIOP-1
NCOMM = NMNH_COMM_WORLD
!
!KCODE to 2D
ALLOCATE(ISIZE_2D(IIU+2*JPHEXT,IJU+2*JPHEXT))
ISIZE_2D(:,:) = -1
ISIZE_2D(1+JPHEXT:IIU+JPHEXT,1+JPHEXT:IJU+JPHEXT) = RESHAPE(KCODE, (/ IIU,IJU /) )
! tasks to whole domaine
ALLOCATE(ISIZE_2D_ALL(IIMAX+2*JPHEXT,IJMAX+2*JPHEXT))
CALL GATHER_XYFIELD(ISIZE_2D,ISIZE_2D_ALL,ISIOP,NMNH_COMM_WORLD)
DEALLOCATE(ISIZE_2D)
ALLOCATE(ISIZE_TOT(IDIM_FULL))
ISIZE_TOT = RESHAPE(ISIZE_2D_ALL(1+JPHEXT:IIMAX+JPHEXT,1+JPHEXT:IJMAX+JPHEXT),(/IDIM_FULL/))
DEALLOCATE(ISIZE_2D_ALL)
#endif
!
ELSEIF (IOLD==1) THEN
!
IDIM_FULL = U%NDIM_FULL
INEAR_NBR = KNEAR_NBR
!
ALLOCATE(ISIZE_TOT(IDIM_FULL))
CALL GATHER_AND_WRITE_MPI(KCODE,ISIZE_TOT)
!
ENDIF
!...known by all tasks
IF (NPROC>1) THEN
#if defined(SFX_MPI) || defined(SFX_MNH)
CALL MPI_BCAST(ISIZE_TOT,IDIM_FULL*KIND(ISIZE_TOT)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
#endif
ENDIF
IF (ALL(ISIZE_TOT/=0)) THEN
IF (IOLD==2) THEN
#ifdef SFX_MNH
NRANK = IRANK_SAVE
NPROC = IPROC_SAVE
NPIO = IPIO_SAVE
NCOMM = ICOMM_SAVE
#endif
ENDIF
DEALLOCATE(ISIZE_TOT)
CALL DR_HOOK('INTERPOL_NPTS_1',1,ZHOOK_HANDLE)
RETURN
ENDIF
IL1 = SIZE(PFIELD,1)
IL2 = SIZE(PFIELD,2)
IF (IOLD==2) THEN
!
#ifdef SFX_MNH
! NNUM
ALLOCATE(INUM_1D(U%NDIM_FULL))
DO JI = 1,U%NDIM_FULL
INUM_1D(JI) = JI
ENDDO
! INUM_1D to 2D
ALLOCATE(INUM_2D(IIU+2*JPHEXT,IJU+2*JPHEXT))
INUM_2D = 0
INUM_2D(1+JPHEXT:IIU+JPHEXT,1+JPHEXT:IJU+JPHEXT) = RESHAPE(INUM_1D, (/ IIU,IJU /) )
DEALLOCATE(INUM_1D)
! tasks to whole domaine
ALLOCATE(INUM_2D_ALL(IIMAX+2*JPHEXT,IJMAX+2*JPHEXT))
CALL GATHER_XYFIELD(INUM_2D,INUM_2D_ALL,ISIOP,NMNH_COMM_WORLD)
DEALLOCATE(INUM_2D)
ALLOCATE(INUM_TOT(IDIM_FULL))
INUM_TOT = RESHAPE(INUM_2D_ALL(1+JPHEXT:IIMAX+JPHEXT,1+JPHEXT:IJMAX+JPHEXT),(/IDIM_FULL/))
DEALLOCATE(INUM_2D_ALL)
!
! NINDEX
ALLOCATE(IINDEX_1D(U%NDIM_FULL))
IINDEX_1D(:) = ISP - 1
! INDEX_1D to 2D
ALLOCATE(IINDEX_2D(IIU+2*JPHEXT,IJU+2*JPHEXT))
IINDEX_2D(:,:) = 0
IINDEX_2D(1+JPHEXT:IIU+JPHEXT,1+JPHEXT:IJU+JPHEXT) = RESHAPE(IINDEX_1D, (/ IIU,IJU /) )
DEALLOCATE(IINDEX_1D)
! tasks to whole domaine
ALLOCATE(IINDEX_2D_ALL(IIMAX+2*JPHEXT,IJMAX+2*JPHEXT))
CALL GATHER_XYFIELD(IINDEX_2D,IINDEX_2D_ALL,ISIOP,NMNH_COMM_WORLD)
DEALLOCATE(IINDEX_2D)
ALLOCATE(IINDEX_TOT(IDIM_FULL))
IINDEX_TOT = RESHAPE(IINDEX_2D_ALL(1+JPHEXT:IIMAX+JPHEXT,1+JPHEXT:IJMAX+JPHEXT),(/IDIM_FULL/))
DEALLOCATE(IINDEX_2D_ALL)
!
!PX, PY: coordinates of all the points
!ZX, ZY: coordinates of the points in this task
ALLOCATE(ZCOORD_2D(IIU+2*JPHEXT,IJU+2*JPHEXT))
ALLOCATE(ZCOORD_2D_ALL(IIMAX+2*JPHEXT,IJMAX+2*JPHEXT))
!
ZCOORD_2D = 0.
ZCOORD_2D(1+JPHEXT:IIU+JPHEXT,1+JPHEXT:IJU+JPHEXT) = RESHAPE(PX, (/ IIU,IJU /) )
CALL GATHER_XYFIELD(ZCOORD_2D,ZCOORD_2D_ALL,ISIOP,NMNH_COMM_WORLD)
ALLOCATE(ZX(IDIM_FULL))
ZX = RESHAPE(ZCOORD_2D_ALL(1+JPHEXT:IIMAX+JPHEXT,1+JPHEXT:IJMAX+JPHEXT),(/IDIM_FULL/))
!
ZCOORD_2D = 0.
ZCOORD_2D(1+JPHEXT:IIU+JPHEXT,1+JPHEXT:IJU+JPHEXT) = RESHAPE(PY, (/ IIU,IJU /) )
CALL GATHER_XYFIELD(ZCOORD_2D,ZCOORD_2D_ALL,ISIOP,NMNH_COMM_WORLD)
ALLOCATE(ZY(IDIM_FULL))
ZY = RESHAPE(ZCOORD_2D_ALL(1+JPHEXT:IIMAX+JPHEXT,1+JPHEXT:IJMAX+JPHEXT),(/IDIM_FULL/))
!
DEALLOCATE(ZCOORD_2D,ZCOORD_2D_ALL)
!
IF (NPROC>1) THEN
CALL MPI_BCAST(INUM_TOT,IDIM_FULL*KIND(INUM_TOT)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
CALL MPI_BCAST(IINDEX_TOT,IDIM_FULL*KIND(IINDEX_TOT)/4,MPI_INTEGER,NPIO,NCOMM,INFOMPI)
WAUTELET Philippe
committed
CALL MPI_BCAST(ZX,IDIM_FULL*KIND(ZX)/4,MPI_REAL,NPIO,NCOMM,INFOMPI)
CALL MPI_BCAST(ZY,IDIM_FULL*KIND(ZY)/4,MPI_REAL,NPIO,NCOMM,INFOMPI)
#endif
!
ELSEIF (IOLD==1) THEN
!
!PX, PY: coordinates of all the points
!ZX, ZY: coordinates of the points in this task
ALLOCATE(ZX(IL1),ZY(IL1))
CALL READ_AND_SEND_MPI(PX,ZX)
CALL READ_AND_SEND_MPI(PY,ZY)
!
IF (INEAR_NBR/=IDIM_FULL) THEN
IF (.NOT.ASSOCIATED(UG%NNEAR)) THEN
ALLOCATE(UG%NNEAR(IL1,INEAR_NBR))
!seach near meshes in the complete grid (xgrid_full_par) for this task
CALL GET_NEAR_MESHES(UG%G%CGRID,UG%NGRID_FULL_PAR,IDIM_FULL,UG%XGRID_FULL_PAR,INEAR_NBR,UG%NNEAR)
ENDIF
ENDIF
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_1',1,ZHOOK_HANDLE)
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_2',0,ZHOOK_HANDLE)
!
ALLOCATE(IINDEX(INEAR_NBR))
IINDEX (:) = 0
!
IF (INEAR_NBR==IDIM_FULL) THEN
ICOUNT = 0
DO JL=1,IDIM_FULL
!is the neareast point available to interpolation
IF (ISIZE_TOT(JL)>0) THEN
ICOUNT = ICOUNT+1
IINDEX(ICOUNT) = JL
END DO
!did we found enough points for interpolate
IF (ICOUNT>=1) THEN
INPTS = MIN(KNPTS,ICOUNT)
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
WHERE(KCODE(:)==0) KCODE(:) = -4
END IF
ALLOCATE(ZDIST (ICOUNT))
ELSE
ALLOCATE(ZDIST (INEAR_NBR))
ENDIF
!
ZDIST (:) = 0.
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_2',1,ZHOOK_HANDLE)
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_3',0,ZHOOK_HANDLE)
!
!indexes of points used for interpolation, for each point to interpolate
ALLOCATE(ININD (IP,KNPTS))
ININD(:,:) = 0
!distances of the points used for interpolation
ALLOCATE(ZNDIST(IP,0:KNPTS))
ZNDIST (:,1:KNPTS) = 1.E20
ZNDIST (:,0) = 0.
!
ICPT=0
!loop on points for this task
DO JL=1,IL1
!
!does this point need to be interpolated?
IF (KCODE(JL)/=0) CYCLE
!
IF (INEAR_NBR/=IDIM_FULL) THEN
DO JK=1,INEAR_NBR
IF (UG%NNEAR(JL,JK)>0) THEN
!is the neareast point available to interpolation
IF (ISIZE_TOT(UG%NNEAR(JL,JK))>0) THEN
!did we found enough points for interpolate
IF (ICOUNT>=KNPTS) THEN
INPTS = KNPTS
ELSE
KCODE(JL) = -4
CYCLE
END IF
!
!one point more to interpolate
ICPT = ICPT + 1
!
IF (IOLD==2) THEN
#ifdef SFX_MNH
ZDIST(1:ICOUNT) = (ZX(IINDEX(1:ICOUNT))-PX(JL))**2 + (ZY(IINDEX(1:ICOUNT))-PY(JL))**2
#endif
ZDIST(1:ICOUNT) = (PX(IINDEX(1:ICOUNT))-ZX(JL))**2 + (PY(IINDEX(1:ICOUNT))-ZY(JL))**2
ENDIF
IF (ZDIST(JP)>ZNDIST(ICPT,INPTS)) CYCLE
IF ( ZDIST(JP)>ZNDIST(ICPT,JK-1) ) THEN
IF ( JK<INPTS ) THEN
DO JKK = INPTS,JK+1,-1
ZNDIST(ICPT,JKK) = ZNDIST(ICPT,JKK-1)
ININD (ICPT,JKK) = ININD (ICPT,JKK-1)
!distances and indexes of points used to interpolate are saved
ZNDIST(ICPT,JK) = ZDIST (JP)
ININD (ICPT,JK) = IINDEX(JP)
!
EXIT
!
ENDIF
!
ENDDO
!
ENDDO
!
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
ENDDO
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_3',1,ZHOOK_HANDLE)
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_4',0,ZHOOK_HANDLE)
!
DEALLOCATE(IINDEX,ZX,ZY,ISIZE_TOT,ZDIST)
!
ZNDIST(:,:) = SQRT(ZNDIST(:,:))
!
ALLOCATE(ISIZE(0:NPROC-1))
!numbers of points to interpolated are gathered
IF (NPROC>1) THEN
#if defined(SFX_MPI) || defined(SFX_MNH)
CALL MPI_ALLGATHER(ICPT,KIND(ICPT)/4,MPI_INTEGER,&
ISIZE,KIND(ISIZE)/4,MPI_INTEGER,NCOMM,INFOMPI)
#endif
ELSE
ISIZE(:) = ICPT
ENDIF
!
!this array contains, for each point to interpolate,
!the correspondance between the task where is the point to use for interpolation
!(NINDEX) and its index in this task (NNUM)
ALLOCATE(ININD0(MAXVAL(ISIZE),KNPTS,0:NPROC-1))
ININD0(:,:,:) = 0
!
!number of points needed to interpolated
DO JL=1,KNPTS
!number of points to interpolated
DO JP=1,ICPT
!index of the point needed in the whole grid
JK = ININD(JP,JL)
!inind0 contains the task and the index in this task for this point
IF (IOLD==2) THEN
#ifdef SFX_MNH
IF (JK/=0) ININD0(JP,JL,IINDEX_TOT(JK)) = INUM_TOT(JK)
#endif
ELSEIF (IOLD==1) THEN
IF (JK/=0) ININD0(JP,JL,NINDEX(JK)) = NNUM(JK)
ENDIF
ENDDO
ENDDO
!
!for each task, points needed and located in this task are gathered
ALLOCATE(ININD_ALL(MAXVAL(ISIZE),KNPTS,0:NPROC-1))
!
IF (NPROC>1) THEN
!for each task
DO JP=0,NPROC-1
#if defined(SFX_MPI) || defined(SFX_MNH)
!inind_all receives from all tasks the points they need that are
!located in it
CALL MPI_GATHER(ININD0(:,:,JP),MAXVAL(ISIZE)*KNPTS*KIND(ININD0)/4,MPI_INTEGER,&
ININD_ALL,MAXVAL(ISIZE)*KNPTS*KIND(ININD_ALL)/4,MPI_INTEGER,&
JP,NCOMM,INFOMPI)
#endif
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
ELSE
ININD_ALL(:,:,:) = ININD0(:,:,:)
ENDIF
!
DEALLOCATE(ININD0)
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_4',1,ZHOOK_HANDLE)
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_5',0,ZHOOK_HANDLE)
!
!zfield contains the values of the points needed located in this task
!(ie values for indexes of ININD_ALL)
ALLOCATE(ZFIELD(MAXVAL(ISIZE),KNPTS,SIZE(PFIELD,2),0:NPROC-1))
ZFIELD(:,:,:,:) = XUNDEF
DO JP=0,NPROC-1
DO JK=1,MAXVAL(ISIZE)
DO JL=1,KNPTS
IF (ININD_ALL(JK,JL,JP)/=0) THEN
!pfield in only on this task
ZFIELD(JK,JL,:,JP) = PFIELD(ININD_ALL(JK,JL,JP),:)
ENDIF
ENDDO
ENDDO
ENDDO
!
DEALLOCATE(ININD_ALL)
!
!ZFIELD2 gathers values needed for this task from all other tasks
!(inverse operation than before)
ALLOCATE(ZFIELD2(ICPT,KNPTS,SIZE(PFIELD,2),0:NPROC-1))
IF (NPROC>1) THEN
DO JP=0,NPROC-1
#if defined(SFX_MPI) || defined(SFX_MNH)
CALL MPI_GATHER(ZFIELD(1:ISIZE(JP),:,:,JP),SIZE(ZFIELD(1:ISIZE(JP),:,:,JP))*KIND(ZFIELD)/4,MPI_REAL,&
ZFIELD2,ISIZE(JP)*KNPTS*SIZE(PFIELD,2)*KIND(ZFIELD2)/4,MPI_REAL,JP,NCOMM,INFOMPI)
#endif
ENDDO
ELSE
ZFIELD2(:,:,:,:) = ZFIELD(:,:,:,:)
ENDIF
!
DEALLOCATE(ZFIELD)
!
!zfield3 contains the values of the points needed for interpolation, gathered from all tasks
ALLOCATE(ZFIELD3(ICPT,KNPTS,SIZE(PFIELD,2)))
DO JP=0,NPROC-1
WHERE (ZFIELD2(:,:,:,JP)/=XUNDEF) ZFIELD3(:,:,:) = ZFIELD2(:,:,:,JP)
ENDDO
DEALLOCATE(ZFIELD2)
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_5',1,ZHOOK_HANDLE)
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_6',0,ZHOOK_HANDLE)
!
!values of the points used for interpolation
ALLOCATE(ZNVAL(IP,IL2))
ZNVAL(:,:) = XUNDEF
!
!znval contains the averaged values for the knpts points
ZNVAL(:,:) = 0.
DO JL=1,ICPT
ZSUM = 0.
DO JP=1,KNPTS
IF (ININD(JL,JP)/=0) THEN
ZNVAL(JL,:) = ZNVAL(JL,:) + ZFIELD3(JL,JP,:)/ZNDIST(JL,JP)
ZSUM = ZSUM + 1./ZNDIST(JL,JP)
ENDIF
ENDDO
IF (ZSUM/=0.) ZNVAL(JL,:) = ZNVAL(JL,:) / ZSUM
ENDDO
!
DEALLOCATE(ININD, ZNDIST, ZFIELD3)
!
!
!finally, pfield contains the interpolated values!
ICPT=0
DO JL=1,IL1
IF (KCODE(JL)/=0) CYCLE
ICPT = ICPT + 1
PFIELD(JL,:) = ZNVAL(ICPT,:)
ENDDO
!
DEALLOCATE(ZNVAL)
!
DEALLOCATE(ISIZE)
!
IF (IOLD==2) THEN
#ifdef SFX_MNH
NRANK = IRANK_SAVE
NPROC = IPROC_SAVE
NPIO = IPIO_SAVE
NCOMM = ICOMM_SAVE
DEALLOCATE(IINDEX_TOT,INUM_TOT)
#endif
ENDIF
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS_6',1,ZHOOK_HANDLE)
!-------------------------------------------------------------------------------
!
END SUBROUTINE INTERPOL_NPTS