ecearth3/sources/lpjg/offline/lpjg_forcing_ifs.F90

program LPJGFORC

! reading from NetCDF files and passing the data to OASIS
! units for the climatic parameters:
! temperature (K)
! precipitation (m s-1)
! radiation (W m-2)

! Previous update: 10 Dec 2015 - add 2 new fields for low veg. - see "ecev3"
! Now reads AMIP output.
! update: Nov. 2017 (Etienne Tourigny) - update indentation, add namelist for several parameters
!   and sync with latest changes in ESM branch (from r4521)
! update: Feb. 2020 (Etienne Tourigny) - dummy tm5 coupling - co2 ppm + fluxes, need to set co2 ppm value in namelist to activate this (284 for PI), or set to 0 to deactivate
! update: May  2020 (Etienne Tourigny) - moved tm5 dummy coupling to ccycle_coupling

  use mod_oasis
  use mpi

  implicit none

!  integer,parameter :: MAXGRID = 35718
!  integer,parameter :: MAXGRID = 88838
#if defined(IFS_RES_T159)
  integer,parameter :: MAXGRID = 35718 ! T159
#elif defined(IFS_RES_T255)
  integer,parameter :: MAXGRID = 88838 ! T255
#else
  integer,parameter :: MAXGRID = 88838 ! T255
#endif
  integer,parameter :: t_max_year = 365
  integer,parameter :: t_max_leap = 366

!    integer,parameter :: num_years = 10
!    integer,parameter :: start_year = 1850
!    integer,parameter :: num_loops = 10
  integer :: num_years,start_year,num_loops
!  character(len=*),parameter::nc_file_path="./ifs_data/"
  character(len=1024) :: nc_file_path
  character(len=*),parameter::nc_file_tail="_dayavg.nc"
  namelist /NAMLPJGFORC/ num_years,start_year,num_loops,nc_file_path
  character(len=*),parameter :: NamelistFileName  = "namelist.lpjg_forcing"
  integer,parameter :: NAMLPJGFORCU = 999
  logical :: LPJGFORCnotFree

  character(len=6),parameter :: model_name    = "LPJGFORC"

  double precision :: minhlai,maxhlai,minllai,maxllai
  double precision :: meanhlai,meanllai,meanvegl,meanvegh
  double precision :: maxcvh,maxcvl,maxtvl,maxtvh,cover

  ! PM, was SPATEMPE, SPAPRECI, SPARADIA, SPATOLAI

  character(len=128),parameter :: f_temp_name = "T2MVeg"
  integer                    :: f_temp_id   = -77
  double precision           :: f_temp_data(MAXGRID,1)

  character(len=128),parameter :: f_prec_name = "TPVeg"
  integer                    :: f_prec_id   = -77
  double precision           :: f_prec_data(MAXGRID,1)
  double precision           :: f_precs_data(MAXGRID,1)
  double precision           :: f_precc_data(MAXGRID,1)
  double precision           :: f_prect_data(MAXGRID,1)

!    character(len=128),parameter :: f_vegl_name = "CVLVeg"
!    integer                    :: f_vegl_id   = -77
!    double precision           :: f_vegl_data(MAXGRID,1)

!    character(len=128),parameter :: f_vegltype_name = "TVLVeg"
!    integer                    :: f_vegltype_id   = -77
!    double precision           :: f_vegltype_data(MAXGRID,1)

!    character(len=128),parameter :: f_vegh_name = "CVHVeg"
!    integer                    :: f_vegh_id   = -77
!    double precision           :: f_vegh_data(MAXGRID,1)

  character(len=128),parameter :: f_snoc_name = "SDVeg"
  integer                    :: f_snoc_id   = -77
  double precision           :: f_snoc_data(MAXGRID,1)

  character(len=128),parameter :: f_snod_name = "SDensVeg"
  integer                    :: f_snod_id   = -77
  double precision           :: f_snod_data(MAXGRID,1)

  character(len=128),parameter :: f_st1l_name = "SoilTVeg.L001"
  integer                    :: f_st1l_id   = -77
  double precision           :: f_st1l_data(MAXGRID,1)

  character(len=128),parameter :: f_st2l_name = "SoilTVeg.L002"
  integer                    :: f_st2l_id   = -77
  double precision           :: f_st2l_data(MAXGRID,1)

  character(len=128),parameter :: f_st3l_name = "SoilTVeg.L003"
  integer                    :: f_st3l_id   = -77
  double precision           :: f_st3l_data(MAXGRID,1)

  character(len=128),parameter :: f_st4l_name = "SoilTVeg.L004"
  integer                    :: f_st4l_id   = -77
  double precision           :: f_st4l_data(MAXGRID,1)

  character(len=128),parameter :: f_sm1l_name = "SoilMVeg.L001"
  integer                    :: f_sm1l_id   = -77
  double precision           :: f_sm1l_data(MAXGRID,1)

  character(len=128),parameter :: f_sm2l_name = "SoilMVeg.L002"
  integer                    :: f_sm2l_id   = -77
  double precision           :: f_sm2l_data(MAXGRID,1)

  character(len=128),parameter :: f_sm3l_name = "SoilMVeg.L003"
  integer                    :: f_sm3l_id   = -77
  double precision           :: f_sm3l_data(MAXGRID,1)

  character(len=128),parameter :: f_sm4l_name = "SoilMVeg.L004"
  integer                    :: f_sm4l_id   = -77
  double precision           :: f_sm4l_data(MAXGRID,1)

  character(len=128),parameter :: f_rads_name = "SSRVeg"
  integer                    :: f_rads_id   = -77
  double precision           :: f_rads_data(MAXGRID,1)

  character(len=128),parameter :: f_radl_name = "SLRVeg"
  integer                    :: f_radl_id   = -77
  double precision           :: f_radl_data(MAXGRID,1)

  character(len=128),parameter :: f_llai_name = "LAILVeg"
  integer                    :: f_llai_id   = -77
  double precision           :: f_llai_data(MAXGRID,1)

  character(len=128),parameter :: f_hlai_name = "LAIHVeg"
  integer                    :: f_hlai_id   = -77
  double precision           :: f_hlai_data(MAXGRID,1)

  character(len=128),parameter :: f_typh_name = "TypeHVeg"
  integer                    :: f_typh_id   = -77
  double precision           :: f_typh_data(MAXGRID,1)	

  character(len=128),parameter :: f_vggh_name = "FracHVeg"
  integer                    :: f_vggh_id   = -77
  double precision           :: f_vggh_data(MAXGRID,1)

