da_balance_equation_adj.inc
References to this file elsewhere.
1 subroutine da_balance_equation_adj(grid, xbx, phi_b, u, v)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Adjoint of da_balance_equation
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type(domain), intent(inout) :: grid
10
11 type (xbx_type),intent(in) :: xbx ! Header & non-gridded vars.
12 real, intent(in) :: phi_b(ims:ime,jms:jme,kms:kme) ! Balanced mass increment.
13 real, intent(inout) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp. (dot pts)
14 real, intent(inout) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp. (dot pts)
15
16 integer :: i, j, k ! Loop counters.
17 integer :: is, ie ! 1st dim. end points.
18 integer :: js, je ! 2nd dim. end points.
19
20 real, dimension(ims:ime,jms:jme) :: coefx, & ! Multiplicative coefficient.
21 coefy, & ! Multiplicative coefficient.
22 term_x, & ! Balance eqn x term
23 term_y ! Balance eqn y term
24
25 real, dimension(ims:ime,jms:jme,kms:kme) :: del2phi_b ! Del**2 phi_b/M**2
26 real :: coeff1, coeff2 ! Multiplicative coefficient.
27
28 if (trace_use) call da_trace_entry("da_balance_equation_adj")
29
30 !---------------------------------------------------------------------------
31 ! [1.0] Initialise:
32 !---------------------------------------------------------------------------
33
34 ! Computation to check for edge of domain:
35 is = its-1; ie = ite+1; js = jts-1; je = jte+1
36 if (.not. global .and. its == ids) is = ids+1
37 if (.not. global .and. ite == ide) ie = ide-1
38 if (jts == jds ) js = jds+1
39 if (jte == jde ) je = jde-1
40
41 if (fg_format == fg_format_kma_global) then
42 coefx(is:ie,js:je) = grid%xb%coefx(is:ie,js:je)
43 coefy(is:ie,js:je) = grid%xb%coefy(is:ie,js:je)
44 else if (fg_format == fg_format_wrf) then
45 coefx(is:ie,js:je) = grid%xb%coefz(is:ie,js:je)
46 coefy(is:ie,js:je) = coefx(is:ie,js:je)
47 else
48 write (unit=message(1),fmt='(A,I3)') ' Wrong FG_FORMAT = ',fg_format
49 call da_error(__FILE__,__LINE__,message(1:1))
50 end if
51
52 ! [1.1] Multiplicative coefficent for conversion RHS->Del**2 phi_b/M**2:
53
54 del2phi_b(:,:,:) = 0.0
55
56 !---------------------------------------------------------------------------
57 ! [3.0] Solve Del**2 phi_b = RHS for phi_b:
58 !---------------------------------------------------------------------------
59
60 call da_solve_poissoneqn_fst_adj(grid,xbx, phi_b, del2phi_b)
61
62 !---------------------------------------------------------------------------
63 ! [2.0] Calculate RHS of balance equation in gridpt space:
64 !---------------------------------------------------------------------------
65
66 do k = kts, kte
67 ! [2.4] Del**2 Phi_b boundary conditions (null as zero boundary conditions):
68
69 ! [2.3] Take divergence to get Del**2 phi_b/M**2:
70
71 term_x(ims:ime,jms:jme) = 0.0
72 term_y(ims:ime,jms:jme) = 0.0
73
74 do j = je, js, -1
75 do i = ie, is, -1
76 coeff1 = coefx(i,j) * del2phi_b(i,j,k)
77 coeff2 = coefy(i,j) * del2phi_b(i,j,k)
78
79 term_x(i+1,j) = term_x(i+1,j) - coeff1
80 term_x(i-1,j) = term_x(i-1,j) + coeff1
81 term_y(i,j+1) = term_y(i,j+1) - coeff2
82 term_y(i,j-1) = term_y(i,j-1) + coeff2
83 end do
84 end do
85
86 ! [2.2] Include cyclostrophic terms in balance eqn if requested:
87
88 if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
89 call da_balance_cycloterm_adj (grid%xb%rho(:,:,k),grid%xb%u(:,:,k),&
90 grid%xb%v(:,:,k), u(:,:,k), v(:,:,k), grid%xb%coefx(:,:), grid%xb%coefy(:,:),&
91 term_x(:,:), term_y(:,:))
92 end if
93
94
95 ! [2.1] Calculate geostrophic terms in balance eqn:
96
97 if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
98 call da_balance_geoterm_adj (grid%xb%cori, grid%xb%rho(:,:,k), term_x, term_y, &
99 u(:,:,k), v(:,:,k))
100 end if
101 end do
102
103 if (trace_use) call da_trace_exit("da_balance_equation_adj")
104
105 end subroutine da_balance_equation_adj
106
107