diff --git a/source/model/lake.f90 b/source/model/lake.f90 index 9a79c2ce5fea5af806987344244fea1135efa228..ba60c978da654990a2fb5deeb6316d735e6e4c1e 100644 --- a/source/model/lake.f90 +++ b/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! diff --git a/source/model/numerics_mod.f90 b/source/model/numerics_mod.f90 index 8f9f0179d6fcccacb5fe7389432b234d7ed67f02..d1ce565191d0bd7ffe3a07aca77443acba62e2ff 100644 --- a/source/model/numerics_mod.f90 +++ b/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