! source file: /Users/oschlies/UVIC/master/source/embm/solve.F
      subroutine solve (n)

!=======================================================================
!     solve for tracer distribution after diffusion

!     input:
!       n    = tracer number

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

      implicit none

      include "param.h"
      include "solve.h"
      include "atm.h"
      include "cembm.h"
      include "csbc.h"
      include "grdvar.h"
      include "coord.h"
      include "levind.h"
      include "ice.h"

      integer i, ii, ierr, j, jj, k, n

      logical done

      real afw, afe, afn, atc, ate, atn, ats, atnc, atsc, atw, b, dt
      real dtss, dfw, dfe, dfn, fa, fb, fc, fd, fe, ff, fg, fh, tmp, x
      real forc(imt,jmt)

      dtss = dts

!-----------------------------------------------------------------------
!     set the forcing for each tracer
!-----------------------------------------------------------------------

      if (n .eq. isat) then

!       temperature

        fa = dtss/(cpatm*rhoatm*sht)
        fb = dtss*vlocn/(cpatm*rhoatm*sht)
        fc = dtss*slice/(cpatm*rhoatm*sht) - fb
        fd = dtss*flice/(cpatm*rhoatm*sht)
        fe = scatter*(1. + pass)
        do j=2,jmtm1
          do i=2,imtm1
            forc(i,j) = fa*(solins(i,j)*sbc(i,j,iaca)*fe
     &                - dnswr(i,j)*scatter - outlwr(i,j)
     &                + uplwr(i,j) + upsens(i,j))
!           latent heat from total precipitation as water
            forc(i,j) = forc(i,j) + precip(i,j)*fb
            if (psno(i,j) .gt. 0.0) then
!             correct for latent heat from snow
              forc(i,j) = forc(i,j) + fc*psno(i,j)
            elseif (land_map(i,j) .eq. 0) then
!             minus heat to melt snow over land if land stores no heat
              forc(i,j) = forc(i,j) + fd*psno(i,j)
            endif
          enddo
        enddo

      elseif (n .eq. ishum) then

!       humidity

        fa = dtss/(rhoatm*shq)
        do j=2,jmtm1
          do i=2,imtm1
            forc(i,j) = fa*evap(i,j)
          enddo
        enddo

      else

!       other tracers

        do j=2,jmtm1
          do i=2,imtm1
            forc(i,j) = c0
          enddo
        enddo

      endif
      call embmbc (forc)

!-----------------------------------------------------------------------
!     calculate new coefficients if required
!-----------------------------------------------------------------------

      if (newcoef(lf,n)) call coef (n)

!-----------------------------------------------------------------------
!     shuffle in time
!-----------------------------------------------------------------------

      do j=1,jmt
        do i=1,imt
          tmp = at(i,j,2,n)
          at(i,j,2,n) = at(i,j,lf,n) + forc(i,j)
          at(i,j,1,n) = tmp
        enddo
      enddo

!-----------------------------------------------------------------------
!     load rhs into the solver array
!-----------------------------------------------------------------------

      k = 0
      do j=2,jmtm1
        do i=2,imtm1
          k = k + 1
          bv(k) = at(i,j,2,n)
          xv(k) = at(i,j,1,n)
        enddo
      enddo

!-----------------------------------------------------------------------
!     solve for tracer
!-----------------------------------------------------------------------

      call mgrid (xv, ap(1,lf,n), an(1,lf,n), as(1,lf,n), ae(1,lf,n)
     &,           aw(1,lf,n), bv, 1, iimtm2, 1, jjmtm2, iimtm2, jjmtm2
     &,           itin(n), levelin, epsin(n), itout(n), levelout
     &,           epsout(n))

      newcoef(lf,n) = .false.
      if (epsout(n) .gt. epsin(n)) write(*,*)
     &  '==> Warning:  atmospheric solver not converging in ',
     &  itout(n),' iterations ( eps = ',epsout(n), ' > ',epsin(n),' )'

!-----------------------------------------------------------------------
!     copy new solution from left hand side
!-----------------------------------------------------------------------

      k = 0
      do j=2,jmtm1
        do i=2,imtm1
          k = k + 1
          at(i,j,2,n) = xv(k)
        enddo
      enddo

