SUBROUTINE ncodaclim(latv, lonv, iyrv, imonv, idayv, ohcv, d20v, d26v, sstv, ierr)

! Version 2.0.0
! Last updated: 10 Mar 2020
! 
! Modifications: 
!       3/10/20 (SS) removed ioper flag
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------

      IMPLICIT NONE
      
      
!     Input variables

      REAL               :: latv		! query latitude (deg N)
      REAL               :: lonv		! query longitude (deg E)
      INTEGER            :: iyrv, imonv, idayv	! query year, month, day

!     Note: lonv can be 0 to 360 or 
!                       0 to 180 for east longitude, and -180 to 0 for west longitude

!     Output variables (NCODA climatology values for query lat/lon/date)

      REAL               :: ohcv		! climatological OHC		
      REAL               :: d20v		! climatological depth 20C
      REAL               :: d26v		! climatological depth 26C
      REAL               :: sstv		! climatological SST
      INTEGER            :: ierr		! error flag
     						!   0 = normal
						!   1 = error open/read climo file 
						!   2 = invalid date
						!   3 = invalid lon/lat
      
      
!     Internal variables

      REAL,    PARAMETER :: rmiss = -999.0
      INTEGER            :: i, j, k, ii, jj, kk
      
      INTEGER,                                SAVE :: nx, ny
      REAL,                                   SAVE :: dlon, dlat
      REAL,                                   SAVE :: lon0, lat0
      REAL,                                   SAVE :: rmissv
      
      REAL,    DIMENSION(:),     ALLOCATABLE, SAVE :: glon, glat    
      REAL,    DIMENSION(:,:,:), ALLOCATABLE, SAVE :: cohc, cd20, cd26, csst
      INTEGER,                                SAVE :: iread = 0      

      INTEGER                                      :: imon1, imon2
      REAL                                         :: wt1, wt2
      REAL                                         :: lonvt
    
      INTEGER, DIMENSION(12),                 SAVE :: dmon  


!     I/O

      INTEGER, PARAMETER :: luinp = 12, luout = 13
      CHARACTER(LEN=256) :: fninp, fnohc, fnd20, fnd26, fnsst  
      CHARACTER(LEN=256) :: coef_location
          

!     NCODA climatology file names

      DATA fnohc /'climOHC26C.dat'/
      DATA fnd20 /'climTOPO20.dat'/
      DATA fnd26 /'climTOPO26.dat'/
      DATA fnsst /'climSURF_T.dat'/

!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------

!     Length of each month in days, accounting for leap year

      dmon(1)  = 31
      dmon(2)  = 28
      dmon(3)  = 31
      dmon(4)  = 30
      dmon(5)  = 31
      dmon(6)  = 30
      dmon(7)  = 31
      dmon(8)  = 31
      dmon(9)  = 30
      dmon(10) = 31
      dmon(11) = 30
      dmon(12) = 31
      
      IF (MOD(iyrv, 4) .EQ. 0) dmon(2) = 29
      
      
!     Initialize ohcv as a missing value      
      
      ohcv = rmiss
      d20v = rmiss
      d26v = rmiss
      sstv = rmiss


!     Check that query date is valid

      IF (imonv .LT. 1 .OR. imonv .GT. 12 .OR. idayv .LT. 1 .OR. idayv .GT. dmon(imonv)) THEN
      
        ierr = 2
	
	RETURN
	
      ENDIF
            

!     Only read in climatological data the first time
  
      IF (iread .NE. 1) THEN
  
        CLOSE(luinp)
	
!       Read clim OHC data

        call getenv( "SHIPS_COEF", coef_location )
        fninp = trim( coef_location ) // trim( fnohc )
        print*,fninp

	OPEN(UNIT = luinp, FILE = fninp, STATUS = 'old', ERR = 900)
      
        READ(luinp,'(26X,2(I7,2F7.2),F7.1)', END = 910, ERR = 910) nx, lon0, dlon, ny, lat0, dlat, rmissv
	
	IF (ALLOCATED(cohc)) DEALLOCATE(cohc)
	
	ALLOCATE(cohc(nx,ny,12))
	
	cohc = rmiss
      
        DO k = 1, 12
      
          READ(luinp,'()', END = 910, ERR = 910)
	 
          DO j = ny, 1, -1
	
            READ(luinp,'(1440(f7.1))', END = 910, ERR = 910) (cohc(i,j,k),i=1,nx)
        
	  ENDDO
	  
        ENDDO
      
        CLOSE(luinp)
	
	
