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 !---------------------------------------------------------------------------
29 ! [1.0] Initialise:
30 !---------------------------------------------------------------------------
31
32 one_third = 1.0 / 3.0
33 div = 0.0
34
35 !---------------------------------------------------------------------------
36 ! Computation to check for edge of domain:
37 !---------------------------------------------------------------------------
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 (.not.global) inv_2ds = 0.5 / xb%ds
45
46 !---------------------------------------------------------------------------
47 ! [2.0] Calculate divergence:
48 !---------------------------------------------------------------------------
49
50 do k = kts, kte
51 ! [2.1] Compute fd divergence at interior points:
52
53 if (global) then
54 do j = js, je
55 do i = is, ie
56 div(i,j,k) = xb%coefx(i,j) * (u(i+1,j,k) - u(i-1,j,k)) + &
57 xb%coefy(i,j) * (v(i,j+1,k) - v(i,j-1,k))
58 end do
59 end do
60
61 ! if (global) cycle
62 else
63
64 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)
65 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)
66
67 ! [2.1] Compute fd divergence at interior points:
68
69 do j = js, je
70 do i = is, ie
71 coeff = xb%map_factor(i,j) * xb%map_factor(i,j) * inv_2ds
72 div(i,j,k) = (um(i+1,j) - um(i-1,j) + vm(i,j+1) - vm(i,j-1)) * coeff
73 end do
74 end do
75
76 ! [2.2] Impose zero divergence gradient condition at boundaries:
77
78 ! [2.2.1] Bottom boundaries:
79
80 if (its == ids) then
81 i = its
82 do j = jts, jte
83 div(i,j,k) = one_third * (4.0 * div(i+1,j,k) - div(i+2,j,k))
84 end do
85 end if
86
87 ! [2.2.2] Top boundaries:
88
89 if (ite == ide) then
90 i = ite
91 do j = jts, jte
92 div(i,j,k) = one_third * (4.0 * div(i-1,j,k) - div(i-2,j,k))
93 end do
94 end if
95
96 ! [2.2.3] Left boundaries:
97
98 if (jts == jds) then
99 j = jts
100 do i = its, ite
101 div(i,j,k) = one_third * (4.0 * div(i,j+1,k) - div(i,j+2,k))
102 end do
103 end if
104
105 ! [2.2.4] right boundaries:
106
107 if (jte == jde) then
108 j = jte
109 do i = its, ite
110 div(i,j,k) = one_third * (4.0 * div(i,j-1,k) - div(i,j-2,k))
111 end do
112 end if
113 end if
114 end do
115
116 if (global) then
117 call da_set_boundary_3d(div)
118 end if
119
120 end subroutine da_uv_to_divergence
121
122