!-----------------------------------------------------------------------
!     set boundary conditions
!-----------------------------------------------------------------------

      call embmbc (at(1,1,2,n))

      return
      end

      subroutine coef (n)

!=======================================================================
!     compute matrix coefficients

!     input:
!       n    = tracer number

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

      implicit none

      include "param.h"
      include "solve.h"
      include "grdvar.h"
      include "cembm.h"
      include "atm.h"
      include "csbc.h"

      integer i, ide, ii, ielm, iord, j, jdn, jj, jord, n, iwx, iwy

      real acej, acnj, adde, addn, adds, addw, cc, ce, cew, cmax, cn
      real cns, cs, cw, dcej, dcnj, fe, fn, fs, fw, ue, un, uw, ve
      real vn, vs, vw

      ielm = 0
      iord = 0

      iwx = 0
      iwy = 0
      if (n .eq. 1) then
        iwx = iwxt
        iwy = iwyt
      endif
      if (n .eq. 2) then
        iwx = iwxq
        iwy = iwyq
      endif

      if (iwx .gt. 0 .and. iwy .gt. 0) then
        call embmbc (sbc(1,1,iwx))
        call embmbc (sbc(1,1,iwy))
      endif

      cmax = 3.9e10

      do j=2,jmtm1
        jj = j - 1
        jdn = j
        do i=2,imtm1
          ii = i - 1
          ide = i

!-----------------------------------------------------------------------
!         set coefficients for implicit diffusion
!-----------------------------------------------------------------------

          cs = dn(i,j-1,n)
          cn = dn(i,jdn,n)
          cw = de(i-1,j,n)
          ce = de(ide,j,n)

!-----------------------------------------------------------------------
!         closed north/south boundary conditions for diffusion
!-----------------------------------------------------------------------

          if (j .eq. 2) cs = c0
          if (j .eq. jmtm1) cn = c0

          cs =-dts*cs*dsgrd(j)
          cn =-dts*cn*dngrd(j)
          cw =-dts*cw*cstr(j)*cstr(j)*dwgrd(i)
          ce =-dts*ce*cstr(j)*cstr(j)*degrd(i)

          cc = 1.0 - cs - cn - cw - ce

          if (n .eq. 1 .or. n .eq. 2) then

!-----------------------------------------------------------------------
!         set coefficients for up-stream advection
!-----------------------------------------------------------------------

            vs = (sbc(i-1,j-1,iwy) + sbc(i,j-1,iwy))
            vn = (sbc(i-1,j,iwy) + sbc(i,j,iwy))
            uw = (sbc(i-1,j-1,iwx) + sbc(i-1,j,iwx))
            ue = (sbc(i,j-1,iwx) + sbc(i,j,iwx))

            fs = p5*(c1 + sign(c1,vs))
            fn = p5*(c1 + sign(c1,vn))
            fw = p5*(c1 + sign(c1,uw))
            fe = p5*(c1 + sign(c1,ue))

!-----------------------------------------------------------------------
!           closed north/south boundary conditions for advection
!-----------------------------------------------------------------------

            if (j .eq. 2) vs = c0
            if (j .eq. jmtm1) vn = c0

            cs = cs - dts*fs*vs*asgrd(j)
            cn = cn + dts*(c1-fn)*vn*angrd(j)
            cw = cw - dts*fw*uw*cstr(j)*azgrd(i)
            ce = ce + dts*(c1-fe)*ue*cstr(j)*azgrd(i)
            cc = cc + dts*(fn*vn*angrd(j)-(c1-fs)*vs*asgrd(j)
     &         + (fe*ue - (c1-fw)*uw)*cstr(j)*azgrd(i))
          endif

          iord = iord + 1

!-----------------------------------------------------------------------
!         load the coefficients for the multigrid solver
!-----------------------------------------------------------------------

          ap(iord,lf,n) = cc
          an(iord,lf,n) = -cn
          as(iord,lf,n) = -cs
          ae(iord,lf,n) = -ce
          aw(iord,lf,n) = -cw

        enddo
      enddo

      return
      end