module veg_drv 6,1
!@sum veg_drv contains variables and routines for vegetation driver
!@auth I. Alienov, N. Kiang

      use model_com, only : im,jm
      implicit none
      private
      save

      public init_vegetation,reset_veg_to_defaults
     &     ,veg_set_cell, veg_save_cell,upd_gh

      real*8,public :: cosday,sinday

      contains



      subroutine init_vegetation(redogh,istart) 1,17
!@sum reads vegetation arrays and initializes the vegetation
      use filemanager
      use param
      use constant, only : twopi,one
      use DOMAIN_DECOMP, only : GRID, GET
      use model_com, only : fearth,jeq,jyear
      use veg_com !, only : vdata,Cint,Qfol
      use ghycom, only : ngm
      use vegetation, only : cond_scheme,crops_yr
      use ghycom, only : dz_ij
      use surf_albedo, only: albvnh  !nyk !!! not used? i.a.

      implicit none

      integer, intent(in) :: istart
      logical, intent(in) :: redogh
!---  local vars
      integer iu_veg
      integer i, j, k
      real*8 dif,frdn,frup,pearth,phase,scs0,scsim,scsre,sfv,sla0
      real*8 almass0, almassre, almassim  !nyk
!      real*8 sla0f, slimf, slref !nyk
      real*8 slim,slre,svh,z
      real*8 snm,snf  ! temporary sums (adf)
      real*8 cwc_sum
      integer iv, l
      logical veg_data_missing
      real*8 fv,dz(ngm)
      integer n
!---  parameters for different types of vegetation
!---          tundr  grass  shrub  trees  decid evrgr  rainf crops
      real*8, parameter :: alamax(8) =
     $     (/ 1.5d0, 2.0d0, 2.5d0, 4.0d0, 6.0d0,10.0d0,8.0d0,4.5d0/)
      real*8, parameter :: alamin(8) =
     $     (/ 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 8.0d0,6.0d0,1.0d0/)

      real*8, parameter :: aroot(8) =
     $     (/ 12.5d0, 0.9d0, 0.8d0,0.25d0,0.25d0,0.25d0,1.1d0,0.9d0/)
      real*8, parameter :: broot(8) =
     $     (/  1.0d0, 0.9d0, 0.4d0,2.00d0,2.00d0,2.00d0,0.4d0,0.9d0/)
      real*8, parameter :: rsar(8) =
     $     (/100d0, 100d0, 200d0, 200d0, 200d0, 300d0,250d0, 125d0/)
      real*8, parameter :: vhght(8) =
     $     (/0.1d0, 1.5d0,   5d0,  15d0,  20d0,  30d0, 25d0,1.75d0/)
! Mean canopy nitrogen (nmv; g/m2) and Rubisco factors (nfv) for each
! vegetation type (adf)
      real*8, parameter :: nmv(8) =
     $     (/1.6d0,0.82d0,2.38d0,1.03d0,1.25d0,2.9d0,2.7d0,0.82d0/)
      real*8, parameter :: nfv(8) =
     $     (/1.4d0,1.5d0 ,1.3d0 ,1.3d0 ,1.5d0 ,0.9d0,1.1d0,1.5d0 /)
      integer, parameter :: laday(8) =
     $     (/ 196,  196,  196,  196,  196,  196,  196,  196/)
      real*8, parameter :: can_w_coef(8) =
     &     (/ 1.d-4, 1.d-4, 1.d-4, 1.d-4, 1.d-4, 1.d-4, 1.d-4, 1.d-4 /)

! Specific leaf areas (sleafa, kg[C]/m2) (adf, nyk)
! Values below 1/(m2/kg) to get kg/m2 for multiplying.
! Sources: White, M.A., et.al. (2000), Earth Interactions, 4:1-85.
!          Leonardos,E.D.,et.al.(2003), Physiologia Plantarum, 117:521+.
!               From winter wheat grown at 20 C (20 m2/kg[dry mass])
!                                      and  5 C (13 m2/kg[dry mass])
!          Francesco Tubiello, personal communication, crop 18-20 m2/kg.
      real*8, parameter :: sleafa(8) =
     $     1./(/30.5d0,49.0d0,30.5d0,40.5d0,32.0d0,8.2d0,32.0d0,18.0d0/)

