read_ndvi.F90

Go to the documentation of this file.

!
!-------------------------------------------------------------------------------
subroutine open_single_td( sib, time, filename)
!-------------------------------------------------------------------------------
! opens single point ascii file of time dependent parameters
! and reads in times for ndvi composite periods
!
! Modifications:
!   Kevin Schaefer created routine (8/1/05)
!-------------------------------------------------------------------------------
!
use kinds
use sibtype
use timetype
use sib_io_module
use sib_const_module
use sib_bc_module
!
implicit none
!
! input/output variables
type(sib_t), dimension(subcount), intent(inout) :: sib
type(time_struct) :: time
!
! local variables
integer(kind=int_kind) :: ntest1    ! compare nsib value of file to simulation
integer(kind=int_kind) :: i,j,k,l,m,n  ! indeces
character*100 filename  ! filename used to read in ndvi data
!
! close time dependent parameter file
    close(32)
!
! open time dependent parameter file
    open(32, file=trim(filename), form='formatted', status='old') 
!
! read number of sib points  
    read(32,*) ntest1
    if(ntest1/=nsib)  stop ' param file nsib no match with model nsib'
!
! read number composite periods
    read(32,*) nper
    if(nper>npermax)  stop ' too many composite periods in param file'
!
! read composite period start and stop times
    read(32,*)
    do i = 1,nper

!itb_modis
!      read(32,*) time%ndvi_start(i), time%ndvi_stop(i), j
      read(32,*) time%modis_start(i), time%modis_stop(i), j
    enddo
!
! read header line before time dep parameters
    read(32,*)
!
end subroutine open_single_td
!
!-------------------------------------------------------------------------------
subroutine read_single_td_param(sib, time)
!-------------------------------------------------------------------------------
! reads in time dependent parameters (ndvi, etc.) for single point from ascii file
!
! Modifications:
!   Kevin Schaefer created routine (8/1/05)
!-------------------------------------------------------------------------------
!
use kinds
use sibtype
use timetype
use sib_io_module
use sib_const_module
use sib_bc_module
!
implicit none
!
! input/output variables
type(sib_t), dimension(subcount), intent(inout) :: sib
type(time_struct), intent(in) :: time
!
! read in time dependent parameters (assume file at time%bc_recnum)
    read(32,*)
!itb_modis
!    read(32,*) sib%param%ndvi3
!    read(32,*) sib%param%NDVI_time3
    read(32,*) sib%param%mlai3
    read(32,*) sib%param%mfpar3
    read(32,*) sib%param%modis_time3
!itb_modis
    read(32,*) sib%param%d13cresp3
    read(32,*) sib%param%physfrac3(1)
    read(32,*) sib%param%physfrac3(2)
    read(32,*) sib%param%physfrac3(3)
    read(32,*) sib%param%physfrac3(4)
    read(32,*) sib%param%physfrac3(5)
    
!itb_modis
!    sib%param%ndvi_period3 = time%bc_recnum
    sib%param%modis_period3 = time%bc_recnum
!itb_modis

!
end subroutine read_single_td_param
!
!-------------------------------------------------------------------------------
subroutine open_global_td( sib, time, filename )
!-------------------------------------------------------------------------------
! opens netcdf file containing time dependent paramter data (ndvi, etc.)
! and reads in times for ndvi composite periods
!
! Modifications:
!  Kevin Schaefer created routine from new_bc (8/4/05)
!-------------------------------------------------------------------------------
!
use sibtype
use timetype
use sib_io_module
use sib_const_module
use sib_bc_module
!
#ifdef PGF
use netcdf
use typeSizes
#endif
!
implicit none
!
! input and output variables
type(sib_t), dimension(subcount), intent(inout) :: sib
type(time_struct) :: time
character*100 filename  ! filename used to read in ndvi data

! local variables
integer(kind=int_kind) :: ntest1   ! compare nsib value of file to simulation
integer(kind=int_kind) :: dimid    ! netcdf dimension id number
integer(kind=int_kind) :: var_id   ! netcdf variable ID number
integer(kind=int_kind) :: status   ! status of netcdf call
character(len=10) :: name          ! not sure
!
! close currently open file
    if(param_id/=0) then
      status = nf90_close ( param_id )
      if (status /= nf90_noerr) call handle_err (status)
    endif
!
! Open the parameter file

print*,'open_global_td: filename=',trim(filename)

    status = nf90_open ( trim(filename), nf90_nowrite, param_id )
    if (status /= nf90_noerr) call handle_err(status)
