da_hydrostaticp_to_rho_lin.inc

1 subroutine da_hydrostaticp_to_rho_lin(grid, p, rho ) 
2 
3    !---------------------------------------------------------------------------
4    !  Purpose: Calculates density increment from pressure increment fiteld.
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(domain), intent(in)  :: grid
16    real,         intent(in)  :: p(ims:ime,jms:jme,kms:kme)   ! Pressure inc. (cross pts)
17    real,         intent(out) :: rho(ims:ime,jms:jme,kms:kme) ! Density inc. (cross pts)
18 
19    integer :: i, j, k      ! Loop counters.
20    real    :: delta1       ! Height difference.
21    real    :: delta2       ! Height difference.
22    real    :: dPdz         ! Vertical pressure gradient.
23 
24    if (trace_use) call da_trace_entry("da_hydrostaticp_to_rho_lin")
25 
26    !---------------------------------------------------------------------------
27    ! [2.0] Calculate density increments at all levels except top/bottom:
28    !---------------------------------------------------------------------------
29 
30    do k = kts+1, kte-1 
31       rho(its:ite,jts:jte,k) =  -( p(its:ite,jts:jte,k+1) - p(its:ite,jts:jte,k-1) ) / &
32                              ( ( grid%xb % h(its:ite,jts:jte,k+1) - &
33                              grid%xb % h(its:ite,jts:jte,k-1) ) * gravity )
34    end do                                  
35 
36    !---------------------------------------------------------------------------
37    ! [3.0] Calculate density increment on bottom level:
38    !---------------------------------------------------------------------------
39 
40    k = kts
41    do j = jts, jte
42       do i = its, ite
43          ! dP~/dz by backwards one-sided 2nd order finite differences:
44          
45          delta1 = grid%xb % h(i,j,k+1) - grid%xb % h(i,j,k)
46          delta2 = grid%xb % h(i,j,k+2) - grid%xb % h(i,j,k)
47          dPdz = -( delta1 + delta2 ) * p(i,j,k)  / ( delta1 * delta2 ) + &
48                  ( delta2 / delta1 ) * p(i,j,k+1)/ ( delta2 - delta1 ) - &
49                  ( delta1 / delta2 ) * p(i,j,k+2)/ ( delta2 - delta1 )
50                       
51          ! Put together to get density increment at top level:
52          rho(i,j,k) = -dPdz / gravity
53       end do
54    end do
55                        
56    !---------------------------------------------------------------------------
57    ! [4.0] Calculate density increment on top level:
58    !---------------------------------------------------------------------------
59 
60    k = kte
61    do j = jts, jte
62       do i = its, ite
63          ! dP~/dz by forwards one-sided 2nd order finite differences:
64          
65          delta1 = grid%xb % h(i,j,k) - grid%xb % h(i,j,k-1)
66          delta2 = grid%xb % h(i,j,k) - grid%xb % h(i,j,k-2)
67          
68          dPdz = ( delta1 + delta2 ) * p(i,j,k)   / ( delta1 * delta2 ) - &
69                 ( delta2 / delta1 ) * p(i,j,k-1) / ( delta2 - delta1 ) + &
70                 ( delta1 / delta2 ) * p(i,j,k-2) / ( delta2 - delta1 )
71 
72          ! Put together to get density increment at bottom level:
73          rho(i,j,k) = -dPdz / gravity
74       end do
75    end do
76 
77    if (trace_use) call da_trace_exit("da_hydrostaticp_to_rho_lin")
78    
79 end subroutine da_hydrostaticp_to_rho_lin
80 
81