#define ITYPE_RST 0
#define MACHINE_WORD 4
#define c_str(s) ('s\0')
#define n_a(s,i) ('s'//numbers(i:i)//'\0')

      subroutine init_rstrt (nxp, nyp, nz, npt, zin, iseg)
c---------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      character*100 str
      common /run/   nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast
      common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100)
      dimension   zin(1), zz(50), iseg(1)

      equivalence (inf(20),kd_sta),
     *            (inf(21),kd_xy),  (inf(22),kd_z),  (inf(23),kd_time), 
     *            (inf(24),kd_seg), (inf(25),kd_uv), (inf(26),kd_temp), 
     *            (inf(27),kd_salt),(inf(28),kd_tr), (inf(29),kd_phi), 
     *            (inf(30),kd_mask),(inf(32),kd_psi), (inf(33),kd_ice), 
     *            (inf(100),kd_end)

      call opda (1, 1, fbo(1:n_out))
      call opda (0, 1, fbo0(1:n_out0))

      do i = 1, 100
         inf(i) = 0
         rnf(i) = 0.
      enddo

      iword   = MACHINE_WORD
      
      kd = iword * 200
      str = 'Restart file for '//finp(1:mlen(finp)) 
      call wrda (1, 1, kd, len(str), str)
      call wrda (0, 1, kd, len(str), str)

      kd = iword * 200 + 100
      kd_sta = kd
      kd_xy  = kd
      kd_keep = kd_xy
      call wrda (1, iword, kd, nxp, xm) 
      call wrda (1, iword, kd, nyp, ym) 
      kd_xy = kd_keep
      call wrda (0, iword, kd, nxp, xm) 
      call wrda (0, iword, kd, nyp, ym) 

      kd_z = kd
      kd_keep = kd_z
      do i = 1, nz
         zz(i) = -zin(i)
      enddo
      call wrda (1, iword, kd, nz, zz) 
      kd_z = kd_keep
      call wrda (0, iword, kd, nz, zz) 

      call segm_from_iox (npt, nseg, iox, iseg)
      kd_seg = kd
      kd_keep = kd_seg
      call wrda (1, iword, kd, 2*nseg, iseg) 
      kd_seg = kd_keep
      call wrda (0, iword, kd, 2*nseg, iseg) 
      
      kd_uv = kd

      nstart     = 1
      nlaststart = 1
      nlast      = nsteps
      nscpu = 0
      nswll = 0

      inf(1) = 120
      inf(2) = nxp
      inf(3) = nyp
      inf(4) = nz
      inf(5) = 1
      inf(6) = npt
      inf(7) = nseg
      inf(8) = ntrac_tot

      inf(11) = nstart

      inf(15) = nlaststart
      inf(16) = nsteps
      inf(17) = nlast

      rnf(1) = delt
      rnf(2) = enso_start
      rnf(3) = enso_scale
      rnf(4) = rnf(2)

      budget_sp = 0.
      budget_ri = 0.
      budget_re = 0.
      budget_eb = 0.
      budget_pr = 0.
      coef_bud = coef_precip

      return
      end

      subroutine read_rstrt (nxp, nyp, nz, npt)
c--------------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_bio.h'
      include 'diffiso.h'
      common /run/   nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast
      common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100)
      common /mean_comm/ nmcount, nmcount_bio, nmcount_mosf, nmcount_conv
      character*4 mname(12)
      save mname

      equivalence (inf(20),kd_sta),
     *            (inf(21),kd_xy),  (inf(22),kd_z),  (inf(23),kd_time), 
     *            (inf(24),kd_seg), (inf(25),kd_uv), (inf(26),kd_temp), 
     *            (inf(27),kd_salt),(inf(28),kd_tr), (inf(29),kd_phi), 
     *            (inf(30),kd_mask),(inf(32),kd_psi), (inf(33),kd_ice), 
     *            (inf(100),kd_end)

      data mname /'Jan.','Feb.','Mar.','Apr.','May ','Jun.',
     *            'Jul.','Aug.','Sep.','Oct.','Nov.','Dec.'/

      iword = MACHINE_WORD

      call opda (1, 0, fbi(1:n_in))
      kd = 0
      call rdda (1, iword, kd, 100, inf) 
      kd = irdda (iword, 100, rnf) 

      ntype = idig(inf(1),1)
      if (ntype .ne. ITYPE_RST) then
         write (6, *) 'The file <' ,fbi(1:n_in),'> is not a RESTART-TYPE one !!'
         stop
      endif
      
      if (nxp.ne.inf(2) .or. nyp.ne.inf(3) .or.
     *     nz.ne.inf(4) .or. npt.ne.inf(6) ) then
         write(6,*)'The model dimensions are different from the restart data!'
         stop
      endif

      nptz   = inf(4)*inf(6)
      mtra   = inf(8)
      nptzp   = inf(4)*inf(6) +inf(6)

      nruns  = inf(12)
      nmcount      = 0
      nmcount_bio  = 0
      nmcount_mosf = 0
      nmcount_conv = 0

      budget_sp = 0.
      budget_ri = 0.
      budget_eb = 0.
      budget_pr = 0.
      coef_bud  = coef_precip
      budget_re = 0.
      if (irest.eq.1) coef_bud  = rnf(19)
      f0 = rnf(5)

      if (irest.ne.4) then
         nstart = inf(11) + 1
         nscpu  = inf(13)
         nswll  = inf(14)
         if (irest .eq. 1) then
            nlaststart = inf(15)
         else
            nlaststart = nstart
         endif
         nlast = nlaststart + nsteps - 1
         enso_start = rnf(2)
         rnf(4) = rnf(2) + rnf(3)*nstart
      else
         nstart     = 1
         nlaststart = 1
         nlast      = nsteps
         nscpu = 0
         nswll = 0
         inf(11) = nstart
         
         inf(15) = nlaststart
         inf(16) = nsteps
         inf(17) = nlast
         
         rnf(2) = enso_start
         rnf(4) = rnf(2)
      endif

      if (irest .lt. 3 .and. irest. ge. 0) then
         ekf1 = rnf(10)
         epf1 = rnf(11)
         hcf1 = rnf(12)
         wcf1 = rnf(13)
         vlf1 = rnf(14)
         do k = 1, inf(4)+1
            hsave(k) = rnf(20+k)
         enddo
      endif

      kd = kd_uv
      call rdda (1, iword, kd, nptz, u) 
      kd = irdda (iword, nptz, v) 
      kd = irdda (iword, nptz, uc) 
      kd = irdda (iword, nptz, vc) 
      kd = irdda (iword, nptz, h) 
      
      if (kd_temp .ne. 0) then
         kd = kd_temp
         call rdda (1, iword, kd, nptz, t) 
      endif
      if (ibaro.ne.0.and.kd_psi .ne. 0) then
         kd = kd_psi
         call rdda (1, iword, kd, npt, psi) 
      endif
      if (use_ice .and. kd_ice .ne. 0) then
         kd = kd_ice
         call rdda (1, iword, kd, npt, trs(i_hice)) 
         kd = irdda (iword, npt, trs(i_cice)) 
         kd = irdda (iword, npt, trs(i_thice)) 
      endif
      if (kd_salt .ne. 0) then
         kd = kd_salt
         call rdda (1, iword, kd, nptz, sal) 
         kd = irdda (iword, nptz, dens) 
      endif
      if ((use_trac .or.use_bio) .and. kd_tr .ne.   0) then
         kd = kd_tr
         call rdda (1, iword, kd, nptz*mtra, tr) 
      endif
     
      call clda(1)

      call opda (1, 1, fbo(1:n_out))
      call opda (0, 1, fbo0(1:n_out0))

      inf(8)  = ntrac_tot
      inf(12) = nruns + 1
      inf(15) = nlaststart
      inf(16) = nsteps
      inf(17) = nlast

      rnf(1) = delt

      call enso2date (rnf(2), id, im, iy)
      write (iout, 101) 'Model Starting Date:        ', mname(im), id, iy
      call enso2date (real(rnf(2)+enso_scale*nlast), id, im, iy)
      write (iout, 101) 'Scheduled End of the Run:   ', mname(im), id, iy
      call enso2date (real(rnf(2)+enso_scale*nstep), id, im, iy)
      write (iout, 101) 'Current Time:   ' , mname(im), id, iy

      call enso2res (real(rnf(3)*inf(17)), id, im, iy)
      write (iout, 201) 'Scheduled Length of the Run:', iy, im, id, inf(17)
      call enso2res (real(rnf(3)*inf(18)), id, im, iy)
      write (iout, 201) 'Elapsed Model Time:         ', iy, im, id, inf(18)

      write (iout, *) 'Number of restarts:', inf(12)

  101 format (a30, a4, i2, ',' , i5)
  201 format (a30, i5, ' Years', i3, ' months', i3, ' days. (', i7, ' steps.)') 

      return
      end

      subroutine close_rstrt 
c---------------------------------
      call clda(1)
      return
      end

      subroutine dump_rstrt
