da_get_height.inc

1 subroutine da_get_height( input_file, dim1, dim2, dim3, height)
2 
3    !---------------------------------------------------------------------------
4    ! Purpose: Calculates T, RH from input WRF file.
5    !---------------------------------------------------------------------------
6    
7    implicit none
8    
9    character(len=200), intent(in) :: input_file       ! NETCDF file nane.
10    integer, intent(in)   :: dim1, dim2, dim3          ! Dimensions.
11    real, intent(out)     :: height(1:dim1,1:dim2,1:dim3) ! Height.
12 
13    character(len=10)    :: var                       ! Variable to search for.
14    integer               :: k                         ! Loop counter.
15    real                  :: gravity_inv               ! 1/gravity.
16    real                  :: phb(1:dim1,1:dim2)        ! Base state geopotential.
17    real                  :: ph(1:dim1,1:dim2)         ! Perturbation geopotential.
18    real                  :: phfull(1:dim1,1:dim2,1:dim3+1)! Geopotential.
19 
20    gravity_inv = 1.0 / gravity
21 
22    do k = 1, dim3+1
23       var = "PHB"
24       call da_get_field( input_file, var, 3, dim1, dim2, dim3+1, k, phb)
25       var = "PH"
26       call da_get_field( input_file, var, 3, dim1, dim2, dim3+1, k, ph)
27 
28       phfull(:,:,k) = phb + ph ! Calculate full geopotential on full(w) model levels:
29    end do
30 
31    do k = 1, dim3
32       height(:,:,k) = 0.5 *( phfull(:,:,k+1) + phfull(:,:,k)) * gravity_inv
33    end do
34 
35 end subroutine da_get_height
36 
37