program source
c
c Program which defines the dust sources as a function 
c of orography and vegetation
c
      parameter (nx=360,ny=180,nveg=34)
c
      character *1   FLG(nx,ny)
      character *5  veget(nveg)
      character *5  data(nx)
c
      dimension nel1(nx,ny),nel2(nx,ny)
      dimension vdlg(nx,ny),dxy(nx,ny),rlate(ny)
      dimension ylat(ny)
      dimension ifrac_vu(34)
c
c    Value  Land Cover Class
c    -----  -------------------------------------------
c    0      water
c    1      broadleaf evergreen forest 
c    2      coniferous evergreen forest and woodland
c    3      high latitude deciduous forest and woodland
c    4      tundra
c    5      mixed coniferous forest and woodland
c    6      wooded grassland
c    7      grassland
c    8      bare ground
c    9      shrubs and bare ground
c    10     cultivated crops
c    11     broadleaf deciduous forest and woodland
c    12     data unavailable
C
C Define latitiudes of data grid pints
C
      do j=1,180
        ylat(j)=-89.5+(j-1)
      enddo
C
C Read DeFries vegetation classes on a 1x1 grid (GISS)
C http://www.geog.umd.edu/landcover/global-cover.html
C
      open(1,file='landcover.txt',status='old')
      do j=1,ny
        do i=1,nx
          read(1,*) nel2(i,ny-j+1)
          if (nel2(i,ny-j+1).eq.8) then
            nel2(i,ny-j+1)=10
          else
            if (nel2(i,ny-j+1).eq.9) then
              nel2(i,ny-j+1)=2
            else
              nel2(i,ny-j+1)=0
            endif
          endif
        enddo
      enddo
      close(1)
C
C Read orography distribution on a 1x1 grid (GISS)
C http://www.giss.nasa.gov/data/
C NB: There is also the NAVY dataset on 10'x10'
      OPEN(1,FILE='orography',status='old')
      READ(1,1010) ((NEL1(I,J),FLG(I,J),I=1,360),J=1,180)
      CLOSE(1)
      DO j=1,ny
        DO i=1,nx/2
          vdlg(i,j)=AMAX1(0.,FLOAT(NEL1(i+nx/2,j)))
          vdlg(i+nx/2,j)=AMAX1(0.,FLOAT(NEL1(i,j)))
        enddo
      enddo
c
c Surface
      pi=3.141592
      re=6371e3
      do j=1,ny
        rlate(j)=(-90.+(j-1))*pi/180.
      enddo
      do i=1,nx
        do j=2,ny
          dxy(i,j)=2.*pi*re**2*(sin(rlate(j))-sin(rlate(j-1)))/nx
        enddo
        dxy(i,1)=dxy(i,ny)
      enddo
c
      call filter(ylat,nel2,vdlg)
c
      totsur=0.
      do j=1,ny
        do i=1,nx
          totsur=totsur+vdlg(i,j)*dxy(i,j)
        enddo
      enddo
      write(*,*) 'Tot surface=',totsur
      open(1,file='source',form='unformatted')
      write(1) ((vdlg(i,j),i=1,nx),j=1,ny)
      close(1)
c
 1010 FORMAT(12(I5,A1))
  101 format(10(3x,a5))
c
      stop
      end
c
      subroutine filter(ylat,nel2,vdlg)
c
c Weight each grid point as a function of its relative
c altitude and other parameters on a surrounding domain 
c of 2*ip (x-direction) and 2*jp (y-direction)
c
      parameter (nx=360,ny=180)
      dimension vdlg0(nx,ny)
      dimension vdlg(nx,ny),ylat(ny)
      dimension nel2(nx,ny)
c
c Initialize
c
      do j=1,ny
        do i=1,nx
          vdlg0(i,j)=vdlg(i,j)
          vdlg(i,j)=0.
        enddo
      enddo
      niter=1
c
c Perform niter iteration to enlarge the surrounding domain
c which will reinforce the weigth of large bassins compare
c to small depressions
c
      do it=1,niter
c
c Fix the number of grid points defining the size of the
c surrounding domain
c
        ip=6+it
        jp=6+it
c
c Loop over the grid points of the global domain
c Warning: Siberia is cut off (to avoid extend the domain
c on both sides of the dateline)
c
        do 20 j1=jp,ny-jp+1
        do 20 i1=ip,nx-ip+1
          zmin=1.e30
          zmax=-1.e30
          av=0.
          iav=0
c
c Loop over the grid points of the surrounding domain
c
          do i2=i1-ip,i1+ip
           do j2=j1-jp,j1+jp
             iav=iav+1
             av=av+vdlg0(i2,j2)
c
c Local altitude extrema
c
             zmax=amax1(zmax,vdlg0(i2,j2))
             zmin=amin1(zmin,vdlg0(i2,j2))
           enddo
          enddo 
          av=av/iav
          delta=amin1(8000.,zmax-zmin)
          deltak=delta/1000.
          if (delta.gt.0.and.vdlg0(i1,j1).gt.0.) then
c
c Weight depending on the gradient
c
            alpha=(8.-0.25*deltak)/8.
            alpha=1
c
c Weight depending on the relative altitude
c
            delta=(zmax-vdlg0(i1,j1))/(zmax-zmin)
            delta=delta**3
c
c Weigth depending on the vegetation
c
            epsilon=nel2(i1,j1)/10.
c
c Weigth depending on the latitude
c
            if (abs(ylat(j1)).gt.60.) then
              phi=0.
            else
              phi=1.
            endif
c
c Potential surface available for dust uplifting (source)
c
            vdlg(i1,j1)=vdlg(i1,j1)+
     >       alpha*delta*epsilon*phi/niter
c
          endif 
   20   continue
      enddo
      return
      end