MODULE module_recoverh
!------------------------------------------------------------------------------!
! Recover observation height, when missing, based on pressure.
!
! Y.-R. GUO, September 2000
!------------------------------------------------------------------------------!
USE module_type
USE module_func
USE module_per_type
USE module_mm5
INCLUDE 'missing.inc'
CONTAINS
!------------------------------------------------------------------------------
SUBROUTINE recover_height_from_pressure(max_number_of_obs , &
obs , number_of_obs, print_hp_recover)
! This routine recovers the missing heights based on the pressure
! under the hydrostatic assumption, or the model reference state.
! for the multi-level OBS data (SOUND, AIREP, etc.).
!
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 ) :: print_hp_recover
TYPE (measurement), POINTER :: current
INTEGER :: iunit
INTEGER :: qc_flag
INTEGER :: i, j, nlevel, k, &
k_start, k_top
CHARACTER (LEN = 80) :: filename
CHARACTER (LEN = 80) :: proc_name = &
"recover_height_from_pressure"
LOGICAL :: connected, correct, failed
INTEGER :: io_error
TYPE (field) , dimension(9000) :: hh
REAL , dimension(9000) :: pp, tt, qq
INCLUDE 'platform_interface.inc'
!------------------------------------------------------------------------------!
WRITE (0,'(A)') &
'------------------------------------------------------------------------------'
WRITE (UNIT = 0, FMT = '(A,/)') 'HEIGHT RECOVERED FROM P, T, Q,..:'
!
! Open diagnostic file
IF (print_hp_recover) THEN
filename = 'obs_recover_height.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
IF (print_hp_recover) &
WRITE (UNIT = IUNIT, FMT = '(/A/)') &
'HEIGHT RECOVERED FROM PRESSURE FOR MULTI-LEVEL OBS DATA:'
failed = .false.
! 1. ESTIMATE H
! ==============
j = 0
! 1.1 Loop over obs
! -------------
loop_all: &
DO i = 1, number_of_obs
IF ((obs (i) % info % discard) .OR. .NOT. ASSOCIATED &
(obs (i) % surface)) THEN
CYCLE loop_all
ENDIF
! 2. SINGLE LEVEL OBS
! ====================
surface:&
IF ((ASSOCIATED (obs (i) % surface)) .AND. &
(.NOT.ASSOCIATED (obs (i) % surface % next))) THEN
! IF height is missing, pressure should be present
IF (eps_equal (&
obs (i) % surface % meas % height % data, missing_r, 1.)) THEN
obs (i) % surface % meas % height % data = ref_height &
(obs (i) % surface % meas % pressure % data)
obs (i) % surface % meas % height % qc = Reference_atmosphere
obs (i) % surface % meas % height % data = NINT &
(obs (i) % surface % meas % height % data + .5)
obs (i) % surface % meas % height % data = MAX &
(obs (i) % surface % meas % height % data, 0.)
IF (print_hp_recover) THEN
WRITE (UNIT = iunit,FMT = '(/,A,A5,1X,A23,2F9.3)') &
"Recover 1 level station id = ", &
TRIM (obs (i) % location % id ) , &
TRIM (obs (i) % location % name), &
obs (i) % location % latitude, &
obs (i) % location % longitude
WRITE (UNIT = iunit, FMT = '(2(A,I5),A)') &
"Use reference state to infer height (", &
INT (obs (i) % surface % meas % height % data), &
"m) from pressure (",&
INT (obs (i) % surface % meas % pressure % data/100.),"hPa)."
ENDIF
ENDIF
! Recover station elevation for single level obs
! IF (eps_equal (&
! obs (i) % info % elevation, missing_r, 1.)) THEN
! obs (i) % info % elevation = &
! obs (i) % surface % meas % height % data
! ENDIF
CYCLE loop_all
ENDIF surface
! 3. MULTI LEVEL OBS
! ==================
! Y.-R. Guo (10/25/2005):.......
call reorder(obs(i), i, 'pressure', failed)
if (failed) then
obs(i) % info % discard = .true.
cycle loop_all
endif
! ..........................................
! 3.1 Get the OBS profile and count the number of levels
! --------------------------------------------------
nlevel = 0
correct = .FALSE.
current => obs(i)%surface
count_level_1:&
DO WHILE (ASSOCIATED (current))
nlevel = nlevel + 1
hh (nlevel) = current%meas%height
pp (nlevel) = current%meas%pressure%data
tt (nlevel) = current%meas%temperature%data
qq (nlevel) = current%meas%qv%data
IF (eps_equal(current%meas%height%data, missing_r, 1.)) THEN
correct = .TRUE.
ENDIF
current => current%next
ENDDO count_level_1
! 3.2 If all levels have height, go to next station
! ---------------------------------------------
IF (.not.correct) CYCLE loop_all
! 3.2 Otherwise recover missing height for upper-level
! ------------------------------------------------
levels:&
IF (nlevel <= 1) THEN
IF (print_hp_recover) THEN
WRITE (UNIT = iunit , FMT = '(A,A5,1X,A23,2F9.3)') &
"No level found for sound id= " , &
TRIM (obs (i) % location % id ) , &
TRIM (obs (i) % location % name), &
obs (i) % location % latitude, &
obs (i) % location % longitude
ENDIF
STOP 'in recover_height.F90'
ELSE IF (nlevel > 1) THEN levels
CALL recover_h_from_ptq (pp, tt, qq, hh, nlevel,k_start,k_top)
IF (k_start >= 1 .or. k_top <= nlevel) THEN
IF (print_hp_recover) &
WRITE (UNIT = iunit, FMT = '(/,A,A5,1X,A23,2F9.3)') &
"Recover upperair station id = ", &
TRIM (obs (i) % location % id ) , &
TRIM (obs (i) % location % name), &
obs (i) % location % latitude, &
obs (i) % location % longitude
ENDIF
ENDIF levels
! 3. CORRECT OBS
! ===============
k = 0
current => obs(i)%surface
correct_levels: &
DO WHILE (ASSOCIATED (current))
k = k + 1
! IF (eps_equal(current%meas%height%data, missing_r, 1.)) THEN
IF (print_hp_recover) THEN
WRITE (UNIT = iunit, FMT = '(2(A,I5),A)') &
"Height missing set to ",INT (hh (k) % data), &
"m, pressure = ",INT (pp(k)/100.),"hpa."
ENDIF
current%meas%height % data = CEILING (hh(k)%data)
current%meas%height % qc = hh(k)%qc
! ENDIF
current => current%next
ENDDO correct_levels
call reorder(obs(i), i, 'pressure', failed)
if (failed) then
obs(i) % info % discard = .true.
cycle loop_all
endif
! 3.4 Go to next station
! ------------------
ENDDO loop_all
IF (print_hp_recover) CLOSE (IUNIT)
END SUBROUTINE recover_height_from_pressure
!-----------------------------------------------------------------------
SUBROUTINE recover_h_from_ptq(P,T,Q,HGT,KXS,K_START,K_TOP)
!-----------------------------------------------------------------------
! To recover the missing heights for the MULTI-LEVEL (KXS>1) OBS data
!
! (1) To compute the heights based on pressure(P), and available
! temperature(T) and specific humidity(Q) under the hydrostatic
! assumption.
!
! Using the available observed heights to calibrate the computed
! heights, then these calibrated and computed heights are used
! where the observed heights are missing.
!
! (2) When the hydrostatic assumption can not be applied at the bottom
! (K_START > 1) or top (K_TOP < KXS) part of the atmosphere, the
! model reference state is used to derive the missing heights based
! on the observed pressure with the calibration from the available
! observed heights.
!
! Note:
!
! In case (1), the derived heights have very high accuracy since only
! a high accurate "hydrostatic assumption" used here.
!
! In case (2), the model reference state parameters are used. These
! parameters are consistent with the model BKG, and independent of
! the model domain settings. Of course, the model reference state
! atmosphere is not a real atmosphere, but it's consistent with the model.
! With the observed heights calibrated, these derived heights
! also included some information from OBS (pressure and height).
!
! This subroutine can not be applied to the SIGLE-LEVEL(KXS=1) data.
! For those OBS, the missing heights can only be derived from the
! observed pressure and the model BKG and reference state (subroutine
! recover_hpt).
!
! Yong-Run Guo
! 11/18/00
!
!-------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(in) :: KXS
REAL , DIMENSION(KXS), INTENT(in) :: T, Q, P
TYPE (field), DIMENSION(KXS), INTENT(inout) :: HGT
INTEGER, INTENT(out) :: k_start,k_top
INTEGER :: k, kk, kwk, L, K0, K1, K2
REAL :: height0, height1, diff_hp, &
TM1, TM2, ALNP, DDH, DDP, DHH, &
dh1, dh2, dp1, dp2, aa, bb, cc
INTEGER,dimension(9000) :: KP
REAL ,DIMENSION(9000) :: TWK, QWK, PWK, HWK, DHGT
LOGICAL :: Vert_ok
TYPE (field), DIMENSION(KXS) :: HGT0, HGT1
include 'constants.inc'
! -------------------------------------------------------------------
HGT0 = HGT
! .. get data at the first levels where both pressure and height
! available
!
! Find the first level with height, pressure and temperature
K_top = kxs+1
K_start = 0
first_level: &
DO k = 1, kxs
IF (.NOT. eps_equal(P (k) , missing_r, 1.) .AND. &
.NOT. eps_equal(HGT(k)%data, missing_r, 1.) .AND. &
.NOT. eps_equal(T (k) , missing_r, 1.)) THEN
K_start = k
EXIT first_level
ENDIF
ENDDO first_level
!
! Obs without temperature and/or height are processed here
!
IF (K_start == 0) THEN
DO k = 1, kxs
IF (.NOT. eps_equal(P (k) , missing_r, 1.) .AND. &
eps_equal(HGT(k)%data, missing_r, 1.)) THEN
HGT(k)%data = ref_height (P (k) )
HGT(k)%qc = reference_atmosphere
ENDIF
ENDDO
! To avoid the duplicated heights between the computed and
! observed heights:
1999 Vert_ok = .True.
DO k = 2, kxs
if ((P(k) < P(K-1)) .and. (HGT(k)%data > HGT(k-1)%data)) then
cycle
else
Vert_ok = .False.
if (HGT(k)%qc <= 0) then
! HGT at level k-1 is computed:
if (k > 2) then
AA = P(k )-P(k-2)
BB = P(k-1)-P(k-2)
CC = AA - BB
HGT(k-1)%data = (HGT(k)%data*BB + HGT(k-2)%data * CC) / AA
! Y.-R. Guo (10/25/2005), Y.-R. Guo (01/16/2006):
if (HGT(k-1)%qc>0) HGT(k-1)%qc = - HGT(k-1)%qc
! print '(" Computed k=",i3," pp,hh,qc:",2f11.2,i8)', k-1,p(k-1),HGT(k-1)%data, HGT(k-1)%qc
else
if ( k <= kxs-1 ) then
AA = P(k+1) - P(k)
BB = P(k-1) - P(k)
CC = AA - BB
HGT(k-1)%data = (HGT(k+1)%data*BB + HGT(k)%data * CC) / AA
else
! Y.-R. Guo (1/31/2008): must be processed separately when k >= kxs
AA = HGT(k)%data - ref_height (P (k) )
HGT(k-1)%data = HGT(k-1)%data + AA
endif
! Y.-R. Guo (10/25/2005), Y.-R. Guo (01/16/2006):
if (HGT(k-1)%qc>0) HGT(k-1)%qc = - HGT(k-1)%qc
endif
else
! HGT at level k is computed
if ( k <= kxs-1 ) then
AA = P(k+1) - P(k-1)
BB = P(k) - P(k-1)
CC = AA - BB
HGT(k)%data = (HGT(k+1)%data*BB + HGT(k-1)%data * CC) / AA
else
! Y.-R. Guo (1/31/2008): must be processed separately when k >= kxs
AA = HGT(k-1)%data - ref_height (P (k-1) )
HGT(k)%data = HGT(k)%data + AA
endif
! Y.-R. Guo (10/25/2005), Y.-R. Guo (01/16/2006):
if (HGT(k-1)%qc>0) HGT(k-1)%qc = - HGT(k-1)%qc
endif
endif
ENDDO
if (.Not. Vert_ok) goto 1999
K_start = 1
K_top = kxs
! Y.-R. Guo (10/25/2005):
DO k = 1, kxs
HGT(k)%qc = abs(HGT(k)%qc)
! print '(i3," p,h,qc:",2f12.2,i8)', k, p(k), HGT(k)%data, HGT(k)%qc
ENDDO
RETURN
ENDIf
! .. if k_start > 1, using the model reference state to get
! the height at any level k below the level k_start,
! keep the height difference between level k_start and level k
! same as between the observed height and the model reference height
! at level k_start
IF (k_start > 1) THEN
! Model reference height
height1 = Ref_height (p(k_start))
DO k = k_start-1, 1, -1
! Missing height
IF (eps_equal (hgt(k)%data, missing_r, 1.)) THEN
IF (p(k)-p(k_start) > 20000.) THEN
! too far below (200mb) the level where the OBS height available
HGT(k)%data = Ref_height(P(k))
HGT(k)%qc = reference_atmosphere
ELSE
HGT(k)%data = HGT(k_start)%data - height1 &
+ Ref_height(P(k))
HGT(k)%qc = reference_OBS_scaled
ENDIF
ENDIF
ENDDO
ENDIF
! Use the hydrostatic equation to correct the heigt
kwk = 0
temp_search: &
DO k = k_start, KXS
IF (.NOT.eps_equal(T(k), missing_r, 1.)) THEN
kwk = kwk+1
pwk (kwk) = P(k)
twk (kwk) = T(k)
qwk (kwk) = q(k)
ENDIF
ENDDO temp_search
HWK(1) = HGT(k_start)%data * G
! ... INTEGRATE HYDROSTATIC EQN
!
hydro_int: &
DO K=2,KWK
ALNP = ALOG (PWK(K-1) /PWK(K)) * gasr
IF (.NOT.eps_equal(QWK(k), missing_r, 1.)) THEN
TM2 = TWK(K )*(1.+0.608*QWK(K ))
ELSE
TM2 = TWK(K )
ENDIF
IF (.NOT.eps_equal(QWK(k-1), missing_r, 1.)) THEN
TM1 = TWK(K-1)*(1.+0.608*QWK(K-1))
ELSE
TM1 = TWK(K-1)
ENDIF
HWK(K) = HWK(K-1) + .5*(TM1+TM2) * ALNP
ENDDO hydro_int
!
! .. CALIBRATION OF THE HWK BASED ON THE AVAILABLE HGT:
K0 = 1
KP(1) = 1
DHGT(1) = 0
calibration: &
DO K = 1,KXS
IF (eps_equal(HGT(K)%data, missing_r, 1.)) then
! nothing
ELSE
DO KK = 1,KWK
IF (P(K).EQ.PWK(KK)) THEN
K0 = K0+1
KP(K0) = KK
DHGT(K0) = HWK(KK)/G - HGT(K)%data
CYCLE calibration
ENDIF
ENDDO
ENDIF
ENDDO calibration
! .. levels between k_start(K0=1) and K0-1
!
DO L = 1,K0-1
K1 = KP(L)
K2 = KP(L+1)
DO KK = K1,K2-1
DDH = DHGT(L+1) - DHGT(L)
DDP = ALOG(PWK(K2)/PWK(K1))
DHH = DHGT(L) + ALOG(PWK(KK)/PWK(K1))*DDH/DDP
HWK(KK) = HWK(KK)/G - DHH
END DO
END DO
!
! .. Above the level KP(K0):
DO K = KP(K0),KWK
HWK(K) = HWK(K)/G - DHGT(K0)
END DO
! WRITE(0,'(I3,2X,4E15.5)') &
! (L,PWK(L),TWK(L),QWK(L),HWK(L),L=1,KWK)
!
! .. FILL BACK THE HEIGHTS AT THE LEVELS WHERE TEMPERATURE AVAILABLE:
height1 = Ref_height(Pwk(kwk))
above_k_start: DO K = k_start,KXS
if (abs(P(K) - PWK(KWK)) < 0.01) k_top = k
! if (eps_equal(HGT(k)%data, missing_r, 1.)) then
IF (P(K) >= PWK(KWK)) then
DO KK = 1,KWK
IF (P(K).EQ.PWK(KK)) THEN
HGT(K)%data = HWK(KK)
ELSE IF (KK.LT.KWK .AND. &
P(K).LT.PWK(KK) .AND. P(K).GT.PWK(KK+1)) THEN
!
! .. Get the interpolated heights at the levels temperature unavailable:
ALNP = ALOG (P(K)/PWK(KK)) / ALOG(PWK(KK+1)/PWK(KK))
HGT(K)%data = HWK(KK) + ALNP*(HWK(KK+1)-HWK(KK))
ENDIF
ENDDO
HGT(k)%qc = Hydrostatic_recover
ELSE
if ((PWK(KWK)-P(K)) > 10000.) then
! too far above (100mb) the level where the OBS height available
HGT(k)%data = Ref_height(P(k))
HGT(k)%qc = reference_atmosphere
else
HGT(k)%data = Hwk(kwk) - height1 &
+ Ref_height(P(k))
HGT(k)%qc = reference_OBS_scaled
endif
ENDIF
! endif
ENDDO above_k_start
!
! Non-hydrastatic adjustment:
!
! In case of only P and H observed without T, TD, the hydrostatic
! height may be different from the observed H. We must keep this
! observed H, and adjust the calculated heights at the adjacent
! levels to avoid inconsistency between P and H when having
! the dense observed levels.
!
! Yong-Run Guo 12/06/2001
! ------------------------------------------------------------
! .. Find the starting level with the observed height:
! Fixed the bug on 07/08/2004
k0 = 1
do k = 1,kxs
if (HGT0(k) % qc == 0) then
k0 = k
exit
endif
enddo
! ------------------------------------------------------------
k1 = 0
! .. Keep the calculated height:
HGT1 = HGT
do k = k0, kxs
! .. Find the starting level (k1+1) without the observed height:
if (HGT0(k) % qc /= 0 .and. k1 == 0) then
k1 = k-1
else if (HGT0(k) % qc == 0 .and. &
abs(HGT0(k)%data - HGT(k)%data) <= 0.10*HGT(k)%data ) then
HGT1(k) = HGT0(k)
endif
! .. Find the ending level (k2-1) without the observed height:
if (HGT0(k) % qc == 0 .and. k1 > 0 .and. k > k1+1) then
k2 = k
if (abs(HGT0(k2)%data - HGT(k2)%data) > 0.10*HGT(k2)%data) then
! ....When the difference between the observed height and retrieved height
! greater than 10% of Hydro. Retrv. height, ignore the observed height,
! and no adjustment done:
HGT0(k) = HGT1(k)
else
HGT1(k) = HGT0(k)
! .. Adjust the calculated hydristatic heights from level k1+1 to k2-1:
dp2 = p(k2) - p(k1)
dh1 = HGT0(k1)%data - HGT(k1)%data
dh2 = HGT0(k2)%data - HGT(k2)%data
do kk = k1+1, k2-1
dp1 = p(kk) - p(k1)
HGT1(kk)%data = HGT(kk)%data + dh1*(1-dp1/dp2) + dh2*dp1/dp2
HGT1(kk)%qc = HGT(kk)%qc
HGT1(kk)%error= HGT(kk)%error
enddo
endif
k1 = 0
k2 = 0
endif
enddo
! .. Fill back to HGT with HGT1:
HGT = HGT1
END subroutine recover_h_from_ptq
!------------------------------------------------------------------------------!
SUBROUTINE reorder(obs, i, order_component, failed)
!------------------------------------------------------------------------------!
IMPLICIT NONE
INTEGER, INTENT (in) :: i
CHARACTER*(*), INTENT (in) :: order_component
TYPE (report), INTENT (inout) :: obs
logical :: failed
INCLUDE 'missing.inc'
TYPE (measurement), pointer :: current, new, tmp, pre
!
integer :: num, ii, num_loop, nlevels
logical :: need_check
failed = .false.
ii = 0
current => obs % surface
count_levels: do while (associated(current))
ii = ii + 1
current => current%next
enddo count_levels
nlevels = ii
current => obs%surface
!--Missing Check:
num = 0
ii = 0
missing_check: do while (associated(current))
num = num + 1
if(eps_equal(comp_field(current,order_component), &
missing_r, 1.0)) then
! remove the missing data link:
num=num-1
end if
current => current%next
end do missing_check
if (num /= nlevels) then
write(0,'(/I5,A,A,A/3X,A,A,A,A,A,2I5,A,f10.3,A,f10.3)') I, &
' There are ',order_component,&
' at several levels missing. Reordering can not be done.', &
' FM=',obs%info%platform(4:5),' ID=',obs%location%id(1:5), &
' nlevels, num:', nlevels, num, &
' lat=', obs%location%latitude, ' long=', obs%location%longitude
! STOP 'in_reorder'
failed = .true.
endif
! Re-set the number of levels to OBS
obs % info % levels = nlevels
if (nlevels <= 1) return
ii = 0
num_loop = 0
need_check = .true.
put2order: do while(need_check)
num_loop = num_loop + 1
head_check: do while(need_check)
current => obs%surface
!--Head Check
if(comp_field(current,order_component) < &
comp_field(current%next,order_component)) then
!--------Cut the first on off and put it in pointer tmp:
tmp => current
current => current%next
nullify(tmp%next)
new => current
obs%surface => new
!--------if first pressure is not good, through tmp away
if(eps_equal(comp_field(current,order_component), &
missing_r, 1.0)) then
num=num-1
! Don't allow happened
STOP 'in reorder_missing data'
! cycle head_check
end if
!--------Now we need to put pointer tmp to a place.
current => obs%surface
new => obs%surface%next
nullify(current%next)
allocate(current%next)
current%next => tmp
tmp%next => new
end if
need_check = .false.
end do head_check
if(num < 3) exit put2order
! write(0,'(/)')
! current => obs%surface
! ii=0
! do while (associated(current))
! ii = ii + 1
! write(0,'(A,I3,A,I3,3f12.2)') 'num_loop=',num_loop, &
! ' I,P,H,T:', ii, &
! current%meas%pressure%data, current%meas%height%data, &
! current%meas%temperature%data
! current => current%next
! end do
!-----Put others in order
pre => obs%surface
current => obs%surface%next
new => obs%surface%next%next
do while (associated(new))
if(comp_field(current,order_component) < &
comp_field(current%next,order_component)) then
tmp => new%next
nullify(pre%next)
nullify(current%next)
nullify(new%next)
current%next => tmp
new%next => current
pre%next => new
need_check = .true.
exit
end if
pre => pre%next
new => new%next
current => current%next
end do
! write(0,'(/)')
! current => obs%surface
! ii=0
! do while (associated(current))
! ii = ii + 1
! write(0,'(A,I3,A,I3,3f12.2)') 'num_loop=',num_loop, &
! ' I,P,H,T:', ii, &
! current%meas%pressure%data, current%meas%height%data, &
! current%meas%temperature%data
! current => current%next
! end do
end do put2order
! write(0,'(//)')
! current => obs%surface
! ii=0
! do while (associated(current))
! ii = ii + 1
! write(0,'(A,I3,3f12.2)') 'I,P,H,T:', ii, &
! current%meas%pressure%data, current%meas%height%data, &
! current%meas%temperature%data
! current => current%next
! end do
END subroutine reorder
!------------------------------------------------------------------------------!
FUNCTION comp_field(current, order_component) result(xxx)
!------------------------------------------------------------------------------!
type (measurement), pointer :: current
character*(*), intent(in) :: order_component
real :: xxx
SELECT CASE (order_component)
CASE ('pressure')
xxx = current%meas%pressure%data
CASE ('height')
xxx = current%meas%height%data
CASE DEFAULT
WRITE(0,'(A,A,A)') 'order_component=',order_component, &
' is not defined correctly'
STOP 'in_reorder'
END SELECT
END function comp_field
!------------------------------------------------------------------------------!
END MODULE module_recoverh