! ecev3
  character(len=128),parameter :: f_typl_name = "TypeLVeg"
  integer                    :: f_typl_id   = -77
  double precision           :: f_typl_data(MAXGRID,1)	

! ecev3
  character(len=128),parameter :: f_vggl_name = "FracLVeg"
  integer                    :: f_vggl_id   = -77
  double precision           :: f_vggl_data(MAXGRID,1)

  integer :: comp_id = -77
  integer :: part_id = -77
  integer :: t,t_step,yr,lyr,t_max,t_step_full,cell,loopyr
  integer :: ierror

  integer::ix,iy

! VARIABLE FIELDS
  character(len=256) :: nc_temp_file
  character(len=256) :: nc_precs_file
  character(len=256) :: nc_precc_file
  character(len=256) :: nc_prect_file
  character(len=256) :: nc_snoc_file
  character(len=256) :: nc_snod_file
  character(len=256) :: nc_st1l_file
  character(len=256) :: nc_st2l_file
  character(len=256) :: nc_st3l_file
  character(len=256) :: nc_st4l_file
  character(len=256) :: nc_sm1l_file
  character(len=256) :: nc_sm2l_file
  character(len=256) :: nc_sm3l_file
  character(len=256) :: nc_sm4l_file
  character(len=256) :: nc_rads_file
  character(len=256) :: nc_radl_file

! FIXED FIELDS - T159!!!!

!      character(len=*),parameter::nc_vegltype_file="/nobackup/rossby15/sm_paumi/preind/tvl.nc"
!      character(len=*),parameter::nc_veghtype_file="/nobackup/rossby15/sm_paumi/preind/tvh.nc"

  logical :: isleapyear = .false.
  character yearstr*4

! VARIABLE NAMES
  character(len=*),parameter::nc_temp_name="var167"
  character(len=*),parameter::nc_precs_name="var142"
  character(len=*),parameter::nc_precc_name="var143"
  character(len=*),parameter::nc_prect_name="var228"
  ! We send fixed values to guess below if do not have the SNOC and SNOD data
  character(len=*),parameter::nc_snoc_name="var141"
  character(len=*),parameter::nc_snod_name="var033"
  character(len=*),parameter::nc_snod_varname="var33"
  character(len=*),parameter::nc_st1l_name="var139"
  character(len=*),parameter::nc_st2l_name="var170"
  ! ST3L and ST4L - we do not have this data, so we'll send ST2L values instead - see below.
  ! TODO use data from IFS which is now available
  character(len=*),parameter::nc_st3l_name="var183"
  character(len=*),parameter::nc_st4l_name="var236"
  character(len=*),parameter::nc_sm1l_name="var039" ! AMIP files have this
  character(len=*),parameter::nc_sm2l_name="var040"
  character(len=*),parameter::nc_sm3l_name="var041"
  character(len=*),parameter::nc_sm4l_name="var042"
  character(len=*),parameter::nc_sm1l_varname="var39"
  character(len=*),parameter::nc_sm2l_varname="var40"
  character(len=*),parameter::nc_sm3l_varname="var41"
  character(len=*),parameter::nc_sm4l_varname="var42"
  character(len=*),parameter::nc_rads_name="var176"
  character(len=*),parameter::nc_radl_name="var177"
!    character(len=*),parameter::nc_vegl_name="tvl"
!    character(len=*),parameter::nc_vegh_name="tvh"
  logical,parameter::ignore_missing_stl=.true. ! ignore missing files for st1l and st4l vars
  logical,parameter::ignore_missing_sml=.true. ! ignore missing files for sm?l vars
  logical,parameter::ignore_missing_sno=.true. ! ignore missing files for snoc and snod
  logical,parameter::convert_fluxes=.false. ! should we convert fluxes to daily accumulated values?

  logical::file_exists,file2_exists,file3_exists

