KUL-RSDA / AquaCrop

AquaCrop source code endorsed by FAO
https://www.fao.org/aquacrop/en/
Other
56 stars 17 forks source link

several bugs when working with monthly climate data #334

Closed mcguirepatr closed 1 year ago

mcguirepatr commented 1 year ago

Hi AquaCroppers: I can't get my monthly climate driving data working with the standalone Windows version of AquaCrop. It seg faults. I was able to get it working with the GitHub version on Linux, but I had to fix 3 different bugs, which I show below. Can you confirm that these are bugs? Can you merge these changes with the GitHub version of AquaCrop? Can you also fix it in the standalone Windows version? I would like to also be able to use the Windows version with monthly data.

I have tried several different monthly climate data sets, including 1 year for Lima, Peru, (provided with the standalone Windows version) normally just using the Precip PLU and Temperature Tnx files, without the ETo files. But I have also tried it when also including the ETo data. And I have tried it with 130 years of data that I processed. Patrick McGuire University of Reading, UK

svn diff
===================================================================
--- run.f90     (revision 1468)
+++ run.f90     (working copy)
@@ -5707,6 +5727,7 @@
     ! 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())
@@ -5909,7 +5961,7 @@
                     if (RunningDay > GetTminDataSet_DayNr(31)) then
                         TminDataSet_temp = GetTminDataSet()
                         TmaxDataSet_temp = GetTmaxDataSet()
