subroutine assim(iday,imo,iyr,ndmo)
c
c  Blends satellite and model chlorophyll data together.  Satellite
c  represents truth field, although weight functions are definable. 
c
      save ifst,wtmon,iday1,ifle
#include "gloparam.F.h"
#include "ostate.F.h"
#include "definebio.h"
#include "combarr.h"
#include "radpth.h"
      parameter(ndmomx=31)
      real satchl(im,jm)
      real chlmod(im,jm)
      real chls(im,jm)
      real flap(im,jm),flapb(im,jm)
      real cblend(im,jm)
      common /bsatmon2/ nachl(im,jm)
      common /bsat/ chlsat(im,jm,ndmomx)
      common /bsatmon/ achlmon(im,jm),sumonchl(im,jm)
      common /blonlat/ xlonb(im),ylatb(jm)
      common /bicec/ picec(len)
      common /bwt/ wtmod(im,jm)
      real Pmod(nchl)
      data flagl /999.0E9/
      data flagm /0.0/
      data epsilon /0.0001/
      data ifst /0/
      data ifle /0/
c
      if (ifst .eq. 0)then
       ifst = 1
       wtmon = 0.25  !weight of monthly mean
       iday1 = iday
       rlat = alat0 - 0.5*dlat
       do j = 1,jm
        rlat = rlat + dlat
        ylatb(j) = rlat
       enddo
       rlon = along0 - 0.5*dlong
       do i = 1,im
        rlon = rlon + dlong
        xlonb(i) = rlon
        if (xlonb(i) .gt. 180.0)xlonb(i) = xlonb(i) - 360.0
       enddo
       achlmon = 0.0
       sumonchl = 0.0
       nachl = 0.0
       chlsat = 0.0
       open(81,file='satchl.dat',status='old',form='unformatted',err=19)
       ifle = 1
       do id = iday,iday+ndmo-1
        read(81,end=18)id2,iyr2
        if (iyr2 .ne. iyr)then
         close(81)
         write(6,*)'No Assimilation Data for Day, Year ',iday,iyr
         go to 20
        else
         read(81,end=18)satchl
         do j = 1,jm
          do i = 1,im
           chlsat(i,j,id2-iday1+1) = satchl(i,j)
           if (satchl(i,j) .gt. flagm .and. satchl(i,j) .lt. flagl)then
            sumonchl(i,j) = sumonchl(i,j) + satchl(i,j)
            nachl(i,j) = nachl(i,j)+1
           endif
          enddo
         enddo
        endif
       enddo
18     close(81)
       do j = 1,jm
        do i = 1,im
         if (nachl(i,j) .gt. 0)then
          achlmon(i,j) = sumonchl(i,j)/float(nachl(i,j))
         else
          achlmon(i,j) = flagl
         endif
        enddo
       enddo
c   Eliminate outliers
       do j = 1,jm
        do i = 1,im
         do id = 1,ndmo
          if (chlsat(i,j,id) .lt. flagl)then
           if (chlsat(i,j,id) .gt. 2.0*achlmon(i,j))then
            sumonchl(i,j) = sumonchl(i,j) - chlsat(i,j,id)
            nachl(i,j) = nachl(i,j)-1
            chlsat(i,j,id) = flagm
           endif
          endif
         enddo
        enddo
       enddo
       chlmod = 0.0 
       k = 1
       m = indext1
       do n = nnut+1,ntyp-nzoo
        nw = 0
        do j = 1,jm   
         do i = 1,im
          if (Depth(i,j) .gt. 200.0)then
           nw = nw+1
           chlmod(i,j) = chlmod(i,j) + P(nw,k,m,n)
          endif
         enddo    !i
        enddo     !j
       enddo       !n
       call assimwts(imo,chlmod)
      endif
c
      if (ifle .eq. 0)go to 20
c  Satellite data
      id = iday-iday1+1
      rmins = 9999.0
      rmaxs = -9999.0
      nsea = 0
      do j = 1,jm
       do i = 1,im
        if (chlsat(i,j,id) .gt. flagm .and. 
     *      chlsat(i,j,id) .lt. flagl)then
         rmins = min(rmins,chlsat(i,j,id))
         rmaxs = max(rmaxs,chlsat(i,j,id))
         nsea = nsea+1
        endif
       enddo
      enddo
      if (nsea .le. 0)then
       write(6,*)'No Assimilation Data for Day, Year ',iday,iyr
       go to 20
      else
       write(6,*)'Assimilated data available'
       write(6,*)'Min, max, no. Assimilation data = ',rmins,rmaxs,nsea
      endif
