c     Subroutine for mapping CLM 1D output to 3D (elevationclass, lat, lon)
c     Written by Raymond Sellevold, for questions e-mail:
c     R.Sellevold-1@tudelft.nl
c     This subroutine can be used with NCL, for more info:
c     https://www.ncl.ucar.edu/Document/Tools/WRAPIT.shtml
c     ===============

C NCLFORTSTART
      subroutine pft2latlo(nt,nk,ltype,jxy,ixy,varpft,nla,nlo,nl,col,ou)
      integer nt, nk, nla, nlo, nl, col
      integer ltype(nk), jxy(nk), ixy(nk)
      dimension ou(nlo,nla,nl,nt), varpft(nk,nt)
      integer t, k
C NCLEND
      mylatidx = 0
      mylonidx = 0
      lev = 0

      do t = 1,nt
       do k = 1,nk
        if (ltype(k) == col) then
         if ((mylatidx.eq.jxy(k)).and.(mylonidx.eq.ixy(k))) then
          lev = lev+1
         else
          mylatidx = jxy(k)
          mylonidx = ixy(k)
          lev = 1
         end if
         ou(ixy(k),jxy(k),lev,t) = varpft(k,t)
         !ou(t,lev,jxy(k),ixy(k)) = varpft(t,k)
        end if
       enddo
      enddo

      return
      end



C NCLFORTSTART
      subroutine bilin_interp(cism_lat,cism_lon,clm_lonx,clm_lat,v,a)
      integer y, x, cism_ny, cism_nx
      integer clm_point_x1, clm_point_x2
      integer clm_point_y1, clm_point_y2
      real*8 cism_lon_x, cism_lat_y
      real*8 clm_lon_x1, clm_lon_x2, clm_lat_y1, clm_lat_y2
      real*8 dist, distA, distB, distC, distD
      real*8 valA, valB, valC, valD
      dimension dists(4), vars(4)

      dimension cism_lat(416,704), cism_lon(416,704)
      dimension clm_lonx(288), clm_lat(192)
      dimension v(288,192) 
      dimension a(416,704)
C NCLEND

      cism_ny = 704
      cism_nx = 416

      do x=1,cism_nx
       do y=1,cism_ny
c        First find the latitudes and longitudes of the CISM points
         cism_lon_x = cism_lon(x,y)
         cism_lat_y = cism_lat(x,y)

c        Find the indexes of the closest CLM longitude and latitude
         clm_point_x1 = minloc(abs(clm_lonx-cism_lon_x),1)
         clm_point_y1 = minloc(abs(clm_lat-cism_lat_y),1)

c        Use this information to locate x1 and y1
         if (clm_lonx(clm_point_x1).gt.cism_lon_x) then
          clm_point_x1 = clm_point_x1 - 1
         end if
         if (clm_lat(clm_point_y1).gt.cism_lat_y) then
          clm_point_y1 = clm_point_y1 - 1
         end if

c        Find x2 and y2
         clm_point_x2 = clm_point_x1 + 1
         clm_point_y2 = clm_point_y1 + 1

c        Assign longitudes and latitudes to x1,y1,x2,y2
         clm_lon_x1 = clm_lonx(clm_point_x1)
         clm_lon_x2 = clm_lonx(clm_point_x2)
         clm_lat_y1 = clm_lat(clm_point_y1)
         clm_lat_y2 = clm_lat(clm_point_y2)

c        Calculate the distances from CISM(x,y) to each CLM point
         dist = ((clm_lon_x2-clm_lon_x1)*(clm_lat_y2-clm_lat_y1))
         distA = ((clm_lon_x2-cism_lon_x)*(clm_lat_y2-cism_lat_y))
         distA = distA / dist
         distB = ((cism_lon_x-clm_lon_x1)*(clm_lat_y2-cism_lat_y))
         distB = distB / dist
         distC = ((clm_lon_x2-cism_lon_x)*(cism_lat_y-clm_lat_y1))
         distC = distC / dist
         distD = ((cism_lon_x-clm_lon_x1)*(cism_lat_y-clm_lat_y1))
         distD = distD / dist

c        Find the corresponding CLM values
         valA = v(clm_point_x1,clm_point_y1)
         valB = v(clm_point_x2,clm_point_y1)
         valC = v(clm_point_x1,clm_point_y2)
         valD = v(clm_point_x2,clm_point_y2)

c        Multiply with the distances
         valA = valA * distA
         valB = valB * distB
         valC = valC * distC
         valD = valD * distD

c        Mask not-a-number values and collect the values in an array named a
         if (abs(valA).gt. 1E+20) then
          valA = -9999.0
          distA = -9999.0
         endif
         if (abs(valB).gt. 1E+20) then
          valB = -9999.0
          distB = -9999.0
         endif
         if (abs(valC).gt. 1E+20) then
          valC = -9999.0
          distC = -9999.0
         endif
         if (abs(valD).gt. 1E+20) then
          valD = -9999.0
          distD = -9999.0
         endif
         vars = (/ valA, valB, valC, valD /)
         vars = pack(vars, vars.ne.-9999.0)
         dists = (/ distA, distB, distC, distD /)
         dists = pack(dists, dists.ne.-9999.0)

         if (sum(dists).eq.0) then
          a(x,y) = 1e+36
         else
          a(x,y) = sum(vars) / sum(dists)
         end if	
       enddo
      enddo

      return
      end


C NCLFORTSTART
      subroutine vert_interp(point,topo,values,valout)
      integer minp, upper, lower
      real*8 dist_upper, dist_lower
      real*8 point
      real valout
      dimension topo(10), values(10)
C NCLEND

c     Find the lower limit
      minp = minloc(abs(topo-point),1)

c     Decide if this point is the upper or lower limit of interpolation
      if (topo(minp).gt.point) then
       upper = minp
       lower = minp - 1
      else
       upper = minp + 1
       lower = minp
      end if

c     Check that the upper limit is higher than the lower limit
      if (topo(upper).gt.topo(lower)) then
       dist_upper = topo(upper) - point
       dist_lower = point - topo(lower)

       val_upper = values(upper) * dist_lower
       val_lower = values(lower) * dist_upper

       valout = (val_upper + val_lower) / (dist_upper + dist_lower)
      else
       valout = values(minp)
      end if

      return
      end


C NCLFORTSTART
      subroutine vert_interp2(topohr,topo,data_bilin,data_out)
      integer minp, upper, lower, x, y
      real*8 dist_upper, dist_lower
      real topohr(704), topo(704,10)
      real data_bilin(704,10), data_out(704)
C NCLEND

      do y=1,704

c       Find the lower limit
        minp = minloc(abs(topo(y,:)-topohr(y)),1)

c       Decide if this point is the upper or lower limit of interpolation
        if (topo(y,minp).gt.topohr(y)) then
         upper = minp
         lower = minp - 1
        else
         upper = minp + 1
         lower = minp
        end if

c       Check that the upper limit is higher than the lower limit
        if (topo(y,upper).gt.topo(y,lower)) then
         dist_upper = topo(y,upper) - topohr(y)
         dist_lower = topohr(y) - topo(y,lower)

         val_upper = data_bilin(y,upper) * dist_lower
         val_lower = data_bilin(y,lower) * dist_upper

         data_out(y) = (val_upper+val_lower)/(dist_upper+dist_lower)
        else
         data_out(y) = data_bilin(y,minp)
        end if
      enddo

      return
      end