! FIELD IDS
  integer::nc_temp_fileid,nc_precs_fileid,nc_rads_fileid,nc_radl_fileid,nc_vegl_fileid,nc_vegh_fileid,nc_precc_fileid,nc_prect_fileid;
  integer::nc_temp_varid,nc_precs_varid,nc_rads_varid,nc_radl_varid,nc_vegl_varid,nc_vegh_varid,nc_precc_varid,nc_prect_varid;
!    integer::nc_vegltype_varid; ! We don't have a T255 file for vegltype yet

  integer::nc_snoc_fileid,nc_snod_fileid;
  integer::nc_snoc_varid,nc_snod_varid;
  integer::nc_st1l_fileid,nc_st2l_fileid,nc_st3l_fileid,nc_st4l_fileid;
  integer::nc_st1l_varid,nc_st2l_varid,nc_st3l_varid,nc_st4l_varid;

  integer::nc_sm1l_fileid,nc_sm2l_fileid,nc_sm3l_fileid,nc_sm4l_fileid;
  integer::nc_sm1l_varid,nc_sm2l_varid,nc_sm3l_varid,nc_sm4l_varid;

  integer::localcomm, cplcomm, icpl;

! *** START ***

  write (*,'(A)') "*II* lpjg_forcing: Hello"
  write (*,*) "*II* lpjg_forcing: MAXGRID=",MAXGRID

  ! Read RunLengthSec,TimeStepSec from namelist
  inquire(NAMLPJGFORCU,opened=LPJGFORCnotFree)
  if (LPJGFORCnotFree) call ERROR('Namelist LPJGFORC not free',1)
  open(NAMLPJGFORCU,file=NamelistFileName,status='OLD')
  read(NAMLPJGFORCU,nml=NAMLPJGFORC)
  close(NAMLPJGFORCU)

  write (*,'(A)') "*II* lpjg_forcing: Now initialising Lpjg_forcing using oasis_..."

  call oasis_init_comp(comp_id,model_name,ierror)
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_init_comp returned  ierror=",ierror

!    call oasis_get_localcomm(localcomm, ierror)
!    write (*,'(A,I3)') "*II* lpjg_forcing: get_localcomm  ierror=",ierror
!    write (*,'(A,I12)') "*II* lpjg_forcing: oasis_get_localcomm returned localcomm =",localcomm

  icpl = 1
  call oasis_create_couplcomm(icpl,MPI_COMM_SELF, cplcomm, ierror)
