source: trunk/LMDZ.PLUTO/libf/phypluto/check_fields.F90 @ 3380

module check_fields_mod

  real,save :: default_temp_min=5.  ! minimum reasonable temperature (K)
  real,save :: default_temp_max=5000. ! maximum reasonable temperature (K)
!$OMP THREADPRIVATE(default_temp_min,default_temp_max)

  real,save :: default_wind_max=5000. ! maximum reasonable wind magnitude (m/s)
!$OMP THREADPRIVATE(default_wind_max)

  real,save :: default_ps_min=80.  ! minimum reasonable surface pressure (Pa)
  real,save :: default_ps_max=3000000. ! maximum reasonable surface pressure (Pa)
!$OMP THREADPRIVATE(default_ps_min,default_ps_max)

contains

subroutine check_physics_fields(message,temp,u,v,pplev,q)
use dimphy, only: klon, klev
use tracer_h, only: nqtot
use ioipsl_getin_p_mod, only: getin_p

implicit none
character(len=*),intent(in):: message ! text message for outputs
real,intent(in) :: temp(klon,klev)
real,intent(in) :: u(klon,klev) ! zonal wind (m/s)
real,intent(in) :: v(klon,klev) ! meridional wind (m/s)
real,intent(in) :: pplev(klon,klev+1) ! pressure at level interfaces (Pa)
real,intent(in) :: q(klon,klev,nqtot) ! tracer mixing ratio (.../kg_air)

character(len=50) :: name="check_physics_fields"
logical :: ok_t,ok_w,ok_ps,ok_q
logical, save :: firstcall = .true.
!$OMP THREADPRIVATE(firstcall)

! 0. Get user defaults
if (firstcall) then
   call getin_p("check_temp_min", default_temp_min)
   call getin_p("check_temp_max", default_temp_max)
   call getin_p("check_wind_max", default_wind_max)
   call getin_p("check_ps_min", default_ps_min)
   call getin_p("check_ps_max", default_ps_max)
   firstcall = .false.
endif

! 1. Initialisations
ok_q=.true.

! 2. Check temperature, winds and surface pressure
call check_temperature(message,temp,ok_t)
call check_winds(message,u,v,ok_w)
call check_ps(message,pplev(:,1),ok_ps)
if (nqtot>=1) then
  call check_tracers(message,q,ok_q)
endif