!
! check number of points
    status = nf90_inq_dimid ( param_id, 'nsib', dimid )
    if (status /= nf90_noerr) call handle_err(status)
    status = nf90_inquire_dimension ( param_id, dimid, name, ntest1 )
    if (status /= nf90_noerr) call handle_err(status)
    if(ntest1 /= nsib) stop ' open: file sib_bc no match with model for nsib'
!
! read in global attributes that are passed on to output files
!itb_modis...this is still called 'ndvi_source' in the file I have...
    status = nf90_get_att( param_id, nf90_global, 'ndvi_source', ndvi_source )
    status = nf90_get_att( param_id, nf90_global, 'c4_source', c4_source )
    status = nf90_get_att( param_id, nf90_global, 'd13cresp_source', d13cresp_source )



!itb_modis...


!
! get ndvi id number
!    status = nf90_inq_varid ( param_id, 'ndvi', ndvi_id )
!    if (status /= nf90_noerr) call handle_err(status)
!print*,'open_global_td: netcdf ids=',param_id,ndvi_id


!
! get ndvi time id number
!    status = nf90_inq_varid ( param_id, 'ndvi_time', ndvi_time_id )
!    if (status /= nf90_noerr) call handle_err(status)

!itb_modis; get modis id info
    status = nf90_inq_varid ( param_id, 'lai', mlai_id )
    if (status /= nf90_noerr) call handle_err(status)
    status = nf90_inq_varid ( param_id, 'fpar', mfpar_id )
    if (status /= nf90_noerr) call handle_err(status)

! get modis time id number
    status = nf90_inq_varid ( param_id, 'modis_time', modis_time_id )
    if (status /= nf90_noerr) call handle_err(status)

!print*,'open_global_td: netcdf ids=',param_id,mlai_id,mfpar_id

!itb_modis





!
! get d13cresp id number
    status = nf90_inq_varid ( param_id, 'd13cresp', d13_id )
    if (status /= nf90_noerr) call handle_err(status)
!
! get physfrac id number
    status = nf90_inq_varid ( param_id, 'physfrac', phys_id )
    if (status /= nf90_noerr) call handle_err(status)
!
! get number of composite periods
    status = nf90_inq_varid ( param_id, 'mapsyear', var_id )
    if (status /= nf90_noerr) call handle_err(status)
    status = nf90_get_var ( param_id, var_id, nper )
    if (status /= nf90_noerr) call handle_err (status)

!itb_modis
!
! get composite period start times
!    status = nf90_inq_varid ( param_id, 'ndvi_start', var_id )
!    if (status /= nf90_noerr) call handle_err(status)
!    status = nf90_get_var ( param_id, var_id, time%ndvi_start(1:nper) )
!    if (status /= nf90_noerr) call handle_err (status)
!
! get composite period stop times
!    status = nf90_inq_varid ( param_id, 'ndvi_stop', var_id )
!    if (status /= nf90_noerr) call handle_err(status)
!    status = nf90_get_var ( param_id, var_id, time%ndvi_stop(1:nper) )
!    if (status /= nf90_noerr) call handle_err (status)


! get composite period start times
    status = nf90_inq_varid ( param_id, 'modis_start', var_id )
    if (status /= nf90_noerr) call handle_err(status)
    status = nf90_get_var ( param_id, var_id, time%modis_start(1:nper) )
    if (status /= nf90_noerr) call handle_err (status)
!
! get composite period stop times
    status = nf90_inq_varid ( param_id, 'modis_stop', var_id )
    if (status /= nf90_noerr) call handle_err(status)
    status = nf90_get_var ( param_id, var_id, time%modis_stop(1:nper) )
    if (status /= nf90_noerr) call handle_err (status)

!itb_modis

!
end subroutine open_global_td
!
!-------------------------------------------------------------------------------
subroutine read_global_td_param(sib, time)
!-------------------------------------------------------------------------------
! Reads time dependent parameter data (ndvi, etc.) from netcdf file
!
! Modifications:
!  Kevin Schaefer created routine from read_ndvi (8/4/05)
!-------------------------------------------------------------------------------

use sibtype
use timetype
use sib_io_module
use sib_const_module
use sib_bc_module

#ifdef PGF
use netcdf
use typeSizes
#endif

Implicit none

! input/output variables
type(sib_t), dimension(subcount), intent(inout) :: sib
type(time_struct), intent(in) :: time
!
! local variables
integer(kind=int_kind) :: i,j
integer(kind=int_kind) :: ntest1     ! value of nsib read from sib_bc file
integer(kind=int_kind) :: status     ! netcdf status variable
integer(kind=int_kind) :: begin (2)  ! indeces where to start reading from file
integer(kind=int_kind) :: finish (2) ! indeces where to finish reading from file

