PLAPS documentation wrong

[marek...@gmail.com]

Hi,

I found a paper from 2016 pointing out that PLAPS was implemented as mm/m, not mm/km as the documentation states. I find the same in my models, getting extreme precipitation if I assume mm/km. Is this really still an issue that hasn't been fixed? 

It's odd, because I see other people using PLAPS as mm/km, for example Karim Abasspour in his tutorial: https://www.youtube.com/watch?v=Y7-u7Ae68UA (45:15). 

Can anyone confirm how they use PLAPS? I am using rev681 and it's definitely using mm/m. 

Best regards,

Bendik

[Steven Jepsen]

Bendik,

I know that the effective precipitation lapse rate is multiplied by a factor that looks something like 365/(# precipitation days per year), or something like that (see THEORY manual). So if PLAPS is 100 mm/km, and there were 182.5 wet days per year, the effective precipitation lapse rate would be 100 mm/km * 2 = 200 mm/km for that particular day. That is how I think of it. 

Steve 

[Steven Jepsen]

Take two!... Guess it pays to check the documentation before posting. PLAPS is in mm/km. The eqn. on p. 91 of the theory manual is:

R_band = R_day + (EL_band - EL_gauge)*(PLAPS)/(days_prec_per_yr * 1000)

where R_band = precipitation in elevation band (mm), R_day = precipitation at gauge (mm), PLAPS = precipitation lapse rate (mm/km), days_prec_per_yr = average # days per year with precipitation, and EL_band and EL_gauge are elevations of elevation band (m) and gauge, respectively. The 1000 in the above equation is needed because PLAPS is in units of mm/km while elevation is in units of m.

Steve

On Thursday, August 27, 2020 at 5:39:20 AM UTC-7, marek...@gmail.com wrote:

[marek...@gmail.com]

Thank you Steve! Now the process makes a lot more sense to me. I should have checked more thoroughly, instead of just reading the file.cio manual. 

An example to show what I expect:

R_day = 40mm

EL_band = 500m

El_gauge = 300m

PLAPS = 20mm/km

Days_prec_per_year = 200

So R_band = 40mm + (500m-300m)*((20mm/km)/(200))*(km/1000m)

R_band = 40mm + 200*20mm/(200*1000) = 40.02mm

However, I still have some serious issues. In the two images I attached, I show my clean simulation without any PLAPS or calibration other than making elevation bands, as well as my results when I set PLAPS = 20 (which, as far as I understand, adds/subtracts 20mm to the yearly precipitation per km the elevation band is above/below the gauge, and scales the daily P accordingly). I intentionally used only one gauge for this example, and the gauge is higher than most of my elevation bands.  It is quite obvious that SWAT does not read the PLAPS correctly, and I think it most likely fails to divide by 1000 somewhere. 

The precipitation on day 1 (for one sub) for example: 0.48mm without plaps, 138mm with plaps = 20. If I set plaps=0.02, the precipitation becomes 0.54mm.

Are you able to change PLAPS as mm/km and get reasonable results?

Best regards,

Bendik

[Steven Jepsen]

Bendik,

How strange. Looks like you are taking a good approach to figure this out.

I've always used gridded data with a "gauge" strategically placed at the mean elevation of each subbasin, so PLAPS has never significantly affected precipitation volume in my subbasins. However, I have noticed some small discrepancies in the past compared to what I expected.

OK, so I get that R_band on day 1 of your example would be 0.50 mm if the 1000 was included in the equation denominator, 20.48 mm if 1000 was not included. That assumed (correctly?) that EL_band - EL_gauge = 200 m, PLAPS = 20 mm/km, and days_prec_per_yr = 200. Do you get this, too? You listed a different value, 138 mm, which would suggest to me that either my assumed values are incorrect or that the missing denominator is not the culprit. However, maybe I assumed incorrect components going into the equation for R_band.

Are you able to confirm each of those equation components in your SWAT model? (i.e., elevation difference, # wet days per year, etc.). If any of those are way off, possibly as a results of a glitch somewhere, that could throw your lapsed precipitation values way off. Maybe also verify that the correct precipitation gauge number, IRGAGE, is listed in the *.sub file of your subbasin. You might also play around with the mean elevation listed in your subbasin *.sub file to see what's going on, starting with the elevation of the gauge and then gradually increasing it. Setting elevations of different bands to the mean elevation of the subbasin would be another simplifying thing to try.

I'd be interested to hear what you find out.

Thanks,

Steve

[marek...@gmail.com]

Hi Steve,

The example was not related to the outputs from my model, so those numbers shouldn't match up. My bad, I see how that was misleading. 

Good idea, I will set the elevation bands in the subbasin to all be the same elevation, then I can compare with some hand-calculations to see if it is simply missing a 1/1000. I think most bands get 0 precipitation (the gauge is higher than 4/5 elevation bands in that sub), but the remaining one gets a tremendous amount. 

I'm not sure if I can verify the number of wet days, but I can calculate that from my precipitation data and do a test like described above. I will let you know once I've done that. 

Bendik

[Steven Jepsen]

Bendik,

Sounds good. I'm not sure about this, but the value of parameter days_prec_per_yr  may come from line 14 of the weather file (*.wgn) of the applicable subbasin, e.g. 000010000.wgn for subbasin #1. Parameter days_prec_per_yr may be calculated as follows:

days_prec_per_yr = sum_{i}[PCPD_{i}]

where PCPD_{i} = average number of wet days during month i.

It's possible that SWAT computes the value of parameter days_prec_per_yr from your actual precipitation forcing data, but I kind of doubt it. 

Steve  

[marek...@gmail.com]

Hi Steve,

I think I found the source of my problems (but not a solution). My colleague read the paper  "Rainfall estimation in SWAT: An alternative method to simulate orographic precipitation" (https://www.sciencedirect.com/science/article/pii/S0022169413008688) more thoroughly and noticed that they also commented on the number of rainy days. 

Apparently SWAT uses a constant of "12" for number of rainy days, and does not divide by 1000 in the equation. With these assumptions, my numbers work out perfectly.

I edited my files to test it:

Elev_gauge = 1330

Elev bands = 1430 (all the same)

R_day = 0.496mm

PLAPS = 20

Assume: days_prec_per_year = 12, and do not divide by 1000

Results from equation: R_band = 0.496+100*(20/12) = 167mm

Results from SWAT: R_band = 166mm

I have no idea where to change this number of rainy days parameter. I can't find it anywhere in my input files, and in the absolute swat values file it says it is in the .wgn files. The .wgn files have different numbers (around 80), and since I am not using and wgens I think they aren't used anyway. What I really don't understand is how this problem can persist when the model is so widely utilized. How come there aren't more people facing this issue?

Best regards,

Bendik

[Steven Jepsen]

Hi Bendik,

The number of rainy days per year is supposed to come from summing monthly average number of rainy days in the the .wgn file, according to p. 93 of the 2009 Theory Manual. I take this to be the case regardless of whether your using WGENX for your met forcings.

Isn't the format of precipitation f5.1? If so, how would it possible to input a value of 0.496 mm? Could precipitation be getting read incorrectly? I'd double-check the format of your precipitation inputs to make sure they comply with the I/O manual. 

I take the PLAPS to mean the change in yearly precipitation per km elevation, with units of mm yr^-1 km^-1. The manual is a little unclear on the year part of this. The values that you stated don't make sense to me. If PLAPS = 20, then annual precipitation should increase by 20 mm per km elevation, right? But you found that precipitation during one day increased 166 mm over 100 m? Something is wrong somewhere, as you've recognized. The only thing I'd suggest is to double-check all of the values going into SWAT, including precipitation format and anything related to time (monthly versus daily values, etc). As a last resort, maybe check what's going on in the source code? Let me know if there is something I can do to help. I'd probably start by simplifying things, as you've done. 

Steve

[marek...@gmail.com]

Hi Steve, 

Good catch about the format of the precipitation. I read the "R_day" from output.sub before applying any PLAPS, since it SHOULD just be the measured precipitation from its chosen gauging station when PLAPS=0. However, it isn't! This might be another clue to what is going wrong. My input pcp is 0.48, but without any PLAPS SWAT gives me 0.496 (in output.sub).

I tried changing the number of rainy days in the .wgn file for the subcatchment, but I get no change in the precipitation output. This supports my theory that SWAT uses a constant 12 for number of rainy days, defined in some mysterious file that I can't find. It is possible that the Fortran code uses a default "1" per month if no weather generators are used...

I have contacted the authors of the previously mentioned paper for advice on how they managed to fix it, and I will contact the SWAT support team now. I would read the source code, but I can't make heads or tails of the Fortran code. 

Bendik

[Steven Jepsen]

Bendik,

I think that the precipitation output in output.sub has been summed over all elevation bands and then normalized by subbasin area, and that it DOES include the effects of precipitation lapse rate. It SHOULD NOT BE the same as the input gage-value of precipitation unless PLAPS=0. Same is true for air temperature.

The input format of precipitation is f5.1, so how you could input a value of 0.48? I ask just to make sure SWAT is reading your precipitation correctly. Maybe you already verified that in your output.sub... 

Where are you getting your gage elevation? I think it should be from the *.pcp forcing file. You can find the precip gauge # in *.sub and then using that gauge # to read the gage elevation off the *.pcp file. 

I'd be surprised if a "12" was being included in the denominator and the factor of 1000 was being left out, unless SWAT redefined precipitation lapse rate and its units since the 2009 theory manual. I'd be happy to check out some of your input files if you'd like a second set of eyes on it, maybe using Google drive?   

Steve

[marek...@gmail.com]

Hi Steve,

That is the problem. PLAPS=0 but the precipitation in output.sub is not exactly the same as in pcp.pcp. It's not even consistent between subbasins, as shown below in the output.sub file. 

You're right though, 0.48 is in my precipitation gauge file, not pcp.pcp. In pcp.pcp it is 000.5 since it was rounded when creating the SWAT files (so it shouldn't be 0.496 in output.sub when PLAPS=0, it should be 0.500). 

Gauge elevation is from pcp.pcp, 1330 for this gauge. 

I would also be surprised at the 12 and 1000 thing if the exact same issue wasn't reported by another user in a published paper. The fact that the days of precipitation numbers in .wgn have no influence on the outputs of the model (tested) is a clear indicator that SWAT reads the number of rainy days from elsewhere, or assumes a number (apparently 12). 

Moreover when I set the elevations in the subbasin lower than the gauge and use PLAPS=20, precipitation in that sub is reduced to 0. This show that SWAT does try to to increase/reduce the precipitation based on elevation, but is off by orders of magnitude. 

I'm not at liberty to share the data, unfortunately. The provider of the data was pretty clear that we can't share it externally. I can post the first few headings of each file, though:

pcp.pcp (I hid the coordinates, but the x's were numbers): 

Station  station1, station2 , station3 , station4 ,

Lati    xx.9 xx.8 xx.6 xx.6

Long    xx.2 xx.4 xx.6 xx.0

Elev     378  972 1330  972

2016001000.0000.0000.5000.3

2016002000.2000.1000.6000.0

2016003049.3058.8029.8017.0

2016004022.5024.0019.8007.1

2016005022.7009.9006.0001.8

2016006050.5055.4028.4017.9

000050000.sub (all elevation bands set to 1430 for testing, but I get the same error regardless of what they are):

 .sub file Subbasin: 5 8/5/2020 12:00:00 AM ArcSWAT 2012.10_2.19

       42.914269    | SUB_KM : Subbasin area [km2]

Climate in subbasin

       xx.808796    | LATITUDE : Latitude of subbasin [degrees]

         1073.88    | ELEV : Elevation of subbasin [m]

               3    | IRGAGE: precip gage data used in subbasin

               3    | ITGAGE: temp gage data used in subbasin

               1    | ISGAGE: solar radiation gage data used in subbasin

               1    | IHGAGE: relative humidity gage data used in subbasin

               3    | IWGAGE: wind speed gage data used in subbasin

000050000.wgn       | WGNFILE: name of weather generator data file

               1    | FCST_REG: Region number used to assign forecast data to the subbasin

Elevation Bands

| ELEVB: Elevation at center of elevation bands [m]

1430.0001430.0001430.0001430.0001430.000   0.000   0.000   0.000   0.000   0.000

| ELEVB_FR: Fraction of subbasin area within elevation band

   0.148   0.252   0.258   0.219   0.116   0.000   0.000   0.000   0.000   0.000

| SNOEB: Initial snow water content in elevation band [mm]

     0.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0

           0.000    | PLAPS : Precipitation lapse rate [mm/km]

           0.000    | TLAPS : Temperature lapse rate [�C/km]

           0.000    | SNO_SUB : Initial snow water content [mm]

output.sub (notice the bolded values, they are all from the same gauge, and PLAPS is still 0 in the .sub files):

SUB GIS MON AREAkm2 PRECIPmm SNOMELTmm PETmm

BIGSUB 1 0 1.69E+02 0.00E+00 0.00E+00 4.22E-01

BIGSUB 2 0 1.47E+02 0.00E+00 0.00E+00 4.21E-01

BIGSUB 3 0 1.17E+03 0.00E+00 0.00E+00 4.23E-01

BIGSUB 4 0 1.13E+03 0.00E+00 0.00E+00 0.00E+00

BIGSUB 5 0 1.43E+02 4.96E-01 0.00E+00 0.00E+00

BIGSUB 6 0 1.11E+03 0.00E+00 0.00E+00 0.00E+00

BIGSUB 7 0 1.78E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 8 0 1.75E+02 4.86E-01 0.00E+00 0.00E+00

BIGSUB 9 0 1.48E+02 2.98E-01 0.00E+00 0.00E+00

BIGSUB 10 0 1.41E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 11 0 1.91E+02 4.78E-01 0.00E+00 0.00E+00

BIGSUB 12 0 1.48E+02 4.98E-01 0.00E+00 0.00E+00

BIGSUB 13 0 1.53E+02 2.97E-01 0.00E+00 0.00E+00

BIGSUB 14 0 1.38E+03 0.00E+00 0.00E+00 0.00E+00

BIGSUB 15 0 1.12E+03 4.94E-01 0.00E+00 0.00E+00

BIGSUB 16 0 1.36E+03 4.99E-01 0.00E+00 0.00E+00

BIGSUB 17 0 1.25E+03 2.99E-01 0.00E+00 0.00E+00

BIGSUB 18 0 1.15E+02 2.98E-01 0.00E+00 0.00E+00

BIGSUB 19 0 1.41E+02 2.89E-01 0.00E+00 0.00E+00

BIGSUB 20 0 1.14E+01 4.87E-01 0.00E+00 0.00E+00

BIGSUB 21 0 1.15E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 22 0 1.49E+01 0.00E+00 0.00E+00 4.25E-01

BIGSUB 23 0 1.48E+02 2.99E-01 0.00E+00 0.00E+00

BIGSUB 24 0 1.77E+02 2.99E-01 0.00E+00 0.00E+00

BIGSUB 25 0 1.26E+02 4.99E-01 0.00E+00 0.00E+00

BIGSUB 26 0 1.19E+02 5.00E-01 0.00E+00 0.00E+00

BIGSUB 27 0 1.12E+03 4.97E-01 0.00E+00 0.00E+00

BIGSUB 28 0 1.21E+02 4.63E-01 0.00E+00 0.00E+00

BIGSUB 29 0 1.55E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 30 0 1.79E+01 0.00E+00 0.00E+00 0.00E+00

BIGSUB 31 0 1.26E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 32 0 1.50E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 33 0 1.26E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 34 0 1.19E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 35 0 1.26E+02 0.00E+00 0.00E+00 4.27E-01

BIGSUB 36 0 1.11E+03 0.00E+00 0.00E+00 4.25E-01

BIGSUB 37 0 1.11E+03 0.00E+00 0.00E+00 4.24E-01

BIGSUB 38 0 1.28E+02 0.00E+00 0.00E+00 4.24E-01

BIGSUB 39 0 1.16E+03 0.00E+00 0.00E+00 4.25E-01

BIGSUB 1 0 2.69E+02 1.99E-01 0.00E+00 4.06E-01

BIGSUB 2 0 2.47E+02 1.93E-01 0.00E+00 4.03E-01

BIGSUB 3 0 2.17E+03 1.98E-01 0.00E+00 4.07E-01

BIGSUB 4 0 2.13E+03 9.87E-02 0.00E+00 0.00E+00

BIGSUB 5 0 2.43E+02 5.96E-01 0.00E+00 0.00E+00

BIGSUB 6 0 2.11E+03 9.94E-02 0.00E+00 0.00E+00

BIGSUB 7 0 2.78E+02 9.61E-02 0.00E+00 0.00E+00

BIGSUB 8 0 2.75E+02 5.84E-01 0.00E+00 0.00E+00

BIGSUB 9 0 2.48E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 10 0 2.41E+02 9.85E-02 0.00E+00 0.00E+00

BIGSUB 11 0 2.91E+02 5.74E-01 0.00E+00 0.00E+00

BIGSUB 12 0 2.48E+02 5.97E-01 0.00E+00 0.00E+00

BIGSUB 13 0 2.53E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 14 0 2.38E+03 9.96E-02 0.00E+00 0.00E+00

BIGSUB 15 0 2.12E+03 5.93E-01 0.00E+00 0.00E+00

BIGSUB 16 0 2.36E+03 5.98E-01 0.00E+00 0.00E+00

BIGSUB 17 0 2.25E+03 0.00E+00 0.00E+00 0.00E+00

BIGSUB 18 0 2.15E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 19 0 2.41E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 20 0 2.14E+01 5.85E-01 0.00E+00 0.00E+00

BIGSUB 21 0 2.15E+02 9.96E-02 0.00E+00 0.00E+00

BIGSUB 22 0 2.49E+01 1.98E-01 0.00E+00 4.09E-01

BIGSUB 23 0 2.48E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 24 0 2.77E+02 0.00E+00 0.00E+00 0.00E+00

BIGSUB 25 0 2.26E+02 5.99E-01 0.00E+00 0.00E+00

BIGSUB 26 0 2.19E+02 5.99E-01 0.00E+00 0.00E+00

BIGSUB 27 0 2.12E+03 5.96E-01 0.00E+00 0.00E+00

BIGSUB 28 0 2.21E+02 5.56E-01 0.00E+00 0.00E+00

BIGSUB 29 0 2.55E+02 9.92E-02 0.00E+00 0.00E+00

BIGSUB 30 0 2.79E+01 9.82E-02 0.00E+00 0.00E+00

BIGSUB 31 0 2.26E+02 9.91E-02 0.00E+00 0.00E+00

BIGSUB 32 0 2.50E+02 9.89E-02 0.00E+00 0.00E+00

BIGSUB 33 0 2.26E+02 9.81E-02 0.00E+00 0.00E+00

BIGSUB 34 0 2.19E+02 9.77E-02 0.00E+00 0.00E+00

BIGSUB 35 0 2.26E+02 2.00E-01 0.00E+00 4.12E-01

BIGSUB 36 0 2.11E+03 1.99E-01 0.00E+00 4.10E-01

BIGSUB 37 0 2.11E+03 2.00E-01 0.00E+00 4.09E-01

BIGSUB 38 0 2.28E+02 2.00E-01 0.00E+00 4.08E-01

BIGSUB 39 0 2.16E+03 1.99E-01 0.00E+00 4.09E-01

000050000.wgn (these values are nonsense, since I have not set up a wgen. And they are obviously not used, since changing them do not affect the outputs. I have tried changing both line 14 and the elevation listed).

 .Wgn file Subbasin: 5 STATION NAME:Sample 8/5/2020 12:00:00 AM ArcSWAT 2012.10_2.19

  LATITUDE =  33.65 LONGITUDE = -95.69

  ELEV [m] = 179.80

  RAIN_YRS =  83.00

 11.90 13.90 18.80 23.60 27.60 32.30 34.80 34.90 31.40 25.70 18.40 13.10

 -0.20  1.60  5.90 11.10 15.70 20.30 22.30 21.90 18.20 11.90  5.60  1.20

  7.36  6.94  6.47  4.74  3.98  3.41  3.27  3.38  4.44  5.32  6.19  6.47

  6.10  5.58  5.59  4.63  3.72  2.76  2.06  2.24  3.96  5.06  5.42  5.51

  66.9  72.3  89.8 120.0 133.2  93.1  81.9  66.4  88.6  87.9  81.3  77.1

 12.70 14.22 14.22 17.78 18.29 19.05 16.51 16.00 22.10 20.32 19.56 16.20

  2.23  2.27  0.91  0.94  1.50  2.78  0.94  1.44  2.27  1.93  2.41  3.00

  0.17  0.18  0.19  0.22  0.21  0.16  0.15  0.13  0.13  0.14  0.13  0.15

  7.12  6.78  7.36  8.15  9.04  6.67  6.20  5.45  5.91  5.86  5.82  6.46

 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00

 16.00 17.53 20.32 27.18 38.35 35.56 34.04 45.21 36.83 25.65 22.86 13.90

  9.42 12.27 16.41 19.40 22.95 25.01 24.50 23.12 19.22 15.54 11.40  9.05

  0.76  2.42  4.85 10.55 16.45 19.86 20.56 20.07 17.15 11.38  5.20  1.66

  4.82  4.74  5.32  5.26  4.55  4.28  3.85  3.75  3.76  4.06  4.45  4.63

[Steven Jepsen]

Bendik,

Darn, I'm stumped! Thanks for posting all of that. 

--Have you tried setting PLAPS to a small value, like 0.00001? Sometimes SWAT replaces 0 with default values.

--I wonder what would happen if you changed ELEV (mean elevation) of subbasin 5 to something other than 1073.88. Would that do anything? Maybe that only has an effect if elevation bands are not specified? Don't know.

--There are rainfall scaling factors, RFINC(mon), in *.sub files. Is it possible that those were mistakenly set to something other than 0?

It's peculiar that the precip values of all of the subbasins using that gauge have been adjusted downward. Maybe that's a clue? 

I might have a look at the source code to see if I can find anything, not that I'm an expert at locating such things. I'll post back here if I find anything of potential interest. 

Steve

[Steven Jepsen]

Bendik,

I went through source code of an old version of SWAT (rev 637) and found the familiar equations, including the 1000 and (# days precipitation per year) in the denominator. I did not find a "12" in the denominator. Not sure why this algorithm would not be taking effect in your model. For what it's worth, here are the line numbers where that code occurs in SWAT rev 637. Hopefully these line numbers are close to the version of SWAT you're using.

Lines of source code in SWAT Rev 637:

-----------------------

pmeas.f: Read precipitation data

138: subp(k) = rmeas(irgage(hru_sub(k)))

-----------------------

readsub.f: Read data from the HRU/subbasin general input file (.sub) 

!! plaps(:) = precipitation lapse ratio adjusted for number of days precipitation per year

189: read (101,*) plaps(i)

451: call readwgn

452: plaps(i) = plaps(i) / pcpdays(i)

-----------------------

readwgn.f: Reads HRU weather generator parameters 

!! pcpd(:) = average number of days of precipitation in the month read from *.wgn file 

158: read (114,5200) (pcpd(mon),mon = 1,12)

260: pcpdays(i) = pcpdays(i) + pcpd(mon)

-----------------------

clicon.f: Adjustment of precipitation for elevation bands 

150: if (pcpsim == 1) call pmeas

297-298: pdif = (elevb(ib,hru_sub(k)) - Real(elevp(irgage(hru_sub(k))))) * plaps(hru_sub(k)) / 1000.

307: pcpband(ib,k) = subp(k) + pdif

329: subp(k) = subp(k) + pcpband(ib,k) * elevb_fr(ib,hru_sub(k))

-----------------------

Steve

[marek...@gmail.com]

Hi Steve,

Thanks again. I figured out one thing at least: My elevation band fractions did not add up to exactly 1 (due to the method I used for generating them it was slightly less in most cases (95%-99%)). This was what caused the odd results when PLAPS=0 in my post above. You were right that there was a hint in the fact that they were all reduced (it was scaled down proportional to how much was missing in the sum of the elevation fractions). This fixes my first issue: now the precipitation in output.sub is correct and consistent when PLAPS=0. It was the error message example at the end of the source code file (readsub.f) that hinted me to this, so you reading the source code did help me, even if not in the intended way. 

I was hoping it would also solve the other issue, but unfortunately it did not. I am still getting results consistent with the equation not dividing by 1000 and pcpdays=12 for every subcatchment when PLAPS=/=0. I I will set up one of the tutorial projects and see how PLAPS behaves in that one, to try to identify what's wrong with my project (assuming it's not the source code that's wrong somewhere).

I checked the source code for rev 681 and it's the same as you posted. It says clearly in readsub.f that it should read the .wgn file and divide PLAPS by pcpdays. Is there any way this code (readwgn.f, 254+) could output 12 every time? 

        !! calculate precipitation-related values

        if (pcpd(mon) <= 0.) pcpd(mon) = .001

        pr_w(3,mon,i) = pcpd(mon) / mdays

        pcp_stat(mon,1,i) = pcpmm(mon) / pcpd(mon)

        if (pcp_stat(mon,3,i) < 0.2) pcp_stat(mon,3,i) = 0.2

        summm_p = summm_p + pcpmm(mon)

        pcpdays(i) = pcpdays(i) + pcpd(mon)

At the beginning of this file, it generates some dummy tables. Is the code for pcpmm supposed to be like that? I don't see it influencing pcpdays, but still curious. This generates a table filled with 0s, as far as I can tell. Maybe something happens while overwriting this. 

      use parm

      character (len=80) :: titldum

      real*8 :: xx, lattan, x1, x2, x3, tav, tmin, tmax, rain_yrs

      real*8 :: summx_t, summn_t, summm_p, sum, rnm2, r6, xlv, pcp

      real*8, dimension (12) :: rainhhmx, rain_hhsm, pcpmm, pcpd

      real*8 :: tmpsoil, sffc, rndm1, dl

      integer :: mon, mdays, j, m1, nda, xrnd

      pcpd = 0.

      rainhhmx = 0.

      pcpmm = 0.

      read (114,5000) titldum

      read (114,5100) wlat(i)

      read (114,5100) welev(i)

      read (114,5100) rain_yrs

      read (114,5200) (tmpmx(mon,i),mon = 1,12)

      read (114,5200) (tmpmn(mon,i),mon = 1,12)

      read (114,5200) (tmpstdmx(mon,i),mon = 1,12)

      read (114,5200) (tmpstdmn(mon,i),mon = 1,12)

      read (114,5201) (pcpmm(mon),mon = 1,12)

      read (114,5200) (pcp_stat(mon,2,i),mon = 1,12)  !pcpstd

      read (114,5200) (pcp_stat(mon,3,i),mon = 1,12)  !pcpskw

      read (114,5200) (pr_w(1,mon,i),mon = 1,12)

      read (114,5200) (pr_w(2,mon,i),mon = 1,12)

      read (114,5200) (pcpd(mon),mon = 1,12)

      read (114,5200) (rainhhmx(mon),mon = 1,12)

      read (114,5200) (solarav(mon,i),mon = 1,12)

      read (114,5200) (dewpt(mon,i),mon = 1,12)

      read (114,5200) (wndav(mon,i),mon = 1,12)

Thank you for all your feedback so far.

Bendik

[Steven Jepsen]

Hi Bendik,

That's good you solved half the problem.

I looked and couldn't find a problem with pcpmm(mon). It seemed decoupled from lines of code affecting lapse rate, at least based on my rusty Fortran.

Let me try to phrase the question we're trying to solve.

You are saying that the following equation seems to be what SWAT is DOING based on SWAT output and PLAPS in *.sub, right?

R_band = R_day + (EL_band - EL_gauge)*(PLAPS)/12 (eqn 1)

I can get an equivalent version of the above equation that works with your *.wgn input by multiplying numerator and denominator by 1e+4:

R_band = R_day + (EL_band - EL_gauge)*(PLAPS*10,000)/(120*1000) (eqn 2)

The denominator of eqn 2 is as it should be. Have you tried varying the number of wet days per month for months INDIVIDUALLY? For example, setting PCPD(1) = 10, PCPD(2) = 20, PCPD(3) = 40, etc? (note: PCPD(i) = # wet days in month i). I would try doing that.

If you do that and still get eqn 1, the only thing I can think to do next would be to insert print statements into the SWAT code and track what SWAT is doing step by step. This would be a sure way to find the problem. Kind of a pain, unless you're set up to compile SWAT. Then it wouldn't be that big a deal.

Steve

[Steven Jepsen]

Forgot to add-- also make sure that you change the sum of the terms (PCPD(1) + PCPD(2) + ... + PCPD(12)) from your *.wgn file. That way the "120" in my eqn 2 gets changed.

Steve

[Suresh Marahatta]

Hi Natalja and SWAT User Friends,

I have well calibrated and validated the model quite well.  However when I made a new model, the SWATOutput of Tablesout file was empty. I repeated several times still no access file is there. Please suggest how to solve this problem.

 Best Regards

--
You received this message because you are subscribed to the Google Groups "SWAT-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to swatuser+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/swatuser/83f4e158-5c4d-4a28-a3bd-20840a98042bo%40googlegroups.com.

--

Suresh Marahatta
Phone : 977-1-4325707
Cell :  98510 06425

Office:977-1-4251921

[Natalja C.]

Hello,

If it is a new model, then the output will be empty. You need to complete the setup and run the model, export the output to the tables. Then they will be populated.

Or I did not really understand your question :) elaborate, please.

All the best,

Natalja

[marek...@gmail.com]

Hi Steve,

Summary for this long post:

1: Changes to pcpd(mon) not influencing output (in other projects as well)

2: Source code has line "if (pcpd(mon) <= 0.) pcpd(mon) = .001" which might be the answer, if SWAT reads pcpd(mon) as 0.

I have no experience with fortran or compiling, so I'm afraid that would be a virtually impossible task for me. It would be a nice way to figure out what goes wrong though.

1:

I just checked PLAPS for one of the tutorial projects and it works almost properly there. This tutorial project utilizes weather generator data, and I think this is part of the problem/solution. Maybe SWAT has issues with the .wgn files when no weather generators are used. 

For the tutorial example, I have a pcp of 001.4 in pcp.pcp. 

The gauge elevation is 1790 m

The subbasin elevation (for all zones) is 2090 m

And plaps is set to 20

The weather generator file is (pcpdays=135.27): 

 .Wgn file Subbasin: 1 STATION NAME:bdrwgn 8/10/2020 12:00:00 AM ArcSWAT 2012.10_7.23

  LATITUDE =  11.60 LONGITUDE =  37.38

  ELEV [m] =1770.00

  RAIN_YRS =  10.00

 26.94 28.63 29.91 30.29 29.62 27.13 24.41 24.44 25.76 26.76 26.89 26.82

  9.09 10.96 13.31 15.43 15.63 14.97 14.53 14.36 13.74 13.73 11.80  9.65

  1.61  1.89  1.81  2.06  1.96  1.79  1.55  1.43  1.26  1.19  1.13  1.22

  2.15  2.35  2.84  2.47  1.79  1.21  1.02  1.06  1.10  1.53  2.11  2.28

   1.8   3.1  13.3  26.3  68.4 204.9 419.6 365.0 199.4  96.3  13.3   2.9

  0.54  1.08  3.11  3.77  6.86 11.56 16.20 13.48 10.23  7.19  2.33  0.77

 11.56 14.81 14.26  5.96  4.94  2.54  2.33  1.87  2.60  3.23  7.22 10.32

  0.01  0.03  0.08  0.11  0.22  0.50  0.83  0.84  0.61  0.23  0.05  0.03

  0.36  0.33  0.41  0.47  0.57  0.79  0.95  0.93  0.79  0.59  0.33  0.19

  0.69  1.13  3.50  4.50  9.69 19.75 29.44 28.63 22.56 11.69  2.69  1.00

  3.95  8.00 28.60 17.80 34.60 35.65 57.05 44.75 35.15 29.10 13.80  5.05

 20.93 22.72 23.61 24.63 23.47 21.58 19.72 19.13 21.29 22.30 21.01 20.40

  0.52  0.45  0.43  0.44  0.51  0.67  0.76  0.77  0.72  0.64  0.57  0.54

  0.53  0.59  0.71  0.80  0.74  0.80  0.66  0.64  0.61  0.65  0.58  0.53

In output.sub, I get a precipitation of 1.44, which is correct according to the equation:

1.4 + 300*20/(135.27*1000)=1.44

HOWEVER, when I change the pcpdays per month in the .wgn file, I get no corresponding change in output.sub. 

My changed .wgn file is (pcpdays=205.27):

 .Wgn file Subbasin: 1 STATION NAME:bdrwgn 8/10/2020 12:00:00 AM ArcSWAT 2012.10_7.23

  LATITUDE =  11.60 LONGITUDE =  37.38

  ELEV [m] =1770.00

  RAIN_YRS =  10.00

 26.94 28.63 29.91 30.29 29.62 27.13 24.41 24.44 25.76 26.76 26.89 26.82

  9.09 10.96 13.31 15.43 15.63 14.97 14.53 14.36 13.74 13.73 11.80  9.65

  1.61  1.89  1.81  2.06  1.96  1.79  1.55  1.43  1.26  1.19  1.13  1.22

  2.15  2.35  2.84  2.47  1.79  1.21  1.02  1.06  1.10  1.53  2.11  2.28

   1.8   3.1  13.3  26.3  68.4 204.9 419.6 365.0 199.4  96.3  13.3   2.9

  0.54  1.08  3.11  3.77  6.86 11.56 16.20 13.48 10.23  7.19  2.33  0.77

 11.56 14.81 14.26  5.96  4.94  2.54  2.33  1.87  2.60  3.23  7.22 10.32

  0.01  0.03  0.08  0.11  0.22  0.50  0.83  0.84  0.61  0.23  0.05  0.03

  0.36  0.33  0.41  0.47  0.57  0.79  0.95  0.93  0.79  0.59  0.33  0.19

 10.69 11.13 13.50 14.50 19.69 19.75 29.44 28.63 22.56 11.69 12.69 11.00

  3.95  8.00 28.60 17.80 34.60 35.65 57.05 44.75 35.15 29.10 13.80  5.05

 20.93 22.72 23.61 24.63 23.47 21.58 19.72 19.13 21.29 22.30 21.01 20.40

  0.52  0.45  0.43  0.44  0.51  0.67  0.76  0.77  0.72  0.64  0.57  0.54

  0.53  0.59  0.71  0.80  0.74  0.80  0.66  0.64  0.61  0.65  0.58  0.53

With this, I expect still get a pcp in output.sub of 1.44, but the equation yields 1.43:

1.4 + 300*20/(205.27*1000)=1.43 

This is the same as for my own project, where changing pcpdays yields no change in the output pcp. 

I think it was a coicidence that I used exactly 120 pcpdays in my example (12*10?), making your equation work. No matter what I set them to, they do not influence the output, neither in my own project or in the example project. The example project clearly does not use a fixed "12", but it does seem to have the original values for pcpdays stored elsewhere...

2:

I found one possible line of code that could be causing this, in readsub.f:

        !! calculate precipitation-related values

        if (pcpd(mon) <= 0.) pcpd(mon) = .001

        pr_w(3,mon,i) = pcpd(mon) / mdays

        pcp_stat(mon,1,i) = pcpmm(mon) / pcpd(mon)

        if (pcp_stat(mon,3,i) < 0.2) pcp_stat(mon,3,i) = 0.2

        summm_p = summm_p + pcpmm(mon)

        pcpdays(i) = pcpdays(i) + pcpd(mon)

If SWAT somehow reads pcpd(mon) as 0, it will replace it with 0.001. If every month is 0, this is equal to 0.012 (12/1000). If this is used in the equation, it is the equivalent of using a fixed "12" and omitting the 1/1000... This could be negating the 1/1000 (delta_elev*PLAPS/((12/1000)*1000)). Since changing pcpd(mon) in the .wgn file does not influence the output pcp, I am curious if SWAT reads these values elsewhere. If no wgen is used, they could be 0 in that location by default. The problem is that there is no indication in the source code regarding this. The only thing I can find when searching for "pcpd" is this, in the readfcst.f file (which should be unrelated, but maybe the parameters get overwritten here?). 

      use parm

      character (len=80) :: titldum

      real*8, dimension (12) :: pcpmm, pcpd

      integer :: mon, mdays, j, fcstregtot

      fcstregtot = 0

      i = 0

      pcpd = 0.

      pcpmm = 0.

      read (109,5000) titldum

      read (109,5100) fcstregtot

      do j = 1, fcstregtot

        read (109,5000) titldum

        read (109,5100) i    !forecast region number

          if (i <= 0 ) i = 1

        read (109,5200) (ftmpmx(mon,i),mon = 1,12)

        read (109,5200) (ftmpmn(mon,i),mon = 1,12)

        read (109,5200) (ftmpstdmx(mon,i),mon = 1,12)

        read (109,5200) (ftmpstdmn(mon,i),mon = 1,12)

        read (109,5200) (pcpmm(mon),mon = 1,12)

        read (109,5200) (fpcp_stat(mon,2,i),mon = 1,12)  !pcpstd

        read (109,5200) (fpcp_stat(mon,3,i),mon = 1,12)  !pcpskw

        read (109,5200) (fpr_w(1,mon,i),mon = 1,12)

        read (109,5200) (fpr_w(2,mon,i),mon = 1,12)

        read (109,5200) (pcpd(mon),mon = 1,12)

!! calculate missing values and additional parameters

      do mon = 1, 12

        mdays = 0

        mdays = ndays(mon+1) - ndays(mon)

        !! calculate values for fpr_w if missing or bad

        if (fpr_w(2,mon,i) <= fpr_w(1,mon,i).or.fpr_w(1,mon,i) <= 0.)   

     &                                                              then

          if (pcpd(mon) < .1) pcpd(mon) = 0.1

          fpr_w(1,mon,i) = .75 * pcpd(mon) / mdays

          fpr_w(2,mon,i) = .25 + fpr_w(1,mon,i)

        else

        !! if fpr_w values good, use calculated pcpd based on these values

        !! using first order Markov chain

        pcpd(mon) = mdays * fpr_w(1,mon,i) /                            

     &                            (1. - fpr_w(2,mon,i) + fpr_w(1,mon,i))

    

        end if

        !! calculate precipitation-related values

        if (pcpd(mon) <= 0.) pcpd(mon) = .001

        fpr_w(3,mon,i) = pcpd(mon) / mdays

        fpcp_stat(mon,1,i) = pcpmm(mon) / pcpd(mon)

        if (fpcp_stat(mon,3,i) < 0.2) fpcp_stat(mon,3,i) = 0.2

      end do

The developers are not getting back to me, unfortunately. 

Bendik

[Steven Jepsen]

Hi Bendik,

Quick reply on point #1. I'll have to look at the code more for point #2, which will take me longer.

1:

"With this, I expect still get a pcp in output.sub of 1.44, but the equation yields 1.43:

1.4 + 300*20/(205.27*1000)=1.43"

Why do you expect a value of 1.44? I get 1.43 with your increase in pcpdays. Did SWAT write 1.43 to your output.sub file? Wouldn't it make sense to use a larger PLAPS value to avoid such small differences (i.e., 1.43 vs 1.44)?

2:

I could be totally wrong about this, but isn't that readfcst.f code only used if you are doing runs in "forecast" mode? Are you working in forecast mode? I've never done that. I recall seeing somewhere a flag for activating forecast mode, maybe in file.cio? Just a guess--don't want to lead you down wrong path.

Possible pragmatic solution: Couldn't you apply a scaling factor your input PLAPS values to give you the lapsed precipitation that you expect? This would in theory be doable by comparing your pcp.pcp data, your output precipitation in output.sub, your elevation bands, and your PLAPS. You would apply a scaling factor to PLAPS which would account for whatever constant value of pcpdays that SWAT was in fact using. This would in a sense be back-solving for pcpdays. Do you think that would that work?

Lastly, let me know if there is something I can do to help. I know a little Fortran and an associate that has successfully compiled SWAT using free gfortran. I could work on this, but it would be on my free time and thus probably a slow process, unfortunately. Given that it involves the forcing of precipitation, it is clearly an important issue to get to the bottom of. 

Steve

[Steven Jepsen]

Bendik,

One more thing came to mind... If you have a forecast file in your project folder (see Ch 13 about *.cst in I/O document), SWAT may be reading weather parameters from that file instead of your *.wgn file. Which of those would depend on the time period of your forecast simulation. Deleting that file may tell SWAT to read parameters from your *.wgn file. A similar problem came up for me when I couldn't tell from where SWAT was reading its nitrate concentration in rainfall. After deleting the ATMO.DAT file, which I must have mistakenly created in ArcSWAT, SWAT reverted to reading from the .bsn file as I was expecting.

Steve

[marek...@gmail.com]

Hi Steve,

1: Sorry there was a typo, I didn't mean to write expect. I meant to write "I still get". I expect to get 1.43 (as the equation shows), but I get 1.44. I ran it again with significantly lower values for pcpd(mon), and no change. I used a total of 21 days with the same data as above. The result should be 1.69, but remains 1.44. I think this is pretty decent proof that SWAT does not read the pcpd(mon) value from the .wgn file. I even DELETED all the numbers in the 14th row in the .wgn file and left it blank, and SWAT still ran fine and gave me the same outputs (the same outputs as if I had not changed the .wgn file at all). They must be stored somewhere else as well, but it's very weird because I have done a search (after setting up indexing on my PC to search content in all SWAT-related file extensions) in all the files in the SWAT-CUP folder where I run SWAT, and no files contain any of the pcpd(mon) values from .wgn, nor the sum of the pcpd(mon) values (other than some random files that happen to have the same number for a concentration or value on a specific day, like discharge).

I also tried changing ALL the .wgn files, not just one specific one and looking at the outputs for one subcatchment, no difference. 

2: I don't use any forecasting, so I also don't think the file should have any influence. It was just that the fcst file was the only additional file I could find that had anything to do with the "pcpd" parameter. 

I think your suggestion for a workaround sounds feasible. I feel like my colleague and I are on the trail of something now with some weather generator stuff and the swat project Microsoft access database, but if that doesn't work out then we'll go for a workaround like this. I will also make another post here to see if anyone knows where SWAT actually reads pcpd(mon), since that is what it seems to boil down to in the end. I don't expect anyone except you and me to read this long exchange.

Thank you for the offer to help, I appreciate it. I don't think we need to go down the "changing the source code" rode for now, unless the other workarounds don't work.

[Steven Jepsen]

Bendik,

1. You've shown well that SWAT is getting its pcpd(mon) values from somewhere other than the *.wgn files. The only other source I can think of where SWAT would get pcpd(mon) is from the forecast input files, so myself I would double-check that a mistake was not made that causes SWAT to read pcpd(mon) from the forecast files rather than the *.wgn files.

2. OK, you don't do forecasting--thanks. Could a mistake have been made in inputs somewhere causing SWAT to read from the forecast file instead of the *.wgn file? If that file name is listed in your file.cio file, then SWAT will read from that file even though you don't intend for that to happen.That is how the "optional" files seem to work with SWAT.

Sounds good. Let me know if I can help. I'm interested in finding out what's causing this problem.

Steve

 

...