SUBROUTINE oifield (xme, xmo, xoi, sigobs, lwave, sigmod, lmax, ndim2)
!**** *OIFIELD* - OPTIMUM INTERPOLATION.
!     PURPOSE.
!     --------
!       TO PRODUCE A MAP OF THE FIELD X, MERGING MEASUREMENT
!       AND MODEL , BY OPTIMUM INTERPOLATION. THE ARRAY XOI
!       AT THE END OF THE SUBROUTINE CONTAINS THE VALUES
!       TO BE USED TO ANALYSE THE SPECTRA, HAVING NEGATIVE RETURN
!       CODES WHERE O.I. PRODUCED NO RESULTS.
!       THE FIRST GUESS FIELD XMO IS NOT MODIFIED
!       IN THIS SUBROUTINE.
!**   INTERFACE.
!     ----------
!       *CALL* *OIFIELD (XME, XMO, XOI, SIGOBS, SIGMOD)*
!         *XME*     REAL     FIELD FROM MEASUREMENTS.
!         *XMO*     REAL     FIELD FROM MODEL (FIRST GUESS).
!         *XOI*     REAL     FIELD FROM O.I.
!         *SIGOBS*  REAL     MEASUREMENT SCATTER.
!         *LWAVE*   LOGICAL  TRUE IF ASSIMILATION FOR WAVE HEIGHTS WITH
!                            ADJUSTED OBSERVATIONAL ERRORS.
!         *SIGMOD*  REAL     MODEL SCATTER.
! ----------------------------------------------------------------------
use params
use gridpar
use map
use testo
implicit none

REAL   , INTENT(IN)    :: xme(ngx,ngy),sigobs, sigmod
REAL   , INTENT(IN OUT):: xmo(ngx,ngy)
REAL   , INTENT(OUT)   :: xoi(ngx,ngy)
LOGICAL, INTENT(IN)    :: lwave
INTEGER, INTENT(IN)    :: lmax, ndim2
!   LMAX = MAXIMUM NUMBER OF GRID STEPS TO SPREAD MEASURED HS AND USTAR
INCLUDE'globals.h'
INTEGER, DIMENSION(NDIM2)             :: IMEAS, JMEAS
REAL   , ALLOCATABLE, DIMENSION (:)   :: V, WORK
INTEGER, ALLOCATABLE, DIMENSION (:)   :: IPVT
REAL   , ALLOCATABLE, DIMENSION (:, :):: XM, D, P
REAL    :: SIGMAOBS (NDIM2)
INTEGER :: LLON(NGY)
REAL    :: COSLON, DELLON, COS2, SIN2, DIST, DIST2, HSMIN
REAL    :: A, B, DOBS, DOBS2, COND, EPS
INTEGER :: I, J, K, I1, I2, IG, IX, IJ, J1, J2, JX, JSN, IOBS, JOBS, NOBS, IER, LLAT
! ----------------------------------------------------------------------
!*    1. INITIALISE SEARCH DISTANCE AND CORRELATION LENGTH.
!     -----------------------------------------------------
IF (XDELLA.LT.0.5) THEN  
!       DIST  = 150000./R
   DIST = 400000. / R  
ELSE  
   DIST = 300000. / R  
ENDIF  
DIST2 = 2. * DIST  
!     DEFINE COEFFICIENT FOR THE ALT DATA ERROR
IF (LWAVE) THEN  
   A = SIGOBS  
   B = - LOG (0.001)  
   HSMIN = 0.01  
ELSE  
   A = 0.  
   B = 0.  
   HSMIN = 0.5  

ENDIF  
!       THE MEASUREMENTS ARE STORED IN THE O.I. FIELD.
!       NEGATIVE VALUES WILL REMAIN WHERE O.I. WILL PRODUCE NO RESULTS
DO J = 1, NGY  
DO I = 1, NGX  
XOI (I, J) = XME (I, J)  
ENDDO  
ENDDO  
DO JSN = 1, NGY  
SIN2 = SINPH (JSN) **2  
COS2 = COSPH (JSN) **2  
COSLON = (COS (DIST2) - SIN2) / COS2  
DELLON = ACOS (COSLON)  
LLON (JSN) = NINT (DEG * DELLON / XDELLO)  
LLON (JSN) = MIN (LMAX, LLON (JSN) )  
ENDDO  
LLAT = NINT (DEG * DIST2 / XDELLA)  
LLAT = MIN (LMAX, LLAT)  
!       LOOP OVER GRID POINTS.
!       ----------------------
DO IG = 1, IGL  

DO IJ = IJS (IG), IJL (IG)  
I = IXLG (IJ, IG)  
J = KXLT (IJ, IG)  
IF (XMO (I, J) .LT.0.01) GOTO 2001  
NOBS = 0  
DO IOBS = 1, NDIM2  
IMEAS (IOBS) = 0  
JMEAS (IOBS) = 0  