c****             tundr grass shrub trees decid evrgr rainf crops
c****
c**** laday(veg type, lat belt) = day of peak lai
c old peak lai:  2nd line is for latitudes < 23.5 deg
c****    1  temperate latitudes
c****    2  non-temperate latitudes
c     data  laday/ 196,  196,  196,  196,  196,  196,  105,  196/
c     data  laday/ 196,  288,  288,  288,  288,  196,  105,  288/
c****
c**** contents of ala(k,i,j),  lai coefficients
c****   1  average leaf area index
c****   2  real amplitude of leaf area index
c****   3  imaginary amplitude of leaf area index
c****
c**** contents of almass(i,j),  leaf mass
c****      leaf mass = ala*sleafa = kg[C]/ground area
c****
c**** contents of acs(k,i,j),  cs coefficients
c****   1  average stomatal conductance
c****   2  real amplitude of stomatal conductance
c****   3  imaginary amplitude of stomatal conductance
c****
!@dbparam ghy_default_data if == 1 reset all GHY data to defaults
!@+ (do not read it from files)
      integer :: ghy_default_data = 0
      integer :: northsouth  !1=south, 2=north hemisphere

      INTEGER :: I_0, I_1, J_1, J_0
      INTEGER :: J_0S, J_1S, J_0STG, J_1STG
      LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S,
     &               J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG,
     &               HAVE_SOUTH_POLE = HAVE_SOUTH_POLE,
     &               HAVE_NORTH_POLE = HAVE_NORTH_POLE)

c**** read rundeck parameters
      call sync_param( "cond_scheme", cond_scheme)  !nyk 5/1/03
      call sync_param( "crops_yr", crops_yr)

c**** read land surface parameters or use defaults
      call openunit("VEG",iu_VEG,.true.,.true.)
      do k=1,10                 !  11 ????
        call readt (iu_VEG,0,vdata(1,1,K),im*jm,vdata(1,1,k),1)
      end do
c**** zero-out vdata(11) until it is properly read in
      vdata(:,:,11) = 0.
      call closeunit(iu_VEG)
C**** Update vegetation file if necessary (i.e. crops_yr =0 or >0)
      if(crops_yr.eq.0) call updveg(jyear,.false.)
      if(crops_yr.gt.0) call updveg(crops_yr,.false.)

      if (istart.le.2 .or. redogh) then ! initial. foliage arrays (adf)
        Cint(:,:)=0.0127D0      ! internal CO2
        Qfol(:,:)=3.D-6         ! surface mixing ratio
        cnc_ij(:,:) = 0.d0
      end if
      if (istart.le.0) return   ! avoid reading unneeded files

     !!! spgsn=.1d0

c**** check whether ground hydrology data exist at this point.
      veg_data_missing = .false.
      do j=J_0,J_1
        do i=1,im
          if (fearth(i,j).gt.0) then
            if ( sum(vdata(i,j,1:10)).eq.0 ) then
              print *,"No vegetation data: i,j=",i,j,vdata(i,j,1:10)
              veg_data_missing = .true.
            end if
          end if
        enddo
      enddo
      if ( veg_data_missing ) then
        write(6,*) 'Vegetation data is missing at some pts'
        write(6,*) 'If you have a non-standard land mask, please'
        write(6,*) 'consider using extended GH data and rfs file.'
        call stop_model(
     &       'Vegetation data is missing at some cells',255)
      endif


      entry upd_gh  1