!       Read clim d20C data

        call getenv( "SHIPS_COEF", coef_location )
        fninp = trim( coef_location ) // trim( fnd20 )
	
	OPEN(UNIT = luinp, FILE = fninp, STATUS = 'old', ERR = 900)
      
        READ(luinp,'(26X,2(I7,2F7.2),F7.1)', END = 910, ERR = 910) nx, lon0, dlon, ny, lat0, dlat, rmissv
	
	IF (ALLOCATED(cd20)) DEALLOCATE(cd20)
	
	ALLOCATE(cd20(nx,ny,12))
	
	cd20 = rmiss
      
        DO k = 1, 12
      
          READ(luinp,'()', END = 910, ERR = 910)
	 
          DO j = ny, 1, -1
	
            READ(luinp,'(1440(f7.1))', END = 910, ERR = 910) (cd20(i,j,k),i=1,nx)
        
	  ENDDO
	  
        ENDDO
      
        CLOSE(luinp)
	

!       Read clim d26C data

        call getenv( "SHIPS_COEF", coef_location )
        fninp = trim( coef_location ) // trim( fnd26 )
	
	OPEN(UNIT = luinp, FILE = fninp, STATUS = 'old', ERR = 900)
      
        READ(luinp,'(26X,2(I7,2F7.2),F7.1)', END = 910, ERR = 910) nx, lon0, dlon, ny, lat0, dlat, rmissv
	
	IF (ALLOCATED(cd26)) DEALLOCATE(cd26)
	
	ALLOCATE(cd26(nx,ny,12))
	
	cd26 = rmiss
      
        DO k = 1, 12
      
          READ(luinp,'()', END = 910, ERR = 910)
	 
          DO j = ny, 1, -1
	
            READ(luinp,'(1440(f7.1))', END = 910, ERR = 910) (cd26(i,j,k),i=1,nx)
        
	  ENDDO
	  
        ENDDO
      
        CLOSE(luinp)
	
	
!       Read clim SST data

        call getenv( "SHIPS_COEF", coef_location )
        fninp = trim( coef_location ) // trim( fnsst )
	
	OPEN(UNIT = luinp, FILE = fninp, STATUS = 'old', ERR = 900)
      
        READ(luinp,'(26X,2(I7,2F7.2),F7.1)', END = 910, ERR = 910) nx, lon0, dlon, ny, lat0, dlat, rmissv
	
	IF (ALLOCATED(csst)) DEALLOCATE(csst)
	
	ALLOCATE(csst(nx,ny,12))
	
	csst = rmiss
      
        DO k = 1, 12
      
          READ(luinp,'()', END = 910, ERR = 910)
	 
          DO j = ny, 1, -1
	
            READ(luinp,'(1440(f7.1))', END = 910, ERR = 910) (csst(i,j,k),i=1,nx)
        
	  ENDDO
	  
        ENDDO
      
        CLOSE(luinp)		
	
	
	WHERE(cohc .EQ. rmissv) cohc = rmiss
	
	WHERE(cd20 .EQ. rmissv) cd20 = rmiss
	
	WHERE(cd26 .EQ. rmissv) cd26 = rmiss
	
	WHERE(csst .EQ. rmissv) csst = rmiss
	
	iread = 1
      
      
        IF (ALLOCATED(glon)) DEALLOCATE(glon)
	
	IF (ALLOCATED(glat)) DEALLOCATE(glat)
	
	ALLOCATE(glon(nx))
	
	ALLOCATE(glat(ny))      
      
        DO i = 1, nx
      
           glon(i) = lon0 + float(i-1)*dlon
	 
        ENDDO
      
        DO j = 1, ny
      
           glat(j) = lat0 + float(j-1)*dlat
	 
        ENDDO 

      ENDIF
      
      