ENDDO  
!      FIND NUMBER AND POSITION OF OBSERVATIONS IN A BOX
!      OF SIZE 2*LMAX+1 AROUND ANALYSIS POINT I,J.
J1 = J - LLAT  
J1 = MAX (J1, 1)  
J2 = J + LLAT  
J2 = MIN (J2, NGY)  
DO JX = J1, J2  
I1 = MAX (1, I - LLON (JX) )  
I2 = MIN (NGX, I + LLON (JX) )  
DO IX = I1, I2  
IF (XME (IX, JX) .GT.0.) THEN  
!               WRITE(*,*) 'XME OIFIELD:',XME(IX,JX),IX,JX,I,J
   IF (XMO (IX, JX) .GE.0.01) THEN  
      NOBS = NOBS + 1  
      SIGMAOBS (NOBS) = SIGOBS + A * EXP ( - B * (XME (IX, JX) &
       - HSMIN) )
      IMEAS (NOBS) = IX  
      JMEAS (NOBS) = JX  
   ENDIF  
ENDIF  
ENDDO  
IF (IPER.EQ.1) THEN  
   I1 = 1  
   I2 = I + LLON (JX) - NGX  
   DO IX = I1, I2  
   IF (XME (IX, JX) .GT.0.) THEN  
      IF (XMO (IX, JX) .GE.0.01) THEN  
         NOBS = NOBS + 1  
         SIGMAOBS (NOBS) = SIGOBS + A * EXP ( - B * (XME (IX, JX) &
          - HSMIN) )
         IMEAS (NOBS) = IX  
         JMEAS (NOBS) = JX  
      ENDIF  
   ENDIF  
   ENDDO  
   I1 = I - LLON (JX) + NGX  
   I2 = NGX  
   DO IX = I1, I2  
   IF (XME (IX, JX) .GT.0.) THEN  
      IF (XMO (IX, JX) .GE.0.01) THEN  
         NOBS = NOBS + 1  
         SIGMAOBS (NOBS) = SIGOBS + A * EXP ( - B * (XME (IX, JX) &
          - HSMIN) )
         IMEAS (NOBS) = IX  
         JMEAS (NOBS) = JX  
      ENDIF  
   ENDIF  
   ENDDO  
ENDIF  
ENDDO  

IF (NOBS.EQ.0) GOTO 2001  
!         WRITE(*,*) 'NU. OF OBSERVATIONS:',NOBS
ALLOCATE (D (NOBS, NOBS) )  
ALLOCATE (P (NOBS, NOBS) )  

ALLOCATE (XM (NOBS, NOBS) )  
!         DETERMINE OBSERVATION CORRELATION MATRIX  D
DO IOBS = 1, NOBS  
DO JOBS = 1, NOBS  
D (IOBS, JOBS) = 0.  
ENDDO  
ENDDO  
DO IOBS = 1, NOBS  
D (IOBS, IOBS) = (SIGMAOBS (IOBS) / SIGMOD) **2  

ENDDO  
!         DETERMINE MODEL CORRELATION MATRIX AT OBSERVATION
!         POINTS.
DO IOBS = 1, NOBS  
DO JOBS = IOBS + 1, NOBS  
I1 = IMEAS (IOBS)  
I2 = IMEAS (JOBS)  
J1 = JMEAS (IOBS)  
J2 = JMEAS (JOBS)  
DELLON = REAL (I1 - 1) * XDELLO - REAL (I2 - 1) * XDELLO  
COSLON = COS (DELLON * RAD)  
DOBS2 = COSLON * COSPH (J1) * COSPH (J2) + SINPH (J1) * SINPH (J2)  
DOBS = ACOS (DOBS2)  
P (IOBS, JOBS) = EXP ( - DOBS / DIST)  
P (JOBS, IOBS) = P (IOBS, JOBS)  
ENDDO  
P (IOBS, IOBS) = 1.  

ENDDO  
DO IOBS = 1, NOBS  
DO JOBS = 1, NOBS  
XM (IOBS, JOBS) = P (IOBS, JOBS) + D (IOBS, JOBS)  
ENDDO  

ENDDO  
DEALLOCATE (D)  

DEALLOCATE (P)  
!         INVERSE MATRIX XM
COND = 0.  

EPS = 0.  
ALLOCATE (V (NOBS) )  
ALLOCATE (WORK (NOBS) )  
ALLOCATE (IPVT (NOBS) )  
DO IOBS = 1, NOBS  
V (IOBS) = 0.  

ENDDO  
  200 CONTINUE  
CALL SYMINV (XM, NOBS, V, WORK, IPVT, IER)  
IF (IER.EQ.1) THEN  
   COND = 0.1  
   WRITE (IU06, * ) 'OIFIELD: MATRIX XM IS SINGULAR!', &
    'INCREASING MAIN DIAGONAL BY:', COND
   DO K = 1, NOBS  
   XM (K, K) = XM (K, K) + COND  
   ENDDO  
   GOTO 200  

ENDIF  
DEALLOCATE (V)  
DEALLOCATE (WORK)  



DEALLOCATE (IPVT)  
!         PRODUCE FIELD AT POINT I,J ACCORDING TO OPTIMUM INTERPOLATION

CALL ANALYSE (XM, NOBS, XME, XMO, XOI (I, J), IMEAS, JMEAS, I, J, &
 DIST, NDIM2)

DEALLOCATE (XM)  
 2001 CONTINUE  
ENDDO  

ENDDO  
RETURN  
END SUBROUTINE oifield
