Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

wrf_fddaobs_in.F @ 623

Last change on this file since 623 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: 76.8 KB

Line
1	!WRF:MEDIATION_LAYER:IO
2	! ---
3
4	! This obs-nudging FDDA module (RTFDDA) is developed by the
5	! NCAR/RAL/NSAP (National Security Application Programs), under the
6	! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is
7	! acknowledged for releasing this capability for WRF community
8	! research applications.
9	!
10	! The NCAR/RAL RTFDDA module was adapted, and significantly modified
11	! from the obs-nudging module in the standard MM5V3.1 which was originally
12	! developed by PSU (Stauffer and Seaman, 1994).
13	!
14	! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module
15	! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
16	! Nov. 2006
17	!
18	! References:
19	!
20	! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
21	! implementation of obs-nudging-based FDDA into WRF for supporting
22	! ATEC test operations. 2005 WRF user workshop. Paper 10.7.
23	!
24	! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update
25	! on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE
26	! and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7.
27
28	!
29	! Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data
30	! assimilation. J. Appl. Meteor., 33, 416-434.
31	!
32	! http://www.rap.ucar.edu/projects/armyrange/references.html
33	!
34
35	SUBROUTINE wrf_fddaobs_in (grid ,config_flags)
36
37	USE module_domain
38	USE module_configure
39	USE module_model_constants !rovg
40
41	IMPLICIT NONE
42	TYPE(domain) :: grid
43	TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
44	#if ( EM_CORE == 1 )
45
46	! Local variables
47	integer :: ktau ! timestep index corresponding to xtime
48	integer :: krest ! restart timestep
49	integer :: inest ! nest level
50	integer :: infreq ! input frequency
51	integer :: nstlev ! nest level
52	real :: dtmin ! dt in minutes
53	real :: xtime ! forecast time in minutes
54	logical :: iprt_in4dob ! print flag
55
56	INTEGER ids , ide , jds , jde , kds , kde , &
57	ims , ime , jms , jme , kms , kme , &
58	ips , ipe , jps , jpe , kps , kpe
59	INTEGER ij, its, ite, jts, jte
60
61	! Modified to also call in4dob intially, since subr in4dob is no
62	! longer called from subr fddaobs_init. Note that itimestep is now
63	! the model step BEFORE the model integration step, because this
64	! routine is now called by med_before_solve_io.
65	ktau = grid%itimestep ! ktau corresponds to xtime
66	krest = grid%fdob%ktaur
67	inest = grid%grid_id
68	nstlev = grid%fdob%levidn(inest)
69	infreq = grid%obs_ionf(grid%parent_grid_ratio*nstlev)
70	iprt_in4dob = grid%obs_ipf_in4dob
71
72	IF( (ktau.GT.krest.AND.MOD(ktau,infreq).EQ.0) &
73	.OR.(ktau.EQ.krest) ) then
74	! Calculate forecast time.
75	dtmin = grid%dt/60.
76	xtime = grid%xtime
77
78	CALL get_ijk_from_grid ( grid , &
79	ids, ide, jds, jde, kds, kde, &
80	ims, ime, jms, jme, kms, kme, &
81	ips, ipe, jps, jpe, kps, kpe )
82
83	!$OMP PARALLEL DO &
84	!$OMP PRIVATE ( ij )
85
86	DO ij = 1 , grid%num_tiles
87	its = grid%i_start(ij)
88	ite = min(grid%i_end(ij),ide-1)
89	jts = grid%j_start(ij)
90	jte = min(grid%j_end(ij),jde-1)
91
92	CALL in4dob(inest, xtime, ktau, krest, dtmin, &
93	grid%julyr, grid%julday, grid%gmt, & !obsnypatch
94	grid%obs_nudge_opt, grid%obs_nudge_wind, grid%obs_nudge_temp, &
95	grid%obs_nudge_mois, grid%obs_nudge_pstr, grid%obs_coef_wind, &
96	grid%obs_coef_temp, grid%obs_coef_mois, grid%obs_coef_pstr, &
97	grid%obs_rinxy, grid%obs_rinsig, grid%fdob%window, &
98	grid%obs_npfi, grid%obs_ionf, &
99	grid%obs_prt_max, grid%obs_prt_freq, &
100	grid%obs_idynin, &
101	grid%obs_dtramp, grid%fdob, grid%fdob%varobs, &
102	grid%fdob%timeob, grid%fdob%nlevs_ob, grid%fdob%lev_in_ob, &
103	grid%fdob%plfo, grid%fdob%elevob, grid%fdob%rio, &
104	grid%fdob%rjo, grid%fdob%rko, &
105	grid%xlat, grid%xlong, &
106	config_flags%cen_lat, &
107	config_flags%cen_lon, &
108	config_flags%stand_lon, &
109	config_flags%truelat1, config_flags%truelat2, &
110	grid%fdob%known_lat, grid%fdob%known_lon, &
111	config_flags%dx, config_flags%dy, rovg, t0, &
112	grid%fdob%obsprt, &
113	grid%fdob%latprt, grid%fdob%lonprt, &
114	grid%fdob%mlatprt, grid%fdob%mlonprt, &
115	grid%fdob%stnidprt, &
116	ide, jde, &
117	ims, ime, jms, jme, &
118	its, ite, jts, jte, &
119	config_flags%map_proj, &
120	model_config_rec%parent_grid_ratio, &
121	model_config_rec%i_parent_start(inest), &
122	model_config_rec%j_parent_start(inest), &
123	model_config_rec%max_dom, &
124	model_config_rec%nobs_ndg_vars, grid%max_obs, iprt_in4dob)
125	ENDDO
126	!$OMP END PARALLEL DO
127
128	ENDIF
129
130	RETURN
131	#endif
132	END SUBROUTINE wrf_fddaobs_in
133	#if ( EM_CORE == 1 )
134	!------------------------------------------------------------------------------
135	! Begin subroutine in4dob and its subroutines
136	!------------------------------------------------------------------------------
137	SUBROUTINE in4dob(inest, xtime, ktau, ktaur, dtmin, &
138	myear, julday, gmt, & !obsnypatch
139	nudge_opt, iswind, istemp, &
140	ismois, ispstr, giv, &
141	git, giq, gip, &
142	rinxy, rinsig, twindo, &
143	npfi, ionf, prt_max, prt_freq, idynin, &
144	dtramp, fdob, varobs, &
145	timeob, nlevs_ob, lev_in_ob, &
146	plfo, elevob, rio, &
147	rjo, rko, &
148	xlat, xlong, &
149	cen_lat, &
150	cen_lon, &
151	stand_lon, &
152	true_lat1, true_lat2, &
153	known_lat, known_lon, &
154	dxm, dym, rovg, t0, &
155	obs_prt, &
156	lat_prt, lon_prt, &
157	mlat_prt, mlon_prt, &
158	stnid_prt, &
159	e_we, e_sn, &
160	ims, ime, jms, jme, &
161	its, ite, jts, jte, &
162	map_proj, &
163	parent_grid_ratio, &
164	i_parent_start, &
165	j_parent_start, &
166	maxdom, &
167	nndgv, niobf, iprt)
168
169	USE module_domain
170	USE module_model_constants, ONLY : rcp
171	USE module_date_time , ONLY : geth_idts
172	USE module_llxy
173
174	IMPLICIT NONE
175
176	! THIS IS SUBROUTINE READS AN OBSERVATION DATA FILE AND
177	! SELECTS ONLY THOSE VALUES OBSERVED AT TIMES THAT FALL
178	! WITHIN A TIME WINDOW (TWINDO) CENTERED ABOUT THE CURRENT
179	! FORECAST TIME (XTIME). THE INCOMING OBS FILES MUST BE
180	! IN CHRONOLOGICAL ORDER.
181	!
182	! NOTE: This routine was originally designed for MM5, which uses
183	! a nonstandard (I,J) coordinate system. For WRF, I is the
184	! east-west running coordinate, and J is the south-north
185	! running coordinate. So "J-slab" here is west-east in
186	! extent, not south-north as for MM5. RIO and RJO have
187	! the opposite orientation here as for MM5. -ajb 06/10/2004
188
189	! NOTE - IN4DOB IS CALLED ONLY FOR THE COARSE MESH TIMES
190	! - VAROBS(IVAR,N) HOLDS THE OBSERVATIONS.
191	! IVAR=1 OBS U
192	! IVAR=2 OBS V
193	! IVAR=3 OBS T
194	! IVAR=4 OBS Q
195	! IVAR=5 OBS Pressure
196	! IVAR=6 OBS Height
197
198	INTEGER, intent(in) :: niobf ! maximum number of observations
199	INTEGER, intent(in) :: nndgv ! number of nudge variables
200	INTEGER, intent(in) :: INEST ! nest level
201	REAL, intent(in) :: xtime ! model time in minutes
202	INTEGER, intent(in) :: KTAU ! current timestep
203	INTEGER, intent(in) :: KTAUR ! restart timestep
204	REAL, intent(in) :: dtmin ! dt in minutes
205	INTEGER, intent(in) :: myear ! model year !obsnypatch
206	INTEGER, intent(in) :: julday ! Julian day
207	REAL, intent(in) :: gmt ! Model GMT at start of run
208	INTEGER, intent(in) :: nudge_opt ! obs-nudge flag for this nest
209	INTEGER, intent(in) :: iswind ! nudge flag for wind
210	INTEGER, intent(in) :: istemp ! nudge flag for temperature
211	INTEGER, intent(in) :: ismois ! nudge flag for moisture
212	INTEGER, intent(in) :: ispstr ! nudge flag for pressure (obsolete)
213	REAL, intent(in) :: giv ! coefficient for wind
214	REAL, intent(in) :: git ! coefficient for temperature
215	REAL, intent(in) :: giq ! coefficient for moisture
216	REAL, intent(in) :: gip ! coefficient for pressure
217	REAL, intent(in) :: rinxy ! horizontal radius of influence (km)
218	REAL, intent(in) :: rinsig ! vertical radius of influence (on sigma)
219	REAL, intent(inout) :: twindo ! (time window)/2 (min) for nudging
220	INTEGER, intent(in) :: npfi ! coarse-grid time-step frequency for diagnostics
221	INTEGER, intent(in) :: ionf ! coarse-grid time-step frequency for obs-nudging calcs
222	INTEGER, intent(in) :: prt_max ! max number of entries of obs for diagnostic printout
223	INTEGER, intent(in) :: prt_freq ! frequency (in obs index) for diagnostic printout.
224	INTEGER, intent(in) :: idynin ! for dynamic initialization using a ramp-down function
225	REAL, intent(in) :: dtramp ! time period in minutes for ramping
226	TYPE(fdob_type), intent(inout) :: fdob ! derived data type for obs data
227	REAL, intent(inout) :: varobs(nndgv,niobf) ! observational values in each variable
228	REAL, intent(inout) :: timeob(niobf) ! model times for each observation (hours)
229	REAL, intent(inout) :: nlevs_ob(niobf) ! numbers of levels in sounding obs
230	REAL, intent(inout) :: lev_in_ob(niobf) ! level in sounding-type obs
231	REAL, intent(inout) :: plfo(niobf) ! index for type of obs-platform
232	REAL, intent(inout) :: elevob(niobf) ! elevations of observations (meters)
233	REAL, intent(inout) :: rio(niobf) ! west-east grid coordinate (non-staggered grid)
234	REAL, intent(inout) :: rjo(niobf) ! south-north grid coordinate (non-staggered grid)
235	REAL, intent(inout) :: rko(niobf) ! vertical grid coordinate
236	REAL, DIMENSION( ims:ime, jms:jme ), &
237	INTENT(IN ) :: xlat, xlong ! lat/lon on mass-pt grid (for diagnostics only)
238	REAL, intent(in) :: cen_lat ! center latitude for map projection
239	REAL, intent(in) :: cen_lon ! center longitude for map projection
240	REAL, intent(in) :: stand_lon ! standard longitude for map projection
241	REAL, intent(in) :: true_lat1 ! truelat1 for map projection
242	REAL, intent(in) :: true_lat2 ! truelat2 for map projection
243	REAL, intent(in) :: known_lat ! latitude of domain origin point (i,j)=(1,1)
244	REAL, intent(in) :: known_lon ! longigude of domain origin point (i,j)=(1,1)
245	REAL, intent(in) :: dxm ! grid size in x (meters)
246	REAL, intent(in) :: dym ! grid size in y (meters)
247	REAL, intent(in) :: rovg ! constant rho over g
248	REAL, intent(in) :: t0 ! background temperature
249	INTEGER, intent(inout) :: obs_prt(prt_max) ! For printout of obs index
250	REAL, intent(inout) :: lat_prt(prt_max) ! For printout of obs latitude
251	REAL, intent(inout) :: lon_prt(prt_max) ! For printout of obs longitude
252	REAL, intent(inout) :: mlat_prt(prt_max) ! For printout of model lat at obs (ri,rj)
253	REAL, intent(inout) :: mlon_prt(prt_max) ! For printout of model lon at obs (ri,rj)
254	INTEGER, intent(inout) :: stnid_prt(40,prt_max) ! For printout of model lon at obs (ri,rj)
255	INTEGER, intent(in) :: e_we ! max grid index in south-north coordinate
256	INTEGER, intent(in) :: e_sn ! max grid index in west-east coordinate
257	INTEGER, intent(in) :: ims ! grid memory start index (west-east dim)
258	INTEGER, intent(in) :: ime ! grid memory end index (west-east dim)
259	INTEGER, intent(in) :: jms ! grid memory start index (south-north dim)
260	INTEGER, intent(in) :: jme ! grid memory end index (south-north dim)
261	INTEGER, intent(in) :: its ! grid tile start index (west-east dim)
262	INTEGER, intent(in) :: ite ! grid tile end index (west-east dim)
263	INTEGER, intent(in) :: jts ! grid tile start index (south-north dim)
264	INTEGER, intent(in) :: jte ! grid tile end index (south-north dim)
265	INTEGER, intent(in) :: map_proj ! map projection index
266	INTEGER, intent(in) :: parent_grid_ratio ! parent to nest grid ration
267	INTEGER, intent(in) :: i_parent_start ! starting i coordinate in parent domain
268	INTEGER, intent(in) :: j_parent_start ! starting j coordinate in parent domain
269	INTEGER, intent(in) :: maxdom ! maximum number of domains
270	LOGICAL, intent(in) :: iprt ! print flag
271
272	!*** DECLARATIONS FOR IMPLICIT NONE
273	integer :: n, ndum, nopen, nvol, idate, imm, iss
274	integer :: nlast ! last obs in list of valid obs from prev window
275	integer :: nsta ! number of stations held in timeobs array
276	integer :: nstaw ! # stations within the actual time window
277	integer :: nprev ! previous n in obs read loop (for printout only)
278	integer :: meas_count, imc, njend, njc, njcc, julob, kn
279	real :: hourob, rjulob
280	real :: xhour, tback, tforwd, rjdate1, timanl1, rtimob
281	real :: rj, ri, elevation, pressure_data
282	real :: pressure_qc, height_data, height_qc, temperature_data
283	real :: temperature_qc, u_met_data, u_met_qc, v_met_data
284	real :: v_met_qc, rh_data, rh_qc, r_data, slp_data, slp_qc
285	real :: ref_pres_data, ref_pres_qc, psfc_data, psfc_qc
286	real :: precip_data, precip_qc, tbar, twdop
287	real*8 :: tempob
288	INTEGER, EXTERNAL :: nvals_le_limit ! for finding #obs with timeobs <= tforwd
289
290	! Local variables
291	TYPE (PROJ_INFO) :: obs_proj ! Structure for obs projection information.
292	character*14 date_char
293	character*19 obs_date !obsnypatch
294	integer idts !obsnypatch
295	character*40 platform,source,id,namef
296	character*2 fonc
297	character(len=200) :: msg ! Argument to wrf_message
298	real latitude,longitude
299	logical :: newpass ! New pass flag (used for printout)
300	logical is_sound,bogus
301	LOGICAL OPENED,exist
302	integer :: ieof(5),ifon(5)
303	data ieof/0,0,0,0,0/
304	data ifon/0,0,0,0,0/
305	integer :: nmove, nvola
306	integer :: iyear, itimob !obsnypatch
307	integer :: errcnt
308	DATA NMOVE/0/,NVOLA/61/
309
310	! if(ieof(inest).eq.2.and.fdob%nstat.eq.0)then
311	! IF (iprt) print *,'returning from in4dob'
312	! return
313	! endif
314	! IF (iprt) print *,'start in4dob ',inest,xtime
315	IF(nudge_opt.NE.1)RETURN
316
317	! Initialize obs for printout
318	obs_prt = -999
319	newpass = .true.
320	errcnt = 0
321
322	! if start time, or restart time, set obs array to missing value
323	IF(KTAU.EQ.0.OR.KTAU.EQ.KTAUR) THEN
324	DO N=1,NIOBF
325	TIMEOB(N)=99999.
326	ENDDO
327	fdob%xtime_at_rest = xtime !yliu 20080127
328	ENDIF
329	! set number of obs=0 if at start or restart
330	IF(KTAU.EQ.KTAUR)fdob%NSTAT=0
331	NSTA=fdob%NSTAT
332
333	XHOUR=XTIME/60.
334	XHOUR=AMAX1(XHOUR,0.0)
335
336	! DEFINE THE MAX LIMITS OF THE WINDOW
337	TBACK=XHOUR-TWINDO
338	TFORWD=XHOUR+TWINDO
339
340	IF (iprt) then
341	write(msg,fmt='(2(a,f8.3),a,i2)') &
342	'OBS NUDGING: Reading new obs for time window TBACK = ', &
343	tback,' TFORWD = ',tforwd,' for grid = ',inest
344	call wrf_message(msg)
345	ENDIF
346
347	! For obs that have become invalid because they are too old for the current time
348	! window, mark with 99999 to set up for removal from the rolling valid-obs list.
349
350	IF(NSTA.NE.0) THEN
351	NDUM=0
352	t_window : DO N=1,NSTA+1
353	IF((TIMEOB(N)-TBACK).LT.0) THEN
354	TIMEOB(N)=99999.
355	ENDIF
356	IF(TIMEOB(N).LT.9.E4) EXIT t_window
357	NDUM=N
358	ENDDO t_window
359
360	IF (iprt .and. ndum>0) THEN
361	write(msg,fmt='(a,i5,2a)') 'OBS NUDGING: ',ndum,' previously read obs ', &
362	'are now too old for the current window and have been removed.'
363	call wrf_message(msg)
364	ENDIF
365
366	! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
367	! IF (iprt) print *,'ndum at 20=',ndum
368	NDUM=ABS(NDUM)
369	NMOVE=NIOBF-NDUM
370	IF(NMOVE.GT.0 .AND. NDUM.NE.0 ) THEN
371	DO N=1,NMOVE
372	do KN = 1,nndgv
373	VAROBS(KN,N)=VAROBS(KN,N+NDUM)
374	enddo
375	! RIO is the west-east coordinate. RJO is south-north. (ajb)
376	RJO(N)=RJO(N+NDUM)
377	RIO(N)=RIO(N+NDUM)
378	RKO(N)=RKO(N+NDUM)
379	TIMEOB(N)=TIMEOB(N+NDUM)
380	nlevs_ob(n)=nlevs_ob(n+ndum)
381	lev_in_ob(n)=lev_in_ob(n+ndum)
382	plfo(n)=plfo(n+ndum)
383	elevob(n)=elevob(n+ndum)
384	ENDDO
385	ENDIF
386	NOPEN=NMOVE+1
387	IF(NOPEN.LE.NIOBF) THEN
388	DO N=NOPEN,NIOBF
389	do KN = 1,nndgv
390	VAROBS(KN,N)=99999.
391	enddo
392	RIO(N)=99999.
393	RJO(N)=99999.
394	RKO(N)=99999.
395	TIMEOB(N)=99999.
396	nlevs_ob(n)=99999.
397	lev_in_ob(n)=99999.
398	plfo(n)=99999.
399	elevob(n)=99999.
400	ENDDO
401	ENDIF
402	ENDIF
403
404	! Compute map projection info.
405	call set_projection(obs_proj, map_proj, cen_lat, cen_lon, &
406	true_lat1, true_lat2, stand_lon, &
407	known_lat, known_lon, &
408	e_we, e_sn, dxm, dym )
409
410	! FIND THE LAST OBS IN THE LIST
411	NLAST=0
412	last_ob : DO N=1,NIOBF
413	! print *,'nlast,n,timeob(n)=',nlast,n,timeob(n)
414	IF(TIMEOB(N).GT.9.E4) EXIT last_ob
415	NLAST=N
416	ENDDO last_ob
417
418	! print *,'in4dob, after 90 ',nlast,ktau,ktaur,nsta
419	! open file if at beginning or restart
420	IF(KTAU.EQ.0.OR.KTAU.EQ.KTAUR) THEN
421	fdob%RTLAST=-999.
422	INQUIRE (NVOLA+INEST-1,OPENED=OPENED)
423	IF (.NOT. OPENED) THEN
424	ifon(inest)=1
425	write(fonc(1:2),'(i2)')ifon(inest)
426	if(fonc(1:1).eq.' ')fonc(1:1)='0'
427	INQUIRE (file='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2) &
428	,EXIST=exist)
429	if(exist)then
430	IF (iprt) THEN
431	write(msg,*) 'opening first fdda obs file, fonc=', &
432	fonc,' inest=',inest
433	call wrf_message(msg)
434	write(msg,*) 'ifon=',ifon(inest)
435	call wrf_message(msg)
436	ENDIF
437	OPEN(NVOLA+INEST-1, &
438	FILE='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2), &
439	FORM='FORMATTED',STATUS='OLD')
440	else
441	! no first file to open
442	IF (iprt) call wrf_message("there are no fdda obs files to open")
443	return
444	endif
445
446	ENDIF
447	ENDIF !end if(KTAU.EQ.0.OR.KTAU.EQ.KTAUR)
448	! print *,'at jc check1'
449
450	!**********************************************************************
451	! -------------- BIG 100 LOOP OVER N --------------
452	!**********************************************************************
453	! NOW CHECK TO SEE IF EXTRA DATA MUST BE READ IN FROM THE
454	! DATA FILE. CONTINUE READING UNTIL THE REACHING THE EOF
455	! (DATA TIME IS NEGATIVE) OR FIRST TIME PAST TFORWD. THE
456	! LAST OBS CURRENTLY AVAILABLE IS IN N=NMOVE.
457
458	N=NLAST
459	IF(N.EQ.0)GOTO 110
460
461	1001 continue
462
463	! ieof=2 means no more files
464
465	IF(IEOF(inest).GT.1) then
466	GOTO 130
467	endif
468
469	100 CONTINUE
470	!ajb 20070116 bugfix for zero array index. N=0 if first obs is not in the domain.
471	IF(N.ne.0) THEN
472	IF(TIMEOB(N).GT.TFORWD.and.timeob(n).lt.99999.) THEN
473	GOTO 130
474	ENDIF
475	ENDIF
476
477	! OBSFILE: no more data in the obsfile
478	! AJB note: This is where we would implement multi-file reading.
479	if(ieof(inest).eq.1 )then
480	ieof(inest)=2
481	goto 130
482	endif
483
484	!**********************************************************************
485	! -------------- 110 SUBLOOP OVER N --------------
486	!**********************************************************************
487	110 continue
488	! THE TIME OF THE MOST RECENTLY ACQUIRED OBS IS .LE. TFORWD,
489	! SO CONTINUE READING
490
491	IF(N.GT.NIOBF-1)GOTO 120
492	! REPLACE NVOLA WITH LUN 70, AND USE NVOLA AS A FILE COUNTER
493	NVOL=NVOLA+INEST-1
494	IF(fdob%IEODI.EQ.1)GOTO 111
495	read(nvol,101,end=111,err=111)date_char
496	101 FORMAT(1x,a14)
497
498	n=n+1
499
500	! Convert the form of the observation date for geth_idts.
501	call fmt_date(date_char, obs_date)
502
503	! Compute the time period in seconds from the model reference
504	! date (fdob%sdate) until the observation date.
505
506	call geth_idts(obs_date, fdob%sdate(1:19), idts)
507
508	! Convert time in seconds to hours.
509	! In case of restart, correct for new sdate.
510	idts = idts + nint(fdob%xtime_at_rest*60.) ! yliu 20080127
511
512	rtimob =float(idts)/3600.
513	timeob(n)=rtimob
514
515	! print *,'read in ob ',n,timeob(n),rtimob
516	IF(IDYNIN.EQ.1.AND.TIMEOB(N)*60..GT.fdob%DATEND) THEN
517	IF (iprt) THEN
518	write(msg,*) ' IN4DOB: FOR INEST = ',INEST,' AT XTIME = ',XTIME, &
519	' TIMEOB = ',TIMEOB(N)*60.,' AND DATEND = ',fdob%DATEND,' :'
520	call wrf_message(msg)
521	write(msg,*) ' END-OF-DATA FLAG SET FOR OBS-NUDGING', &
522	' DYNAMIC INITIALIZATION'
523	call wrf_message(msg)
524	ENDIF
525	fdob%IEODI=1
526	TIMEOB(N)=99999.
527	rtimob=timeob(n)
528	ENDIF
529	read(nvol,102)latitude,longitude
530	102 FORMAT(2x,2(f9.4,1x))
531
532	! if(ifon.eq.4)print *,'ifon=4',latitude,longitude
533	! this works only for lc projection
534	! yliu: add llxy for all 3 projection
535
536	!ajb Arguments ri and rj have been switched from MM5 orientation.
537
538	CALL latlon_to_ij(obs_proj, latitude, longitude, ri, rj)
539
540	!ajb ri and rj are referenced to the non-staggered grid (not mass-pt!).
541	! (For MM5, they were referenced to the dot grid.)
542
543	ri = ri + .5 !ajb Adjust from mass-pt to non-staggered grid.
544	rj = rj + .5 !ajb Adjust from mass-pt to non-staggered grid.
545
546	rio(n)=ri
547	rjo(n)=rj
548
549	read(nvol,1021)id,namef
550	1021 FORMAT(2x,2(a40,3x))
551	read(nvol,103)platform,source,elevation,is_sound,bogus,meas_count
552	103 FORMAT( 2x,2(a16,2x),f8.0,2x,2(l4,2x),i5)
553
554	! write(6,*) '----- OBS description ----- '
555	! write(6,*) 'platform,source,elevation,is_sound,bogus,meas_count:'
556	! write(6,*) platform,source,elevation,is_sound,bogus,meas_count
557
558	! yliu
559	elevob(n)=elevation
560	! jc
561	! change platform from synop to profiler when needed
562	if(namef(2:9).eq.'PROFILER')platform(7:14)='PROFILER'
563	! yliu
564	if(namef(2:6).eq.'ACARS')platform(7:11)='ACARS'
565	if(namef(1:7).eq.'SATWNDS') platform(1:11)='SATWNDS '
566	if(namef(1:8).eq.'CLASS DA')platform(7:10)='TEMP'
567	! yliu end
568
569	rko(n)=-99.
570	!yliu 20050706
571	! if((platform(7:11).eq.'METAR').or.(platform(7:11).eq.'SPECI').or.
572	! 1 (platform(7:10).eq.'SHIP').or.(platform(7:11).eq.'SYNOP').or.
573	! 1 (platform(1:4).eq.'SAMS'))
574	! 1 rko(n)=1.0
575	if(.NOT. is_sound) rko(n)=1.0
576	!yliu 20050706 end
577
578	! plfo is inFORMATion on what platform. May use this later in adjusting weights
579	plfo(n)=99.
580	if(platform(7:11).eq.'METAR')plfo(n)=1.
581	if(platform(7:11).eq.'SPECI')plfo(n)=2.
582	if(platform(7:10).eq.'SHIP')plfo(n)=3.
583	if(platform(7:11).eq.'SYNOP')plfo(n)=4.
584	if(platform(7:10).eq.'TEMP')plfo(n)=5.
585	if(platform(7:11).eq.'PILOT')plfo(n)=6.
586	if(platform(1:7).eq.'SATWNDS')plfo(n)=7.
587	if(platform(7:11).eq.'SATWI')plfo(n)=7.
588	if(platform(1:4).eq.'SAMS')plfo(n)=8.
589	if(platform(7:14).eq.'PROFILER')plfo(n)=9.
590	! yliu: ACARS->SATWINDS
591	if(platform(7:11).eq.'ACARS')plfo(n)=7.
592	! yliu: end
593	if(plfo(n).eq.99.) then
594	IF (iprt) then
595	write(msg,*) 'n=',n,' unknown ob of type ',platform
596	call wrf_message(msg)
597	ENDIF
598	endif
599
600	!======================================================================
601	!======================================================================
602	! THIS PART READS SOUNDING INFO
603	IF(is_sound)THEN
604	nlevs_ob(n)=real(meas_count)
605	lev_in_ob(n)=1.
606	do imc=1,meas_count
607	! write(6,*) '0 inest = ',inest,' n = ',n
608	! the sounding has one header, many levels. This part puts it into
609	! "individual" observations. There's no other way for nudob to deal
610	! with it.
611	if(imc.gt.1)then ! sub-loop over N
612	n=n+1
613	if(n.gt.niobf)goto 120
614	nlevs_ob(n)=real(meas_count)
615	lev_in_ob(n)=real(imc)
616	timeob(n)=rtimob
617	rio(n)=ri
618	rjo(n)=rj
619	rko(n)=-99.
620	plfo(n)=plfo(n-imc+1)
621	elevob(n)=elevation
622	endif
623
624	read(nvol,104)pressure_data,pressure_qc, &
625	height_data,height_qc, &
626	temperature_data,temperature_qc, &
627	u_met_data,u_met_qc, &
628	v_met_data,v_met_qc, &
629	rh_data,rh_qc
630	104 FORMAT( 1x,6(f11.3,1x,f11.3,1x))
631
632	! yliu: Ensemble - add disturbance to upr obs
633	! if(plfo(n).eq.5.or.plfo(n).eq.6.or.plfo(n).eq.9) then FORE07E08
634	! if(imc.eq.1) then FORE07E08
635	! call srand(n)
636	! t_rand =- (rand(2)-0.5)*6
637	! call srand(n+100000)
638	! u_rand =- (rand(2)-0.5)*6
639	! call srand(n+200000)
640	! v_rand =- (rand(2)-0.5)*6
641	! endif FORE07E08
642	! if(temperature_qc.ge.0..and.temperature_qc.lt.30000..and.
643	! & temperature_data .gt. -88880.0 )
644	! & temperature_data = temperature_data + t_rand
645	! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
646	! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
647	! make sure at least 1 of the components is .ne.0
648	! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
649	! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
650	! u_met_data = u_met_data + u_rand
651	! v_met_data = v_met_data + v_rand
652	! endif
653	! endif FORE07E08
654	! yliu: Ens test - end
655
656
657	! jc
658	! hardwire to switch -777777. qc to 0. here temporarily
659	! -777777. is a sounding level that no qc was done on.
660
661	if(temperature_qc.eq.-777777.)temperature_qc=0.
662	if(pressure_qc.eq.-777777.)pressure_qc=0.
663	if(height_qc.eq.-777777.)height_qc=0.
664	if(u_met_qc.eq.-777777.)u_met_qc=0.
665	if(v_met_qc.eq.-777777.)v_met_qc=0.
666	if(rh_qc.eq.-777777.)rh_qc=0.
667	if(temperature_data.eq.-888888.)temperature_qc=-888888.
668	if(pressure_data.eq.-888888.)pressure_qc=-888888.
669	if(height_data.eq.-888888.)height_qc=-888888.
670	if(u_met_data.eq.-888888.)u_met_qc=-888888.
671	if(v_met_data.eq.-888888.)v_met_qc=-888888.
672	if(rh_data.eq.-888888.)rh_qc=-888888.
673
674	! jc
675	! Hardwire so that only use winds in pilot obs (no winds from temp) and
676	! only use temperatures and rh in temp obs (no temps from pilot obs)
677	! Do this because temp and pilot are treated as 2 platforms, but pilot
678	! has most of the winds, and temp has most of the temps. If use both,
679	! the second will smooth the effect of the first. Usually temps come in after
680	! pilots. pilots usually don't have any temps, but temp obs do have some
681	! winds usually.
682	! plfo=5 is TEMP ob, range sounding is an exception
683	!yliu start -- comment out to test with merged PILOT and TEMP and
684	! do not use obs interpolated by little_r
685	! if(plfo(n).eq.5. .and. namef(1:8).ne.'CLASS DA')then
686	! u_met_data=-888888.
687	! v_met_data=-888888.
688	! u_met_qc=-888888.
689	! v_met_qc=-888888.
690	! endif
691	if(plfo(n).eq.5..and.(u_met_qc.eq.256..or.v_met_qc.eq.256.))then
692	u_met_data=-888888.
693	v_met_data=-888888.
694	u_met_qc=-888888.
695	v_met_qc=-888888.
696	endif
697	!yliu end
698	! plfo=6 is PILOT ob
699	if(plfo(n).eq.6.)then
700	temperature_data=-888888.
701	rh_data=-888888.
702	temperature_qc=-888888.
703	rh_qc=-888888.
704	endif
705
706	!ajb Store temperature for WRF
707	! NOTE: The conversion to potential temperature, performed later in subroutine
708	! errob, requires good pressure data, either directly or via good height data.
709	! So here, in addition to checking for good temperature data, we must also
710	! do a check for good pressure or height.
711	if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
712
713	if( (pressure_qc.ge.0..and.pressure_qc.lt.30000.) .or. &
714	(height_qc .ge.0..and.height_qc .lt.30000.) ) then
715
716	varobs(3,n) = temperature_data
717	else
718	varobs(3,n)=-888888.
719	endif
720
721	else
722	varobs(3,n)=-888888.
723	endif
724
725	!ajb Store obs height
726	if(height_qc.ge.0..and.height_qc.lt.30000.)then
727	varobs(6,n)=height_data
728	else
729	varobs(6,n)=-888888.
730	endif
731
732	if(pressure_qc.ge.0..and.pressure_qc.lt.30000.)then
733	! if(pressure_qc.ge.0.)then
734	varobs(5,n)=pressure_data
735	else
736	varobs(5,n)=-888888.
737	IF (iprt) THEN
738	if(varobs(6,n).eq.-888888.000) then
739	if (errcnt.le.10) then
740	write(msg,) '** PROBLEM: sounding, p and ht undefined',latitude,longitude
741	call wrf_message(msg)
742	errcnt = errcnt + 1
743	if (errcnt.gt.10) call wrf_message("MAX of 10 warnings issued.")
744	endif
745	endif
746	ENDIF
747	endif
748	if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
749	! don't use data above 80 mb
750	if((varobs(5,n).gt.0.).and.(varobs(5,n).le.8.))then
751	u_met_data=-888888.
752	v_met_data=-888888.
753	u_met_qc=-888888.
754	v_met_qc=-888888.
755	temperature_data=-888888.
756	temperature_qc=-888888.
757	rh_data=-888888.
758	rh_qc=-888888.
759	endif
760
761
762	! Store horizontal wind components for WRF
763	if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
764	(v_met_qc.ge.0..and.v_met_qc.lt.30000.).and. &
765	! make sure at least 1 of the components is .ne.0
766	(u_met_data.ne.0..or.v_met_data.ne.0.))then
767
768	! If Earth-relative wind vector, need to rotate it to grid-relative coords
769	if(u_met_qc.eq.129. .and. v_met_qc.eq.129.) then
770	CALL rotate_vector(longitude,u_met_data,v_met_data, &
771	obs_proj,map_proj)
772	endif
773	varobs(1,n)=u_met_data
774	varobs(2,n)=v_met_data
775	else
776	varobs(1,n)=-888888.
777	varobs(2,n)=-888888.
778	endif
779
780	r_data=-888888.
781
782	if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
783	if((pressure_qc.ge.0.).and.(temperature_qc.ge.0.).and. &
784	(pressure_qc.lt.30000.).and.(temperature_qc.lt.30000.))then
785	call rh2r(rh_data,temperature_data,pressure_data*.01, &
786	r_data,0) ! yliu, change last arg from 1 to 0
787	else
788	! print *,'rh, but no t or p to convert',temperature_qc, &
789	! pressure_qc,n
790	r_data=-888888.
791	endif
792	endif
793	varobs(4,n)=r_data
794	enddo ! end do imc=1,meas_count
795	! print *,'--- sdng n=',n,nlevs_ob(n),lev_in_ob(n),timeob(n)
796	! read in non-sounding obs
797
798	ELSEIF(.NOT.is_sound)THEN
799	nlevs_ob(n)=1.
800	lev_in_ob(n)=1.
801	read(nvol,105)slp_data,slp_qc, &
802	ref_pres_data,ref_pres_qc, &
803	height_data,height_qc, &
804	temperature_data,temperature_qc, &
805	u_met_data,u_met_qc, &
806	v_met_data,v_met_qc, &
807	rh_data,rh_qc, &
808	psfc_data,psfc_qc, &
809	precip_data,precip_qc
810	105 FORMAT( 1x,9(f11.3,1x,f11.3,1x))
811
812	! Ensemble: add disturbance to sfc obs
813	! call srand(n)
814	! t_rand =+ (rand(2)-0.5)*5
815	! call srand(n+100000)
816	! u_rand =+ (rand(2)-0.5)*5
817	! call srand(n+200000)
818	! v_rand =+ (rand(2)-0.5)*5
819	! if(temperature_qc.ge.0..and.temperature_qc.lt.30000. .and.
820	! & temperature_data .gt. -88880.0 )
821	! & temperature_data = temperature_data + t_rand
822	! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
823	! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
824	! make sure at least 1 of the components is .ne.0
825	! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
826	! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
827	! u_met_data = u_met_data + u_rand
828	! v_met_data = v_met_data + v_rand
829	! endif
830	! yliu: Ens test - end
831
832	!Lilis
833
834	! calculate psfc if slp is there
835	if((psfc_qc.lt.0.).and.(slp_qc.ge.0..and.slp_qc.lt.30000.).and. &
836	(temperature_qc.ge.0..and.temperature_qc.lt.30000.).and. &
837	(slp_data.gt.90000.))then
838	tbar=temperature_data+0.5elevation.0065
839	psfc_data=slp_dataexp(-elevation/(rovgtbar))
840	varobs(5,n)=psfc_data
841	psfc_qc=0.
842	endif
843
844	!c No Very rough estimate of psfc from sfc elevation if UUtah ob and elev>1000m
845	! estimate psfc from temp and elevation
846	! Do not know sfc pressure in model at this point.
847	! if((psfc_qc.lt.0.).and.(elevation.gt.1000.).and.
848	! 1 (temperature_qc.ge.0..and.temperature_qc.lt.30000.)
849	! 1 .and.(platform(7:16).eq.'SYNOP PRET'))then
850	if((psfc_qc.lt.0.).and. &
851	(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
852	tbar=temperature_data+0.5elevation.0065
853	psfc_data=100000.exp(-elevation/(rovgtbar))
854	varobs(5,n)=psfc_data
855	psfc_qc=0.
856	endif
857
858	if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and.(psfc_data.gt.70000. &
859	.and.psfc_data.lt.105000.))then
860	varobs(5,n)=psfc_data
861	else
862	varobs(5,n)=-888888.
863	endif
864
865	if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
866
867	!Lilie
868	!ajb Store temperature for WRF
869	if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
870
871	if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and. &
872	(psfc_data.gt.70000. .and.psfc_data.lt.105000.))then
873
874	varobs(3,n) = temperature_data
875	else
876	varobs(3,n)=-888888.
877	endif
878	else
879	varobs(3,n)=-888888.
880	endif
881
882	! Store horizontal wind components for WRF
883	if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
884	(v_met_qc.ge.0..and.v_met_qc.lt.30000.).and. &
885	! make sure at least 1 of the components is .ne.0
886	(u_met_data.ne.0..or.v_met_data.ne.0.))then
887
888	! If Earth-relative wind vector, need to rotate it to grid-relative coords
889	if(u_met_qc.eq.129. .and. v_met_qc.eq.129.) then
890	CALL rotate_vector(longitude,u_met_data,v_met_data, &
891	obs_proj,map_proj)
892	endif
893	varobs(1,n)=u_met_data
894	varobs(2,n)=v_met_data
895	else
896	varobs(1,n)=-888888.
897	varobs(2,n)=-888888.
898	endif
899
900	! jc
901	! if a ship ob has rh<70%, then throw out
902
903	if(plfo(n).eq.3..and.rh_qc.ge.0..and.rh_data.lt.70.)then
904	rh_qc=-888888.
905	rh_data=-888888.
906	endif
907	!
908	r_data=-888888.
909	if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
910	if((psfc_qc.ge.0..and.psfc_qc.lt.30000.) &
911	.and.(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
912	! rh_data=amin1(rh_data,96.) ! yliu: do not allow surface to be saturated
913	call rh2r(rh_data,temperature_data,psfc_data*.01, &
914	r_data,0) ! yliu, change last arg from 1 to 0
915	else
916	! print *,'rh, but no t or p',temperature_data,
917	! 1 psfc_data,n
918	r_data=-888888.
919	endif
920	endif
921	varobs(4,n)=r_data
922	ELSE
923	IF (iprt) THEN
924	call wrf_message(" ====== ")
925	call wrf_message(" NO Data Found ")
926	ENDIF
927	ENDIF !end if(is_sound)
928	! END OF SFC OBS INPUT SECTION
929	!======================================================================
930	!======================================================================
931	! check if ob time is too early (only applies to beginning)
932	IF(RTIMOB.LT.TBACK-TWINDO)then
933	IF (iprt) call wrf_message("ob too early")
934	n=n-1
935	GOTO 110
936	ENDIF
937
938	! check if this ob is a duplicate
939	! this check has to be before other checks
940	njend=n-1
941	if(is_sound)njend=n-meas_count
942	do njc=1,njend
943	! Check that time, lat, lon, and platform all match exactly.
944	! Platforms 1-4 (surface obs) can match with each other. Otherwise,
945	! platforms have to match exactly.
946	if( (timeob(n).eq.timeob(njc)) .and. &
947	(rio(n).eq.rio(njc)) .and. &
948	(rjo(n).eq.rjo(njc)) .and. &
949	(plfo(njc).ne.99.) ) then
950	!yliu: if two sfc obs are departed less than 1km, consider they are redundant
951	! (abs(rio(n)-rio(njc))*dscg.gt.1000.) &
952	! .or. (abs(rjo(n)-rjo(njc))*dscg.gt.1000.) &
953	! .or. (plfo(njc).eq.99.) )goto 801
954	!yliu end
955	! If platforms different, and either > 4, jump out
956	if( ( (plfo(n).le.4.).and.(plfo(njc).le.4.) ) .or. &
957	(plfo(n).eq.plfo(njc)) ) then
958
959	! if not a sounding, and levels are the same then replace first occurrence
960	if((.not.is_sound).and.(rko(njc).eq.rko(n))) then
961	! print *,'dup single ob-replace ',n,inest,
962	! plfo(n),plfo(njc)
963	! this is the sfc ob replacement part
964	do KN = 1,nndgv
965	VAROBS(KN,njc)=VAROBS(KN,n)
966	enddo
967	! don't need to switch these because they're the same
968	! RIO(njc)=RIO(n)
969	! RJO(njc)=RJO(n)
970	! RKO(njc)=RKO(n)
971	! TIMEOB(njc)=TIMEOB(n)
972	! nlevs_ob(njc)=nlevs_ob(n)
973	! lev_in_ob(njc)=lev_in_ob(n)
974	! plfo(njc)=plfo(n)
975	! end sfc ob replacement part
976
977	n=n-1
978	goto 100
979	! It's harder to fix the soundings, since the number of levels may be different
980	! The easiest thing to do is to just set the first occurrence to all missing, and
981	! keep the second occurrence, or vice versa.
982	! For temp or profiler keep the second, for pilot keep the one with more levs
983	! This is for a temp or prof sounding, equal to same
984	! also if a pilot, but second one has more obs
985	elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
986	( (plfo(njc).eq.5.).or.(plfo(njc).eq.9.).or. &
987	( (plfo(njc).eq.6.).and. &
988	(nlevs_ob(n).ge.nlevs_ob(njc)) ) ) )then
989	IF (iprt) THEN
990	write(msg,*) 'duplicate sounding - eliminate first occurrence', &
991	n,inest,meas_count,nlevs_ob(njc), &
992	latitude,longitude,plfo(njc)
993	call wrf_message(msg)
994	ENDIF
995	if(lev_in_ob(njc).ne.1.) then
996	IF (iprt) THEN
997	write(msg,) 'problem ****** - dup sndg ', &
998	lev_in_ob(njc),nlevs_ob(njc)
999	call wrf_message(msg)
1000	ENDIF
1001	endif
1002	! n=n-meas_count
1003	! set the first sounding ob to missing
1004	do njcc=njc,njc+nint(nlevs_ob(njc))-1
1005	do KN = 1,nndgv
1006	VAROBS(KN,njcc)=-888888.
1007	enddo
1008	plfo(njcc)=99.
1009	enddo
1010	goto 100
1011	! if a pilot, but first one has more obs
1012	elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
1013	(plfo(njc).eq.6.).and. &
1014	(nlevs_ob(n).lt.nlevs_ob(njc)) )then
1015	IF (iprt) THEN
1016	write(msg,*) &
1017	'duplicate pilot sounding - eliminate second occurrence', &
1018	n,inest,meas_count,nlevs_ob(njc), &
1019	latitude,longitude,plfo(njc)
1020	call wrf_message(msg)
1021	ENDIF
1022	if(lev_in_ob(njc).ne.1.) then
1023	IF (iprt) THEN
1024	write(msg,) 'problem ****** - dup sndg ', &
1025	lev_in_ob(njc),nlevs_ob(njc)
1026	call wrf_message(msg)
1027	ENDIF
1028	endif
1029	n=n-meas_count
1030
1031	!ajb Reset timeob for discarded indices.
1032	do imc = n+1, n+meas_count
1033	timeob(imc) = 99999.
1034	enddo
1035	goto 100
1036	! This is for a single-level satellite upper air ob - replace first
1037	elseif( (is_sound).and. &
1038	(nlevs_ob(njc).eq.1.).and. &
1039	(nlevs_ob(n).eq.1.).and. &
1040	(varobs(5,njc).eq.varobs(5,n)).and. &
1041	(plfo(njc).eq.7.).and.(plfo(n).eq.7.) ) then
1042	IF (iprt) then
1043	write(msg,*) &
1044	'duplicate single lev sat-wind ob - replace first',n, &
1045	inest,meas_count,varobs(5,n)
1046	call wrf_message(msg)
1047	ENDIF
1048	! this is the single ua ob replacement part
1049	do KN = 1,nndgv
1050	VAROBS(KN,njc)=VAROBS(KN,n)
1051	enddo
1052	! don't need to switch these because they're the same
1053	! RIO(njc)=RIO(n)
1054	! RJO(njc)=RJO(n)
1055	! RKO(njc)=RKO(n)
1056	! TIMEOB(njc)=TIMEOB(n)
1057	! nlevs_ob(njc)=nlevs_ob(n)
1058	! lev_in_ob(njc)=lev_in_ob(n)
1059	! plfo(njc)=plfo(n)
1060	! end single ua ob replacement part
1061	n=n-1
1062	goto 100
1063	else
1064	! IF (iprt) THEN
1065	! write(msg,*) 'duplicate location, but no match otherwise',n,njc, &
1066	! plfo(n),varobs(5,n),nlevs_ob(n),lev_in_ob(n), &
1067	! plfo(njc),varobs(5,njc),nlevs_ob(njc),lev_in_ob(njc)
1068	! call wrf_message(msg)
1069	! ENDIF
1070	endif
1071	endif
1072	endif
1073	! end of njc do loop
1074	enddo
1075
1076	! check if ob is a sams ob that came in via UUtah - discard
1077	if( plfo(n).eq.4..and.(platform(7:16).eq.'SYNOP PRET').and. &
1078	(id(7:15).eq.'METNET= 3') )then
1079	! print *,'elim metnet=3',latitude,longitude,rtimob
1080	n=n-1
1081	goto 100
1082	endif
1083
1084	! check if ob is in the domain
1085	if( (ri.lt.2.).or.(ri.gt.real(e_we-1)).or.(rj.lt.2.).or. &
1086	(rj.gt.real(e_sn-1)) ) then
1087
1088	n=n-meas_count
1089	!ajb Reset timeob for discarded indices.
1090	do imc = n+1, n+meas_count
1091	timeob(imc) = 99999.
1092	enddo
1093	goto 100
1094	endif
1095
1096	IF(TIMEOB(N).LT.fdob%RTLAST) THEN
1097	IF (iprt) THEN
1098	call wrf_message("2 OBS ARE NOT IN CHRONOLOGICAL ORDER")
1099	call wrf_message("NEW YEAR?")
1100	write(msg,*) 'timeob,rtlast,n=',timeob(n),fdob%rtlast,n
1101	call wrf_message(msg)
1102	ENDIF
1103	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 111' )
1104	ELSE
1105	fdob%RTLAST=TIMEOB(N)
1106	ENDIF
1107	! Save obs and model latitude and longitude for printout
1108	CALL collect_obs_info(newpass,inest,n,latitude,longitude, &
1109	nlast,nprev,niobf,id,stnid_prt, &
1110	rio,rjo,prt_max,prt_freq,xlat,xlong, &
1111	obs_prt,lat_prt,lon_prt,mlat_prt,mlon_prt, &
1112	e_we,e_sn,ims,ime,jms,jme,its,ite,jts,jte)
1113	GOTO 100
1114	111 CONTINUE
1115	!**********************************************************************
1116	! -------------- END BIG 100 LOOP OVER N --------------
1117	!**********************************************************************
1118
1119	if (iprt) then
1120	write(msg,5403) NVOL,XTIME
1121	call wrf_message(msg)
1122	endif
1123	IEOF(inest)=1
1124
1125	close(NVOLA+INEST-1)
1126	IF (iprt) then
1127	write(msg,*) 'closed fdda file for inest=',inest,nsta
1128	call wrf_message(msg)
1129	ENDIF
1130
1131	! AJB note: Go back and check for more files. (Multi-file implementation)
1132	goto 1001
1133
1134	120 CONTINUE
1135	! THE OBSERVATION ARRAYS ARE FULL AND THE MOST RECENTLY
1136	! ACQUIRED OBS STILL HAS TIMEOB .LE. TFORWD. SO START
1137	! DECREASING THE SIZE OF THE WINDOW
1138	! get here if too many obs
1139	IF (iprt) THEN
1140	write(msg,121) N,NIOBF
1141	call wrf_message(msg)
1142	ENDIF
1143	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 122' )
1144
1145	130 CONTINUE
1146	! READ CYCLE IS COMPLETED. DETERMINE THE NUMBER OF OBS IN
1147	! THE CURRENT WINDOW
1148	!
1149	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1150	! BUT FIRST, WHEN KTAU.EQ.0 (OR IN GENERAL, KTAUR), DISCARD THE
1151	! "OLD" OBS FIRST...
1152
1153	! Get here if at end of file, or if obs time is beyond what we need right now.
1154	! On startup, we report the index of the last obs read.
1155	! For restarts, we need to remove any old obs and then repack obs list.
1156	IF(KTAU.EQ.KTAUR)THEN
1157	NSTA=0
1158	keep_obs : DO N=1,NIOBF
1159	! try to keep all obs, but just don't use yet
1160	! (don't want to throw away last obs read in - especially if
1161	! its a sounding, in which case it looks like many obs)
1162	IF(TIMEOB(N).GT.9.e4) EXIT keep_obs
1163	if(timeob(n).gt.tforwd) then
1164	if(iprt) then
1165	write(msg,950) inest
1166	call wrf_message(msg)
1167	write(msg,951) n,timeob(n),tforwd
1168	call wrf_message(msg)
1169	endif
1170	950 FORMAT('Saving index of first ob after end of current time window ', &
1171	'for nest = ', i3,':')
1172	951 FORMAT(' ob index = ',i8,', time of ob = ',f8.4, &
1173	' hrs, end of time window = ',f8.4,' hrs')
1174	endif
1175	NSTA=N
1176	ENDDO keep_obs
1177
1178	NDUM=0
1179	! make time=99999. if ob is too old
1180	! print *,'tback,nsta=',tback,nsta
1181	old_obs : DO N=1,NSTA+1
1182	IF((TIMEOB(N)-TBACK).LT.0)THEN
1183	TIMEOB(N)=99999.
1184	ENDIF
1185	! print *,'n,ndum,timeob=',n,ndum,timeob(n)
1186	IF(TIMEOB(N).LT.9.E4) EXIT old_obs
1187	NDUM=N
1188	ENDDO old_obs
1189
1190	! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
1191	IF (iprt .and. ktaur > 0) THEN
1192	write(msg,fmt='(a,i5,a)') 'OBS NUDGING: Upon restart, skipped over ',ndum, &
1193	' obs that are now too old for the current obs window.'
1194	call wrf_message(msg)
1195	ENDIF
1196
1197	NDUM=ABS(NDUM)
1198	NMOVE=NIOBF-NDUM
1199	IF( NMOVE.GT.0 .AND. NDUM.NE.0) THEN
1200	DO N=1,NMOVE
1201	do KN = 1,nndgv
1202	VAROBS(KN,N)=VAROBS(KN,N+NDUM)
1203	enddo
1204	RJO(N)=RJO(N+NDUM)
1205	RIO(N)=RIO(N+NDUM)
1206	RKO(N)=RKO(N+NDUM)
1207	TIMEOB(N)=TIMEOB(N+NDUM)
1208	nlevs_ob(n)=nlevs_ob(n+ndum)
1209	lev_in_ob(n)=lev_in_ob(n+ndum)
1210	plfo(n)=plfo(n+ndum)
1211	ENDDO
1212	ENDIF
1213	! moved obs up. now fill remaining space with 99999.
1214	NOPEN=NMOVE+1
1215	IF(NOPEN.LE.NIOBF) THEN
1216	DO N=NOPEN,NIOBF
1217	do KN = 1,nndgv
1218	VAROBS(KN,N)=99999.
1219	enddo
1220	RIO(N)=99999.
1221	RJO(N)=99999.
1222	RKO(N)=99999.
1223	TIMEOB(N)=99999.
1224	ENDDO
1225	ENDIF
1226	ENDIF
1227	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1228	NSTA=0
1229	! print *,'nsta at restart setting is ',nsta
1230	! recalculate nsta after moving things around
1231	recalc : DO N=1,NIOBF
1232	! try to save all obs - don't throw away latest read in
1233	IF(TIMEOB(N).GT.9.e4) EXIT recalc
1234	NSTA=N
1235	! nsta=n-1 ! yliu test
1236	ENDDO recalc
1237
1238	! Find the number of stations that are actually within the time window.
1239	nstaw = nvals_le_limit(nsta, timeob, tforwd)
1240
1241	IF (iprt) then
1242	write(msg,160) KTAU,XTIME,NSTAW
1243	call wrf_message(msg)
1244	ENDIF
1245	IF(KTAU.EQ.KTAUR)THEN
1246	IF(nudge_opt.EQ.1)THEN
1247	TWDOP=TWINDO*60.
1248	IF (iprt) THEN
1249	write(msg,1449) INEST,RINXY,RINSIG,TWDOP
1250	call wrf_message(msg)
1251	IF(ISWIND.EQ.1) then
1252	write(msg,1450) GIV
1253	call wrf_message(msg)
1254	ELSE
1255	write(msg,1455) INEST
1256	call wrf_message("")
1257	call wrf_message(msg)
1258	call wrf_message("")
1259	ENDIF
1260	IF(ISTEMP.EQ.1) then
1261	write(msg,1451) GIT
1262	call wrf_message(msg)
1263	ELSE
1264	write(msg,1456) INEST
1265	call wrf_message("")
1266	call wrf_message(msg)
1267	ENDIF
1268	IF(ISMOIS.EQ.1) then
1269	call wrf_message("")
1270	write(msg,1452) GIQ
1271	call wrf_message(msg)
1272	ELSE
1273	write(msg,1457) INEST
1274	call wrf_message("")
1275	call wrf_message(msg)
1276	call wrf_message("")
1277	ENDIF
1278	ENDIF
1279	ENDIF
1280	ENDIF
1281	IF(KTAU.EQ.KTAUR)THEN
1282	IF(fdob%IWTSIG.NE.1)THEN
1283	IF (iprt) THEN
1284	write(msg,555)
1285	call wrf_message(msg)
1286	write(msg,556) fdob%RINFMNRINXY,fdob%RINFMXRINXY,fdob%PFREE*10.
1287	call wrf_message(msg)
1288	ENDIF
1289	IF(fdob%RINFMN.GT.fdob%RINFMX) then
1290	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 556' )
1291	ENDIF
1292	! IS MINIMUM GREATER THAN MAXIMUM?
1293
1294	IF (iprt) then
1295	write(msg,557) fdob%DPSMX*10.,fdob%DCON
1296	call wrf_message(msg)
1297	ENDIF
1298	IF(fdob%DPSMX.GT.10.) then
1299	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 557' )
1300	ENDIF
1301	ENDIF
1302	ENDIF
1303
1304	IF(KTAU.EQ.KTAUR)THEN
1305	IF (iprt) then
1306	write(msg,601) INEST,IONF
1307	call wrf_message(msg)
1308	call wrf_message("")
1309	ENDIF
1310	ENDIF
1311	fdob%NSTAT=NSTA
1312	fdob%NSTAW=NSTAW
1313
1314	555 FORMAT(1X,' ABOVE THE SURFACE LAYER, OBS NUDGING IS PERFORMED', &
1315	' ON PRESSURE LEVELS,')
1316	556 FORMAT(1X,' WHERE RINXY VARIES LINEARLY FROM ',E11.3,' KM AT', &
1317	' THE SURFACE TO ',E11.3,' KM AT ',F7.2,' MB AND ABOVE')
1318	557 FORMAT(1X,' IN THE SURFACE LAYER, WXY IS A FUNCTION OF ', &
1319	'DPSMX = ',F7.2,' MB WITH DCON = ',E11.3, &
1320	' - SEE SUBROUTINE NUDOB')
1321	601 FORMAT('FOR EFFICIENCY, THE OBS NUDGING FREQUENCY ', &
1322	'FOR MESH #',I2,' IS ',1I2,' CGM TIMESTEPS ')
1323	121 FORMAT(' WARNING: NOBS = ',I4,' IS GREATER THAN NIOBF = ', &
1324	I4,': INCREASE PARAMETER NIOBF')
1325	5403 FORMAT(1H0,'-------------EOF REACHED FOR NVOL = ',I3, &
1326	' AND XTIME = ',F10.2,'-------------------')
1327	160 FORMAT('****** CALL IN4DOB AT KTAU = ',I5,' AND XTIME = ', &
1328	F10.2,': NSTA = ',I7,' ******')
1329	1449 FORMAT('*****NUDGING INDIVIDUAL OBS ON MESH #',I2, &
1330	' WITH RINXY = ', &
1331	E11.3,' KM, RINSIG = ',E11.3,' AND TWINDO (HALF-PERIOD) = ', &
1332	E11.3,' MIN')
1333	1450 FORMAT(1X,'NUDGING IND. OBS WINDS WITH GIV = ',E11.3)
1334	1451 FORMAT(1X,'NUDGING IND. OBS TEMPERATURE WITH GIT = ',E11.3)
1335	1452 FORMAT(1X,'NUDGING IND. OBS MOISTURE WITH GIQ = ',E11.3)
1336	1455 FORMAT(1X,'*** OBS WIND NUDGING FOR MESH ',I2,' IS TURNED OFF!!')
1337	1456 FORMAT(1X,'*** OBS TEMPERATURE NUDGING FOR MESH ',I2,' IS TURNED OFF!!')
1338	1457 FORMAT(1X,'*** OBS MOISTURE NUDGING FOR MESH ',I2,' IS TURNED OFF!!')
1339
1340	RETURN
1341	END SUBROUTINE in4dob
1342
1343	SUBROUTINE julgmt(mdate,julgmtn,timanl,julday,gmt,ind)
1344	! CONVERT MDATE YYMMDDHH TO JULGMT (JULIAN DAY * 100. +GMT)
1345	! AND TO TIMANL (TIME IN MINUTES WITH RESPECT TO MODEL TIME)
1346	! IF IND=0 INPUT MDATE, OUTPUT JULGMTN AND TIMANL
1347	! IF IND=1 INPUT TIMANL, OUTPUT JULGMTN
1348	! IF IND=2 INPUT JULGMTN, OUTPUT TIMANL
1349	INTEGER, intent(in) :: MDATE
1350	REAL, intent(out) :: JULGMTN
1351	REAL, intent(out) :: TIMANL
1352	INTEGER, intent(in) :: JULDAY
1353	REAL, intent(in) :: GMT
1354	INTEGER, intent(in) :: IND
1355
1356	!*** DECLARATIONS FOR IMPLICIT NONE
1357	real :: MO(12), rjulanl, houranl, rhr
1358
1359	integer :: iyr, idate1, imo, idy, ihr, my1, my2, my3, ileap
1360	integer :: juldayn, juldanl, idymax, mm
1361
1362
1363	IF(IND.EQ.2)GOTO 150
1364	IYR=INT(MDATE/1000000.+0.001)
1365	IDATE1=MDATE-IYR*1000000
1366	IMO=INT(IDATE1/10000.+0.001)
1367	IDY=INT((IDATE1-IMO*10000.)/100.+0.001)
1368	IHR=IDATE1-IMO10000-IDY100
1369	MO(1)=31
1370	MO(2)=28
1371	! IS THE YEAR A LEAP YEAR? (IN THIS CENTURY)
1372	IYR=IYR+1900
1373	MY1=MOD(IYR,4)
1374	MY2=MOD(IYR,100)
1375	MY3=MOD(IYR,400)
1376	ILEAP=0
1377	! jc
1378	! IF(MY1.EQ.0.AND.MY2.NE.0.OR.MY3.EQ.0)THEN
1379	IF(MY1.EQ.0)THEN
1380	ILEAP=1
1381	MO(2)=29
1382	ENDIF
1383	IF(IND.EQ.1)GOTO 200
1384	MO(3)=31
1385	MO(4)=30
1386	MO(5)=31
1387	MO(6)=30
1388	MO(7)=31
1389	MO(8)=31
1390	MO(9)=30
1391	MO(10)=31
1392	MO(11)=30
1393	MO(12)=31
1394	JULDAYN=0
1395	DO 100 MM=1,IMO-1
1396	JULDAYN=JULDAYN+MO(MM)
1397	100 CONTINUE
1398
1399	IF(IHR.GE.24)THEN
1400	IDY=IDY+1
1401	IHR=IHR-24
1402	ENDIF
1403	JULGMTN=(JULDAYN+IDY)*100.+IHR
1404	! CONVERT JULGMT TO TIMANL WRT MODEL TIME IN MINUTES (XTIME)
1405	150 CONTINUE
1406	JULDANL=INT(JULGMTN/100.+0.000001)
1407	RJULANL=FLOAT(JULDANL)*100.
1408	HOURANL=JULGMTN-RJULANL
1409	TIMANL=(FLOAT(JULDANL-JULDAY)24.-GMT+HOURANL)60.
1410	RETURN
1411	200 CONTINUE
1412	RHR=GMT+TIMANL/60.+0.000001
1413	IDY=JULDAY
1414	IDYMAX=365+ILEAP
1415	300 IF(RHR.GE.24.0)THEN
1416	RHR=RHR-24.0
1417	IDY=IDY+1
1418	GOTO 300
1419	ENDIF
1420	IF(IDY.GT.IDYMAX)IDY=IDY-IDYMAX
1421	JULGMTN=FLOAT(IDY)*100.+RHR
1422	RETURN
1423	END SUBROUTINE julgmt
1424
1425	SUBROUTINE rh2r(rh,t,p,r,iice)
1426
1427	! convert rh to r
1428	! if iice=1, use saturation with respect to ice
1429	! rh is 0-100.
1430	! r is g/g
1431	! t is K
1432	! p is mb
1433	!
1434	REAL, intent(in) :: rh
1435	REAL, intent(in) :: t
1436	REAL, intent(in) :: p
1437	REAL, intent(out) :: r
1438	INTEGER, intent(in) :: iice
1439
1440	!*** DECLARATIONS FOR IMPLICIT NONE
1441	real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
1442	real esat, rsat
1443
1444	eps=0.62197
1445	e0=6.1078
1446	eslcon1=17.2693882
1447	eslcon2=35.86
1448	esicon1=21.8745584
1449	esicon2=7.66
1450	t0=260.
1451
1452	! print *,'rh2r input=',rh,t,p
1453	rh1=rh*.01
1454
1455	if(iice.eq.1.and.t.le.t0)then
1456	esat=e0exp(esicon1(t-273.16)/(t-esicon2))
1457	else
1458	esat=e0exp(eslcon1(t-273.16)/(t-eslcon2))
1459	endif
1460	rsat=eps*esat/(p-esat)
1461	! print *,'rsat,esat=',rsat,esat
1462	r=rh1*rsat
1463
1464	! print *,'rh2r rh,t,p,r=',rh1,t,p,r
1465
1466	return
1467	END SUBROUTINE rh2r
1468
1469	SUBROUTINE rh2rb(rh,t,p,r,iice)
1470
1471	! convert rh to r
1472	! if iice=1, use daturation with respect to ice
1473	! rh is 0-100.
1474	! r is g/g
1475	! t is K
1476	! p is mb
1477
1478	REAL, intent(in) :: rh
1479	REAL, intent(in) :: t
1480	REAL, intent(in) :: p
1481	REAL, intent(out) :: r
1482	INTEGER, intent(in) :: iice
1483
1484	!*** DECLARATIONS FOR IMPLICIT NONE
1485	real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
1486	real esat, rsat
1487	character(len=200) :: msg ! Argument to wrf_message
1488
1489	eps=0.622
1490	e0=6.112
1491	eslcon1=17.67
1492	eslcon2=29.65
1493	esicon1=22.514
1494	esicon2=6.15e3
1495	t0=273.15
1496
1497	write(msg,*) 'rh2r input=',rh,t,p
1498	call wrf_message(msg)
1499	rh1=rh*.01
1500
1501	if(iice.eq.1.and.t.le.t0)then
1502	esat=e0*exp(esicon1-esicon2/t)
1503	else
1504	esat=e0exp(eslcon1(t-t0)/(t-eslcon2))
1505	endif
1506	rsat=eps*esat/(p-esat)
1507	! print *,'rsat,esat=',rsat,esat
1508	r=rh1epsrsat/(eps+rsat*(1.-rh1))
1509
1510	write(msg,*) 'rh2r rh,t,p,r=',rh1,t,p,r
1511	call wrf_message(msg)
1512	rh1=rh*.01
1513
1514	return
1515	END SUBROUTINE rh2rb
1516
1517	SUBROUTINE set_projection (obs_proj, map_proj, cen_lat, cen_lon, &
1518	true_lat1, true_lat2, stand_lon, &
1519	known_lat, known_lon, &
1520	e_we, e_sn, dxm, dym )
1521
1522	USE module_llxy
1523
1524	!*************************************************************************
1525	! Purpose: Set map projection information which will be used to map the
1526	! observation (lat,lon) location to its corresponding (x,y)
1527	! location on the WRF (coarse) grid. using the selected map
1528	! projection (e.g., Lambert, Mercator, Polar Stereo, etc).
1529	!*************************************************************************
1530
1531	IMPLICIT NONE
1532
1533	TYPE(PROJ_INFO), intent(out) :: obs_proj ! structure for obs projection info.
1534	INTEGER, intent(in) :: map_proj ! map projection index
1535	REAL, intent(in) :: cen_lat ! center latitude for map projection
1536	REAL, intent(in) :: cen_lon ! center longiture for map projection
1537	REAL, intent(in) :: true_lat1 ! truelat1 for map projection
1538	REAL, intent(in) :: true_lat2 ! truelat2 for map projection
1539	REAL, intent(in) :: stand_lon ! standard longitude for map projection
1540	INTEGER, intent(in) :: e_we ! max grid index in south-north coordinate
1541	INTEGER, intent(in) :: e_sn ! max grid index in west-east coordinate
1542	REAL, intent(in) :: known_lat ! latitude of domain origin point (i,j)=(1,1)
1543	REAL, intent(in) :: known_lon ! longigude of domain origin point (i,j)=(1,1)
1544	REAL, intent(in) :: dxm ! grid size in x (meters)
1545	REAL, intent(in) :: dym ! grid size in y (meters)
1546
1547	! Set up map transformation structure
1548	CALL map_init(obs_proj)
1549
1550	! Mercator
1551	IF (map_proj == PROJ_MERC) THEN
1552	CALL map_set(PROJ_MERC, obs_proj, &
1553	truelat1 = true_lat1, &
1554	lat1 = known_lat, &
1555	lon1 = known_lon, &
1556	knowni = 1., &
1557	knownj = 1., &
1558	dx = dxm)
1559
1560	! Lambert conformal
1561	ELSE IF (map_proj == PROJ_LC) THEN
1562	CALL map_set(PROJ_LC, obs_proj, &
1563	truelat1 = true_lat1, &
1564	truelat2 = true_lat2, &
1565	stdlon = stand_lon, &
1566	lat1 = known_lat, &
1567	lon1 = known_lon, &
1568	knowni = 1., &
1569	knownj = 1., &
1570	dx = dxm)
1571
1572	! Polar stereographic
1573	ELSE IF (map_proj == PROJ_PS) THEN
1574	CALL map_set(PROJ_PS, obs_proj, &
1575	truelat1 = true_lat1, &
1576	stdlon = stand_lon, &
1577	lat1 = known_lat, &
1578	lon1 = known_lon, &
1579	knowni = 1., &
1580	knownj = 1., &
1581	dx = dxm)
1582	! Cassini (global ARW)
1583	ELSE IF (map_proj == PROJ_CASSINI) THEN
1584	CALL map_set(PROJ_CASSINI, obs_proj, &
1585	latinc = dym360.0/(2.0EARTH_RADIUS_M*PI), &
1586	loninc = dxm360.0/(2.0EARTH_RADIUS_M*PI), &
1587	lat1 = known_lat, &
1588	lon1 = known_lon, &
1589	! We still need to get POLE_LAT and POLE_LON metadata variables before
1590	! this will work for rotated poles.
1591	lat0 = 90.0, &
1592	lon0 = 0.0, &
1593	knowni = 1., &
1594	knownj = 1., &
1595	stdlon = stand_lon)
1596
1597	! Rotated latitude-longitude
1598	ELSE IF (map_proj == PROJ_ROTLL) THEN
1599	CALL map_set(PROJ_ROTLL, obs_proj, &
1600	! I have no idea how this should work for NMM nested domains
1601	ixdim = e_we-1, &
1602	jydim = e_sn-1, &
1603	phi = real(e_sn-2)*dym/2.0, &
1604	lambda = real(e_we-2)*dxm, &
1605	lat1 = cen_lat, &
1606	lon1 = cen_lon, &
1607	latinc = dym, &
1608	loninc = dxm, &
1609	stagger = HH)
1610
1611	END IF
1612
1613	END SUBROUTINE set_projection
1614
1615	SUBROUTINE fmt_date(idate,odate) !obsnypatch
1616
1617	!*************************************************************************
1618	! Purpose: Re-format a character date string from YYYYMMDDHHmmss form
1619	! to YYYY-MM-DD_HH:mm:ss form.
1620	! INPUT:
1621	! IDATE - Date string as YYYYMMDDHHmmss
1622	! OUTPUT:
1623	! ODATE - Date string as YYYY-MM-DD_HH:mm:ss
1624	!*************************************************************************
1625
1626	IMPLICIT NONE
1627
1628	CHARACTER*14, intent(in) :: idate ! input date string
1629	CHARACTER*19, intent(out) :: odate ! output date string
1630
1631	odate(1:19) = "0000-00-00_00:00:00"
1632	odate(1:4) = idate(1:4) ! Year
1633	odate(6:7) = idate(5:6) ! Month
1634	odate(9:10) = idate(7:8) ! Day
1635	odate(12:13) = idate(9:10) ! Hours
1636	odate(15:16) = idate(11:12) ! Minutes
1637	odate(18:19) = idate(13:14) ! Seconds
1638
1639	RETURN
1640	END SUBROUTINE fmt_date
1641
1642	INTEGER FUNCTION nvals_le_limit(isize, values, limit)
1643	!------------------------------------------------------------------------------
1644	! PURPOSE: Return the number of values in a (real) non-decreasing array that
1645	! are less than or equal to the specified upper limit.
1646	! NOTE: It is important that the array is non-decreasing!
1647	!
1648	!------------------------------------------------------------------------------
1649	IMPLICIT NONE
1650
1651	INTEGER, INTENT(IN) :: isize ! Size of input array
1652	REAL, INTENT(IN) :: values(isize) ! Input array of reals
1653	REAL, INTENT(IN) :: limit ! Upper limit
1654
1655	! Local variables
1656	integer :: n
1657
1658	! Search the array from largest to smallest values.
1659	find_nvals: DO n = isize, 1, -1
1660	if(values(n).le.limit) EXIT find_nvals
1661	ENDDO find_nvals
1662	nvals_le_limit = n
1663
1664	RETURN
1665	END FUNCTION nvals_le_limit
1666
1667	SUBROUTINE collect_obs_info(newpass,inest,n,latitude,longitude, &
1668	nlast,nprev,niobf,station_id,stnid, &
1669	rio,rjo,prt_max,prt_freq,xlat,xlong, &
1670	obs, lat,lon, mlat,mlon, &
1671	e_we,e_sn,ims,ime,jms,jme,its,ite,jts,jte)
1672	!*************************************************************************
1673	! Purpose: Collect the obs index, obs latitude, obs longitude, obs station
1674	! id, and model latitude and longitude values for print
1675	! diagnostics. Note that THIS SUBROUTINE IS CALLED INTERATIVELY
1676	! FROM IN4DOB, WITHIN THE OBS READ LOOP that reads new obser-
1677	! vations needed for the new time window. Flag newpass is true
1678	! the first time collect_obs_info is called from the read-loop
1679	! for a new time window. So for each pass of IN4DOB, newpass is
1680	! true the first time IN4DOB calls collec_obs_info.
1681
1682	! OBS (soundings) contain multiple obs levels. So on each sub-
1683	! sequent call of collect_obs_info for a specific pass of IN4DOB,
1684	! n will jump by the number of levels in the sounding.
1685	!
1686	! Here, nlast refers to the index of the last valid-time obs
1687	! that was read in during the last pass of IN4DOB (after the old
1688	! obs were removed). This way we can properly start storing
1689	! obs information for the new obs that are being read on this
1690	! pass of IN4DOB, beginning with the first newly read obs for
1691	! this pass of IN4DOB.
1692	!
1693	! Note that nprev is needed to properly handle soundings. On
1694	! each pass, n is stored into nprev, and on each subsequent
1695	! pass of collect_obs_info, a loop is performed from nprev+1 to n.
1696	!*************************************************************************
1697
1698	IMPLICIT NONE
1699
1700	LOGICAL, intent(inout) :: newpass ! New pass flag
1701	INTEGER, intent(in) :: inest ! nest index
1702	INTEGER, intent(in) :: n ! Observation index
1703	REAL, intent(in) :: latitude ! Latitude of obs
1704	REAL, intent(in) :: longitude ! Latitude of obs
1705	INTEGER, intent(in) :: nlast ! Last obs of valid obs, prev window
1706	INTEGER, intent(inout) :: nprev ! Previous obs in new window read seq
1707	INTEGER, intent(in) :: niobf ! Maximum number of observations
1708	CHARACTER*15, intent(in) :: station_id ! First 15 chars of station id for obs n
1709	INTEGER, intent(in) :: prt_max ! Max no. of obs for diagnostic printout
1710	INTEGER, intent(inout) :: stnid(40,prt_max) ! Station ids for diagnostic printout
1711	REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
1712	REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
1713	INTEGER, intent(in) :: prt_freq ! Frequency for diagnostic printout
1714	REAL, DIMENSION( ims:ime, jms:jme ), &
1715	intent(in ) :: xlat, xlong ! Lat/lon on mass-pt grid
1716	INTEGER, intent(inout) :: obs(prt_max) ! Obs index for printout
1717	REAL, intent(inout) :: lat(prt_max) ! Obs latitude for printout
1718	REAL, intent(inout) :: lon(prt_max) ! Obs longitude for printout
1719	REAL, intent(inout) :: mlat(prt_max) ! Model latitude at (rio,rjo) for printout
1720	REAL, intent(inout) :: mlon(prt_max) ! Model longitude at (rio,rjo) for printout
1721	INTEGER, intent(in) :: e_we ! Max grid index in south-north
1722	INTEGER, intent(in) :: e_sn ! Max grid index in west-east
1723	INTEGER, intent(in) :: ims ! Grid mem start (west-east)
1724	INTEGER, intent(in) :: ime ! Grid mem end (west-east)
1725	INTEGER, intent(in) :: jms ! Grid mem start (south-north)
1726	INTEGER, intent(in) :: jme ! Grid mem end (south-north)
1727	INTEGER, intent(in) :: its ! Grid tile start (west-east)
1728	INTEGER, intent(in) :: ite ! Grid tile end (west-east)
1729	INTEGER, intent(in) :: jts ! Grid tile start (south-north)
1730	INTEGER, intent(in) :: jte ! Grid tile end (south-north)
1731
1732	! Local variables
1733	integer i ! Loop counter over station id character
1734	integer nn ! Loop counter over obs index
1735	integer ndx,ndxp ! Index into printout arrays (ndx and prev ndx)
1736	real :: ri, rj ! Mass-pt coord of obs
1737	integer :: ril, rjl ! Mass-pt integer coord immed sw of obs
1738	integer :: iend, jend ! Upper i, j index for interpolation
1739	real :: dxob, dyob ! Grid fractions for interp
1740	logical :: llsave ! Save lat/lon values if true
1741	character(len=200) :: msg ! Argument to wrf_message
1742
1743	if(newpass) then
1744	newpass = .false.
1745	nprev = nlast ! Reset in case old obs have been discarded from prev window
1746	endif
1747
1748	! Start iteration only if we have not yet stored prt_max number of obs for printing.
1749	! Note: The loop below could represent multiple levels in a sounding, so we
1750	! go ahead and start the loop if the beginning index (ndx) is prt_max or
1751	! less, and then exit the loop if ndx exceeds prt_max.
1752	if(prt_freq.gt.0) then
1753	ndx = (n-nlast-1)/prt_freq + 1
1754	ndxp = (nprev-nlast-1)/prt_freq + 1
1755	else
1756	write(msg,*) 'STOP! OBS NAMELIST INPUT obs_prt_freq MUST BE GREATER THAN ZERO.'
1757	call wrf_message(msg)
1758	write(msg,*) 'THE NAMELIST VALUE IS',prt_freq,' FOR NEST ',inest
1759	call wrf_message(msg)
1760	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP' )
1761	endif
1762
1763	! write(6,'5(a,i5),a,a15') 'n = ',n,' nlast = ',nlast,' ndx = ',ndx, &
1764	! ' nprev = ',nprev,' ndxp = ',ndxp, &
1765	! ' station id = ',station_id
1766
1767	if(ndxp .lt. prt_max) then
1768
1769	MODCHK : do nn = nprev+1, n
1770	llsave = .false.
1771
1772	! if( mod(nn-1,prt_freq) .eq. 0 ) then
1773	if( mod(nn-nlast-1,prt_freq) .eq. 0 ) then
1774	ndx = (nn-nlast-1)/prt_freq + 1
1775	if(ndx.gt.prt_max) EXIT MODCHK ! Limit printout to prt_max entries
1776	llsave = .true.
1777	endif
1778	if(llsave) then
1779
1780	! Collect obs index and latitude and longitude.
1781	obs(ndx) = nn
1782	lat(ndx) = latitude
1783	lon(ndx) = longitude
1784
1785	! Collect first 15 chars of obs station id (in integer format).
1786	do i = 1,15
1787	stnid(i,ndx) = ichar(station_id(i:i))
1788	enddo
1789
1790	! Compute and collect the model latitude and longitude at the obs point.
1791	CALL get_model_latlon(nn,niobf,rio,rjo,xlat,xlong,e_we,e_sn, &
1792	ims,ime,jms,jme,its,ite,jts,jte, &
1793	mlat(ndx),mlon(ndx))
1794	endif !end if(llsave)
1795	enddo MODCHK
1796
1797	endif !end if(ndx .le. prt_max)
1798
1799	! Save index of previous obs in read loop.
1800	nprev = n
1801
1802	END SUBROUTINE collect_obs_info
1803
1804	SUBROUTINE get_model_latlon(n,niobf,rio,rjo,xlat,xlong,e_we,e_sn, &
1805	ims,ime,jms,jme,its,ite,jts,jte, &
1806	mlat,mlon)
1807	!*************************************************************************
1808	! Purpose: Use bilinear interpolation to compute the model latitude and
1809	! longitude at the observation point.
1810	!*************************************************************************
1811
1812	IMPLICIT NONE
1813
1814	INTEGER, intent(in) :: n ! Observation index
1815	INTEGER, intent(in) :: niobf ! Maximum number of observations
1816	REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
1817	REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
1818	REAL, DIMENSION( ims:ime, jms:jme ), &
1819	intent(in ) :: xlat, xlong ! Lat/lon on mass-pt grid
1820	INTEGER, intent(in) :: e_we ! Max grid index in south-north
1821	INTEGER, intent(in) :: e_sn ! Max grid index in west-east
1822	INTEGER, intent(in) :: ims ! Grid mem start (west-east)
1823	INTEGER, intent(in) :: ime ! Grid mem end (west-east)
1824	INTEGER, intent(in) :: jms ! Grid mem start (south-north)
1825	INTEGER, intent(in) :: jme ! Grid mem end (south-north)
1826	INTEGER, intent(in) :: its ! Grid tile start (west-east)
1827	INTEGER, intent(in) :: ite ! Grid tile end (west-east)
1828	INTEGER, intent(in) :: jts ! Grid tile start (south-north)
1829	INTEGER, intent(in) :: jte ! Grid tile end (south-north)
1830	REAL, intent(out) :: mlat ! Model latitude at obs point
1831	REAL, intent(out) :: mlon ! Model longitude at obs point
1832
1833	! Local variables
1834	integer ndx ! Index into save arrays
1835	real :: ri, rj ! Mass-pt coord of obs
1836	integer :: ril, rjl ! Mass-pt integer coord immed sw of obs
1837	integer :: iend, jend ! Upper i, j index for interpolation
1838	real :: dxob, dyob ! Grid fractions for interp
1839
1840	! Compute model latitude and longitude if point on tile.
1841	ri = rio(n) - .5 ! mass-pt west-east obs grid coord
1842	rj = rjo(n) - .5 ! mass-pt south-north obs grid coord
1843	ril = int(ri)
1844	rjl = int(rj)
1845	dxob = ri - float(ril)
1846	dyob = rj - float(rjl)
1847	iend = min(ite+1,e_we-2)
1848	jend = min(jte+1,e_sn-2)
1849	mlat = -999
1850	mlon = -999
1851
1852	if(ri.ge.its .and. ri.lt.iend .and. rj.ge.jts .and. rj.lt.jend) then
1853
1854	! bilinear interpolation
1855	mlat = ((1.-dyob)((1.-dxob)xlat(ril,rjl)+ &
1856	dxob*xlat(ril+1,rjl) &
1857	)+dyob((1.-dxob)xlat(ril,rjl+1)+ &
1858	dxob*xlat(ril+1,rjl+1)))
1859
1860	mlon = ((1.-dyob)((1.-dxob)xlong(ril,rjl)+ &
1861	dxob*xlong(ril+1,rjl) &
1862	)+dyob((1.-dxob)xlong(ril,rjl+1)+ &
1863	dxob*xlong(ril+1,rjl+1)))
1864	endif
1865
1866	END SUBROUTINE get_model_latlon
1867
1868	SUBROUTINE rotate_vector(lon,u,v,obs_proj,map_proj)
1869
1870	USE module_llxy
1871
1872	!*************************************************************************
1873	! Purpose: Rotate a single Earth-relative wind vector to a grid-relative
1874	! wind vector.
1875	!*************************************************************************
1876
1877	IMPLICIT NONE
1878
1879	REAL, intent(in) :: lon ! Longitude (deg)
1880	REAL, intent(inout) :: u ! U-component of wind vector
1881	REAL, intent(inout) :: v ! V-component of wind vector
1882	TYPE(PROJ_INFO),intent(in) :: obs_proj ! Structure for obs projection
1883	INTEGER, intent(in) :: map_proj ! Map projection index
1884
1885	! Local variables
1886	real diff, alpha
1887	double precision udbl, vdbl
1888
1889	! Only rotate winds for Lambert conformal or polar stereographic
1890	if (map_proj == PROJ_LC .or. map_proj == PROJ_PS) then
1891
1892	diff = obs_proj%stdlon - lon
1893	if (diff > 180.) then
1894	diff = diff - 360.
1895	else if (diff < -180.) then
1896	diff = diff + 360.
1897	end if
1898
1899	! Calculate the rotation angle, alpha, in radians
1900	if (map_proj == PROJ_LC) then
1901	alpha = diff * obs_proj%cone * rad_per_deg * obs_proj%hemi
1902	else
1903	alpha = diff * rad_per_deg * obs_proj%hemi
1904	end if
1905
1906	udbl = vsin(alpha) + ucos(alpha)
1907	vdbl = vcos(alpha) - usin(alpha)
1908	u = udbl
1909	v = vdbl
1910
1911	endif
1912	END SUBROUTINE rotate_vector
1913
1914	#endif
1915	!-----------------------------------------------------------------------
1916	! End subroutines for in4dob
1917	!-----------------------------------------------------------------------

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

Download in other formats:

Original Format