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 do k = kts, kte
53 ! [2.1] Compute fd divergence at interior points:
54
55 if (global) then
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
63 ! if (global) cycle
64 else
65
66 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)
67 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)
68
69 ! [2.1] Compute fd divergence at interior points:
70
71 do j = js, je
72 do i = is, ie
73 coeff = xb%map_factor(i,j) * xb%map_factor(i,j) * inv_2ds
74 div(i,j,k) = (um(i+1,j) - um(i-1,j) + vm(i,j+1) - vm(i,j-1)) * coeff
75 end do
76 end do
77
78 ! [2.2] Impose zero divergence gradient condition at boundaries:
79
80 ! [2.2.1] Bottom boundaries:
81
82 if (its == ids) then
83 i = its
84 do j = jts, jte
85 div(i,j,k) = one_third * (4.0 * div(i+1,j,k) - div(i+2,j,k))
86 end do
87 end if
88
89 ! [2.2.2] Top boundaries:
90
91 if (ite == ide) then
92 i = ite
93 do j = jts, jte
94 div(i,j,k) = one_third * (4.0 * div(i-1,j,k) - div(i-2,j,k))
95 end do
96 end if
97
98 ! [2.2.3] Left boundaries:
99
100 if (jts == jds) then
101 j = jts
102 do i = its, ite
103 div(i,j,k) = one_third * (4.0 * div(i,j+1,k) - div(i,j+2,k))
104 end do
105 end if
106
107 ! [2.2.4] right boundaries:
108
109 if (jte == jde) then
110 j = jte
111 do i = its, ite
112 div(i,j,k) = one_third * (4.0 * div(i,j-1,k) - div(i,j-2,k))
113 end do
114 end if
115 end if
116 end do
117
118 if (global) then
119 call da_set_boundary_3d(div)
120 end if
121
122 if (trace_use) call da_trace_exit("da_uv_to_divergence")
123
124 end subroutine da_uv_to_divergence
125
126