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