Index: trunk/LMDZ.GENERIC/libf/phystd/check_fields.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/check_fields.F90	(revision 2663)
+++ trunk/LMDZ.GENERIC/libf/phystd/check_fields.F90	(revision 2663)
@@ -0,0 +1,269 @@
+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_min)
+   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
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 2662)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 2663)
@@ -13,5 +13,6 @@
                   flxw,                    &
                   pdu,pdv,pdt,pdq,pdpsrf)
- 
+
+      use ioipsl_getin_p_mod, only: getin_p 
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
       use watercommon_h, only : RLVTT, Psat_water,epsi,su_watercycle, RV, T_h2o_ice_liq
@@ -61,4 +62,5 @@
                               waterrain, global1d, szangle
       use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
+      use check_fields_mod, only: check_physics_fields
       use conc_mod, only: rnew, cpnew, ini_conc_mod
       use phys_state_var_mod
@@ -439,4 +441,9 @@
 #endif
 
+      ! flags to trigger extra sanity checks
+      logical, save ::  check_physics_inputs=.false.
+      logical, save ::  check_physics_outputs=.false.
+!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)     
+
       ! Misc
       character*2 :: str2
@@ -451,4 +458,7 @@
 ! --------------------------------
       if (firstcall) then
+        call getin_p("check_physics_inputs", check_physics_inputs)
+        call getin_p("check_physics_outputs", check_physics_outputs)
+
         ! Allocate saved arrays (except for 1D model, where this has already
         ! been done)
@@ -714,4 +724,9 @@
          write(*,*) "physiq: end of firstcall"
       endif ! end of 'firstcall'
+
+      if (check_physics_inputs) then
+         !check the validity of input fields coming from the dynamics
+         call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
+      endif
 
 ! ------------------------------------------------------
@@ -2407,4 +2422,9 @@
 #endif
       
+      if (check_physics_outputs) then
+         ! Check the validity of updated fields at the end of the physics step
+         call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
+      endif
+
       icount=icount+1