!itb_modis
!real(kind=real_kind), dimension(nsib) :: ndvi      ! temp ndvi variable
!real(kind=real_kind), dimension(nsib) :: ndvi_time ! temp ndvi time variable

real(kind=real_kind), dimension(nsib) :: mlai       ! temp lai variable
real(kind=real_kind), dimension(nsib) :: mfpar       ! temp fpar variable
real(kind=real_kind), dimension(nsib) :: modis_time ! temp MODIS time variable

!itb_modis

real(kind=real_kind), dimension(nsib,physmax)  :: frac ! temp physfrac variabla
real(kind=real_kind), dimension(nsib)  :: d13  ! temp d12cresp variable
character(len=10) :: name
!
! set start and stop points of data to read
    begin (1)  = time%bc_recnum 
    begin (2)  = 1 
    finish (1) = 1 
    finish (2) = nsib 

!itb_modis
!
! get ndvi data
!print*,'read_ndvi:',param_id,ndvi_id,begin,finish
!    status = nf90_get_var ( param_id, ndvi_id, ndvi, begin, finish )
!    if (status /= nf90_noerr) call handle_err (status,'read_ndvi',1)
!
! get ndvi times
!    status = nf90_get_var ( param_id, ndvi_time_id, ndvi_time, begin, finish )
!    if (status /= nf90_noerr) call handle_err (status,'read_ndvi',2)

! get MODIS data
!print*,'read_modis:',param_id,mlai_id,mfpar_id,begin,finish
    status = nf90_get_var ( param_id, mlai_id, mlai, begin, finish )
    if (status /= nf90_noerr) call handle_err (status,'read_modis',1)

    status = nf90_get_var ( param_id, mfpar_id, mfpar, begin, finish )
    if (status /= nf90_noerr) call handle_err (status,'read_modis',2)
!
! get modis times
    status = nf90_get_var ( param_id, modis_time_id, modis_time, begin, finish )
    if (status /= nf90_noerr) call handle_err (status,'read_modis',3)


!itb_modis




!
! get d13cresp data
    status = nf90_get_var ( param_id, d13_id, d13 )
    if (status /= nf90_noerr) call handle_err (status,'read_ndvi',3)
!
! get physfrac data
    status = nf90_get_var ( param_id, phys_id, frac )
    if (status /= nf90_noerr) call handle_err (status,'read_ndvi',4)
!
! copy only points in subdomain
    do i = 1, subcount

!itb_modis
!        sib(i)%param%ndvi3 = ndvi(subset(i))
!        sib(i)%param%ndvi_time3 = ndvi_time(subset(i))
!        sib(i)%param%ndvi_period3 = time%bc_recnum

        sib(i)%param%mlai3 = mlai(subset(i))
        sib(i)%param%mfpar3 = mfpar(subset(i))
!print*,'READ_NDVI:',sib(i)%param%mlai3,sib(i)%param%mfpar3
        sib(i)%param%modis_time3 = modis_time(subset(i))
        sib(i)%param%modis_period3 = time%bc_recnum
!itb_modis

        
        sib(i)%param%d13cresp = d13(subset(i))
        do j = 1, physmax
            sib(i)%param%physfrac3(j) = frac(subset(i),j)
        enddo
    enddo
!
end subroutine read_global_td_param
!
!-------------------------------------------------------------------------------
subroutine calculate_td_param (sib, lat)
!-------------------------------------------------------------------------------
! Calls mapper and calculates the time dependent parameters
!
! MODIFICATIONS:
!  Owen Leonard created routine (8/10/2001)
!  Kevin Schaefer changed to single point (8/1/05)
!-------------------------------------------------------------------------------
!
use sibtype
use sib_const_module
use sib_bc_module
use kinds
!
implicit none
!
! input/output variables
type(sib_t), intent(inout) :: sib
    real(kind=real_kind) lat ! (deg) latitude of point
!
! local variables
integer(kind=int_kind) :: i,j,k,l,m,n   ! indeces
type time_dep_var
    real(kind=real_kind) :: fpar       ! canopy absorbed fraction of par
    real(kind=real_kind) :: lai        ! leaf-area index
    real(kind=real_kind) :: green      ! canopy greeness fraction of lai
    real(kind=real_kind) :: zo         ! canopy roughness coeff 
    real(kind=real_kind) :: zp_disp    ! zero plane displacement
    real(kind=real_kind) :: rbc        ! rb coefficient (c1)
    real(kind=real_kind) :: rdc        ! rc coefficient (c2)
    real(kind=real_kind) :: gmudmu     ! time-mean leaf projection