c
c  Create im,jm oriented model data
      chlmod = 0.0
      k = 1
      m = indext1
      do n = nnut+1,ntyp-nzoo
       nw = 0
       do j = 1,jm
        do i = 1,im
         if (Depth(i,j) .gt. 200.0)then
          nw = nw+1
          chlmod(i,j) = chlmod(i,j) + P(nw,k,m,n)
         endif
        enddo    !i
       enddo     !j
      enddo       !n

c  Convert to log; only use Assimilation points which are co-located
c  with model points
      chls = 0.0
      do j = 1,jm
       do i = 1,im
        if (chlmod(i,j) .gt. flagm .and. chlmod(i,j) .lt. flagl)then
         chlmod(i,j) = log10(chlmod(i,j)) + 3.0
        else
         chlmod(i,j) = 0.0
        endif
        if (chlsat(i,j,id) .gt. flagm .and. 
     *      chlsat(i,j,id) .lt. flagl .and. 
     *      chlmod(i,j) .lt. flagl)then
         chlsmonwt = (1.0-wtmon)*chlsat(i,j,id) + wtmon*achlmon(i,j)
         chls(i,j) = log10(chlsmonwt) + 3.0
         chls(i,j) = (1.0-wtmod(i,j))*chls(i,j) + wtmod(i,j)*chlmod(i,j)
        endif
       enddo
      enddo
c  Exclude points
      nw = 0
      do j = 1,jm
       do i = 1,im
        if (Depth(i,j) .gt. 200.0)nw = nw+1
c   ice exclusion
        if (picec(nw) .gt. 0.0)chls(i,j) = 0.0
c        chlstmp = 0.0
c        if (chls(i,j) .gt. flagm)chlstmp = 10.0**(chls(i,j)-3.0)
c        chlmodtmp = 0.0
c        if (chlmod(i,j) .gt. flagm)chlmodtmp = 10.0**(chlmod(i,j)-3.0)
c        if (chlstmp .gt. 1.0)then
c         write(6,*)'Assim: id,i,j,nw = ',id,i,j,nw
c         write(6,10)'Assim: chlsmonwt,chlmod = ',chlstmp,chlmodtmp,
c     *nachl(i,j)
c        endif
       enddo
      enddo
c
c  Laplacian of the satellite data
      call laplace(im,jm,chlmod,flap)
c
c  Blended field
      do j = 1,jm
       do i = 1,im
        if (chls(i,j) .gt. flagm)then
         cblend(i,j) = chls(i,j)
        else
         cblend(i,j) = chlmod(i,j)
        endif
       enddo
      enddo
c
c  Perform blending
      call blendsub(im,jm,epsilon,flap,chls,flapb,cblend)
c
c  Remove log
      do j = 1,jm
       do i = 1,im
        cblend(i,j) = 10.0**(cblend(i,j)-3.0)
       enddo
      enddo
c
c  Re-initialize
      m = indext1
      m2 = indext3
      k = 1
      nw = 0
      do j = 1,jm
       do i = 1,im
        if (Depth(i,j) .gt. 200.0)then
         nw = nw+1
         tot = 0.0
         do n = nnut+1,ntyp-nzoo
          tot = tot + P(nw,k,m,n)
         enddo
         tot = max(tot,1.0E-36)
         if (picec(nw) .le. 0.0)then   !do not assim over ice
          do n = nnut+1,ntyp-nzoo
           Pmod(n-nnut) = P(nw,k,m,n)
           Pmod(n-nnut) = max(Pmod(n-nnut),1.0E-5)
           P(nw,k,m,n) = cblend(i,j)/tot*P(nw,k,m,n)
           adj = P(nw,k,m,n)/Pmod(n-nnut)
           P(nw,k,m2,n) = P(nw,k,m2,n)*adj
          enddo
         endif
        endif
       enddo
      enddo
c
      go to 20
19    write(6,*)'No Assimilation Data for Month, Year ',imo,iyr
      close(81)
20    continue
10    format(a27,2f15.5,i5)
      return
      end
c
c **********************************************************************
      subroutine laplace(igrd,jgrd,data,flap)
c
c  Computes Laplacian of input data file.  Outputs Laplacian data
c  file for use in CRAM (blending).  Does not include delta x and 
c  delta y terms because they cancel out in CRAM method.
c
      real data(igrd,jgrd), flap(igrd,jgrd)
