Support sub-daily forcing time steps and units - #44
Conversation
martynpclark
left a comment
martynpclark
left a comment
There was a problem hiding this comment.
This is looking very good. I have a number of comments to improve code readability and functionality.
| ! Syntax: length per time | ||
| ipos = index(units_lc, " per ") | ||
|
|
||
| if (ipos <= 1 .or. ipos + 4 > len_trim(units_lc)) then | ||
| ierr = 20 | ||
| message = trim(message)// & | ||
| "unsupported flux units '"//trim(cunits)// & | ||
| "': expected length/time or length per time" | ||
| return | ||
| end if |
There was a problem hiding this comment.
This is an unusual time unit. Fine it is supported.
| end if | ||
| end do | ||
|
|
||
| ! Current supported syntax: length/time or length per time |
There was a problem hiding this comment.
I think this would be better as a case statement. You can identify unique sub-strings and process accordingly. The most commonly used time units are likely L/T and L T-1. The best approach may be to first search for sub-strings, i.e., "/" and "-1" (and if you want "per", and then search for delimiters, i.e., "/", " ", and " per " for the three supported options. Then you can extract sub-strings for length and time. For generality it is probable best to have the sub-string extraction based on the length of the delimiter and avoid magic numbers in the code.
| ierr = 20 | ||
| message = trim(message)// & | ||
| "unsupported flux units '"//trim(cunits)// & | ||
| "': expected length/time" |
There was a problem hiding this comment.
We expect other flux units also.
| if (ipos > 0) then | ||
|
|
||
| ! Syntax: length/time | ||
| if (ipos <= 1 .or. ipos >= len_trim(units_lc)) then |
There was a problem hiding this comment.
I do not understand why this erroe check is here since it is within the ipos>0 if statement. Also, can ipos ever be greater than len_trim(units_lc)?
| "centimetre", "centimetres") | ||
| length_to_mm = 10._wp | ||
|
|
||
| case ("dm", "decimeter", "decimeters", & |
There was a problem hiding this comment.
Are these units common? If not, then suggest delete this case to ensure we exit with an error.
| case (iPRECIP); gForce_3d(:,:,1:numtim)%ppt = gTemp(:,:,:) * amult_ppt | ||
| case (iTEMP) ; gForce_3d(:,:,1:numtim)%temp = gTemp(:,:,:) | ||
| case (iPET) ; gForce_3d(:,:,1:numtim)%pet = gTemp(:,:,:) | ||
| case (iQOBS) ; aValid( :,:,1:numtim)%obsq = gTemp(:,:,:) ! TODO: check dimensions (works for nx=1, ny=1) |
There was a problem hiding this comment.
We need to retain this TODO comment somewhere.
|
|
||
| subroutine get_forcing_varids(ncid, info, ierr, message) | ||
| use multiforce, only: amult_ppt, amult_pet, amult_q |
There was a problem hiding this comment.
We need to move away from using shared variables in modules. This would be better stored in the info data structure.
| ! -------------------------------------------------------------------------------------- | ||
|
|
||
| subroutine get_forcing_varids(ncid, info, ierr, message) |
There was a problem hiding this comment.
I think we need to change the name of the subroutine because it now does a lot more -- suggest get_forcing_metadata
| ! Verify that the forcing time axis is increasing and regularly spaced. | ||
| do i = 3, size(jdate) | ||
| dt_current = (time_steps(i) - time_steps(i-1)) * scale_to_days | ||
|
|
||
| if (dt_current <= 0._wp) then | ||
| ierr = 1 | ||
| message = trim(message)//"forcing time axis must be strictly increasing" | ||
| return | ||
| endif | ||
|
|
||
| if (abs(dt_current - deltim_days) > tolerance) then | ||
| ierr = 1 | ||
| message = trim(message)//"forcing time steps are not equally spaced" | ||
| return | ||
| endif | ||
| enddo |
There was a problem hiding this comment.
This will increase compute time. Can we wrap this in a flag, e.g.,
logical(lgt), parameter :: do_timeCheck=.false.
and then
if(timeCheck)then
newcode
endif
It is possible that the extra time is minimal, and, if so, the flag can be =.true.
|
|
||
| ! Compute the number of routing bins from the maximum routing horizon and the forcing timestep, both expressed in days. | ||
| NTDH = CEILING(TDH_MAX / info%time%deltim_days, KIND=I4B) | ||
|
|
||
| ! Allocate the runoff fractions assigned to future routing bins. | ||
| ALLOCATE(DPARAM%FRAC_FUTURE(NTDH), STAT=ISTAT) | ||
| IF (ISTAT /= 0) THEN | ||
| ERR = 20 | ||
| MESSAGE = TRIM(MESSAGE)//"cannot allocate DPARAM%FRAC_FUTURE" | ||
| RETURN | ||
| ENDIF | ||
|
|
||
| ! Allocate the routing queue using the same number of bins. | ||
| ALLOCATE(FUTURE(NTDH), STAT=ISTAT) | ||
| IF (ISTAT /= 0) THEN | ||
| DEALLOCATE(DPARAM%FRAC_FUTURE) | ||
| ERR = 20 | ||
| MESSAGE = TRIM(MESSAGE)//"cannot allocate FUTURE" | ||
| RETURN | ||
| ENDIF | ||
|
|
||
| ! Initialise the routing queue before the first model evaluation. | ||
| FUTURE = 0._WP | ||
|
|
There was a problem hiding this comment.
Can we put these allocate statements in q_timedelay? This will require passing the info structures through par_derive. We could also elevate the call to q_timedelay here.
There was a problem hiding this comment.
Making this comment because the intent is to keep the setup routines short and readable.
2 participants
Summary
This PR adds support for sub-daily forcing time steps and automatic conversion of input fluxes units.
Main changes
mm/dayrepresentationSupported units
Length:
Time:
Both
length/timeandlength per timesyntaxes are supported (case-insensitive, with common aliases such ashour,hours,hr, etc.).Testing