da_uvprho_to_w_lin.inc
References to this file elsewhere.
1 subroutine da_uvprho_to_w_lin(grid)
2
3 !------------------------------------------------------------------------------
4 ! Purpose: Calculates vertical velocity increments from Richardson's Eq.
5 !
6 ! Method: Richardson's Eq., which
7 ! combines continuity Eq., thermodynamic Eq. and hrdrostatic Eq.
8 !------------------------------------------------------------------------------
9
10 implicit none
11
12 type (domain), intent(inout) :: grid
13
14 integer :: is, ie ! 1st dim. end points.
15 integer :: js, je ! 2nd dim. end points.
16
17 integer :: I,J,K
18
19 real :: urho(ims:ime,jms:jme,kms:kme)
20 real :: vrho(ims:ime,jms:jme,kms:kme)
21 real :: div(ims:ime,jms:jme,kms:kme)
22 real :: wz(ims:ime,jms:jme,kms:kme)
23 real :: term3
24
25 if (trace_use) call da_trace_entry("da_uvprho_to_w_lin")
26
27 ! Computation to check for edge of domain:
28 is = its
29 ie = ite
30 js = jts
31 je = jte
32 if (its == ids) is = ids+1
33 if (ite == ide) ie = ide-1
34 if (jts == jds) js = jds+1
35 if (jte == jde) je = jde-1
36
37 WZ(:,:,:) = 0.0
38 ! Term 1.1: perturbed pressure advection along the basic wind
39 do K=kts,kte
40 do J=js,je
41 do I=is,ie
42 WZ(I,J,K)=WZ(I,J,K)-grid%xb%u(I,J,K)*(grid%xa%p(I+1,J,K)-grid%xa%p(I-1,J,K))* &
43 grid%xb%coefx(I,J)
44 WZ(I,J,K)=WZ(I,J,K)-grid%xb%v(I,J,K)*(grid%xa%p(I,J+1,K)-grid%xa%p(I,J-1,K))* &
45 grid%xb%coefy(I,J)
46 end do
47 end do
48 end do
49
50 ! Term 1.2: Basic pressure advection along the perturbed wind
51
52 do K=kts,kte
53 do J=js,je
54 do I=is,ie
55 WZ(I,J,K)=WZ(I,J,K)-grid%xa%u(I,J,K)*(grid%xb%p(I+1,J,K)-grid%xb%p(I-1,J,K))* &
56 grid%xb%coefx(I,J)
57 WZ(I,J,K)=WZ(I,J,K)-grid%xa%v(I,J,K)*(grid%xb%p(I,J+1,K)-grid%xb%p(I,J-1,K))* &
58 grid%xb%coefy(I,J)
59 end do
60 end do
61 end do
62
63 ! Dealing the laterial boundary because of the advection.
64 ! boundary too simple? (It is the same as fill in interpf, fill can be used)
65
66 if (its == ids) then
67 do K=kts,kte
68 do J=js,je
69 WZ(its,J,K)=WZ(its+1,J,K)
70 end do
71 end do
72 end if
73
74 if (ite == ide) then
75 do K=kts,kte
76 do J=js,je
77 WZ(ite,J,K)=WZ(ite-1,J,K)
78 end do
79 end do
80 end if
81
82 if (jts == jds) then
83 do K=kts,kte
84 do I=its, ite
85 WZ(I,jts,K)=WZ(I,jts+1,K)
86 end do
87 end do
88 end if
89
90 if (jte == jde) then
91 do K=kts,kte
92 do I=its, ite
93 WZ(I,jte,K)=WZ(I,jte-1,K)
94 end do
95 end do
96 end if
97
98 ! Term 2.1: Divergence term from perturbed wind
99
100 call da_uv_to_divergence(grid%xb, grid%xa%u, grid%xa%v, DIV)
101
102 WZ(its:ite,jts:jte,kts:kte)=WZ(its:ite,jts:jte,kts:kte)-GAMMA*grid%xb%p(its:ite,jts:jte,kts:kte)*DIV(its:ite,jts:jte,kts:kte)
103
104 ! Term 2.2: Divergence term from basic wind
105
106 call da_uv_to_divergence(grid%xb, grid%xb%u, grid%xb%v, DIV)
107
108 WZ(its:ite,jts:jte,kts:kte)=WZ(its:ite,jts:jte,kts:kte)-GAMMA*grid%xa%p(its:ite,jts:jte,kts:kte)*DIV(its:ite,jts:jte,kts:kte)
109
110 ! Computation to check for edge of domain:
111 is = its-1; ie = ite+1; js = jts-1; je = jte+1
112 if (its == ids) is = ids; if (ite == ide) ie = ide
113 if (jts == jds) js = jds; if (jte == jde) je = jde
114
115 ! Term 3.1: Vertical integration of the perturbed mass divergence
116
117 URHO(is:ie,js:je,kts:kte)=grid%xb%rho(is:ie,js:je,kts:kte)*grid%xa%u(is:ie,js:je,kts:kte)
118 VRHO(is:ie,js:je,kts:kte)=grid%xb%rho(is:ie,js:je,kts:kte)*grid%xa%v(is:ie,js:je,kts:kte)
119
120 call da_uv_to_divergence(grid%xb, URHO, VRHO, DIV)
121
122 do J=jts,jte
123 do I=its,ite
124 TERM3=0.0
125
126 do K=kte-1,kts,-1
127 TERM3=TERM3+GRAVITY*(DIV(I,J,K+1)+DIV(I,J,K))*0.5 *(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
128 WZ(I,J,K)=WZ(I,J,K)+TERM3
129 end do
130 end do
131 end do
132
133 ! Term 3.2: Vertical integration of the basic mass divergence
134
135 URHO(is:ie,js:je,kts:kte)=grid%xa%rho(is:ie,js:je,kts:kte)*grid%xb%u(is:ie,js:je,kts:kte)
136 VRHO(is:ie,js:je,kts:kte)=grid%xa%rho(is:ie,js:je,kts:kte)*grid%xb%v(is:ie,js:je,kts:kte)
137
138 call da_uv_to_divergence(grid%xb, URHO, VRHO, DIV)
139
140 do J=jts,jte
141 do I=its,ite
142 TERM3=0.0
143
144 do K=kte-1,kts,-1
145 TERM3=TERM3+GRAVITY*(DIV(I,J,K+1)+DIV(I,J,K))*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
146 WZ(I,J,K)=WZ(I,J,K)+TERM3
147 end do
148 end do
149 end do
150
151 ! Term 4: Derivative of basic vertical velocity with respect to z.
152
153 do J=jts,jte
154 do I=its,ite
155 do K=kts,kte
156 WZ(I,J,K)=WZ(I,J,K)-GAMMA*grid%xa%p(I,J,K)*(grid%xb%w(I,J,K+1)-grid%xb%w(I,J,K))/ &
157 (grid%xb%hf(I,J,K+1)-grid%xb%hf(I,J,K))
158 end do
159 end do
160 end do
161
162 ! Divide by constant
163
164 WZ(its:ite,jts:jte,kts:kte)=WZ(its:ite,jts:jte,kts:kte)/(GAMMA*grid%xb%p(its:ite,jts:jte,kts:kte))
165
166 ! integration to calculate the vertical velocity
167
168 call da_w_adjustment_lin(grid%xb,grid%xa%w,WZ)
169
170 do J=jts,jte
171 do I=its,ite
172 grid%xa%w(I,J,kte+1)=0.0
173 do K=kte,kts,-1
174 grid%xa%w(I,J,K)=grid%xa%w(I,J,K+1) + WZ(I,J,K)*(grid%xb%hf(I,J,K)-grid%xb%hf(I,J,K+1))
175 end do
176 end do
177 end do
178
179 if (trace_use) call da_trace_exit("da_uvprho_to_w_lin")
180
181 end subroutine da_uvprho_to_w_lin
182
183