c
c  Compute Laplacian (ignore deltax,y terms in denominator -- will
c  cancel out in the matrix solution)
c   Laplacian
      do j = 2,jgrd-1
       do i = 2,igrd-1
        term1 = data(i,j-1)
        term2 = data(i-1,j)
        term3 = -4.0*data(i,j)
        term4 = data(i+1,j)
        term5 = data(i,j+1)
        flap(i,j) = term1 + term2 + term3 + term4 + term5
       enddo
      enddo
      j = 1
      do i = 2,igrd-1
       term1 = 0.0
       term2 = data(i-1,j)
       term3 = -4.0*data(i,j)
       term4 = data(i+1,j)
       term5 = data(i,j+1)
       flap(i,j) = term1 + term2 + term3 + term4 + term5
      enddo
      j = jgrd
      do i = 2,igrd-1
       term1 = data(i,j-1)
       term2 = data(i-1,j)
       term3 = -4.0*data(i,j)
       term4 = data(i+1,j)
       term5 = 0.0
       flap(i,j) = term1 + term2 + term3 + term4 + term5
      enddo
      i = 1
      do j = 2,jgrd-1
       term1 = data(i,j-1)
       term2 = data(igrd,j)
       term3 = -4.0*data(i,j)
       term4 = data(i+1,j)
       term5 = data(i,j+1)
       flap(i,j) = term1 + term2 + term3 + term4 + term5
      enddo
      i = igrd
      do j = 2,jgrd-1
       term1 = data(i,j-1)
       term2 = data(i-1,j)
       term3 = -4.0*data(i,j)
       term4 = data(1,j)
       term5 = data(i,j+1)
       flap(i,j) = term1 + term2 + term3 + term4 + term5
      enddo
      j = 1
      i = 1
      term1 = 0.0
      term2 = data(igrd,j)
      term3 = -4.0*data(i,j)
      term4 = data(i+1,j)
      term5 = data(i,j+1)
      flap(i,j) = term1 + term2 + term3 + term4 + term5
      j = 1
      i = igrd
      term1 = 0.0
      term2 = data(i-1,j)
      term3 = -4.0*data(i,j)
      term4 = data(igrd,j)
      term5 = data(i,j+1)
      flap(i,j) = term1 + term2 + term3 + term4 + term5
      j = jgrd
      i = 1
      term1 = data(i,j-1)
      term2 = data(igrd,j)
      term3 = -4.0*data(i,j)
      term4 = data(i+1,j)
      term5 = 0.0
      flap(i,j) = term1 + term2 + term3 + term4 + term5
      j = jgrd
      i = igrd
      term1 = data(i,j-1)
      term2 = data(i-1,j)
      term3 = -4.0*data(i,j)
      term4 = data(igrd,j)
      term5 = 0.0
      flap(i,j) = term1 + term2 + term3 + term4 + term5
c
      return
      end
c
c *********************************************************************
      subroutine blendsub(igrd,jgrd,epsilon,flap,chls,flapb,cblend)
c
c  Subroutine that actually performs the blended analysis.
c
      real chls(igrd,jgrd)
      real flap(igrd,jgrd)
      real flapb(igrd,jgrd)
      real cblend(igrd,jgrd)
c
c
c  Use iteration to converge blended product to pre-specified limit
      alpha = 0.1
      rmax1 = 9999.0
      m = 0
      do while (rmax1 .gt. epsilon)
       m = m+1
       call laplace(igrd,jgrd,cblend,flapb)
       rmax = -9999.0
       rmin = 9999.0
       do j = 1,jgrd
        do i = 1,igrd
         if (chls(i,j) .gt. 0.0)then
          cblend(i,j) = chls(i,j)
         else
          R = flapb(i,j) - flap(i,j)
          rmax = max(rmax,R)
          rmin = min(rmin,cblend(i,j))
          alphaR = alpha*R
          if (R .lt. 0.0 .and. abs(alphaR) .ge. cblend(i,j))then
           cblend(i,j) = cblend(i,j)
          else
           cblend(i,j) = cblend(i,j) + alpha*R
          endif
         endif
        enddo
       enddo
       rmax1 = rmax
c       if (mod(m,200) .eq. 0)write(6,*)'m, Max/Min R = ',m, rmax,rmin
      enddo
      write(6,*)'Final no iter., Max/Min R = ',m, rmax,rmin
c
      return
      end
c
c ***********************************************************************
      subroutine assimwts(imo,chlmod)
c
c  Computes weighting functions for model vs. data for assimilation.
c  Performed regionally and monthly.
c
#include "gloparam.F.h"
      integer icst(im,jm)
      real chlmod(im,jm)
      common /bwt/ wtmod(im,jm)
      common /blonlat/ xlonb(im),ylatb(jm)
      common /bsatmon/ achlmon(im,jm),sumonchl(im,jm)
      data flag /999.0E9/
c
      wtmod = 0.10
