da_balance_equation_lin.inc
References to this file elsewhere.
1 subroutine da_balance_equation_lin( xb, xbx, xp, 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 (xb_type), intent(in) :: xb ! First guess structure.
15 type (xbx_type),intent(in) :: xbx ! Header & non-gridded vars.
16 type (xpose_type),intent(inout) :: xp ! Dimensions and xpose buffers.
17
18 real, intent(in) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp.
19 real, intent(in) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp.
20 real, intent(out) :: phi_b(ims:ime,jms:jme,kms:kme) ! Balanced mass increment.
21
22 integer :: i, j, k ! Loop counters.
23 integer :: is, ie ! 1st dim. end points.
24 integer :: js, je ! 2nd dim. end points.
25
26 real, dimension(ims:ime,jms:jme) :: coefx, & ! Multiplicative coefficient.
27 coefy, & ! Multiplicative coefficient.
28 term_x, & ! Balance eqn x term
29 term_y ! Balance eqn y term
30
31 real, dimension(ims:ime,jms:jme,kms:kme) :: del2phi_b ! Del**2 phi_b/M**2
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 = xb%coefx
46 coefy = xb%coefy
47 else if( fg_format == fg_format_wrf) then
48 coefx = 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
59 term_x(ims:ime,jms:jme) = 0.0
60 term_y(ims:ime,jms:jme) = 0.0
61
62 !---------------------------------------------------------------------------
63 ! [2.0] Calculate RHS of balance equation in gridpt space:
64 !---------------------------------------------------------------------------
65
66 ! [2.1] Include geostrophic terms in balance eqn if requested:
67
68 if (balance_type == balance_geo .OR. &
69 balance_type == balance_geocyc ) then
70 call da_balance_geoterm_lin( xb % cori, xb % rho(:,:,k), &
71 u(:,:,k), v(:,:,k), &
72 term_x, term_y)
73 end if
74
75 ! [2.2] Include cyclostrophic terms in balance eqn if requested:
76
77 if (balance_type == balance_cyc .OR. &
78 balance_type == balance_geocyc ) then
79 call da_balance_cycloterm_lin( xb % rho(:,:,k), xb % u(:,:,k), &
80 xb % v(:,:,k), u(:,:,k), v(:,:,k),&
81 xb % coefx(:,:), xb % coefy(:,:), &
82 term_x(:,:), term_y(:,:))
83 end if
84
85 ! [2.3] Take divergence to get Del**2 phi_b/M**2:
86
87 do j = js, je
88 do i = is, ie
89 del2phi_b(i,j,k) = -coefx(i,j)*( term_x(i+1,j) - term_x(i-1,j)) &
90 -coefy(i,j)*( term_y(i,j+1) - term_y(i,j-1))
91 end do
92 end do
93
94 ! [2.4] Del**2 Phi_b boundary conditions:
95
96 if (.not. global .and. its == ids ) del2phi_b(its,js:je,k) = 0.0
97 if (.not. global .and. ite == ide ) del2phi_b(ite,js:je,k) = 0.0
98 if (jts == jds ) del2phi_b(is:ie,jts,k) = 0.0
99 if (jte == jde ) del2phi_b(is:ie,jte,k) = 0.0
100
101 if (.not. global .and. (its == ids .and. jts == jds) ) del2phi_b(its,jts,k) = 0.0
102 if (.not. global .and. (its == ids .and. jte == jde) ) del2phi_b(its,jte,k) = 0.0
103 if (.not. global .and. (ite == ide .and. jts == jds) ) del2phi_b(ite,jts,k) = 0.0
104 if (.not. global .and. (ite == ide .and. jte == jde) ) del2phi_b(ite,jte,k) = 0.0
105
106 end do
107
108 !------------------------------------------------------------------------------
109 ! [3.0] Solve Del**2 phi_b = RHS for phi_b:
110 !------------------------------------------------------------------------------
111
112 call da_solve_poissoneqn_fst( xbx, xp, del2phi_b, phi_b)
113
114 end subroutine da_balance_equation_lin
115
116