diff --git a/src/global.f90 b/src/global.f90 index 3e1938b1..0a8d6c27 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -5415,8 +5415,10 @@ real(dp) function SeasonalSumOfKcPot(TheDaysToCCini, TheGDDaysToCCini, L0, L12, ! 3. Calculate Sum do Dayi = 1, L1234 ! 3.1 calculate growing degrees for the day + !WRITE(*,*) Dayi,L1234,L123,GDDL1234 if (GetTemperatureFile() /= '(None)') then read(fhandle, *) Tndayi, Txdayi + !WRITE(*,*) Tndayi, Txdayi GDDi = DegreesDay(Tbase, Tupper, Tndayi, Txdayi, & GetSimulParam_GDDMethod()) else diff --git a/src/project_input.f90 b/src/project_input.f90 index 7adf8eac..c19b77e2 100644 --- a/src/project_input.f90 +++ b/src/project_input.f90 @@ -159,6 +159,8 @@ subroutine initialize_project_input(filename, NrRuns) integer :: i, NrRuns_local +! write(*,*)'IPI',NrRuns + if (present(NrRuns)) then NrRuns_local = NrRuns else @@ -169,6 +171,7 @@ subroutine initialize_project_input(filename, NrRuns) do i = 1, NrRuns_local call ProjectInput(i)%read_project_file(filename, i) end do +! write(*,*)'IPI_local',NrRuns_local end subroutine initialize_project_input @@ -237,16 +240,26 @@ subroutine read_project_file(self, filename, NrRun) iostat=rc) read(fhandle, '(a)', iostat=rc) buffer self%Description = trim(buffer) +! write(*,*)'RPF ',self%Description read(fhandle, *, iostat=rc) self%VersionNr ! AquaCrop version Nr +! write(*,*)'RPF ',self%VersionNr +! write(*,*)'RPF ',NrRun if (NrRun > 1) then ! Skip sections belonging to previous runs do Runi = 1, (NrRun - 1) do i = 1, 47 - read(fhandle, *, iostat=rc) ! 5 + 42 lines with files + read(fhandle, *, iostat=rc) buffer ! 5 + 42 lines with files !buffer added +! if (i.eq.1) then +! write(*,*)'RPFSkip ',trim(buffer) +! end if end do end do end if + if (NrRun > 127) then !PCM + read(fhandle, *, iostat=rc) buffer ! 5 + 42 lines with files !buffer added +! write(*,*)'RPFSkip ',trim(buffer) + end if ! 0. Year of cultivation and Simulation and Cropping period read(fhandle, *, iostat=rc) self%Simulation_YearSeason @@ -258,8 +271,10 @@ subroutine read_project_file(self, filename, NrRun) ! 1. Climate read(fhandle, '(a)', iostat=rc) buffer self%Climate_Info = trim(buffer) +! write(*,*)'RPF ',self%Climate_Info read(fhandle, *, iostat=rc) buffer self%Climate_Filename = trim(buffer) +! write(*,*)'RPF ',self%Climate_Filename read(fhandle, *, iostat=rc) buffer self%Climate_Directory = trim(buffer) @@ -268,6 +283,7 @@ subroutine read_project_file(self, filename, NrRun) self%Temperature_Info = trim(buffer) read(fhandle, *, iostat=rc) buffer self%Temperature_Filename = trim(buffer) +! write(*,*)'RPF ',self%Temperature_Filename read(fhandle, *, iostat=rc) buffer self%Temperature_Directory = trim(buffer) diff --git a/src/run.f90 b/src/run.f90 index 4feeef1b..b419c69d 100644 --- a/src/run.f90 +++ b/src/run.f90 @@ -566,6 +566,7 @@ module ac_run integer(int8) :: StageCode integer(int32) :: PreviousDayNr logical :: NoYear +logical :: debug = .FALSE. character(len=:), allocatable :: fEval_filename @@ -3490,7 +3491,7 @@ subroutine CheckForPrint(TheProjectFile) SaltIn = GetSumWaBal_SaltIn() - GetPreviousSum_SaltIn() SaltOut = GetSumWaBal_SaltOut() - GetPreviousSum_SaltOut() CRsalt = GetSumWaBal_CRsalt() - GetPreviousSum_CRsalt() - call WriteTheResults(int(undef_int,kind=int8), DayN, MonthN, YearN, DayN, MonthN, & + call WriteTheResults(int(undef_int,kind=int32), DayN, MonthN, YearN, DayN, MonthN, & YearN, GetRain(), GetETo(), GetGDDayi(), GetIrrigation(), & GetInfiltrated(), GetRunoff(), GetDrain(), & GetCRwater(), GetEact(), GetEpot(), GetTact(), & @@ -3964,6 +3965,9 @@ subroutine RelationshipsForFertilityAndSaltStress() integer(int8) :: BioTop, BioLow real(dp) :: StrTop, StrLow + if(debug)then + write(*,*)'BB0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 1. Soil fertility call SetFracBiomassPotSF(1._dp) @@ -4028,6 +4032,9 @@ subroutine RelationshipsForFertilityAndSaltStress() call SetFracBiomassPotSF(GetFracBiomassPotSF()/100._dp) end if + if(debug)then + write(*,*)'BB GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 2. soil salinity (Coeffb0Salt,Coeffb1Salt,Coeffb2Salt : CCx/KsSto - Salt stress) if (GetSimulation_SalinityConsidered() .eqv. .true.) then call CCxSaltStressRelationship(GetCrop_DaysToCCini(), & @@ -4115,7 +4122,7 @@ end subroutine DetermineGrowthStage subroutine WriteTitleDailyResults(TheProjectType, TheNrRun) integer(intenum), intent(in) :: TheProjectType - integer(int8), intent(in) :: TheNrRun + integer(int32), intent(in) :: TheNrRun character(len=1025) :: Str1, Str2, tempstring real(dp) :: NodeD, Zprof @@ -4125,7 +4132,7 @@ subroutine WriteTitleDailyResults(TheProjectType, TheNrRun) call fDaily_write('') if (TheProjectType == typeproject_TypePRM) then write(Str1, '(i4)') TheNrRun - call fDaily_write(' Run:'// trim(Str1)) + call fDaily_write('# Run:'// trim(Str1)) end if ! B. thickness of soil profile and root zone @@ -4213,10 +4220,10 @@ subroutine WriteTitleDailyResults(TheProjectType, TheNrRun) if (GetOut5CompWC()) then call fDaily_write(trim(' WC01'), .false.) do Compi = 2, (GetNrCompartments()-1) - write(Str1, '(i2)') Compi + write(Str1, '(i0.2)') Compi call fDaily_write(' WC'// trim(Str1), .false.) end do - write(Str1,'(i2)') GetNrCompartments() + write(Str1,'(i0.2)') GetNrCompartments() if ((GetOut6CompEC()) .or. (GetOut7Clim())) then call fDaily_write(' WC'// trim(Str1), .false.) else @@ -4227,10 +4234,10 @@ subroutine WriteTitleDailyResults(TheProjectType, TheNrRun) if (GetOut6CompEC()) then call fDaily_write(trim(' ECe01'), .false.) do Compi = 2, (GetNrCompartments()-1) - write(Str1, '(i2)') Compi + write(Str1, '(i0.2)') Compi call fDaily_write(' ECe'// trim(Str1), .false.) end do - write(Str1, '(i2)') GetNrCompartments() + write(Str1, '(i0.2)') GetNrCompartments() if (GetOut7Clim()) then call fDaily_write(' ECe'// trim(Str1), .false.) else @@ -4247,13 +4254,13 @@ subroutine WriteTitleDailyResults(TheProjectType, TheNrRun) if (GetOut1Wabal()) then if ((GetOut2Crop()) .or. (GetOut3Prof()) .or. (GetOut4Salt()) & .or. (GetOut5CompWC()) .or. (GetOut6CompEC()) .or. (GetOut7Clim())) then - write(tempstring, '(2a)') ' mm mm mm mm' // & + write(tempstring, '(2a)') '# mm mm mm mm' // & ' mm mm mm mm m ', & ' mm mm % mm mm %' // & ' mm mm %' call fDaily_write(trim(tempstring), .false.) else - write(tempstring, '(2a)') ' mm mm mm mm' // & + write(tempstring, '(2a)') '# mm mm mm mm' // & ' mm mm mm mm m ', & ' mm mm % mm mm %' // & ' mm mm %' @@ -4353,7 +4360,7 @@ end subroutine WriteTitleDailyResults subroutine FinalizeRun2(NrRun, TheProjectType) - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun integer(intenum), intent(in) :: TheProjectType call CloseClimateFiles() @@ -4369,7 +4376,7 @@ subroutine FinalizeRun2(NrRun, TheProjectType) subroutine CloseEvalDataPerformEvaluation(NrRun) - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun character(len=:), allocatable :: totalnameEvalStat character(len=1024) :: StrNr @@ -4544,38 +4551,43 @@ end subroutine OpenIrrigationFile subroutine WriteTitlePart1MultResults(TheProjectType, TheNrRun) integer(intEnum), intent(in) :: TheProjectType - integer(int8), intent(in) :: TheNrRun + integer(int32), intent(in) :: TheNrRun integer(int32) :: Dayi, Monthi, Yeari real(dp) :: Nr character(len=1024) :: tempstring ! A. Run number - call fHarvest_write('') - if (TheProjectType == typeproject_TypePRM) then + if(TheNrRun.EQ.1) then + call fHarvest_write('') + if (TheProjectType == typeproject_TypePRM) then write(tempstring, '(i4)') TheNrRun - call fHarvest_write(' Run:' // trim(tempstring)) - end if + call fHarvest_write('# Run:' // trim(tempstring)) + end if ! B. Title - call fHarvest_write(' Nr Day Month Year DAP Interval Biomass ' // & - 'Sum(B) Dry-Yield Sum(Y) Fresh-Yield Sum(Y)') - call fHarvest_write(' days ton/ha ' // & - 'ton/ha ton/ha ton/ha ton/ha ton/ha') + !call fHarvest_write(' Nr Day Month Year DAP Interval Biomass ' // & + ! 'Sum(B) Dry-Yield Sum(Y) Fresh-Yield Sum(Y)') + call fHarvest_write(' Nr Day Month Year ' // & + ' Biomass-Sum(B) Dry-Yield-Sum(Y) Fresh-Yield-Sum(Y)') + call fHarvest_write('# ' // & + ' ton/ha ton/ha ton/ha') ! C. start crop cycle - call DetermineDate(GetCrop_Day1(), Dayi, Monthi, Yeari) - call SetNoYear(Yeari == 1901) - if (GetNoYear()) then + call DetermineDate(GetCrop_Day1(), Dayi, Monthi, Yeari) + call SetNoYear(Yeari == 1901) + if (GetNoYear()) then if (Dayi == 0) then Dayi = 1 end if Yeari = 9999 + end if + Nr = 0._dp + !write(tempstring, '(i6, 3i6, f34.3, 2f20.3)') 0, Dayi, Monthi, Yeari, & + ! Nr, Nr, Nr + !write(tempstring, '(i6, 3i6, f34.3)') 0, Dayi, Monthi, Yeari + !call fHarvest_write(trim(tempstring)) end if - Nr = 0._dp - write(tempstring, '(i6, 3i6, f34.3, 2f20.3)') 0, Dayi, Monthi, Yeari, & - Nr, Nr, Nr - call fHarvest_write(trim(tempstring)) end subroutine WriteTitlePart1MultResults @@ -4585,7 +4597,7 @@ subroutine WriteTheResults(ANumber, Day1, Month1, Year1, DayN, MonthN, & TrxPer, SalInPer, SalOutPer, & SalCRPer, BiomassPer, BUnlimPer, BmobPer, BstoPer, & TheProjectFile) - integer(int8), intent(in) :: ANumber + integer(int32), intent(in) :: ANumber integer(int32), intent(in) :: Day1 integer(int32), intent(in) :: Month1 integer(int32), intent(in) :: Year1 @@ -4767,6 +4779,15 @@ subroutine InitializeSimulationRunPart1() integer(int32) :: Crop_DaysToFullCanopySF_temp logical :: WaterTableInProfile_temp + if(debug)then + write(*,*)'AA0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if + if(GetCrop_DaysToHarvest().gt.364)THEN + write(*,*)'AA1 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + write(*,*)'growing season too short for GDD: returning' + return + end if + ! 1. Adjustments at start ! 1.1 Adjust soil water and salt content if water table IN soil profile WaterTableInProfile_temp = GetWaterTableInProfile() @@ -4855,6 +4876,9 @@ subroutine InitializeSimulationRunPart1() call SetStressTot_Sto(0._dp) call SetStressTot_Weed(0._dp) + if(debug)then + write(*,*)'AA1 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 6. Soil fertility stress ! Coefficients for soil fertility - biomass relationship ! AND for Soil salinity - CCx/KsSto relationship @@ -4897,6 +4921,10 @@ subroutine InitializeSimulationRunPart1() end if end if + if(debug)then + write(*,*)'GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if + ! Maximum sum Kc (for reduction WP in season if soil fertility stress) call SetSumKcTop(SeasonalSumOfKcPot(GetCrop_DaysToCCini(), & GetCrop_GDDaysToCCini(), GetCrop_DaysToGermination(), & @@ -5449,7 +5477,7 @@ end subroutine InitializeSimulationRunPart2 subroutine CreateEvalData(NrRun) - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun integer(int32) :: dayi, monthi, yeari, integer_temp character(len=:), allocatable :: TempString @@ -5508,13 +5536,13 @@ subroutine CreateEvalData(NrRun) call fEval_open(GetfEval_filename(), 'w') write(tempstring2, '(a)') GetAquaCropDescriptionWithTimeStamp() call fEval_write(trim(tempstring2)) - call fEval_write('Evaluation of simulation results - Data') + call fEval_write('# Evaluation of simulation results - Data') write(TempString, '(f5.2)') GetZeval() - call fEval_write(' ' // & + call fEval_write('# ' // & ' for soil depth: ' // trim(TempString) // ' m') call fEval_write(' Day Month Year DAP Stage CCsim CCobs CCstd Bsim ' // & 'Bobs Bstd SWCsim SWCobs SWstd') - call fEval_write(' % % % ton/ha ' // & + call fEval_write('# % % % ton/ha ' // & 'ton/ha ton/ha mm mm mm') end subroutine CreateEvalData @@ -5542,7 +5570,7 @@ subroutine OpenOutputRun(TheProjectType) ' SaltIn SaltOut SaltUp SaltProf' // & ' Cycle SaltStr FertStr WeedStr TempStr ExpStr StoStr' // & ' BioMass Brelative HI Y(dry) Y(fresh) WPet Bin Bout DayN MonthN YearN') - call fRun_write(' mm mm degC.day ppm' // & + call fRun_write('# mm mm degC.day ppm' // & ' mm mm mm mm mm mm % mm mm %' // & ' ton/ha ton/ha ton/ha ton/ha' // & ' days % % % % % % ' // & @@ -5585,7 +5613,7 @@ subroutine OpenPart1MultResults(TheProjectType) call fHarvest_open(GetfHarvest_filename(), 'w') write(tempstring, '(a)') GetAquaCropDescriptionWithTimeStamp() call fHarvest_write(trim(tempstring)) - call fHarvest_write('Biomass and Yield at Multiple cuttings') + call fHarvest_write('# Biomass and Yield at Multiple cuttings') end subroutine OpenPart1MultResults @@ -5603,6 +5631,9 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) type(rep_DayEventDbl), dimension(31) :: TminDataSet_temp, TmaxDataSet_temp real(dp) :: Tmin_temp, Tmax_temp type(rep_DayEventDbl), dimension(31) :: EToDataSet_temp, RainDataSet_temp + if(debug .eqv. .TRUE.)then + write(*,*)'RR0 FromSimDay, ToSimDay=',FromSimDay, ToSimDay + end if ! 1. ETo file if (GetEToFile() /= '(None)') then @@ -5707,6 +5738,7 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) ! 2. Rain File if (GetRainFile() /= '(None)') then totalname = GetRainFilefull() + RunningDay = FromSimDay if (FileExists(totalname)) then ! open file and find first day of simulation period select case (GetRainRecord_DataType()) @@ -5736,8 +5768,14 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) end do call SetRain(GetRainDataSet_Param(i)) case(datatype_Monthly) + if(debug .eqv. .TRUE.)then + write(*,*)'RR1 FromSimDay, RunningDay=',FromSimDay, RunningDay + end if RainDataSet_temp = GetRainDataSet() call GetMonthlyRainDataSet(FromSimDay, RainDataSet_temp) + if(debug .eqv. .TRUE.)then + write(*,*)'RR1b FromSimDay, RunningDay=',FromSimDay, RunningDay + end if call SetRainDataSet(RainDataSet_temp) i = 1 do while (GetRainDataSet_DayNr(i) /= FromSimDay) @@ -5752,7 +5790,13 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) action='write') write(fRainS, '(f10.4)') GetRain() ! next days of simulation period + if(debug .eqv. .TRUE.)then + write(*,*)'RR2 FromSimDay, ToSimDay=',FromSimDay, ToSimDay + end if do RunningDay = (FromSimDay + 1), ToSimDay + if(debug)then + write(*,*)'RR2 RunningDay=',RunningDay + end if select case (GetRainRecord_DataType()) case(datatype_daily) if (rc == iostat_end) then @@ -5791,8 +5835,18 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) i = 1 do while (GetRainDataSet_DayNr(i) /= RunningDay) i = i+1 + if(debug)then + write(*,*)'RR2a RunningDay,i=',RunningDay,i + write(*,*)'RR2a GRDS=',GetRainDataSet_DayNr(i) + end if end do + if(debug)then + write(*,*)'RR2b RunningDay,i=',RunningDay,i + end if call SetRain(GetRainDataSet_Param(i)) + if(debug)then + write(*,*)'RR2c RunningDay,i=',RunningDay,i + end if end select write(fRainS, '(f10.4)') GetRain() end do @@ -5803,6 +5857,9 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) close(fRainS) end if end if + if(debug)then + write(*,*)'RR2d FromSimDay, ToSimDay=',FromSimDay, ToSimDay + end if ! 3. Temperature file if (GetTemperatureFile() /= '(None)') then @@ -5864,7 +5921,13 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) action='write') write(fTempS, '(2f10.4)') GetTmin(), GetTmax() ! next days of simulation period + if(debug)then + write(*,*)'RR3 FromSimDay, ToSimDay=',FromSimDay, ToSimDay + end if do RunningDay = (FromSimDay + 1), ToSimDay + if(debug)then + write(*,*)'RR3 RunningDay=',RunningDay + end if select case (GetTemperatureRecord_Datatype()) case(datatype_Daily) if (rc == iostat_end) then @@ -5919,8 +5982,17 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) do while (GetTminDataSet_DayNr(i) /= RunningDay) i = i+1 end do + if(debug)then + write(*,*)'RR3b RunningDay,i=',RunningDay,i + write(*,*)'RR3a GRDS=',GetTminDataSet_DayNr(i) + write(*,*)'RR3c Tmin=',GetTminDataSet_Param(i) + write(*,*)'RR3c Tmax=',GetTmaxDataSet_Param(i) + end if call SetTmin(GetTminDataSet_Param(i)) call SetTmax(GetTmaxDataSet_Param(i)) + if(debug)then + write(*,*)'RR3d' + end if end select write(fTempS, '(2f10.4)') GetTmin(), GetTmax() end do @@ -5931,6 +6003,10 @@ subroutine CreateDailyClimFiles(FromSimDay, ToSimDay) close(fTempS) end if end if + if(debug)then + write(*,*)'RR4 FromSimDay, ToSimDay=',FromSimDay, ToSimDay + write(*,*)'RR4 RunningDay=',RunningDay + endif end subroutine CreateDailyClimFiles @@ -6062,7 +6138,7 @@ end subroutine FinalizeSimulation subroutine WriteSimPeriod(NrRun, TheProjectFile) - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun character(len=*), intent(in) :: TheProjectFile integer(int32) :: Day1, Month1, Year1, DayN, MonthN, YearN @@ -6120,7 +6196,7 @@ subroutine WriteIntermediatePeriod(TheProjectFile) BstoPer = GetSimulation_Storage_Btotal() - GetPreviousBsto() ! write - call WriteTheResults(int(undef_int, kind=int8), Day1, Month1, Year1, DayN, & + call WriteTheResults(int(undef_int, kind=int32), Day1, Month1, Year1, DayN, & MonthN, YearN, RPer, EToPer, GDDPer, IrriPer, InfiltPer, & ROPer, DrainPer, CRwPer, EPer, ExPer, TrPer, TrWPer, & TrxPer, SalInPer, SalOutPer, SalCRPer, BiomassPer, & @@ -6592,7 +6668,7 @@ subroutine InitializeRunPart1(NrRun, TheProjectType) !! Loads the run input from the project file !! Initializes parameters and states !! Calls InitializeSimulationRunPart1 - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun integer(intEnum), intent(in) :: TheProjectType type(rep_sum) :: SumWaBal_temp, PreviousSum_temp @@ -6601,8 +6677,14 @@ subroutine InitializeRunPart1(NrRun, TheProjectType) ! Do nothing return end if + if(debug)then + write(*,*)'CC0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call LoadSimulationRunProject(int(NrRun, kind=int32)) + if(debug)then + write(*,*)'CC1 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call AdjustCompartments() SumWaBal_temp = GetSumWaBal() @@ -6611,6 +6693,9 @@ subroutine InitializeRunPart1(NrRun, TheProjectType) PreviousSum_temp = GetPreviousSum() call ResetPreviousSum(PreviousSum_temp) call SetPreviousSum(PreviousSum_temp) + if(debug)then + write(*,*)'CC2 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call InitializeSimulationRunPart1() contains @@ -6673,12 +6758,12 @@ subroutine InitializeRunPart2(NrRun, TheProjectType) !! Part2 (after reading the climate) of the run initialization !! Calls InitializeSimulationRunPart2 !! Initializes write out for the run - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun integer(intEnum), intent(in) :: TheProjectType call InitializeSimulationRunPart2() - if (GetOutDaily()) then + if (GetOutDaily().AND.(NrRun.EQ.1)) then call WriteTitleDailyResults(TheProjectType, NrRun) end if @@ -6712,7 +6797,7 @@ subroutine RecordHarvest(NrCut, DayInSeason) if (NrCut == 9999) then ! last line at end of season - write(tempstring, '(4i6, f34.3)') NrCut, Dayi, Monthi, Yeari, & + write(tempstring, '(4i6, f20.3)') NrCut, Dayi, Monthi, Yeari, & GetSumWaBal_Biomass() call fHarvest_write(trim(tempstring), .false.) if (GetCrop_DryMatter() == undef_int) then @@ -7388,7 +7473,7 @@ end subroutine SetGDDVariablesNextDay subroutine FinalizeRun1(NrRun, TheProjectFile, TheProjectType) - integer(int8), intent(in) :: NrRun + integer(int32), intent(in) :: NrRun character(len=*), intent(in) :: TheProjectFile integer(intEnum), intent(in) :: TheProjectType @@ -7766,10 +7851,13 @@ subroutine FileManagement() RepeatToDay = GetSimulation_ToDayNr() loop: do + !write(*,*)'FM1 ',GetDayNri()-1,RepeatToDay call AdvanceOneTimeStep(WPi) call ReadClimateNextDay() call SetGDDVariablesNextDay() + !write(*,*)'FM2 ',GetDayNri()-1,RepeatToDay if ((GetDayNri()-1) == RepeatToDay) exit loop + if (GetDayNri() == 364) exit loop end do loop end subroutine FileManagement @@ -7778,9 +7866,12 @@ subroutine RunSimulation(TheProjectFile_, TheProjectType) character(len=*), intent(in) :: TheProjectFile_ integer(intEnum), intent(in) :: TheProjectType - integer(int8) :: NrRun + integer(int32) :: NrRun integer(int32) :: NrRuns + if(debug)then + write(*,*)'DD0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call InitializeSimulation(TheProjectFile_, TheProjectType) select case (TheProjectType) @@ -7794,6 +7885,20 @@ subroutine RunSimulation(TheProjectFile_, TheProjectType) call InitializeRunPart1(NrRun, TheProjectType) call InitializeClimate() call InitializeRunPart2(NrRun, TheProjectType) + + if(debug)then + write(*,*)'DD1 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if + if(GetCrop_DaysToHarvest().gt.364)THEN + write(*,*)'DD1b GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + write(*,*)'growing season too short for GDD: skipping year' + !call SetDayNri(GetSimulation_ToDayNr()) + !call RecordHarvest((9999), & + ! (GetDayNri() - GetCrop_Day1()+1)) ! last line at end of season + !call WriteSimPeriod(NrRun, TheProjectFile) + !cycle + call SetDayNri(GetSimulation_FromDayNr()) + end if call FileManagement() call FinalizeRun1(NrRun, GetTheProjectFile(), TheProjectType) call FinalizeRun2(NrRun, TheProjectType) diff --git a/src/tempprocessing.f90 b/src/tempprocessing.f90 index 409134cc..ae0ed20f 100644 --- a/src/tempprocessing.f90 +++ b/src/tempprocessing.f90 @@ -607,6 +607,7 @@ subroutine GetMonthlyTemperatureDataSet(DayNri, TminDataSet, TmaxDataSet) real(dp) :: aOver3Max, bOver2Max, cMax call DetermineDate(DayNri, Dayi, Monthi, Yeari) +! write(*,*)'TTT0,DayN,Dayi,Monthi,Yeari',DayNri,Dayi,Monthi,Yeari call GetSetofThreeMonths(Monthi, Yeari, & C1Min, C2Min, C3Min, C1Max, C2Max, C3Max, X1, X2, X3, t1) @@ -683,9 +684,11 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & Yfile = GetTemperatureRecord_FromY() end if OK3 = .false. +! write(*,*)'TTT1,Monthi,Yeari',Monthi,Yeari ! 2. IF 3 or less records if (GetTemperatureRecord_NrObs() <= 3) then +! write(*,*)'TTT2' call ReadMonth(C1Min, C1Max, fhandle, rc) X1 = n1 select case (GetTemperatureRecord_NrObs()) @@ -743,20 +746,24 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & ! 3. If first observation if ((.not. OK3) .and. ((Monthi == Mfile) .and. (Yeari == Yfile))) then +! write(*,*)'TTT3' t1 = 0 call ReadMonth(C1Min, C1Max, fhandle, rc) +! write(*,*)'TTT3a,Mfile,Y,Obsi,TMx=',Mfile,Yfile,Obsi,C1Max/30.0 X1 = n1 Mfile = Mfile + 1 if (Mfile > 12) then call AdjustMONTHandYEAR(Mfile, Yfile) end if call ReadMonth(C2Min, C2Max, fhandle, rc) +! write(*,*)'TTT3a,Mfile,Y,Obsi,TMx=',Mfile,Yfile,Obsi,C2Max/30.0 X2 = X1 + n2 Mfile = Mfile + 1 if (Mfile > 12) then call AdjustMONTHandYEAR(Mfile, Yfile) end if call ReadMonth(C3Min, C3Max, fhandle, rc) +! write(*,*)'TTT3a,Mfile,Y,Obsi,TMx=',Mfile,Yfile,Obsi,C3Max/30.0 X3 = X2 + n3 OK3 = .true. end if @@ -765,6 +772,7 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & if ((.not. OK3) .and. (Monthi == GetTemperatureRecord_ToM())) then if ((GetTemperatureRecord_FromY() == 1901) & .or. (Yeari == GetTemperatureRecord_ToY())) then +! write(*,*)'TTT4' do Nri = 1, (GetTemperatureRecord_NrObs()-3) read(fhandle, *, iostat=rc) Mfile = Mfile + 1 @@ -793,6 +801,7 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & ! 5. IF not previous cases if (.not. OK3) then +! write(*,*)'TTT5,Monthi,Yeari=',Monthi,Yeari Obsi = 1 do while (.not. OK3) if ((Monthi == Mfile) .and. (Yeari == Yfile)) then @@ -802,7 +811,7 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & if (Mfile > 12) then call AdjustMONTHandYEAR(Mfile, Yfile) end if - Obsi = Obsi + 1 + Obsi = Obsi + 1 end if end do Mfile = GetTemperatureRecord_FromM() @@ -814,6 +823,7 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & end if end do call ReadMonth(C1Min, C1Max, fhandle, rc) +! write(*,*)'TTT5a,Mfile,Y,Obsi,TMx=',Mfile,Yfile,Obsi,C1Max/30.0 X1 = n1 t1 = X1 Mfile = Mfile + 1 @@ -821,15 +831,18 @@ subroutine GetSetofThreeMonths(Monthi, Yeari, & call AdjustMONTHandYEAR(Mfile, Yfile) end if call ReadMonth(C2Min, C2Max, fhandle, rc) +! write(*,*)'TTT5b,Mfile,Y,Obsi,Tmx=',Mfile,Yfile,Obsi,C2Max/30.0 X2 = X1 + n2 Mfile = Mfile + 1 if (Mfile > 12) then call AdjustMONTHandYEAR(Mfile, Yfile) end if call ReadMonth(C3Min, C3Max, fhandle, rc) +! write(*,*)'TTT5c,Mfile,Y,Obsi,Tmx=',Mfile,Yfile,Obsi,C3Max/30.0 X3 = X2 + n3 end if +! write(*,*)'TTT6' close(fhandle) end subroutine GetSetofThreeMonths @@ -1048,12 +1061,16 @@ integer(int32) function SumCalendarDays(ValGDDays, FirstDayCrop, Tbase, Tupper,& type(rep_DayEventDbl), dimension(31) :: TminDataSet, TmaxDataSet logical :: AdjustDayNri real(dp) :: TDayMin_loc, TDayMax_loc + logical :: debug = .FALSE. TDayMin_loc = TDayMin TDayMax_loc = TDayMax NrCdays = 0 + if(debug)then + write(*,*)'SCD DtoH()=',NrCDays + end if if (ValGDDays > 0) then if (GetTemperatureFile() == '(None)') then ! given average Tmin and Tmax @@ -1082,6 +1099,9 @@ integer(int32) function SumCalendarDays(ValGDDays, FirstDayCrop, Tbase, Tupper,& select case (GetTemperatureRecord_DataType()) case (datatype_daily) + if(debug)then + write(*,*)'SCD0 DtoH()=',NrCDays + end if ! Tmin and Tmax arrays contain the TemperatureFilefull data i = DayNri - GetTemperatureRecord_FromDayNr() + 1 TDayMin_loc = Tmin(i) @@ -1155,6 +1175,9 @@ integer(int32) function SumCalendarDays(ValGDDays, FirstDayCrop, Tbase, Tupper,& end if case(datatype_monthly) + if(debug)then + write(*,*)'SCD1 DtoH()=',NrCDays + end if call GetMonthlyTemperatureDataSet(DayNri, & TminDataSet, TmaxDataSet) i = 1 @@ -1168,6 +1191,10 @@ integer(int32) function SumCalendarDays(ValGDDays, FirstDayCrop, Tbase, Tupper,& NrCDays = NrCDays + 1 RemainingGDDays = RemainingGDDays - DayGDD DayNri = DayNri + 1 + if(debug)then + write(*,*)'SCD2 DtoH()=',NrCDays + write(*,*)'SCD2b ',RemainingGDDays,DayNri,GetTemperatureRecord_ToDayNr(),AdjustDayNri + end if do while ((RemainingGDDays > 0) & .and. ((DayNri < GetTemperatureRecord_ToDayNr()) & .or. AdjustDayNri)) @@ -1186,16 +1213,28 @@ integer(int32) function SumCalendarDays(ValGDDays, FirstDayCrop, Tbase, Tupper,& NrCDays = NrCDays + 1 RemainingGDDays = RemainingGDDays - DayGDD DayNri = DayNri + 1 + if(debug)then + write(*,*)'SCD2c ',RemainingGDDays,DayGDD,DayNri,Tbase,Tupper,TdayMin_loc,TDayMax_loc + end if end do + if(debug)then + write(*,*)'SCD3 DtoH()=',NrCDays + end if if (RemainingGDDays > 0) then NrCDays = undef_int end if + if(debug)then + write(*,*)'SCD4 DtoH()=',NrCDays + end if end select else NrCDays = undef_int endif endif endif + if(debug)then + write(*,*)'SCD5 DtoH()=',NrCDays + end if SumCalendarDays = NrCDays end function SumCalendarDays @@ -1371,6 +1410,7 @@ subroutine AdjustCalendarDays(PlantDayNr, InfoCropType,& real(dp) :: tmp_NoTempFileTMin, tmp_NoTempFileTMax integer :: ExtraDays, ExtraGDDays + logical :: debug = .FALSE. tmp_NoTempFileTMin = NoTempFileTMin tmp_NoTempFileTMax = NoTempFileTMax @@ -1398,6 +1438,9 @@ subroutine AdjustCalendarDays(PlantDayNr, InfoCropType,& Tbase, Tupper, tmp_NoTempFileTMin, tmp_NoTempFileTMax) DHarvest = SumCalendarDays(GDDHarvest, PlantDayNr,& Tbase, Tupper, tmp_NoTempFileTMin, tmp_NoTempFileTMax) + if(debug)then + write(*,*)'ZZ0 DHarvest=',DHarvest + end if end if DLZmax = SumCalendarDays(GDDLZmax, PlantDayNr,& @@ -1433,6 +1476,9 @@ subroutine AdjustCalendarDays(PlantDayNr, InfoCropType,& Succes = .false. end if + if(debug)then + write(*,*)'ZZ1 DHarvest=',DHarvest + end if if (Succes) then CGC = (real(GDDL12, kind=dp)/real(D12, kind=dp)) * GDDCGC call GDDCDCToCDC(PlantDayNr, D123, GDDL123, GDDHarvest,& @@ -1440,6 +1486,9 @@ subroutine AdjustCalendarDays(PlantDayNr, InfoCropType,& call DetermineLengthGrowthStages(CCo, CCx, CDC, D0, DHarvest,& IsCGCGiven, TheDaysToCCini, & ThePlanting, D123, StLength, D12, CGC) + if(debug)then + write(*,*)'ZZ2 DHarvest=',DHarvest + end if if ((InfoCropType == subkind_Grain) .or. (InfoCropType == subkind_Tuber)) then dHIdt = real(HIndex, kind=dp)/real(LHImax, kind=dp) end if @@ -1482,6 +1531,7 @@ subroutine AdjustCalendarCrop(FirstCropDay) real(dp) :: Crop_CGC_temp real(dp) :: Crop_CDC_temp real(dp) :: Crop_dHIdt_temp + logical :: debug = .FALSE. CGCisGiven = .true. @@ -1500,6 +1550,9 @@ subroutine AdjustCalendarCrop(FirstCropDay) Crop_DaysToFlowering_temp = GetCrop_DaysToFlowering() Crop_LengthFlowering_temp = GetCrop_LengthFlowering() Crop_DaysToSenescence_temp = GetCrop_DaysToSenescence() + if(debug)then + write(*,*)'ACC0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if Crop_DaysToHarvest_temp = GetCrop_DaysToHarvest() Crop_DaysToMaxRooting_temp = GetCrop_DaysToMaxRooting() Crop_DaysToHIo_temp = GetCrop_DaysToHIo() @@ -1956,7 +2009,11 @@ subroutine LoadSimulationRunProject(NrRun) integer(int32) :: ZiAqua_temp type(rep_clim) :: etorecord_tmp, rainrecord_tmp real(dp) :: ECiAqua_temp, SurfaceStorage_temp + logical :: debug = .FALSE. + if(debug)then + write(*,*)'WW0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 0. Year of cultivation and Simulation and Cropping period call SetSimulation_YearSeason(ProjectInput(NrRun)%Simulation_YearSeason) call SetCrop_Day1(ProjectInput(NrRun)%Crop_Day1) @@ -1976,6 +2033,9 @@ subroutine LoadSimulationRunProject(NrRun) call SetClimateDescription(trim(TempString)) close(fClim) end if + if(debug)then + write(*,*)'WW1 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 1.1 Temperature call SetTemperatureFile(ProjectInput(NrRun)%Temperature_Filename) @@ -2002,6 +2062,9 @@ subroutine LoadSimulationRunProject(NrRun) call SetTemperatureRecord(temperature_record) end if + if(debug)then + write(*,*)'WW2 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 1.2 ETo call SetEToFile(ProjectInput(NrRun)%ETo_Filename) if ((GetEToFile() == '(None)') .or. & @@ -2018,6 +2081,9 @@ subroutine LoadSimulationRunProject(NrRun) call SetEToRecord(etorecord_tmp) end if + if(debug)then + write(*,*)'WW3 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 1.3 Rain call SetRainFile(ProjectInput(NrRun)%Rain_Filename) if ((GetRainFile() == '(None)') .or. & @@ -2035,6 +2101,9 @@ subroutine LoadSimulationRunProject(NrRun) call SetRainRecord(rainrecord_tmp) end if + if(debug)then + write(*,*)'WW4 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 1.4 CO2 call SetCO2File(ProjectInput(NrRun)%CO2_Filename) if (GetCO2File() /= '(None)') then @@ -2049,6 +2118,9 @@ subroutine LoadSimulationRunProject(NrRun) end if call AdjustOnsetSearchPeriod() ! Set initial StartSearch and StopSearchDayNr + if(debug)then + write(*,*)'WW5 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 2. Calendar call SetCalendarFile(trim(ProjectInput(NrRun)%Calendar_Filename)) if (GetCalendarFile() == '(None)') then @@ -2061,12 +2133,24 @@ subroutine LoadSimulationRunProject(NrRun) call SetCalendarDescription(CalendarDescriptionLocal) end if + if(debug)then + write(*,*)'WW5d GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 3. Crop call SetSimulation_LinkCropToSimPeriod(.true.) + if(debug)then + write(*,*)'WW5e GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call SetCropFile(ProjectInput(NrRun)%Crop_Filename) call SetCropFilefull(ProjectInput(NrRun)%Crop_Directory // GetCropFile()) + if(debug)then + write(*,*)'WW5f GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call LoadCrop(GetCropFilefull()) + if(debug)then + write(*,*)'WW6 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! Adjust crop parameters of Perennials if (GetCrop_subkind() == subkind_Forage) then ! adjust crop characteristics to the Year (Seeding/Planting or @@ -2096,19 +2180,37 @@ subroutine LoadSimulationRunProject(NrRun) Crop_DaysToHarvest_temp = GetCrop_DaysToHarvest() Crop_GDDaysToSenescence_temp = GetCrop_GDDaysToSenescence() Crop_GDDaysToHarvest_temp = GetCrop_GDDaysToHarvest() + if(debug)then + write(*,*)'WW6.0 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call AdjustCropFileParameters(GetCropFileSet(),& GetCrop_DaysToHarvest(), GetCrop_Day1(), & GetCrop_ModeCycle(), GetCrop_Tbase(), GetCrop_Tupper(),& Crop_DaysToSenescence_temp, Crop_DaysToHarvest_temp,& Crop_GDDaysToSenescence_temp, Crop_GDDaysToHarvest_temp) call SetCrop_DaysToSenescence(Crop_DaysToSenescence_temp) + if(debug)then + write(*,*)'WW6,1 Crop_DaysToHarvest_temp=',Crop_DaysToHarvest_temp + end if call SetCrop_DaysToHarvest(Crop_DaysToHarvest_temp) + if(debug)then + write(*,*)'WW6.2 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call SetCrop_GDDaysToSenescence(Crop_GDDaysToSenescence_temp) call SetCrop_GDDaysToHarvest(Crop_GDDaysToHarvest_temp) end if + if(debug)then + write(*,*)'WW6a GetCrop_DaysToHarv,GetCrop_Day1=',GetCrop_DaysToHarvest(),GetCrop_Day1() + end if call AdjustCalendarCrop(GetCrop_Day1()) + if(debug)then + write(*,*)'WW6a2 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if call CompleteCropDescription + if(debug)then + write(*,*)'WW6a3 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! Onset.Off := true; if (GetClimFile() == '(None)') then Crop_Day1_temp = GetCrop_Day1() @@ -2120,6 +2222,9 @@ subroutine LoadSimulationRunProject(NrRun) else call SetCrop_DayN(GetCrop_Day1() + GetCrop_DaysToHarvest() - 1) end if + if(debug)then + write(*,*)'WW6b GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! adjusting ClimRecord.'TO' for undefined year with 365 days if ((GetClimFile() /= '(None)') .and. (GetClimRecord_FromY() == 1901) & @@ -2128,6 +2233,9 @@ subroutine LoadSimulationRunProject(NrRun) end if ! adjusting simulation period call AdjustSimPeriod + if(debug)then + write(*,*)'WW6c GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 4. Irrigation call SetIrriFile(ProjectInput(NrRun)%Irrigation_Filename) @@ -2140,6 +2248,9 @@ subroutine LoadSimulationRunProject(NrRun) // GetIrriFile()) call LoadIrriScheduleInfo(GetIrriFileFull()) end if + if(debug)then + write(*,*)'WW6d GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 5. Field Management call SetManFile(ProjectInput(NrRun)%Management_Filename) @@ -2167,6 +2278,9 @@ subroutine LoadSimulationRunProject(NrRun) call SetSimulation_EffectStress_RedCCX(RedCCX_temp) end if + if(debug)then + write(*,*)'WW7 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 6. Soil Profile call SetProfFile(ProjectInput(NrRun)%Soil_Filename) if (GetProfFile() == '(External)') then @@ -2192,6 +2306,9 @@ subroutine LoadSimulationRunProject(NrRun) ! Loading the groundwater is done after loading the soil profile (see ! 9.) end if + if(debug)then + write(*,*)'WW8 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 8. Set simulation period call SetSimulation_FromDayNr(ProjectInput(NrRun)%Simulation_DayNr1) @@ -2318,6 +2435,9 @@ subroutine LoadSimulationRunProject(NrRun) call ResetSWCToFC() end if + if(debug)then + write(*,*)'WW9 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if ! 11. Off-season conditions call SetOffSeasonFile(ProjectInput(NrRun)%OffSeason_Filename) if (GetOffSeasonFile() == '(None)') then @@ -2341,6 +2461,9 @@ subroutine LoadSimulationRunProject(NrRun) call GetFileDescription(GetObservationsFileFull(), observations_descr) call SetObservationsDescription(observations_descr) end if + if(debug)then + write(*,*)'WW10 GetCrop_DaysToHarvest()=',GetCrop_DaysToHarvest() + end if contains @@ -3006,6 +3129,7 @@ real(dp) function BiomassRatio(TempDaysToCCini, TempGDDaysToCCini,& real(dp) :: RatDGDD, SumBPot, SumBSF integer(int32) :: tSwitch, DaysYieldFormation + write(*,*)'AA tempproc GetCrop_DaysToHarvest()=',TempHarvest ! 1. Initialize ! 1 - a. Maximum sum Kc SumKcTop = SeasonalSumOfKcPot(TempDaysToCCini, TempGDDaysToCCini,& @@ -3144,6 +3268,7 @@ subroutine StressBiomassRelationship(TheDaysToCCini, TheGDDaysToCCini,& call SetSimulation_DelayedDays(0) ! required for CalculateETpot L12SF = L12 ! to calculate SumKcTop (no stress) GDDL12SF = GDDL12 ! to calculate SumKcTop (no stress) + write(*,*)'BB tempproc GetCrop_DaysToHarvest()=',L1234 ! Maximum sum Kc (no stress) SumKcTop = SeasonalSumOfKcPot(TheDaysToCCini, TheGDDaysToCCini,& L0, L12, L123, L1234, GDDL0, GDDL12, GDDL123, GDDL1234,& @@ -3356,6 +3481,7 @@ subroutine CCxSaltStressRelationship(TheDaysToCCini, TheGDDaysToCCini,& call SetSimulation_DelayedDays(0) ! required for CalculateETpot GDDL12SS = GDDL12 ! to calculate SumKcTop (no stress) BNor100 = real(undef_int, kind=dp) + write(*,*)'CC tempproc GetCrop_DaysToHarvest()=',L1234 ! Maximum sum Kc (no stress) SumKcTop = SeasonalSumOfKcPot(TheDaysToCCini, TheGDDaysToCCini,& L0, L12, L123, L1234, GDDL0, GDDL12, GDDL123, GDDL1234,& diff --git a/src/utils.f90 b/src/utils.f90 index e3eb690e..e215712f 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -54,7 +54,7 @@ function GetAquaCropDescriptionWithTimeStamp() result(str) WRITE(timestr,8)d(5),d(6),d(7) 10 FORMAT(I2.2, '-', I2.2, '-', I4.4) 8 FORMAT(I2.2, ':', I2.2, ':', I2.2) - str = GetAquaCropDescription() // ' - Output created on (date) : ' // & + str = '# ' // GetAquaCropDescription() // ' - Output created on (date) : ' // & datestr // ' at (time) : ' // timestr end function GetAquaCropDescriptionWithTimeStamp