c****
c**** set the global arrays  ala, acs, afb, afr, anm, anf
c****
      ala(:,:,:)=0.
      alaf(:,:,:,:)=0.          !nyk lai by veg functional type
      alaif(:,:,:)=0.
      acs(:,:,:)=0.
      avh(:,:)=0.
      afb(:,:)=0.
      afr(:,:,:)=0.
      acs(1,:,:)=.01d0
      anm(:,:)=0.0D0            ! Global mean canopy N array (adf)
      anf(:,:)=0.0D0            ! Global Rubisco factors (adf)
      almass(:,:,:)=0.0D0       ! Global leaf mass at a time, nyk
!rar  aalbveg(:,:)=0.08D0 ! no need, it is set in daily_earth
      can_w_capacity(:,:) = 0.d0

      do j=J_0,J_1
        if(j.le.jm/2) then
          northsouth=1.d0  !southern hemisphere
        else
          northsouth=2.d0  !northern hemisphere
        end if
        do i=1,im
          pearth=fearth(i,j)
          afb(i,j)=vdata(i,j,1)+vdata(i,j,10)
          if(afb(i,j).gt..999) afb(i,j)=1.
          if(pearth.le.0..or.afb(i,j).ge.1.) cycle
c**** calculate lai, cs coefficicents
          sfv=0.d0
          sla0=0.     !For calculating sla
          slre=0.     !For calculating sla
          slim=0.     !For calculating sla
          !sla0f=0.     !For calculating alaf
          !slref=0.     !For calculating alaf
          !slimf=0.     !For calculating alaf
          almass0=0.  !nyk For specific leaf area
          almassre=0. !nyk For specific leaf area
          almassim=0. !nyk For specific leaf area
          scs0=0.
          scsre=0.
          scsim=0.
          svh=0.
          snm=0. ! adf
          snf=0. ! adf
          cwc_sum = 0.d0
          do iv=1,8
            phase=twopi*laday(iv)/365.
            if(j.lt.jeq) phase=phase+twopi/2.
            fv=vdata(i,j,iv+1)
            sfv=sfv+fv
            svh=svh+fv*vhght(iv)
            snm=snm+fv*nmv(iv) ! adf
            snf=snf+fv*nfv(iv) ! adf
            dif=(alamax(iv) - alamin(iv))
            sla0=sla0+fv*(alamax(iv) + alamin(iv))
            slre=slre+fv*dif*cos(phase)
            slim=slim+fv*dif*sin(phase)
            !nyk-------------
            !alaf and alaif
            alaf(1,iv,i,j) = 0.5*(alamax(iv) + alamin(iv))
            alaf(2,iv,i,j) = 0.5*dif*cos(phase)
            alaf(3,iv,i,j) = 0.5*dif*sin(phase)
            alaif(iv,i,j)= alaf(1,iv,i,j)+
     $           cosday*alaf(2,iv,i,j)+sinday*alaf(3,iv,i,j) !save ij
           !nyk-------------
            !almaxmin = almaxmin + sleafa(iv)*fv*(alamax(iv)-alamin(iv))
            almass0 = almass0 + sleafa(iv)*fv*(alamax(iv) - alamin(iv))
            !almass0 = almass0 + sleafa(iv)*fv*(alamax(iv) + alamin(iv))
            !almassre = almassre + sleafa(iv)*fv*dif*cos(phase)
            !almassim = almassim + sleafa(iv)*fv*dif*sin(phase)
            !----------------
            scs0=scs0+fv*(alamax(iv) + alamin(iv))/rsar(iv)
            scsre=scsre+fv*dif*cos(phase)/rsar(iv)
            scsim=scsim+fv*dif*sin(phase)/rsar(iv)
            cwc_sum = cwc_sum + fv*can_w_coef(iv)
          end do
          ala(1,i,j)=.5/sfv*sla0
          ala(2,i,j)=.5/sfv*slre
          ala(3,i,j)=.5/sfv*slim

          acs(1,i,j)=.5/sfv*scs0
          acs(2,i,j)=.5/sfv*scsre
          acs(3,i,j)=.5/sfv*scsim
          avh(i,j)=svh/sfv
          can_w_capacity(i,j) = cwc_sum/sfv
          anm(i,j)=snm/sfv ! adf
          anf(i,j)=snf/sfv ! adf
          almass(1,i,j) = almass0   !This just computes total growth for
          almass(2,i,j) = 0.   ! year via difference between max and min
          almass(3,i,j) = 0.
          !almass(1,i,j)= 0.5/sfv*almass0 !nyk
          !almass(2,i,j)= 0.5/sfv*almassre !nyk
          !almass(3,i,j)= 0.5/sfv*almassim !nyk

