Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_fdda_psufddagd.F @ 1509

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 48.1 KB

Line
1	!wrf:model_layer:physics
2	!
3	!
4	!
5	MODULE module_fdda_psufddagd
6
7	CONTAINS
8	!
9	!-------------------------------------------------------------------
10	!
11	SUBROUTINE fddagd(itimestep,dx,dt,xtime, &
12	id,analysis_interval, end_fdda_hour, &
13	if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
14	if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
15	guv, gt, gq, if_ramping, dtramp_min, &
16
17	grid_sfdda, &
18	analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
19	rinblw, &
20
21	u3d,v3d,th3d,t3d, &
22	qv3d, &
23	p3d,pi3d, &
24	u_ndg_old,v_ndg_old,t_ndg_old,q_ndg_old,mu_ndg_old, &
25	u_ndg_new,v_ndg_new,t_ndg_new,q_ndg_new,mu_ndg_new, &
26	u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
27	rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
28	u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
29	rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
30	RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
31	pblh, ht, regime, znt, z, z_at_w, &
32	ids,ide, jds,jde, kds,kde, &
33	ims,ime, jms,jme, kms,kme, &
34	its,ite, jts,jte, kts,kte )
35
36	!-------------------------------------------------------------------
37	implicit none
38	!-------------------------------------------------------------------
39	!
40	! This code was implemented by Aijun Deng (Penn State). The 3-D analysis nudging was comleted
41	! and released in December 2006. The surface analysis nudging capability was added and
42	! released in March 2009 with WRFV3.1.
43	!
44	!-- u3d 3d u-velocity staggered on u points
45	!-- v3d 3d v-velocity staggered on v points
46	!-- th3d 3d potential temperature (k)
47	!-- t3d temperature (k)
48	!-- qv3d 3d water vapor mixing ratio (kg/kg)
49	!-- p3d 3d pressure (pa)
50	!-- pi3d 3d exner function (dimensionless)
51	!-- rundgdten staggered u tendency due to
52	! fdda grid nudging (m/s/s)
53	!-- rvndgdten staggered v tendency due to
54	! fdda grid nudging (m/s/s)
55	!-- rthndgdten theta tendency due to
56	! fdda grid nudging (K/s)
57	!-- rqvndgdten qv tendency due to
58	! fdda grid nudging (kg/kg/s)
59	!-- rmundgdten mu tendency due to
60	! fdda grid nudging (Pa/s)
61	!-- ids start index for i in domain
62	!-- ide end index for i in domain
63	!-- jds start index for j in domain
64	!-- jde end index for j in domain
65	!-- kds start index for k in domain
66	!-- kde end index for k in domain
67	!-- ims start index for i in memory
68	!-- ime end index for i in memory
69	!-- jms start index for j in memory
70	!-- jme end index for j in memory
71	!-- kms start index for k in memory
72	!-- kme end index for k in memory
73	!-- its start index for i in tile
74	!-- ite end index for i in tile
75	!-- jts start index for j in tile
76	!-- jte end index for j in tile
77	!-- kts start index for k in tile
78	!-- kte end index for k in tile
79	!-------------------------------------------------------------------
80	!
81	INTEGER, INTENT(IN) :: itimestep, analysis_interval, end_fdda_hour
82	INTEGER, INTENT(IN) :: analysis_interval_sfc, end_fdda_hour_sfc
83	INTEGER, INTENT(IN) :: grid_sfdda
84
85	INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
86	if_no_pbl_nudging_q
87	INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_q
88	INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_q
89	INTEGER, INTENT(IN) :: if_ramping
90
91	INTEGER , INTENT(IN) :: id
92	REAL, INTENT(IN) :: DT, dx, xtime, dtramp_min
93
94	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
95	ims,ime, jms,jme, kms,kme, &
96	its,ite, jts,jte, kts,kte
97
98	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
99	INTENT(IN) :: qv3d, &
100	p3d, &
101	pi3d, &
102	th3d, &
103	t3d, &
104	z, &
105	z_at_w
106
107	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
108	INTENT(INOUT) :: rundgdten, &
109	rvndgdten, &
110	rthndgdten, &
111	rqvndgdten
112
113	REAL, DIMENSION( ims:ime, jms:jme ), &
114	INTENT(INOUT) :: rmundgdten
115
116	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
117	INTENT(IN) :: u_ndg_old, &
118	v_ndg_old, &
119	t_ndg_old, &
120	q_ndg_old, &
121	u_ndg_new, &
122	v_ndg_new, &
123	t_ndg_new, &
124	q_ndg_new
125
126	REAL, DIMENSION( ims:ime, jms:jme ), &
127	INTENT(IN) :: u10_ndg_old, &
128	v10_ndg_old, &
129	t2_ndg_old, &
130	th2_ndg_old, &
131	q2_ndg_old, &
132	rh_ndg_old, &
133	psl_ndg_old, &
134	ps_ndg_old, &
135	u10_ndg_new, &
136	v10_ndg_new, &
137	t2_ndg_new, &
138	th2_ndg_new, &
139	q2_ndg_new, &
140	rh_ndg_new, &
141	psl_ndg_new, &
142	ps_ndg_new
143
144	REAL, DIMENSION( ims:ime, jms:jme ), &
145	INTENT(IN) :: tob_ndg_old, &
146	tob_ndg_new
147
148	REAL, DIMENSION( ims:ime, jms:jme ), &
149	INTENT(INOUT) :: mu_ndg_old, &
150	mu_ndg_new
151
152	REAL, DIMENSION( ims:ime, jms:jme ), &
153	INTENT(IN) :: odis_ndg_old, odis_ndg_new
154
155	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
156	INTENT(IN) :: u3d, &
157	v3d
158
159	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
160	ht, &
161	regime, &
162	znt
163
164	REAL, INTENT(IN) :: guv, gt, gq
165	REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc, rinblw
166
167	INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
168	REAL :: xtime_old, xtime_new, coef, val_analysis
169	INTEGER :: kpbl, dbg_level
170
171	REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
172	REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ) :: wpbl ! 1: u, 2: v, 3: t, 4: q
173	REAL, DIMENSION( kts:kte, 4 ) :: wzfac ! 1: u, 2: v, 3: t, 4: q
174
175	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
176
177	CHARACTER (LEN=256) :: message
178	INTEGER :: int4
179
180	int4 = 1 ! 1: temporal interpolation. else: target nudging toward *_ndg_new values
181
182	actual_end_fdda_min = end_fdda_hour*60.0
183	IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
184	actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
185	IF( xtime > actual_end_fdda_min ) THEN
186	! If xtime is greater than the end time, no need to calculate tendencies. Just set the tnedencies
187	! to zero to turn off nudging and return.
188	DO j = jts, jte
189	DO k = kts, kte
190	DO i = its, ite
191	RUNDGDTEN(i,k,j) = 0.0
192	RVNDGDTEN(i,k,j) = 0.0
193	RTHNDGDTEN(i,k,j) = 0.0
194	RQVNDGDTEN(i,k,j) = 0.0
195	IF( k .EQ. kts ) RMUNDGDTEN(i,j) = 0.
196	ENDDO
197	ENDDO
198	ENDDO
199	RETURN
200	ENDIF
201
202	IF( analysis_interval <= 0 )CALL wrf_error_fatal('In grid FDDA, gfdda_interval_m must be > 0')
203	xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
204	xtime_new = xtime_old + analysis_interval
205	IF( int4 == 1 ) THEN
206	coef = (xtime-xtime_old)/(xtime_new-xtime_old)
207	ELSE
208	coef = 1.0 ! Nudging toward a target value (*_ndg_new values)
209	ENDIF
210
211	IF ( wrf_dm_on_monitor()) THEN
212
213	CALL get_wrf_debug_level( dbg_level )
214
215	IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
216
217	IF( xtime < end_fdda_hour*60.0 ) THEN
218	WRITE(message,'(a,i1,a,f10.3,a)') &
219	'D0',id,' 3-D analysis nudging reads new data at time = ', xtime, ' min.'
220	CALL wrf_message( TRIM(message) )
221	WRITE(message,'(a,i1,a,2f8.2,a)') &
222	'D0',id,' 3-D analysis nudging bracketing times = ', xtime_old, xtime_new, ' min.'
223	CALL wrf_message( TRIM(message) )
224	ENDIF
225
226	actual_end_fdda_min = end_fdda_hour*60.0
227	IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
228	actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
229
230	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
231	! Find the mid point of the tile and print out the sample values
232	i0 = (ite-its)/2+its
233	j0 = (jte-jts)/2+jts
234
235	IF( guv > 0.0 ) THEN
236	DO k = kts, kte
237	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
238	' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
239	' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
240	CALL wrf_message( TRIM(message) )
241	ENDDO
242	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
243	' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
244	' mu_ndg_old=', mu_ndg_old(i0,j0), ' mu_ndg_new=', mu_ndg_new(i0,j0)
245	CALL wrf_message( TRIM(message) )
246	DO k = kts, kte
247	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
248	' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
249	' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
250	CALL wrf_message( TRIM(message) )
251	ENDDO
252	ENDIF
253
254	IF( gt > 0.0 ) THEN
255	DO k = kts, kte
256	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
257	' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
258	' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
259	CALL wrf_message( TRIM(message) )
260	ENDDO
261	ENDIF
262
263	IF( gq > 0.0 ) THEN
264	DO k = kts, kte
265	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
266	' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
267	' q_ndg_old=', q_ndg_old(i0,k,j0), ' q_ndg_new=', q_ndg_new(i0,k,j0)
268	CALL wrf_message( TRIM(message) )
269	ENDDO
270	ENDIF
271
272	IF( int4 == 1 ) then
273	WRITE(message,'(a,i1,a)') ' D0',id, &
274	' 3-D nudging towards the temporally interpolated analysis'
275	ELSE
276	WRITE(message,'(a,i1,a)') ' D0',id, &
277	' 3-D nudging towards the target analysis'
278	ENDIF
279
280	ENDIF
281	ENDIF
282	ENDIF
283
284	jtsv=MAX0(jts,jds+1)
285	itsu=MAX0(its,ids+1)
286
287	jtf=MIN0(jte,jde-1)
288	ktf=MIN0(kte,kde-1)
289	itf=MIN0(ite,ide-1)
290	!
291	! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t,
292	! if_no_pbl_nudging_q swithes) are set to 1, compute the weighting function, wpbl(:,k,:,:),
293	! based on the PBL depth. wpbl = 1 above the PBL and wpbl = 0 in the PBL. If all
294	! the switche are set to zero, wpbl = 1 (default value).
295	!
296	wpbl(:,:,:,:) = 1.0
297
298	IF( if_no_pbl_nudging_uv == 1 .OR. grid_sfdda == 1 ) THEN
299
300	DO j=jts,jtf
301	DO i=itsu,itf
302
303	kpbl = 1
304	zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
305
306	loop_ku: DO k=kts,ktf
307	! zagl = 0.5 * ( z(i-1,k,j)-ht(i-1,j) + z(i,k,j)-ht(i,j) )
308	zagl_bot = 0.5 * ( z_at_w(i-1,k, j)-ht(i-1,j) + z_at_w(i,k, j)-ht(i,j) )
309	zagl_top = 0.5 * ( z_at_w(i-1,k+1,j)-ht(i-1,j) + z_at_w(i,k+1,j)-ht(i,j) )
310	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
311	kpbl = k
312	EXIT loop_ku
313	ENDIF
314	ENDDO loop_ku
315
316	DO k=kts,ktf
317	IF( k <= kpbl ) wpbl(i, k, j, 1) = 0.0
318	IF( k == kpbl+1 ) wpbl(i, k, j, 1) = 0.1
319	IF( k > kpbl+1 ) wpbl(i, k, j, 1) = 1.0
320	ENDDO
321
322	ENDDO
323	ENDDO
324
325	DO i=its,itf
326	DO j=jtsv,jtf
327
328	kpbl = 1
329	zpbl = 0.5 * ( pblh(i,j-1) + pblh(i,j) )
330
331	loop_kv: DO k=kts,ktf
332	! zagl = 0.5 * ( z(i,k,j-1)-ht(i,j-1) + z(i,k,j)-ht(i,j) )
333	zagl_bot = 0.5 * ( z_at_w(i,k, j-1)-ht(i,j-1) + z_at_w(i,k, j)-ht(i,j) )
334	zagl_top = 0.5 * ( z_at_w(i,k+1,j-1)-ht(i,j-1) + z_at_w(i,k+1,j)-ht(i,j) )
335	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
336	kpbl = k
337	EXIT loop_kv
338	ENDIF
339	ENDDO loop_kv
340
341	DO k=kts,ktf
342	IF( k <= kpbl ) wpbl(i, k, j, 2) = 0.0
343	IF( k == kpbl+1 ) wpbl(i, k, j, 2) = 0.1
344	IF( k > kpbl+1 ) wpbl(i, k, j, 2) = 1.0
345	ENDDO
346
347	ENDDO
348	ENDDO
349
350	ENDIF
351
352	IF( if_no_pbl_nudging_t == 1 .OR. grid_sfdda == 1 ) THEN
353
354	DO j=jts,jtf
355	DO i=its,itf
356
357	kpbl = 1
358	zpbl = pblh(i,j)
359
360	loop_kt: DO k=kts,ktf
361	! zagl = z(i,k,j)-ht(i,j)
362	zagl_bot = z_at_w(i,k, j)-ht(i,j)
363	zagl_top = z_at_w(i,k+1,j)-ht(i,j)
364	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
365	kpbl = k
366	EXIT loop_kt
367	ENDIF
368	ENDDO loop_kt
369
370	DO k=kts,ktf
371	IF( k <= kpbl ) wpbl(i, k, j, 3) = 0.0
372	IF( k == kpbl+1 ) wpbl(i, k, j, 3) = 0.1
373	IF( k > kpbl+1 ) wpbl(i, k, j, 3) = 1.0
374	ENDDO
375
376	ENDDO
377	ENDDO
378
379	ENDIF
380
381	IF( if_no_pbl_nudging_q == 1 .OR. grid_sfdda == 1 ) THEN
382
383	DO j=jts,jtf
384	DO i=its,itf
385
386	kpbl = 1
387	zpbl = pblh(i,j)
388
389	loop_kq: DO k=kts,ktf
390	! zagl = z(i,k,j)-ht(i,j)
391	zagl_bot = z_at_w(i,k, j)-ht(i,j)
392	zagl_top = z_at_w(i,k+1,j)-ht(i,j)
393	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
394	kpbl = k
395	EXIT loop_kq
396	ENDIF
397	ENDDO loop_kq
398
399	DO k=kts,ktf
400	IF( k <= kpbl ) wpbl(i, k, j, 4) = 0.0
401	IF( k == kpbl+1 ) wpbl(i, k, j, 4) = 0.1
402	IF( k > kpbl+1 ) wpbl(i, k, j, 4) = 1.0
403	ENDDO
404
405	ENDDO
406	ENDDO
407
408	ENDIF
409	!
410	! If the user-defined namelist switches (if_zfac_uv, if_zfac_t,
411	! if_zfac_q swithes) are set to 1, compute the weighting function, wzfac(k,:),
412	! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_q) below which analysis
413	! nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x). If all
414	! the switche are set to zero, wzfac = 1 (default value).
415	!
416	wzfac(:,:) = 1.0
417
418	IF( if_zfac_uv == 1 ) THEN
419
420	DO j=jts,jtf
421	DO i=itsu,itf
422	DO k=kts,ktf
423	IF( k <= k_zfac_uv ) wzfac(k, 1:2) = 0.0
424	IF( k == k_zfac_uv+1 ) wzfac(k, 1:2) = 0.1
425	IF( k > k_zfac_uv+1 ) wzfac(k, 1:2) = 1.0
426	ENDDO
427	ENDDO
428	ENDDO
429
430	ENDIF
431
432	IF( if_zfac_t == 1 ) THEN
433
434	DO j=jts,jtf
435	DO i=itsu,itf
436	DO k=kts,ktf
437	IF( k <= k_zfac_t ) wzfac(k, 3) = 0.0
438	IF( k == k_zfac_t+1 ) wzfac(k, 3) = 0.1
439	IF( k > k_zfac_t+1 ) wzfac(k, 3) = 1.0
440	ENDDO
441	ENDDO
442	ENDDO
443
444	ENDIF
445
446	IF( if_zfac_q == 1 ) THEN
447
448	DO j=jts,jtf
449	DO i=itsu,itf
450	DO k=kts,ktf
451	IF( k <= k_zfac_q ) wzfac(k, 4) = 0.0
452	IF( k == k_zfac_q+1 ) wzfac(k, 4) = 0.1
453	IF( k > k_zfac_q+1 ) wzfac(k, 4) = 1.0
454	ENDDO
455	ENDDO
456	ENDDO
457
458	ENDIF
459	!
460	! If if_ramping and dtramp_min are defined by user, comput a time weighting function, tfac,
461	! for analysis nudging so that at the end of the nudging period (which has to be at a
462	! analysis time) we ramp down the nudging coefficient, based on the use-defined sign of dtramp_min.
463	!
464	! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at
465	! end_fdda_hour-ABS(dtramp_min).
466	!
467	! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at
468	! end_fdda_hour+ABS(dtramp_min). In this case, the obs values are extrapolated using
469	! the obs tendency saved from the previous FDDA wondow. More specifically for extrapolation,
470	! coef (see codes below) is recalculated to reflect extrapolation during the ramping period.
471	!
472	tfac = 1.0
473
474	IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
475
476	IF( dtramp_min <= 0.0 ) THEN
477	actual_end_fdda_min = end_fdda_hour*60.0
478	ELSE
479	actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
480	ENDIF
481
482	IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN
483	tfac = 1.0
484	ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
485	tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
486	IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/(analysis_interval*1.0)
487	ELSE
488	tfac = 0.0
489	ENDIF
490
491	ENDIF
492	!
493	! Surface Analysis Nudging
494	!
495	IF( grid_sfdda == 1 ) THEN
496	CALL SFDDAGD(itimestep,dx,dt,xtime, id, &
497	analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
498	rinblw, &
499	u3d,v3d,th3d,t3d, &
500	qv3d, &
501	p3d,pi3d, &
502	u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
503	rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
504	u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
505	rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
506	RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
507	pblh, ht, regime, znt, z, z_at_w, &
508	ids,ide, jds,jde, kds,kde, &
509	ims,ime, jms,jme, kms,kme, &
510	its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
511	actual_end_fdda_min, tfac )
512	ENDIF
513	!
514	! Compute 3-D nudging tendencies for u, v, t and q
515	!
516	DO j=jts,jtf
517	DO k=kts,ktf
518	DO i=itsu,itf
519	val_analysis = u_ndg_old(i,k,j) ( 1.0 - coef ) + u_ndg_new(i,k,j) coef
520	RUNDGDTEN(i,k,j) = RUNDGDTEN(i,k,j) + guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
521	( val_analysis - u3d(i,k,j) )
522	ENDDO
523	ENDDO
524	ENDDO
525
526	DO j=jtsv,jtf
527	DO k=kts,ktf
528	DO i=its,itf
529	val_analysis = v_ndg_old(i,k,j) ( 1.0 - coef ) + v_ndg_new(i,k,j) coef
530	RVNDGDTEN(i,k,j) = RVNDGDTEN(i,k,j) + guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
531	( val_analysis - v3d(i,k,j) )
532	ENDDO
533	ENDDO
534	ENDDO
535
536	DO j=jts,jtf
537	DO k=kts,ktf
538	DO i=its,itf
539	val_analysis = t_ndg_old(i,k,j) ( 1.0 - coef ) + t_ndg_new(i,k,j) coef
540	RTHNDGDTEN(i,k,j) = RTHNDGDTEN(i,k,j) + gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
541	( val_analysis - th3d(i,k,j) + 300.0 )
542
543	val_analysis = q_ndg_old(i,k,j) ( 1.0 - coef ) + q_ndg_new(i,k,j) coef
544	RQVNDGDTEN(i,k,j) = RQVNDGDTEN(i,k,j) + gq * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
545	( val_analysis - qv3d(i,k,j) )
546	ENDDO
547	ENDDO
548	ENDDO
549
550	END SUBROUTINE fddagd
551
552	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
553	SUBROUTINE sfddagd(itimestep,dx,dt,xtime, &
554	id, analysis_interval_sfc, end_fdda_hour_sfc, &
555	guv_sfc, gt_sfc, gq_sfc, rinblw, &
556	u3d,v3d,th3d,t3d, &
557	qv3d, &
558	p3d,pi3d, &
559	u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
560	rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
561	u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
562	rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
563	RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN, &
564	pblh, ht, regime, znt, z, z_at_w, &
565	ids,ide, jds,jde, kds,kde, &
566	ims,ime, jms,jme, kms,kme, &
567	its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
568	actual_end_fdda_min, tfac)
569
570	!-------------------------------------------------------------------
571	USE module_model_constants
572
573	implicit none
574
575	!-------------------------------------------------------------------
576	!
577	! This code was implemented by Aijun Deng (Penn State). The 3-D analysis nudging was comleted
578	! and released in December 2006. The surface analysis nudging capability was added and
579	! released in March 2009 with WRFV3.1.
580	!
581	!-- u3d 3d u-velocity staggered on u points
582	!-- v3d 3d v-velocity staggered on v points
583	!-- th3d 3d potential temperature (k)
584	!-- t3d temperature (k)
585	!-- qv3d 3d water vapor mixing ratio (kg/kg)
586	!-- p3d 3d pressure (pa)
587	!-- pi3d 3d exner function (dimensionless)
588	!-- rundgdten staggered u tendency due to
589	! fdda grid nudging (m/s/s)
590	!-- rvndgdten staggered v tendency due to
591	! fdda grid nudging (m/s/s)
592	!-- rthndgdten theta tendency due to
593	! fdda grid nudging (K/s)
594	!-- rqvndgdten qv tendency due to
595	! fdda grid nudging (kg/kg/s)
596	!-- rmundgdten mu tendency due to
597	! fdda grid nudging (Pa/s)
598	!-- ids start index for i in domain
599	!-- ide end index for i in domain
600	!-- jds start index for j in domain
601	!-- jde end index for j in domain
602	!-- kds start index for k in domain
603	!-- kde end index for k in domain
604	!-- ims start index for i in memory
605	!-- ime end index for i in memory
606	!-- jms start index for j in memory
607	!-- jme end index for j in memory
608	!-- kms start index for k in memory
609	!-- kme end index for k in memory
610	!-- its start index for i in tile
611	!-- ite end index for i in tile
612	!-- jts start index for j in tile
613	!-- jte end index for j in tile
614	!-- kts start index for k in tile
615	!-- kte end index for k in tile
616	!-------------------------------------------------------------------
617	!
618	INTEGER, INTENT(IN) :: itimestep, analysis_interval_sfc, end_fdda_hour_sfc
619
620	INTEGER , INTENT(IN) :: id
621	REAL, INTENT(IN) :: dx,DT, xtime
622
623	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
624	ims,ime, jms,jme, kms,kme, &
625	its,ite, jts,jte, kts,kte
626
627	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
628	INTENT(IN) :: qv3d, &
629	p3d, &
630	pi3d, &
631	th3d, &
632	t3d, &
633	z, &
634	z_at_w
635
636	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
637	INTENT(INOUT) :: rundgdten, &
638	rvndgdten, &
639	rthndgdten, &
640	rqvndgdten
641
642	REAL, DIMENSION( ims:ime, jms:jme ), &
643	INTENT(INOUT) :: rmundgdten
644
645	REAL, DIMENSION( ims:ime, jms:jme ), &
646	INTENT(IN) :: u10_ndg_old, &
647	v10_ndg_old, &
648	t2_ndg_old, &
649	th2_ndg_old, &
650	q2_ndg_old, &
651	rh_ndg_old, &
652	psl_ndg_old, &
653	ps_ndg_old, &
654	u10_ndg_new, &
655	v10_ndg_new, &
656	t2_ndg_new, &
657	th2_ndg_new, &
658	q2_ndg_new, &
659	rh_ndg_new, &
660	psl_ndg_new, &
661	ps_ndg_new
662
663	REAL, DIMENSION( ims:ime, jms:jme ), &
664	INTENT(IN) :: tob_ndg_old, &
665	tob_ndg_new
666
667	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
668	INTENT(IN) :: u3d, &
669	v3d
670
671	REAL, DIMENSION( ims:ime, jms:jme ), &
672	INTENT(IN) :: odis_ndg_old, odis_ndg_new
673
674	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
675	ht, &
676	regime, &
677	znt
678
679	REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc, rinblw
680
681	INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, j0
682	REAL :: xtime_old_sfc, xtime_new_sfc, coef, val_analysis, es
683	INTEGER :: kpbl, dbg_level
684
685	REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
686
687	REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ), &
688	INTENT(IN) :: wpbl ! 1: u, 2: v, 3: t, 4: q
689	REAL, DIMENSION( kts:kte, 4 ), &
690	INTENT(IN) :: wzfac ! 1: u, 2: v, 3: t, 4: q
691	REAL, DIMENSION( its:ite, jts:jte) :: wndcor_u, wndcor_v
692	REAL, DIMENSION( its-2:ite+1, jts-2:jte+1) :: blw_old, blw_new
693	REAL, DIMENSION( its:ite, kts:kte, jts:jte) :: qsat
694
695	REAL :: m, b=1.8, blw, rindx, x
696
697	REAL :: difz, wr14, wr1z, wr24, wr2z, wndfac, reg, znt0
698	INTEGER, INTENT(IN) :: if_ramping
699	REAL, INTENT(IN) :: dtramp_min
700
701	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
702
703	CHARACTER (LEN=256) :: message
704	INTEGER :: iwinds, idd, iqsat, int4
705
706	iwinds = 1 ! 1: Scale the surface wind analysis to the lowest model level,
707	! if the first model half-layer is greater than 10 meters
708	! and we are in the free convection regime (REGIME=4.0). else: no
709	idd = 1 ! 1: Obs data density correction is applied. else: no
710	iqsat = 1 ! 1: Remove super saturation. eles: no
711	int4 = 1 ! 1: temporal ionterpolation. else: target nudging toward *_ndg_new values
712
713	IF( analysis_interval_sfc <= 0 )CALL wrf_error_fatal('In grid sfc FDDA, sgfdda_interval_m must be > 0')
714	xtime_old_sfc = FLOOR(xtime/analysis_interval_sfc) * analysis_interval_sfc * 1.0
715	xtime_new_sfc = xtime_old_sfc + analysis_interval_sfc
716	IF( int4 == 1 ) THEN
717	coef = (xtime-xtime_old_sfc)/(xtime_new_sfc-xtime_old_sfc) ! Temporal interpolation
718	ELSE
719	coef = 1.0 ! Nudging toward a target value (*_ndg_new values)
720	ENDIF
721
722	IF ( wrf_dm_on_monitor()) THEN
723
724	CALL get_wrf_debug_level( dbg_level )
725
726	IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
727
728	IF( xtime < end_fdda_hour_sfc*60.0 ) THEN
729	WRITE(message,'(a,i1,a,f10.3,a)') &
730	'D0',id,' surface analysis nudging reads new data at time = ', xtime, ' min.'
731	CALL wrf_message( TRIM(message) )
732	WRITE(message,'(a,i1,a,2f8.2,a)') &
733	'D0',id,' surface analysis nudging bracketing times = ', xtime_old_sfc, xtime_new_sfc, ' min.'
734	CALL wrf_message( TRIM(message) )
735	ENDIF
736
737	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
738	! Find the mid point of the tile and print out the sample values
739	i0 = (ite-its)/2+its
740	j0 = (jte-jts)/2+jts
741
742	IF( guv_sfc > 0.0 ) THEN
743	WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
744	' D0',id,' sample surface analysis values at i,j=', i0, j0, &
745	' u10_ndg_old=', u10_ndg_old(i0,j0), ' u10_ndg_new=', u10_ndg_new(i0,j0)
746	CALL wrf_message( TRIM(message) )
747	WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
748	' D0',id,' sample surface analysis values at i,j=', i0, j0, &
749	' v10_ndg_old=', v10_ndg_old(i0,j0), ' v10_ndg_new=', v10_ndg_new(i0,j0)
750	CALL wrf_message( TRIM(message) )
751	ENDIF
752
753	IF( gt_sfc > 0.0 ) THEN
754	WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
755	' D0',id,' sample surface analysis values at i,j=', i0, j0, &
756	' th2_ndg_old=', th2_ndg_old(i0,j0), ' th2_ndg_new=', th2_ndg_new(i0,j0)
757	CALL wrf_message( TRIM(message) )
758	ENDIF
759
760	IF( gq_sfc > 0.0 ) THEN
761	WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
762	' D0',id,' sample surface analysis values at i,j=', i0, j0, &
763	' q2_ndg_old=', q2_ndg_old(i0,j0), ' q2_ndg_new=', q2_ndg_new(i0,j0)
764	CALL wrf_message( TRIM(message) )
765	ENDIF
766
767	IF( iwinds == 1 ) &
768	WRITE(message,'(a,i1,a)') ' D0',id, &
769	' surface wind analysis s scaled to the lowest model level, if dz1 > 10m and REGIME=4.'
770
771	IF( idd == 1 ) &
772	WRITE(message,'(a,i1,a)') ' D0',id, &
773	' obs data density is used for additional weighting function'
774
775	IF( iqsat == 1 ) &
776	WRITE(message,'(a,i1,a)') ' D0',id, &
777	' super saturation is not allowed for q analysis'
778
779	IF( int4 == 1 ) then
780	WRITE(message,'(a,i1,a)') ' D0',id, &
781	' surface nudging towards the temporally interpolated analysis'
782	ELSE
783	WRITE(message,'(a,i1,a)') ' D0',id, &
784	' surface nudging towards the target analysis'
785	ENDIF
786
787	ENDIF
788	ENDIF
789	ENDIF
790
791
792	jtsv=MAX0(jts,jds+1)
793	itsu=MAX0(its,ids+1)
794
795	jtf=MIN0(jte,jde-1)
796	ktf=MIN0(kte,kde-1)
797	itf=MIN0(ite,ide-1)
798
799	!
800	! Compute the vertical weighting function to scale the surface wind analysis to
801	! the lowest model level, if the first model half-layer is greater
802	! than 10 meters and we are in the free convection regime (REGIME=4.0).
803	!
804	IF( iwinds == 1 ) THEN
805	wndcor_u(:,:) = 1.0
806	DO j=jts,jtf
807	DO i=itsu,itf
808	reg = 0.5 * ( regime(i-1, j) + regime(i, j) )
809	difz = 0.5 * ( z(i-1,1,j) - ht(i-1,j) &
810	+ z(i, 1,j) - ht(i, j) )
811	IF( reg > 3.5 .AND. difz > 10.0 ) THEN
812	znt0 = 0.5 * ( znt(i-1, j) + znt(i, j) )
813	IF( znt0 <= 0.2) THEN
814	wndcor_u(i,j) = 1.0+0.320znt0*0.2
815	ELSE
816	wndcor_u(i,j) = 1.169+0.315*znt0
817	ENDIF
818
819	wr14 = log(40.0/0.05)
820	wr1z = log(difz/0.05)
821	wr24 = log(40.0/1.0)
822	wr2z = log(difz/1.0)
823	wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
824	wndcor_u(i,j) = wndfac*wndcor_u(i,j)
825	ENDIF
826	ENDDO
827	ENDDO
828
829	IF ( wrf_dm_on_monitor()) THEN
830	IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
831	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
832	i0 = (ite-its)/2+its
833	j0 = (jte-jts)/2+jts
834	WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
835	' D0',id,' sample wndcor_u values at i,j=', i0, j0, &
836	' wndcor_u=', wndcor_u(i0,j0)
837	CALL wrf_message( TRIM(message) )
838	ENDIF
839	ENDIF
840	ENDIF
841	ELSE
842	wndcor_u(:,:) = 1.0
843	ENDIF
844
845	IF( iwinds == 1 ) THEN
846	wndcor_v(:,:) = 1.0
847	DO j=jtsv,jtf
848	DO i=its,itf
849	reg = 0.5 * ( regime(i, j-1) + regime(i, j) )
850	difz = 0.5 * ( z(i,1,j-1) - ht(i,j-1) &
851	+ z(i,1,j ) - ht(i,j ) )
852	IF( reg > 3.5 .AND. difz > 10.0 ) THEN
853	znt0 = 0.5 * ( znt(i, j-1) + znt(i, j) )
854	IF( znt0 <= 0.2) THEN
855	wndcor_v(i,j) = 1.0+0.320znt0*0.2
856	ELSE
857	wndcor_v(i,j) = 1.169+0.315*znt0
858	ENDIF
859
860	wr14 = log(40.0/0.05)
861	wr1z = log(difz/0.05)
862	wr24 = log(40.0/1.0)
863	wr2z = log(difz/1.0)
864	wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
865	wndcor_v(i,j) = wndfac*wndcor_v(i,j)
866	ENDIF
867	ENDDO
868	ENDDO
869	IF ( wrf_dm_on_monitor()) THEN
870	IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
871	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
872	i0 = (ite-its)/2+its
873	j0 = (jte-jts)/2+jts
874	WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
875	' D0',id,' sample wndcor_v values at i,j=', i0, j0, &
876	' wndcor_v=', wndcor_v(i0,j0)
877	CALL wrf_message( TRIM(message) )
878	ENDIF
879	ENDIF
880	ENDIF
881	ELSE
882	wndcor_v(:,:) = 1.0
883	ENDIF
884	!
885	! Compute saturation mixing ratio so that nudging to a super-saturated state
886	! is not allowed.
887	!
888	IF( iqsat == 1 ) THEN
889	DO j=jts,jtf
890	DO k=kts,ktf
891	DO i=its,itf
892	es = SVP1EXP(SVP2(t3d(i,k,j)-SVPT0)/(t3d(i,k,j)-SVP3)) * 10.0 ! mb
893	qsat(i,k,j) = EP_2*es/(p3d(i,k,j)/100.0-es)
894	ENDDO
895	ENDDO
896	ENDDO
897
898	IF ( wrf_dm_on_monitor()) THEN
899	IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
900	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
901	i0 = (ite-its)/2+its
902	j0 = (jte-jts)/2+jts
903	DO k = kts, kte
904	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
905	' D0',id,' sample moisture values (g/kg) at i,k,j=', i0, k, j0, &
906	' qv3d=', qv3d(i0,k,j0)1000.0, ' qsat=', qsat(i0,k,j0)1000.0
907	CALL wrf_message( TRIM(message) )
908	ENDDO
909	ENDIF
910	ENDIF
911	ENDIF
912	ENDIF
913	!
914	! Obs data density weighting.
915	!
916	IF( idd == 1 ) THEN
917
918	IF( rinblw < 0.001 ) THEN
919	IF ( wrf_dm_on_monitor()) THEN
920	WRITE(message,'(a)') 'Error in rinblw, please specify a reasonable value ***'
921	CALL wrf_message( TRIM(message) )
922	ENDIF
923	CALL wrf_error_fatal('In grid FDDA')
924	ENDIF
925
926	rindx = rinblw*1000.0/dx
927	m = -0.8*2.0/rindx
928
929	DO j=MAX(jts-2,jds),MIN(jtf+1,jde-1)
930	DO i=MAX(its-2,ids),MIN(itf+1,ide-1)
931	IF( odis_ndg_old(i,j) < 0.5*rinblw ) THEN
932	blw_old(i,j) = 1.0
933	ELSE
934	x = min( odis_ndg_old(i,j)*1000./dx, rindx )
935	blw_old(i,j) = m * x + b
936	ENDIF
937
938	IF( odis_ndg_new(i,j) < 0.5*rinblw ) THEN
939	blw_new(i,j) = 1.0
940	ELSE
941	x = min( odis_ndg_new(i,j)*1000./dx, rindx )
942	blw_new(i,j) = m * x + b
943	ENDIF
944	ENDDO
945	ENDDO
946
947	! Smoother applies one point outside the tile, but one point in from boundaries
948	CALL smther(blw_old, its-2,itf+1, jts-2,jtf+1, 1, &
949	MAX(its-2,ids+1), MIN(ite+1,ide-2), MAX(jts-2,jds+1), MIN(jte+1,jde-2))
950	CALL smther(blw_new, its-2,itf+1, jts-2,jtf+1, 1, &
951	MAX(its-2,ids+1), MIN(ite+1,ide-2), MAX(jts-2,jds+1), MIN(jte+1,jde-2))
952
953	WHERE ( blw_old > 1.0)
954	blw_old = 1.0
955	END WHERE
956	WHERE ( blw_new > 1.0)
957	blw_new = 1.0
958	END WHERE
959	WHERE ( blw_old < 0.0)
960	blw_old = 0.0
961	END WHERE
962	WHERE ( blw_new < 0.0)
963	blw_new = 0.0
964	END WHERE
965
966	IF ( wrf_dm_on_monitor()) THEN
967	IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
968	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
969	i0 = (ite-its)/2+its
970	j0 = (jte-jts)/2+jts
971	WRITE(message,'(a,i1,a,2i4,4(a,f10.4))') &
972	' D0',id,' sample blw values at i,j=', i0, j0, &
973	' odis_ndg_old=', odis_ndg_old(i0,j0), ' km odis_ndg_new=', odis_ndg_new(i0,j0), &
974	' km blw_old=', blw_old(i0,j0), ' blw_new=', blw_new(i0,j0)
975	CALL wrf_message( TRIM(message) )
976	ENDIF
977	ENDIF
978	ENDIF
979
980	ENDIF
981	!
982	! TFAC for surface analysis nudging
983	!
984	IF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min &
985	.AND. dtramp_min > 0.0 .AND. if_ramping == 1 ) &
986	coef = (xtime-xtime_old_sfc+analysis_interval_sfc)/(analysis_interval_sfc*1.0)
987	! print*, 'coef =', xtime_old_sfc, xtime, xtime_new_sfc, coef
988	!
989	! Compute surface analysis nudging tendencies for u, v, t and q
990	!
991	DO j=jts,jtf
992	DO k=kts,ktf
993	DO i=itsu,itf
994	IF( idd == 1 ) THEN
995	blw = 0.5* (blw_old(i-1,j)+blw_old(i,j)) * ( 1.0 - coef ) &
996	+ 0.5* (blw_new(i-1,j)+blw_new(i,j)) * coef
997	ELSE
998	blw = 1.0
999	ENDIF
1000	val_analysis = u10_ndg_old(i,j) ( 1.0 - coef ) + u10_ndg_new(i,j) coef
1001	val_analysis = val_analysis * wndcor_u(i,j)
1002	RUNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,1)) * wzfac(k,1) * tfac * blw * &
1003	( val_analysis - u3d(i,1,j) )
1004	ENDDO
1005	ENDDO
1006	ENDDO
1007
1008	DO j=jtsv,jtf
1009	DO k=kts,ktf
1010	DO i=its,itf
1011	IF( idd == 1 ) THEN
1012	blw = 0.5* (blw_old(i,j-1)+blw_old(i,j)) * ( 1.0 - coef ) &
1013	+ 0.5* (blw_new(i,j-1)+blw_new(i,j)) * coef
1014	ELSE
1015	blw = 1.0
1016	ENDIF
1017	val_analysis = v10_ndg_old(i,j) ( 1.0 - coef ) + v10_ndg_new(i,j) coef
1018	val_analysis = val_analysis * wndcor_v(i,j)
1019	RVNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,2)) * wzfac(k,2) * tfac * blw * &
1020	( val_analysis - v3d(i,1,j) )
1021	ENDDO
1022	ENDDO
1023	ENDDO
1024
1025	DO j=jts,jtf
1026	DO k=kts,ktf
1027	DO i=its,itf
1028	IF( idd == 1 ) THEN
1029	blw = blw_old(i,j) * ( 1.0 - coef ) + blw_new(i,j) * coef
1030	ELSE
1031	blw = 1.0
1032	ENDIF
1033	val_analysis = th2_ndg_old(i,j) ( 1.0 - coef ) + th2_ndg_new(i,j) coef
1034	RTHNDGDTEN(i,k,j) = gt_sfc * (1.0-wpbl(i,k,j,3)) * wzfac(k,3) * tfac * blw * &
1035	( val_analysis - th3d(i,1,j))
1036
1037	val_analysis = q2_ndg_old(i,j) ( 1.0 - coef ) + q2_ndg_new(i,j) coef
1038	IF( iqsat == 1 .AND. val_analysis > qsat(i,k,j) ) val_analysis = qsat(i,k,j)
1039	RQVNDGDTEN(i,k,j) = gq_sfc * (1.0-wpbl(i,k,j,4)) * wzfac(k,4) * tfac * blw * &
1040	( val_analysis - qv3d(i,k,j) )
1041	ENDDO
1042	ENDDO
1043	ENDDO
1044
1045	END SUBROUTINE sfddagd
1046
1047	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1048
1049	SUBROUTINE fddagdinit(id,rundgdten,rvndgdten,rthndgdten,rqvndgdten,rmundgdten,&
1050	run_hours, &
1051	if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
1052	if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
1053	guv, gt, gq, if_ramping, dtramp_min, end_fdda_hour, &
1054	grid_sfdda, guv_sfc, gt_sfc, gq_sfc, &
1055	restart, allowed_to_read, &
1056	ids, ide, jds, jde, kds, kde, &
1057	ims, ime, jms, jme, kms, kme, &
1058	its, ite, jts, jte, kts, kte )
1059	!-------------------------------------------------------------------
1060	IMPLICIT NONE
1061	!-------------------------------------------------------------------
1062	!
1063	INTEGER , INTENT(IN) :: id
1064	LOGICAL, INTENT(IN) :: restart, allowed_to_read
1065	INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1066	ims, ime, jms, jme, kms, kme, &
1067	its, ite, jts, jte, kts, kte
1068	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
1069	rundgdten, &
1070	rvndgdten, &
1071	rthndgdten, &
1072	rqvndgdten
1073	INTEGER, INTENT(IN) :: run_hours
1074	INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
1075	if_no_pbl_nudging_q, end_fdda_hour
1076	INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_q
1077	INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_q
1078	INTEGER, INTENT(IN) :: if_ramping, grid_sfdda
1079	REAL, INTENT(IN) :: dtramp_min
1080	REAL, INTENT(IN) :: guv, gt, gq
1081	REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc
1082	REAL :: actual_end_fdda_min
1083
1084	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: rmundgdten
1085	INTEGER :: i, j, k
1086
1087	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
1088
1089	CHARACTER (LEN=256) :: message
1090
1091	IF ( wrf_dm_on_monitor() ) THEN
1092
1093	IF( guv > 0.0 ) THEN
1094	WRITE(message,'(a,i1,a,e12.4)') &
1095	'D0',id,' 3-D analysis nudging for wind is applied and Guv= ', guv
1096	CALL wrf_message(TRIM(message))
1097	ELSE IF( guv < 0.0 ) THEN
1098	CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
1099	ELSE
1100	WRITE(message,'(a,i1,a,e12.4)') &
1101	'D0',id,' 3-D analysis nudging for wind is not applied and Guv= ', guv
1102	CALL wrf_message(TRIM(message))
1103	ENDIF
1104
1105	IF( gt > 0.0 ) THEN
1106	WRITE(message,'(a,i1,a,e12.4)') &
1107	'D0',id,' 3-D analysis nudging for temperature is applied and Gt= ', gt
1108	CALL wrf_message(TRIM(message))
1109	ELSE IF( gt < 0.0 ) THEN
1110	CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
1111	ELSE
1112	WRITE(message,'(a,i1,a,e12.4)') &
1113	'D0',id,' 3-D analysis nudging for temperature is not applied and Gt= ', gt
1114	CALL wrf_message(TRIM(message))
1115	ENDIF
1116
1117	IF( gq > 0.0 ) THEN
1118	WRITE(message,'(a,i1,a,e12.4)') &
1119	'D0',id,' 3-D analysis nudging for water vapor mixing ratio is applied and Gq= ', gq
1120	CALL wrf_message(TRIM(message))
1121	ELSE IF( gq < 0.0 ) THEN
1122	CALL wrf_error_fatal('In grid FDDA, Gq must be positive.')
1123	ELSE
1124	WRITE(message,'(a,i1,a,e12.4)') &
1125	'D0',id,' 3-D analysis nudging for water vapor mixing ratio is not applied and Gq= ', gq
1126	CALL wrf_message(TRIM(message))
1127	ENDIF
1128
1129	IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
1130	WRITE(message,'(a,i1,a)') &
1131	'D0',id,' 3-D analysis nudging for wind is turned off within the PBL.'
1132	CALL wrf_message(TRIM(message))
1133	ENDIF
1134
1135	IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
1136	WRITE(message,'(a,i1,a)') &
1137	'D0',id,' 3-D analysis nudging for temperature is turned off within the PBL.'
1138	CALL wrf_message(TRIM(message))
1139	ENDIF
1140
1141	IF( gq > 0.0 .AND. if_no_pbl_nudging_q == 1 ) THEN
1142	WRITE(message,'(a,i1,a)') &
1143	'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off within the PBL.'
1144	CALL wrf_message(TRIM(message))
1145	ENDIF
1146
1147	IF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
1148	WRITE(message,'(a,i1,a,i3)') &
1149	'D0',id,' 3-D analysis nudging for wind is turned off below layer', k_zfac_uv
1150	CALL wrf_message(TRIM(message))
1151	ENDIF
1152
1153	IF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
1154	WRITE(message,'(a,i1,a,i3)') &
1155	'D0',id,' 3-D analysis nudging for temperature is turned off below layer', k_zfac_t
1156	CALL wrf_message(TRIM(message))
1157	ENDIF
1158
1159	IF( gq > 0.0 .AND. if_zfac_q == 1 ) THEN
1160	WRITE(message,'(a,i1,a,i3)') &
1161	'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off below layer', &
1162	k_zfac_q
1163	CALL wrf_message(TRIM(message))
1164	ENDIF
1165
1166	IF( grid_sfdda ==1 ) THEN
1167	IF( guv_sfc > 0.0 ) THEN
1168	WRITE(message,'(a,i1,a,e12.4)') &
1169	'D0',id,' surface analysis nudging for wind is applied and Guv_sfc= ', guv_sfc
1170	CALL wrf_message(TRIM(message))
1171	ELSE IF( guv_sfc < 0.0 ) THEN
1172	CALL wrf_error_fatal('In grid FDDA, Guv_sfc must be positive.')
1173	ELSE
1174	WRITE(message,'(a,i1,a,e12.4)') &
1175	'D0',id,' surface analysis nudging for wind is not applied and Guv_sfc= ', guv_sfc
1176	CALL wrf_message(TRIM(message))
1177	ENDIF
1178
1179	IF( gt_sfc > 0.0 ) THEN
1180	WRITE(message,'(a,i1,a,e12.4)') &
1181	'D0',id,' surface analysis nudging for temperature is applied and Gt_sfc= ', gt_sfc
1182	CALL wrf_message(TRIM(message))
1183	ELSE IF( gt_sfc < 0.0 ) THEN
1184	CALL wrf_error_fatal('In grid FDDA, Gt_sfc must be positive.')
1185	ELSE
1186	WRITE(message,'(a,i1,a,e12.4)') &
1187	'D0',id,' surafce analysis nudging for temperature is not applied and Gt_sfc= ', gt_sfc
1188	CALL wrf_message(TRIM(message))
1189	ENDIF
1190
1191	IF( gq_sfc > 0.0 ) THEN
1192	WRITE(message,'(a,i1,a,e12.4)') &
1193	'D0',id,' surface analysis nudging for water vapor mixing ratio is applied and Gq_sfc= ', gq_sfc
1194	CALL wrf_message(TRIM(message))
1195	ELSE IF( gq_sfc < 0.0 ) THEN
1196	CALL wrf_error_fatal('In grid FDDA, Gq_sfc must be positive.')
1197	ELSE
1198	WRITE(message,'(a,i1,a,e12.4)') &
1199	'D0',id,' surface analysis nudging for water vapor mixing ratio is not applied and Gq_sfc= ', gq_sfc
1200	CALL wrf_message(TRIM(message))
1201	ENDIF
1202
1203	ENDIF
1204
1205	IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
1206	IF( dtramp_min <= 0.0 ) THEN
1207	actual_end_fdda_min = end_fdda_hour*60.0
1208	ELSE
1209	actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
1210	ENDIF
1211
1212	IF( actual_end_fdda_min <= run_hours*60. ) THEN
1213	WRITE(message,'(a,i1,a)') &
1214	'D0',id,' analysis nudging is ramped down near the end of the nudging period,'
1215	CALL wrf_message(TRIM(message))
1216
1217	WRITE(message,'(a,f6.2,a,f6.2,a)') &
1218	' starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0, &
1219	'h, ending at ', actual_end_fdda_min/60.0,'h.'
1220	CALL wrf_message(TRIM(message))
1221	ENDIF
1222	ENDIF
1223
1224	ENDIF
1225
1226	IF(.not.restart) THEN
1227	DO j = jts,jte
1228	DO k = kts,kte
1229	DO i = its,ite
1230	rundgdten(i,k,j) = 0.
1231	rvndgdten(i,k,j) = 0.
1232	rthndgdten(i,k,j) = 0.
1233	rqvndgdten(i,k,j) = 0.
1234	if(k.eq.kts) rmundgdten(i,j) = 0.
1235	ENDDO
1236	ENDDO
1237	ENDDO
1238	ENDIF
1239
1240	END SUBROUTINE fddagdinit
1241
1242	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1243	SUBROUTINE smther(slab, idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd)
1244	!
1245	! PURPOSE: SPATIALLY SMOOTH DATA IN SLAB TO DAMPEN SHORT
1246	! WAVELENGTH COMPONENTS.
1247	!
1248	! Implemented based on the same smoothing subroutine used in MM5, with modifications to
1249	! remove the extra loop that causes unneccessary desmoothing. Aijun Deng (Penn State),
1250	! December 2008
1251	!
1252	IMPLICIT NONE
1253
1254	INTEGER :: idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd
1255	INTEGER :: i, j, k, kp
1256	REAL :: asv, aplus, cell
1257	REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB
1258	REAL, DIMENSION(2) :: XNU
1259
1260	IF(NPASS.EQ.0)RETURN
1261	XNU(1)=0.50
1262	XNU(2)=-0.52
1263	DO K=1,NPASS
1264	! FIRST, SMOOTH IN THE J DIRECTION
1265	DO J=JST,JND
1266	ASV=SLAB(IST-1,J)
1267	DO I=IST,IND
1268	APLUS=SLAB(I+1,J)
1269	CELL=SLAB(I,J)
1270	SLAB(I,J)=SLAB(I,J)+XNU(K)*((ASV+APLUS)/2.0-SLAB(I,J))
1271	ASV=CELL
1272	ENDDO
1273	ENDDO
1274
1275	! NOW, SMOOTH IN THE I DIRECTION
1276	DO I=IST,IND
1277	ASV=SLAB(I,JST-1)
1278	DO J=JST,JND
1279	APLUS=SLAB(I,J+1)
1280	CELL=SLAB(I,J)
1281	SLAB(I,J)=SLAB(I,J)+XNU(K)*((ASV+APLUS)/2.0-SLAB(I,J))
1282	ASV=CELL
1283	ENDDO
1284	ENDDO
1285
1286	ENDDO
1287
1288	END SUBROUTINE smther
1289
1290	!-------------------------------------------------------------------
1291	END MODULE module_fdda_psufddagd

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

Download in other formats:

Original Format