!     Find the indices of the climatology grid point closest to lonv, latv      
      lonvt = lonv
      if (lonvt .lt. 0.0) lonvt = lonvt + 360.0
      
      ii = NINT((lonvt - lon0)/dlon) + 1
      jj = NINT((latv  - lat0)/dlat) + 1
      
      
!     Account for fact climatology grid is cyclic in the x direction
      
      IF (ii .EQ. 0) ii = nx
      
      
!     Error processing 
 
      IF (ii .LT. 1 .OR. ii .GT. nx) THEN
      
        ierr = 3
      
!        WRITE(6, '(A,F7.1)') '**Invalid input lon: ', lonv
	
	RETURN
	
      ENDIF

      IF (jj .LT. 1 .OR. jj .GT. ny) THEN
      
        ierr = 3
	
!	WRITE(6, '(A,F7.1)') '**Invalid input lat: ', latv
	
	RETURN
	
      ENDIF
      
      
!     Interpolate between monthly climatology values

      IF (idayv .LE. INT(dmon(imonv)/2)) THEN
      
        IF (imonv .NE. 1) THEN
	
	  imon1 = imonv - 1
	  
	ELSE
	
	  imon1 = 12
	  
	ENDIF
	
	imon2 = imonv
	
	wt1 = (REAL(dmon(imon2))/2.0 - REAL(idayv))
	wt1 = wt1 / (0.5*(REAL(dmon(imon1) + dmon(imon2))))
	
	wt2 = (REAL(dmon(imon1)) - REAL(dmon(imon1))/2.0 + REAL(idayv))
	wt2 = wt2 / (0.5*(REAL(dmon(imon1) + dmon(imon2))))
	
      ELSE
      
        IF (imonv .NE. 12) THEN
	
	  imon2 = imonv + 1
	  
	ELSE
	
	  imon2 = 1
	  
	ENDIF
	
	imon1 = imonv
	
	wt1 = (REAL(dmon(imon1)) - REAL(idayv) + REAL(dmon(imon2))/2.0)
	wt1 = wt1 / (0.5*(REAL(dmon(imon1) + dmon(imon2))))
	
	wt2 = (REAL(idayv) - REAL(dmon(imon1))/2.0)
	wt2 = wt2 / (0.5*(REAL(dmon(imon1) + dmon(imon2))))
	  
      ENDIF  
      
      
!     Compute interpolated climatological OHC, D20C, D26C, and SST at the closest grid point

      IF (cohc(ii,jj,imon1) .NE. rmiss .AND. cohc(ii,jj,imon2) .NE. rmiss) THEN
      
        ohcv = wt1*cohc(ii,jj,imon1) + wt2*cohc(ii,jj,imon2)
	
      ENDIF

      IF (cd20(ii,jj,imon1) .NE. rmiss .AND. cd20(ii,jj,imon2) .NE. rmiss) THEN
      
        d20v = wt1*cd20(ii,jj,imon1) + wt2*cd20(ii,jj,imon2)
	
      ENDIF
      
      IF (cd26(ii,jj,imon1) .NE. rmiss .AND. cd26(ii,jj,imon2) .NE. rmiss) THEN
      
        d26v = wt1*cd26(ii,jj,imon1) + wt2*cd26(ii,jj,imon2)
	
      ENDIF
      
      IF (csst(ii,jj,imon1) .NE. rmiss .AND. csst(ii,jj,imon2) .NE. rmiss) THEN
      
        sstv = wt1*csst(ii,jj,imon1) + wt2*csst(ii,jj,imon2)
	
      ENDIF
      
      ierr = 0
    
      
      RETURN
    
    
      900 CONTINUE
    
!      WRITE(6, '(A,A)') '**Error opening climatology file: ', TRIM(fninp)
      
      ierr = 1  
    
      RETURN
      
          
      910 CONTINUE
    
!      WRITE(6, '(A,A)') '**Error reading climatology file: ', TRIM(fninp) 
      
      ierr = 1 
    
      RETURN
      
    
END SUBROUTINE ncodaclim