| 1 | MODULE module_compute_geop |
|---|
| 2 | |
|---|
| 3 | CONTAINS |
|---|
| 4 | SUBROUTINE compute_500mb_height ( ph, phb, p, pb, & |
|---|
| 5 | height, track_level, & |
|---|
| 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 | INTEGER , INTENT(IN) :: track_level |
|---|
| 29 | |
|---|
| 30 | ! local variables |
|---|
| 31 | |
|---|
| 32 | integer :: i,j,k |
|---|
| 33 | real, dimension(kms:kme) :: pressure,geopotential |
|---|
| 34 | real :: interp_output |
|---|
| 35 | real :: track_level_p |
|---|
| 36 | |
|---|
| 37 | ! slow version of code, we'll call interp routine for each column |
|---|
| 38 | |
|---|
| 39 | track_level_p = float(track_level) |
|---|
| 40 | |
|---|
| 41 | do j = jts, min(jde-1,jte) |
|---|
| 42 | do i = its, min(ide-1,ite) |
|---|
| 43 | |
|---|
| 44 | do k=kds,kde |
|---|
| 45 | pressure(k) = p(i,k,j) + pb(i,k,j) |
|---|
| 46 | geopotential(k) = 0.5*( ph(i,k ,j)+phb(i,k ,j) & |
|---|
| 47 | +ph(i,k+1,j)+phb(i,k+1,j) ) |
|---|
| 48 | enddo |
|---|
| 49 | |
|---|
| 50 | ! call interp_p( geopotential, pressure, 70000., interp_output, & |
|---|
| 51 | call interp_p( geopotential, pressure, track_level_p, interp_output, & |
|---|
| 52 | kds,kde-1,kms,kme, i,j ) |
|---|
| 53 | |
|---|
| 54 | height(i,j) = interp_output/9.81 ! 500 mb height in meters |
|---|
| 55 | |
|---|
| 56 | enddo |
|---|
| 57 | enddo |
|---|
| 58 | |
|---|
| 59 | end subroutine compute_500mb_height |
|---|
| 60 | |
|---|
| 61 | !-------- |
|---|
| 62 | |
|---|
| 63 | subroutine interp_p(a,p,p_loc,a_interp,ks,ke,kms,kme,i,j) |
|---|
| 64 | implicit none |
|---|
| 65 | |
|---|
| 66 | integer, intent(in) :: ks,ke,kms,kme,i,j |
|---|
| 67 | real, dimension(kms:kme), intent(in) :: a,p |
|---|
| 68 | real, intent(in) :: p_loc |
|---|
| 69 | real, intent(out) :: a_interp |
|---|
| 70 | |
|---|
| 71 | !--- local variables |
|---|
| 72 | |
|---|
| 73 | integer :: kp, pk, k |
|---|
| 74 | real :: wght1, wght2, dp, pressure |
|---|
| 75 | character*256 mess |
|---|
| 76 | |
|---|
| 77 | !cys_change: set high value at below-ground point |
|---|
| 78 | if (p(ks).lt.p_loc) then |
|---|
| 79 | a_interp=9.81e5 |
|---|
| 80 | else |
|---|
| 81 | |
|---|
| 82 | kp = ks+1 |
|---|
| 83 | pk = p(kp) |
|---|
| 84 | pressure = p_loc |
|---|
| 85 | do while( pk .gt. pressure ) |
|---|
| 86 | |
|---|
| 87 | kp = kp+1 |
|---|
| 88 | |
|---|
| 89 | if(kp .gt. ke) then |
|---|
| 90 | write(mess,*) ' interp too high: pressure, p(ke), i, j = ',pressure,p(ke),i,j |
|---|
| 91 | write(0,*)'p: ',p |
|---|
| 92 | call wrf_error_fatal( mess ) |
|---|
| 93 | end if |
|---|
| 94 | |
|---|
| 95 | pk = p(kp) |
|---|
| 96 | |
|---|
| 97 | enddo |
|---|
| 98 | |
|---|
| 99 | dp = p(kp-1) - p(kp) |
|---|
| 100 | wght2 = (p(kp-1)-pressure)/dp |
|---|
| 101 | wght1 = 1.-wght2 |
|---|
| 102 | |
|---|
| 103 | a_interp = wght1*a(kp-1) + wght2*a(kp) |
|---|
| 104 | |
|---|
| 105 | endif !cys_change |
|---|
| 106 | |
|---|
| 107 | end subroutine interp_p |
|---|
| 108 | |
|---|
| 109 | END MODULE module_compute_geop |
|---|