-                        call GetMonthlyTemperatureDataSet(FromSimDay, &
+                        call GetMonthlyTemperatureDataSet(RunningDay, &
                                                           TminDataSet_temp, &
                                                           TmaxDataSet_temp)
                         call SetTminDataSet(TminDataSet_temp)

===================================================================
--- climprocessing.f90  (revision 1468)
+++ climprocessing.f90  (working copy)
@@ -311,10 +323,14 @@
             Mfile = Mfile + 1
             if (Mfile > 12) then
                 call AdjustMONTHandYEAR(Mfile, Yfile)
+            else
+                read(fETo, *) C3
+                C3 = C3 * ni
+                X3 = X2 + n3
             end if
-            read(fETo, *) C3
-            C3 = C3 * ni
-            X3 = X2 + n3
         end if

         close(fETo)
@@ -749,10 +784,17 @@
                     end if
                     Obsi = Obsi + 1
                 end if
                 if (OK3) exit loop
             end do loop
             Mfile = GetRainRecord_FromM()
-            do Nri = 1, (Obsi-2)
+            do Nri = 1, (Obsi-3)
                 read(fRain, *)
                 Mfile = Mfile + 1
                 if (Mfile > 12) then
mcguirepatr commented 1 year ago

As requested, here are some settings which would trigger some of the errors. The PLU, ETo, and Tnx files are attached as Txt files.

      1         : Year number of cultivation (Seeding/planting year)
    366         : First day of simulation period - 1 Jan 1902
   1460         : Last day of simulation period - 31 December 1905
    420         : First day of cropping period 
    460         : Last day of cropping period 
-- 1. Climate (CLI) file
   (None)
   (None)
   1.1 Temperature (Tnx or TMP) file
   Lima4.Tnx
   '../../../GUI_AC71/AquaCropV71No17072
   1.2 Reference ET (ETo) file
   Lima4.ETo
   '../../../GUI_AC71/AquaCropV71No17072023/DATA/'
   1.3 Rain (PLU) file
   Lima4.PLU
   '../../../GUI_AC71/AquaCropV71No17072023/DATA/'
   1.4 Atmospheric CO2 concentration (CO2) file
   MaunaLoa.CO2
   './SIMUL/'

Lima4_Tnx.txt Lima4_PLU.txt Lima4_ETo.txt

KUL-RSDA commented 1 year ago

@mcguirepatr We much appreciate your bug report. We confirm that it's a bug originating from the conversion from Pascal to Fortran. The monthly input data feature was unfortunately not in our set of test cases for ensuring zero-diff.
We opened a pull request that fixes this bug #335 We are running some last tests and will merge the PR in the coming days. Then we will replace the various compiled standalone programs in the release v7.1 together with a post-release bug fix note. The download links at the FAO website will then directly refer to the updated v7.1 files.

mcguirepatr commented 1 year ago

Thanks for examining my bug report, and for the rapid turnaround. Patrick McGuire

mbechtold commented 1 year ago

All executables (linux, windows, macos) have been replaced with the monthly climate reading bug fix. Thanks again @mcguirepatr for pointing out this issue!

mcguirepatr commented 1 year ago

Hi mbechtold: Thanks for fixing this. I can confirm that it's mostly working on linux with the GitHub version. But when I use only 12 months of monthly Tnx & PLU driving data with either the linux/GitHub version or the windows standalone version, then it crashes at either line 763, 764, or 733 in climprocessing.f90 . Can you fix this?

For example, line 763 is read(fRain, *) C1

It seems that you didn't need to update the Windows GUI version. I have now figured out why my monthly driving data wasn't working with the Windows GUI version: I had neglected to convert the PLU & Tnx ASCII text files from linux to windows format. I have now done that with the linux command unix2dos, and things are better now, and I can visualize the driving data in the AquaCrop GUI in windows. Patrick

mbechtold commented 1 year ago

Hi @mcguirepatr sorry that it apparently does not work smoothly for you. We were not yet able to reproduce your issue. It works for us even when there is only 12 month of input climate data. Can you check whether this example (only 12 month data) works for you as it did for us: 12month_example.zip If yes, can you share your monthly files to see what causes the problem? Thx.

mcguirepatr commented 1 year ago

Thanks, mbechtold I tried with your setup, and it works fine. But when I change this in LIST/TunisMonth.PRM:

                          81         : First day of simulation period - 21 May 2014
                         205        : Last day of simulation period - 31 October 2014

to this:

                            1         : First day of simulation period - 1 January 2014
                         365        : Last day of simulation period - 31 Decober 2014

then it crashes in line 315 of climprocessing.f90

Maybe some of the following changes are the right thing to do?:

#line 293:
            do Nri = 1, (Obsi-2)
#new line 293:
            do Nri = 1, (Obsi-3) #only skip the first 9 months, not 10 months (in the case of a 12-month run)
#line 729:
             do Nri = 1, (GetRainRecord_NrObs()-2)
#new line 729
             do Nri = 1, (GetRainRecord_NrObs()-3) 
#line 755:
             do Nri = 1, (Obsi-2)
#new line 755:
             do Nri = 1, (Obsi-3)

Patrick

mbechtold commented 1 year ago

Thanks to you again @mcguirepatr We solved the issue related to a simulation time that is ending exactly with the length of the climate data. The problem was sitting somewhere else, see PR #338 The 'main' branch is updated, and the linux and windows executables are replaced (macos still following).

mcguirepatr commented 1 year ago

Thanks. I just did an 'svn update' of my local copy of the code from the github code. The change that you made in global.f90 now allows me to simulate all 365 days from a 12 month dataset of driving data (using the TunisMonth.PRM data set).

mcguirepatr commented 1 year ago

Hi Michel et al. @mbechtold et al.: In the climprocessing.f90 code that I referred you to (see below for a copy), it looks like, in the case of Obsi=12 (is that my case here of 12 months of data?), that it skips the first 12-2=10 observations, and then reads the next 3. So I am a little confused, since it reads 13 times, which is greater than the length of the file. See the code below. I am worried that the month-reading is one month off (ahead or behind). Can you explain this so that I understand better? Thanks, Patrick

            do Nri = 1, (Obsi-2)
                read(fRain, *)
                Mfile = Mfile + 1
                if (Mfile > 12) then
                    call AdjustMONTHandYEAR(Mfile, Yfile)
                end if
            end do
            read(fRain, *) C1
            read(fRain, *) C2
            read(fRain, *) C3
mbechtold commented 1 year ago

@mcguirepatr The "Mfile>12" is not related to your 12 months of data case but to the end of a year. We fixed the 12-month issue in another part of the code (see response above). The code you share here is in a function called GetSetofThreeMonths. The reader was programmed many years ago. It's not the most intuitive code, but there was likely a good reason for how it was programmed. We converted code almost a 1:1 conversion from Pascal to Fortran, to ensure step-by-step the zero-difference in simulation output. Now with all having in Fortran, we have the intention to restructure some code in the coming years if time and funding allow. For now, however, I can confirm that the current code reads in the 12-month data record without a shift like you were worried about. You can verify yourself by comparing the PLU file with the 'rain variable' in the output in OUTP/PRMday.OUT

mcguirepatr commented 1 year ago

Thanks for the feedback and for the suggestion for how to check for a shift. I followed your suggestion, and I made a quick plot. I think I plotted it approximately correctly. And I don't see any shift. The code seems to be doing it properly. See the plot below. Patrick

Screenshot 2023-10-16 at 20 59 21
mcguirepatr commented 1 year ago

I do wonder now, however, why (for example) the daily-precip data points for days 50-60 are so much higher than any scaled monthly data points during that time period. I thought the daily precip was just interpolated from the monthly data, so it is a bit puzzling how the daily precip data can stick out so much from the monthly trend. Or maybe I am misplotting this?

mcguirepatr commented 1 year ago

Or is that time period from days 50-60 the time period of the short month at the end of February, so my scaling is a bit off?