c**** calculate root fraction afr averaged over vegetation types
          do n=1,ngm
            dz(n)=dz_ij(i,j,n)
            if(dz(n).le.0.) go to 320
          end do
 320      n=n-1
          do iv=1,8
            fv=vdata(i,j,iv+1)
            z=0.
            frup=0.
            do l=1,n
              z=z+dz(l)
              frdn=aroot(iv)*z**broot(iv)
              frdn=min(frdn,one)
              if(l.eq.n)frdn=1.
              afr(l,i,j) = afr(l,i,j) + fv*(frdn-frup)
              frup=frdn
            end do
          end do
          do l=1,n
            afr(l,i,j) = afr(l,i,j)/(1.-afb(i,j))
          end do
        end do
      end do

      ! write(98,*) afr

c**** recompute ground hydrology data if necessary (new soils data)
      if (redogh) then
        Cint(:,:)=0.0127D0      ! Internal foliage CO2(adf)
        Qfol(:,:)=3.D-6         ! Foliage surface mixing ratio (adf)
        cnc_ij(:,:) = 0.d0
        print *, 'reset vegetation prognostic variables to defaults'
      end if

      return
      end subroutine init_vegetation



      subroutine reset_veg_to_defaults( reset_prognostic ) 1,3
      use veg_com, only: vdata
      use DOMAIN_DECOMP, only : GRID, GET
      !use ghycom
      logical, intent(in) :: reset_prognostic
      integer i,j

      INTEGER :: I_0, I_1, J_1, J_0
      INTEGER :: J_0S, J_1S, J_0STG, J_1STG
      LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S,
     &               J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG,
     &               HAVE_SOUTH_POLE = HAVE_SOUTH_POLE,
     &               HAVE_NORTH_POLE = HAVE_NORTH_POLE)

      do j=J_0,J_1
      do i=1,im

      vdata(i,j,1:11)= (/  0.00000000d+00,  0.00000000d+00,
     &     0.00000000d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.62451148d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.37548852d+00,  0.00000000d+00,  0.00000000d+00 /)
      enddo
      enddo

      end subroutine reset_veg_to_defaults



      subroutine veg_set_cell (i0,j0,ghy_data_only) 2,5
!@sum resets the vegetation module to a new cell i0,j0
      use ghycom, only : ngm
      use sle001, only : fr,snowm,ws,shc,shw
      use vegetation, only : alaie,rs,nm,nf,alai,vh
     &     ,alait,vfraction
     &     ,fdir,parinc,vegalbedo,sbeta,Ci,Qf ! added by adf, nyk
     &     ,cond_scheme         !nyk
     &     ,cnc
      use radncb, only : cosz1
     &    ,FSRDIR,SRVISSURF  !adf, nyk
      use veg_com, only :
     &     Cint,Qfol           ! added by adf
     &     ,cnc_ij
     &     ,aalbveg    ! nyk
     &     ,afr,avh,anm,anf,ala,alaif,alaf,vdata,acs


      implicit none

      integer, intent(in) :: i0,j0
!@var ghy_data_only set only ghy data (no call to vegetatin expected)
      logical, optional :: ghy_data_only
      real*8, parameter :: spgsn=.1d0
      integer l
      real*8 aa
!      real*8 aalbveg0, sfv  !nyk
!      integer northsouth,iv  !nyk
!      real*8 alaic,vh,shtpr,alai  ! adf
      real*8 alaic

