module_compute_geop.F

1 MODULE module_compute_geop
2 
3 CONTAINS
4   SUBROUTINE compute_500mb_height  ( ph, phb, p, pb,                  &
5                                    height,                          &
6                                    ids, ide, jds, jde, kds, kde,    &
7                                    ims, ime, jms, jme, kms, kme,    &
8                                    its, ite, jts, jte, kts, kte    )
9 
10    IMPLICIT NONE
11 
12 
13    !  Input data.
14 
15    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
16                                     ims, ime, jms, jme, kms, kme, &
17                                     its, ite, jts, jte, kts, kte
18 
19    REAL , DIMENSION(  ims:ime , kms:kme, jms:jme ) ,                      &
20                                                INTENT(IN   ) ::           &
21                                                                  ph,      &
22                                                                  phb,     &
23                                                                  pb,      &
24                                                                  p
25 
26    REAL , DIMENSION( ims:ime , jms:jme ) ,    INTENT(  OUT) :: height
27 
28 !  local variables
29 
30    integer :: i,j,k
31    real, dimension(kms:kme) :: pressure,geopotential
32    real :: interp_output
33 
34 !  slow version of code, we'll call interp routine for each column
35 
36    do j = jts, min(jde-1,jte)
37    do i = its, min(ide-1,ite)
38 
39       do k=kds,kde
40         pressure(k) = p(i,k,j) + pb(i,k,j)
41         geopotential(k) = 0.5*( ph(i,k  ,j)+phb(i,k  ,j)  &
42                                +ph(i,k+1,j)+phb(i,k+1,j) )
43       enddo
44 
45       call interp_p( geopotential, pressure, 50000., interp_output,  &
46                      kds,kde-1,kms,kme, i,j )
47 
48       height(i,j) = interp_output/9.81  !  500 mb height in meters
49 
50    enddo
51    enddo
52 
53    end subroutine compute_500mb_height
54 
55 !--------
56 
57   subroutine interp_p(a,p,p_loc,a_interp,ks,ke,kms,kme,i,j)
58   implicit none
59 
60   integer, intent(in) :: ks,ke,kms,kme,i,j
61   real, dimension(kms:kme), intent(in) :: a,p
62   real, intent(in)  :: p_loc
63   real, intent(out) :: a_interp
64 
65 !---  local variables
66 
67   integer :: kp, pk, k
68   real    :: wght1, wght2, dp, pressure
69   character*256 mess
70 
71     kp = ks+1
72     pk = p(kp)
73     pressure = p_loc
74     do while( pk .gt. pressure )
75 
76       kp = kp+1
77 
78       if(kp .gt. ke) then
79         write(mess,*) ' interp too high: pressure, p(ke), i, j = ',pressure,p(ke),i,j
80         write(0,*)'p: ',p
81         call wrf_error_fatal( mess )
82       end if
83  
84       pk = p(kp)
85 
86     enddo
87 
88     dp = p(kp-1) - p(kp)
89     wght2 = (p(kp-1)-pressure)/dp
90     wght1 = 1.-wght2
91 
92     a_interp = wght1*a(kp-1) + wght2*a(kp)
93 
94     end subroutine interp_p
95 
96 END MODULE module_compute_geop