c------------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      character*100 str
      character*4 mname(12)
      common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100)
      save mname

      data mname /'Jan.','Feb.','Mar.','Apr.','May ','Jun.',
     *            'Jul.','Aug.','Sep.','Oct.','Nov.','Dec.'/

      iword = MACHINE_WORD
      
      call opda (1, 0, fbi(1:n_in))
      kd = 0
      call rdda (1, iword, kd, 100, inf) 
      kd = irdda (iword, 100, rnf)       
      kd = irdda (1, 100, str)       

      call clda(1)

      write (6, *) 'Dump of the data file <',fbi(1:n_in),'>:'
      write (6, *) 'LABEL: ', str
      write (6, *) 'TYPE', inf(1)
      write (6, *) 'NX =', inf(2)
      write (6, *) 'NY =', inf(3)
      write (6, *) 'NZ =', inf(4)
      write (6, *) 'NT =', inf(5)
      write (6, *) 'NPACK =', inf(6)
      write (6, *) 'NSEGM =', inf(7)
      write (6, *) 'NTRACERS =', inf(8)

      call enso2date (rnf(2), id, im, iy)
      write (6, 101) 'Model Starting Date:        ', mname(im), id, iy
      call enso2date (real(rnf(2)+rnf(3)*inf(17)), id, im, iy)
      write (6, 101) 'Scheduled End of the Run:   ', mname(im), id, iy
      call enso2date (real(rnf(2)+rnf(3)*inf(18)), id, im, iy)
      write (6, 101) 'Current End of the Run:   ' , mname(im), id, iy

      call enso2res (real(rnf(3)*inf(17)), id, im, iy)
      write (6, 201) 'Scheduled Length of the Run:', iy, im, id, inf(17)
      call enso2res (real(rnf(3)*inf(18)), id, im, iy)
      write (6, 201) 'Elapsed Model Time:         ', iy, im, id, inf(18)

      write (6, *) 'Number of restarts:', inf(12)

      call enso2date (rnf(4), id, im, iy)
      write (6, 101) 'Last Re-start:              ', mname(im), id, iy
      call enso2date (real(rnf(2) + rnf(3)*inf(15)), id, im, iy)
      write (6, 101) 'Last New Run:               ', mname(im), id, iy
      call enso2date (real(rnf(2) + rnf(3)*inf(11)), id, im, iy)
      write (6, 101) 'Last Save:                  ', mname(im), id, iy

      i = inf(13)
      if (i/3600 .eq. 1) then
         write (6, 301) 'Total CPU  time:', i/3600, mod(i,3600)/60, mod(i,60) 
      else
         write (6, 302) 'Total CPU  time:', i/3600, mod(i,3600)/60, mod(i,60) 
      endif
      i = inf(14)
      if (i/3600 .eq. 1) then
         write (6, 301) 'Total WALL time:', i/3600, mod(i,3600)/60, mod(i,60) 
      else
         write (6, 302) 'Total WALL time:', i/3600, mod(i,3600)/60, mod(i,60) 
      endif

  101 format (a30, a4, i2, ',' , i5)
  201 format (a30, i5, ' Years', i3, ' months', i3, ' days. (', i7, ' steps.)') 
  301 format (a20, i5, ' hour ', i3, ' min', i3, ' sec.')
  302 format (a20, i5, ' hours', i3, ' min', i3, ' sec.') 

      return
      end

      subroutine segm_from_iox (npt, nseg, iox, iseg)
c------------------------------------------------------
      dimension iox(1), iseg(2,1)
      
      ista   = iox(1)
      inext  = ista + 1
      nseg   = 0
      icount = 0

      do i = 2, npt
         icurr = iox(i)
         icount = icount + 1

         if (inext .ne. icurr) then
            nseg = nseg + 1
            iseg(1,nseg) = ista
            iseg(2,nseg) = ista + icount - 1
            ista  = icurr
            icount = 0
         endif

         inext = icurr + 1
      enddo

      nseg = nseg + 1
      iseg(1,nseg) = ista
      iseg(2,nseg) = ista + icount 

      return
      end

c------------------------------------------------------------
      subroutine model_input(npt,mbox)
c------------------------------------------------------------
c.....read model parameters and grid from file.
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_para.h'
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_bio.h'
      include 'biology.h'
      include 'amlice.h'
      include 'diffiso.h'
      include 'comm_topo.h'

      character*50 ubeg,vbeg,hbeg,tbeg,uout,vout,hout,tout,eout,
     +             blk,wndx,wndy,cwndx,cwndy,cloud,ccloud,wout
      common/files/ubeg,vbeg,hbeg,tbeg,uout,vout,hout,tout,eout,wout
      common/run/  nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast
      common/param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20),
     +              itherm,mlc,limp
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch
      common/winds/mtau,matau,tausc,atau,froude
      common/wnfils/wndx,wndy,cwndx,cwndy,cloud,ccloud
      common /vert/  zin(MAXNZ+1), hin(MAXNZ), t_in(MAXNZ+1), s_in(MAXNZ+1),
     *               bint(MAXNZ), cint(MAXNZ), dzin(MAXNZ+1), sigma(MAXNZ),
     *               facz(MAXNZ), ttr_fac(MAXNZ)
      common/strech/xs(MAXXS),alpha(MAXXS),beta(MAXXS)
      common /errors/ ioerr, nstep
      common /new_filt/  MAXFO, nbxk, nbyk, nfxk, nfpxk, nfyk



c----------------------------------------------------------------------------
      common /pbl_param/ pbl_pnu,pbl_delta,pbl_pml,pbl_depth,pbl_betav,
     *  pbl_grad,nstep_pbl, ipbl_advec, ipbl_jsta,ipbl_jend, 
     *  pbl_south,pbl_north, pbl_wmin
c----------------------------------------------------------------------------
      character*80       f_bar
      common /baro_files/ n_bar, f_bar
      integer iflag(12), RM12_NN 
      real    aflag(12)
      common /y12m_input/ RM12_NN, aflag, iflag 
      common /baro_input/ n_def_cor, mod_scheme, mod_solver, BAR_DEPMIN,
     *                    BAR_DSINK, ibar_key, nbaro, rayl, nonlin_baro

      character*1 BC_I(0:9)
      common /masks/BC_I

      logical use_per_island, use_stan_island, use_sunk
      common /baro_island/
     *                alons_min(10),alons_max(10),alats_min(10),alats_max(10)
     *               ,alon1_min,alon1_max,alat1_min,alat1_max
     *               ,per_lat,use_sunk,use_per_island,use_stan_island
     *               ,itrans,ftrans

      common /new_island/n_islands,in_rhs(0:9),ityp_island(0:9),
     *               alon_island(9),alat_island(9),btrans(0:9),
     *               n_poly(9),alon_poly(100,9),alat_poly(100,9),
     *               i_poly(100,9),j_poly(100,9)
c--------------------------------------------------------------

      logical sink_region
      common /bathymetry/sink_region,sr_depth, alon_sr, blon_sr, 
     *                   alat_sr, blat_sr

      logical use_hi
      common /order/ use_hi
      common/friction/ ibfric,b_fric, sw_fric
c----------------------------------------------------------------------------

      real inp_flt, inp_days
      logical inp_def
      dimension flt(100)
      save         flt
      character*80 fbdir

      external time, ctime
      integer time
      character*24 ctime
      character*30 date,pwd
      character*10 name
