! source file: /Users/oschlies/UVIC/master/testcase/updates/setembm.F
      subroutine setembm

!=======================================================================
!     initialize the energy-moisture balance model

!     based on code by: A. Fanning and M. Eby
!=======================================================================

      implicit none

      include "param.h"
      include "calendar.h"
      include "solve.h"
      include "switch.h"
      include "coord.h"
      include "grdvar.h"
      include "cembm.h"
      include "atm.h"
      include "insolation.h"
      include "ice.h"
      include "evp.h"
      include "riv.h"
      include "tmngr.h"
      include "levind.h"
      include "csbc.h"
      include "scalar.h"
      include "veg.h"
      real rveg(imt,jmt)

      character(120) :: fname, new_file_name
      character(3) :: a3

      integer i, ie, ii, iou, is, j, je, jj, js, jz, k, m, n, nsolve
      integer nu, nsum, iice, ib(10), ic(10)

      logical exists, inqvardef

      real dlam, dphi, dtatms, dte, dyz, eccice, grarea, saltmax
      real si, ssh, t1, tair, yz_max, yz_min, wz, calday
      real zrel, C2K, c100
      real dmsk(imt,jmt)

      namelist /tsteps/ dtts, dtuv, dtsf, dtatm, dtatms, namix, segtim
      namelist /embm/   rlapse, rfactor, scatter, rhmax
      namelist /carbon/ co2ccn, co2yri, co2yrf, co2rate, co2for, co2_yr
     &,                 c14_yr, dc14ccn, c14prod
      namelist /paleo/  orbit_yr, pyear, eccen, obliq, mvelp
      namelist /ice/    niats, nivts, dampice, tsno, hsno_max, ice_yr
     &,                  ice_calb, sno_calb
      namelist /veg/    veg_alb, veg_rl, veg_rs, veg_smd, crops_yr
     &,                 icrops, idesert, iice
      namelist /solar/  solarconst, solar_yr

      is = 1
      ie = imt
      js = 1
      je = jmt

      c100 = 100.
      C2K = 273.15

!-----------------------------------------------------------------------
!     define default parameters
!-----------------------------------------------------------------------

!     set defaults namelist embm
      rlapse     = 5.e-5
      rfactor    = 0.2
      scatter    = 0.23
      rhmax      = 0.85

!     set defaults namelist carbon
      co2ccn     = 280
      co2yri     = 1.e20
      co2yrf     = 1.e20
      co2rate    = 0.
      co2for     = 5.35e03
      co2_yr     = 2.e20
      c14_yr     = 2.e20
      dc14ccn    = 0.
      c14prod    = 0.

!     set defaults namelist paleo
      pyear      = 2.e20
      orbit_yr   = 2.e20
!     default user orbit is for 1950
      eccen      = 1.672393E-02
      obliq      = 23.44627
      mvelp      = 102.0391

!     set defaults namelist ice
      niats      = 1
      nivts      = 1
      dampice    = 5.
      tsno       = 0.
      hsno_max   = 1000.
      ice_yr     = 2.e20
      ice_calb   = 0.35
      sno_calb   = 0.25