c
c   Coastal border points
      icst = 0
      do j = 1,jm
       jst = max(j-1,1)
       jnd = min(j+1,jm)
       do i = 1,im
        if (chlmod(i,j) .gt. 0.0)then
         ist = max(i-1,1)
         ind = min(i+1,im)
         do jj = jst,jnd
          do ii = ist,ind
           if (chlmod(ii,jj) .le. 0.0)then
            icst(i,j) = 1
           endif
          enddo
         enddo
        endif
       enddo
      enddo
c
      do j = 1,jm
       do i = 1,im
c   Antarctic uses global mean except in Ross Sea and Mackenzie Bay
        if (ylatb(j) .lt. -55.0)then
         wtmod(i,j) = 0.25
        endif
c  Argentina coast
        if (ylatb(j) .lt. -40.0 .and. ylatb(j) .gt. -50.0)then
         if (xlonb(i) .gt. -60.0 .and. xlonb(i) .le. -50.0)then
          if (icst(i,j) .eq. 1)wtmod(i,j) = 0.75
         endif
        endif
c   North and Equatorial Indian are contminated by dust
        if (ylatb(j) .gt. -11.0)then
         if (xlonb(i) .gt. 30.0 .and. xlonb(i) .le. 115.0)then
          if (achlmon(i,j) .gt. 1.0)then
           wtmod(i,j) = 0.85
          else if (achlmon(i,j) .gt. 0.5)then
           wtmod(i,j) = 0.75
          else
           wtmod(i,j) = 0.5
           if (icst(i,j) .eq. 1)wtmod(i,j) = 0.75
          endif
         endif
        endif
c   Congo River high chl is mostly CDOM
        if (ylatb(j) .gt. -20.0 .and. ylatb(j) .lt. 10.0)then
         if (xlonb(i) .ge. 0.0 .and. xlonb(i) .lt. 30.0)then
          wtmod(i,j) = 0.5
          if (achlmon(i,j) .gt. 2.0)wtmod(i,j) = 0.75
          if (icst(i,j) .eq. 1)wtmod(i,j) = 0.75
         endif
        endif
c   Offshore Mauritania is contaminated by dust
        if (ylatb(j) .gt. 0.0 .and. ylatb(j) .lt. 30.0 .and.
     *      xlonb(i) .gt. -30.0 .and. xlonb(i) .lt. 10.0)then
         wtmod(i,j) = 0.5
         if (achlmon(i,j) .gt. 1.0)wtmod(i,j) = 0.85
         if (icst(i,j) .eq. 1)wtmod(i,j) = 0.90
        endif
c   Also Amazon plume is mostly CDOM
        if (ylatb(j) .gt. -5.0 .and. ylatb(j) .lt. 15.0 .and.
     *      xlonb(i) .gt. -60.0 .and. xlonb(i) .lt. -40.0)then
         wtmod(i,j) = 0.5
         if (achlmon(i,j) .gt. 2.0)wtmod(i,j) = 0.75
        endif
c   North Pacific
        if (imo .eq. 5 .or. imo .eq. 6)then
         if (ylatb(j) .gt. 50.0)then
          if (xlonb(i) .ge. 160.0)then
           wtmod(i,j) = 0.25
           if (icst(i,j) .eq. 1)wtmod(i,j) = 0.75
          endif
          if (xlonb(i) .ge. -180.0 .and. xlonb(i) .lt. -120.0)then
           wtmod(i,j) = 0.25
          endif
         endif
        endif
c   Sea of Okhotsk
        if (imo.eq.4 .or. imo.eq.5 .or. imo.eq.6)then
         if (ylatb(j) .gt. 40.0)then
          if (xlonb(i) .ge. 140.0 .and. xlonb(i) .lt. 155.0)then
           wtmod(i,j) = 0.75
           if (icst(i,j) .eq. 1)wtmod(i,j) = 0.75
          endif
         endif
         if (ylatb(j) .gt. 50.0)then
          if (xlonb(i) .ge. 150.0 .and. xlonb(i) .lt. 160.0)then
           wtmod(i,j) = 0.75
          endif
         endif
        endif
cc  North Central Pacific
c        if (imo .eq. 4 .or. imo .eq. 5)then
c         if (ylatb(j) .gt. 30.0 .and. ylatb(j) .lt. 50.0)then
c          if (xlonb(i) .ge. 140.0)then
c           wtmod(i,j) = 0.25
c          endif
c          if (xlonb(i) .ge. -180.0)then
c           wtmod(i,j) = 0.90
c           if (icst(i,j) .eq. 1)wtmod(i,j) = 0.75
c          endif
c         endif
c        endif
c  Make sure coastal regions are at least the minimum
        if (icst(i,j) .eq. 1)then
         if (wtmod(i,j) .lt. 0.5)wtmod(i,j) = 0.5
        endif
 
       enddo
      enddo
c
      return
      end