!Reads RTOF's fields and set up interpolation 
!in POM's domain.
!Author: Biju Thomas  on 2016/01/12

MODULE readrtfinpmd
USE invdistinpmd
IMPLICIT NONE
  INTEGER, PARAMETER :: idm = 4500, jdm = 3298, ngdf = 12
  INTEGER, PARAMETER :: jdm1 = jdm-1, nstride = 50
  REAL, PARAMETER :: maskin = 1.2676506E+30, maskout = 1.2676506E+30
  
   CONTAINS
   SUBROUTINE readrtfinp(im, jm, kdm, lont, latt, lonu,latu, lonv,latv,    & 
                         zrtof, tpom, spom, upom, vpom, hpom)
     INTEGER, INTENT(IN) :: im, jm, kdm
     REAL, DIMENSION(idm, jdm) :: plon, plat, plon0, ulon, ulat, ulon0
     REAL, DIMENSION(idm, jdm) :: vlon, vlat, vlon0, hrf
     REAL, DIMENSION(idm, jdm, kdm) :: trf, srf, urf, vrf
     REAL, ALLOCATABLE :: f2d(:,:), f3d(:,:,:)   
     REAL, INTENT(IN), DIMENSION(im, jm) :: lont, latt, lonu, latu, lonv, latv
     REAL, INTENT(IN), DIMENSION(kdm) :: zrtof 
     REAL, INTENT(OUT), DIMENSION(im, jm) :: hpom
     REAL, INTENT(OUT), DIMENSION(im, jm, kdm) :: tpom, spom, upom, vpom
     REAL :: x, y, w11, w12, w21, w22
     REAL :: rflagged 
     INTEGER :: i1, j1, i2, j2
     INTEGER :: irs, ire, jrs, jre, irs0, ire0
     INTEGER :: recl
     INTEGER :: i, j, k
     IF ( .NOT. ALLOCATED(f2d)) ALLOCATE(f2d(idm,jdm))
     IF ( .NOT. ALLOCATED(f3d)) ALLOCATE(f3d(idm,jdm,kdm))
     INQUIRE (IOLENGTH = recl) 1.0
     recl = recl*idm*jdm
     OPEN(10,FILE='lonlatrtf.dat', STATUS='OLD', RECL=recl, ACCESS='DIRECT')
     READ(10,REC=1)f2d
     READ(10,REC=2)plat
     OPEN(11,FILE='temprtf.dat', FORM='UNFORMATTED')
     READ(11)f3d
     CLOSE(11)
     DO j = 1, jdm1
       DO i = 1, idm
         DO k=1, kdm
           trf(i,j,k) = f3d(i,j,k)
         ENDDO
         plon(i,j) = f2d(i,j)
       ENDDO
     ENDDO
     DO i = 1, idm
       DO k=1, kdm
          trf(i,jdm,k) = f3d(idm-i+1,jdm,k)
       ENDDO
       plon(i,jdm) = f2d(idm-i+1,jdm)
     ENDDO

     READ(10,REC=3)f2d
     READ(10,REC=4)f2d
   
     OPEN(12,FILE='salnrtf.dat', FORM='UNFORMATTED')
     READ(12)f3d
     CLOSE(12)
     DO j = 1, jdm1
       DO i = 1, idm
         DO k=1, kdm
           srf(i,j,k) = f3d(i,j,k)
         ENDDO
       ENDDO
     ENDDO
     DO i = 1, idm
       DO k=1, kdm
         srf(i,jdm,k) = f3d(idm-i+1,jdm,k)
       ENDDO
     ENDDO

     READ(10,REC=5)f2d
     READ(10,REC=6)ulat
     OPEN(13,FILE='uwndrtf.dat', FORM='UNFORMATTED')
     READ(13)f3d
     CLOSE(13)
     DO j = 1, jdm1
       DO i = 1, idm
         DO k=1, kdm
           urf(i,j,k) = f3d(i,j,k)
         ENDDO
       ulon(i,j) = f2d(i,j)
       ENDDO
     ENDDO
     DO i = 1, idm
       DO k=1, kdm
         urf(i,jdm,k) = f3d(idm-i+1,jdm,k)
       ENDDO
       ulon(i,jdm) = f2d(idm-i+1,jdm)
     ENDDO

     READ(10,REC=7)f2d
     READ(10,REC=8)vlat
     CLOSE(10)
     OPEN(14,FILE='vwndrtf.dat', FORM='UNFORMATTED')
     READ(14)f3d
     CLOSE(14)
     DO j = 1, jdm1
       DO i = 1, idm
         DO k=1, kdm
           vrf(i,j,k) = f3d(i,j,k)
         ENDDO
         vlon(i,j) = f2d(i,j)
       ENDDO
     ENDDO
     DO i = 1, idm
       DO k=1, kdm
         vrf(i,jdm,k) = f3d(idm-i+1,jdm,k)
       ENDDO
       vlon(i,jdm) = f2d(idm-i+1,jdm)
     ENDDO
   
     OPEN(10, FILE='depthrtf.dat', STATUS='OLD', RECL=recl, ACCESS='DIRECT')
     READ(10,REC=1) f2d
     DO j = 1, jdm1
       DO i = 1, idm
         hrf(i,j) = f2d(i,j)
       ENDDO
     ENDDO
     DO i = 1, idm
       hrf(i,jdm) = f2d(idm-i+1,jdm)
     ENDDO
     DEALLOCATE (f3d, f2d)
     DO k = 1, kdm
       DO j = 1, jdm
         DO i = 1, idm
           IF( zrtof(k) >  hrf(i,j) ) THEN
              trf(i,j,k) = maskin
              srf(i,j,k) = maskin
              urf(i,j,k) = maskin
              vrf(i,j,k) = maskin  
           ENDIF
         ENDDO
       ENDDO
     ENDDO

     irs0 = 1
     ire0 = idm
     jrs = 1
     jre = jdm

     DO j = 1,jm
     irs = irs0
     ire = ire0
     DO i = 1,im
       CALL nearpts(idm, jdm, plon, plat, lont(i,j), latt(i,j),  &
                    irs, ire, jrs, jre, trf(:,:,1), i1, j1, i2, j2, maskin )        

       CALL weights(i,j, idm, jdm, plon, plat, lont(i,j), latt(i,j),  &
                    trf(:,:,1), i1, j1, i2, j2, w11, w12, w21, w22, maskin)    

       CALL invdistinp(maskin,maskout,hrf,hpom(i,j),   &
                       i1, j1, i2, j2, w11, w12, w21, w22,idm,jdm)
       DO k =1, kdm
         CALL check_weights(maskin, trf(i1,j1,k), trf(i1,j2,k),   &
               trf(i2,j1,k),  trf(i2,j2,k), w11, w12, w21, w22)
         CALL invdistinp(maskin,maskout,trf(:,:,k),tpom(i,j,k),   &
                         i1, j1, i2, j2, w11, w12, w21, w22,idm,jdm)
         CALL invdistinp(maskin,maskout,srf(:,:,k),spom(i,j,k),   &
                         i1, j1, i2, j2, w11, w12, w21, w22,idm,jdm)
       ENDDO
       IF ( i == 1 ) THEN
         irs0 = MAX(1, i1-nstride)
         ire0 = MIN(idm, i2+nstride)
         irs = irs0
         ire = ire0
       ELSE
         irs = MAX(1, i1-nstride)
         ire = MIN(idm, i2+nstride)
       ENDIF
       jrs = MAX(1, j1-nstride)
       jre = MIN(jdm, j2+nstride)
     END DO
     END DO

     irs0 = 1
     ire0 = idm
     jrs = 1
     jre = jdm

     DO j = 1,jm
      irs = irs0
      ire = ire0
      DO i = 1,im
       CALL nearpts(idm, jdm, ulon, ulat, lonu(i,j), latu(i,j),  &
                    irs, ire, jrs, jre, trf(:,:,1), i1, j1, i2, j2, maskin )

       CALL weights(i,j, idm, jdm, ulon, ulat, lonu(i,j), latu(i,j),  &
                    trf(:,:,1), i1, j1, i2, j2, w11, w12, w21, w22, maskin)
       DO k =1, kdm
         CALL check_weights(maskin, urf(i1,j1,k), urf(i1,j2,k),   &
               urf(i2,j1,k), urf(i2,j2,k), w11, w12, w21, w22)
         CALL invdistinp(maskin,maskout,urf(:,:,k),upom(i,j,k),   &
                         i1, j1, i2, j2, w11, w12, w21, w22,idm,jdm)
       ENDDO
       IF ( i == 1 ) THEN
         irs0 = MAX(1, i1-nstride)
         ire0 = MIN(idm, i2+nstride)
         irs = irs0
         ire = ire0
       ELSE
         irs = MAX(1, i1-nstride)
         ire = MIN(idm, i2+nstride)
       ENDIF
       jrs = MAX(1, j1-nstride)
       jre = MIN(jdm, j2+nstride)

     END DO
     END DO

     irs0 = 1
     ire0 = idm
     jrs = 1
     jre = jdm

     DO j = 1,jm
       irs = irs0
       ire = ire0
       DO i = 1,im
        CALL nearpts(idm, jdm, vlon, vlat, lonv(i,j), latv(i,j),  &
                     irs, ire, jrs, jre, trf(:,:,1), i1, j1, i2, j2, maskin )

        CALL weights(i,j, idm, jdm, vlon, vlat, lonv(i,j), latv(i,j),  &
                     trf(:,:,1), i1, j1, i2, j2, w11, w12, w21, w22, maskin)
        DO k =1, kdm
          CALL check_weights(maskin, vrf(i1,j1,k), vrf(i1,j2,k),   &
               vrf(i2,j1,k), vrf(i2,j2,k), w11, w12, w21, w22)
          CALL invdistinp(maskin,maskout,vrf(:,:,k),vpom(i,j,k),   &
                          i1, j1, i2, j2, w11, w12, w21, w22,idm,jdm)
        ENDDO
        IF ( i == 1 ) THEN
          irs0 = MAX(1, i1-nstride)
          ire0 = MIN(idm, i2+nstride)
          irs = irs0
          ire = ire0
        ELSE
          irs = MAX(1, i1-nstride)
          ire = MIN(idm, i2+nstride)
        ENDIF
       jrs = MAX(1, j1-nstride)
       jre = MIN(jdm, j2+nstride)
      END DO
    END DO
    
    rflagged = maskout
    IF ( ABS(selmax(rflagged, im, jm, hpom))  > maskout*0.0001 ) &
      PRINT*, 'H INTP messed'
    DO k =1, kdm
      IF ( ABS(selmax(rflagged, im, jm, tpom(:,:,k)))  > maskout*0.0001 ) &
        PRINT*, 'T INTP L/S messed : ', k
      IF ( ABS(selmax(rflagged, im, jm, spom(:,:,k)))  > maskout*0.0001 ) &
        PRINT*, 'S INTP L/S messed : ', k
    ENDDO
       
  END SUBROUTINE readrtfinp

  FUNCTION selmax(maskval, nx, ny, fin) RESULT(smax)
    INTEGER, INTENT(IN) :: nx, ny
    REAL, INTENT(IN) :: maskval
    REAL, INTENT(IN) :: fin(nx,ny)
    REAL :: smax
    INTEGER :: i, j

    smax = -ABS(maskval)
    DO j = 1, ny
    DO i = 1, nx
      IF ( fin(i,j) /= maskval ) THEN
         IF ( fin(i,j) > smax ) smax = fin(i,j)
      ENDIF
    ENDDO
    ENDDO
  END FUNCTION selmax

END MODULE readrtfinpmd