!    call oasis_create_couplcomm(icpl,localcomm, cplcomm, ierror)
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_create_couplcomm  ierror=",ierror
  write (*,'(A,I12)') "*II* lpjg_forcing: oasis_create_couplcomm returned cplcomm =",cplcomm

  call oasis_def_partition(part_id,(/ 0,0,MAXGRID*1 /),ierror)
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_partition returned part_id =",part_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_partition returned  ierror =",ierror

  call oasis_def_var( f_temp_id, f_temp_name,     &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_temp_id =",f_temp_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_prec_id, f_prec_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_prec_id =",f_prec_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

!    call oasis_def_var( f_vegl_id, f_vegl_name,   &
!       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
!    write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_vegl_id =",f_vegl_id
!    write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

!    call oasis_def_var( f_vegh_id, f_vegh_name,   &
!       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
!    write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_vegh_id =",f_vegh_id
!    write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

!   call oasis_def_var( f_vegltype_id, f_vegltype_name,   &
!       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
!    write (*,'(A,I3)') "*II* lpjg_forcing: f_vegltype_id =",f_vegltype_id
!    write (*,'(A,I3)') "*II* lpjg_forcing: returned ierror=",ierror

  call oasis_def_var( f_snoc_id, f_snoc_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_snoc_id =",f_snoc_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_snod_id, f_snod_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_snod_id =",f_snod_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_st1l_id, f_st1l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_st1l_id =",f_st1l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_st2l_id, f_st2l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_st2l_id =",f_st2l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_st3l_id, f_st3l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_st3l_id =",f_st3l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_st4l_id, f_st4l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_st4l_id =",f_st4l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_sm1l_id, f_sm1l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_sm1l_id =",f_sm1l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_sm2l_id, f_sm2l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_sm2l_id =",f_sm2l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_sm3l_id, f_sm3l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_sm3l_id =",f_sm3l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_sm4l_id, f_sm4l_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_sm4l_id =",f_sm4l_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_rads_id, f_rads_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_rads_id =",f_rads_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  call oasis_def_var( f_radl_id, f_radl_name,   &
       part_id, (/ 2,1 /), PRISM_Out, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_radl_id =",f_radl_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

! LLAI
  call oasis_def_var( f_llai_id, f_llai_name,   &
       part_id, (/ 2,1 /), PRISM_In, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_llai_id =",f_llai_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

! HLAI
  call oasis_def_var( f_hlai_id, f_hlai_name,   &
       part_id, (/ 2,1 /), PRISM_In, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_hlai_id =",f_hlai_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

! TYPH
  call oasis_def_var( f_typh_id, f_typh_name,   &
       part_id, (/ 2,1 /), PRISM_In, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_typh_id =",f_typh_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

! VEGH
  call oasis_def_var( f_vggh_id, f_vggh_name,   &
       part_id, (/ 2,1 /), PRISM_In, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_vggh_id =",f_vggh_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

! ecev3 TYPL
  call oasis_def_var( f_typl_id, f_typl_name,   &
       part_id, (/ 2,1 /), PRISM_In, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_typl_id =",f_typl_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

! ecev3 VEGL
  call oasis_def_var( f_vggl_id, f_vggl_name,   &
       part_id, (/ 2,1 /), PRISM_In, (/ 1,MAXGRID,1,1 /), PRISM_Real, ierror )
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned f_vggl_id =",f_vggl_id
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_def_var returned    ierror =",ierror

  write (*,'(A)') "*II* lpjg_forcing: before call oasis_enddef(ierror)"
  call oasis_enddef(ierror)
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_enddef returned ierror =",ierror

! FIXED FIELDS

!    write (*,'(A)') "*II* lpjg_forcing: Reading vegl and vegh"

! VEGL
  ! call read_fixedfiels(MAXGRID,1,nc_vegl_file,nc_vegl_name,&
  ! f_vegl_data,nc_vegl_fileid,nc_vegl_varid)

! VEGH
  ! call read_fixedfiels(MAXGRID,1,nc_vegh_file,nc_vegh_name,&
  ! f_vegh_data,nc_vegh_fileid,nc_vegh_varid)

  write (*,'(A)') "*II* lpjg_forcing: Beginning time loop"

  t_step_full = 0

! Repetitions of forcing
  do lyr = 1, num_loops

     ! TODO fix various write call format, e.g. this one does not print lyr
     write (*,*) "*III* lpjg_forcing: Beginning loop #",lyr

! Year loop
     do loopyr = start_year, start_year+num_years-1

        write (*,*) "*III* lpjg_forcing: Beginning year",loopyr

        yr = loopyr
        ! ET disabled May 2019
        !if (loopyr < 1870) yr = loopyr + 20
        !if (loopyr < 1870 .and. loopyr > 1859) yr = loopyr + 10

! Leap year?
        isleapyear = .false.
        if (mod(yr,4) .eq. 0)   isleapyear = .true.
        if (mod(yr,100) .eq. 0) isleapyear = .false.
        if (mod(yr,400) .eq. 0) isleapyear = .true.

! Create file names

        write(unit=yearstr, fmt='(I4)') yr

        write (*,'(A)') "netcdf files to read:"

        nc_temp_file = trim(nc_file_path)//"/"//nc_temp_name//"_"//yearstr//nc_file_tail
        write (*,'(A)') nc_temp_file

        nc_precs_file = trim(nc_file_path)//"/"//nc_precs_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_precs_file

        nc_precc_file = trim(nc_file_path)//"/"//nc_precc_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_precc_file

        nc_prect_file = trim(nc_file_path)//"/"//nc_prect_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_prect_file

        nc_snoc_file=trim(nc_file_path)//"/"//nc_snoc_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_snoc_file

        nc_snod_file=trim(nc_file_path)//"/"//nc_snod_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_snod_file

        nc_st1l_file=trim(nc_file_path)//"/"//nc_st1l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_st1l_file

        nc_st2l_file=trim(nc_file_path)//"/"//nc_st2l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_st2l_file

        nc_st3l_file=trim(nc_file_path)//"/"//nc_st3l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_st3l_file

        nc_st4l_file=trim(nc_file_path)//"/"//nc_st4l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_st4l_file

        nc_sm1l_file=trim(nc_file_path)//"/"//nc_sm1l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_sm1l_file

        nc_sm2l_file = trim(nc_file_path)//"/"//nc_sm2l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_sm2l_file

        nc_sm3l_file = trim(nc_file_path)//"/"//nc_sm3l_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_sm3l_file

        nc_sm4l_file = trim(nc_file_path)//"/"//nc_sm4l_name //"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_sm4l_file

        nc_rads_file = trim(nc_file_path)//"/"//nc_rads_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_rads_file

        nc_radl_file = trim(nc_file_path)//"/"//nc_radl_name//"_"//yearstr//trim(nc_file_tail)
        write (*,'(A)') nc_radl_file

        if (isleapyear .eq. .true.) then
           t_max = t_max_leap
        else
           t_max = t_max_year
        endif

        ! Day/6-hr loop for each year

        write (*,'(A,I5)') "*II* lpjg_forcing: t_max =",t_max

        do t_step=0,t_max-1
           write (*,*) "*III* lpjg_forcing: Beginning dayloop ",t_step

           t = t_step_full*86400

           ! reading climate data from NetCDF file

           ! PM - changed to 6,12 from 2,6
           write (*,'(A,I6,A,I12)') "*II* lpjg_forcing before read_ifstest: Time step t = ",t_step," - time t=",t

           ! PM to _pi version
           call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_temp_file,nc_temp_name,&
                f_temp_data,nc_temp_fileid,nc_temp_varid)

           ! liquid precipitaion can come from either 142 and 143 (from IFS) or 228 (from OSM)
           ! snowfall is not used anymore since it is included in other vars
           ! TODO check real ifs experiment with 142 143 and 228
           inquire(file=nc_precs_file, exist=file_exists)
           inquire(file=nc_precc_file, exist=file2_exists)
           inquire(file=nc_prect_file, exist=file3_exists)

           if (file_exists .and. file2_exists) then ! precs+precc
              call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_precs_file,nc_precs_name,&
                   f_precs_data,nc_precs_fileid,nc_precs_varid)
              call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_precc_file,nc_precc_name,&
                   f_precc_data,nc_precc_fileid,nc_precc_varid)
              f_prect_data(:,:) = 0
              f_prec_data = f_precs_data + f_precc_data
           else if (file3_exists) then ! prect
              f_precs_data(:,:) = 0
              f_precc_data(:,:) = 0
              call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_prect_file,nc_prect_name,&
                   f_prect_data,nc_prect_fileid,nc_prect_varid)
              f_prec_data = f_prect_data
           else
              call ERROR('ERROR! cannot file forcing file(s) for precipitation, looked for '&
                   //trim(nc_precs_file)//' '//trim(nc_precc_file)//' '//trim(nc_prect_file), 1)
           end if

           ! PM - changed to 6,12 from 2,6
           write (*,'(A,I6,A,I12)') "*II* lpjg_forcing after precip: Time step t = ",t_step," - time t=",t

           ! soil layers 1 and 4 are unused in LPJG, so send soil layers 2 and 3 if the input files are not found
           call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_st2l_file,nc_st2l_name,&
                f_st2l_data,nc_st2l_fileid,nc_st2l_varid,.false.)

           call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_st3l_file,nc_st3l_name,&
                f_st3l_data,nc_st3l_fileid,nc_st3l_varid,.false.)

           ! snow
           inquire(file=nc_snoc_file, exist=file_exists)
           if (file_exists .or. (.not.ignore_missing_sno) ) then
              call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_snoc_file,nc_snoc_name,&
                   f_snoc_data,nc_snoc_fileid,nc_snoc_varid,ignore_missing_sno)
           else
              if (t_step .eq. 0) write (*,*) "*II* lpjg_forcing: input file ",trim(nc_snoc_file)," is not found, setting snoc values to 0.0"
              f_snoc_data(:,:) = 0.0
           endif
           inquire(file=nc_snod_file, exist=file_exists)
           if (file_exists .or. (.not.ignore_missing_sno) ) then
              call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_snod_file,nc_snod_varname,&
                   f_snod_data,nc_snod_fileid,nc_snod_varid,ignore_missing_sno)
           else
              if (t_step .eq. 0) write (*,*) "*II* lpjg_forcing: input file ",trim(nc_snod_file)," is not found, setting snod values to 330.0"
              f_snod_data(:,:) = 330.0
           endif

           inquire(file=nc_st1l_file, exist=file_exists)
           if (file_exists .or. (.not.ignore_missing_stl) ) then
              call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_st1l_file,nc_st1l_name,&
                   f_st1l_data,nc_st1l_fileid,nc_st1l_varid,ignore_missing_stl)
           else
              if (t_step .eq. 0) write (*,*) "*II* lpjg_forcing: input file ",trim(nc_st1l_file)," is not found, using data from ",trim(nc_st2l_file)
              f_st1l_data(:,:) = f_st2l_data(:,:)
           endif

           inquire(file=nc_st4l_file, exist=file_exists)
           if (file_exists .or. (.not.ignore_missing_stl) ) then
              call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_st4l_file,nc_st4l_name,&
                   f_st4l_data,nc_st4l_fileid,nc_st4l_varid,ignore_missing_stl)
           else
              if (t_step .eq. 0) write (*,*) "*II* lpjg_forcing: input file ",trim(nc_st4l_file)," is not found, using data from ",trim(nc_st3l_file)
              f_st4l_data(:,:) = f_st3l_data(:,:)
           endif

           call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_sm1l_file,nc_sm1l_varname,&
                f_sm1l_data,nc_sm1l_fileid,nc_sm1l_varid,ignore_missing_sml)

           call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_sm2l_file,nc_sm2l_varname,&
                f_sm2l_data,nc_sm2l_fileid,nc_sm2l_varid,ignore_missing_sml)

           call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_sm3l_file,nc_sm3l_varname,&
                f_sm3l_data,nc_sm3l_fileid,nc_sm3l_varid,ignore_missing_sml)

           call read_ifstest_pi_depth(MAXGRID,1,t_step,t_max,nc_sm4l_file,nc_sm4l_varname,&
                f_sm4l_data,nc_sm4l_fileid,nc_sm4l_varid,ignore_missing_sml)

           ! PM - changed to 6,12 from 2,6
           write (*,'(A,I6,A,I12)') "*II* lpjg_forcing after SM files: Time step t = ",t_step," - time t=",t

           call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_rads_file,nc_rads_name,&
                f_rads_data,nc_rads_fileid,nc_rads_varid)

           call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_radl_file,nc_radl_name,&
                f_radl_data,nc_radl_fileid,nc_radl_varid)

           ! convert fluxes to daily accumulated values? This assumes 6h post-processing timestep
           if ( convert_fluxes ) then
           do cell = 1, MAXGRID
              ! Total precip
              ! AMIP data have units m/6h average, so *4 *
              ! 1000 to get mm/day
              ! Now * 1000 and / 1000 to go to kg m-2 day-1,
              ! and / 86400.0 to get to kg m-2 s-1
              ! PM_Apr2012 - convert m to kg m-2 s-1 to mimic ECE - we'll convert back in LPJG
              !f_prec_data(cell,1) = 4 * 1000.0 *(f_precs_data(cell,1) +  f_precc_data(cell,1) +  f_prect_data(cell,1)
              f_prec_data(cell,1)=f_prec_data(cell,1) * 4 * 1000.0 / 86400.0
              ! PM_Apr2012 - change J m-2 to W m-2 to be consistent with ECE. We'll remove the / in LPJG.
              f_rads_data(cell,1) = 4* f_rads_data(cell,1) / 86400.0
              f_radl_data(cell,1) = 4* f_radl_data(cell,1) / 86400.0
           enddo
           endif

!                        do cell = 1, MAXGRID
!                                ! Enforce vegltype as grassland
!                                f_vegltype_data(cell,1) = 2;
!                        enddo

           ! PM - changed to 6,12 from 2,6
           write (*,'(A,I6,A,I12)') "*II* lpjg_forcing after ifstest: Time step t = ",t_step," - time t=",t

           ! -----------------------------------------------------------------------------------
           ! *** PUT variables
           ! -----------------------------------------------------------------------------------

           write (*,'(A,I6,A,I12)') "*II* lpjg_forcingPUT : Time step t = ",t_step," - time t=",t

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_temp_id =",f_temp_id
           call oasis_put(f_temp_id,t,f_temp_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_rads_id =",f_rads_id
           call oasis_put(f_rads_id,t,f_rads_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_radl_id =",f_radl_id
           call oasis_put(f_radl_id,t,f_radl_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_sm1l_id =",f_sm1l_id
           call oasis_put(f_sm1l_id,t,f_sm1l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_sm2l_id =",f_sm2l_id
           call oasis_put(f_sm2l_id,t,f_sm2l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_sm3l_id =",f_sm3l_id
           call oasis_put(f_sm3l_id,t,f_sm3l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_sm4l_id =",f_sm4l_id
           call oasis_put(f_sm4l_id,t,f_sm4l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_st1l_id =",f_st1l_id
           call oasis_put(f_st1l_id,t,f_st1l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_st2l_id =",f_st2l_id
           call oasis_put(f_st2l_id,t,f_st2l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_st3l_id =",f_st3l_id
           call oasis_put(f_st3l_id,t,f_st3l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_st4l_id =",f_st4l_id
           call oasis_put(f_st4l_id,t,f_st4l_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           ! Snow mass/unit surface (kg/m2)
           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_snoc_id =",f_snoc_id
           call oasis_put(f_snoc_id,t,f_snoc_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           ! Snow density (kg/m3)
           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_snod_id =",f_snod_id
           call oasis_put(f_snod_id,t,f_snod_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           ! Veg low
!                        write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_vegl_id =",f_vegl_id
!                        call oasis_put(f_vegl_id,t,f_vegl_data,ierror)
!                        write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           ! Veg type low
!                        write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_vegltype_id =",f_vegltype_id
!                        call oasis_put(f_vegltype_id,t,f_vegltype_data,ierror)
!                        write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           ! Veg high
!                        write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_vegh_id =",f_vegh_id
!                        call oasis_put(f_vegh_id,t,f_vegh_data,ierror)
!                        write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_put with f_prec_id =",f_prec_id
           call oasis_put(f_prec_id,t,f_prec_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_put returned ierror =",ierror

           ! -----------------------------------------------------------------------------------
           ! *** GET variables
           ! -----------------------------------------------------------------------------------

           write (*,'(A,I6,A,I12)') "*II* lpjg_forcingGET : Time step t = ",t_step," - time t=",t

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_get with f_llai_id =",f_llai_id
           call oasis_get(f_llai_id,t,f_llai_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_get returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_get with f_hlai_id =",f_hlai_id
           call oasis_get(f_hlai_id,t,f_hlai_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_get returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_get with f_typh_id =",f_typh_id
           call oasis_get(f_typh_id,t,f_typh_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_get returned ierror =",ierror

           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_get with f_vggh_id =",f_vggh_id
           call oasis_get(f_vggh_id,t,f_vggh_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_get returned ierror =",ierror

! ecev3
           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_get with f_typl_id =",f_typl_id
           call oasis_get(f_typl_id,t,f_typl_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_get returned ierror =",ierror

! ecev3
           write (*,'(A,I3)') "*II* lpjg_forcing: calling oasis_get with f_vggl_id =",f_vggl_id
           call oasis_get(f_vggl_id,t,f_vggl_data,ierror)
           write (*,'(A,I3)') "*II* lpjg_forcing: oasis_get returned ierror =",ierror

           ! checking whether the incoming LAI data are ok:
           write(*,*)"*II* lpjg_forcing: LLAI obtained from OASIS for test gridcell (g=8000) is ",f_llai_data(8000,1)
           write(*,*)"*II* lpjg_forcing: HLAI obtained from OASIS for test gridcell (g=8000) is ",f_hlai_data(8000,1)
           write(*,*)"*II* lpjg_forcing: TYPH obtained from OASIS for test gridcell (g=8000) is ",f_typh_data(8000,1)
           write(*,*)"*II* lpjg_forcing: VEGH obtained from OASIS for test gridcell (g=8000) is ",f_vggh_data(8000,1)
           write(*,*)"*II* lpjg_forcing: TYPL obtained from OASIS for test gridcell (g=8000) is ",f_typl_data(8000,1)
           write(*,*)"*II* lpjg_forcing: VEGL obtained from OASIS for test gridcell (g=8000) is ",f_vggl_data(8000,1)

           write (*,'(A,I6,A,I12)') "*II* lpjg_forcing DIAGS LAI check: Time step t = ",t_step," - time t=",t

           ! Check for min/max LAI
           minhlai=0
           maxhlai=0
           minllai=0
           maxllai=0
           meanhlai=0
           meanllai=0
           meanvegl=0
           meanvegh=0
           maxcvh=-1
           maxcvl=-1
           maxtvh=-1
           maxtvl=-1
           cover=0.0

           do cell = 1, MAXGRID

              if (f_llai_data(cell,1)>maxllai)then
                 maxllai=f_llai_data(cell,1)
              endif

              meanllai = meanllai + f_llai_data(cell,1)/MAXGRID

              if (f_llai_data(cell,1)<minllai)then
                 minllai=f_llai_data(cell,1)
              endif

              if (f_hlai_data(cell,1)>maxhlai)then
                 maxhlai=f_hlai_data(cell,1)
              endif

              if (f_hlai_data(cell,1)<minhlai)then
                 minhlai=f_hlai_data(cell,1)
              endif

              meanhlai = meanhlai + f_hlai_data(cell,1)/MAXGRID

              if (f_vggl_data(cell,1)>maxcvl)then
                 maxcvl=f_vggl_data(cell,1)
              endif

              if (f_vggh_data(cell,1)>maxcvh)then
                 maxcvh=f_vggh_data(cell,1)
              endif

!              meanvegl = meanvegl + f_vegl_data(cell,1)/MAXGRID
!              meanvegh = meanvegh + f_vegh_data(cell,1)/MAXGRID

              if (f_typl_data(cell,1)>maxtvl)then
                 maxtvl=f_typl_data(cell,1)
              endif

              if (f_typh_data(cell,1)>maxtvh)then
                 maxtvh=f_typh_data(cell,1)
              endif

              cover=f_vggl_data(cell,1)+f_vggh_data(cell,1)
              if (cover>1.001) then
                 write(*,*)"*II* lpjg_forcing: COVER ERROR!",cover
              endif

           enddo

           write(*,*)"*II* lpjg_forcing DIAGS: MAX HLAI :",maxhlai
           write(*,*)"*II* lpjg_forcing: MIN HLAI :",minhlai
           write(*,*)"*II* lpjg_forcing: MAX LLAI :",maxllai
           write(*,*)"*II* lpjg_forcing: MIN LLAI :",minllai
           write(*,*)"*II* lpjg_forcing: MEAN LLAI :",meanllai
           write(*,*)"*II* lpjg_forcing: MEAN HLAI :",meanhlai
!                        write(*,*)"*II* lpjg_forcing: MEAN VEGL :",meanvegl
!                        write(*,*)"*II* lpjg_forcing: MEAN VEGH :",meanvegh
!                        write(*,*)"*II* lpjg_forcing: MAX VEGL :",maxcvl
!                        write(*,*)"*II* lpjg_forcing: MAX VEGH :",maxcvh
           write(*,*)"*II* lpjg_forcing: MAX TYPL :",maxtvl
           write(*,*)"*II* lpjg_forcing: MAX TYPH :",maxtvh

           write (*,*) "*III* lpjg_forcing: Finished dayloop!",t_step, t_max

           ! Increase the full simulation counter
           t_step_full = t_step_full + 1

           ! End of dayloop
        enddo

        write (*,*) "*II* lpjg_forcing: End of year!",loopyr

        ! End of year loop
     enddo

     write (*,'(A)') "*II* lpjg_forcing: Finished time loop!"

     ! End of repetition loop
  enddo

  write (*,'(A)') "*II* lpjg_forcing: oasis_terminate..."
  call oasis_terminate(ierror)
  write (*,'(A,I3)') "*II* lpjg_forcing: oasis_terminate returned ierror =",ierror

end program LPJGFORC


! PM - for pre-Industrial nc files.
! now reads new_ files
! ---------------------------------------------------------------------------------------------------------------
subroutine read_ifstest_pi(NX,NY,t_step,t_max,filename,parname, &
                        data_in,fileid,varid)
! NetCDF
    use netcdf

    implicit none

    integer::NX,NY,t_step,t_max
    character(len=*)::filename
    character(len=*)::parname
    double precision::data_in(NX,NY)
    integer::fileid
    integer::varid
    logical::file_exists

    ! for older dpic_ files:
    !integer::start(2),count(2)
    !start=(/1,t_step+1/)
    !count=(/NX,NY/)

    ! for new_ files:
    integer::start(3),count(3)
    start=(/1,1,t_step+1/)
    count=(/NX,NY,1/)

    !write (*,*) "*II* lpjg_forcing: reading variable ",trim(parname)," from input file ",trim(filename)

    ! check that files exists
    inquire(file=filename, exist=file_exists)
    if ( .not. file_exists) then
      call ERROR("*II* lpjg_forcing: could not read variable "//trim(parname)//" from input file "//trim(filename),1)
    endif

    ! read file properties in first occasion
    if(t_step==0)then
      call check(nf90_open(filename,NF90_NOWRITE,fileid))
      call check(nf90_inq_varid(fileid,parname,varid))
    endif

    ! read field for given timestep
    call check(nf90_get_var(fileid,varid,data_in,start=start,count=count))

    ! close NetCDF file
    if(t_step==t_max-1)then
      call check(nf90_close(fileid))
    endif

end

! PM - for pre-Industrial nc files with DEPTH
! ---------------------------------------------------------------------------------------------------------------
subroutine read_ifstest_pi_depth(NX,NY,t_step,t_max,filename,parname, &
                        data_in,fileid,varid,ignore_file_missing)
! NetCDF
    use netcdf

    implicit none

    integer::NX,NY,t_step,t_max
    character(len=*)::filename
    character(len=*)::parname
    double precision::data_in(NX,NY)
    integer::fileid
    integer::varid
    logical::ignore_file_missing
    logical::file_exists
!    double precision::data_in(NX,NY)
!    integer::fileid
!    integer::varid

    !integer::start(3),count(3)
    !start=(/1,1,t_step+1/)
    !count=(/NX,NY,1/)

    ! new_ files
    integer::start(4),count(4)
    start=(/1,1,1,t_step+1/)
    count=(/NX,NY,1,1/)

    !write (*,*) "*II* lpjg_forcing: reading variable ",trim(parname)," from input file ",trim(filename)

    ! check that files exists
    inquire(file=filename, exist=file_exists)
    if ( .not. file_exists) then
      if (ignore_file_missing) then
        if(t_step==0) write (*,*) "*II* lpjg_forcing: could not read variable ",&
             trim(parname)," from input file ",trim(filename)," continuing since ignore_file_missing=true"
        data_in(:,:) = 0
        return
      else
        call ERROR("*II* lpjg_forcing: could not read variable "//trim(parname)//" from input file "//trim(filename),1)
      endif
    endif

    ! read file properties in first occasion
    if(t_step==0)then
      call check(nf90_open(filename,NF90_NOWRITE,fileid))
      call check(nf90_inq_varid(fileid,parname,varid))
    endif

    ! read field for given timestep
    call check(nf90_get_var(fileid,varid,data_in,start=start,count=count))

    ! close NetCDF file
    if(t_step==t_max-1)then
      call check(nf90_close(fileid))
    endif

end


subroutine ERROR(msg,status)
    character(len=*),intent(in) :: msg
    integer,intent(in) :: status

    write (*,'("*EE*",A)') msg
    call EXIT(status)
end subroutine ERROR


! PM - for fixed nc files.
! ---------------------------------------------------------------------------------------------------------------
subroutine read_fixedfiels(NX,NY,filename,parname,data_in,fileid,varid)

! NetCDF
    use netcdf

    implicit none

    integer::NX,NY
    character(len=*)::filename
    character(len=*)::parname
    double precision::data_in(NX,NY)
    integer::fileid
    integer::varid

    integer::start(2),count(2)

    start=(/1,1/)
    count=(/NX,NY/)

    !write (*,*) "*II* lpjg_forcing: reading variable ",trim(parname)," from input file ",trim(filename)

    ! read file properties in first occasion
    call check(nf90_open(filename,NF90_NOWRITE,fileid))
    call check(nf90_inq_varid(fileid,parname,varid))

    ! read field
    call check(nf90_get_var(fileid,varid,data_in,start=start,count=count))

    ! close NetCDF file
    call check(nf90_close(fileid))

end



! ---------------------------------------------------------------------------------------------------------------
subroutine read_ifstest(NX,NY,t_step,t_max,filename,parname, &
                        data_in,fileid,varid)
! NetCDF
    use netcdf

    implicit none

    integer::NX,NY,t_step,t_max
    character(len=*)::filename
    character(len=*)::parname
    double precision::data_in(NX,NY)
    integer::fileid
    integer::varid

    integer::start(3),count(3)

    start=(/1,1,t_step+1/)
    count=(/NX,NY,1/)

    ! read file properties in first occasion
    if(t_step==0)then
      call check(nf90_open(filename,NF90_NOWRITE,fileid))	
      call check(nf90_inq_varid(fileid,parname,varid))
    endif

    ! read field for given timestep
    call check(nf90_get_var(fileid,varid,data_in,start=start,count=count))

    ! close NetCDF file
    if(t_step==t_max-1)then
      call check(nf90_close(fileid))
    endif

end

! ---------------------------------------------------------------------------------------------------------------
subroutine check(status)

    use netcdf

    integer, intent ( in) :: status

    if(status /= nf90_noerr) then
      print *, trim(nf90_strerror(status))
!      stop 1
      call abort
    end if
  end subroutine check
! ---------------------------------------------------------------------------------------------------------------