!     set defaults namelist veg
      veg_alb(1:7) = (/0.17, 0.17, 0.22, 0.22, 0.22, 0.30, 0.60/)
      veg_rl(1:7) = (/2.9, 0.91, 0.11, 0.053, 0.039, 0.039, 0.020/)
      veg_rs(1:7) = (/0.75, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
      veg_smd(1:7) = (/2.0, 3.0, 0.1, 0.1, 0.5, 0.01, 0.01/)
      crops_yr = 2.e20
      icrops = 3
      idesert = 6
      iice = 7

!     set defaults namelist solar
      solarconst = 1.368e6
      solar_yr = 2.e20

!     set other stuff
      cdatm     = 1.e-3
      cpatm     = 1.004e7
      sht       = 8.4e5
      shq       = 1.8e5
      shc       = 7.9e5
      rhoatm    = 1.250e-3
      esatm     = 4.6e-05
      pcfactor  = 0.
      cssh      = 3.8011e-3
      cfc11ccnn = 0.
      cfc11ccns = 0.
      cfc12ccnn = 0.
      cfc12ccns = 0.
      dc14ccnn  = 0.
      dc14ccne  = 0.
      dc14ccns  = 0.

      rhoocn    = 1.035
      esocn     = 5.4e-5
      vlocn     = 2.501e10

      cdice     = 5.5e-3
      rhoice    = 0.913
      rhosno    = 0.330
      esice     = 5.347e-5
      slice     = 2.835e10
      flice     = 3.34e9
      condice   = 2.1656e5

      soilmax   = 15.
      eslnd     = 5.347e-5
      cs_alb    = 0.08

      nivc      = 1
      dtatms    = 1800.
      ns        = 30

      dalt_v    = 3.3e-3
      dalt_o    = 1.4e-3
      dalt_i    = 1.4e-3

!-----------------------------------------------------------------------
!     read namelist files
!-----------------------------------------------------------------------
      call getunit (iou, 'control.in', 'f s r')
      read  (iou, tsteps, end=100)
100   continue
      write (stdout, tsteps)
      call relunit (iou)

      call getunit (iou, 'control.in', 'f s r')
      read  (iou, embm, end=101)
101   continue
      write (stdout,embm)
      call relunit (iou)

      call getunit (iou, 'control.in', 'f s r')
      read  (iou, carbon, end=102)
102   continue
      write (stdout,carbon)
      call relunit (iou)
      co2ccni = co2ccn

      call getunit (iou, 'control.in', 'f s r')
      read  (iou, paleo, end=103)
103   continue
      write (stdout,paleo)
      call relunit (iou)

      call getunit (iou, 'control.in', 'f s r')
      read  (iou, ice, end=104)
104   continue
      write (stdout,ice)
      call relunit (iou)

      call getunit (iou, 'control.in', 'f s r')
      read  (iou, veg, end=105)
105   continue
      write (stdout,veg)
      call relunit (iou)

      call getunit (iou, 'control.in', 'f s r')
      read  (iou, solar, end=106)
106   continue
      write (stdout,solar)
      call relunit (iou)

!     ensure pass is between zero and one.
      pass =  min(max((1. - scatter), 0.), 1.)

!     sort out forcing years
      if (pyear  .gt. 1.e20) pyear = 1850.
      if (co2_yr .gt. 1.e20) co2_yr = pyear
      if (c14_yr .gt. 1.e20) c14_yr = pyear
      if (orbit_yr .gt. 1.e20) orbit_yr = pyear
      if (crops_yr .gt. 1.e20) crops_yr = pyear
      if (ice_yr .gt. 1.e20) ice_yr = pyear
      if (solar_yr .gt. 1.e20) solar_yr = pyear

      if (mod(timavgint, segtim) .gt. 1.e-6 .and. timavgint .gt. 0.)
     &  then
        t1 = nint(timavgint/segtim)*segtim
        if (t1 .lt. segtim) t1 = segtim
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: "timavgint" does not contain an integral number'
     &,   '              of coupling time steps "segtim".              '
        write (stdout,*) '              (changed "timavgint" from '
     &, timavgint,' days to ', t1,' days to insure this condition)'
        timavgint = t1
      endif
      if (timavgint .eq. 0.) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: averaging interval "timavgint" = 0. implies no '
     &,   '             averaging when "time_averages" is enabled      '
      endif
      if (timavgint .gt. timavgper) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: the interval "timavgint" exceeds the averaging '
     &,   '             period "timavgper" for option "time_averages"  '
      endif
      if (timavgint .lt. timavgper) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: averaging period "timavgper" exceeds interval  '
     &,   '             "timavgint". Setting timavgper = timavgint     '
        timavgper = timavgint
      endif
      if (timavgper .eq. 0.) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: the averaging period "timavgper" is zero. The  '
     &,   '             average will be over only one time step!       '
      endif
      write (stdout,'(/,1x,a,f10.2,a,/,1x,a,f10.2,a)')
     &  '==> Time averages will be written every ', timavgint, ' days, '
     &, '    with an averaging period of         ', timavgper, ' days. '

!-----------------------------------------------------------------------
!     calculate the expansion coefficients for Berger's solution for
!     the year of the initial conditions
!-----------------------------------------------------------------------

      call orbit (orbit_yr, eccen, obliq, mvelp, lambm0)
      write (stdout,*)
      write (stdout,*) 'Initial Orbital Parameters:'
      write (stdout,*) '  Orbital Year:', orbit_yr
      write (stdout,*) '  Eccentricity:', eccen
      write (stdout,*) '  Obliquity:   ', obliq
      write (stdout,*) '  Longitude of Perihelion:', mvelp+180.

!-----------------------------------------------------------------------
!     calculate annual average insolation and Coriolis factor
!-----------------------------------------------------------------------

      radian = 360./(2.*pi)
      do j=1,jmt
        do i=1,imt
!         calculate coriolis parameter
          fcor(i,j) = 2.*omega*sin(ulat(i,j)/radian)
        enddo
      enddo

!-----------------------------------------------------------------------
!     read diffusion
!-----------------------------------------------------------------------

      dn(:,:,:) = 1.e10
      de(:,:,:) = 1.e10
      fname = new_file_name ("diff.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
        print*, "Warning => ", trim(fname), " does not exist."
      else
        call openfile (fname, iou)
        ib(:) = 1
        ic(:) = imt
        ic(2) = jmt
        do n=1,nat
          if (n .lt. 1000) write(a3,'(i3)') n
          if (n .lt. 100) write(a3,'(i2)') n
          if (n .lt. 10) write(a3,'(i1)') n
          call getvara ('dn_'//trim(a3), iou, imt*jmt, ib, ic, dn(1,1,n)
     &,     c1, c0)
          call getvara ('de_'//trim(a3), iou, imt*jmt, ib, ic, de(1,1,n)
     &,     c1, c0)
        enddo
        call closefile (iou)
      endif

!-----------------------------------------------------------------------
!     set solver parameters
!-----------------------------------------------------------------------

      nsolve = 0
      itin(1:nat)  = 500        ! max solver iterations
      epsin(1) = 1.e-3
      epsin(2:nat) = 1.e-5
      nsolve = nsolve + 1
      levelin = 20              ! max coarse grid level
      if (nsolve .ne. 1) then
        write(*,*) '==> Error: more or less than one solver defined.'
        write(*,*) '           Use only one of uvic_embm_adi,'
     &,   ' uvic_embm_mgrid, uvic_embm_slap, uvic_embm_essl or'
     &,   ' or uvic_embm_sparskit or uvic_embm_explicit'
        stop '=>setembm'
      endif

!-----------------------------------------------------------------------
!     check latent heats will sum to zero
!-----------------------------------------------------------------------

      if (slice .ne. vlocn + flice) write (stdout,'(/,a)')
     &   '==> Warning: changing latent heat of fusion to conserve heat'
        flice = slice - vlocn

!-----------------------------------------------------------------------
!     calculate grid terms for the atmospheric solver
!-----------------------------------------------------------------------

      do j=2,jmtm1
        dsgrd(j) = csu(j-1)/(dyu(j-1)*cst(j)*dyt(j))
        dngrd(j) = csu(j)/(dyu(j)*cst(j)*dyt(j))
        asgrd(j) = csu(j-1)/(2.*cst(j)*dyt(j))
        angrd(j) = csu(j)/(2.*cst(j)*dyt(j))
      enddo
      do i=2,imtm1
        dwgrd(i) = 1./(dxu(i-1)*dxt(i))
        degrd(i) = 1./(dxu(i)*dxt(i))
        azgrd(i) = 1./(2.*dxt(i))
      enddo

!-----------------------------------------------------------------------
!     set initial conditions or read a restart
!-----------------------------------------------------------------------

      newcoef(:,:) = .true.

      nats = namix
      dayoyr = 1.
      itt = 0
      irstdy = 0
      msrsdy = 0
      at(:,:,:,:) = 0.
      tair = 13.
      at(:,:,:,isat) = tair
      ssh = cssh*exp(17.67*tair/(tair + 243.5))
      rh(:,:) = rhmax
      at(:,:,:,ishum) = rhmax*ssh
cao# if defined atmco2const
cao      at(:,:,:,ico2) = co2ccn
cao# endif
      precip(:,:) = 0.
      aicel(:,:,:) = 0.
      hicel(:,:,:) = 0.
      soilm(:,:,:) = 0.
      surf(:,:) = 0.
      hice(:,:,:) = 0.
      aice(:,:,:) = 0.
      tice(:,:) = 0.
      hsno(:,:,:) = 0.
      uice(:,:) = 0.
      vice(:,:) = 0.
      sbc(:,:,isu) = 0.
      sbc(:,:,isv) = 0.
      sbc(:,:,igu) = 0.
      sbc(:,:,igv) = 0.
      sig11n(:,:) = 0.
      sig11e(:,:) = 0.
      sig11s(:,:) = 0.
      sig11w(:,:) = 0.
      sig22n(:,:) = 0.
      sig22e(:,:) = 0.
      sig22s(:,:) = 0.
      sig22w(:,:) = 0.
      sig12n(:,:) = 0.
      sig12e(:,:) = 0.
      sig12s(:,:) = 0.
      sig12w(:,:) = 0.
      bv(:) = 0.
      xv(:) = 0.

!-----------------------------------------------------------------------
!     set land ice data and tracer grid ocean mask
!-----------------------------------------------------------------------

      call icedata

      if (.not. init) then
        fname = new_file_name ("restart_embm.nc")
        call embm_rest_in (fname, is, ie, js, je)
      endif
      flux(:,:,:) = 0.0

!-----------------------------------------------------------------------
!     read any atmospheric forcing data
!-----------------------------------------------------------------------

      call setatm

!-----------------------------------------------------------------------
!     read land elevations
!-----------------------------------------------------------------------

      elev(:,:) = 0.
      fname = new_file_name ("elev.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
        print*, "Warning => ", trim(fname), " does not exist."
      else
        ib(:) = 1
        ic(:) = imt
        ic(2) = jmt
        call openfile (fname, iou)
        call getvara ('elev', iou, imt*jmt, ib, ic, elev, c100, c0)
        call closefile (iou)
      endif
!     check for negative elevations
      where (elev(:,:) .lt. 0.) elev(:,:) = 0.

!-----------------------------------------------------------------------
!     set velocity grid ocean mask
!-----------------------------------------------------------------------

      umsk(:,:) = 0.
      do j=2,jmtm1
        do i=2,imtm1
          umsk(i,j) = min (tmsk(i,j), tmsk(i+1,j), tmsk(i,j+1)
     &,                    tmsk(i+1,j+1))
        enddo
      enddo
      call embmbc (umsk)
!     remove isolated bays
      do j=2,jmtm1
        do i=2,imtm1
          tmsk(i,j) = max (umsk(i,j), umsk(i-1,j), umsk(i,j-1)
     &,                    umsk(i-1,j-1))
        enddo
      enddo
      call embmbc (tmsk)
      do j=2,jmtm1
        do i=2,imtm1
          umsk(i,j) = min (tmsk(i,j), tmsk(i+1,j), tmsk(i,j+1)
     &,                    tmsk(i+1,j+1))
        enddo
      enddo
      call embmbc (umsk)

      do j=1,jmt
        do i=1,imt
          if (tmsk(i,j) .ge. 0.5) then
            if (hice(i,j,1) .le. 0.) aice(i,j,1) = 0.
            if (hice(i,j,2) .le. 0.) aice(i,j,2) = 0.
          endif
        enddo
      enddo
      call embmbc (aice(1,1,1))
      call embmbc (aice(1,1,2))

!-----------------------------------------------------------------------
!     set the river model
!-----------------------------------------------------------------------

      call rivinit

!-----------------------------------------------------------------------
!     set ocean coalbedo
!-----------------------------------------------------------------------

      do j=1,jmt
        do i=1,imt
          if (kmt(i,j) .gt. 0) then
!           varies from 0.895 at the equator to 0.815 at the pole
            sbc(i,j,isca) = 0.855 + 0.04*cos(abs(tlat(i,j))*2./radian)
          endif
        enddo
      enddo

!-----------------------------------------------------------------------
!     read vegetation class
!-----------------------------------------------------------------------

      rveg(:,:) = 0.
      fname = new_file_name ("veg_type.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
        print*, "Warning => ", trim(fname), " does not exist."
      else
        call openfile (fname, iou)
        call getvara ('veg', iou, imt*jmt, ib, ic, rveg, c1, c0)
        call closefile (iou)
      endif
      do j=1,jmt
        do i=1,imt
          if (kmt(i,j) .gt. 0) then

            iveg(i,j) = iice
          elseif (rveg(i,j) .gt. 0.6 .and. rveg(i,j) .lt. 7.4) then
            iveg(i,j) = nint(rveg(i,j))
          else
            print*, "==> Warning: invalid vegetation type, ", rveg(i,j)
     &,             " at point (i,j), ", i, j
            print*, "    Setting point to vegetation type ice"
            iveg(i,j) = iice
          endif
        enddo
      enddo

      call gvsbc

!----------------------------------------------------------------------
!     initialize elastic viscous plastic variables
!-----------------------------------------------------------------------

      dlam = dxu(int(imt/2))/100.
      dphi = dyu(int(jmt/2))/100.
      diff1 = 0.004
      diff1 = diff1*dlam
      diff2 = diff1*dlam**2
      eccice = 2.
      ecc2 = 1./(eccice**2)
      ecc2m = 2.*(1.-ecc2)
      ecc2p = (1.+ecc2)
      zetamin = 4.e11
      eyc = 0.25
      dte = dtatm/float(ndte)
      dtei = 1./dte
      floor = 1.e-11
      do j=2,jmtm1
        do i=2,imtm1
           xyminevp=(min(cst(j)*dxt(i),dyt(j)))**2
        enddo
      enddo

!-----------------------------------------------------------------------
!     check for even flux averaging
!-----------------------------------------------------------------------

      if (mod(namix,2) .ne. mod(nint(segtim*daylen/dtatm),2)) then
        write(*,*) '==> Error: time steps between mixing and coupling'
        write(*,*) '           in the atmosphere must both be even to'
        write(*,*) '           use uvic_embm_even_fluxes.'
        stop '=>setembm'
      endif

      if (mod(nmix,2) .ne. mod(nint(segtim*daylen/dtocn),2)) then
        write(*,*) '==> Error: time steps between mixing and coupling'
        write(*,*) '           in the ocean must both be even to use'
        write(*,*) '           uvic_embm_even_fluxes.'
        stop '=>setembm'
      endif

      if (mod(namix,2) .ne. mod(nats,2)) then
        write(*,*) '==> Warning: restart was not saved with even flux'
        write(*,*) '             averaging. Starting with a mixing time'
        write(*,*) '             step.'
        nats = namix
      endif

!-----------------------------------------------------------------------
!     check ice velocity calculation
!-----------------------------------------------------------------------

      if (nivts .gt. nint(segtim*daylen/dtatm)) then
        write(*,*) '==> Warning: ice velocities will be calculated'
        write(*,*) '             every coupling time.'
        nivts =  nint(segtim*daylen/dtatm)
      endif
!-----------------------------------------------------------------------
!     zero time average accumulators
!-----------------------------------------------------------------------

      call ta_embm_snap (is, ie, js, je, 0)

!-----------------------------------------------------------------------
!     zero integrated time average accumulators
!-----------------------------------------------------------------------

      call ta_embm_tsi (0)

      return
      end