var/da/da_dynamics/da_hydrostaticp_to_rho_adj.inc

   1 subroutine da_hydrostaticp_to_rho_adj (xb, rho, p)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Adjoint of da_hydrostaticp_to_rho.
   5    !
   6    !  Method:  Standard adjoint coding.
   7    !
   8    !  Assumptions: 1) Hydrostatic pressure.
   9    !               2) Model level stored top down.
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    type(xb_type), intent(in)    :: xb                           ! First guess structure.
  15    real,          intent(in)    :: rho(ims:ime,jms:jme,kms:kme) ! Density inc. (cross pts).
  16    real,          intent(inout) :: p(ims:ime,jms:jme,kms:kme)   ! Pressure inc. (cross pts)
  17
  18    integer               :: i, j, k      ! Loop counters.
  19    real                  :: delta1       ! Height difference.
  20    real                  :: delta2       ! Height difference.
  21    real                  :: dPdz         ! Vertical pressure gradient.
  22    real                  :: temp(its:ite,jts:jte) ! Temporary array.
  23
  24    if (trace_use) call da_trace_entry("da_hydrostaticp_to_rho_adj")
  25
  26    !---------------------------------------------------------------------------
  27    ! [4.0] Calculate density increment on top level:
  28    !---------------------------------------------------------------------------
  29
  30    do j = jts, jte
  31       do i = its, ite
  32
  33          ! Put together to get density increment at bottom level:
  34          dPdz = -rho(i,j,kte) / gravity
  35
  36          ! dP~/dz by forwards one-sided 2nd order finite differences:
  37
  38          delta1 = xb % h(i,j,kte) - xb % h(i,j,kte-1)
  39          delta2 = xb % h(i,j,kte) - xb % h(i,j,kte-2)
  40
  41          p(i,j,k) = p(i,j,kte) + ( delta1 + delta2 ) * dPdz / (delta1 * delta2)
  42          p(i,j,k-1) = p(i,j,kte-1) - (delta2 / delta1) * dPdz / (delta2 - delta1)
  43          p(i,j,k-2) = p(i,j,kte-2) + (delta1 / delta2) * dPdz / (delta2 - delta1)
  44       end do
  45    end do
  46
  47    !---------------------------------------------------------------------------
  48    ! [3.0] Calculate density increment on top level:
  49    !---------------------------------------------------------------------------
  50
  51    do j = jts, jte
  52       do i = its, ite
  53
  54          ! Put together to get density increment at top level:
  55          dPdz = -rho(i,j,kts) / gravity
  56
  57          ! dP~/dz by backwards one-sided 2nd order finite differences:
  58
  59          delta1 = xb % h(i,j,kts+1) - xb % h(i,j,kts)
  60          delta2 = xb % h(i,j,kts+2) - xb % h(i,j,kts)
  61
  62          p(i,j,kts)   = p(i,j,k)     - (delta1 + delta2) * dPdz / (delta1 * delta2)
  63          p(i,j,kts+1) = p(i,j,kts+1) + (delta2 / delta1) * dPdz / (delta2 - delta1)
  64          p(i,j,kts+2) = p(i,j,kts+2) - (delta1 / delta2) * dPdz / (delta2 - delta1)
  65       end do
  66    end do
  67
  68    !---------------------------------------------------------------------------
  69    ! [2.0] Calculate density increments at all levels except top/bottom:
  70    !---------------------------------------------------------------------------
  71
  72    do k = kte-1, kts+1, -1
  73       temp(its:ite,jts:jte) = -rho(its:ite,jts:jte,k) / &
  74          ((xb%h(its:ite,jts:jte,k+1) - xb%h(its:ite,jts:jte,k-1)) * gravity)
  75
  76       p(its:ite,jts:jte,k-1) = p(its:ite,jts:jte,k-1) - temp(its:ite,jts:jte)
  77       p(its:ite,jts:jte,k+1) = p(its:ite,jts:jte,k+1) + temp(its:ite,jts:jte)
  78    end do
  79
  80 end subroutine da_hydrostaticp_to_rho_adj
  81
  82