end type time_dep_var
type(time_dep_var) :: timevar
real(kind=real_kind) :: temptran (2,2)
real(kind=real_kind) :: tempref (2,2)
type(aero_var) :: tempaerovar(50,50)
!
! convert biome number to integer
    i = int(sib%param%biome)
!
! double to single precision conversion
    temptran = sib%param%tran(:,:)
    tempref = sib%param%ref(:,:)
    tempaerovar = aerovar(:,:,i)
!
! calculate new MODIS derived parameters
    call mapper(          &
      lat,&
      sib%param%modis_time1,           &      
      sib%param%modis_time2,           &
      sib%param%mlai1,      &
      sib%param%mlai2,       &
      sib%param%mfpar1,      &
      sib%param%mfpar2,       &
      sib%param%vcover, &
      sib%param%chil,   &
      temptran,         &
      tempref,          &
      morphtab(i),      &
      tempaerovar,      &
      laigrid,          &
      fvcovergrid,      &
      timevar)
!
! move ndvi parameters to sib variable tree
    sib%param%aparc2 = timevar%fpar
    sib%param%zlt2 = timevar%lai
    sib%param%green2 = timevar%green
    sib%param%z0d2 = timevar%zo
    sib%param%zp_disp2 = timevar%zp_disp
    sib%param%rbc2 = timevar%rbc
    sib%param%rdc2 = timevar%rdc
    sib%param%gmudmu2 = timevar%gmudmu 

end subroutine calculate_td_param
!
!------------------------------------------------------------------------------
subroutine need_to_switch (sib,time)
!------------------------------------------------------------------------------
! determines if time dependent parameters need to be switched ( 3-->2, 2-->1 )
!
! MODIFICATIONS:
!  Andrew Philpott created routine (9/28/05)
!------------------------------------------------------------------------------
!
use sibtype
use timetype
use sib_const_module
use sib_bc_module
use kinds
!
implicit none
!
! parameters
type(sib_t), dimension(subcount), intent(inout) :: sib
type(time_struct), intent(in) :: time
!
! local variables
integer i      ! index
real(kind=real_kind) :: ndtime   ! (day) local version of ndvi2 time
integer(kind=int_kind) :: ndp    ! (#) local version of ndvi2 period
real(kind=real_kind) :: conv  ! (day) local real copy of time%sec_per_day
!
conv = real(time%sec_per_day)
do i=1,subcount
!
        ! copy local versions of ndvi2 time and ndvi2 period
!itb_modis
!        ndtime = sib(i)%param%ndvi_time2
        ndtime = sib(i)%param%modis_time2

!        ndp = sib(i)%param%ndvi_period2
        ndp = sib(i)%param%modis_period2
!
        ! two cases: Case A) -- ndvi2 time is at the end of ndvi2 period
        !            in this case we want to switch one timestep prior to 
        !            then because we want to copy the old period 3 values to 
        !            period 2 before overwriting them in the read routines
        !  Note -- checking that two real numbers are equal requires rather 
        !       checking that they are very nearly equal (avoid roundoff diffs)


!itb_modis; I'm just replacing 'ndvi' with 'modis' here, hoping the 
!           conversion is seamless...

!itb        if ( ndtime >= time%ndvi_stop(ndp) - 0.1*time%dtsib/conv .and. &
!itb               ndtime <= time%ndvi_stop(ndp) + 0.1*time%dtsib/conv ) then
!itb            if ( time%real_doy >= ndtime - 1.1*time%dtsib/conv .and. &
!itb                   time%real_doy <= ndtime - 0.9*time%dtsib/conv .and. &
!itb                   ndp<nper   ) then
!itb              call switch_td_param (sib(i))
!itb              call calculate_td_param (sib(i), latsib(subset(i)))
!       print*,'A1 switch',sib%param%ndvi_time1,sib%param%ndvi_time2,sib%param%ndvi_time3
!itb            else if ( time%real_doy >= 365. - 1.1*time%dtsib/conv .and. &
!itb                      time%real_doy <= 365. - 0.9*time%dtsib/conv .and. &
!itb                      ndp==nper   ) then
!itb              call switch_td_param (sib(i))
!itb              call calculate_td_param (sib(i), latsib(subset(i)))
!       print*,'A2 switch',sib%param%ndvi_time1,sib%param%ndvi_time2,sib%param%ndvi_time3                     
!itb            end if
        ! Case B) -- ndvi2 time is not at the end of ndvi2 period
        !            switch td params at the beginning of the first day 
        !            that falls after ndvi2 time                
