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