da_hydrostaticp_to_rho_lin.inc
References to this file elsewhere.
1 subroutine da_hydrostaticp_to_rho_lin( xb, xp, p, rho )
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculates density increment from pressure increment field.
5 !
6 ! Method: Hydrostatic eqn => rho~ = -dP~/dz / g.
7 !
8 ! Assumptions: 1) Hydrostatic pressure.
9 ! 2) Model level stored top down.
10 !
11 !---------------------------------------------------------------------------
12
13 implicit none
14
15 type (xb_type), intent(in) :: xb ! First guess structure.
16 type (xpose_type), intent(in) :: xp ! Dimensions and xpose buffers.
17 real, intent(in) :: p(xp%ims:xp%ime,xp%jms:xp%jme,xp%kms:xp%kme)
18 ! Pressure inc. (cross pts)
19 real, intent(out) :: rho(xp%ims:xp%ime,xp%jms:xp%jme,xp%kms:xp%kme)
20 ! Density inc. (cross pts)
21
22 integer :: is, ie ! 1st dim. end points.
23 integer :: js, je ! 2nd dim. end points.
24 integer :: ks, ke ! 3rd dim. end points.
25 integer :: i, j, k ! Loop counters.
26 real :: delta1 ! Height difference.
27 real :: delta2 ! Height difference.
28 real :: dPdz ! Vertical pressure gradient.
29
30 !---------------------------------------------------------------------------
31 ! [1.0] Initialise:
32 !---------------------------------------------------------------------------
33
34 is = xp%its
35 ie = xp%ite
36 js = xp%jts
37 je = xp%jte
38 ks = xp%kts
39 ke = xp%kte
40
41 !---------------------------------------------------------------------------
42 ! [2.0] Calculate density increments at all levels except top/bottom:
43 !---------------------------------------------------------------------------
44
45 do k = ks+1, ke-1
46 rho(is:ie,js:je,k) = -( p(is:ie,js:je,k+1) - p(is:ie,js:je,k-1) ) / &
47 ( ( xb % h(is:ie,js:je,k+1) - &
48 xb % h(is:ie,js:je,k-1) ) * gravity )
49 end do
50
51 !---------------------------------------------------------------------------
52 ! [3.0] Calculate density increment on bottom level:
53 !---------------------------------------------------------------------------
54
55 k = ks
56 do j = js, je
57 do i = is, ie
58 ! dP~/dz by backwards one-sided 2nd order finie differences:
59
60 delta1 = xb % h(i,j,k+1) - xb % h(i,j,k)
61 delta2 = xb % h(i,j,k+2) - xb % h(i,j,k)
62 dPdz = -( delta1 + delta2 ) * p(i,j,k) / ( delta1 * delta2 ) + &
63 ( delta2 / delta1 ) * p(i,j,k+1)/ ( delta2 - delta1 ) - &
64 ( delta1 / delta2 ) * p(i,j,k+2)/ ( delta2 - delta1 )
65
66 ! Put together to get density increment at top level:
67 rho(i,j,k) = -dPdz / gravity
68 end do
69 end do
70
71 !---------------------------------------------------------------------------
72 ! [4.0] Calculate density increment on top level:
73 !---------------------------------------------------------------------------
74
75 k = ke
76 do j = js, je
77 do i = is, ie
78 ! dP~/dz by forwards one-sided 2nd order finite differences:
79
80 delta1 = xb % h(i,j,k) - xb % h(i,j,k-1)
81 delta2 = xb % h(i,j,k) - xb % h(i,j,k-2)
82
83 dPdz = ( delta1 + delta2 ) * p(i,j,k) / ( delta1 * delta2 ) - &
84 ( delta2 / delta1 ) * p(i,j,k-1) / ( delta2 - delta1 ) + &
85 ( delta1 / delta2 ) * p(i,j,k-2) / ( delta2 - delta1 )
86
87 ! Put together to get density increment at bottom level:
88 rho(i,j,k) = -dPdz / gravity
89 end do
90 end do
91
92 end subroutine da_hydrostaticp_to_rho_lin
93
94