ecearth3/sources/lpjg/offline/sparring_old.F90

program SPARRING

! 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)

! Last update: 10 Dec 2015 - add 2 new fields for low veg. - see "ecev3"
! Now reads AMIP output.

    use mod_oasis
    use mpi
 
    implicit none

!    integer,parameter :: MAXGRID = 35718
    integer,parameter :: MAXGRID = 88838
   
    integer,parameter :: num_years = 10
    integer,parameter :: start_year = 1850

    integer,parameter :: num_loops = 10
	
    integer,parameter :: t_max_year = 365 
    integer,parameter :: t_max_leap = 366
	
    character(len=6),parameter :: model_name    = "SPARRI"

	
    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_snowfall_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_radi_name = "SSRVeg"
    integer                    :: f_radi_id   = -77
    double precision           :: f_radi_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
!    character(len=*),parameter::nc_temp_file="/nobackup/rossby15/sm_paumi/preind/dpic_var167+1849.nc"
!    character(len=*),parameter::nc_prec_file="/nobackup/rossby15/sm_paumi/preind/dpic_var142+1849.nc"
!    character(len=*),parameter::nc_radi_file="/nobackup/rossby15/sm_paumi/preind/dpic_var169+1849.nc"

! VARIABLE FIELDS
	character(len=100) :: nc_temp_file
	character(len=100) :: nc_precs_file
        character(len=100) :: nc_precc_file
        character(len=100) :: nc_snowfall_file
	character(len=100) :: nc_st1l_file
        character(len=100) :: nc_st2l_file
        character(len=100) :: nc_st3l_file
        character(len=100) :: nc_st4l_file
        character(len=100) :: nc_sm1l_file
        character(len=100) :: nc_sm2l_file
        character(len=100) :: nc_sm3l_file
        character(len=100) :: nc_sm4l_file
	character(len=100) :: nc_radi_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"

! STEM FOR VARIABLE FIELDS
!	character(len=*),parameter::nc_file_path="/nobackup/rossby15/sm_paumi/ecev3/processed/"

! NEW, AMIP OUTPUT:
        character(len=*),parameter::nc_file_path="/nobackup/rossby18/sm_wyser/ecearth3-amip/processed/daily/"
        character(len=*),parameter::nc_file_tail="_dayavg.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_snowfall_name="var144"

    ! We do not have the SNOC and SNOD data, so we send fixed values to guess below.
 
    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.
    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_radi_name="var176"	
    character(len=*),parameter::nc_vegl_name="tvl"	
    character(len=*),parameter::nc_vegh_name="tvh"	

	
! FIELD IDS
    integer::nc_temp_fileid,nc_precs_fileid,nc_radi_fileid,nc_vegl_fileid,nc_vegh_fileid,nc_precc_fileid,nc_snowfall_fileid;
    integer::nc_temp_varid,nc_precs_varid,nc_radi_varid,nc_vegl_varid,nc_vegh_varid,nc_precc_varid,nc_snowfall_varid;   
    integer::nc_vegltype_varid; ! We don't have a T255 file for vegltype yet 
 
    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* sparring: Hello" 

    write (*,'(A)') "*II* sparring: Now initialising Sparring using oasis_..."

    call oasis_init_comp(comp_id,model_name,ierror)
    write (*,'(A,I3)') "*II* sparring: oasis_init_comp returned  ierror=",ierror
   
!    call oasis_get_localcomm(localcomm, ierror)
!    write (*,'(A,I3)') "*II* sparring: get_localcomm  ierror=",ierror
!    write (*,'(A,I12)') "*II* sparring: 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* sparring: oasis_create_couplcomm  ierror=",ierror
    write (*,'(A,I12)') "*II* sparring: oasis_create_couplcomm returned cplcomm =",cplcomm

    call oasis_def_partition(part_id,(/ 0,0,MAXGRID*1 /),ierror)
    write (*,'(A,I3)') "*II* sparring: oasis_def_partition returned part_id =",part_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_temp_id =",f_temp_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_prec_id =",f_prec_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_vegl_id =",f_vegl_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_vegh_id =",f_vegh_id
    write (*,'(A,I3)') "*II* sparring: oasis_def_var 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* sparring: oasis_def_var returned f_snoc_id =",f_snoc_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_snod_id =",f_snod_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_st1l_id =",f_st1l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_st2l_id =",f_st2l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_st3l_id =",f_st3l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_st4l_id =",f_st4l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_sm1l_id =",f_sm1l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_sm2l_id =",f_sm2l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_sm3l_id =",f_sm3l_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_sm4l_id =",f_sm4l_id
    write (*,'(A,I3)') "*II* sparring: oasis_def_var returned    ierror =",ierror



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


! VEGLTYPE

 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* sparring: f_vegltype_id =",f_vegltype_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_llai_id =",f_llai_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_hlai_id =",f_hlai_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_typh_id =",f_typh_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_vggh_id =",f_vggh_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_typl_id =",f_typl_id
    write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_def_var returned f_vggl_id =",f_vggl_id
    write (*,'(A,I3)') "*II* sparring: oasis_def_var returned    ierror =",ierror

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

