subroutine AFFORESTATION (LAND_PTS, LAND_INDEX, DZ_SOIL, M
     &,                         FRAC)

#if defined O_mtlm

! This is a climate engineering subroutine for irrigating N. Africa
! and Australia.

      implicit none

      include "size.h"
      include "csbc.h"

! LAND_PTS   = IN Number of land points.
! LAND_INDEX = IN Indices of landpoints.
! GREEN_THETA_SAT = Saturated volumetric soil moisture
!                   concentration resulting from irrigation
! DZ_SOIL = IN Soil layer thickness (m).
! FRAC_VS = fraction of grid cell containing veg


      integer LAND_PTS, LAND_INDEX(POINTS)
      integer I, J, L, N
      real GREEN_THETA_SAT, M(POINTS) , DZ_SOIL
      real FRAC(POINTS,NTYPE)

      GREEN_THETA_SAT = 0.40

# if defined O_set_afforest_veg
!	To force a land (veg) type use the following example
!   as a template. Where:
! NTYPE = Number of surface types.
! Land surface types :
!     1 - Broadleaf Tree
!     2 - Needleleaf Tree
!     3 - C3 Grass
!     4 - C4 Grass
!     5 - Shrub
!     6 - Soil
!
! Example for forcing Broadleaf trees (N=1):
!        FRAC(L,1) = 1.0
!        do N=2,NTYPE
!           FRAC(L,N)= 0.0
!        enddo      
# endif

# if defined O_green_australia
! Irrigate Australia and watch the plants grow 
      do J=29,44
       do I=33,43
        L = land_map(I,J)
        M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
!        print*,'M(L)=',M(L)
       enddo
      enddo
# endif

# if defined O_green_africa
! Irrigate N. Africa and watch the plants grow
! VSAT = 0.458m3 H2O/m3 soil = theta_sat
! RHO_WATER = 1000kg/m3
! ROOTDEP = dz_soil
! First do the Western part
      do J=61,72
       do I=96,101
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
      enddo
! Then most of the Eastern part
      do J=61,72
       do I=1,9
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
      enddo
! Now fill in an area to the south
      J=60
       do I=4,11
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
! Now fill in more Eastern area (Eritrea and Sudan)
      do J=61,69
       do I=10,11
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
      enddo
! Now fill in more Eastern area (Egypt and Sudan)
      do J=60,64
         I=12
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
# if defined O_green_Saudi_Arabia
!
      do J=60,68
       do I=12,14
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
      enddo
!
      do J=60,68
         I=15
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
      do J=61,66
       do I=16,17
          L = land_map(I,J)
          M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
# if defined O_set_afforest_veg
        FRAC(L,6) = 1.0
        do N=1,5
           FRAC(L,N)= 0.0
        enddo
# endif
       enddo
      enddo
# endif
# endif

      return
      end