ICAOHEIGHT.f

      SUBROUTINE icaoheight(MAXWP,     & !input
                            maxwicaoz)   ! output

! SUBPROGRAM:    ICAOHEIGHT      
!   PRGRMMR: CHUANG         ORG: W/NP2      DATE: 09-05-08
!     
! ABSTRACT:
!     THIS ROUTINE CONVERT PRESSURE FIELDS TO ICAO HEIGHT
!

! Description:
!   Convert pressure (Pa) to height (m) using ICAO standard atmosphere
!   adpt code from uk
! Method:
!
!
! Code Description:
!   Language:           Fortran 90
!   Software Standards: UMDP3 v6

use ctlblk_mod, only: jsta, jend, spval, im, ista, iend
use physcons_post, only: con_g, con_rd
IMPLICIT None

! Subroutine Arguments:
!REAL, INTENT(IN) :: SPVAL
REAL, INTENT(IN) :: maxwp(ista:iend,jsta:jend)       !P field for conversion

REAL, INTENT(INOUT) :: maxwicaoz(ista:iend,jsta:jend)   !ICAO height in m
!INTEGER, INTENT(INOUT) :: ErrorStatus

! Local Constants:
REAL, PARAMETER :: g_over_r = con_g / con_rd
REAL, PARAMETER :: lapse_ratel = 6.5e-03  ! For levels below 11,000 gpm
REAL, PARAMETER :: lapse_rateu = -1.0e-03 ! For levels above 11,000 gpm
REAL, PARAMETER :: press_bot = 101325.    ! ICAO std: surface pressure
REAL, PARAMETER :: press_mid = 22632.     !      pressure @ 11,000 gpm
REAL, PARAMETER :: press_top = 5474.87    !      pressure @ 20,000 gpm
REAL, PARAMETER :: temp_bot = 288.15      ! Surface temperature
REAL, PARAMETER :: temp_top = 216.65      ! Temperature of isotherm
REAL, PARAMETER :: gpm1 = 11000.0  ! Ht limit (gpm) for std lower
                                   ! lapse rate
REAL, PARAMETER :: gpm2 = 20000.0  ! Ht (gpm) of top of isothermal layer
REAL, PARAMETER :: zp1 = lapse_ratel/g_over_r ! Exponents used for
REAL, PARAMETER :: zp2 = lapse_rateu/g_over_r ! calculation

! Local Variables:
INTEGER :: i,j      ! loop counters
REAL :: pressure    ! Local pressure

! End of header --------------------------------------------------------

!IF ( ErrorStatus /= StatusOK ) THEN
  ! Previous error - do not proceed
!  GO TO 9999
!END IF

!IF ( ASSOCIATED( IHField % RData ) ) THEN
!  DEALLOCATE( IHField % RData )
!END IF



DO j=jsta,jend
  DO i=ista,iend
    pressure = maxwp(i,j)
    IF ( (pressure <= 1000.) .AND. (pressure >= 0.) ) THEN
      pressure = 1000.
      print*,'lower ICAO pressure to 10 mb'
    END IF
    IF ( pressure > press_bot .and. pressure<spval) THEN
      pressure = press_bot
    END IF

    IF (pressure == spval) THEN
      maxwicaoz(i,j) = spval
    ELSE IF (pressure > press_mid) THEN ! Hts up to 11,000 GPM
      pressure = pressure/press_bot
      pressure = 1.0 - pressure**zp1
      maxwicaoz(i,j) = pressure*temp_bot/lapse_ratel

    ELSE IF (pressure > press_top) THEN ! Hts between 11,000
                                        !     and     20,000 GPM
      pressure = pressure/press_mid
      pressure = -alog(pressure)
      maxwicaoz(i,j) = gpm1 + pressure*temp_top/g_over_r

    ELSE                                ! Hts above 20,000 GPM
      pressure = pressure/press_top
      pressure = 1.0 - pressure**zp2
      maxwicaoz(i,j) = gpm2 + pressure*temp_top/lapse_rateu

    END IF

  ENDDO
ENDDO


!9999 CONTINUE

END SUBROUTINE icaoheight