MODULE module_recoverp USE module_type USE module_func USE module_mm5 INCLUDE 'missing.inc' CONTAINS !------------------------------------------------------------------------------ SUBROUTINE recover_pressure_from_height (max_number_of_obs , obs , & number_of_obs, print_hp_recover) !------------------------------------------------------------------------------ ! ! Since the data merging is based on the observed pressure (see ! module_obs_merge.F90/subroutine link_level), the observed pressure ! is not allowed to be missing. ! ! The pressure missing from the decoded OBS sounding (SOUND, AIREP) data ! is recovered ! (1) from the observed heights available with the hydrostatic assumption ! (dP/dH)_hydr if possible, ! (2) or from the observed heights available with the model reference ! state (dP/dH)_ref when hydrostatic interpolation is impossible. ! ! There are two steps to do this: ! ! (1) To guarantee the original pressure and height are correct, first, ! a sequence check, i.e. ! height at the next level > height at the current level and ! pressure at the next level < pressure at the current level, ! is completed (subroutine hp_sequence_check). ! ! (2) second, a consistency check between the observed pressure and ! height agaist the model reference state is completed ! (subroutine pressure_recover, FUNCTION hp_consist). ! ! (3) To fill the missing pressure by using the observed pressure and ! heights that passed the above checks (subroutine recover_p_from_h). ! ! Yong-Run Guo ! 11/16/00 ! !------------------------------------------------------------------------------ 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 :: qc_flag INTEGER :: i, nlevel CHARACTER (LEN=80) :: filename LOGICAL :: connected, consistent INTEGER :: iunit, io_error CHARACTER (LEN = 32) :: proc_name = & "recover_pressure_from_height" INCLUDE 'platform_interface.inc' !------------------------------------------------------------------------------! WRITE (0,'(A)') & '------------------------------------------------------------------------------' WRITE (UNIT = 0, FMT = '(A,/)') 'PRESSURE RECOVERED FROM HEIGHTS:' IF (print_hp_recover) THEN filename = 'obs_recover_pressure.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 ! 2. ESTIMATE P ! ============== ! 2.1 Print parameters ! ---------------- IF (print_hp_recover) THEN WRITE (UNIT = IUNIT, FMT = '(A,/)') 'PRESSURE RECOVERED FROM HEIGHTS:' WRITE(UNIT = iunit , FMT = '(A,/,3(A,F10.2,A),/,3(A,F10.2,A),/)') & "MODEL REFERENCE STATE: ", & " HTOP = ", htop,"m, ", & " PTOP = ", ptop,"Pa, ", & " PS0 = ", ps0 ,"Pa, ", & " TS0 = ", ts0 ,"K, ", & " TLP = ", tlp ,"K, ",& " DIS = ", DIS(idd),"m." ENDIF ! 2.3 Loop over obs ! ------------- loop_all_1:& DO i = 1, number_of_obs IF ((obs (i) % info % discard) .OR. .NOT. ASSOCIATED & (obs (i) % surface)) THEN CYCLE loop_all_1 ENDIF ! 3. SINGLE LEVEL OBS ! ==================== levels:& IF (obs (i) % info % bogus .OR. (.NOT. obs (i) % info % is_sound)) THEN ! IF pressure is missing, height should be present IF (eps_equal & (obs (i) % surface % meas % pressure % data, missing_r, 1.))THEN obs (i) % surface % meas % pressure % data = floor (ref_pres & (obs (i) % surface % meas % height % data)) obs (i) % surface % meas % pressure % qc = & reference_atmosphere ELSE IF (obs (i) % surface % meas % pressure % qc == 0 .and. & obs (i) % surface % meas % height % qc == 0 .and. & eps_equal & (obs (i) % surface % meas % height % data, 0., 1.)) then ! When both h and p have qc = 0 and h = 0. (at sea level ! because the pressure is MSLP for surface data, check ! the consistency between h and p consistent = hp_consist ( & obs (i) % surface % meas % height % data, & obs (i) % surface % meas % pressure % data, & .true., iunit) if (.not.consistent) obs (i) % info % discard = .true. ENDIF ELSE levels ! 3. MULTI LEVEL OBS ! ================== ! 3.1 Height and pressure sequence check ! ---------------------------------- CALL hp_sequence_check (obs (i), i, nlevel, print_hp_recover, iunit) CALL hp_missing_check (obs(i), i, nlevel) ! 3.2 Single level sounding are removed ! --------------------------------- IF ((nlevel == 1) .AND. eps_equal & (obs (i) % surface % meas % pressure% data, missing_r, 1.0)) THEN obs (i) % info % discard = .TRUE. IF (print_hp_recover) THEN WRITE (UNIT = iunit, FMT = '(A,A,I5,A,A)') & "In recover_pressure_from_height: ","I = ", I," ID = ", & obs (i) % location % id (1:5) WRITE (UNIT = iunit, FMT = '(A,I5,A,F9.0,A,F9.0)') & "nlevel = ", nlevel, & "height = ", obs (i) % surface % meas % height % data, & "pressure = ",obs (i) % surface % meas % pressure % data ENDIF CYCLE loop_all_1 ENDIF ! 3.2 If the actual number of valid levels differs from expected, print it ! -------------------------------------------------------------------- IF (nlevel /= obs (i) % info % levels) THEN IF (print_hp_recover) THEN WRITE (UNIT = iunit, FMT = '(3A,I4,A,I4)') & "ID = ", obs(i)%location%id(1:5), ", expect ", & obs (i) % info % levels, " levels, get ", nlevel ENDIF ENDIF ! 3.3 Pressure is recovered based on height if the pressure is missing ! ---------------------------------------------------------------- CALL recover_p_from_h (obs (i), i, nlevel, print_hp_recover, iunit) ENDIF levels ENDDO loop_all_1 IF (print_hp_recover) CLOSE (iunit) END SUBROUTINE recover_pressure_from_height !------------------------------------------------------------------------------! SUBROUTINE hp_sequence_check (obs, i, nlevel, print_hp_recover, iunit) !------------------------------------------------------------------------------! INTEGER, INTENT (in) :: i, iunit INTEGER, INTENT (out) :: nlevel LOGICAL, INTENT (in) :: print_hp_recover TYPE (report), INTENT (inout) :: obs TYPE (measurement), POINTER :: current REAL :: height1, height2, & press1 , press2 LOGICAL :: first_h, first_p !------------------------------------------------------------------------------ current => obs % surface first_h = .false. first_p = .false. nlevel = 0 loop_hp_level: & DO WHILE (ASSOCIATED (current)) nlevel = nlevel + 1 ! 1. HEIGHT CHECK ! ================ IF (.NOT. eps_equal (current%meas%height%data, missing_r, 1.0)) then IF (.NOT. first_h) then height1 = current%meas%height%data first_h = .true. ELSE height2 = current%meas%height%data IF (height2 <= height1) then current%meas%height%qc = missing current%meas%height%data = missing_r IF (print_hp_recover) THEN WRITE (UNIT = iunit, FMT = '(A,/,4A,A,I3,2(A,f12.2))') & "HEIGHT VIOLATION: ", & " FM = ", obs % info % platform (4:6), & " ID = ", obs % location % id, & " LEVEL = ", nlevel, & " H2 = ", height2, & " H1 = ", height1 ENDIF ELSE height1 = height2 ENDIF ENDIF ENDIF ! 2. PRESSURE CHECK ! ================== IF (.NOT.eps_equal(current%meas%pressure%data, missing_r, 1.0)) then IF (.NOT.first_p) THEN press1 = current%meas%pressure%data first_p = .true. ELSE PRESS2 = current%meas%pressure%data IF (press2 >= press1) THEN current%meas%pressure%data = missing_r current%meas%pressure%qc = missing IF (print_hp_recover) THEN WRITE (UNIT = iunit, FMT = '(A,/,4A,A,I3,2(A,f12.2))') & "PRESSURE VIOLATION: ", & " FM = ", obs % info % platform (4:6), & " ID = ", obs % location % id, & " LEVEL = ", nlevel, & " P2 = ", press1, & " P1 = ", press1 ENDIF ELSE press1 = press2 ENDIF ENDIF ENDIF ! Next level current => current%next ENDDO loop_hp_level END subroutine hp_sequence_check !------------------------------------------------------------------------------! SUBROUTINE recover_p_from_h(obs, i, nlevel, print_hp_recover, iunit) !------------------------------------------------------------------------------! IMPLICIT NONE INTEGER, INTENT (in) :: i, iunit, nlevel LOGICAL, INTENT (in) :: print_hp_recover TYPE (report), INTENT (inout) :: obs TYPE (measurement), POINTER :: current TYPE (field) , dimension(nlevel) :: pp, hh INTEGER :: k, Lb, Le REAL :: elev LOGICAL :: do_it !------------------------------------------------------------------------------! ! 1. GET THE P/H PROFILE ! ======================= k = 0 current => obs % surface loop_level:& DO WHILE (ASSOCIATED (current)) k = k + 1 hh (k) = current % meas % height pp (k) = current % meas % pressure current => current % next ENDDO loop_level ! 2. PRESSURE FILED BASED ON THE HEIGHT IF THE PRESSURE IS MISSING ! ================================================================= ! 2.1 Recover pressure ! ---------------- CALL pressure_recover (nlevel, pp, hh, iunit, LB, LE, do_it) ! 2.2 Print status ! ------------ IF (do_it) THEN IF (print_hp_recover) THEN WRITE(UNIT = iunit,FMT = '(5A)', ADVANCE = 'no') & "== PRESSURE RECOVER DONE: ", & " for FM=", obs % info % platform (4:6), & " ID=", obs % location % id ENDIF IF (print_hp_recover) THEN IF (LB > 1 .OR. LE < nlevel) THEN WRITE(UNIT = iunit,FMT = '(A)') & " Reference pressure may have been used." ELSE WRITE(UNIT = iunit,FMT = '(A)') " " ENDIF ENDIF ! 2. SEND BACK TO OBS ! ==================== current => obs % surface k = 0 back_level: & DO WHILE (ASSOCIATED (current)) k = k + 1 ! current % meas % height % data = CEILING (hh (k) % data) current % meas % height % data = hh (k) % data current % meas % height % qc = hh (k) % qc current % meas % pressure % data = FLOOR (pp (k) % data) current % meas % pressure % qc = pp (k) % qc current => current % next ENDDO back_level ENDIF END SUBROUTINE recover_p_from_h !------------------------------------------------------------------------------! FUNCTION hp_consist (hin,pin,print_out,iunit) RESULT (hout) !------------------------------------------------------------------------------! IMPLICIT NONE REAL, INTENT (in) :: hin, pin LOGICAL, INTENT (in) :: print_out INTEGER, INTENT (in) :: iunit LOGICAL :: hout REAL, parameter :: hmax = 1000. ! a tolerance apart from ! the model reference state ! (50000/pin)*hmax REAL :: h_ref, hdiff !------------------------------------------------------------------------------! if (pin <= 0) then write(unit=iunit, fmt='(A,f12.2)') "Pressure voilation: P=",pin STOP 55555 endif hout = .true. h_ref = Ref_height(pin) hdiff = ABS (hin - h_ref) IF (hdiff > (50000/pin)*hmax) THEN IF (print_out) THEN WRITE (UNIT = iunit, FMT = '(/,A,/,3(A,F12.3,/))') & " Pressure / height inconsistency: ", & " Pressure = ", pin, & " height = ", hin, & " ref_height = ", h_ref ENDIF hout = .false. ENDIF END FUNCTION hp_consist !------------------------------------------------------------------------------! subroutine pressure_recover(level, pp, hh, iunit, LB, LE, do_it) IMPLICIT NONE INTEGER, INTENT(in) :: level, iunit TYPE (field), dimension(level), INTENT(inout) :: hh, pp INTEGER, INTENT(out) :: LB, LE LOGICAL, INTENT(out) :: do_it INTEGER :: k, L, L1, L2, Lstart REAL :: height1, height2, & press1 , press2 , bb, diff_pr LOGICAL :: first_L, second_L, & consist1, consist2 ! ----------------------------------------------------------------- LB = 0 LE = 0 do_it = .false. L1 = 0 L2 = 0 first_L = .false. second_L = .false. consist1 = .true. consist2 = .true. ! if there is ONLY one level OBS: if (level == 1) then return endif loop_level: do k = 1,level if (.NOT.first_L) then ! To find the first level with both height and pressure if (.NOT.eps_equal (pp(k)%data, missing_r, 1.0) .and. & .NOT.eps_equal (hh(k)%data, missing_r, 1.0) ) then L1 = k LB = L1 Lstart = L1 first_L = .true. endif else ! To find the second level with both height and pressure if (.NOT.eps_equal (pp(k)%data, missing_r, 1.0) .and. & .NOT.eps_equal (hh(k)%data, missing_r, 1.0) ) then L2 = k second_L = .true. endif endif if (first_L) then height1 = hh(L1) % data press1 = pp(L1) % data if (Lstart > 1) then ! If the first level is not the level 1 and pp(1) is missing if (eps_equal(pp(1)%data, missing_r, 1.0)) then L2 = L1 second_L = .true. L1 = 1 Lstart = 1 height1 = hh(L1) % data ! Set the correction from reference pressure (h1,h2) plus p2 press1 = Ref_pres(height1) press2 = Ref_pres(hh(L2)%data) diff_pr = press1 - press2 press1 = pp(L2)%data + diff_pr pp(L1)%data = press1 pp(L1)%qc = reference_atmosphere else do L = 1, L2 ! There are pressure missing below the level where both p and h available: if (eps_equal(pp(L)%data, missing_r, 1.0)) then pp(L)%data = Ref_pres(hh(L)%data) pp(L)%qc = reference_atmosphere endif enddo endif do_it = .true. endif ! height/pressure consistency check consist1 = hp_consist (height1,press1,.true.,iunit) if (.NOT.consist1) then WRITE(UNIT=IUNIT, FMT='(A,2F12.2)') & "FAILED IN H/P CONSISTENCY CHECK: H1,P1:",height1,press1 first_L = .FALSE. hh(L1) % data = missing_r hh(L1) % qc = missing consist1 = .true. CYCLE loop_level endif endif if (second_L) then height2 = hh(L2) % data press2 = pp(L2) % data consist2 = hp_consist (height2,press2,.true.,iunit) if (.NOT.consist2) then WRITE(UNIT=IUNIT, FMT='(A,I4,2F12.2)') & "FAILED IN H/P CONSISTENCY CHECK: L2,H2,P2:",L2,height2,press2 second_L = .FALSE. hh(L2) % data = missing_r hh(L2) % qc = missing consist2 = .true. CYCLE loop_level endif ! If L1 and L2 are consecutive if (L2 <= (L1 + 1)) then L1 = L2 second_L = .false. else ! If L1 and L2 are not consecutive BB = (height2 - height1)/ALOG(press2/press1) do L = L1+1, L2-1 ! .. get the interpolated pressure from the height if pressure is missing: if (eps_equal(pp(L)%data, missing_r, 1.0)) then if (eps_equal(hh(L)%data, missing_r, 1.0)) then WRITE(0,'(/A,A)') & "Both pressure and height are missing, ", & "this should never be happened???" STOP 33333 endif pp(L)%data = press1*EXP((hh(L)%data - height1)/BB) pp(L)%qc = 0 do_it = .true. endif enddo L1 = L2 second_L = .false. endif endif enddo loop_level ! If L2 is not equal to level LE = L2 if (level > L2) then do k = L2+1, level if (eps_equal(pp(k)%data, missing_r, 1.0)) then ! Set the reference pressure if (L2 > 0) then press1 = Ref_pres(hh(L2)%data) press2 = Ref_pres(hh(k )%data) diff_pr = press1 - press2 pp(k)%data = pp(L2)%data - diff_pr else pp(k)%data = Ref_pres(hh(k )%data) endif pp(k)%qc = reference_atmosphere do_it = .true. endif enddo endif end subroutine pressure_recover SUBROUTINE hp_missing_check (obs, i, nlevel) !------------------------------------------------------------------------------! INTEGER, INTENT (in) :: i INTEGER, INTENT (out) :: nlevel TYPE (report), INTENT (inout) :: obs TYPE (measurement), POINTER :: current, tmp, new, pre integer :: nn !------------------------------------------------------------------------------ current => obs % surface new => current % next nn = 0 loop_hp_level: & DO WHILE (ASSOCIATED (new)) nn = nn + 1 tmp => new % next ! If both P and H missing at the next level? IF (eps_equal (new % meas % pressure % data, missing_r, 1.) .and. & eps_equal (new % meas % height % data, missing_r, 1.) ) THEN write(0,'("==> Missing P and H:",i3,2f15.6)') nn, & new%meas%pressure%data, new%meas%height%data nullify(current % next) current % next => tmp 1001 continue tmp => tmp%next ! If both P and H missing at the next next level? if (associated (current%next) .and. & eps_equal (current%next%meas%pressure% data, missing_r, 1.) .and.& eps_equal (current%next%meas% height % data, missing_r, 1.) ) THEN nullify(current % next) current % next => tmp goto 1001 endif current => current % next ELSE current => current % next ENDIF new => current % next END DO loop_hp_level ! To count the total number of levels: nn = 0 current => obs % surface do while (associated(current)) nn = nn + 1 ! write(0,'(I3," P,H:",2f15.5)') nn, & ! current%meas%pressure%data, current%meas%height%data current => current%next enddo nlevel = nn if (obs%info%levels /= nlevel ) then obs%info%levels = nlevel write(0, fmt='(a, i3, 2x, 2a, 2x, a, 2f10.3/)') & 'After missing check: Number of levels = ', obs%info%levels, & 'OBS=', obs%info%platform(1:12), & 'LOC=', obs%location%latitude, obs%location%longitude endif END SUBROUTINE hp_missing_check END MODULE module_recoverp