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)maxhlai)then maxhlai=f_hlai_data(cell,1) endif if (f_hlai_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 ! ---------------------------------------------------------------------------------------------------------------