c**** fr: root fraction in layer l  (1=fr(1)+fr(2)+...+fr(n))
      do l=1,ngm
        fr(l)=afr(l,i0,j0)
      end do
c**** vh: vegetation height
      vh=avh(i0,j0)
      nm=anm(i0,j0) ! mean canopy nitrogen (g/m2) (adf)
      nf=anf(i0,j0) ! canopy nitrogen factor (adf)
      snowm=vh*spgsn

c**** alai: leaf area index
      alai=ala(1,i0,j0)+cosday*ala(2,i0,j0)+sinday*ala(3,i0,j0)
      alai=max(alai,1.d0)
      !lai by functional type, not normalized by cover fraction!
      do L=1,8
        alaif(L,i0,j0)=alaf(1,L,i0,j0)+
     &       cosday*alaf(2,L,i0,j0)+sinday*alaf(3,L,i0,j0) !save ij
        alait(L) = alaif(L,i0,j0)  !for VEGETATION.f
        vfraction(L)=vdata(i0,j0,L+1) !for VEGETATION.f
      end do

      alaic=5.0
      alaie=alaic*(1.-exp(-alai/alaic))
c**** rs: minimum stomatal resistance
      rs=alai/(acs(1,i0,j0)+cosday*acs(2,i0,j0)+sinday*acs(3,i0,j0))
c???  cnc=alai/rs   redefined before being used (qsbal,cond)

      !---------------------------------------------------------
      !nyk vegetation albedo.  Only really updated daily, but have to
      !get it initialized somewhere after ALBVNH is calculated.
      !ALBVNH is unfortunately set *after* ground hydr is initialized.
      !albvnh(9,6,2)=albvnh(1+8veg,6bands,2hemi), band 1 is VIS.
!      aalbveg0=0.d0               !nyk
!      sfv=0.d0
!      if (j0.le.jm/2) then
!
!       northsouth=1.d0         !southern hemisphere
!      else
!        northsouth=2.d0         !northern hemisphere
!      end if
!      do iv=1,8
!        aalbveg0 = aalbveg0 + fv*(ALBVNH(iv+1,1,northsouth))
!        sfv = sfv + fv
!      end do
!      aalbveg(i0,j0) = aalbveg0/sfv
!      write (99,*) 'aalbveg ghinij', aalbveg(i0,j0) !nyk
      !---------------------------------------------------------

      ws(0,2)=.0001d0*alai
!!!      ws(0,2)=can_w_capacity(i0,j0)*alai
c shc(0,2) is the heat capacity of the canopy
      aa=ala(1,i0,j0)
      shc(0,2)=(.010d0+.002d0*aa+.001d0*aa**2)*shw

! This is a hack to prevent the program from using uninitialized 
! radiation arrays at the beginning of the run.
! If program returns at this point, it sets only the data needed for
! ground hydrology but not the vegetation which is ok if veg_conductance
! is not called.
      if ( present(ghy_data_only) ) then
        if ( ghy_data_only ) return
      endif

!----------------------------------------------------------------------!
      if (cond_scheme.eq.2) then  !new conductance scheme
! Sine of solar elevation (rad).
        sbeta=cosz1(i0,j0)
! adf Fraction of solar radiation at ground that is direct beam.
        fdir=FSRDIR(i0,j0)
! nyk Calculate incident PAR, photosynthetically active radiation.
!     *SRVISSURF is from SRDVIS:
!     = incident visible solar radiation (dir+dif) on the surface
!     = estimated as 53% of total solar flux density, so wavelength
!       range is UV through ~760 or 770 nm cutoff (not strict)
!     *SRDVIS is normalized to solar zenith = 0.
!     *From integrating solar flux density (TOA) over solar spectrum,
!       PAR(400-700 nm) is ~ 82% of the flux density of SRDVIS.
        parinc=0.82*SRVISSURF(i0,j0)*sbeta
! nyk Get vegetation grid albedo, temporary until canopy scheme in place
        vegalbedo = aalbveg(i0,j0)
