Check for NaN/Infinity values in atmospheric forcing added

source/model/lake.f90

@@ -269,7 +269,7 @@ use TRIBUTARIES, only : TRIBTEMP
 
 use BATHYMSUBR, only : BATHYMETRY
 
-use NUMERICS, only : ACCUMM
+use NUMERICS, only : ACCUMM, VALUE_IS_FINITE
 
 use BL_MOD_LAKE
 
@@ -471,6 +471,7 @@ data nstep_meas /0/
  
 SAVE
   
+
 allocate(a(1:nveclen),b(1:nveclen),c(1:nveclen), &
 &        d(1:nveclen),Temp(1:nveclen))
 
@@ -584,7 +585,6 @@ else
   netrad = netrad1
 endif
 
-!print*, tempair, humair, pressure, uwind, vwind, longwave, shortwave
 
 if (longwave1 == missing_value .and. firstcall) then
   write(*,*) 'Atmospheric radiation is missing in input file and will be calculated &
@@ -601,37 +601,36 @@ if (uwind == 0) uwind = minwind
 if (vwind == 0) vwind = minwind
 wind = sqrt(uwind**2+vwind**2)
 !wind10 = wind*log(10./roughness0)/log(zref/roughness0)
-precip = 0.
 !Control of the input atmospheric forcing
 flag = .false.
-if (tempair>60.or.tempair<-90) then
+if (tempair > 60. .or. tempair < -90. .or. (.not. VALUE_IS_FINITE(tempair) )) then
   print*, 'The air temperature ', tempair, 'is illegal: STOP'
   flag = .true.
-elseif (abs(humair)>1.) then
+elseif (abs(humair) > 1. .or. (.not. VALUE_IS_FINITE(humair) )) then
   print*, 'The air humidity ', humair, 'is illegal: STOP'
   flag = .true.
-elseif (pressure > 110000 .or. pressure < 80000.) then
+elseif (pressure > 110000 .or. pressure < 80000. .or. (.not. VALUE_IS_FINITE(pressure) )) then
   print*, 'The air pressure ', pressure, 'is illegal: STOP'
   flag = .true.
-elseif (abs(uwind)>200.) then
+elseif (abs(uwind) > 200. .or. (.not. VALUE_IS_FINITE(uwind) )) then
   print*, 'The x-component of wind ', uwind, 'is illegal: STOP'
   flag = .true.
-elseif (abs(vwind)>200.) then
+elseif (abs(vwind) > 200. .or. (.not. VALUE_IS_FINITE(vwind) )) then
   print*, 'The y-component of wind ', vwind, 'is illegal: STOP'
   flag = .true.
-elseif (abs(longwave)>1000.) then
+elseif (abs(longwave) > 1000. .or. (.not. VALUE_IS_FINITE(longwave) )) then
   print*, 'The longwave radiation ',longwave,'is illegal: STOP'
   flag = .true.
-elseif (abs(shortwave)>1400.) then
+elseif (abs(shortwave) > 1400. .or. (.not. VALUE_IS_FINITE(shortwave) )) then
   print*, 'The shortwave radiation ', shortwave, 'is illegal: STOP'
   flag = .true.
-elseif (abs(precip)>1.) then
+elseif (abs(precip) > 1. .or. (.not. VALUE_IS_FINITE(precip) )) then
   print*, 'The atmospheric precipitation ',precip, &
   & 'is illegal: STOP'
   flag = .true.
 endif
 if (flag) then
-  write(*,*) 'Temperature = ', tempair, &
+  write(*,*) 'LAKE: atmospheric forcing is incorrect: Temperature = ', tempair, &
   & 'Humidity = ', humair, &
   & 'Pressure = ', pressure, &
   & 'Uwind = ', uwind, &
@@ -2447,7 +2446,7 @@ endif
 
 
 !print*,h1,l1,hs1,ls1,Tw2(1),Sal2(1),Sal2(M+1)
-!if (ix == 10 .and. iy == 10) print*, 'Lake', &
+!if (ix == 1 .and. iy == 1) print*, 'Lake', &
 !& 'T1=',Tw1(1),'T2=',Tw1(2),'h1=',h1,'S=',shortwave, &
 !& 'hh=',hflux,'LE=',Elatent, &
 !& 'eflux=',eflux, 'Longwave=', longwave, &
@@ -2912,6 +2911,7 @@ deallocate(radflux)
 deallocate(a,b,c,d)
 deallocate(Temp)
 
+
 if (firstcall) firstcall = .false.
 
 !FORMATS!

source/model/numerics_mod.f90

@@ -2,8 +2,18 @@ MODULE NUMERICS
 
 use LAKE_DATATYPES, only : ireals, iintegers
 use NUMERIC_PARAMS, only : nveclen
+use IEEE_ARITHMETIC, only : IEEE_IS_NAN, IEEE_IS_FINITE
 
 contains
+
+!> Function VALUE_IS_FINITE checks, whether the value is not infinite nor NaN
+FUNCTION VALUE_IS_FINITE(val)
+implicit none
+real(kind=ireals), intent(in) :: val !>Real value to be checked
+logical :: VALUE_IS_FINITE
+VALUE_IS_FINITE = (.not.(IEEE_IS_NAN(val)) .and. IEEE_IS_FINITE(val))
+END FUNCTION VALUE_IS_FINITE
+
 SUBROUTINE MATRIXSUM(nn,a,b,c,k)
  implicit none