Forcing=2 with IceModel=1 and UseFilteredWind=False and SnowModel=0 reads the wrong number of entries from forcing file

[akesandgren]

With Forcing=2, IceModel=1, UseFilteredWind=False and SnowModel=0 forcing_update sets nval_offset=1 and then calls forcing_read with

call self%read (state%datum, A_s, A_e, A_cur, 5 + nval_offset, state%first_timestep)

This makes forcing_read to read date + 6 values from the forcing file. A properly created forcing file with date + 5 values (Forcing=2) then gets read in such a way that it skips every second line.

The

read (20, *, end=9) tb_start, (A_s(i), i=1, nval)

will read the date, the 5 values and then the first value of the next line. The next read (tb_end) will then read the 3rd line.

This is with GCC/10.3.

The forcing_read needs the cfg%ice_model, cfg%use_filtered_wind values to decide how much to read, or the code that sets nval_offset is wrong:

if (cfg%use_filtered_wind .and. cfg%ice_model == 0) then
      nval_offset = 1
else if (cfg%ice_model == 1 .and. cfg%use_filtered_wind) then
     nval_offset = 2
else if (cfg%ice_model == 1) then
     nval_offset = 1
else
     nval_offset = 0
end if

Should the

else if (cfg%ice_model == 1 .and. cfg%use_filtered_wind) then
  nval_offset = 2
else if (cfg%ice_model == 1) then
  nval_offset = 1

check cfg%snow_model instead of cfg%ice_model ??

[f-baerenbold]

Hello,

Yes indeed, it should use cfg%snow_model instead of cfg%ice_model.

I don't understand why you think nval_offset should be 2 if cfg%snow_model ==1 but no filtered wind. Can you explain? In my opinion, the filtered wind and the snow_model each add one additional column. So if only snow_model is on but not the filtered wind, nval_offset should be 1.

[akesandgren]

Ah, that's just a typo when cut&past-ing the code. Sorry, fixed the original text.

[f-baerenbold]

Ok, now worries.

But I still don't see why one has to change this sequence except replacing cfg%ice_model with cfg%snow_model. Can you explain?

[akesandgren]

I'm not suggesting to change the sequence, I just cut out the part of it where ice should be changed to snow.

[f-baerenbold]

Ok, I see. I'll adjust it on the master.

[akesandgren]

Actually I think there is a problem with your change (because of the original code being even more wrong than I thought): Looking at what the code under forcing=2, 3, or 5 does (example from forcing=2)

            if (cfg%use_filtered_wind) state%Wf = A_cur(6) ! AG 2014
            if (cfg%snow_model == 1 .and. cfg%use_filtered_wind) then
               state%precip = max(A_cur(7),0.0_RK)
            else if (cfg%snow_model == 1) then
               state%precip = max(A_cur(6),0.0_RK)
            end if

the nval_offset should look like:

      if (cfg%use_filtered_wind .and. cfg%snow_model == 1 ) then
         nval_offset = 2
      else if (cfg%use_filtered_wind) then
         nval_offset = 1
      else if (cfg%snow_model == 1) then
         nval_offset = 1
      else
         nval_offset = 0
      end if

or perhaps

nval_offset = 0
if (cfg%use_filtered_wind) then
  nval_offset = nval_offset + 1
end if
if (cfg%snow_model == 1) then
  nval_offset = nval_offset + 1
endif

[akesandgren]

@f-baerenbold Just making sure you see my last comment

[f-baerenbold]

I don't really see what is wrong here. As far as I see, the code currently sets nval_offset = 2 if both filtered wind and snow model are "on", nval_offset = 1 if only one is "on" and 0 if none is "on".

The currently implemented solution is not the most elegant, but I don't see how the result differs from your suggestion. Can you elaborate in which case it would be wrong?

[akesandgren]

You're right, I misread the change. The change you made altered the

if (cfg%use_filtered_wind .and. cfg%ice_model == 0) then

to use "cfg%snow_model==0" instead and that confused me a bit.

The last suggestion I made above is however much simpler to read and does the same thing.

Sorry about the confusion.

[f-baerenbold]

Ok good. I agree that your last solution looks better. I'll implement it.

If you like, you can directly contribute yourself by creating a pull request next time you spot something.

[akesandgren]

I know, I just didn't have time to pull yet another repo down to make PR's :-) I have way too many already...