! Internal foliage CO2 concentration (mol/m3).
        Ci=Cint(i0,j0)
! Foliage surface mixing ratio (kg/kg).
        Qf=Qfol(i0,j0)
        CNC = cnc_ij(i0,j0)
      end if
!----------------------------------------------------------------------!
      return
      end subroutine veg_set_cell



      subroutine veg_save_cell(i,j) 1,2
      !Save vegetation arrays for cell.
      use vegetation, only: cond_scheme, Ci, Qf, CNC
      use veg_com, only:  Cint, Qfol, cnc_ij

      integer, intent(in):: i,j

      if (cond_scheme.eq.2) then      !new conductance scheme only
        Cint(i,j)=Ci                  ! Save last value
        Qfol(i,j)=Qf                  ! Save last value
        cnc_ij(i,j)=CNC
      end if
      end subroutine veg_save_cell



      end module veg_drv

!***********************************************************************


      subroutine updveg (year,reset_veg) 3,10
!@sum  reads appropriate crops data and updates the vegetation file
!@auth R. Ruedy
!@ver  1.0
      USE FILEMANAGER
      USE MODEL_COM, only : im,jm
      USE DOMAIN_DECOMP, only : GRID, GET
      use veg_com, only : vdata
      USE GEOM, only : imaxj
      use veg_drv, only : upd_gh
      implicit none
      integer, intent(in) :: year
      logical, intent(in) :: reset_veg

      real*8 wt,crops         ! temporary vars ! ,crop(im,jm)

      integer :: year1,year2,year_old=-1, iu, i,j,k
      real*8 vdata0(im,jm,11) ! to limit i/o
      real*8  crop1(im,jm)    ! to limit i/o
      real*8  crop2(im,jm)    ! to limit i/o
      save   year1,year2,year_old,vdata0,crop1,crop2,iu      ! to limit i/o

      character*80 title
      real*4 crop4(im,GRID%J_STRT_HALO:GRID%J_STOP_HALO)

      INTEGER :: I_0, I_1, J_1, J_0
      INTEGER :: J_0S, J_1S, J_0STG, J_1STG
      LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S,
     &               J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG,
     &               HAVE_SOUTH_POLE = HAVE_SOUTH_POLE,
     &               HAVE_NORTH_POLE = HAVE_NORTH_POLE)

C**** check whether update is needed
      if (year.eq.year_old) return

C**** first iteration actions:
      if (year_old.lt.0) then
C****     check whether a no-crops vege-file was used
        do j=J_0S,J_1S
        do i=1,im
           if(vdata(i,j,9).gt.0.)
     *     call stop_model('updveg: use no_crops_VEG_file',255)
        end do
        end do
C****     open and read input file
        call openunit('CROPS',iu,.true.,.true.)
        read(iu) title,crop4
        read(title,*) year1
        crop1=crop4 ; crop2=crop4 ; year2=year1
        if (year1.ge.year)          year2=year+1
C****     save orig. (no-crop) vdata to preserve restart-independence
        vdata0 = vdata
      end if

      wt=0.
      do while (year2.lt.year)
         year1 = year2 ; crop1 = crop2
         read (iu,end=10) title,crop4
         read(title,*) year2
         crop2 = crop4
      end do
      wt = (year-year1)/(real(year2-year1,kind=8))
   10 continue
      write(6,*) 'Using crops data from year',year1+wt*(year2-year1)

C**** Modify the vegetation fractions
      do j=J_0,J_1
      do i=1,imaxj(j)
         if (crop1(i,j).ge.0.) then
            crops = crop1(i,j) + wt*(crop2(i,j)-crop1(i,j))
            do k=1,11
              vdata(i,j,k) = vdata0(i,j,k)*(1.-crops)
            end do
              vdata(i,j,9) = crops
         end if
      end do
      end do
      if(reset_veg) call upd_gh

      year_old = year

      return
      end subroutine updveg