| 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 |
|---|
