source: lmdz_wrf/trunk/WRFV3/share/module_compute_geop.F @ 409

Last change on this file since 409 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 2.9 KB
Line 
1MODULE module_compute_geop
2
3CONTAINS
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
109END MODULE module_compute_geop
Note: See TracBrowser for help on using the repository browser.