Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_force_scm.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 25.0 KB

Line
1	MODULE module_force_scm
2
3	! AUTHOR: Josh Hacker (NCAR/RAL)
4	! Forces a single-column (3x3) version of WRF
5
6	CONTAINS
7
8	SUBROUTINE force_scm(itimestep, dt, scm_force, dx, num_force_layers &
9	, scm_th_adv, scm_qv_adv &
10	, scm_ql_adv &
11	, scm_wind_adv, scm_vert_adv &
12	, scm_soilT_force, scm_soilQ_force &
13	, scm_force_th_largescale &
14	, scm_force_qv_largescale &
15	, scm_force_ql_largescale &
16	, scm_force_wind_largescale &
17	, u_base, v_base, z_base &
18	, z_force, z_force_tend &
19	, u_g, v_g &
20	, u_g_tend, v_g_tend &
21	, w_subs, w_subs_tend &
22	, th_upstream_x, th_upstream_x_tend &
23	, th_upstream_y, th_upstream_y_tend &
24	, qv_upstream_x, qv_upstream_x_tend &
25	, qv_upstream_y, qv_upstream_y_tend &
26	, ql_upstream_x, ql_upstream_x_tend &
27	, ql_upstream_y, ql_upstream_y_tend &
28	, u_upstream_x, u_upstream_x_tend &
29	, u_upstream_y, u_upstream_y_tend &
30	, v_upstream_x, v_upstream_x_tend &
31	, v_upstream_y, v_upstream_y_tend &
32	, tau_x, tau_x_tend &
33	, tau_y, tau_y_tend &
34	,th_largescale &
35	,th_largescale_tend &
36	,qv_largescale &
37	,qv_largescale_tend &
38	,ql_largescale &
39	,ql_largescale_tend &
40	,u_largescale &
41	,u_largescale_tend &
42	,v_largescale &
43	,v_largescale_tend &
44	,tau_largescale &
45	,tau_largescale_tend &
46	, num_force_soil_layers, num_soil_layers &
47	, soil_depth_force, zs &
48	, tslb, smois &
49	, t_soil_forcing_val, t_soil_forcing_tend &
50	, q_soil_forcing_val, q_soil_forcing_tend &
51	, tau_soil &
52	, z, z_at_w, th, qv, ql, u, v &
53	, thten, qvten, qlten, uten, vten &
54	, ids, ide, jds, jde, kds, kde &
55	, ims, ime, jms, jme, kms, kme &
56	, ips, ipe, jps, jpe, kps, kpe &
57	, kts, kte &
58	)
59
60	! adds forcing to bl tendencies and also to base state/geostrophic winds.
61
62	USE module_init_utilities, ONLY : interp_0
63	IMPLICIT NONE
64
65
66	INTEGER, INTENT(IN ) :: itimestep
67	INTEGER, INTENT(IN ) :: num_force_layers, scm_force
68	REAL, INTENT(IN ) :: dt,dx
69	LOGICAL, INTENT(IN ) :: scm_th_adv, &
70	scm_qv_adv, &
71	scm_ql_adv, &
72	scm_wind_adv, &
73	scm_vert_adv, &
74	scm_soilT_force, &
75	scm_soilQ_force, &
76	scm_force_th_largescale, &
77	scm_force_qv_largescale, &
78	scm_force_ql_largescale, &
79	scm_force_wind_largescale
80
81	REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN ) :: z, th, qv, ql
82	REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN ) :: u, v
83	REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN ) :: z_at_w
84	REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: thten, qvten
85	REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: qlten
86	REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: uten, vten
87	REAL, DIMENSION( kms:kme ), INTENT(INOUT) :: u_base, v_base
88	REAL, DIMENSION( kms:kme ), INTENT(INOUT) :: z_base
89	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: z_force
90	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_g,v_g
91
92	REAL, DIMENSION(num_force_layers), INTENT (IN) :: z_force_tend
93	REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_g_tend,v_g_tend
94	REAL, DIMENSION(num_force_layers), INTENT (IN) :: w_subs_tend
95	REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_upstream_x_tend
96	REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_upstream_y_tend
97	REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_upstream_x_tend
98	REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_upstream_y_tend
99	REAL, DIMENSION(num_force_layers), INTENT (IN) :: ql_upstream_x_tend
100	REAL, DIMENSION(num_force_layers), INTENT (IN) :: ql_upstream_y_tend
101	REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_upstream_x_tend
102	REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_upstream_y_tend
103	REAL, DIMENSION(num_force_layers), INTENT (IN) :: v_upstream_x_tend
104	REAL, DIMENSION(num_force_layers), INTENT (IN) :: v_upstream_y_tend
105	REAL, DIMENSION(num_force_layers), INTENT (IN) :: tau_x_tend
106	REAL, DIMENSION(num_force_layers), INTENT (IN) :: tau_y_tend
107
108	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_upstream_x
109	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_upstream_y
110	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_upstream_x
111	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_upstream_y
112	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_upstream_x
113	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_upstream_y
114	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_upstream_x
115	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_upstream_y
116	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_upstream_x
117	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_upstream_y
118	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: w_subs
119	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_x
120	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_y
121
122	! WA 1/8/10 for large-scale forcing
123	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_largescale
124	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_largescale_tend
125	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_largescale
126	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_largescale_tend
127	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_largescale
128	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_largescale_tend
129	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_largescale
130	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_largescale_tend
131	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_largescale
132	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_largescale_tend
133	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_largescale
134	REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_largescale_tend
135
136	! WA 1/3/10 For soil forcing
137	INTEGER, INTENT(IN ) :: num_force_soil_layers, num_soil_layers
138	REAL, DIMENSION(ims:ime,num_soil_layers,jms:jme),INTENT(INOUT) :: tslb, smois
139	REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: t_soil_forcing_val
140	REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: t_soil_forcing_tend
141	REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: q_soil_forcing_val
142	REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: q_soil_forcing_tend
143	REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: tau_soil
144	REAL, DIMENSION(num_force_soil_layers), INTENT (IN ) :: soil_depth_force
145	REAL, DIMENSION(num_soil_layers), INTENT (IN ) :: zs
146
147	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
148	ims,ime, jms,jme, kms,kme, &
149	ips,ipe, jps,jpe, kps,kpe, &
150	kts,kte
151
152	! Local
153	INTEGER :: i,j,k
154	LOGICAL :: debug = .false.
155	REAL :: t_x, t_y, qv_x, qv_y, ql_x, ql_y
156	REAL :: u_x, u_y, v_x, v_y
157	REAL, DIMENSION(kms:kme) :: th_adv_tend, qv_adv_tend, ql_adv_tend
158	REAL, DIMENSION(kms:kme) :: u_adv_tend, v_adv_tend
159	REAL, DIMENSION(kms:kme) :: dthdz, dudz, dvdz, dqvdz, dqldz
160	REAL :: w
161	REAL, DIMENSION(kms:kme) :: w_dthdz, w_dudz, w_dvdz, w_dqvdz, w_dqldz
162	REAL, DIMENSION(kms:kme) :: adv_timescale_x, adv_timescale_y
163	CHARACTER*256 :: message
164	! Large-scale forcing WA 1/8/10
165	REAL :: t_ls, qv_ls, ql_ls
166	REAL :: u_ls, v_ls
167	REAL, DIMENSION(kms:kme) :: th_ls_tend, qv_ls_tend, ql_ls_tend
168	REAL, DIMENSION(kms:kme) :: u_ls_tend, v_ls_tend
169	REAL, DIMENSION(kms:kme) :: ls_timescale
170	! Soil forcing WA 1/3/10
171	INTEGER :: ks
172	REAL :: t_soil, q_soil
173	REAL, DIMENSION(num_soil_layers) :: t_soil_tend, q_soil_tend
174	REAL, DIMENSION(num_soil_layers) :: timescale_soil
175
176	IF ( scm_force .EQ. 0 ) return
177
178	! NOTES
179	! z is kts:kte
180	! z_at_w is kms:kme
181
182	! this is a good place for checks on the configuration
183	if ( z_force(1) > z(ids,1,jds) ) then
184	CALL wrf_message("First forcing level must be lower than first WRF half-level")
185	WRITE( message , * ) 'z forcing = ',z_force(1), 'z = ',z(ids,1,jds)
186	! print*,"z forcing = ",z_force(1), "z = ",z(ids,1,jds)
187	CALL wrf_error_fatal( message )
188	endif
189
190	z_force = z_force + dt*z_force_tend
191	u_g = u_g + dt*u_g_tend
192	v_g = v_g + dt*v_g_tend
193	tau_x = tau_x + dt*tau_x_tend
194	tau_y = tau_y + dt*tau_y_tend
195	tau_largescale = tau_largescale + dt*tau_largescale_tend
196
197	if ( scm_th_adv .AND. th_upstream_x(1) > 0.) then
198	th_upstream_x = th_upstream_x + dt*th_upstream_x_tend
199	th_upstream_y = th_upstream_y + dt*th_upstream_y_tend
200	endif
201	if ( scm_qv_adv .AND. qv_upstream_x(1) > 0.) then
202	qv_upstream_x = qv_upstream_x + dt*qv_upstream_x_tend
203	qv_upstream_y = qv_upstream_y + dt*qv_upstream_y_tend
204	endif
205	if ( scm_ql_adv .AND. ql_upstream_x(1) > 0.) then
206	ql_upstream_x = ql_upstream_x + dt*ql_upstream_x_tend
207	ql_upstream_y = ql_upstream_y + dt*ql_upstream_y_tend
208	endif
209	if ( scm_wind_adv .AND. u_upstream_x(1) > -900.) then
210	u_upstream_x = u_upstream_x + dt*u_upstream_x_tend
211	u_upstream_y = u_upstream_y + dt*u_upstream_y_tend
212	v_upstream_x = v_upstream_x + dt*v_upstream_x_tend
213	v_upstream_y = v_upstream_y + dt*v_upstream_y_tend
214	endif
215	if ( scm_vert_adv ) then
216	w_subs = w_subs + dt*w_subs_tend
217	endif
218
219	if ( scm_force_th_largescale .AND. th_largescale(1) > 0.) then
220	th_largescale = th_largescale + dt*th_largescale_tend
221	endif
222	if ( scm_force_qv_largescale .AND. qv_largescale(1) > 0.) then
223	qv_largescale = qv_largescale + dt*qv_largescale_tend
224	endif
225	if ( scm_force_ql_largescale.AND. ql_largescale(1) > 0.) then
226	ql_largescale = ql_largescale + dt*ql_largescale_tend
227	endif
228	if ( scm_force_wind_largescale .AND. u_largescale(1) > -900.) then
229	u_largescale = u_largescale + dt*u_largescale_tend
230	v_largescale = v_largescale + dt*v_largescale_tend
231	endif
232
233	if ( scm_soilT_force ) then
234	t_soil_forcing_val = t_soil_forcing_val + dt*t_soil_forcing_tend
235	endif
236	if ( scm_soilQ_force ) then
237	q_soil_forcing_val = q_soil_forcing_val + dt*q_soil_forcing_tend
238	endif
239
240	! 0 everything in case we don't set it later
241	th_adv_tend = 0.0
242	qv_adv_tend = 0.0
243	ql_adv_tend = 0.0
244	u_adv_tend = 0.0
245	v_adv_tend = 0.0
246	th_ls_tend = 0.0
247	qv_ls_tend = 0.0
248	ql_ls_tend = 0.0
249	u_ls_tend = 0.0
250	v_ls_tend = 0.0
251	w_dthdz = 0.0
252	w_dqvdz = 0.0
253	w_dqldz = 0.0
254	w_dudz = 0.0
255	w_dvdz = 0.0
256	adv_timescale_x = 0.0
257	adv_timescale_y = 0.0
258
259	! now interpolate forcing to model vertical grid
260
261	! if ( debug ) print*,' z u_base v_base '
262	CALL wrf_debug(100,'k z_base u_base v_base')
263	do k = kms,kme-1
264	z_base(k) = z(ids,k,jds)
265	u_base(k) = interp_0(u_g,z_force,z_base(k),num_force_layers)
266	v_base(k) = interp_0(v_g,z_force,z_base(k),num_force_layers)
267	! if ( debug ) print*,z_base(k),u_base(k),v_base(k)
268	WRITE( message, '(i4,3f12.4)' ) k,z_base(k),u_base(k),v_base(k)
269	CALL wrf_debug ( 100, message )
270	enddo
271
272	if ( scm_th_adv .or. scm_qv_adv .or. scm_ql_adv .or. scm_wind_adv ) then
273	if ( scm_th_adv ) CALL wrf_debug ( 100, 'k tau_x tau_y t_ups_x t_ups_y t_m ' )
274	do k = kms,kme-1
275	adv_timescale_x(k) = interp_0(tau_x,z_force,z(ids,k,jds),num_force_layers)
276	adv_timescale_y(k) = interp_0(tau_y,z_force,z(ids,k,jds),num_force_layers)
277	enddo
278	endif
279
280	if ( scm_th_adv ) then
281	if ( th_upstream_x(1) > 0.) then
282	do k = kms,kme-1
283	t_x = interp_0(th_upstream_x,z_force,z(ids,k,jds),num_force_layers)
284	t_y = interp_0(th_upstream_y,z_force,z(ids,k,jds),num_force_layers)
285
286	th_adv_tend(k) = (t_x-th(ids,k,jds))/adv_timescale_x(k) + (t_y-th(ids,k,jds))/adv_timescale_y(k)
287	WRITE( message, '(i4,5f12.4)' ) k,adv_timescale_x(k), adv_timescale_y(k), t_x, t_y, th(ids,k,jds)
288	CALL wrf_debug ( 100, message )
289	enddo
290	else ! WA if upstream is empty, use tendency only not value+tend
291	do k = kms,kme-1
292	t_x = interp_0(dt*th_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
293	t_y = interp_0(dt*th_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
294
295	th_adv_tend(k) = t_x/adv_timescale_x(k) + t_y/adv_timescale_y(k)
296	WRITE( message, '(i4,5f12.4)' ) k,adv_timescale_x(k), adv_timescale_y(k), t_x, t_y, th(ids,k,jds)
297	CALL wrf_debug ( 100, message )
298	enddo
299	endif
300	endif
301	if (minval(tau_x) < 0) then
302	print*,tau_x
303	stop 'TAU_X'
304	endif
305	if (minval(tau_y) < 0) then
306	print*,z_force
307	print*,tau_y
308	stop 'TAU_Y'
309	endif
310
311	if ( scm_qv_adv ) then
312	if ( qv_upstream_x(1) > 0.) then
313	do k = kms,kme-1
314	qv_x = interp_0(qv_upstream_x,z_force,z(ids,k,jds),num_force_layers)
315	qv_y = interp_0(qv_upstream_y,z_force,z(ids,k,jds),num_force_layers)
316
317	qv_adv_tend(k) = (qv_x-qv(ids,k,jds))/adv_timescale_x(k) + (qv_y-qv(ids,k,jds))/adv_timescale_y(k)
318	WRITE( message, * ) 'qv_adv_tend branch 1',k,adv_timescale_x(k), qv_upstream_x(k), adv_timescale_y(k), qv_x, qv_y, qv(ids,k,jds), qv_adv_tend(k)
319	CALL wrf_debug ( 100, message )
320	enddo
321	else ! WA if upstream is empty, use tendency only not value+tend
322	do k = kms,kme-1
323	qv_x = interp_0(dt*qv_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
324	qv_y = interp_0(dt*qv_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
325
326	qv_adv_tend(k) = qv_x/adv_timescale_x(k) + qv_y/adv_timescale_y(k)
327	WRITE( message, * ) 'qv_adv_tend branch 2',k,adv_timescale_x(k), adv_timescale_y(k), qv_upstream_x(k), qv_x, qv_y, qv(ids,k,jds), qv_adv_tend(k)
328	CALL wrf_debug ( 100, message )
329	enddo
330	endif
331	endif
332
333	if ( scm_ql_adv ) then
334	if ( ql_upstream_x(1) > 0.) then
335	do k = kms,kme-1
336	ql_x = interp_0(ql_upstream_x,z_force,z(ids,k,jds),num_force_layers)
337	ql_y = interp_0(ql_upstream_y,z_force,z(ids,k,jds),num_force_layers)
338
339	ql_adv_tend(k) = (ql_x-ql(ids,k,jds))/adv_timescale_x(k) + (ql_y-ql(ids,k,jds))/adv_timescale_y(k)
340	enddo
341	else ! WA if upstream is empty, use tendency only not value+tend
342	do k = kms,kme-1
343	ql_x = interp_0(dt*ql_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
344	ql_y = interp_0(dt*ql_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
345
346	ql_adv_tend(k) = ql_x/adv_timescale_x(k) + ql_y/adv_timescale_y(k)
347	enddo
348	endif
349	endif
350
351	if ( scm_wind_adv ) then
352	if ( u_upstream_x(1) > -900.) then
353	do k = kms,kme-1
354	u_x = interp_0(u_upstream_x,z_force,z(ids,k,jds),num_force_layers)
355	u_y = interp_0(u_upstream_y,z_force,z(ids,k,jds),num_force_layers)
356
357	v_x = interp_0(v_upstream_x,z_force,z(ids,k,jds),num_force_layers)
358	v_y = interp_0(v_upstream_y,z_force,z(ids,k,jds),num_force_layers)
359
360	u_adv_tend(k) = (u_x-u(ids,k,jds))/adv_timescale_x(k) + (u_y-u(ids,k,jds))/adv_timescale_y(k)
361	v_adv_tend(k) = (v_x-v(ids,k,jds))/adv_timescale_x(k) + (v_y-v(ids,k,jds))/adv_timescale_y(k)
362	enddo
363	else ! WA if upstream is empty, use tendency only not value+tend
364	do k = kms,kme-1
365	u_x = interp_0(dt*u_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
366	u_y = interp_0(dt*u_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
367
368	v_x = interp_0(dt*v_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
369	v_y = interp_0(dt*v_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
370
371	u_adv_tend(k) = u_x/adv_timescale_x(k) + u_y/adv_timescale_y(k)
372	v_adv_tend(k) = v_x/adv_timescale_x(k) + v_y/adv_timescale_y(k)
373	enddo
374	endif
375	endif
376
377	! Large scale forcing starts here 1/8/10 WA
378
379	if ( scm_force_th_largescale .or. scm_force_qv_largescale .or. scm_force_ql_largescale .or. scm_force_wind_largescale ) then
380	do k = kms,kme-1
381	ls_timescale(k) = interp_0(tau_largescale,z_force,z(ids,k,jds),num_force_layers)
382	enddo
383	endif
384
385	if ( scm_force_th_largescale ) then
386	if ( th_largescale(1) > 0.) then
387	do k = kms,kme-1
388	t_ls = interp_0(th_largescale,z_force,z(ids,k,jds),num_force_layers)
389	th_ls_tend(k) = (t_ls-th(ids,k,jds))/ls_timescale(k)
390	enddo
391	else ! WA if upstream is empty, use tendency only not value+tend
392	do k = kms,kme-1
393	t_ls = interp_0(dt*th_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
394	th_ls_tend(k) = t_ls/ls_timescale(k)
395	enddo
396	endif
397	endif
398
399	if ( scm_force_qv_largescale ) then
400	if ( qv_largescale(1) > 0.) then
401	do k = kms,kme-1
402	qv_ls = interp_0(qv_largescale,z_force,z(ids,k,jds),num_force_layers)
403	qv_ls_tend(k) = (qv_ls-qv(ids,k,jds))/ls_timescale(k)
404	enddo
405	else ! WA if upstream is empty, use tendency only not value+tend
406	do k = kms,kme-1
407	qv_ls = interp_0(dt*qv_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
408	qv_ls_tend(k) = qv_ls/ls_timescale(k)
409	enddo
410	endif
411	endif
412
413	if ( scm_force_ql_largescale ) then
414	if ( ql_largescale(1) > 0.) then
415	do k = kms,kme-1
416	ql_ls = interp_0(ql_largescale,z_force,z(ids,k,jds),num_force_layers)
417	ql_ls_tend(k) = (ql_ls-ql(ids,k,jds))/ls_timescale(k)
418	enddo
419	else ! WA if upstream is empty, use tendency only not value+tend
420	do k = kms,kme-1
421	ql_ls = interp_0(dt*ql_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
422	ql_ls_tend(k) = ql_ls/ls_timescale(k)
423	enddo
424	endif
425	endif
426
427	if ( scm_force_wind_largescale ) then
428	if ( u_largescale(1) > -900.) then
429	do k = kms,kme-1
430	u_ls = interp_0(u_largescale,z_force,z(ids,k,jds),num_force_layers)
431	v_ls = interp_0(v_largescale,z_force,z(ids,k,jds),num_force_layers)
432	u_ls_tend(k) = (u_ls-u(ids,k,jds))/ls_timescale(k)
433	v_ls_tend(k) = (v_ls-v(ids,k,jds))/ls_timescale(k)
434	enddo
435	else ! WA if upstream is empty, use tendency only not value+tend
436	do k = kms,kme-1
437	u_ls = interp_0(dt*u_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
438	v_ls = interp_0(dt*v_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
439	u_ls_tend(k) = u_ls/ls_timescale(k)
440	v_ls_tend(k) = v_ls/ls_timescale(k)
441	enddo
442	endif
443	endif
444
445	! Now do vertical advection. Note that no large-scale vertical advection
446	! is implemented at this time, may not make sense anyway (WA).
447	! loops are set so that the top and bottom (w=0) are handled correctly
448	! vertical derivatives
449	do k = kms+1,kme-1
450	dthdz(k) = (th(2,k,2)-th(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
451	dqvdz(k) = (qv(2,k,2)-qv(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
452	dqldz(k) = (ql(2,k,2)-ql(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
453	dudz(k) = (u(2,k,2)-u(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
454	dvdz(k) = (v(2,k,2)-v(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
455	enddo
456
457	! w on full levels, then advect
458	if ( scm_vert_adv ) then
459	do k = kms+1,kme-1
460	w = interp_0(w_subs,z_force,z_at_w(ids,k,jds),num_force_layers)
461	w_dthdz(k) = w*dthdz(k)
462	w_dqvdz(k) = w*dqvdz(k)
463	w_dqldz(k) = w*dqldz(k)
464	w_dudz(k) = w*dudz(k)
465	w_dvdz(k) = w*dvdz(k)
466	enddo
467	endif
468
469	! set tendencies for return
470	! vertical advection tendencies need to be interpolated back to half levels
471	CALL wrf_debug ( 100, 'j, k, th_adv_ten, qv_adv_ten, ql_adv_ten, u_adv_ten, v_adv_ten')
472	do j = jms,jme
473	do k = kms,kme-1
474	if(j==1) WRITE( message, * ) k,th_adv_tend(k),qv_adv_tend(k),ql_adv_tend(k), u_adv_tend(k),v_adv_tend(k)
475	if(j==1) CALL wrf_debug ( 100, message )
476	do i = ims,ime
477	thten(i,k,j) = thten(i,k,j) + th_adv_tend(k) + &
478	0.5*(w_dthdz(k) + w_dthdz(k+1)) &
479	+ th_ls_tend(k)
480	qvten(i,k,j) = qvten(i,k,j) + qv_adv_tend(k) + &
481	0.5*(w_dqvdz(k) + w_dqvdz(k+1)) &
482	+ qv_ls_tend(k)
483	qlten(i,k,j) = qlten(i,k,j) + ql_adv_tend(k) + &
484	0.5*(w_dqldz(k) + w_dqldz(k+1)) &
485	+ ql_ls_tend(k)
486	uten(i,k,j) = uten(i,k,j) + u_adv_tend(k) + &
487	0.5*(w_dudz(k) + w_dudz(k+1)) &
488	+ u_ls_tend(k)
489	vten(i,k,j) = vten(i,k,j) + v_adv_tend(k) + &
490	0.5*(w_dvdz(k) + w_dvdz(k+1)) &
491	+ v_ls_tend(k)
492	enddo
493	enddo
494	enddo
495
496	! soil forcing 1/3/10 WA
497	if ( scm_soilT_force ) then
498	do ks = 1,num_soil_layers
499	t_soil = interp_0(t_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
500	timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
501	t_soil_tend(ks) = (t_soil-tslb(ids,ks,jds))/timescale_soil(ks)
502	enddo
503	do j = jms,jme
504	do ks = 1,num_soil_layers
505	do i = ims,ime
506	tslb(ids,ks,jds) = tslb(ids,ks,jds) + t_soil_tend(ks)
507	enddo
508	enddo
509	enddo
510	endif
511	if ( scm_soilQ_force ) then
512	do ks = 1,num_soil_layers
513	q_soil = interp_0(q_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
514	timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
515	q_soil_tend(ks) = (q_soil-smois(ids,ks,jds))/timescale_soil(ks)
516	enddo
517	do j = jms,jme
518	do ks = 1,num_soil_layers
519	do i = ims,ime
520	smois(ids,ks,jds) = smois(ids,ks,jds) + q_soil_tend(ks)
521	enddo
522	enddo
523	enddo
524	endif
525
526	RETURN
527
528	END SUBROUTINE force_scm
529
530	END MODULE module_force_scm

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format