Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_fdda_spnudging.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: 35.8 KB

Line
1	!wrf:model_layer:physics
2	!
3	! Code contributed by Gonzalo Miguez-Macho, Univ. de Santiago de Compostela, Galicia, Spain
4	!
5	MODULE module_fdda_spnudging
6
7	#ifdef DM_PARALLEL
8	USE module_dm , ONLY : ntasks_x, ntasks_y, local_communicator_x, local_communicator_y, data_order_xzy
9	#endif
10	USE module_wrf_error , ONLY : wrf_err_message
11
12	CONTAINS
13	!
14	!-------------------------------------------------------------------
15	!
16	SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fdda_hour, &
17	if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph,&
18	if_zfac_uv, k_zfac_uv, dk_zfac_uv, &
19	if_zfac_t, k_zfac_t, dk_zfac_t, &
20	if_zfac_ph, k_zfac_ph, dk_zfac_ph, &
21	guv, gt, gph, if_ramping, dtramp_min, xwavenum, ywavenum, &
22	u3d,v3d,th3d,ph3d, &
23	u_ndg_old,v_ndg_old,t_ndg_old,ph_ndg_old, &
24	u_ndg_new,v_ndg_new,t_ndg_new,ph_ndg_new, &
25	RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RPHNDGDTEN,&
26	pblh, ht, z, z_at_w, &
27	ids,ide, jds,jde, kds,kde, &
28	ims,ime, jms,jme, kms,kme, &
29	i_start,i_end, j_start,j_end, kts,kte, num_tiles, &
30	ips,ipe,jps,jpe,kps,kpe, &
31	imsx,imex,jmsx,jmex,kmsx,kmex, &
32	ipsx,ipex,jpsx,jpex,kpsx,kpex, &
33	imsy,imey,jmsy,jmey,kmsy,kmey, &
34	ipsy,ipey,jpsy,jpey,kpsy,kpey )
35
36
37
38	USE module_state_description
39	USE module_domain, ONLY : domain
40
41	!-------------------------------------------------------------------
42	implicit none
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	!---ph3d 3d perturbation geopotential
48	!-- rundgdten staggered u tendency due to
49	! spectral nudging (m/s/s)
50	!-- rvndgdten staggered v tendency due to
51	! spectral nudging (m/s/s)
52	!-- rthndgdten theta tendency due to
53	! spectral nudging (K/s)
54	!-- phndgdten ph tendency due to
55	! spectral nudging (m2/s2/s)
56	!-- ids start index for i in domain
57	!-- ide end index for i in domain
58	!-- jds start index for j in domain
59	!-- jde end index for j in domain
60	!-- kds start index for k in domain
61	!-- kde end index for k in domain
62	!-- ims start index for i in memory
63	!-- ime end index for i in memory
64	!-- jms start index for j in memory
65	!-- jme end index for j in memory
66	!-- kms start index for k in memory
67	!-- kme end index for k in memory
68	!-- its start index for i in tile
69	!-- ite end index for i in tile
70	!-- jts start index for j in tile
71	!-- jte end index for j in tile
72	!-- kts start index for k in tile
73	!-- kte end index for k in tile
74	!-------------------------------------------------------------------
75	!
76	TYPE(domain) , TARGET :: grid
77
78	INTEGER, INTENT(IN) :: itimestep, analysis_interval, end_fdda_hour
79
80	INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
81	if_no_pbl_nudging_ph
82	INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_ph
83	INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_ph
84	INTEGER, INTENT(IN) :: dk_zfac_uv, dk_zfac_t, dk_zfac_ph
85	INTEGER, INTENT(IN) :: if_ramping
86	INTEGER, INTENT(IN) :: xwavenum,ywavenum
87
88	INTEGER , INTENT(IN) :: id
89	REAL, INTENT(IN) :: DT, xtime, dtramp_min
90
91	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
92	ims,ime, jms,jme, kms,kme, &
93	kts,kte, num_tiles, &
94	ips,ipe,jps,jpe,kps,kpe, &
95	imsx,imex,jmsx,jmex,kmsx,kmex, &
96	ipsx,ipex,jpsx,jpex,kpsx,kpex, &
97	imsy,imey,jmsy,jmey,kmsy,kmey, &
98	ipsy,ipey,jpsy,jpey,kpsy,kpey
99
100	INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
101	& i_start,i_end,j_start,j_end
102
103	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
104	INTENT(IN) :: ph3d, &
105	th3d, &
106	z, &
107	z_at_w
108
109	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
110	INTENT(INOUT) :: rundgdten, &
111	rvndgdten, &
112	rthndgdten, &
113	rphndgdten
114
115
116	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
117	INTENT(INOUT) :: u_ndg_old, &
118	v_ndg_old, &
119	t_ndg_old, &
120	ph_ndg_old, &
121	u_ndg_new, &
122	v_ndg_new, &
123	t_ndg_new, &
124	ph_ndg_new
125
126	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
127	INTENT(IN) :: u3d, &
128	v3d
129
130	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
131	ht
132
133	REAL, INTENT(IN) :: guv, gt ,gph
134
135	INTEGER :: its,ite, jts,jte,ij
136	INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
137	REAL :: xtime_old, xtime_new, coef, val_analysis
138	INTEGER :: kpbl, dbg_level
139
140	REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
141	REAL, DIMENSION( ips:ipe, kps:kpe, jps:jpe, 4 ) :: wpbl ! 1: u, 2: v, 3: t, 4: ph
142	REAL, DIMENSION( kps:kpe, 4 ) :: wzfac ! 1: u, 2: v, 3: t, 4: ph
143
144	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
145
146	CHARACTER (LEN=256) :: message
147
148	#if ! ( NMM_CORE == 1 )
149
150	DO ij = 1 , num_tiles
151	its=i_start(ij)
152	ite=i_end(ij)
153	jts=j_start(ij)
154	jte=j_end(ij)
155	ENDDO
156	!GMM default values for vertical coefficients, set here
157	wpbl(:,:,:,:) = 1.0
158	wzfac(:,:) = 1.0
159
160	!$OMP PARALLEL DO &
161	!$OMP PRIVATE ( ij, i,j,k ) ! only for part of the calculation.
162	DO ij = 1 , num_tiles
163	DO j = jts, jte
164	DO k = kts, kte
165	DO i = its, ite
166	RUNDGDTEN(i,k,j) = 0.0
167	RVNDGDTEN(i,k,j) = 0.0
168	RTHNDGDTEN(i,k,j) = 0.0
169	RPHNDGDTEN(i,k,j) = 0.0
170	ENDDO
171	ENDDO
172	ENDDO
173	ENDDO
174	!$OMP END PARALLEL DO
175
176	actual_end_fdda_min = end_fdda_hour*60.0
177	IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
178	actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
179	IF( xtime > actual_end_fdda_min ) RETURN
180
181	!$OMP PARALLEL DO &
182	!$OMP PRIVATE ( ij, i,j,k ) ! only for part of the calculation.
183	!GMM Transpose and filtering needs to be done on patch
184
185	DO ij = 1 , num_tiles
186
187	! actual_end_fdda_min = end_fdda_hour*60.0
188	! IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
189	! actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
190
191	xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
192	xtime_new = xtime_old + analysis_interval
193	coef = (xtime-xtime_old)/(xtime_new-xtime_old)
194
195	IF ( wrf_dm_on_monitor()) THEN
196
197	CALL get_wrf_debug_level( dbg_level )
198
199	IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
200
201	IF( xtime < end_fdda_hour*60.0 ) THEN
202	WRITE(message,'(a,i1,a,f10.3,a)') &
203	'D0',id,' Spectral nudging read in new data at time = ', xtime, ' min.'
204	CALL wrf_message( TRIM(message) )
205	WRITE(message,'(a,i1,a,2f8.2,a)') &
206	'D0',id,' Spectral nudging bracketing times = ', xtime_old, xtime_new, ' min.'
207	CALL wrf_message( TRIM(message) )
208	ENDIF
209
210	actual_end_fdda_min = end_fdda_hour*60.0
211	IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
212	actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
213
214	IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
215	! Find the mid point of the tile and print out the sample values
216	i0 = (ite-its)/2+its
217	j0 = (jte-jts)/2+jts
218
219	IF( guv > 0.0 ) THEN
220	DO k = kts, kte
221	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
222	' D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
223	' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
224	CALL wrf_message( TRIM(message) )
225	ENDDO
226	DO k = kts, kte
227	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
228	' D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
229	' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
230	CALL wrf_message( TRIM(message) )
231	ENDDO
232	ENDIF
233
234	IF( gt > 0.0 ) THEN
235	DO k = kts, kte
236	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
237	' D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
238	' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
239	CALL wrf_message( TRIM(message) )
240	ENDDO
241	ENDIF
242
243	IF( gph > 0.0 ) THEN
244	DO k = kts, kte
245	WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
246	' D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
247	' ph_ndg_old=', ph_ndg_old(i0,k,j0), ' ph_ndg_new=', ph_ndg_new(i0,k,j0)
248	CALL wrf_message( TRIM(message) )
249	ENDDO
250	ENDIF
251
252	ENDIF
253	ENDIF
254	ENDIF
255
256	jtsv=MAX0(jts,jds+1)
257	itsu=MAX0(its,ids+1)
258
259	jtf=MIN0(jte,jde-1)
260	ktf=MIN0(kte,kde-1)
261	itf=MIN0(ite,ide-1)
262	!
263	! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t,
264	! if_no_pbl_nudging_ph swithes) are set to 1, compute the weighting function, wpbl(:,k,:,:),
265	! based on the PBL depth. wpbl = 1 above the PBL and wpbl = 0 in the PBL. If all
266	! the switche are set to zero, wpbl = 1 (default value).
267	!
268	!GMM If those are set to zero, then check if the user-defined namelist switches (if_zfac_uv, if_zfac_t,
269	! if_zfac_ph swithes) are set to 1, compute the weighting function, wzfac(k,:),
270	! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_ph)
271	! below which analysis nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x).
272	! The strength of nudging increases linearly from k_zfac to kzfac + dk_zfac
273	! (dk_zfac_uv, dk_zfac_t and kd_zfac_ph are also set in the namelist, default value is 1).
274	!If all switches are set to zero, wzfac = 1 (default value).
275	!
276
277	IF( if_no_pbl_nudging_uv == 1 ) THEN
278
279	DO j=jts,jtf
280	DO i=itsu,itf
281
282	kpbl = 1
283	zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
284
285	loop_ku: DO k=kts,ktf
286	zagl_bot = 0.5 * ( z_at_w(i-1,k, j)-ht(i-1,j) + z_at_w(i,k, j)-ht(i,j) )
287	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) )
288	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
289	kpbl = k
290	EXIT loop_ku
291	ENDIF
292	ENDDO loop_ku
293
294	DO k=kts,ktf
295	wpbl(i, k, j, 1) = max ( min ( float(k-kpbl) / float(dk_zfac_uv) , 1. ) , 0.)
296	ENDDO
297
298	ENDDO
299	ENDDO
300
301	DO i=its,itf
302	DO j=jtsv,jtf
303
304	kpbl = 1
305	zpbl = 0.5 * ( pblh(i,j-1) + pblh(i,j) )
306
307	loop_kv: DO k=kts,ktf
308	zagl_bot = 0.5 * ( z_at_w(i,k, j-1)-ht(i,j-1) + z_at_w(i,k, j)-ht(i,j) )
309	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) )
310	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
311	kpbl = k
312	EXIT loop_kv
313	ENDIF
314	ENDDO loop_kv
315
316	DO k=kts,ktf
317	wpbl(i, k, j, 2) = max ( min ( float(k-kpbl) / float(dk_zfac_uv) , 1. ) , 0.)
318	ENDDO
319
320	ENDDO
321	ENDDO
322
323	ELSEIF( if_zfac_uv == 1 ) THEN
324
325	DO j=jts,jte
326	DO k=kts,ktf
327	DO i=its,ite
328	wpbl(i, k, j, 1) = max ( min ( float(k-k_zfac_uv) / float(dk_zfac_uv) , 1. ) , 0.)
329	wpbl(i, k, j, 2) = wpbl(i, k, j, 1)
330	ENDDO
331	ENDDO
332	ENDDO
333
334	ENDIF
335
336
337	IF( if_no_pbl_nudging_t == 1 ) THEN
338
339	DO j=jts,jtf
340	DO i=its,itf
341
342	kpbl = 1
343	zpbl = pblh(i,j)
344
345	loop_kt: DO k=kts,ktf
346	zagl_bot = z_at_w(i,k, j)-ht(i,j)
347	zagl_top = z_at_w(i,k+1,j)-ht(i,j)
348	IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
349	kpbl = k
350	EXIT loop_kt
351	ENDIF
352	ENDDO loop_kt
353
354	DO k=kts,ktf
355	wpbl(i, k, j, 3) = max ( min ( float(k-kpbl) / float(dk_zfac_t) , 1. ) , 0.)
356	ENDDO
357
358	ENDDO
359	ENDDO
360
361	ELSEIF( if_zfac_t == 1 ) THEN
362
363	DO j=jts,jtf
364	DO k=kts,ktf
365	DO i=its,itf
366	wpbl(i, k, j, 3) = max ( min ( float(k-k_zfac_t) / float(dk_zfac_t) , 1. ) , 0.)
367	ENDDO
368	ENDDO
369	ENDDO
370
371	ENDIF
372
373
374	IF( if_no_pbl_nudging_ph == 1 ) THEN
375
376	DO j=jts,jtf
377	DO i=its,itf
378
379	kpbl = 1
380	zpbl = pblh(i,j)
381
382	loop_kph: DO k=kts,kte
383	zagl = z_at_w(i,k, j)-ht(i,j)
384	IF( zpbl >= zagl) THEN
385	kpbl = k
386	EXIT loop_kph
387	ENDIF
388	ENDDO loop_kph
389
390	DO k=kts,kte
391	wpbl(i, k, j, 4) = max ( min ( float(k-kpbl) / float(dk_zfac_ph) , 1. ) , 0.)
392	ENDDO
393
394	ENDDO
395	ENDDO
396
397	ELSEIF( if_zfac_ph == 1 ) THEN
398
399	DO j=jts,jtf
400	DO k=kts,kte
401	DO i=its,itf
402	wpbl(i, k, j, 4) = max ( min ( float(k-k_zfac_ph) / float(dk_zfac_ph) , 1. ) , 0.)
403	ENDDO
404	ENDDO
405	ENDDO
406
407	ENDIF
408
409
410	! If if_ramping and dtramp_min are defined by user, comput a time weighting function, tfac,
411	! for analysis nudging so that at the end of the nudging period (which has to be at a
412	! analysis time) we ramp down the nudging coefficient, based on the use-defined sign of dtramp_min.
413	!
414	! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at
415	! end_fdda_hour-ABS(dtramp_min).
416	!
417	! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at
418	! end_fdda_hour+ABS(dtramp_min). In this case, the obs values are extrapolated using
419	! the obs tendency saved from the previous FDDA wondow. More specifically for extrapolation,
420	! coef (see codes below) is recalculated to reflect extrapolation during the ramping period.
421	!
422	tfac = 1.0
423
424	IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
425
426	IF( dtramp_min <= 0.0 ) THEN
427	actual_end_fdda_min = end_fdda_hour*60.0
428	ELSE
429	actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
430	ENDIF
431
432	IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN
433	tfac = 1.0
434	ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
435	tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
436	IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval60.0)/(analysis_interval60.0)
437	ELSE
438	tfac = 0.0
439	ENDIF
440
441	ENDIF
442
443	ENDDO
444	!$OMP END PARALLEL DO
445
446	!
447	! Compute 3-D nudging tendencies for u, v
448	!
449	!
450	!GMM Fist calculate differences between model variables and analysis values,
451	!then filter in the x and y direction all wave numbers higher than
452	! xwavenum and ywavenum, as specified in the namelist.
453	!If either xwavenum or ywavenum are not specified,
454	! default values are zero, and spectral nudging is turned off
455	!Then use the filtered differences to calculate the spectral nudging tendencies
456
457	IF(guv > 0. ) then
458
459	!$OMP PARALLEL DO &
460	!$OMP PRIVATE ( ij, i,j,k )
461	DO ij = 1 , num_tiles
462
463	DO j=jts,jte
464	DO k=kts,ktf
465	DO i=its,ite
466	val_analysis = u_ndg_old(i,k,j) ( 1.0 - coef ) + u_ndg_new(i,k,j) coef
467
468	grid%dif_analysis(i,k,j) = val_analysis - u3d(i,k,j)
469
470	ENDDO
471	ENDDO
472	ENDDO
473
474	ENDDO
475	!$OMP END PARALLEL DO
476
477	!Filter
478
479	#ifdef DM_PARALLEL
480	# include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
481
482	CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum, &
483	ids, ide, jds, jde, kds, kde, &
484	imsx, imex, jmsx, jmex, kmsx, kmex, &
485	ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
486
487	# include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
488
489	CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum, &
490	ids, ide, jds, jde, kds, kde, &
491	imsy, imey, jmsy, jmey, kmsy, kmey, &
492	ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
493
494	# include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
495
496
497	#else
498	CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum, &
499	ids, ide, jds, jde, kds, kde, &
500	ims, ime, jms, jme, kms, kme, &
501	ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
502
503	CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum, &
504	ids, ide, jds, jde, kds, kde, &
505	ims, ime, jms, jme, kms, kme, &
506	ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
507	#endif
508
509	! Calculate tendency
510
511	!$OMP PARALLEL DO &
512	!$OMP PRIVATE ( ij, i,j,k )
513	DO ij = 1 , num_tiles
514
515	DO j=jts,jte
516	DO k=kts,ktf
517	DO i=its,ite
518	RUNDGDTEN(i,k,j) = guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
519	grid%dif_analysis(i,k,j)
520	ENDDO
521	ENDDO
522	ENDDO
523
524
525	!
526	! Now V component
527	!
528
529	DO j=jts,jte
530	DO k=kts,ktf
531	DO i=its,ite
532	val_analysis = v_ndg_old(i,k,j) ( 1.0 - coef ) + v_ndg_new(i,k,j) coef
533
534	grid%dif_analysis(i,k,j) = val_analysis - v3d(i,k,j)
535
536	ENDDO
537	ENDDO
538	ENDDO
539
540	ENDDO
541	!$OMP END PARALLEL DO
542
543	!
544	! Filter
545	!
546	#ifdef DM_PARALLEL
547	# include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
548	CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum, &
549	ids, ide, jds, jde, kds, kde, &
550	imsx, imex, jmsx, jmex, kmsx, kmex, &
551	ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
552
553	# include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
554	CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum, &
555	ids, ide, jds, jde, kds, kde-1, &
556	imsy, imey, jmsy, jmey, kmsy, kmey, &
557	ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
558
559	# include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
560
561
562	#else
563	CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum, &
564	ids, ide, jds, jde, kds, kde, &
565	ims, ime, jms, jme, kms, kme, &
566	ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
567
568	CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum, &
569	ids, ide, jds, jde, kds, kde, &
570	ims, ime, jms, jme, kms, kme, &
571	ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
572	#endif
573
574	! Calculate tendency
575
576	!$OMP PARALLEL DO &
577	!$OMP PRIVATE ( ij, i,j,k )
578	DO ij = 1 , num_tiles
579
580	DO j=jts,jte
581	DO k=kts,ktf
582	DO i=its,ite
583	RVNDGDTEN(i,k,j) = guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
584	grid%dif_analysis(i,k,j)
585	ENDDO
586	ENDDO
587	ENDDO
588
589	ENDDO
590	!$OMP END PARALLEL DO
591
592	ENDIF
593
594	IF(gt > 0. ) then
595
596	!$OMP PARALLEL DO &
597	!$OMP PRIVATE ( ij, i,j,k )
598	DO ij = 1 , num_tiles
599
600	DO j=jts,jte
601	DO k=kts,ktf
602	DO i=its,ite
603	val_analysis = t_ndg_old(i,k,j) ( 1.0 - coef ) + t_ndg_new(i,k,j) coef
604
605	grid%dif_analysis(i,k,j) = val_analysis - th3d(i,k,j) + 300.
606
607	ENDDO
608	ENDDO
609	ENDDO
610
611	ENDDO
612	!$OMP END PARALLEL DO
613
614	!Filter
615
616	#ifdef DM_PARALLEL
617	# include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
618
619	CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum, &
620	ids, ide, jds, jde, kds, kde, &
621	imsx, imex, jmsx, jmex, kmsx, kmex, &
622	ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
623
624	# include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
625
626	CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum, &
627	ids, ide, jds, jde, kds, kde, &
628	imsy, imey, jmsy, jmey, kmsy, kmey, &
629	ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
630
631	# include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
632
633
634	#else
635	CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum, &
636	ids, ide, jds, jde, kds, kde, &
637	ims, ime, jms, jme, kms, kme, &
638	ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
639
640	CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum, &
641	ids, ide, jds, jde, kds, kde, &
642	ims, ime, jms, jme, kms, kme, &
643	ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
644	#endif
645
646	! Calculate tendency
647
648	!$OMP PARALLEL DO &
649	!$OMP PRIVATE ( ij, i,j,k )
650	DO ij = 1 , num_tiles
651
652	DO j=jts,jte
653	DO k=kts,ktf
654	DO i=its,ite
655	RTHNDGDTEN(i,k,j) = gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
656	grid%dif_analysis(i,k,j)
657	ENDDO
658	ENDDO
659	ENDDO
660
661	ENDDO
662	!$OMP END PARALLEL DO
663
664	ENDIF
665
666	IF(gph > 0. ) then
667
668	!$OMP PARALLEL DO &
669	!$OMP PRIVATE ( ij, i,j,k )
670	DO ij = 1 , num_tiles
671
672	DO j=jts,jte
673	DO k=kts,kte
674	DO i=its,ite
675	val_analysis = ph_ndg_old(i,k,j) ( 1.0 - coef ) + ph_ndg_new(i,k,j) coef
676
677	grid%dif_analysis(i,k,j) = val_analysis - ph3d(i,k,j)
678
679	ENDDO
680	ENDDO
681	ENDDO
682
683	ENDDO
684	!$OMP END PARALLEL DO
685
686	!Filter
687
688	#ifdef DM_PARALLEL
689	# include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
690
691	CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum, &
692	ids, ide, jds, jde, kds, kde, &
693	imsx, imex, jmsx, jmex, kmsx, kmex, &
694	ipsx, ipex, jpsx, jpex, kpsx, kpex )
695
696	# include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
697
698	CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum, &
699	ids, ide, jds, jde, kds, kde, &
700	imsy, imey, jmsy, jmey, kmsy, kmey, &
701	ipsy, ipey, jpsy, jpey, kpsy, kpey )
702
703	# include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
704
705
706	#else
707	CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum, &
708	ids, ide, jds, jde, kds, kde, &
709	ims, ime, jms, jme, kms, kme, &
710	ips, ipe, jps, jpe, kps, kpe )
711
712	CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum, &
713	ids, ide, jds, jde, kds, kde, &
714	ims, ime, jms, jme, kms, kme, &
715	ips, ipe, jps, jpe, kps, kpe )
716	#endif
717
718	! Calculate tendency
719
720	!$OMP PARALLEL DO &
721	!$OMP PRIVATE ( ij, i,j,k )
722	DO ij = 1 , num_tiles
723
724	DO j=jts,jte
725	DO k=kts,kte
726	DO i=its,ite
727	RPHNDGDTEN(i,k,j) = gph * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
728	grid%dif_analysis(i,k,j)
729	ENDDO
730	ENDDO
731	ENDDO
732
733	ENDDO
734	!$OMP END PARALLEL DO
735
736	ENDIF
737
738	#endif
739
740	END SUBROUTINE spectral_nudging
741
742	!------------------------------------------------------------------------------
743
744	SUBROUTINE spectral_nudging_filter_3dx( f, nwave, &
745	ids, ide, jds, jde, kds, kde, &
746	ims, ime, jms, jme, kms, kme, &
747	its, ite, jts, jte, kts, kte )
748
749	IMPLICIT NONE
750
751	INTEGER , INTENT(IN ) :: nwave
752	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
753	ims, ime, jms, jme, kms, kme, &
754	its, ite, jts, jte, kts, kte
755
756	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: f
757
758	REAL , DIMENSION(1:ide-ids+1,1:kte-kts+1) :: sheet
759	INTEGER :: i, j, j_end, k, nx, ny
760
761	! Variables will stay in domain form since this routine is meaningless
762	! unless tile extent is the same as domain extent in E/W direction, i.e.,
763	! the processor has access to all grid points in E/W direction.
764	! There may be other ways of doing FFTs, but we haven't learned them yet...
765
766	! Check to make sure we have full access to all E/W points
767	IF ((its /= ids) .OR. (ite /= ide)) THEN
768	WRITE ( wrf_err_message , * ) 'module_spectral_nudging: 3d: (its /= ids) or (ite /= ide)',its,ids,ite,ide
769	CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
770	END IF
771
772
773	nx = ide-ids+1
774	ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
775	j_end = MIN(jte, jde-1)
776	IF (j_end == jde-1) j_end = jde
777	DO j = jts, j_end
778
779	DO k=kts,kte
780	DO i=ids,ide-1
781	sheet(i-ids+1,k-kts+1) = f(i,k,j)
782	END DO
783	sheet(ide,k-kts+1) = 0.
784	END DO
785
786	CALL spectral_nudging_filter_fft_2d_ncar(nx,ny,nwave,sheet)
787
788	DO k=kts,kte
789	DO i=ids,ide
790	f(i,k,j) = sheet(i-ids+1,k-kts+1)
791	END DO
792	END DO
793	END DO ! outer j (latitude) loop
794
795	END SUBROUTINE spectral_nudging_filter_3dx
796	!------------------------------------------------------------------------------
797
798	SUBROUTINE spectral_nudging_filter_3dy( f, nwave, &
799	ids, ide, jds, jde, kds, kde, &
800	ims, ime, jms, jme, kms, kme, &
801	its, ite, jts, jte, kts, kte )
802
803	IMPLICIT NONE
804
805	INTEGER , INTENT(IN ) :: nwave
806	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
807	ims, ime, jms, jme, kms, kme, &
808	its, ite, jts, jte, kts, kte
809
810	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: f
811
812	REAL , DIMENSION(1:jde-jds+1,1:kte-kts+1) :: sheet
813	INTEGER :: i, j, i_end, k, nx, ny
814
815	! Variables will stay in domain form since this routine is meaningless
816	! unless tile extent is the same as domain extent in S/N direction, i.e.,
817	! the processor has access to all grid points in S/N direction.
818	! There may be other ways of doing FFTs, but we haven't learned them yet...
819
820	! Check to make sure we have full access to all S/N points
821	IF ((jts /= jds) .OR. (jte /= jde)) THEN
822	WRITE ( wrf_err_message , * ) 'module_spectral_nudging: 3d: (jts /= jds) or (jte /= jde)',jts,ids,ite,ide
823	CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
824	END IF
825
826
827	nx = jde-jds+1
828	ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
829	i_end = MIN(ite, ide-1)
830	IF (i_end == ide-1) i_end = ide
831	DO i = its, i_end
832
833	DO k=kts,kte
834	DO j=jds,jde
835	sheet(j-jds+1,k-kts+1) = f(i,k,j)
836	END DO
837	sheet(jde,k-kts+1) = 0.
838	END DO
839
840	CALL spectral_nudging_filter_fft_2d_ncar(nx,ny,nwave,sheet)
841
842	DO k=kts,kte
843	DO j=jds,jde
844	f(i,k,j) = sheet(j-jds+1,k-kts+1)
845	END DO
846	END DO
847	END DO ! outer i (longitude) loop
848
849	END SUBROUTINE spectral_nudging_filter_3dy
850
851	!------------------------------------------------------------------------------
852
853	SUBROUTINE spectral_nudging_filter_fft_2d_ncar(nx,ny,nwave,fin)
854	IMPLICIT NONE
855	INTEGER , INTENT(IN) :: nx, ny, nwave
856	REAL , DIMENSION(nx,ny), INTENT(INOUT) :: fin
857
858	INTEGER :: i, j
859	REAL, dimension(nx,ny) :: fp
860
861	INTEGER :: lensave, ier, nh, n1
862	INTEGER :: lot, jump, n, inc, lenr, lensav, lenwrk
863	REAL, DIMENSION(nx+15) :: wsave
864	REAL, DIMENSION(nx,ny) :: work
865
866
867
868	! we are following the naming convention of the fftpack5 routines
869
870	n = nx
871	lot = ny
872	lensav = n+15
873	inc = 1
874	lenr = nx*ny
875	jump = nx
876	lenwrk = lenr
877
878	! forward transform
879	! initialize coefficients, place in wsave
880	! (should place this in init and save wsave at program start)
881
882	call rfftmi(n,wsave,lensav,ier)
883	IF(ier /= 0) THEN
884	write(0,*) ' error in rfftmi ',ier
885	END IF
886
887	! do the forward transform
888
889	call rfftmf( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
890	IF(ier /= 0) THEN
891	write(0,*) ' error in rfftmf ',ier
892	END IF
893
894	if(MOD(n,2) == 0) then
895	nh = n/2 - 1
896	else
897	nh = (n-1)/2
898	end if
899
900
901	! filter all waves with wavenumber larger than nwave
902
903	fp = 1.
904
905	DO j=1,ny
906	DO i=nwave+1,nh
907	fp(2*i-1,j) = 0.
908	fp(2*i,j) = 0.
909	ENDDO
910	ENDDO
911
912	DO j=1,ny
913	DO i=1,nx
914	fin(i,j) = fp(i,j)*fin(i,j)
915	ENDDO
916	ENDDO
917
918	! do the backward transform
919
920	call rfftmb( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
921	IF(ier /= 0) THEN
922	write(0,*) ' error in rfftmb ',ier
923	END IF
924
925	END SUBROUTINE spectral_nudging_filter_fft_2d_ncar
926
927	!-------------------------------------------------------------------
928
929	SUBROUTINE fddaspnudginginit(id,rundgdten,rvndgdten,rthndgdten,rphndgdten, &
930	run_hours, &
931	if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph, &
932	if_zfac_uv, k_zfac_uv, dk_zfac_uv, &
933	if_zfac_t, k_zfac_t, dk_zfac_t, &
934	if_zfac_ph, k_zfac_ph, dk_zfac_ph, &
935	guv, gt, gph, if_ramping, dtramp_min, end_fdda_hour, &
936	xwavenum,ywavenum, &
937	restart, allowed_to_read, &
938	ids, ide, jds, jde, kds, kde, &
939	ims, ime, jms, jme, kms, kme, &
940	its, ite, jts, jte, kts, kte )
941	!-------------------------------------------------------------------
942	IMPLICIT NONE
943	!-------------------------------------------------------------------
944	!
945	INTEGER , INTENT(IN) :: id
946	LOGICAL, INTENT(IN) :: restart, allowed_to_read
947	INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
948	ims, ime, jms, jme, kms, kme, &
949	its, ite, jts, jte, kts, kte
950	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
951	rundgdten, &
952	rvndgdten, &
953	rthndgdten, &
954	rphndgdten
955	INTEGER, INTENT(IN) :: run_hours
956	INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
957	if_no_pbl_nudging_ph, end_fdda_hour
958	INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_ph
959	INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_ph
960	INTEGER, INTENT(IN) :: dk_zfac_uv, dk_zfac_t, dk_zfac_ph
961	INTEGER, INTENT(IN) :: if_ramping
962	INTEGER, INTENT(IN) :: xwavenum,ywavenum
963	REAL, INTENT(IN) :: dtramp_min
964	REAL, INTENT(IN) :: guv, gt, gph
965	REAL :: actual_end_fdda_min
966
967	INTEGER :: i, j, k
968
969	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
970
971	CHARACTER (LEN=256) :: message
972
973	IF ( wrf_dm_on_monitor() ) THEN
974
975	IF( guv > 0.0) THEN
976	WRITE(message,'(a,i1,a,e12.4,a,i4,a,i4)') &
977	'D0',id,' Spectral nudging for wind is turned on and Guv= ', guv,' xwave= ',xwavenum,' ywavenum= ',ywavenum
978	CALL wrf_message(TRIM(message))
979	ELSE IF( guv < 0.0 ) THEN
980	CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
981	ELSE
982	WRITE(message,'(a,i1,a,e12.4)') &
983	'D0',id,' Spectral nudging for wind is turned off and Guv= ', guv
984	CALL wrf_message(TRIM(message))
985	ENDIF
986
987	IF( gt > 0.0) THEN
988	WRITE(message,'(a,i1,a,e12.4,a,i4,a,i4)') &
989	'D0',id,' Spectral nudging for temperature is turned on and Gt= ', gt,' xwave= ',xwavenum,' ywavenum= ',ywavenum
990	CALL wrf_message(TRIM(message))
991	ELSE IF( gt < 0.0 ) THEN
992	CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
993	ELSE
994	WRITE(message,'(a,i1,a,e12.4)') &
995	'D0',id,' Spectral nudging for temperature is turned off and Gt= ', gt
996	CALL wrf_message(TRIM(message))
997	ENDIF
998
999	IF( gph > 0.0) THEN
1000	WRITE(message,'(a,i1,a,e12.4,a,i4,a,i4)') &
1001	'D0',id,' Spectral nudging for geopotential is turned on and Gph= ', gph,' xwave= ',xwavenum,' ywavenum= ',ywavenum
1002	CALL wrf_message(TRIM(message))
1003	ELSE IF( gph < 0.0 ) THEN
1004	CALL wrf_error_fatal('In grid FDDA, Gph must be positive.')
1005	ELSE
1006	WRITE(message,'(a,i1,a,e12.4)') &
1007	'D0',id,' Spectral nudging for geopotential is turned off and Gph= ', gph
1008	CALL wrf_message(TRIM(message))
1009	ENDIF
1010
1011	IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
1012	WRITE(message,'(a,i1,a)') &
1013	'D0',id,' Spectral nudging for wind is turned off within the PBL.'
1014	CALL wrf_message(TRIM(message))
1015	IF( dk_zfac_uv < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_uv must be greater or equal than 1.')
1016	ELSEIF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
1017	WRITE(message,'(a,i1,a,i3)') &
1018	'D0',id,' Spectral nudging for wind is turned off below layer', k_zfac_uv
1019	CALL wrf_message(TRIM(message))
1020	IF( dk_zfac_uv < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_uv must be greater or equal than 1.')
1021	ENDIF
1022
1023
1024	IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
1025	WRITE(message,'(a,i1,a)') &
1026	'D0',id,' Spectral nudging for temperature is turned off within the PBL.'
1027	CALL wrf_message(TRIM(message))
1028	IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
1029	ELSEIF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
1030	WRITE(message,'(a,i1,a,i3)') &
1031	'D0',id,' Spectral nudging for temperature is turned off below layer', k_zfac_t
1032	CALL wrf_message(TRIM(message))
1033	IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
1034	ENDIF
1035
1036
1037	IF( gph > 0.0 .AND. if_no_pbl_nudging_ph == 1 ) THEN
1038	WRITE(message,'(a,i1,a)') &
1039	'D0',id,' Spectral nudging for geopotential is turned off within the PBL.'
1040	CALL wrf_message(TRIM(message))
1041	IF( dk_zfac_ph < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_ph must be greater or equal than 1.')
1042	ELSEIF( gph > 0.0 .AND. if_zfac_ph == 1 ) THEN
1043	WRITE(message,'(a,i1,a,i3)') &
1044	'D0',id,' Spectral nudging for geopotential is turned off below layer', &
1045	k_zfac_ph
1046	CALL wrf_message(TRIM(message))
1047	IF( dk_zfac_ph < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_ph must be greater or equal than 1.')
1048	ENDIF
1049
1050	IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
1051	IF( dtramp_min <= 0.0 ) THEN
1052	actual_end_fdda_min = end_fdda_hour*60.0
1053	ELSE
1054	actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
1055	ENDIF
1056
1057	IF( actual_end_fdda_min <= run_hours*60. ) THEN
1058	WRITE(message,'(a,i1,a)') &
1059	'D0',id,' Spectral nudging is ramped down near the end of the nudging period,'
1060	CALL wrf_message(TRIM(message))
1061
1062	WRITE(message,'(a,f6.2,a,f6.2,a)') &
1063	' starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0,&
1064	'h, ending at ', actual_end_fdda_min/60.0,'h.'
1065	CALL wrf_message(TRIM(message))
1066	ENDIF
1067	ENDIF
1068
1069	ENDIF
1070
1071	IF(.not.restart) THEN
1072	DO j = jts,jte
1073	DO k = kts,kte
1074	DO i = its,ite
1075	rundgdten(i,k,j) = 0.
1076	rvndgdten(i,k,j) = 0.
1077	rthndgdten(i,k,j) = 0.
1078	rphndgdten(i,k,j) = 0.
1079	ENDDO
1080	ENDDO
1081	ENDDO
1082	ENDIF
1083
1084	END SUBROUTINE fddaspnudginginit
1085	!-------------------------------------------------------------------
1086
1087	END MODULE module_fdda_spnudging

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

Download in other formats:

Original Format