createtrees Subroutine

    subroutine createtrees
    use modglobal,  only : ltrees,ntrees,tree,cexpnr,ifinput,zh,zf,dzh,dzfi,dzhi,dzf,Qstar,&
                           dec,lad,kb,ke,cp,rhoa,ntree_max,dQdt,tr_A,dy,xh
    use modfields,  only : um,vm,wm,thlm,qt0,svp,up,vp,wp,thlp,qtp,Rn,clai,qc,qa,ladzh,ladzf
    use modmpi,     only : myid,comm3d,mpierr,MY_REAL
    use modsurfdata,only : wtsurf
    use modibmdata ,only : bctfz
    implicit none
    integer :: n,k
    real :: Rq
    character(80) chmess

    if ((ltrees .eqv. .false.) .or. (ntrees==0)) return

      allocate(tree(ntrees,6))

      ! read global trees
      if(myid==0) then
        ! write(*,*) '1, myid, ntrees, ltrees, cexpnr', myid, ntrees, ltrees, cexpnr
        if (ntrees>0) then
          open (ifinput,file='trees.inp.'//cexpnr)
          read (ifinput,'(a80)') chmess
          read (ifinput,'(a80)') chmess
          do n=1,ntrees
            read (ifinput,*) &
                  tree(n,1), &
                  tree(n,2), &
                  tree(n,3), &
                  tree(n,4), &
                  tree(n,5), &
                  tree(n,6) 
          end do

          ! write (6,*) 'Tree number,   il, iu, jl, ju, kl, ku '
          ! do n=1,ntrees
          !   write (1,*) &
          !         n , &
          !         tree(n,1), &
          !         tree(n,2), &
          !         tree(n,3), &
          !         tree(n,4), &
          !         tree(n,5), &
          !         tree(n,6)                 
          ! end do
          close (ifinput)
        end if

      end if ! end if myid==0

      ! call MPI_BCAST(ntrees ,1,MPI_INTEGER ,0,comm3d,mpierr)
      call MPI_BCAST(tree ,6*ntrees,MPI_INTEGER ,0,comm3d,mpierr)

      !! capability to read lad of trees from a lad.inp.xxx file
      !! tg3315 commented - for now assume constant lad or hardcode as seen below
      !if (myid==0) then

      !  allocate(ladt(nlad)) 

      !  if (nlad>0) then
      !    open (ifinput,file='lad.inp.'//cexpnr)
      !    read (ifinput,'(a80)') chmess
      !    do n=1,nlad
      !      read (ifinput,*) &
      !            ladt(n)
      !      end do

      !    write (*,*) 'lad'
      !    do n=1,nlad
      !      write (1,*) &
      !            n , &
      !            ladt(n)
      !    end do
      !  end if

      ! input into leaf area density array (cell faces)
      !lad(tree(1,5):tree(1,6)+1) = ladt !defined at cell faces

      ! initialise tree variables that are constant in time
      ! calculate clai for tallest possible tree
      ! assume no grid stretching in lowest layer with trees
      if (myid==0) then

        ntree_max = maxval(tree(:,6))-minval(tree(:,5))+1

        ! hard code non-uniform lad profiles - temporary
        if (.false.) then
          ladzh(1:ntree_max) = (/  0.1, 0.11, 0.12, 0.14, 0.16, 0.19, 0.22, 0.265, 0.31, 0.335, 0.36, 0.37, 0.38, 0.3725, 0.365, 0.3375, 0.31, 0.2375, 0.165, 0.0875 /) !LAI =2 from Shaw, 1992
          ! assumes lad at ntree_max+1 is 0.!
        else
          ladzh(1:ntree_max+1) = lad !up to ntree_max+1 cos ladzh is at cell faces
        end if

        ! interpolate to find lad at cell centres
        do k=1,ntree_max
          ladzf(k) = 0.5*(ladzh(k)+ladzh(k+1))
        end do

        ! clai at cell faces using lad at cell centre
        do k = 1,ntree_max
          clai(k) = sum(0.5*(ladzh(k:ntree_max)+ladzh(k+1:ntree_max+1))*(dzf(maxval(tree(:,6))-ntree_max+k:maxval(tree(:,6)))))
        end do

        ! Net radiation at cell faces! !W/m^2
        do k = 1,ntree_max+1
          Rn(k) = Qstar*exp(-dec*clai(k))
        end do

        ! Change in radiation over each layer (at cell centres) W/m^2
        do k = 1,ntree_max
          qc(k) = Rn(k+1) - Rn(k)
        end do

        ! incorporate the distributed storage term - updated for zero heat storage in tree
        do k = 1,ntree_max
          qa(k) = qc(k) !(1 - Rq*0.11) * qc(k) - 0.11*Rq*dQdt + Rq*12.3
        end do

        ! write relevant fields
        ! write(*,*) 'ntree_max', ntree_max
        ! write(*,*) 'qa', qa(1:ntree_max)
        ! write(*,*) 'qc', qc(1:ntree_max)
        ! write(*,*) 'clai', clai(1:ntree_max+1)
        ! write(*,*) 'Rn', Rn(1:ntree_max+1)
        ! write(*,*) 'ladzf', ladzf
        ! write(*,*) 'ladzh', ladzh

      end if

      !apply storage to set wtsurf and bctfz as a function of Qstar
      wtsurf = -((1-0.7)*Qstar-0.33*dQdt+38)/(rhoa*cp)
      bctfz  = -((1-0.7)*Qstar-0.33*dQdt+38)/(rhoa*cp)

      ! calc tree cover area - used to apply steady-state boundary condition
      if (myid==0) then
        do n = 1,ntrees
          tr_A = tr_A + (xh(tree(n,2)+1) - xh(tree(n,1)))*dy*(tree(n,4) - tree(n,3) + 1)
        end do
      end if

      ! write updated ground and roof heat fluxes
      write(*,*) 'wtsurf', wtsurf
      write(*,*) 'bctfz', bctfz

      ! broadcast variables
      call MPI_BCAST(ntree_max , 1,MY_REAL ,0,comm3d,mpierr)
      call MPI_BCAST(tr_A , 1,MY_REAL ,0,comm3d,mpierr)
      call MPI_BCAST(Rn , ke+1,MY_REAL ,0,comm3d,mpierr)
      call MPI_BCAST(qc , ke+1,MY_REAL ,0,comm3d,mpierr)
      call MPI_BCAST(clai  , ke+1,MY_REAL ,0,comm3d,mpierr)  
      call MPI_BCAST(qa  , ke+1,MY_REAL ,0,comm3d,mpierr)
      call MPI_BCAST(ladzf  , ke+1,MY_REAL ,0,comm3d,mpierr)
      call MPI_BCAST(ladzh , ke+1,MY_REAL ,0,comm3d,mpierr)

    end subroutine createtrees