if ((.not.ok_t).or.(.not.ok_w).or.(.not.ok_ps).or.(.not.ok_q)) then
  ! Something is wrong, might as well stop here
  if (.not.ok_t) write(*,*) trim(name)//":"//trim(message)//"Bad temperature values"
  if (.not.ok_w) write(*,*) trim(name)//":"//trim(message)//"Bad wind values "
  if (.not.ok_ps) write(*,*) trim(name)//":"//trim(message)//"Bad surface pressure values"
  if (.not.ok_q) write(*,*) trim(name)//":"//trim(message)//"Bad tracers values"
  call abort_physic(trim(name),trim(message)//" Invalid field values",1)
endif

end subroutine check_physics_fields


subroutine check_temperature(message,temp,ok,temp_min,temp_max)
use dimphy, only: klon, klev
implicit none
character(len=*),intent(in):: message ! text message for outputs
real,intent(in) :: temp(klon,klev)
logical,intent(out) :: ok ! returns .true. if everything OK, .false. otherwise
real,intent(in),optional :: temp_min ! user provided acceptable minimum
real,intent(in),optional :: temp_max ! user provided acceptable maximum

character(len=50) :: name="check_temperature"
real :: tmin,tmax
integer :: i,k

! 0. Check optional inputs
if (present(temp_min)) then
  tmin=temp_min
else
  tmin=default_temp_min
endif

if (present(temp_max)) then
  tmax=temp_max
else
  tmax=default_temp_max
endif

! 1. initializations
ok=.true.

! 2. do the checking
do i=1,klon
  do k=1,klev
    ! look for NaN
    if (temp(i,k).ne.temp(i,k)) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" temp(i,k)=",temp(i,k),&
                " for i=",i," k=",k
    endif
    ! check for temperatures too low
    if (temp(i,k).lt.tmin) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" temp(i,k)=",temp(i,k),&
      " for i=",i," k=",k," <",tmin
    endif
    ! check for temperatures too high
    if (temp(i,k).gt.tmax) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" temp(i,k)=",temp(i,k),&
      " for i=",i," k=",k," >",tmax
    endif
  enddo
enddo

end subroutine check_temperature

subroutine check_winds(message,u,v,ok,wind_max)
use dimphy, only: klon, klev
implicit none
character(len=*),intent(in):: message ! text message for outputs
real,intent(in) :: u(klon,klev) ! zonal wind (m/s)
real,intent(in) :: v(klon,klev) ! meridional wind (m/s)
logical,intent(out) :: ok ! returns .true. if everything OK, .false. otherwise
real,intent(in),optional :: wind_max ! user provided acceptable maximum magnitude

character(len=50) :: name="check_winds"
real :: wmax
integer :: i,k

! 0. Check optional inputs

if (present(wind_max)) then
  wmax=wind_max
else
  wmax=default_wind_max
endif

! 1. initializations
ok=.true.

! 2. do the checking
do i=1,klon
  do k=1,klev
    ! look for NaN
    if (u(i,k).ne.u(i,k)) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" u(i,k)=",u(i,k),&
      " for i=",i," k=",k
    endif
    if (v(i,k).ne.v(i,k)) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" v(i,k)=",v(i,k),&
      " for i=",i," k=",k
    endif
    ! check for magnitudes too high
    if (abs(u(i,k)).gt.wmax) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" u(i,k)=",u(i,k),&
      " for i=",i," k=",k," >",wmax
    endif
    if (abs(v(i,k)).gt.wmax) then
      ok=.false.
      write(*,*)trim(message)//" "//trim(name)//" v(i,k)=",v(i,k),&
      " for i=",i," k=",k," >",wmax
    endif
  enddo
enddo

end subroutine check_winds

subroutine check_ps(message,ps,ok,ps_min,ps_max)
use dimphy, only: klon
implicit none
character(len=*),intent(in):: message ! text message for outputs
real,intent(in) :: ps(klon)
logical,intent(out) :: ok ! returns .true. if everything OK, .false. otherwise
real,intent(in),optional :: ps_min ! user provided acceptable minimum
real,intent(in),optional :: ps_max ! user provided acceptable maximum

character(len=50) :: name="check_ps"
real :: pmin,pmax
integer :: i

! 0. Check optional inputs
if (present(ps_min)) then
  pmin=ps_min
else
  pmin=default_ps_min
endif

if (present(ps_max)) then
  pmax=ps_max
else
  pmax=default_ps_max
endif

! 1. initializations
ok=.true.

! 2. do the checking
do i=1,klon
  ! look for NaN
  if (ps(i).ne.ps(i)) then
    ok=.false.
    write(*,*)trim(message)//" "//trim(name)//" ps(i)=",ps(i)," for i=",i
  endif
  ! check for pressures too low
  if (ps(i).lt.pmin) then
    ok=.false.
    write(*,*)trim(message)//" "//trim(name)//" ps(i)=",ps(i)," for i=",i,&
    " <",pmin
  endif
  ! check for pressures too high
  if (ps(i).gt.pmax) then
    ok=.false.
    write(*,*)trim(message)//" "//trim(name)//" ps(i)=",ps(i)," for i=",i,&
    " >",pmax
  endif
enddo

end subroutine check_ps

subroutine check_tracers(message,q,ok)
use dimphy, only: klon, klev
use tracer_h, only: nqtot,noms
implicit none
character(len=*),intent(in):: message ! text message for outputs
real,intent(in) :: q(klon,klev,nqtot) ! tracer mixing ratio
logical,intent(inout) :: ok ! set to .false. if something is wrong

character(len=50) :: name="check_tracers"
integer :: i,k,iq
integer :: nb_bad_neg ! number of problematic points (negative values)
integer :: nb_bad_nan ! number of problematic points (NaNs)

! 1. initializations
nb_bad_nan=0

! 2. do the checking
do iq=1,nqtot
 nb_bad_neg=0 ! initialize
 do i=1,klon
  do k=1,klev
    ! look for NaN
    if (q(i,k,iq).ne.q(i,k,iq)) then
      ok=.false.
      nb_bad_nan=nb_bad_nan+1
    endif
    ! look for negative values
    if (q(i,k,iq).lt.0) then
!      ok=.false.
      nb_bad_neg=nb_bad_neg+1
    endif
  enddo
 enddo
 if (nb_bad_neg>0) then
   write(*,*)trim(message)//" "//trim(name)//" tracer "//trim(noms(iq))//&
   " contains ",nb_bad_neg," negative values!"
 endif
enddo !of do iq=1,nqtot

end subroutine check_tracers

end module check_fields_mod