c----------------------------------------------------------------------------
      call inp_file(finp)
      call inp_vrnt(c_str(Variant), inp_int(c_str(Use_variant), -1))
      if (inp_int(c_str(Trace),0) .ne. 0) 
     *     call inp_trace(finp(1:mlen(finp))//'.tr\0')

      call version (date,pwd,name)
      iret = inp_str(c_str(Compile_date), date//'\0', date)
      iret = inp_str(c_str(Compile_dir), pwd//'\0', pwd)
      iret = inp_str(c_str(Compile_who), name//'\0', name)

      n_debug = inp_int(c_str(Debug),0)

      n_dir= inp_str(c_str(Output_dir),c_str(output), fbdir)
 
      n_tios= inp_str(c_str(Base_file),
     *     fbdir(1:n_dir)//'/'//finp(1:mlen(finp))//'\0',fbt)
      n_out = inp_str(c_str(Save_file), fbt(1:n_tios)//'.save\0', fbo)
      n_out0 = inp_str(c_str(Save_file0), fbt(1:n_tios)//'.bak\0', fbo0)
      n_in  = inp_str(c_str(Restart_file), fbo(1:n_out)//'\0', fbi)

      n_tau = inp_str(c_str(Wind_file),  c_str(wind_data), fbtau)

      n_tem = inp_str(c_str(Temp_file) , c_str(temp_data), fbtem)
      n_sal = inp_str(c_str(Salt_file),  c_str(salt_data), fbsal)
      n_psi = inp_str(c_str(Psi_file),   c_str(psi_data), fbpsi)

      n_sst = inp_str(c_str(SST_file),   c_str(sst_data), fbsst)
      n_sss = inp_str(c_str(SSS_file),   c_str(sss_data), fbsss)

      n_cld = inp_str(c_str(Cloud_file), c_str(cld_data), fbcld)
      n_slr = inp_str(c_str(Solar_file), c_str(slr_data), fbslr)
      n_prp = inp_str(c_str(Precip_file),    c_str(prp_data), fprp)

      n_hcl = inp_str(c_str(MxlTcl_file),c_str(mltc_data), fbhcl)

      n_wsp = inp_str(c_str(Wndspd_file),c_str(wndspd_data), fwsp)
      use_wndspd = fwsp(1:n_wsp).ne.'NONE'

      n_uwd = inp_str(c_str(Uwnd_file),c_str(uwnd_data), fuwd)
      n_vwd = inp_str(c_str(Vwnd_file),c_str(vwnd_data), fvwd)
      n_ah  = inp_str(c_str(Airhum_file),c_str(ahum_data), fah)
      n_at  = inp_str(c_str(Airtem_file),c_str(atem_data), fat)

      n_dep = inp_str(c_str(Bath_file),  c_str(bath_data), fbdep)
      n_sp = inp_str(c_str(Sponge_file), c_str(sponge_data), fbsp)
      n_riv = inp_str(c_str(River_file),  c_str(river_data), fbriv)

      n_grad = inp_str(c_str(Grad_file),  c_str(grad_data), fbgrad)

      if ( inp_def(c_str(Map_file)) ) 
     *     n_map = inp_str(c_str(Map_file), c_str(map_data), fbmap)

      irest = inp_int(c_str(Restart), 0)        ! 0 - start; 1 - restart

      lev_err = inp_int(c_str(Error_level), 1)

      ncyc  = inp_int(c_str(Ncycles),  4)

      delt = inp_days (c_str(Time_step), 1./24.)
      t_fac = inp_flt(c_str(Time_fac), 1.)
      bt_fac = inp_flt(c_str(Time_fac_baro), t_fac)
      stpd = 1./delt

      tr_zmin = inp_flt(c_str(Time_trace_zmin), 10000.)
      tr_zmax = inp_flt(c_str(Time_trace_zmax), 0.)
      tr_fac  = inp_flt(c_str(Time_trace_fac), 1.)

      mtmp = stpd - ncyc
      if (mtmp .lt. 0) mtmp = 0
      nsteps = nint(stpd * inp_days(c_str(Run_time),   365.)) + mtmp
      nskip  = nint(stpd * inp_days(c_str(Save_step),    5.))
      nskip = max(4,nskip - mod(nskip,4))
      nergy  = nint(stpd * inp_days(c_str(Energy_step), 1.))

      call inp_date (c_str(Starting_date), 1,1,2000, imon,iday,iyear)
      enso_start = date2enso (iday, imon, iyear)
      enso_scale = delt * 12./365.

      if (irest .eq. 0) then
         open (unit = iout, file = fout)
         if (lev_err .ge. 1) then
            write (iout, *) 'Run with control file <',finp(1:mlen(finp)),'>'
            write (iout, *) 'STARTED: ', ctime(time())
         endif
      else
         open (unit = iout, file = fout, access='APPEND')
         write (iout, *)'---------------------------------------------------' 
         if (lev_err .ge. 1) then
            write (iout, *) 'RE-STARTED: (code=',irest,' at:', ctime(time())
         endif
      endif


      land  = 0
      iglob = inp_int(c_str(Periodic), 0)
      iperx = (iglob.eq.1)
      mbc   = inp_int(c_str(BC_type),  2)                  
      if (mbc.eq.1.and.ibaro.ne.0) then
         print*,'using no-slip b.c. and barotropic, retry with free-slip b.c.'
      endif

      bc_time = inp_days (c_str(BC_coef), 30.) 
      bc_coef = 0.
      if (bc_time.ne.0.) bc_coef = delt/bc_time

      mseg  = inp_int(c_str(Discretization),  2) 
      use_hi  = (mseg .eq. 4)

      NXP = inp_int(c_str(NX), 10)
      NYP = inp_int(c_str(NY), 10)
      NZ  = inp_int(c_str(NZ), 2)
      nsig= inp_int(c_str(Nsigma), 0)
      if (nsig.eq.1.or.nsig.eq.2.or.nsig.gt.nz.or.nsig.lt.0)then
         print*,'number of sigma layers must be 0 or between 3 and nz'
         stop
      endif

      if(max0(nxp,nyp).gt.MAXSID) call wspace('MAXSID', max0(nxp,nyp))
      if(NZ .gt. MAXNZ)           call wspace('MAXNZ', NZ)

      alat = inp_flt(c_str(South), 0.)
      blat = inp_flt(c_str(North), real(nyp-1))
      alon = inp_flt(c_str(West), 0.)
      blon = inp_flt(c_str(East), real(nxp-1))

      mgrid = inp_int(c_str(Grid_type), 1)
      mbox  = inp_int(c_str(Box_type), 1)
      
      nsx = inp_rarr(c_str(X_stretch_pnts), 0, xs,  xs)
      nsx = inp_rarr(c_str(X_alpha),        0, alpha, alpha)
      nsx = inp_rarr(c_str(X_beta),         0, beta,  beta)

      nystrch = inp_int(c_str(Y_stretch_type), 1)
      if (nystrch.eq.1) then
         nsy = inp_rarr(c_str(Y_stretch_pnts), 0, xs(nsx+1),    xs(nsx+1))
         nsy = inp_rarr(c_str(Y_alpha),        0, alpha(nsx+1), alpha(nsx+1))
         nsy = inp_rarr(c_str(Y_beta),         0, beta(nsx+1),  beta(nsx+1))
         if(nsx+nsy.gt.maxxs)         call wspace('MAXXS', nsx+nsy)
      endif

      ipole = inp_int(c_str(Pole_shift), 0)
      if (ipole .eq. 1) then
         pole_alp = inp_flt(c_str(Pole_alpha), 0.)
         pole_bet = inp_flt(c_str(Pole_beta), -90.)
         pole_gam = inp_flt(c_str(Pole_gamma), 0.)
         call comp_rotma (pole_alp, pole_bet, pole_gam)
      endif

      b_fric  = inp_flt(c_str(Bottom_friction),  .001)
      sw_fric  = inp_flt(c_str(Sidewall_friction),  b_fric)
      ibfric  = inp_int(c_str(Bottom_draglaw), 1)

      dep_min = inp_flt(c_str(Bath_min), 100.)
      ibaro = inp_int(c_str(Baro_solv), 0)
      if (ibaro .ne. 0) call baro_init

      nonlin= inp_int(c_str(Nonlin), 1)
      if (nonlin.ne.1) then
         print*, 'Nonlin disabled (ie, always nonlinear), please do not use'
         stop
      endif

      itemp  = inp_int(c_str(Thermo),   1)
      if (itemp.eq.0) then
         print*, 'Thermo = 0 disabled'
         stop
      endif

      idens  = inp_int(c_str(Density), 1)

      isalt  = inp_int(c_str(Salinity), 4)
      isolrp = inp_int(c_str(Solar_penrad), 0)
      if (isolrp .eq. 1) solr_gamma = inp_flt(c_str(Solar_gamma), .333333)

      icl_h    = inp_int(c_str(Clim_H), 0)   ! 0-init from hin;1-Cnst;2-Vary;
      icl_psi    = inp_int(c_str(Clim_psi),0) ! 0-don't use;1-Cnst;2-Vary;

      if (icl_h .ne. 0.and.nsig.eq.0) then
         print*,'must use Nsigma > 2 if Clim_H <> 0'
         stop
      endif
      icl_htop = inp_int(c_str(Clim_htop), 0)  ! 0-Const; 1-Vary;
      
      icl_ts   = inp_int(c_str(Clim_TS),  0) ! 0-init from Tin;1-Cnst;2-Vary;
      nz_fix   = inp_int(c_str(Clim_nz),  nz+1) 
      icl_fix  = inp_int(c_str(Clim_nz_type),  1) 
      if (icl_h.eq.2 .and. icl_ts.ne.2) 
     *     write (ioerr, *) 'Warning: check Clim_H vs. Clim_TS...'
      if (icl_ts.eq.0 .and. initt.gt.1) 
     *     write (ioerr, *) 'Warning: to initialize temp, use icl_ts <> 0'
      if (icl_ts.eq.0 .and. inits.gt.0) 
     *     write (ioerr, *) 'Warning: to initialize salt, use icl_ts <> 0'
      icl_rlx  = inp_int(c_str(Clim_relax), 0)    ! 0-no relax; 1-relax N-S;
      ksponge  = inp_int(c_str(Sponge_width), 5) 
      krelax   = inp_int(c_str(Relax_width), 5) 

      clm_time = inp_days (c_str(Sponge_time), 30.) 
      clm_time_v = inp_days (c_str(Sponge_time_v), 10000.) 
      clm_time_tr = inp_days (c_str(Sponge_time_tr), 0.) 
      clm_time_trb = inp_days (c_str(Sponge_time_trb), 0.) 

      clm_time_psi = inp_days (c_str(Sponge_time_psi), 30.) 

      clm_coef = ncyc*delt/clm_time
      clm_coef_v = clm_time/clm_time_v
      
      clm_coef_tr  = 0.
      clm_coef_trb = 0.
      if (clm_time_tr.ne.0) clm_coef_tr = ncyc*delt/clm_time_tr
      if (clm_time_trb.ne.0) clm_coef_trb = ncyc*delt/clm_time_trb

      clm_psi = clm_coef* nbaro*clm_time/clm_time_psi

      clm_no   = inp_flt(c_str(Clim_nlat),  90.)  ! North relaxation latitude.
      clm_so   = inp_flt(c_str(Clim_slat), -90.)  ! South relaxation latitude.

      imix  = inp_int(c_str(Mixing),   0)
      if (imix.eq.8) then    ! kpp mixing
         kpp_type = inp_int(c_str(KPP_type), 0)
         if (kpp_type.eq.0) print*,'Using windspeed version of KPP mixing'
         if (kpp_type.eq.1) print*,'Using windstress version of KPP mixing'
      endif
         
      limp = inp_int(c_str(Mix_step),  3*ncyc)
      cm_mix = inp_flt(c_str(Mix_cm), 1.25)  ! "m" turbulence coefficient
      cn_mix = inp_flt(c_str(Mix_cn), 0.17)  ! "n" turbulence coefficient
      cn_mix = 1. - cn_mix
      hmin_mix = inp_flt(c_str(Mix_hmin), 10.)  ! H1 - min
      hmax_mix = inp_flt(c_str(Mix_hmax), 100.) ! H1 - max

      ric1_mix = inp_flt(c_str(Mix_ric1), .65) ! Ri for k = 1
      ric2_mix = inp_flt(c_str(Mix_ric2), .25) ! Ri for k = 2:NZ


      iuse_gam = inp_int(c_str(Mix_usegam), 0)   ! 1 - use gammas in jpmix 

      gam1_mix = inp_flt(c_str(Mix_gam1), 1.)  ! Gamma for k = 1
      gam2_mix = inp_flt(c_str(Mix_gam2), 1.)  ! Gamma for k = 2:NZ
      
      iwnd_mix = inp_int(c_str(Mix_wnd),  0)   ! 0:use tau; 1: use windspeeds

      mix_wtop = inp_int(c_str(Mix_wtop),   2)     

      if (imix .eq. 4.and.nsig .eq. 0) then
         print*, 'must set Nsigma>1 to use variable depth mixed layer'
         stop
      endif

      ntrac = inp_int(c_str(Tracers),  0)
      iice  = inp_int(c_str(Ice), 0)
      mice  = inp_int(c_str(Ice_dynamics), 0)
      if (iice.eq.0) mice = 0

      use_trac = (ntrac .ne. 0)
      use_ice  = (iice .ne. 0)
      use_dyice_old = (mice .eq. 1) ! crude ice dynamics
      use_dyice = (mice .eq. 2)     ! Bruno's ice dynamics
      ibio_type  = inp_int(c_str(Bio_type),0) 
      use_bio  = ibio_type .ne. 0
      use_bio_old  = ibio_type .eq. 2
      use_river = (inp_int(c_str(Rivers),  0) .eq. 1)
      coef_river = inp_flt(c_str(River_coef),1.)
      dmin_river = inp_flt(c_str(River_dmin),5.)

      salt_fixed = inp_flt(c_str(Salt_fixed),36.)

      ntrac_bio = 0
      if (use_bio) ntrac_bio  = inp_int(c_str(Bio_ntrac),NTRBIOdef)
      ntrac_tot  = ntrac + ntrac_bio

      ibio_key   = inp_int(c_str(Bio_err), 0)

      if (ibio_key.ne.0)
     *        n_bio = inp_str(c_str(Bio_errfile), 
     *                        finp(1:mlen(finp))//'.bio\0', f_bio)
      nstep_bio = max(1,nint(stpd * inp_days(c_str(Bio_step), delt)))
      nstep_bio_print = nint(stpd * inp_days(c_str(Bio_print), 1.))
      nstep_propag  = nint(stpd * inp_days(c_str(Bio_step_propag), 1.))

      if (ntrac_tot .gt. ntrac_max) then
         print*, 'need to increase NTRAC_MAX in comm_new.h'
         stop
      endif

      if (use_bio) call bio_init

      do i = 1, ntrac
         tr_in(i) = 100.
      enddo
      if (use_bio) then
         inittr(ntrac+ndet) = 2
         inittr(ntrac+nzoo) = 2
         tr_in(ntrac+ndet) = 0.1
         tr_in(ntrac+nzoo) = 0.0143

         inittr(ntrac+nnut) = 4
         inittr(ntrac+nphy) = 4
      endif

      nret = inp_iarr (c_str(Trac_init), ntrac_tot, inittr,  inittr)
      nret = inp_rarr (c_str(Trac_init_value), ntrac_tot, tr_in,  tr_in)

      if (irest.gt.0) then
         do i = 1, ntrac_tot
            inittr(i) = 0
         enddo
      endif


      if (use_bio.or.use_trac) then
         do i = 1, ntrac
            ftrnm_def(i) = "TR"
         enddo
         ftrnm_def(ntrac+ndet) = "DET"
         ftrnm_def(ntrac+nzoo) = "ZOO"
         ftrnm_def(ntrac+nphy) = "PHY"
         ftrnm_def(ntrac+nnut) = "NUT"
      endif


      num1 = inp_sarr(c_str(Tr_name),ntrac_tot,ftrnm_def,80,name_tr,ftrnm)
      num2 = inp_sarr(c_str(Tr_clim_file),0,fbtr,80,n_tr,fbtr)
      if_tr_adv   = inp_int(c_str(Tr_adv_type),1)

      if (if_tr_adv.eq.0) print*,'Tracer advection turned off'

      if (use_bio) then
         itra = ntrac + nnut
      n_tr(itra)=inp_str(c_str(Nitrate_file),c_str(nitrate_data),fbtr(itra))
      endif

c      write(*,*) " "
c      write(*,*) " ftrnm= ", ftrnjm
c      write(*,*) " fbtr= ", fbtr
c      write(*,*) " tr_atm= ", tr_atm
c      write(*,*) " number1= ", num1
c      write(*,*) " number2= ", num2

      isod  = inp_int(c_str(Iso_diffusion), 0)
      use_trdiff  = (isod .ne. 0)
      isod  = inp_int(c_str(Diffiso_scl), 0)
      use_trdiff  = use_trdiff.or.(isod .ne. 0)
      isod  = inp_int(c_str(Diffiso_vel), 0)
      use_modiff  = (isod .ne. 0)
      use_diffiso  = use_trdiff.or.use_modiff
      if (use_modiff.and.nsig.ne.0) then
         print*, 'cannot use non-shapiro diffusion options with sigma layers'
      endif

      iv_top  = inp_int(c_str(Vert_Top),  1)
      if (iv_top.ne.1.and.nsig.eq.0) then
         print*, 'must use Sigma layers for non-constant depth mixed layer'
      endif

      iv_bot  = inp_int(c_str(Vert_Bot),  99) !! 1-CNST;2-Sigm;3-W=0;4-No Motion 
      if (iv_bot.ne.99) then
         print*, 'Vert_Bot ignored, use Nsigma to determine last layer type'
      endif

      iv_bump = inp_int(c_str(Vert_Bump), 0) !! 0-Sigm;2-Sm.bmp;3-Lg.bmp;4-Exp. 
      if (iv_bump.ne.0) then
         print*, 'no interface bumps allowed, reset Vert_Bump = 0'
         stop
      endif
      
      iv_fix  = inp_int(c_str(Vert_Fix),  1) !! 0-Free Surface; 1-Fixed Depth
      if (iv_fix.ne.1) then
         print*, 'reduced gravity runs disabled, reset Vert_Fix = 1'
         stop
      endif

      iv_sys  = inp_int(c_str(Vert_Sys),  99) !! disabled
      if (iv_sys.ne.99) then
         print*, 'Vert_Sys ignored, H is filtered only in sigma layers'
      endif

      mbot_bc = inp_int(c_str(BC_bot_temp), 1) !! 0-constant temp; 1 - d/dz=0
      if (mbot_bc.ne.1) then
         print*, 'do you REALLY want constant bottom temp with bathymetry?'
         print*, 'reconsider using BC_bot_temp = 1'
         stop
      endif


      nordu  = inp_int(c_str(Shap_vel_order), 4)
      mshapu = inp_int(c_str(Shap_vel_type),  3)
      if (mshapu.lt.0) then
         print*,'Shap_vel_type reset to 0'
         mshapu = 0
      endif
      nshapu = inp_int(c_str(Shap_vel_step),  3*ncyc)
      dshapu = inp_flt(c_str(Shap_vel_damp), 1.0)
      dshapv = inp_flt(c_str(Shap_vel_dampv), dshapu)

      nordh  = inp_int(c_str(Shap_scl_order), 4)
      mshaph = inp_int(c_str(Shap_scl_type),  0)
      nshaph = inp_int(c_str(Shap_scl_step),  3*ncyc)
      dshaph = inp_flt(c_str(Shap_scl_damp), 1.0)

      nordp  = inp_int(c_str(Shap_psi_order), 4)
      mshapp = inp_int(c_str(Shap_psi_type),  0)
      nshapp = inp_int(c_str(Shap_psi_step),  1)
      dshapp = inp_flt(c_str(Shap_psi_damp), 1.0)

      nord  = max(nordu,nordh,nordp,nordd)
      MAXFO = nord

      mtau  = inp_int(c_str(Wind_forc),  0)
      if (mtau .eq. 5) itau_cos = inp_int(c_str(Wind_cos), 0)
      flt(1) = 1.
      flt(2) = 1.
      nret =  inp_rarr (c_str(Wind_tauxy), 2, flt, flt)
      tausc = flt(1)
      atau  = flt(2)

      initt  = inp_int(c_str(Temp_init), 0)
      temp_coef  = inp_flt(c_str(Temp_coef), 0.)
      inits  = inp_int(c_str(Salt_init), 0)
      if (inits.ne.0.and.inits.ne.3) then
         print*,'must use Salt_init = 0 or 3'
         stop
      endif

      froude = inp_flt(c_str(Temp_froude), 0.01)

      initq  = inp_int(c_str(Heat_forc), 0)
      initep = inp_int(c_str(EP_forc), 0)
 
      if (initep.eq.8) then
         coef_precip   = inp_flt(c_str(EP_coef_precip),  1.)
         v_sponge      = inp_flt(c_str(EP_sponge),  0.)
         use_budget    = inp_int(c_str(EP_budget),  0) .ne. 0
         tscl_bud      = inp_days (c_str(EP_bud_time), 0) 
         if (tscl_bud.ne.0) tscl_bud = ncyc*delt/tscl_bud
         print*,'budget_coef (tscl_bud),coef ',tscl_bud,coef_precip
      endif

      iep_key   = inp_int(c_str(EP_err),    0)
      if (iep_key .ne. 0)
     *     n_ep      = inp_str(c_str(EP_errfile), 
     *     finp(1:mlen(finp))//'.ep\0', f_ep)

      if (iep_key.gt.0) then
         if (irest .eq. 1. .or. irest .eq. 2) then
            open (unit = IEP_OUT, file = f_ep(1:n_ep), access='APPEND')
         else
            open (unit = IEP_OUT, file = f_ep(1:n_ep))
            write(IEP_OUT,*)'                              Salt Budget'
            write(IEP_OUT,*)
            write(IEP_OUT,*)' coef_pr        evap-brne    precip      sponge      river       relax       total'
            write(IEP_OUT,*)
         endif
      endif

      igas  = inp_int(c_str(Gas_exchange), 0)
      use_wnsp  = (igas .ne. 0)
      if (use_ice.or.initq.eq.8) use_wnsp = .true.

      if (use_ice.and.initq.ne.8) then
         print*,'Using AML-ICE, so I am ignoring your Heat_forc parameter'
      endif

      if (initq.ne.8 .and. iice.ne.1)iwnd_mix = 0 ! valid only if PBL is "ON"

      qcon     = inp_flt (c_str(Rho_CP), 4.12e6)
      tdef = 30.
      if (use_ice) tdef = 0.
      rlx_time_sst = 86400.*inp_days (c_str(Rlx_time_sst), tdef) 
      rlx_time_sss = 86400.*inp_days (c_str(Rlx_time_sss), 30.) 

      itank = inp_int(c_str(Tank_type), 0)
      if (itank.ne.0) cint_fac = inp_int(c_str(Tank_cint_fac), 0.)
      if (itank.ne.0) bint_fac = inp_int(c_str(Tank_bint_fac), 0.)

      if (initq .ge. 8 .or. use_ice) then
         nstep_pbl  = nint(stpd * inp_days(c_str(PBL_step), 5.))
         ipbl_advec = inp_int(c_str(PBL_advec), 1)
         pbl_wmin   = inp_flt(c_str(PBL_wmin),  4.)
         pbl_pnu    = inp_flt(c_str(PBL_pnu),   0.4e+7)
         pbl_delta  = inp_flt(c_str(PBL_delta), 0.25 )
         pbl_pml    = inp_flt(c_str(PBL_pml),   6000.)
         pbl_depth  = inp_flt(c_str(PBL_depth), 600. )
         pbl_betav  = inp_flt(c_str(PBL_betav), 0.17 )
         pbl_grad   = inp_flt(c_str(PBL_grad),    -2.)
         pbl_south  = inp_flt(c_str(PBL_south), -200.)
         pbl_north  = inp_flt(c_str(PBL_north),  200.)
         pbl_eddy   = inp_flt(c_str(PBL_eddy), pbl_eddydef)
         pbl_eddy_nlat  = inp_flt(c_str(PBL_eddy_nlat),pbl_eddy_nlat_def)
         pbl_eddy_nwidth= inp_flt(c_str(PBL_eddy_nwidth),pbl_eddy_nwidth_def)
         pbl_eddy_slat  = inp_flt(c_str(PBL_eddy_slat),pbl_eddy_slat_def)
         pbl_eddy_swidth= inp_flt(c_str(PBL_eddy_swidth),pbl_eddy_swidth_def)
         pbl_eddy_xbox= inp_flt(c_str(PBL_eddy_xbox),pbl_eddy_xbox_def)
         pbl_eddy_ybox= inp_flt(c_str(PBL_eddy_ybox),pbl_eddy_ybox_def)
         c0_fac   = inp_flt(c_str(PBL_c0_fac), c0_facdef)
      endif

      mshapd = 0

      initkmap = inp_int(c_str(Diffiso_map), 0) 
      if (use_diffiso) then
c             default values are given diffiso.h
         n_dif = inp_str(c_str(Dif_file),  c_str(dif_data), fbdif)
         diffiso_alpha = inp_flt(c_str(Diffiso_alpha), diffiso_alphadef) 
         if (diffiso_alpha.ne.0) then
            nordd  = inp_int(c_str(Shap_dens_order), 4)
            mshapd = inp_int(c_str(Shap_dens_type),  -1)
            nshapd = inp_int(c_str(Shap_dens_step),  1)
            dshapd = inp_flt(c_str(Shap_dens_damp), 1.0)
         endif

         diffiso_order = inp_int(c_str(Diffiso_order), 1) 
         nslope  = inp_int(c_str(Diffiso_step), 1) 
         nkmap  = 4*nint(stpd * inp_days(c_str(Kmap_step), 30.)/4.)
         diffiso_eps   = inp_flt(c_str(Diffiso_eps), diffiso_epsdef) 
         diffiso_coef  = inp_flt(c_str(Diffiso_coef), diffiso_coefdef) 
         diff_coef_tr  = diffiso_coef
         diffiso_coef  = inp_flt(c_str(Diffiso_coef_mo), diff_coef_tr) 
         diff_coef_mo  = diffiso_coef
         diff_dz_mo    = inp_flt(c_str(Diffiso_dz_mo), 10000.) 
         diffiso_slmax = inp_flt(c_str(Diffiso_slmax), diffiso_slmaxdef) 
         slred=inp_flt(c_str(Diffiso_slred),diffiso_slreddef) 
         sigzmin=inp_flt(c_str(Diffiso_sigzmin),sigzmindef) 
         if (sigzmin.ge.0) then
            print*,'Diffiso_sigzmin should be negative'
            stop
         endif
         delta_rho = inp_flt(c_str(Diffiso_delta_rho), delta_rhodef) 

         f0 = 0.
         tscl = inp_days (c_str(Diffiso_time), diffiso_tscldef)
         if (tscl.gt.delt/4.) then
            tscl = (delt/tscl)/4.
         else
            tscl = 1.
         endif

         diffiso_cadv=inp_flt(c_str(Diffiso_cadv),diff_coef_tr) 
c           use_diff_cadv = diffiso_alpha.ne.0 .and. diffiso_cadv.ne.0
         use_diff_cadv = diffiso_alpha.ne.0
         if (use_diff_cadv) then
            facz_cnst = inp_flt(c_str(Diffiso_facz_cnst), facz_cnstdef) 
            do k = 1, nz
               facz(k) = facz_cnst
            enddo
            nret = inp_rarr(c_str(Diffiso_facz),nz,facz,facz) 
            psi_rel=inp_flt(c_str(Diffiso_psi_relax),psi_relaxdef) 
            print*,'using Gent-McWilliams mixing ...'
         endif
         if (use_trdiff.ne.0) print*,'diff_coef_tr = ',diff_coef_tr
         if (use_modiff.ne.0) print*,'diff_coef_mo = ',diff_coef_mo
      endif

      
      diff_coef_dd = inp_flt(c_str(Div_damp_coef),0.) 
      isteer = inp_flt(c_str(Div_steering),0) 
      use_div_damp = diff_coef_dd.ne.0.

      use_tvd_vert = inp_def(c_str(TVD_vert))
      use_tvd_horz = inp_def(c_str(TVD_horz))

      cnst_upwind_z = inp_flt(c_str(Upwind_cnst_z),cnst_upwinddef) 
      cnst_upwind_z_ts = inp_flt(c_str(Upwind_cnst_ts),cnst_upwind_z) 
      cnst_upwind_z_tr = inp_flt(c_str(Upwind_cnst_tr),cnst_upwind_z) 
      cnst_upwind_xy = inp_flt(c_str(Upwind_cnst_xy),cnst_upwinddef) 
      cnst_upwind_xy_ts = inp_flt(c_str(Upwind_cnst_xy_ts),cnst_upwind_xy) 
      cnst_upwind_xy_tr = inp_flt(c_str(Upwind_cnst_xy_tr),cnst_upwind_xy) 
      if (cnst_upwind_z.lt.1. .or. cnst_upwind_xy.lt.1.or.
     *    cnst_upwind_z_tr.lt.1. .or. cnst_upwind_xy_tr.lt.1.or.
     *    cnst_upwind_z_ts.lt.1. .or. cnst_upwind_xy_ts.lt.1) then
         print*,'do not use Upwind_cnst less than one'
         stop
      endif
      cupi    = 1./(cnst_upwind_z + 1.)
      cupi_ts = 1./(cnst_upwind_z_ts + 1.)
      cupi_tr = 1./(cnst_upwind_z_tr + 1.)
      hupi    = 1./(cnst_upwind_xy + 1.)
      hupi_ts = 1./(cnst_upwind_xy_ts + 1.)
      hupi_tr = 1./(cnst_upwind_xy_tr + 1.)

      initb = inp_int(c_str(Bathymetry), 0)
      if (initb.ge.3) then
         i_ridge_min = inp_int(c_str(Ridge_min),0)
         i_ridge_max = inp_int(c_str(Ridge_max),nxp+1)
      endif

      initg = inp_int(c_str(Bath_grad), 2)
      grad_fac = inp_flt(c_str(Bath_grad_fac), 1.)
      
      initbt = inp_int(c_str(Bath_type), 0)
      depth_fac = inp_flt(c_str(Bath_fac), .2)

      nzz = inp_rarr(c_str(Z_profile), nz, zin, zin)
      nzh = inp_rarr(c_str(H_profile), nz, hin, hin)

      if (nz .gt. nzz .and. nz .gt. nzh) then
         call perror1('Z_profile or H_profile: not enough terms...',-1)
         stop
      endif

      zin(1) = hin(1)/2.
      z_bot  = hin(1)
      do k = 2, nz
         zin(k) = - zin(k-1) + 2.*z_bot
         z_bot = z_bot + hin(k)
      enddo
      zin(nz+1) = z_bot
      dzin(1) = zin(1)
      do j = 2, nz
         dzin(j) = (zin(j)-zin(j-1))/2.
      enddo
      dzin(nz+1) = zin(nz+1)-zin(nz)

c  Z_profile specification overrides H_profile 

      zdum = inp_flt(c_str(Z_bot), 0)

      if (zdum.eq.0.and.z_bot.eq.0) then
         print*,'must specifiy Z_bot as well as Z_profile'
         stop
      elseif (zdum.gt.0) then
c        reread  Z_profile
         ndum =  inp_rarr(c_str(Z_profile), nz, zin, zin)
         
         zin(1) = zin(2)/3.     ! enforce this condition
         zin(nz+1) = zdum
         
         dzin(1) = zin(1)
         do j = 2, nz
            dzin(j) = (zin(j)-zin(j-1))/2.
            if (dzin(j).lt.0.)then
               print*,'negative layer depths, check Z_profile'
               stop
            endif
         enddo
         dzin(nz+1) = zin(nz+1)-zin(nz)
         if (dzin(nz+1).lt.0.)then
            print*,'negative depth of last layer, check Z_bot'
            stop
         endif
         
         hin(1) = 2.*dzin(1)
         do j = 1, nz
            hin(j) = dzin(j) + dzin(j+1)
         enddo
      endif

      do j = 1, nz
         if (hin(j).le.0) then
            print*,'a layer depth is negative, check Z_profile or H_profile'
            stop
         endif
      enddo

      z_begin = zin(nsig)

      sigma(1) = -1./3.
      sigma(2) = sigma(1)
      do j = 3, nsig
         sigma(j) = dzin(j)/(z_begin-3.*dzin(1))
      enddo
      do j = nsig + 1, nz
         sigma(j) = 0.
      enddo

      if (irest .eq. 0) then
         write(iout,*)'Z_Profile=',(zin(j),j=1,nz)
         write(iout,*)'H_Profile=',(hin(j),j=1,nz)
         write(iout,*)'W_Profile=',(zin(j)+dzin(j+1),j=1,nz)
      endif

      call flush (iout)

      dep = 0.
      do j = 1, nz
         dep = dep + hin(j)
         if (dep.le.tr_zmin) then
            ttr_fac(j)=  1.
         elseif (dep.ge.tr_zmax) then
            ttr_fac(j)=  tr_fac
         else
            ttr_fac(j)=  1.+ (tr_fac - 1.)* (dep- tr_zmin) /(tr_zmax-tr_zmin)
         endif
         if (ttr_fac(j).ne.1.) print*,j,dep,ttr_fac(j)
      enddo

      trans_coef_sst = 0.
      trans_coef_sss = 0.
      
      if (rlx_time_sst.ne.0) trans_coef_sst = (zin(1)+zin(2))/2./rlx_time_sst
      if (rlx_time_sss.ne.0) trans_coef_sss = (zin(1)+zin(2))/2./rlx_time_sss

      print*,'Haney surface coefs, (T,S)',trans_coef_sst,trans_coef_sss

      t_in(nz+1) = inp_flt(c_str(Temp_bot), 0.)
      burger     = inp_flt(c_str(Temp_burger), 0.)

      if (burger.ne.0.) then
         del_rho = - burger**2 * 0.14158
         print*,'temperature profile:'
         do k = 1, nz
            z = zin(k)
            t_in(k) = t_in(nz+1) + 3.81847*(z-4500.)/4500. * del_rho
            print*,k,zin(k),t_in(k)
         enddo
      else
         i = inp_rarr(c_str(T_profile), nz, t_in, t_in)
         if  (initt.le.1.and.i .lt. nz) then
            call perror1('T_profile: not enough terms...',-1)
         endif
      endif
      
      TATM = inp_flt(c_str(Temp_atm), 30.)

      i = inp_rarr(c_str(S_profile), nz, s_in, s_in)
      if  (inits.le.1.and.i .lt. nz) then
         call perror1('S_profile: not enough terms...',-1)
      endif

      s_in(nz+1) = inp_flt(c_str(Salt_bot), 36.)

      SATM = inp_flt(c_str(Salt_atm), 35.4)

      TEMP_BOT = t_in(nz+1)
      SALT_BOT = s_in(nz+1)
      POTND_BOT=inp_flt(c_str(Pdens_bot),pdens_pnt(TEMP_BOT,SALT_BOT))
      SITUD_BOT=inp_flt(c_str(Sdens_bot),sdens_pnt(TEMP_BOT,SALT_BOT,zin(nz+1)))

      rk_bot = inp_flt (c_str(Vert_kbot), rk_botdef)

c.....Bi & Ci are now *real* "Nu" & "Ka" (not scaled by depth):      
      bi = inp_flt (c_str(Bint), 1.e-3)
      ci = inp_flt (c_str(Cint), 1.e-4)
      do i = 1, nz
         bint(i) = bi
         cint(i) = ci
      enddo

      nret=inp_str(c_str(Grid_label),'Data: '//finp(1:mlen(finp))//'\0',label)
      call mem_alloc (p_mask, NXP*NYP*NZ, 1, 'mask')

      if (n_map  .ne. 0) call inp_file(fbmap(1:n_map)//'\0')
      call read_mask ('Grid_mask', nxp, nyp, nxyc, mask)

      call mem_alloc (p_iox, nxyc, 1, 'iox')
      call mem_alloc (p_inx, nxyc, 1, 'inx')
      call mem_alloc (p_iny, nxyc, 1, 'iny')

      npt  = nxyc
      npt1 = 1
      npt2 = 1 + npt 
      npt3 = 1 + 2*npt 
      npt4 = 1 + 3*npt 

      first_step = .true.

      call hfx_pert_init

      call wind_stir_init

      call  nao_init

      ntrac_sur = 0

      if (use_ice) then
         ntrac_ice  = inp_int(c_str(Ice_ntrac),3)
         call ice_init(nx,ny,ntrac_ice,npt)
         ntrac_sur = ntrac_sur + ntrac_ice
      endif

      if (lev_err .ge. 1) then
         write (iout, *) 'NX =', nxp
         write (iout, *) 'NY =', nyp
         write (iout, *) 'NZ =', nz
         write (iout, *) 'NPACK =', nxyc
         write (iout, *) 'NTRACERS_KEITH =', ntrac
         write (iout, *) 'NTRACERS_BIO =', ntrac_bio
         if (ibaro .eq. 0) then
            write (iout, *) 'BARO OFF'
         else
            write (iout, *) 'BARO ON'
         endif
         call flush(iout)
      endif

      return
      end

      subroutine read_mask (tag, nx, ny, npack, mask)
c-----------------------------------------------------
      character*(*) tag
      dimension mask(nx, 1)
      character*1  ch, buff(1000)
      character*465 number
      logical inp_def

      equivalence (ch, buff(1)), (number, buff(2))

      if ( inp_def(tag//'\0') ) then
         ix = 1
         jy = ny
         ic0 = ichar('0')
         do while ( inp_wnxt(buff) .gt. 0)
            if     (ch.eq.'0' .or. ch.eq.'1' .or. ch.eq.'2') then
               do i = 1, nx
                  mask(i,jy) = ichar(buff(i)) - ic0
               enddo
               jy = jy - 1
            elseif (ch.eq.'w'.or.ch.eq.'z'.or.ch.eq.'x'.or.ch.eq.'s'
     *               .or.ch.eq.'b'.or.ch.eq.'t') then
               read (number, *) kx 
               ix = ix + kx
               if (ix - 1 .gt. NX) goto 200
               if (ch.eq.'w') ich = 1 
               if (ch.eq.'x') ich = 2 
               if (ch.eq.'s') ich = 3 
               if (ch.eq.'t') ich = 4 
               if (ch.eq.'b') ich = 5 
               if (ch.eq.'z') ich = 0 
               do i = ix-kx, ix-1
                  mask(i,jy) = ich
               enddo
               if (ix-1 .eq. NX) then
                  ix = 1
                  jy = jy - 1
               endif
            elseif (ch .eq. 'r') then
               read (number, *) ky
               jy = jy - ky + 2
               if (jy .lt. 1) goto 200
               do j = jy+ky-2, jy, -1
                  do i = 1, nx
                     mask(i,j) = mask(i,j+1)
                  enddo
               enddo
               jy = jy - 1
            else
               goto 200
            endif
            
            if (jy .eq. 0) goto 100
         enddo


  100    npack = 0

         do j = 1, ny
            do i = 1, nx
               if (mask(i,j) .ne. 0) npack = npack + 1
            enddo
         enddo

      else
         npack = nx*ny
         do j = 1, ny
            do i = 1, nx
               mask(i,j) = 1
            enddo
         enddo
      endif

      return

  200 write (6, *) '!!!read_mask: wrong mask data, i,j:', ix,jy
      stop
      end
     
c------------------------------------------------------------
      subroutine bio_init
c------------------------------------------------------------
      include 'comm_para.h'
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_bio.h'
      include 'biology.h'

      real inp_flt, inp_days
      logical inp_def
      dimension flt(100)

c      n_bio = inp_str(c_str(Nitrate_file) , c_str(nitrate_data), fbnitrate)

      JPBIO      = inp_int(c_str(Bio_JPBIO),JPBIOdef)
      NFLBIO     = inp_int(c_str(Bio_NFLBIO),NFLBIOdef)

      NDET       = inp_int(c_str(Bio_NDET),NDETdef)
      NZOO       = inp_int(c_str(Bio_NZOO),NZOOdef)
      NPHY       = inp_int(c_str(Bio_NPHY),NPHYdef)
      NNUT       = inp_int(c_str(Bio_NNUT),NNUTdef)

      APMIN  = inp_flt(c_str(Bio_APMIN),APMINdef)
      AZMIN  = inp_flt(c_str(Bio_AZMIN),AZMINdef)
      ANMIN  = inp_flt(c_str(Bio_ANMIN),ANMINdef)
      ADMIN  = inp_flt(c_str(Bio_ADMIN),ADMINdef)
      REDF   = inp_flt(c_str(Bio_REDF),REDFdef)
      SLOPET = inp_flt(c_str(Bio_SLOPET),SLOPETdef)
      TOPTP  = inp_flt(c_str(Bio_TOPTP),TOPTPdef)
      AKNUT  = inp_flt(c_str(Bio_AKNUT),AKNUTdef)
      RCCHL  = inp_flt(c_str(Bio_RCCHL),RCCHLdef)
      RGAMMA = inp_flt(c_str(Bio_RGAMMA),RGAMMAdef)
      TOPTGZ = inp_flt(c_str(Bio_TOPTGZ),TOPTGZdef)
      TMAXGZ = inp_flt(c_str(Bio_TMAXGZ),TMAXGZdef)
      RGZ    = inp_flt(c_str(Bio_RGZ),RGZdef)
      RPPZ   = inp_flt(c_str(Bio_RPPZ),RPPZdef)
      TAUS   = inp_flt(c_str(Bio_TAUS),TAUSdef)
      AKS    = inp_flt(c_str(Bio_AKS),AKSdef)
      FILMAX = inp_flt(c_str(Bio_FILMAX),FILMAXdef)
      RPNAZ  = inp_flt(c_str(Bio_RPNAZ),RPNAZdef)
      RDNAZ  = inp_flt(c_str(Bio_RDNAZ),RDNAZdef)
      EGGZOO = inp_flt(c_str(Bio_EGGZOO),EGGZOOdef)
      TAUZN  = inp_flt(c_str(Bio_TAUZN),TAUZNdef)
      TMMAXP = inp_flt(c_str(Bio_TMMAXP),TMMAXPdef)
      TMMINP = inp_flt(c_str(Bio_TMMINP),TMMINPdef)
      TMMAXZ = inp_flt(c_str(Bio_TMMAXZ),TMMAXZdef)
      TMMINZ = inp_flt(c_str(Bio_TMMINZ),TMMINZdef)
      ANUMIN = inp_flt(c_str(Bio_ANUMIN),ANUMINdef)
      AFDMIN = inp_flt(c_str(Bio_AFDMIN),AFDMINdef)
      TAUDN  = inp_flt(c_str(Bio_TAUDN),TAUDNdef)
      VSED   = inp_flt(c_str(Bio_VSED),VSEDdef)
      AKI    = inp_flt(c_str(Bio_AKI),AKIdef)
      TMAXR  = inp_flt(c_str(Bio_TMAXR),TMAXRdef)
      TMINR  = inp_flt(c_str(Bio_TMINR),TMINRdef)
      ALPH   = inp_flt(c_str(Bio_ALPH),ALPHdef)
      TMUMAX = inp_flt(c_str(Bio_TMUMAX),TMUMAXdef)

c
c       NAMOPT   OPTICAL PARAMETERS 
c
c  XKG0  green absorption coefficient of water
c  XKR0  red absorption coefficent of water
c  XKGP  green absorption coefficient of chl
c  XKRP  red absorption coefficient of chl
c  XLG   green chl exposant for absorption
c  XLR   red chl exposant for absorption
c  RPIG  chlahla+pheo ratio
c  XKD_A  coefficient to take diffused light into account
c  XPAR2PUR_R coeff to transform par to pur in red
c  XPAR2PUR_G coeff to transform par to pur in green
c
      XKG0 = inp_flt(c_str(Bio_XKG0),XKG0def)
      XKR0 = inp_flt(c_str(Bio_XKR0),XKR0def)
      XKGP = inp_flt(c_str(Bio_XKGP),XKGPdef)
      XKRP = inp_flt(c_str(Bio_XKRP),XKRPdef)
      XLG  = inp_flt(c_str(Bio_XLG),XLGdef)
      XLR  = inp_flt(c_str(Bio_XLR),XLRdef)
      RPIG = inp_flt(c_str(Bio_RPIG),RPIGdef)
      XKD_A = inp_flt(c_str(Bio_XKD_A),XKD_Adef)
      XPAR2PUR_R = inp_flt(c_str(Bio_XPAR2PUR_R),XPAR2PUR_Rdef)
      XPAR2PUR_G = inp_flt(c_str(Bio_XPAR2PUR_G),XPAR2PUR_Gdef)
c
c
c  RAU0  density of reference
c
      RAU0   = inp_flt(c_str(Bio_RAU0),RAU0def)
      DELTAR = inp_flt(c_str(Bio_DELTAR),DELTARdef)

      return
      end

c------------------------------------------------------------
      subroutine ice_init(nx,ny,ntrac_ice,npt)
c------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_para.h'
      include 'comm_amlice.h'
      include 'amlice.h'
      include 'dyice.h'
      include 'icedyn.h'

      real inp_flt, inp_days
      integer inp_rarr
      logical inp_def
      dimension flt(100)

      do i = 1, ntrac_ice
         ftrsnm_def(i) = "ICE_TR"
      enddo

      n_cice  = 1
      n_hice  = 2
      n_thice = 3
      ftrsnm_def(n_cice) = "CICE"
      ftrsnm_def(n_hice) = "HICE"
      ftrsnm_def(n_thice)= "THICE"
      i_cice  = npt*(n_cice - 1) + 1
      i_hice  = npt*(n_hice - 1) + 1
      i_thice = npt*(n_thice- 1) + 1

      iret= inp_sarr(c_str(Trs_name),ntrac_ice,ftrsnm_def,80,name_trs,ftrsnm)

c             default values are given amlice.h
      albedoocean = inp_flt(c_str(ICE_albedoocean), albedooceandef)
      albedoice   = inp_flt(c_str(ICE_albedoice  ), albedoicedef) 
      albedof     = inp_flt(c_str(ICE_albedof    ), albedofdef) 
      tfreeze     = inp_flt(c_str(ICE_tfreeze    ), tfreezedef) 
      cicemax     = inp_flt(c_str(ICE_cicemax    ), cicemaxdef)
      hsnow       = inp_flt(c_str(ICE_hsnow      ), hsnowdef)
      sice        = inp_flt(c_str(ICE_sice       ), sicedef)
      itermax     = inp_int(c_str(ICE_itermax    ), itermaxdef)
      ssticemax   = inp_flt(c_str(ICE_ssticemax  ), ssticemaxdef)
      hicemin     = inp_flt(c_str(ICE_hicemin    ), hicemindef)
      tksnow      = inp_flt(c_str(ICE_tksnow     ), tksnowdef)
      tkice       = inp_flt(c_str(ICE_tkice      ), tkicedef)
      tkocean     = inp_flt(c_str(ICE_tkocean    ), tkoceandef)
      hq          = inp_flt(c_str(ICE_hq         ), hqdef)
      hf          = inp_flt(c_str(ICE_hf         ), hfdef)
      if (use_dyice) then

         dyice_lat   = inp_flt(c_str(Dyice_lat),dyice_lat_def)

         Cdair       = inp_flt(c_str(Dyice_Cdair),Cdair_def)
         Cdwater     = inp_flt(c_str(Dyice_Cdwater),Cdwater_def)
         Pstar       = inp_flt(c_str(Dyice_Pstar),Pstar_def)
         C           = inp_flt(c_str(Dyice_C),C_def)
         phi         = inp_flt(c_str(Dyice_phi),phi_def)
         delta       = inp_flt(c_str(Dyice_delta),delta_def)
         etamax      = inp_flt(c_str(Dyice_etamax),etamax_def)
         theta_a     = inp_flt(c_str(Dyice_theta_a),theta_a_def)
         theta_w     = inp_flt(c_str(Dyice_theta_w),theta_w_def)
         do k = 1, ntrac_ice
            Ktracer(k)   = 0.0     ! diff coeff for h (ice thick) [m2/s]
         enddo
         nic  = inp_rarr(c_str(Dyice_Ktracer),0,Ktracer,Ktracer)
         BndyCond = inp_int(c_str(Dyice_BndyCond),1)
         DragLaw  = inp_int(c_str(Dyice_DragLaw),2)
         Rheology = inp_int(c_str(Dyice_Rheology),3)
         InitialGuess = inp_int(c_str(Dyice_InitialGuess),2)
      endif
      if (use_dyice_old) then
         dyice_diff   = inp_flt(c_str(Ice_dyice_diff),dyice_diffdef)
      endif
      
      return
      end

c---------------------------------------------------------
      subroutine baro_init  
c---------------------------------------------------------
      include 'comm_para.h'
      include 'comm_new.h'
      include 'comm_data.h'

      common/run/  nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast
      common/param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20),
     +              itherm,mlc,limp
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch
      common /errors/ ioerr, nstep

c--------------------------------------------------------------
      character*80       f_bar
      common /baro_files/ n_bar, f_bar
      integer iflag(12), RM12_NN 
      real    aflag(12)
      common /y12m_input/ RM12_NN, aflag, iflag 
      common /baro_input/ n_def_cor, mod_scheme, mod_solver, BAR_DEPMIN,
     *                    BAR_DSINK, ibar_key, nbaro, rayl, nonlin_baro

      character*1 BC_I(0:9)
      common /masks/BC_I

      logical use_per_island, use_stan_island, use_sunk
      common /baro_island/
     *                alons_min(10),alons_max(10),alats_min(10),alats_max(10)
     *               ,alon1_min,alon1_max,alat1_min,alat1_max
     *               ,per_lat,use_sunk,use_per_island,use_stan_island
     *               ,itrans,ftrans

      common /new_island/n_islands,in_rhs(0:9),ityp_island(0:9),
     *               alon_island(9),alat_island(9),btrans(0:9),
     *               n_poly(9),alon_poly(100,9),alat_poly(100,9),
     *               i_poly(100,9),j_poly(100,9)
c--------------------------------------------------------------

      logical sink_region
      common /bathymetry/sink_region,sr_depth, alon_sr, blon_sr, alat_sr, blat_sr
      common/friction/ ibfric,b_fric, sw_fric
      integer inp_rarr
      real inp_flt, inp_days, xf(100)
      logical inp_def

      character*10 numbers
      numbers = '1234567890'

      nbaro      = inp_int(c_str(Baro_step),        6)
      n_def_cor  = inp_int(c_str(Baro_defcor_step), 0)
      mod_scheme = inp_int(c_str(Baro_scheme),      9)
      mod_solver = inp_int(c_str(Baro_solver),      1)
      nonlin_baro= inp_int(c_str(Baro_nonlin),      1)
      if (nonlin_baro.ne.1) then
         print*,'Baro_nonlin 0 no longer allowed'
         stop
      endif

      BAR_DSINK  = inp_flt(c_str(Baro_depsink), dep_min - 1.)
      BAR_DEPMIN  = inp_flt(c_str(Baro_depmin),  dep_min)

      rayl       = inp_flt(c_str(Baro_rayl),   0.)
      ibar_key   = inp_int(c_str(Baro_err),    1)
      adelt      = inp_flt(c_str(Island_relax),   0.)
      if (ibar_key .ne. 0)
     *     n_bar      = inp_str(c_str(Baro_errfile), 
     *     finp(1:mlen(finp))//'.bar\0', f_bar)
      
      use_per_island = inp_def(c_str(Antarctica)) .or.
     *     inp_def(c_str(Channel))
      if (use_per_island) then
         ityp_island(0) = 1   ! no sinking of periodic island allowed
         if (iglob.eq.0) then
            print*,'cannot have periodic island in non-periodic domain,' 
            print*,'  - check Periodic and Periodic_island settings'
            stop
         endif
         if (inp_def(c_str(Antarctica))) then
            per_lat = inp_flt(c_str(Antarctica_lat),-60)
            atrans  = inp_flt(c_str(Antarctica_transport),0.)
            if (atrans.ne.0) ityp_island(0) = 3
         elseif (inp_def(c_str(Channel))) then
            per_lat = inp_flt(c_str(Channel_lat),(alat+blat)/2.)
            atrans = inp_flt(c_str(Channel_transport),0.)
            itrans = inp_int(c_str(Channel_transport_type),0)
            ftrans = inp_days(c_str(Channel_transport_freq),1.)
            if (atrans.ne.0) ityp_island(0) = 3
         endif
         btrans(0) = atrans
      endif

      nis  = inp_rarr(c_str(New_island),0,aflag,aflag)
      n_islands = aflag(1)
      if (n_islands.gt.9) then
         print*,'since you are using more than 9 islands'
         print*,'you must allocate more islands in barotropic.h'
         stop
      endif
      if (n_islands.gt.0.and.nis.ne.3*n_islands+1) then
         print*,'you must specify a long,lat and island type for each island'
         stop
      endif

      use_sunk = .false.
      use_stan_island = .false.

      do kis=1,n_islands
         ityp_island(kis) = aflag(3*(kis-1)+2)
         if (ityp_island(kis).eq.2) use_sunk = .true.
         if (ityp_island(kis).eq.1.or.ityp_island(kis).eq.3) 
     *                              use_stan_island = .true.
         alon_island(kis) = aflag(3*(kis-1)+3)
         alat_island(kis) = aflag(3*(kis-1)+4)

         BC_I(kis) = numbers(kis:kis)
         nis  = inp_rarr(n_a(Island_poly,kis),0,xf,xf)
         btrans(kis) = inp_flt(n_a(Island_tran,kis),0.)
         n_poly(kis) = nis/2
         if (nis/2 .ne. float(nis)/2.) then
            print*,'check long,lat points for New_island_poly',kis
            stop
         endif
         do i = 1, n_poly(kis)
            alon_poly(i,kis) = xf(2*(i-1)+1)
            alat_poly(i,kis) = xf(2*i)
         enddo
      enddo

      sink_region = .false.
      nsr  = inp_rarr(c_str(Sink_region),0,aflag,aflag)
      if (nsr.gt.0) then  
         sink_region = .true.
         if (nsr.ne.5) then  
            print*,'For Sink_region, must specify 
     *           long(min,max),lat(min,max),depth'
            stop
         endif
         sr_depth = aflag(5)
         alon_sr = aflag(1)
         blon_sr = aflag(2)
         alat_sr = aflag(3)
         blat_sr = aflag(4)
      endif

c------------FROM ".y12m"--------------------------------------------
c     RM12_NN    : size of workspace in numbers of NONZ
c     iflag(11)  : number of iterations
c     iflag    : see y12m documentation
c     aflag    : see y12m documentation
c     default values:

      do i = 1, 12
         iflag(i) = 0
         aflag(i) = 0
      enddo

      RM12_NN   = inp_int(c_str(Baro_nnonz), 4)
      iflag(11) = inp_int(c_str(Baro_niter), 100)

      iflag(1)  =   1
      iflag(2)  =   3
      iflag(3)  =   1
      iflag(5)  =   2
      aflag(1)  = 2.
      aflag(2)  = 1.e-4
      aflag(3)  = 1.e6
      aflag(4)  = 1.e-12
      aflag(5)  = -1000.

      nret = inp_iarr (c_str(Baro_iflag), 5, iflag,  iflag)
      nret = inp_rarr (c_str(Baro_aflag), 5, aflag,  aflag)

      return
      end
