da_hydrostaticp_to_rho_adj.inc

References to this file elsewhere.
1 subroutine da_hydrostaticp_to_rho_adj( xb, xp, 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    type (xpose_type), intent(in)  :: xp       ! Dimensions and xpose buffers.
16    real, intent(in)    :: rho(xp%ims:xp%ime,xp%jms:xp%jme,xp%kms:xp%kme) 
17                                                    ! Density inc. (cross pts).
18    real, intent(inout) :: p(xp%ims:xp%ime,xp%jms:xp%jme,xp%kms:xp%kme)   
19                                                    ! Pressure inc. (cross pts)
20 
21    integer               :: is, ie       ! 1st dim. end points.
22    integer               :: js, je       ! 2nd dim. end points.
23    integer               :: ks, ke       ! 3rd dim. end points.
24    integer               :: i, j, k      ! Loop counters.
25    real                  :: delta1       ! Height difference.
26    real                  :: delta2       ! Height difference.
27    real                  :: dPdz         ! Vertical pressure gradient.
28    real                  :: temp(xp%its:xp%ite,xp%jts:xp%jte) ! Temporary array.
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    ! [4.0] Calculate density increment on top level:
43    !---------------------------------------------------------------------------
44 
45    k = ke
46    do j = js, je
47       do i = is, ie
48       
49          ! Put together to get density increment at bottom level:
50          dPdz = -rho(i,j,k) / gravity      
51       
52          ! dP~/dz by forwards one-sided 2nd order finite differences:
53          
54          delta1 = xb % h(i,j,k) - xb % h(i,j,k-1)
55          delta2 = xb % h(i,j,k) - xb % h(i,j,k-2)
56          
57          p(i,j,k) = p(i,j,k) + ( delta1 + delta2 ) * dPdz / &
58                                ( delta1 * delta2 )
59          p(i,j,k-1) = p(i,j,k-1) - ( delta2 / delta1 ) * dPdz / &
60                                    ( delta2 - delta1 )
61          p(i,j,k-2) = p(i,j,k-2) + ( delta1 / delta2 ) * dPdz / &
62                                    ( delta2 - delta1 )
63       end do
64    end do
65 
66    !---------------------------------------------------------------------------
67    ! [3.0] Calculate density increment on top level:
68    !---------------------------------------------------------------------------
69 
70    k = ks
71    do j = js, je
72       do i = is, ie
73 
74          ! Put together to get density increment at top level:
75          dPdz = -rho(i,j,k) / gravity
76 
77          ! dP~/dz by backwards one-sided 2nd order finite differences:
78 
79          delta1 = xb % h(i,j,k+1) - xb % h(i,j,k)
80          delta2 = xb % h(i,j,k+2) - xb % h(i,j,k)
81 
82          p(i,j,k) = p(i,j,k) - ( delta1 + delta2 ) * dPdz / &
83                                ( delta1 * delta2 )
84          p(i,j,k+1) = p(i,j,k+1) + ( delta2 / delta1 ) * dPdz / &
85                                    ( delta2 - delta1 )
86          p(i,j,k+2) = p(i,j,k+2) - ( delta1 / delta2 ) * dPdz / &
87                       ( delta2 - delta1 )
88                       
89       end do
90    end do  
91    
92    !---------------------------------------------------------------------------
93    ! [2.0] Calculate density increments at all levels except top/bottom:
94    !---------------------------------------------------------------------------
95 
96    do k = ke-1, ks+1, -1 
97       temp(is:ie,js:je) = -rho(is:ie,js:je,k) / &
98                         ( ( xb % h(is:ie,js:je,k+1) - xb % h(is:ie,js:je,k-1) ) * &
99                           gravity )
100 
101       p(is:ie,js:je,k-1) = p(is:ie,js:je,k-1) - temp(is:ie,js:je)
102       p(is:ie,js:je,k+1) = p(is:ie,js:je,k+1) + temp(is:ie,js:je)
103    end do                                  
104 
105 end subroutine da_hydrostaticp_to_rho_adj
106 
107