! FIXED FIELDS

    write (*,'(A)') "*II* sparring: 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* sparring: Beginning time loop"


	t_step_full = 0

! Repetitions of forcing	
    do lyr = 1, num_loops
 
! Year loop	
    do loopyr = start_year, start_year+num_years-1

    yr = loopyr
    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

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

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

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

        nc_snowfall_file=nc_file_path//nc_snowfall_name//"_"//yearstr//nc_file_tail
        write (*,'(A)') nc_snowfall_file

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

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

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

        !nc_st4l_file = nc_file_path // nc_st4l_name // "_" // yearstr // ".nc"
        !write (*,'(A)') nc_st4l_file

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

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

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

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

	nc_radi_file = nc_file_path // nc_radi_name//"_"//yearstr//nc_file_tail
	write (*,'(A)') nc_radi_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* sparring: t_max =",t_max
			   
		do t_step=0,t_max-1

			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* sparring 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)
			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)

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


                        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)

                        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)

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


                        !call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_st3l_file,nc_st3l_name,&
                        !                                          f_st3l_data,nc_st3l_fileid,nc_st3l_varid)

                        !call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_st4l_file,nc_st4l_name,&
                        !                                          f_st4l_data,nc_st4l_fileid,nc_st4l_varid)

                        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)

                        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)

                        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)

                        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)

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



                        call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_snowfall_file,nc_snowfall_name,&
                                                                  f_snowfall_data,nc_snowfall_fileid,nc_snowfall_varid)


                        call read_ifstest_pi(MAXGRID,1,t_step,t_max,nc_radi_file,nc_radi_name,&
                                                                  f_radi_data,nc_radi_fileid,nc_radi_varid)


			! Sum precipitation and snow components
			! Replace ST3L and ST4L with values from ST2L
			do cell = 1, MAXGRID
				! Total precip includes snow here
				f_prec_data(cell,1) = 4 * 1000.0 *(f_precs_data(cell,1) +  f_precc_data(cell,1) +  f_snowfall_data(cell,1))
                                ! 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)=f_prec_data(cell,1)/86400.0

				f_st3l_data(cell,1) = f_st2l_data(cell,1)
				f_st4l_data(cell,1) = f_st2l_data(cell,1)
				! Snow
				f_snoc_data(cell,1) = 0.0
				f_snod_data(cell,1) = 330.0
				! PM_Apr2012 - change J m-2 to W m-2 to be consistent with ECE. We'll remove the / in LPJG.
				f_radi_data(cell,1) = 4* f_radi_data(cell,1) / 86400.0 
			enddo

                        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* sparring after ifstest: Time step t = ",t_step," - time t=",t

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

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

			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror
	  
                        write (*,'(A,I3)') "*II* sparring: calling oasis_put with f_radi_id =",f_radi_id
                        call oasis_put(f_radi_id,t,f_radi_data,ierror)
                        write (*,'(A,I3)') "*II* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

			! Snow mass/unit surface (kg/m2)
                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        ! Snow density (kg/m3)
                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        ! Veg low
                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror
                        
                        ! Veg type low
                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror



                        ! Veg high
                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

                        write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_put returned ierror =",ierror

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

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

			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_get returned ierror =",ierror

			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_get returned ierror =",ierror
	  
			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_get returned ierror =",ierror
			
			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_get returned ierror =",ierror
			
! ecev3
			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_get returned ierror =",ierror
	
! ecev3		
			write (*,'(A,I3)') "*II* sparring: 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* sparring: oasis_get returned ierror =",ierror


			! checking whether the incoming LAI data are ok:
			write(*,*)"*II* sparring: LLAI obtained from OASIS for test gridcell (g=8000) is ",f_llai_data(8000,1)
			write(*,*)"*II* sparring: HLAI obtained from OASIS for test gridcell (g=8000) is ",f_hlai_data(8000,1)
			write(*,*)"*II* sparring: TYPH obtained from OASIS for test gridcell (g=8000) is ",f_typh_data(8000,1)
 			write(*,*)"*II* sparring: VEGH obtained from OASIS for test gridcell (g=8000) is ",f_vggh_data(8000,1)
			write(*,*)"*II* sparring: TYPL obtained from OASIS for test gridcell (g=8000) is ",f_typl_data(8000,1)
 			write(*,*)"*II* sparring: VEGL obtained from OASIS for test gridcell (g=8000) is ",f_vggl_data(8000,1)
  
           		write (*,'(A,I6,A,I12)') "*II* sparring 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* sparring: COVER ERROR!",cover
                                endif

			enddo

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


    write (*,*) "*II* sparring: Finished dayloop!",t_step, t_max
                        ! Increase the full simulation counter
                        t_step_full = t_step_full + 1
		
		! End of dayloop
		enddo
    write (*,*) "*II* sparring: End of year!",loopyr

! End of year loop	
    enddo

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

! End of repetition loop
enddo 

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

end program SPARRING


! 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

    ! 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/)

    ! 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)
! 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
!    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/)

    ! 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 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/)

    ! 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 "Stopped"
    end if
  end subroutine check  
! ---------------------------------------------------------------------------------------------------------------