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