c  This program reads in layer fields and outputs level fields

      parameter (imax=162)
      parameter (jmax=110)

      parameter (maxlev=401)
      dimension vert(maxlev)

      dimension rntk(4),rntkm(4)
      dimension phyk(4),phykm(4)
      dimension zook(4),zookm(4)
      dimension detk(4),detkm(4)

      dimension div(4)

C  Open input file (ASCII)
      open(10,file='layercode_data_71.asc',status='old',readonly)

C  Open output file (binary)
      open(11,file='levelcode.dat',status='new',form='unformatted')

      pi=4.*abs(atan(1.))
      rearth=6370.e5
      dypdeg=pi/180.*rearth
      ymin=-29.*dypdeg
      ymax=+25.*dypdeg
      xmin=0.*dypdeg
      xmax=+80.*dypdeg
      dx=(xmax-xmin)/(imax-2)
      dy=(ymax-ymin)/(jmax-2)

5     continue

c  Units:
c  djd is Julian day
c  w12a, w23a, w34a (entrainment velocities): cm/s
c  qsol (solar radiation multiplied by PAR): watts/m2
c  h1b,h2b,h3b,h4b: meters
c  All biological variables: micromole N/Kg
c  div: s^-1

      read(10,*,end=900)djd
      read(10,*)w12a,w23a,w34a,sol,qsol
      read(10,*)h1b,h2b,h3b,h4b
      read(10,*)rntk,rntkm
      read(10,*)phyk,phykm
      read(10,*)zook,zookm
      read(10,*)detk,detkm
      read(10,*)div

c  Interpolate velocity level fields from convergences
      do m=1,maxlev
c      Find the desired layer
       dep=-(m-1)
       dlay=(m-1)
       if(dlay.le.h1b) then
         lay=1
       elseif(dlay.le.(h1b+h2b)) then
         lay=2
       elseif(dlay.le.(h1b+h2b+h3b)) then
         lay=3
       elseif(dlay.le.(h1b+h2b+h3b+h4b)) then
         lay=4
       else
         lay=5
       endif

       if(lay.eq.1) then
         vert(m)=-dep*div(1)
       elseif(lay.eq.2) then
         vert(m)=h1b*div(1)-(dep+h1b)*div(2)
       elseif(lay.eq.3) then
         vert(m)=h1b*div(1)+h2b*div(2)-(dep+(h1b+h2b))*div(3)
       elseif(lay.eq.4) then
         vert(m)=h1b*div(1)+h2b*div(2)+h3b*div(3)-(dep+(h1b+h2b+h3b))
     .           *div(4)
         vert4=vert(m)
       elseif(lay.eq.5) then
         vert(m)=h1b*div(1)+h2b*div(2)+h3b*div(3)+h4b*div(4)
       endif
      enddo

      write(11)djd,qsol,h1b
      write(11)vert

      goto 5
900   close(10)
      close(11)

      end