Commit 362cabf1 authored by kbk's avatar kbk
Browse files

optional sanity checks on velocities, temperature and salinity fields

parent 3fd03b25
!$Id: m3d.F90,v 1.32 2006-08-25 09:00:19 kbk Exp $
!$Id: m3d.F90,v 1.33 2006-12-15 09:57:50 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -57,6 +57,8 @@
logical :: bdy3d=.false.
integer :: bdyfmt_3d,bdyramp_3d
character(len=PATH_MAX) :: bdyfile_3d
integer :: vel_check=0
REALTYPE :: min_vel=-4.,max_vel=4.
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
......@@ -116,7 +118,8 @@
bdy3d,bdyfmt_3d,bdyramp_3d,bdyfile_3d, &
vel_hor_adv,vel_ver_adv,vel_adv_split, &
calc_temp,calc_salt, &
avmback,avhback,ip_method
avmback,avhback,ip_method, &
vel_check,min_vel,max_vel
!
!EOP
!-------------------------------------------------------------------------
......@@ -211,6 +214,17 @@
case default
end select
LEVEL2 'vel_check=',vel_check
if (vel_check .ne. 0) then
LEVEL3 'doing sanity checks on velocities every ',vel_check,' timesteps'
if (vel_check .gt. 0) then
LEVEL3 'out-of-bound values result in termination of program'
end if
if (vel_check .lt. 0) then
LEVEL3 'out-of-bound values result in warnings only'
end if
end if
dt = M*timestep
#ifdef CONSTANT_VISCOSITY
......@@ -383,12 +397,12 @@
huo=hun
hvo=hvn
if (ufirst) then
call uu_momentum_3d(bdy3d)
call vv_momentum_3d(bdy3d)
call uu_momentum_3d(n,bdy3d)
call vv_momentum_3d(n,bdy3d)
ufirst=.false.
else
call vv_momentum_3d(bdy3d)
call uu_momentum_3d(bdy3d)
call vv_momentum_3d(n,bdy3d)
call uu_momentum_3d(n,bdy3d)
ufirst=.true.
end if
#ifndef MUDFLAT
......
!$Id: salinity.F90,v 1.22 2006-08-25 09:00:19 kbk Exp $
!$Id: salinity.F90,v 1.23 2006-12-15 09:57:50 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -23,7 +23,7 @@
#endif
use domain, only: iimin,jjmin,iimax,jjmax,kmax
use domain, only: H,az
use variables_3d, only: S,hn,adv_schemes
use variables_3d, only: S,hn,adv_schemes,kmin
use halo_zones, only: update_3d_halo,wait_halo,D_TAG
IMPLICIT NONE
!
......@@ -40,6 +40,8 @@
integer :: salt_hor_adv=1,salt_ver_adv=1
integer :: salt_adv_split=0
REALTYPE :: salt_AH=-_ONE_
integer :: salt_check=0
REALTYPE :: min_salt=0.,max_salt=40.
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
......@@ -100,10 +102,11 @@
integer, parameter :: nmax=100
REALTYPE :: zlev(nmax),prof(nmax)
integer :: salt_field_no=1
NAMELIST /salt/ &
salt_method,salt_const,salt_file, &
salt_format,salt_name,salt_field_no, &
salt_hor_adv,salt_ver_adv,salt_adv_split,salt_AH
NAMELIST /salt/ &
salt_method,salt_const,salt_file, &
salt_format,salt_name,salt_field_no, &
salt_hor_adv,salt_ver_adv,salt_adv_split,salt_AH, &
salt_check,min_salt,max_salt
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -301,6 +304,18 @@ salt_field_no=1
call wait_halo(D_TAG)
call mirror_bdy_3d(S,D_TAG)
LEVEL3 'salt_check=',salt_check
if (salt_check .ne. 0) then
LEVEL4 'doing sanity check on salinity'
if (salt_check .gt. 0) then
LEVEL4 'out-of-bound values result in termination of program'
end if
if (salt_check .lt. 0) then
LEVEL4 'out-of-bound values result in warnings only'
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving init_salinity()'
write(debug,*)
......@@ -368,6 +383,7 @@ salt_field_no=1
REALTYPE :: delxu(I2DFIELD),delxv(I2DFIELD)
REALTYPE :: delyu(I2DFIELD),delyv(I2DFIELD)
REALTYPE :: area_inv(I2DFIELD)
integer :: status
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -489,6 +505,22 @@ salt_field_no=1
end do
#endif
if (salt_check .ne. 0 .and. mod(n,abs(salt_check)) .eq. 0) then
call check_3d_fields(imin,jmin,imax,jmax,az, &
iimin,jjmin,iimax,jjmax,kmax, &
kmin,S,min_salt,max_salt,status)
if (status .gt. 0) then
if (salt_check .gt. 0) then
call getm_error("do_salinity()", &
"out-of-bound values encountered")
end if
if (salt_check .lt. 0) then
LEVEL1 'do_salinity(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving do_salinity()'
write(debug,*)
......
!$Id: temperature.F90,v 1.18 2006-08-25 09:00:19 kbk Exp $
!$Id: temperature.F90,v 1.19 2006-12-15 09:57:50 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -20,7 +20,7 @@
use exceptions
use domain, only: imin,jmin,imax,jmax,H,az
use domain, only: iimin,jjmin,iimax,jjmax,kmax
use variables_3d, only: T,hn,adv_schemes
use variables_3d, only: T,hn,adv_schemes,kmin
use halo_zones, only: update_3d_halo,wait_halo,D_TAG
IMPLICIT NONE
!
......@@ -37,6 +37,8 @@
integer :: temp_hor_adv=1,temp_ver_adv=1
integer :: temp_adv_split=0
REALTYPE :: temp_AH=-1.
integer :: temp_check=0
REALTYPE :: min_temp=-2.,max_temp=35.
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
......@@ -85,9 +87,10 @@
REALTYPE :: zlev(nmax),prof(nmax)
integer :: temp_field_no=1
NAMELIST /temp/ &
temp_method,temp_const,temp_file, &
temp_format,temp_name,temp_field_no, &
temp_hor_adv,temp_ver_adv,temp_adv_split,temp_AH
temp_method,temp_const,temp_file, &
temp_format,temp_name,temp_field_no, &
temp_hor_adv,temp_ver_adv,temp_adv_split,temp_AH, &
temp_check,min_temp,max_temp
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -157,7 +160,7 @@ temp_field_no=1
call getm_error("init_3d()", &
"temp_adv_split=0: temp_ver_adv not valid (2-6)")
end select
LEVEL2 "1D split --> full u, full v, full w"
LEVEL3 "1D split --> full u, full v, full w"
case (1)
select case (temp_hor_adv)
case (2,3,4,5,6)
......@@ -171,7 +174,7 @@ temp_field_no=1
call getm_error("init_3d()", &
"temp_adv_split=1: temp_ver_adv not valid (2-6)")
end select
LEVEL2 "1D split --> half u, half v, full w, half v, half u"
LEVEL3 "1D split --> half u, half v, full w, half v, half u"
case (2)
select case (temp_hor_adv)
case (2,7)
......@@ -185,7 +188,7 @@ temp_field_no=1
call getm_error("init_3d()", &
"temp_adv_split=2: temp_ver_adv not valid (2-6)")
end select
LEVEL2 "2D-hor, 1D-vert split --> full uv, full w"
LEVEL3 "2D-hor, 1D-vert split --> full uv, full w"
case default
end select
......@@ -193,6 +196,17 @@ temp_field_no=1
call wait_halo(D_TAG)
call mirror_bdy_3d(T,D_TAG)
LEVEL3 'temp_check=',temp_check
if (temp_check .ne. 0) then
LEVEL4 'doing sanity check on temperature'
if (temp_check .gt. 0) then
LEVEL4 'out-of-bound values result in termination of program'
end if
if (temp_check .lt. 0) then
LEVEL4 'out-of-bound values result in warnings only'
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving init_temperature()'
write(debug,*)
......@@ -264,6 +278,7 @@ temp_field_no=1
REALTYPE :: area_inv(I2DFIELD)
REALTYPE :: swr_loc,shf_loc
REALTYPE :: zz,rad(0:1000),A=0.58,g1=0.35,g2=23.0
integer :: status
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -362,6 +377,23 @@ temp_field_no=1
end do
end do
if (temp_check .ne. 0 .and. mod(n,abs(temp_check)) .eq. 0) then
call check_3d_fields(imin,jmin,imax,jmax,az, &
iimin,jjmin,iimax,jjmax,kmax, &
kmin,T,min_temp,max_temp,status)
if (status .gt. 0) then
if (temp_check .gt. 0) then
call getm_error("do_temperature()", &
"out-of-bound values encountered")
end if
if (temp_check .lt. 0) then
LEVEL1 'do_temperature(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving do_temperature()'
write(debug,*)
......
!$Id: uu_momentum_3d.F90,v 1.11 2006-08-25 09:00:19 kbk Exp $
!$Id: uu_momentum_3d.F90,v 1.12 2006-12-15 09:57:50 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -6,7 +6,7 @@
! !ROUTINE: uu_momentum_3d - $x$-momentum eq.\ \label{sec-uu-momentum-3d}
!
! !INTERFACE:
subroutine uu_momentum_3d(bdy3d)
subroutine uu_momentum_3d(n,bdy3d)
!
! !DESCRIPTION:
!
......@@ -41,7 +41,9 @@
! is activated), the result for $j=2$ is copied to $j=3$.
!
! !USES:
use exceptions
use parameters, only: g,avmmol,rho_0
use domain, only: imin,jmin,imax,jmax
use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HU,min_depth
use domain, only: dry_u,coru,au,av,az,ax
#if defined CURVILINEAR || defined SPHERICAL
......@@ -62,9 +64,11 @@
#endif
use halo_zones, only: update_3d_halo,wait_halo,U_TAG
use meteo, only: tausx,airp
use m3d, only: vel_check,min_vel,max_vel
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: n
logical, intent(in) :: bdy3d
!
! !INPUT/OUTPUT PARAMETERS:
......@@ -84,7 +88,7 @@
REALTYPE :: zp,zm,zx,ResInt,Diff,Vloc
REALTYPE :: gamma=g*rho_0
REALTYPE :: cord_curv=_ZERO_
integer :: status
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -216,6 +220,24 @@ end do
call wait_halo(U_TAG)
call mirror_bdy_3d(uu,U_TAG)
if (vel_check .ne. 0 .and. mod(n,abs(vel_check)) .eq. 0) then
call check_3d_fields(imin,jmin,imax,jmax,au, &
iimin,jjmin,iimax,jjmax,kmax, &
kumin,uu,min_vel,max_vel,status)
if (status .gt. 0) then
if (vel_check .gt. 0) then
call getm_error("uu_momentum_3d()", &
"out-of-bound values encountered")
end if
if (vel_check .lt. 0) then
LEVEL1 'do_salinity(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving uu_momentum_3d()'
write(debug,*)
......
!$Id: vv_momentum_3d.F90,v 1.14 2006-08-25 09:00:19 kbk Exp $
!$Id: vv_momentum_3d.F90,v 1.15 2006-12-15 09:57:50 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -6,7 +6,7 @@
! !ROUTINE: vv_momentum_3d - $y$-momentum eq.\ \label{sec-vv-momentum-3d}
!
! !INTERFACE:
subroutine vv_momentum_3d(bdy3d)
subroutine vv_momentum_3d(n,bdy3d)
!
! !DESCRIPTION:
!
......@@ -44,7 +44,9 @@
! prescribed, which has to be hard-coded as local variable.
!
! !USES:
use exceptions
use parameters, only: g,avmmol,rho_0
use domain, only: imin,jmin,imax,jmax
use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HV,min_depth
use domain, only: dry_v,corv,au,av,az,ax
#if defined CURVILINEAR || defined SPHERICAL
......@@ -68,9 +70,11 @@
#endif
use halo_zones, only: update_3d_halo,wait_halo,V_TAG
use meteo, only: tausy,airp
use m3d, only: vel_check,min_vel,max_vel
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: n
logical, intent(in) :: bdy3d
!
! !INPUT/OUTPUT PARAMETERS:
......@@ -93,6 +97,7 @@
#ifdef XZ_PLUME_TEST
REALTYPE :: yslope=0.001
#endif
integer :: status
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -232,6 +237,23 @@
call wait_halo(V_TAG)
call mirror_bdy_3d(vv,V_TAG)
if (vel_check .ne. 0 .and. mod(n,abs(vel_check)) .eq. 0) then
call check_3d_fields(imin,jmin,imax,jmax,av, &
iimin,jjmin,iimax,jjmax,kmax, &
kvmin,vv,min_vel,max_vel,status)
if (status .gt. 0) then
if (vel_check .gt. 0) then
call getm_error("vv_momentum_3d()", &
"out-of-bound values encountered")
end if
if (vel_check .lt. 0) then
LEVEL1 'do_salinity(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving vv_momentum_3d()'
write(debug,*)
......
#$Id: Makefile,v 1.7 2006-02-07 07:16:22 kbk Exp $
#$Id: Makefile,v 1.8 2006-12-15 09:57:48 kbk Exp $
#
# Makefile to build utilities written in Fortran90 - libfutils.a
#
......@@ -13,7 +13,7 @@ LIBSRC = ver_interpol.F90 kbk_interpol.F90 tridiagonal.F90 pos.F90 \
cnv_2d.F90 cnv_3d.F90 eta_mask.F90 col_interpol.F90 \
to_2d_vel.F90 to_2d_u.F90 to_2d_v.F90 \
to_3d_vel.F90 to_3d_uu.F90 to_3d_vv.F90 \
tow.F90 c2x.F90
tow.F90 c2x.F90 check_3d_fields.F90
SRC = $(MODSRC) $(LIBSRC)
......@@ -47,7 +47,8 @@ ${LIB}(to_3d_vel.o) \
${LIB}(to_3d_uu.o) \
${LIB}(to_3d_vv.o) \
${LIB}(tow.o) \
${LIB}(c2x.o)
${LIB}(c2x.o) \
${LIB}(check_3d_fields.o)
all: modules objects
......
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: check_3d_fields() - check 3D-fields for having sane values
!
! !INTERFACE:
subroutine check_3d_fields(imin,jmin,imax,jmax,mask, &
iimin,jjmin,iimax,jjmax,kmax, &
kmin,f,minval,maxval,status)
!
! !DESCRIPTION:
! This routine scans over a 3D field and checks that values are within
! specified bounds. The routine is mainly intended for debugging
! purposes - but can with little overhead also be applied for production
! runs. If status is returned as a non-zero number out-of-bound values
! have been found. It is up to the calling routine to take any actions.
! The out-of-bound values are written to stderr.
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: imin,jmin,imax,jmax
integer, intent(in) :: mask(E2DFIELD)
integer, intent(in) :: iimin,jjmin,iimax,jjmax,kmax
integer, intent(in) :: kmin(I2DFIELD)
REALTYPE, intent(in) :: f(I3DFIELD)
REALTYPE, intent(in) :: minval,maxval
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
integer, intent(out) :: status
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding
!
! $Log: check_3d_fields.F90,v $
! Revision 1.1 2006-12-15 09:57:48 kbk
! optional sanity checks on velocities, temperature and salinity fields
!
!
! !LOCAL VARIABLES:
integer :: i,j,k
!EOP
!-----------------------------------------------------------------------
!BOC
status=0
do k=0,kmax
do j=jjmin,jjmax
do i=iimin,iimax
if ( mask(i,j) .gt. 0 .and. k .ge. kmin(i,j) ) then
if(f(i,j,k) .lt. minval .or. f(i,j,k) .gt. maxval) then
LEVEL1 'check_3d_fields() ',i,j,k,f(i,j,k)
status=status+1
end if
end if
end do
end do
end do
return
end subroutine check_3d_fields
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2006 - Karsten Bolding, Hans Burchard
!-----------------------------------------------------------------------
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment