MODULE module_qc
!-----------------------------------------------------------------------------!
! Vertical sounding checks: correct spikes and super-adiabatic rate
! Height checks: flag level above the model lid
!
! HISTORY:
!
! D. GILL, April 1998
! F. VANDENBERGHE, May 2000
! Y.-G. GUO, September 2000
!
! 01/13/2003 - Updated for Profiler obs. S. R. H. Rizvi
!
! 02/04/2003 - Updated for Buoy obs. S. R. H. Rizvi
!
! 02/11/2003 - Reviewed and modified for Profiler
! and Buoy obs. Y.-R. Guo
! 06/30/2006 - Updated for AIRS retrievals Syed RH Rizvi
! 11/09/2006 - Updated for GPS RO Y.-R. Guo
!------------------------------------------------------------------------------
INTEGER :: nsynops_qc (2), nshipss_qc (2), nmetars_qc (2), &
npilots_qc (2), nsounds_qc (2), nsatems_qc (2), &
nsatobs_qc (2), naireps_qc (2), ngpspws_qc (2), &
nssmt1s_qc (2), nssmt2s_qc (2), nssmis_qc (2), &
ntovss_qc (2), nothers_qc (2), namdars_qc (2), &
nqscats_qc (2), nbuoyss_qc (2), nprofls_qc (2), &
ngpsztd_qc (2), ngpsref_qc (2), ngpseph_qc (2), &
nboguss_qc (2), nairss_qc(2) , ntamdar_qc (2)
CONTAINS
!------------------------------------------------------------------------------
SUBROUTINE proc_qc1 ( max_number_of_obs , obs , number_of_obs, &
qc_test_vert_consistency , qc_test_convective_adj , &
print_qc_vert , print_qc_conv)
! Driver routine for QC
!
! This routine is a driver routine for quality control ( or QC )
! of observations. It includes:
! 1. vertical sounding checks
!
! This procedure for quality control was applied before those
! calculations related to the model background fields:
! surface_correction, innovation, etc.
!
! Modified by Yong-Run Guo
! 11/09/00
!
USE module_type
USE module_func
USE module_per_type
IMPLICIT NONE
INTEGER, INTENT ( IN ) :: max_number_of_obs
TYPE (report), DIMENSION (max_number_of_obs):: obs
INTEGER , INTENT ( IN ) :: number_of_obs
LOGICAL , INTENT ( IN ) :: qc_test_vert_consistency, &
qc_test_convective_adj
LOGICAL , INTENT ( IN ) :: print_qc_vert, &
print_qc_conv
TYPE (measurement), POINTER :: current
REAL :: err, err_max
INTEGER :: qc_flag
INTEGER :: i
LOGICAL :: corrected, failed
CHARACTER (LEN=80) :: filename
LOGICAL :: connected
INTEGER :: iunit, io_error
CHARACTER (LEN = 32) :: proc_name = 'proc_qc: '
INCLUDE 'platform_interface.inc'
INCLUDE 'missing.inc'
!------------------------------------------------------------------------------!
! WRITE (0,'(A)') &
!'------------------------------------------------------------------------------'
! WRITE ( UNIT = 0, FMT = '(A)') 'DATA CHECK AND QC:'
! Open diagnostic file
IF (print_qc_vert .OR. print_qc_conv) THEN
filename = 'obs_qc1.diag'
iunit = 999
INQUIRE (UNIT = iunit, OPENED = connected )
IF (connected) CLOSE (iunit)
OPEN (UNIT = iunit , FILE = filename , FORM = 'FORMATTED' , &
ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error )
IF (io_error .NE. 0) THEN
CALL error_handler (proc_name, &
"Unable to open output diagnostic file. ", filename, .TRUE.)
! ELSE
! WRITE (UNIT = 0, FMT = '(A,A,/)') &
! "Diagnostics in file ", TRIM (filename)
ENDIF
ENDIF
! The first quality control (QC) that can be performed is with
! data that is vertically stacked. These reports have the temperature,
! speed and direction compared against reasonable benchmarks.
! Bad values have flags set, though switching the sign of the temperature
! (in degrees C) is allowed. This test is not performed if the entire
! report was discarded, and no error checks are performed on bogus
! data types. Since the vertical consistency check and the dry
! convective adjustment are for vertically stacked data, we can ignore
! these test for data with only one level.
consistency_check: &
IF ( qc_test_vert_consistency ) THEN
WRITE (0,'(A)') &
'------------------------------------------------------------------------------'
WRITE (UNIT = 0, FMT = '(A,/)') 'VERTICAL CONSISTENCY TEST QC:'
IF (print_qc_vert) THEN
WRITE (UNIT = 0, FMT = '(A,A,/)') &
"Diagnostics in file ", TRIM (filename)
ENDIF
loop_all_1 : DO i = 1, number_of_obs
! For GPS Refractivity (FM=116) and GPS Excess Phase (FM=118), no check applied:
if ((obs(i)%info%platform(4:6) == '116') &
.or. (obs(i)%info%platform(4:6) == '118')) cycle loop_all_1
valid_ob_1 : IF ( (.NOT. obs(i)%info%discard ) .AND. &
(.NOT. obs(i)%info%bogus ) .AND. &
(obs(i)%info%is_sound ) .AND. &
(.NOT. obs(i)%location%name=='OSSEs').AND. &
(ASSOCIATED ( obs(i)%surface ) ) ) THEN
failed = .FALSE.
CALL vert_cons_check ( obs ( i ) , i , print_qc_vert, iunit, failed )
IF (failed) THEN
READ (obs (i) % info % platform (4:6), '(I3)') fm
CALL fm_decoder (fm, platform, bogus=nboguss_qc(1), &
synop=nsynops_qc(1), ship =nshipss_qc(1), &
metar=nmetars_qc(1), pilot=npilots_qc(1), &
sound=nsounds_qc(1), satem=nsatems_qc(1), &
satob=nsatobs_qc(1), airep=naireps_qc(1), &
gpspw=ngpspws_qc(1), gpszd=ngpsztd_qc(1), &
gpsrf=ngpsref_qc(1), ssmt1=nssmt1s_qc(1), &
ssmt2=nssmt2s_qc(1), gpsep=ngpseph_qc(1), &
ssmi =nssmis_qc (1), &
tovs =ntovss_qc (1), other=nothers_qc(1), &
amdar=namdars_qc(1), qscat=nqscats_qc(1), &
profl=nprofls_qc(1), buoy =nbuoyss_qc(1), &
airs =nairss_qc(1) , tamdar=ntamdar_qc(1) )
ENDIF
END IF valid_ob_1
END DO loop_all_1
END IF consistency_check
! The second QC check that is available to be performed on the
! the observations is the convective adjustment. This is used
! only for the upper level temperature reports.
convective_check: &
IF (qc_test_convective_adj) THEN
WRITE (0,'(A)') &
'------------------------------------------------------------------------------'
WRITE (UNIT = 0, FMT = '(A,/)') 'CONVECTIVE ADJUSTEMENT TEST QC:'
IF (print_qc_conv) THEN
WRITE (UNIT = 0, FMT = '(A,A,/)') &
"Diagnostics in file ", TRIM (filename)
ENDIF
loop_all_2 : DO i = 1, number_of_obs
! For GPS Refractivity (FM=116) and GPS Excess Phase (FM=118), no check applied:
if ((obs(i)%info%platform(4:6) == '116') &
.or.(obs(i)%info%platform(4:6) == '118')) cycle loop_all_2
valid_ob_2 : IF ( (.NOT. obs(i)%info%discard ) .AND. &
(.NOT. obs(i)%info%bogus ) .AND. &
(obs(i)%info%is_sound ) .AND. &
(.NOT. obs(i)%location%name=='OSSEs').AND. &
(ASSOCIATED ( obs(i)%surface ) ) ) THEN
corrected = .FALSE.
CALL dry_convective_adjustment (obs (i), i, &
print_qc_conv, iunit, corrected)
IF (corrected) THEN
READ (obs (i) % info % platform (4:6), '(I3)') fm
CALL fm_decoder (fm, platform, bogus=nboguss_qc(2), &
synop=nsynops_qc(2), ship =nshipss_qc(2), &
metar=nmetars_qc(2), pilot=npilots_qc(2), &
sound=nsounds_qc(2), satem=nsatems_qc(2), &
satob=nsatobs_qc(2), airep=naireps_qc(2), &
gpspw=ngpspws_qc(2), gpszd=ngpsztd_qc(2), &
gpsrf=ngpsref_qc(2), gpsep=ngpseph_qc(2), &
ssmt1=nssmt1s_qc(2), &
ssmt2=nssmt2s_qc(2), ssmi =nssmis_qc (2), &
tovs =ntovss_qc (2), other=nothers_qc(2), &
amdar=namdars_qc(2), qscat=nqscats_qc(2), &
profl=nprofls_qc(2), buoy =nbuoyss_qc(2), &
airs =nairss_qc(2) , tamdar=ntamdar_qc(2) )
ENDIF
END IF valid_ob_2
END DO loop_all_2
END IF convective_check
! Close diagnostic file
IF (print_qc_vert .OR. print_qc_conv) THEN
CLOSE (iunit)
ENDIF
END SUBROUTINE proc_qc1
!------------------------------------------------------------------------------
SUBROUTINE modify_qc_flag ( qc_flag , err , factored_difference , test_number )
IMPLICIT NONE
INTEGER :: test_number
INTEGER :: qc_flag
REAL :: err, factored_difference
INTEGER :: qc_small, qc_large
! If the difference between observation and analysis at obs location
! (err) exceeds errmx, mark those obs as flagged by the error checking
! routines.
IF ( ABS ( err ) .GT. factored_difference ) THEN
qc_small = mod ( qc_flag , test_number )
qc_large = ( qc_flag / ( test_number * 2 ) ) * 2
qc_flag = qc_large*test_number + qc_small + test_number
ENDIF
END SUBROUTINE modify_qc_flag
!------------------------------------------------------------------------------
SUBROUTINE dry_convective_adjustment (obs , counter , print_dry , iunit, &
corrected)
! This subroutine performs the dry convective adjustment for a sounding
! and prints out the points and levels where the adjustment have been
! invoked. The adjustment removes any super-adiabatic lapse rate in the
! sounding by making use of the conservation of dry static energy.
USE module_type
USE module_func
IMPLICIT NONE
INTEGER :: iunit
INTEGER :: counter, nlevel
INTEGER :: k, n1, kk, klop,& ! parameters for do loops
numb, num1, i, &
kst, & ! kst is the starting level of adjustment
ked ! ked is the ending level of adjustment
REAL :: esum, & ! sum of static energy between 2 adjacent levels
eav, & ! mean static energy between 2 adjacent levels
psfc ! surface pressure
REAL , DIMENSION ( : ) , ALLOCATABLE :: p , & ! pressure
t , & ! temperature
td, & ! dew point
dp, & ! dew point depression
z , & ! height
zt, & ! height calculated
! from temperature
e ! dry static energy
REAL , DIMENSION ( : ) , ALLOCATABLE :: theta ! potential temperature
INTEGER :: index ! number of vertical levels
TYPE ( measurement ) , POINTER :: current
TYPE ( report ) :: obs
CHARACTER ( LEN = 120) :: message
LOGICAL :: print_dry
LOGICAL :: corrected
REAL :: es, qs
INCLUDE 'constants.inc'
! INCLUDE 'error.inc'
INCLUDE 'missing.inc'
! INTERFACE
! INCLUDE 'error.int'
! END INTERFACE
! Try to find out the surface level first
psfc = missing_r
NULLIFY ( current )
nlevel = 0
current => obs%surface
loop_psfc: DO WHILE ( ASSOCIATED ( current ) )
! IF (eps_equal (current%meas%height%data, obs%info%elevation, 0.1)) THEN
IF (.NOT. eps_equal (current%meas%height%data , missing_r, 1.0) .AND. &
.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.0) ) THEN
psfc = current%meas%pressure%data
!
! Write out the stations with (h < elevation):
!
if ( print_dry ) then
if (current%meas%height%data < obs%info%elevation) &
write(iunit,'(/A,A,2X,A,2X,A/A,2I6,3f12.2)') &
"Platform, ID, NAME: ", &
obs % info % platform, &
obs % location % id(1:10), &
obs % location % name, &
"PSFC found ==> I, nlevel, psfc, height, elevation:",&
counter, nlevel, psfc, current % meas % height%data, &
obs % info % elevation
end if !print
EXIT loop_psfc
ELSE
! Do nothing
END IF
current => current%next
END DO loop_psfc
! Find out the number of vertical levels, then ALLOCATE arrays.
NULLIFY ( current )
index = 1
current => obs%surface
DO WHILE ( ASSOCIATED ( current ) )
IF ((.NOT. eps_equal (current%meas%pressure%data, missing_r, 1. )).AND. &
(.NOT. eps_equal (psfc , missing_r, 1. )).AND. &
(current%meas%pressure%data .LE. psfc ) .AND. &
(.NOT. eps_equal (current%meas%temperature%data, missing_r, 1.)))THEN
index = index + 1
ELSE
! Do nothing, missing or bad data
END IF
current => current%next
END DO
index = index - 1
IF ((index .LE. 3) .OR. (eps_equal (psfc, missing_r, 1.))) THEN
! The vertical levels are less than 3, do nothing
ELSE ! Make dry adjustment for temperature
ALLOCATE ( P ( index ) )
ALLOCATE ( t ( index ) )
ALLOCATE ( td( index ) )
ALLOCATE ( dp( index ) )
ALLOCATE ( z ( index ) )
ALLOCATE ( zt( index ) )
ALLOCATE ( e ( index ) )
ALLOCATE ( theta ( index ) )
NULLIFY ( current )
! Read data set into array and perform quality check.
index = 1
current => obs%surface
DO WHILE ( ASSOCIATED ( current ) )
IF ((.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.)).AND.&
( current%meas%pressure%data .LE. psfc) .AND. &
(.NOT. eps_equal (current%meas%temperature%data, missing_r, 1.)))&
THEN
p ( index ) = current%meas%pressure%data
t ( index ) = current%meas%temperature%data
td( index ) = current%meas%dew_point%data
IF (.NOT. eps_equal (current%meas%dew_point%data, missing_r, 1.)) &
THEN
dp (index) = t (index) - td (index)
ELSE
dp (index) = missing_r
ENDIF
z (index) = current%meas%height%data
index = index + 1
ELSE
! Do nothing, missing or bad data
END IF
current => current%next
END DO
index = index - 1
n1 = index - 1
IF (.NOT. eps_equal ( z (1) , missing_r, 1.)) THEN
! We have the first level height, go ahead make dry adjustment.
zt (1) = z (1)
! Compute potential temperature.
DO k = 1, index
theta ( k ) = t ( k ) * ( 100000.0 / p ( k ) ) ** rcp
END DO
loop_adjust: DO k = 1 , index - 2
IF (theta (k) .LE. theta (k+1)) THEN
! Stable layer, do nothing.
CYCLE loop_adjust
ELSE
! Compute the height z at all levels by taking
! a mean temperature of the layer.
DO kk = 2 , n1
zt (kk) = zt (kk - 1) + 0.5 * (t (kk - 1) + t (kk)) * &
(gasr/g) * LOG (p (kk -1) / p (kk))
END DO
! Compute the dry static energy.
DO kk = 1 , n1
e (kk) = g * zt (kk) + cp * t (kk)
END DO
esum = e (k) + e (k + 1)
eav = esum * 0.5
! Downward checking the layer in which superadiabat exists.
numb = 2
loop_down: DO klop = k - 1 , 1 , -1
IF (klop .LT. 1) THEN
EXIT loop_down
ELSE
IF (eav .GE. e (klop) ) THEN
EXIT loop_down
ELSE
esum = esum + e ( klop )
numb = numb + 1
eav = esum / REAL ( numb )
END IF
END IF
END DO loop_down
num1 = numb
! Upward checking the layer in which superadiabat exists.
loop_up: DO klop = k + 2 , index , 1
IF (klop .GT. n1) THEN
EXIT loop_up
ELSE
IF (eav .LE. e (klop)) THEN
EXIT loop_up
ELSE
esum = esum + e (klop)
numb = numb + 1
eav = esum / REAL (numb)
END IF
END IF
END DO loop_up
! Define the starting and ending levels.
kst = k + 2 - num1
ked = kst + numb -1
corrected = .TRUE.
IF (print_dry) THEN
WRITE (message, FMT = '(&
& "Dry convective adjustment", &
& " ID: ",A5,", NAME: ",A8, &
& ", from k = ",I2, " to ",I2, &
& ", from P = ",F6.0," to ",F6.0," Pa." )') &
& obs%location%id, obs%location%name, &
& kst, ked, p(kst)/100., p(ked)/100.
WRITE (iunit,'(A)') message
ENDIF
! Evaluates the new temperature profile according to
! the mean static energy of the layer.
DO kk = kst , ked
t (kk) = ((eav - g * zt (kst)) / cp ) * &
(p (kk) / p (kst)) ** rcp
END DO
END IF
END DO loop_adjust
ELSE
! We do not have the first level height, do nothing.
END IF
current => obs%surface
looptt: DO WHILE ( ASSOCIATED ( current ) )
DO i = 1, index
IF ((eps_equal (p (i), current%meas%pressure%data, 0.1)) .AND. &
(.NOT. eps_equal (t (i), current%meas%temperature%data, 0.1))) THEN
IF (print_dry) THEN
WRITE (message , FMT = &
& '("T dry conv adj: ", &
& " ID = ",A5,", NAME = ",A8, &
& ", T old = " ,F5.1," (K),", &
& ", T new = " ,F5.1," (K),", &
& ", P = ", F6.1," (hPa)")') &
& obs%location%id, obs%location%name, &
& current%meas%temperature%data, t(i), p(i)/100.
WRITE (iunit,'(A)') message
ENDIF
! Correct t, td, rh and qv
IF (.NOT. eps_equal &
(current%meas%temperature%data,missing_r, 1.)) THEN
current%meas%temperature%data = t (i)
ENDIF
IF (.NOT. eps_equal &
(current%meas%dew_point%data, missing_r,1.)) THEN
current%meas%dew_point%data = t ( i ) - dp ( i)
current%meas%rh%data = 100. * exp (5418.12 &
* (1./t(i) - 1./(t(i)-dp(i))))
IF (.NOT.eps_equal (current%meas%qv%data, missing_r, 1.)) THEN
es = 6.112 * EXP (17.67 * (t (i) - 273.15) &
/ (t (i) - 273.15+243.5))
qs = 0.622 * es /(p (i)/100. -es)
current%meas%qv%data = MAX (0.01*current%meas%rh%data*qs,&
0.);
ENDIF
ENDIF
! Update QC
IF (current%meas%temperature%qc .NE. missing ) THEN
current%meas%temperature%qc = convective_adjustment &
+ current%meas%temperature%qc
END IF
IF (current%meas%dew_point%qc .NE. missing ) THEN
current%meas%dew_point%qc = convective_adjustment &
+ current%meas%dew_point%qc
END IF
IF (current%meas%rh%qc .NE. missing ) THEN
current%meas%rh%qc = convective_adjustment &
+ current%meas%rh%qc
END IF
IF (current%meas%qv%qc .NE. missing ) THEN
current%meas%qv%qc = convective_adjustment &
+ current%meas%qv%qc
END IF
ELSE
! Do nothing
END IF
END DO
current => current%next
END DO looptt
DEALLOCATE ( P )
DEALLOCATE ( t )
DEALLOCATE ( td )
DEALLOCATE ( dp )
DEALLOCATE ( z )
DEALLOCATE ( zt )
DEALLOCATE ( e )
DEALLOCATE ( theta )
NULLIFY ( current )
END IF
END SUBROUTINE dry_convective_adjustment
!-------------------------------------------------------------------------------
SUBROUTINE vert_cons_check ( obs , counter , print_vert, iunit, failed )
! This subroutine is called to perform vertical
! consistency check for a single sounding.
USE module_type
USE module_func
IMPLICIT NONE
INTEGER :: i, index, iflag , counter, nlevel
REAL , DIMENSION ( : ), ALLOCATABLE :: p, t, td, ws, wd, u, v, rh, qv, &
work1, work2
INTEGER, DIMENSION ( : ), ALLOCATABLE :: tqc, tdqc, wsqc, wdqc, uqc, vqc, &
rhqc, qvqc, iindex
TYPE ( measurement ), POINTER :: current
TYPE ( report ) :: obs
REAL :: d1, d2, che1, che2, th, th1, th2, &
psfc, depression
REAL :: es, qs
CHARACTER ( LEN = 120 ) :: proc_name = "vert_cons_check"
CHARACTER ( LEN = 120 ) :: message
LOGICAL :: print_vert
REAL :: p1 , p2 , h1 , h2
LOGICAL :: found
LOGICAL :: failed
LOGICAL :: fatal, listing
INTEGER :: iunit
! INCLUDE 'error.inc'
INCLUDE 'constants.inc'
INCLUDE 'missing.inc'
! INTERFACE
! INCLUDE 'error.int'
! END INTERFACE
! Check for spikes in wind curve.
! Both the lapse rate check and potential temperature check will be
! performed.
! Is there a duplicate surface value?
! We can check the pressure and height values.
! If we find two (or more) levels with the height equal to the surface
! elevation,then we set all of the information in the duplicate level to
! missing.
NULLIFY ( current )
current => obs%surface
IF (ASSOCIATED ( current ) ) THEN
p1 = current%meas%pressure%data
h1 = current%meas%height%data
END IF
current => current%next
found = .FALSE.
DO WHILE ( ASSOCIATED ( current ) )
IF (eps_equal (current%meas%height%data, obs%info%elevation, 0.1) &
.AND. .NOT.eps_equal (obs%info%elevation, missing_r, 0.1) ) THEN
p2 = current%meas%pressure%data
h2 = current%meas%height%data
IF (( eps_equal (h1 , h2 , 0.1 )) .AND. &
(.NOT. eps_equal (p1 , p2 , 0.1 ))) THEN
WRITE (message, FMT = '(" Duplicate surface found at ",A8,A8)') &
TRIM (obs%location%id), TRIM (obs%location%name)
fatal = .true.
listing = .false.
! To discard the OBS:
obs%info % discard = .TRUE.
CALL error_handler (proc_name, message, "",fatal)
current%meas%pressure%data = missing_r
current%meas%height%data = missing_r
current%meas%temperature%data = missing_r
current%meas%dew_point%data = missing_r
current%meas%speed%data = missing_r
current%meas%direction%data = missing_r
current%meas%rh%data = missing_r
current%meas%qv%data = missing_r
current%meas%u%data = missing_r
current%meas%v%data = missing_r
current%meas%pressure%qc = missing
current%meas%height%qc = missing
current%meas%temperature%qc = missing
current%meas%dew_point%qc = missing
current%meas%speed%qc = missing
current%meas%direction%qc = missing
current%meas%rh%qc = missing
current%meas%qv%qc = missing
current%meas%u%qc = missing
current%meas%v%qc = missing
END IF
END IF
current => current%next
END DO
! Try to find out the surface level first
psfc = missing_r
NULLIFY (current)
nlevel = 0
current => obs%surface
loop_psfc: DO WHILE (ASSOCIATED (current))
nlevel = nlevel + 1
! IF (eps_equal (current%meas%height%data, obs%info%elevation, 0.1 ) ) THEN
IF (.NOT. eps_equal (current%meas%height%data , missing_r, 1.0) .AND. &
.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.0) ) THEN
psfc = current%meas%pressure%data
!
! Write out the stations with (h < elevation):
!
if ( print_vert ) then
if (current%meas%height%data < obs%info%elevation) &
write(iunit,'(/A,A,2X,A,2X,A/A,2I6,3f12.2)') &
"Platform, ID, NAME: ", &
obs % info % platform, &
obs % location % id(1:10), &
obs % location % name, &
"PSFC found ==> I, nlevel, psfc, height, elevation:",&
counter, nlevel, psfc, current % meas % height%data, &
obs % info % elevation
end if !print
EXIT loop_psfc
ELSE
! Do nothing
END IF
current => current%next
END DO loop_psfc
! Find out the vertical level number, then ALLOCATE arrays
NULLIFY ( current )
index = 1
current => obs%surface
DO WHILE ( ASSOCIATED ( current ) )
IF ((.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.)) .AND. &
(.NOT. eps_equal (psfc, missing_r, 1.)) .AND. &
(current%meas%pressure%data .LE. psfc) .AND. &
(.NOT. eps_equal (current%meas%temperature%data, missing_r, 1.))) THEN
index = index + 1
ELSE
! Do nothing, missing or bad data
END IF
current => current%next
END DO
index = index - 1
IF ((index .LE. 3 ) .OR. (eps_equal (psfc, missing_r, 1.))) THEN
! The vertical levels are less than 3, do nothing
ELSE ! Check superadiabatic and inversion in temperature curve
ALLOCATE ( P ( index ) )
ALLOCATE ( t ( index ) )
ALLOCATE ( td ( index ) )
ALLOCATE ( rh ( index ) )
ALLOCATE ( qv ( index ) )
ALLOCATE ( tqc ( index ) )
ALLOCATE ( tdqc ( index ) )
ALLOCATE ( rhqc ( index ) )
ALLOCATE ( qvqc ( index ) )
ALLOCATE ( iindex ( index ) )
NULLIFY ( current )
! Read data set into array and perform quality check
index = 1
current => obs%surface
DO WHILE ( ASSOCIATED ( current ) )
IF ((.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.)).AND.&
(current%meas%pressure%data .LE. psfc ) .AND. &
(.NOT. eps_equal (current%meas%temperature%data, missing_r, 1.))) &
THEN
p (index) = current%meas%pressure%data
t (index) = current%meas%temperature%data
td (index) = current%meas%dew_point%data
rh (index) = current%meas%rh%data
qv (index) = current%meas%qv%data
tqc (index) = 0
tdqc (index) = 0
rhqc (index) = 0
qvqc (index) = 0
index = index + 1
ELSE
! Do nothing, missing or bad data
END IF
current => current%next
END DO
index = index - 1
loop1: DO i = 1, index
iindex ( i ) = 0
END DO loop1
! First level
che1 = LOG (t (2) / t (1)) / LOG (p (2) / p (1))
che2 = LOG (t (3) / t (2)) / LOG (p (3) / p (2))
IF (((che1 .LT. -1.5 ) .OR. (che1 .GT. 0.9)) .AND. &
((che2 .GT. -1.5 ) .AND. (che2 .LT. 0.9))) iindex (1) = 1
! IF (((che1 .LT. -1.5 ) .OR. (che1 .GT. 0.9)) .AND. &
! ((che2 .LT. -1.5 ) .OR. (che2 .GT. 0.9))) iindex (1) = 0
! Last level
che1 = LOG (t (index - 1) / t (index - 2)) / LOG (p (index - 1) &
/ p (index - 2))
che2 = LOG (t (index ) / t (index - 1)) / LOG (p (index ) &
/ p (index - 1))
IF (((che2 .LT. -0.9) .OR. (che2 .GT. 0.9)) .AND. &
((che1 .GT. -0.9) .AND. (che1 .LT. 0.9))) iindex (index) = 1
! IF (((che2 .LT. -0.9) .OR. (che2 .GT. 0.9)) .AND. &
! ((che1 .LT. -0.9) .OR. (che1 .GT. 0.9))) iindex (index) = 0
loop2: DO i = 2, index - 1
che1 = LOG (t (i + 1) / t (i )) / LOG (p (i + 1) / p (i ))
che2 = LOG (t (i ) / t (i - 1)) / LOG (p (i ) / p (i - 1 ))
th1 = t (i - 1) * (100000.0 / p (i - 1)) ** RCP
th = t (i ) * (100000.0 / p (i )) ** RCP
th2 = t (i + 1) * (100000.0 / p (i + 1)) ** RCP
IF ((((che1 .LT. -0.9) .OR. (che1 .GT. 0.9)) .AND. &
((che2 .LT. -0.9) .OR. (che2 .GT. 0.9))) .OR. &
((th .LT. (th1 - 10.0)) .AND. (th .LT. (th2 - 10.0))) .OR. &
((th .GT. (th1 + 10.0)) .AND. (th .GT. (th2 + 10.0)))) &
iindex ( i ) = 1
IF ((i .GE. 3) .AND. ((che1 .GT. -0.6) .AND. (che1 .LT. 0.6)) .AND. &
((che2 .GT. -0.6 ) .AND. (che2 .LT. 0.6)) .AND. &
((che1 * che2) .LT. 0.0)) THEN
iindex (i) = 2
ENDIF
END DO loop2
loop3: DO i = 2, index - 1
IF (iindex (i) .EQ. 0) THEN
! Do nothing
ELSE
! Change the sign of temperature ( in Celsius )
che1 = LOG (t (i + 1) / (-t (i) + 2. * 273.15)) &
/ LOG (p (i + 1) / p (i))
che2 = LOG ((-t (i) + 2. * 273.15) / t (i - 1)) &
/ LOG (p (i) / p ( i - 1))
th1 = t (i - 1) * (100000.0 / p (i - 1)) ** RCP
th = (-t ( i) + 2. * 273.15 ) * (100000.0 / p (i)) ** RCP
th2 = t (i + 1) * (100000.0 / p (i + 1)) ** RCP
IF (iindex ( i ) .EQ. 1) THEN
IF ((((che1 .LT. -0.6) .OR. (che1 .GT. 0.6)) .OR. &
((che2 .LT. -0.6) .OR. ( che2 .GT. 0.6))) .OR. &
((th .GT. ( th1 + 10.0)) .AND. (th .GT. (th2 + 10.0)))) &
THEN
! Do nothing
ELSE
tqc ( i ) = wrong_t_sign
tdqc ( i ) = wrong_t_sign
rhqc ( i ) = wrong_t_sign
qvqc ( i ) = wrong_t_sign
IF ( .NOT. eps_equal ( td ( i ) , missing_r , 1. ) ) THEN
depression = t(i) - td(i)
END IF
t ( i ) = -t ( i ) + 2. * 273.15
IF ( .NOT. eps_equal (td (i), missing_r, 1.)) THEN
td (i) = t (i) - depression
rh (i) = 100. * exp ( 5418.12 * ( 1./t(i) - 1./td(i)))
END IF
IF (.NOT. eps_equal (qv (i), missing_r, 1.)) THEN
es = 6.112 * EXP (17.67 * (t (i) - 273.15) &
/ (t (i) - 273.15+243.5))
qs = 0.622 * es / (p (i)/100. - es)
qv (i) = MAX (0.01*rh (i)*qs, 0.);
ENDIF
IF (print_vert) THEN
WRITE (message , FMT = '( &
& "Failed super-adiabatic check: ", &
& " ID = ",A5,","," NAME = " ,A8,", Modified", &
& " T = " ,F5.1," (K)",","," P = ",F6.1, " (hPa)", &
& " I =",I3," iindex=",I4)') &
& TRIM (obs%location%id), TRIM (obs%location%name), &
& t (i), p (i)/100., i, iindex(i)
WRITE (iunit, '(A)') message
ENDIF
iindex ( i ) = 0
END IF
ELSE
IF (((che1 .GT. -0.6 ) .AND. (che1 .LT. 0.6)) .AND. &
((che2 .GT. -0.6 ) .AND. (che2 .LT. 0.6)) .AND. &
((che1 * che2) .GT. 0.0)) THEN
tqc ( i ) = wrong_t_sign
tdqc ( i ) = wrong_t_sign
rhqc ( i ) = wrong_t_sign
qvqc ( i ) = wrong_t_sign
IF (.NOT. eps_equal ( td ( i ) , missing_r , 1. ) ) THEN
depression = t(i) - td(i)
END IF
t ( i ) = -t ( i ) + 2. * 273.15
IF (.NOT. eps_equal ( td ( i ) , missing_r , 1. ) ) THEN
td (i) = t(i) - depression
rh (i) = 100. * exp ( 5418.12 * ( 1./t(i) - 1./td(i)))
END IF
IF (print_vert) THEN
WRITE (message , FMT = '( &
& "Failed super-adiabatic check: ", &
& " ID = ",A5,","," NAME = " ,A8,", Modified", &
& " T = " ,F5.1," (K)",","," P = ",F6.1, " (hPa)", &
& " I =",I3," iindex=",I4)') &
& TRIM (obs%location%id), TRIM (obs%location%name), &
& t (i), p (i)/100., i, iindex(i)
WRITE (iunit, '(A)') message
ENDIF
iindex ( i ) = 0
ELSE
iindex (i) = 0
END IF
END IF
END IF
END DO loop3
loop4: DO i = 1, index
IF (iindex (i) .EQ. 1) THEN
IF (print_vert) THEN
WRITE (message , FMT = '( &
& "Failed super-adiabatic check: ", &
& " ID = ",A5,","," NAME = " ,A8,",", &
& " T = " ,F5.1," (K)",","," P = ",F6.1, " (hPa)", &
& " I =",I3," iindex=",I4)') &
& TRIM (obs%location%id), TRIM (obs%location%name), &
& t (i), p (i)/100., i, iindex(i)
WRITE (iunit,'(A)') message
ENDIF
tqc (i) = t_fail_supa_inver
tdqc (i) = t_fail_supa_inver
rhqc (i) = t_fail_supa_inver
qvqc (i) = t_fail_supa_inver
! t (i) = missing_r
! td (i) = missing_r
! rh (i) = missing_r
! qv (i) = missing_r
failed = .TRUE.
END IF
END DO loop4
NULLIFY (current)
! Put the "Quality Flag" and changed data into "obs" data
current => obs%surface
loopthermal : DO WHILE (ASSOCIATED (current))
DO i = 1, index
IF (eps_equal (p (i), current%meas%pressure%data, 0.1)) THEN
current%meas%temperature%data = t (i)
current%meas%dew_point%data = td (i)
current%meas%rh%data = rh (i)
current%meas%qv%data = qv (i)
IF (current%meas%temperature%qc .NE. missing ) THEN
current%meas%temperature%qc = current%meas%temperature%qc &
+ tqc ( i )
END IF
IF (current%meas%dew_point%qc .NE. missing ) THEN
current%meas%dew_point%qc = current%meas%dew_point%qc &
+ tdqc ( i )
END IF
IF (current%meas%rh%qc .NE. missing ) THEN
current%meas%rh%qc = rhqc ( i ) &
+ current%meas%rh%qc
END IF
IF (current%meas%qv%qc .NE. missing ) THEN
current%meas%qv%qc = qvqc (i ) &
+ current%meas%qv%qc
END IF
ELSE
! Do nothing
END IF
END DO
current => current%next
END DO loopthermal
DEALLOCATE ( p )
DEALLOCATE ( t )
DEALLOCATE ( td )
DEALLOCATE ( rh )
DEALLOCATE ( qv )
DEALLOCATE ( tqc )
DEALLOCATE ( tdqc )
DEALLOCATE ( rhqc )
DEALLOCATE ( qvqc )
DEALLOCATE ( iindex )
NULLIFY ( current )
END IF
! End checking thermal data
! Check for spikes in hodograph and mark them in wind data,
! if wind speed is too strong or wind direction is wrong.
! Find out the vertical level number, then ALLOCATE arrays
index = 1
current => obs%surface
DO WHILE (ASSOCIATED ( current ))
IF ((.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.)) .AND.&
(.NOT. eps_equal (psfc, missing_r, 1.)) .AND.&
(current%meas%pressure%data .LE. psfc ) .AND.&
(.NOT. eps_equal (current%meas%speed%data, missing_r, 1.)) .AND.&
(.NOT. eps_equal (current%meas%direction%data,missing_r, 1.))) THEN
index = index + 1
ELSE
! Do nothing, missing or bad data
END IF
current => current%next
END DO
index = index - 1
IF ((index .LE. 3) .OR. (eps_equal (psfc , missing_r, 1.))) THEN
! The vertical levels are less than 3, do nothing
ELSE ! Check wind data
ALLOCATE ( p ( index ) )
ALLOCATE ( u ( index ) )
ALLOCATE ( v ( index ) )
ALLOCATE ( ws ( index ) )
ALLOCATE ( wd ( index ) )
ALLOCATE ( wsqc ( index ) )
ALLOCATE ( wdqc ( index ) )
ALLOCATE ( uqc ( index ) )
ALLOCATE ( vqc ( index ) )
ALLOCATE ( work1 ( index ) )
ALLOCATE ( work2 ( index ) )
NULLIFY ( current )
! Read data set into array and perform quality check
index = 1
current => obs%surface
DO WHILE (ASSOCIATED (current))
IF ((.NOT. eps_equal (current%meas%pressure%data, missing_r, 1.)).AND.&
(current%meas%pressure%data .LE. psfc) .AND.&
(.NOT. eps_equal (current%meas%speed%data, missing_r, 1.)).AND.&
(.NOT. eps_equal (current%meas%direction%data,missing_r, 1.))) THEN
p (index) = current%meas%pressure%data
ws (index) = current%meas%speed%data
wd (index) = current%meas%direction%data
u (index) = current%meas%u%data
v (index) = current%meas%v%data
wsqc (index) = 0
wdqc (index) = 0
uqc (index) = 0
vqc (index) = 0
index = index + 1
ELSE
! Do nothing, missing or bad data
END IF
current => current%next
END DO
index = index - 1
NULLIFY (current)
loop5: DO i = 2, index - 1
d1 = wd (i - 1)
d2 = wd (i + 1)
IF (( d2 - d1) .GT. 180.) d1 = d1 + 360.
IF (( d2 - d1) .LT. -180.) d2 = d2 + 360.
work1 (i) = wd (i) - 0.5 * (d1 + d2) + 0.0001
IF (work1 (i) .LT. -180. ) work1 (i) = work1 (i) + 360.
IF (work1 (i) .GT. 180. ) work1 (i) = work1 (i) - 360.
work2 (i) = ws (i) - 0.5 * (ws (i - 1) + ws (i + 1)) + 0.0001
END DO loop5
loop6: DO i = 3, index - 2
! The wind speed at this level is reasonable compared to up and down levels;
! The wind direction at the current level is reasonable compared to up
! and down levels.
IF (((ABS (work2 (i)) .LT. 50.) .OR. &
(work2 (i - 1) * work2 (i + 1) .LE. 0.0) .OR. &
(MAX (work2 (i + 1) / work2 (i - 1), &
work2 (i - 1) / work2 (i + 1)) .GT. 3.)) .AND. &
(((wd (i) .LE. 360.) .AND. &
(ws (i) * ABS (work1 (i)) .LT. 1600.)) .OR. &
((work1 (i - 1) * work1 (i + 1)) .LE. 0.0 ) .OR. &
(MAX (work1 (i + 1) / work1 (i - 1), &
work1 (i - 1) / work1 (i + 1)) .GT. 3.))) THEN
! wind data is good, do nothing
ELSE
IF (print_vert) THEN
WRITE (message, FMT = '(&
& "Spikes in wind detected: ID = ",A5,", NAME=" ,A8, &
& " Speed = ",F6.1," (m/s), Direction = ",F5.1," (degrees)")')&
& obs%location%id , obs%location%name , ws(i) , wd(i)
WRITE (iunit,'(A)') message
ENDIF
wsqc (i) = wrong_wind_data
wdqc (i) = wrong_wind_data
uqc (i) = wrong_wind_data
vqc (i) = wrong_wind_data
! ws (i) = missing_r
! wd (i) = missing_r
! u (i) = missing_r
! v (i) = missing_r
END IF
END DO loop6
! Put the "Quality Flag" into "obs" data
current => obs%surface
loopwind : DO WHILE (ASSOCIATED (current))
DO i = 1, index
IF (eps_equal (p (i), current%meas%pressure%data, 0.1)) THEN
current%meas%speed%data = ws (i)
current%meas%direction%data = wd (i)
current%meas%u%data = u (i)
current%meas%v%data = v (i)
IF (current%meas%speed%qc .NE. missing) THEN
current%meas%speed%qc = wsqc (i) &
+ current%meas%speed%qc
END IF
IF (current%meas%direction%qc .NE. missing) THEN
current%meas%direction%qc = wdqc (i) &
+ current%meas%direction%qc
END IF
IF ((current%meas%u%qc .NE. missing) .AND. &
(current%meas%v%qc .NE. missing)) THEN
current%meas%u%qc = uqc (i) &
+ current%meas%u%qc
current%meas%v%qc = vqc (i) &
+ current%meas%v%qc
END IF
ELSE
! Do nothing
END IF
END DO
current => current%next
END DO loopwind
DEALLOCATE (p )
DEALLOCATE (u )
DEALLOCATE (v )
DEALLOCATE (ws )
DEALLOCATE (wd )
DEALLOCATE (wsqc )
DEALLOCATE (wdqc )
DEALLOCATE (uqc )
DEALLOCATE (vqc )
DEALLOCATE (work1 )
DEALLOCATE (work2 )
NULLIFY (current)
END IF
! End checking wind data
END SUBROUTINE vert_cons_check
!-------------------------------------------------------------------------------
SUBROUTINE proc_qc2 (max_number_of_obs, obs , number_of_obs, &
qc_test_above_lid, print_height_qc)
!------------------------------------------------------------------------------
USE module_type
USE module_func
USE module_mm5
IMPLICIT NONE
INTEGER, INTENT (in) :: max_number_of_obs
TYPE (report), DIMENSION (max_number_of_obs) :: obs
INTEGER , INTENT (in) :: number_of_obs
LOGICAL, INTENT (in) :: qc_test_above_lid
LOGICAL, INTENT (in) :: print_height_qc
REAL :: hlid
TYPE (measurement), POINTER :: current
INTEGER :: i
LOGICAL :: connected
LOGICAL :: found = .FALSE.
INTEGER :: iunit, io_error
CHARACTER (LEN = 32) :: proc_name = 'proc_qc2: '
CHARACTER (LEN = 80) :: title, fmt_id, fmt_level
CHARACTER (LEN = 80) :: filename
INCLUDE 'missing.inc'
!-------------------------------------------------------------------------------
IF (.NOT. qc_test_above_lid) RETURN
! 1. OPEN DIAGNOSTIC FILE
! =======================
WRITE (0,'(A)') &
'------------------------------------------------------------------------------'
WRITE (UNIT = 0, FMT = '(A,/)') 'CHECK OBS HEIGHT...:'
IF (print_height_qc) THEN
filename = 'obs_qc2.diag'
iunit = 999
INQUIRE (UNIT = iunit, OPENED = connected)
IF (connected) CLOSE (iunit)
OPEN (UNIT = iunit , FILE = filename , FORM = 'FORMATTED' , &
ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error )
IF (io_error .NE. 0) THEN
CALL error_handler (proc_name,&
"Unable to open output diagnostic file. " , filename, .TRUE.)
ELSE
WRITE (UNIT = 0, FMT = '(A,A)') &
"Diagnostics in file ", TRIM (filename)
ENDIF
fmt_level = '(/,A,2(F10.2,2X),A,I8)'
WRITE (UNIT = iunit , FMT = '(A67)', ADVANCE = 'no') filename
ENDIF
! 1. CHECK IF OBSERVATIONS ARE WITHIN THE MODEL VERTICAL DOMAIN
! ==============================================================
! 1.0 Model Lid height
! -----------------
hlid = ref_height (ptop)
loop_all_obs: &
DO i = 1, number_of_obs
! 1.1 Check only valid obs and not gpspw
! ----------------------------------
valid_obs:&
IF ((.NOT. obs (i) % info % discard) .AND. &
(ASSOCIATED (obs (i) % surface))) THEN
! 1.1 Print obs id
! ------------
IF (print_height_qc) THEN
IF (.NOT. found) THEN
fmt_id = '(TL67,A20,A5,1X,A23,2F9.3)'
ELSE
fmt_id = '(//,A20,A5,1X,A23,2F9.3)'
ENDIF
WRITE (UNIT = iunit , FMT = TRIM (fmt_id), ADVANCE = 'no') &
"Found Name and ID = " , &
TRIM (obs (i) % location % id ) , &
TRIM (obs (i) % location % name), &
obs (i) % location % latitude, &
obs (i) % location % longitude
!SDH010706 found = .FALSE.
found = .TRUE.
ENDIF
! 2. Loop over upper levels
! ==========================
current => obs (i) % surface
all_levels:&
DO WHILE (ASSOCIATED (current))
! 2.1 Upper bound
! -----------
! Use Ptop check instead of hlid check (Y.-R. Guo 07/13/2004):
IF (current % meas % pressure % data < ptop) THEN
current % meas % pressure % qc = above_model_lid
if ( print_height_qc ) then
WRITE (UNIT = iunit , FMT = TRIM (fmt_level), ADVANCE = 'no') &
"Model ptop & obs pressure = ", ptop, &
current % meas % pressure % data, &
" qc = ", current % meas % pressure % qc
end if !print
found = .TRUE.
ENDIF
! IF (current % meas % height % data >= hlid) THEN
! current % meas % height % qc = above_model_lid
! WRITE (UNIT = iunit , FMT = TRIM (fmt_level), ADVANCE = 'no') &
! "Model lid & obs height = ", hlid, &
! current % meas % height % data, &
! " qc = ", current % meas % height % qc
! found = .TRUE.
! ENDIF
! 2.2 Next level
! ----------
current => current % next
ENDDO all_levels
ENDIF valid_obs
ENDDO loop_all_obs
IF (print_height_qc) WRITE (0,'(A)') ""
IF (print_height_qc) CLOSE (iunit)
END SUBROUTINE proc_qc2
!-------------------------------------------------------------------------------
END MODULE module_qc