var/da/da_dynamics/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