!itb        else if (time%new_day) then
!itb          if (time%real_doy>=ndtime .and. &
!itb              (ndp>1.or.time%real_doy<time%ndvi_start(nper))) then
!itb            call switch_td_param (sib(i))
!itb            call calculate_td_param (sib(i), latsib(subset(i)))
!       print*,'B1 switch',sib%param%ndvi_time1,sib%param%ndvi_time2,sib%param%ndvi_time3
!itb          end if
!itb    end if

!print*,'need_to_switch:',ndtime,ndp
        if ( ndtime >= time%modis_stop(ndp) - 0.1*time%dtsib/conv .and. &
               ndtime <= time%modis_stop(ndp) + 0.1*time%dtsib/conv ) then
            if ( time%real_doy >= ndtime - 1.1*time%dtsib/conv .and. &
                   time%real_doy <= ndtime - 0.9*time%dtsib/conv .and. &
                   ndp<nper   ) then
              call switch_td_param (sib(i))
              call calculate_td_param (sib(i), latsib(subset(i)))
!       print*,'A1 switch',sib%param%modis_time1,sib%param%modis_time2,sib%param%modis_time3
            else if ( time%real_doy >= 365. - 1.1*time%dtsib/conv .and. &
                      time%real_doy <= 365. - 0.9*time%dtsib/conv .and. &
                      ndp==nper   ) then
              call switch_td_param (sib(i))
              call calculate_td_param (sib(i), latsib(subset(i)))
!       print*,'A2 switch',sib%param%modis_time1,sib%param%modis_time2,sib%param%modis_time3                  
            end if
        ! Case B) -- ndvi2 time is not at the end of ndvi2 period
        !            switch td params at the beginning of the first day 
        !            that falls after ndvi2 time                
        else if (time%new_day) then
          if (time%real_doy>=ndtime .and. &
              (ndp>1.or.time%real_doy<time%modis_start(nper))) then
            call switch_td_param (sib(i))
            call calculate_td_param (sib(i), latsib(subset(i)))
!       print*,'B1 switch',sib%param%modis_time1,sib%param%modis_time2,sib%param%modis_time3
          end if
        end if
!
end do
!
end subroutine need_to_switch                                   
                                    
!-------------------------------------------------------------------------------
subroutine switch_td_param (sib)
!-------------------------------------------------------------------------------
! switches time dependent parameters
!
! MODIFICATIONS:
!  Kevin Schaefer created routine (8/1/05)
!-------------------------------------------------------------------------------
!
use sibtype
use sib_const_module
use sib_bc_module
use kinds
!
implicit none
!
! input/output variables
type(sib_t), intent(inout) :: sib
!
! local variables
integer(kind=int_kind) :: i,j,k,l,m,n   ! indeces
!
! switch time dependent parameters
!itb_modis

!        sib%param%ndvi1        = sib%param%ndvi2
!        sib%param%ndvi_time1   = sib%param%ndvi_time2

        sib%param%mlai1        = sib%param%mlai2
        sib%param%mfpar1        = sib%param%mfpar2
        sib%param%modis_time1   = sib%param%modis_time2

!        sib%param%ndvi_period1 = sib%param%ndvi_period2        
        sib%param%modis_period1 = sib%param%modis_period2
!itb_modis
        
        sib%param%aparc1       = sib%param%aparc2
        sib%param%zlt1         = sib%param%zlt2
        sib%param%green1       = sib%param%green2
        sib%param%z0d1         = sib%param%z0d2
        sib%param%zp_disp1     = sib%param%zp_disp2
        sib%param%rbc1         = sib%param%rbc2
        sib%param%rdc1         = sib%param%rdc2
        sib%param%gmudmu1      = sib%param%gmudmu2
        sib%param%d13cresp1    = sib%param%d13cresp2
        do i = 1, physmax
            sib%param%physfrac1(i) = sib%param%physfrac2(i)
        enddo
!
!itb_modis...
! move in new ndvi data
!        sib%param%ndvi2       = sib%param%ndvi3
!        sib%param%ndvi_time2   = sib%param%ndvi_time3
!        sib%param%ndvi_period2 = sib%param%ndvi_period3        
        sib%param%mlai2       = sib%param%mlai3
        sib%param%mfpar2       = sib%param%mfpar3
        sib%param%modis_time2   = sib%param%modis_time3
        sib%param%modis_period2 = sib%param%modis_period3
!itb_modis

        
        sib%param%d13cresp2    = sib%param%d13cresp3
        do i = 1, physmax
            sib%param%physfrac2(i) = sib%param%physfrac3(i)
        enddo 

!print*,'switch:',sib%param%ndvi_time1,sib%param%ndvi_time2,sib%param%ndvi_time3

end subroutine switch_td_param
!