da_uv_to_divergence.inc
References to this file elsewhere.
1 subroutine da_uv_to_divergence(xb, u, v, div)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculate divergence on a co-ordinate surface, given an input
5 ! wind field.
6 !
7 ! d U d V
8 ! Div = m^2 *[---(---) + ---(---) ]
9 ! dx m dy M
10 !---------------------------------------------------------------------------
11
12 implicit none
13
14 type (xb_type), intent(in) :: xb ! First guess structure.
15 real, intent(in) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp.
16 real, intent(in) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp.
17 real, intent(inout):: div(ims:ime,jms:jme,kms:kme) ! Divergence.
18
19 integer :: i, j, k ! Loop counters.
20 integer :: is, ie ! 1st dim. end points.
21 integer :: js, je ! 2nd dim. end points.
22 real :: one_third ! 1/3.
23
24 real :: coeff, inv_2ds
25 real :: um(ims:ime,jms:jme) ! Temp. storage of u/m.
26 real :: vm(ims:ime,jms:jme) ! Temp. storage of v/m.
27
28 if (trace_use) call da_trace_entry("da_uv_to_divergence")
29
30 !---------------------------------------------------------------------------
31 ! [1.0] Initialise:
32 !---------------------------------------------------------------------------
33
34 one_third = 1.0 / 3.0
35 div = 0.0
36
37 !---------------------------------------------------------------------------
38 ! Computation to check for edge of domain:
39 !---------------------------------------------------------------------------
40
41 is = its; ie = ite; js = jts; je = jte
42 if (.not. global .and. its == ids) is = ids+1
43 if (.not. global .and. ite == ide) ie = ide-1
44 if (jts == jds) js = jds+1; if (jte == jde) je = jde-1
45
46 if (.not.global) inv_2ds = 0.5 / xb%ds
47
48 !---------------------------------------------------------------------------
49 ! [2.0] Calculate divergence:
50 !---------------------------------------------------------------------------
51
52 if (global) then
53 do k = kts, kte
54 ! [2.1] Compute fd divergence at interior points:
55
56 do j = js, je
57 do i = is, ie
58 div(i,j,k) = xb%coefx(i,j) * (u(i+1,j,k) - u(i-1,j,k)) + &
59 xb%coefy(i,j) * (v(i,j+1,k) - v(i,j-1,k))
60 end do
61 end do
62 end do
63 call da_set_boundary_3d(div)
64 else
65 do k = kts, kte
66
67 um(is-1:ie+1,js-1:je+1) = u(is-1:ie+1,js-1:je+1,k) / xb%map_factor(is-1:ie+1,js-1:je+1)
68 vm(is-1:ie+1,js-1:je+1) = v(is-1:ie+1,js-1:je+1,k) / xb%map_factor(is-1:ie+1,js-1:je+1)
69
70 ! [2.1] Compute fd divergence at interior points:
71
72 do j = js, je
73 do i = is, ie
74 coeff = xb%map_factor(i,j) * xb%map_factor(i,j) * inv_2ds
75 div(i,j,k) = (um(i+1,j) - um(i-1,j) + vm(i,j+1) - vm(i,j-1)) * coeff
76 end do
77 end do
78
79 ! [2.2] Impose zero divergence gradient condition at boundaries:
80
81 ! [2.2.1] Bottom boundaries:
82
83 if (its == ids) then
84 i = its
85 do j = jts, jte
86 div(i,j,k) = one_third * (4.0 * div(i+1,j,k) - div(i+2,j,k))
87 end do
88 end if
89
90 ! [2.2.2] Top boundaries:
91
92 if (ite == ide) then
93 i = ite
94 do j = jts, jte
95 div(i,j,k) = one_third * (4.0 * div(i-1,j,k) - div(i-2,j,k))
96 end do
97 end if
98
99 ! [2.2.3] Left boundaries:
100
101 if (jts == jds) then
102 j = jts
103 do i = its, ite
104 div(i,j,k) = one_third * (4.0 * div(i,j+1,k) - div(i,j+2,k))
105 end do
106 end if
107
108 ! [2.2.4] right boundaries:
109
110 if (jte == jde) then
111 j = jte
112 do i = its, ite
113 div(i,j,k) = one_third * (4.0 * div(i,j-1,k) - div(i,j-2,k))
114 end do
115 end if
116 end do
117 end if
118
119 if (trace_use) call da_trace_exit("da_uv_to_divergence")
120
121 end subroutine da_uv_to_divergence
122
123