Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

wrf_fddaobs_in.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 61.5 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
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 = dtmin*grid%itimestep
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	CALL in4dob(inest, xtime, ktau, krest, dtmin, grid%julday, grid%gmt, &
84	grid%obs_nudge_opt, grid%obs_nudge_wind, grid%obs_nudge_temp, &
85	grid%obs_nudge_mois, grid%obs_nudge_pstr, grid%obs_coef_wind, &
86	grid%obs_coef_temp, grid%obs_coef_mois, grid%obs_coef_pstr, &
87	grid%obs_rinxy, grid%obs_rinsig, grid%fdob%window, &
88	grid%obs_npfi, grid%obs_ionf, grid%obs_nobs_prt, &
89	grid%obs_idynin, &
90	grid%obs_dtramp, grid%fdob, grid%fdob%varobs, &
91	grid%fdob%timeob, grid%fdob%nlevs_ob, grid%fdob%lev_in_ob, &
92	grid%fdob%plfo, grid%fdob%elevob, grid%fdob%rio, &
93	grid%fdob%rjo, grid%fdob%rko, &
94	config_flags%cen_lat, &
95	config_flags%cen_lon, &
96	config_flags%stand_lon, &
97	config_flags%truelat1, config_flags%truelat2, &
98	grid%fdob%known_lat, grid%fdob%known_lon, &
99	config_flags%dx, config_flags%dy, rovg, t0, &
100	ide, jde, &
101	grid%fdob%sn_maxcg, grid%fdob%we_maxcg, config_flags%map_proj, &
102	model_config_rec%parent_grid_ratio, &
103	model_config_rec%i_parent_start(inest), &
104	model_config_rec%j_parent_start(inest), &
105	model_config_rec%nobs_ndg_vars, grid%max_obs, iprt_in4dob)
106	ENDIF
107
108	RETURN
109	#endif
110	END SUBROUTINE wrf_fddaobs_in
111	#if ( EM_CORE == 1 )
112	!------------------------------------------------------------------------------
113	! Begin subroutine in4dob and its subroutines
114	!------------------------------------------------------------------------------
115	SUBROUTINE in4dob(inest, xtime, ktau, ktaur, dtmin, julday, gmt, &
116	nudge_opt, iswind, istemp, &
117	ismois, ispstr, giv, &
118	git, giq, gip, &
119	rinxy, rinsig, twindo, &
120	npfi, ionf, nobs_prt, idynin, &
121	dtramp, fdob, varobs, &
122	timeob, nlevs_ob, lev_in_ob, &
123	plfo, elevob, rio, &
124	rjo, rko, &
125	cen_lat, &
126	cen_lon, &
127	stand_lon, &
128	true_lat1, true_lat2, &
129	known_lat, known_lon, &
130	dxm, dym, rovg, t0, e_we, e_sn, &
131	sn_maxcg, we_maxcg, map_proj, &
132	parent_grid_ratio, &
133	i_parent_start, &
134	j_parent_start, &
135	nndgv, niobf, iprt)
136
137	USE module_domain
138	USE module_model_constants, ONLY : rcp
139	USE module_llxy
140
141	IMPLICIT NONE
142
143	! THIS IS SUBROUTINE READS AN OBSERVATION DATA FILE AND
144	! SELECTS ONLY THOSE VALUES OBSERVED AT TIMES THAT FALL
145	! WITHIN A TIME WINDOW (TWINDO) CENTERED ABOUT THE CURRENT
146	! FORECAST TIME (XTIME). THE INCOMING OBS FILES MUST BE
147	! IN CHRONOLOGICAL ORDER.
148	!
149	! NOTE: This routine was originally designed for MM5, which uses
150	! a nonstandard (I,J) coordinate system. For WRF, I is the
151	! east-west running coordinate, and J is the south-north
152	! running coordinate. So "J-slab" here is west-east in
153	! extent, not south-north as for MM5. RIO and RJO have
154	! the opposite orientation here as for MM5. -ajb 06/10/2004
155
156	! NOTE - IN4DOB IS CALLED ONLY FOR THE COARSE MESH TIMES IN4DOB.10
157	! - VAROBS(IVAR,N) HOLDS THE OBSERVATIONS. IN4DOB.11
158	! IVAR=1 UOBS IN4DOB.12
159	! IVAR=2 VOBS IN4DOB.13
160	! IVAR=3 TOBS IN4DOB.14
161	! IVAR=4 QOBS IN4DOB.15
162	! IVAR=5 PSOBS (CROSS) IN4DOB.16
163
164	INTEGER, intent(in) :: niobf ! maximum number of observations
165	INTEGER, intent(in) :: nndgv ! number of nudge variables
166	INTEGER, intent(in) :: INEST ! nest level
167	REAL, intent(in) :: xtime ! model time in minutes
168	INTEGER, intent(in) :: KTAU ! current timestep
169	INTEGER, intent(in) :: KTAUR ! restart timestep
170	REAL, intent(in) :: dtmin ! dt in minutes
171	INTEGER, intent(in) :: julday ! Julian day
172	REAL, intent(in) :: gmt ! Greenwich Mean Time
173	INTEGER, intent(in) :: nudge_opt ! obs-nudge flag for this nest
174	INTEGER, intent(in) :: iswind ! nudge flag for wind
175	INTEGER, intent(in) :: istemp ! nudge flag for temperature
176	INTEGER, intent(in) :: ismois ! nudge flag for moisture
177	INTEGER, intent(in) :: ispstr ! nudge flag for pressure
178	REAL, intent(in) :: giv ! coefficient for wind
179	REAL, intent(in) :: git ! coefficient for temperature
180	REAL, intent(in) :: giq ! coefficient for moisture
181	REAL, intent(in) :: gip ! coefficient for pressure
182	REAL, intent(in) :: rinxy ! horizontal radius of influence (km)
183	REAL, intent(in) :: rinsig ! vertical radius of influence (on sigma)
184	REAL, intent(in) :: twindo ! (time window)/2 (min) for nudging
185	INTEGER, intent(in) :: npfi ! coarse-grid time-step frequency for diagnostics
186	INTEGER, intent(in) :: ionf ! coarse-grid time-step frequency for obs-nudging calcs
187	INTEGER, intent(in) :: nobs_prt ! Number of current obs to print grid information for.
188	INTEGER, intent(in) :: idynin ! for dynamic initialization using a ramp-down function
189	REAL, intent(in) :: dtramp ! time period in minutes for ramping
190	TYPE(fdob_type), intent(inout) :: fdob ! derived data type for obs data
191	REAL, intent(inout) :: varobs(nndgv,niobf) ! observational values in each variable
192	REAL, intent(inout) :: timeob(niobf) ! model times for each observation (hours)
193	REAL, intent(inout) :: nlevs_ob(niobf) ! numbers of levels in sounding obs
194	REAL, intent(inout) :: lev_in_ob(niobf) ! level in sounding-type obs
195	REAL, intent(inout) :: plfo(niobf) ! index for type of obs-platform
196	REAL, intent(inout) :: elevob(niobf) ! elevations of observations (meters)
197	REAL, intent(inout) :: rio(niobf) ! west-east grid coordinate (non-staggered grid)
198	REAL, intent(inout) :: rjo(niobf) ! south-north grid coordinate (non-staggered grid)
199	REAL, intent(inout) :: rko(niobf) ! vertical grid coordinate
200	REAL, intent(in) :: cen_lat ! center latitude for map projection
201	REAL, intent(in) :: cen_lon ! center longiture for map projection
202	REAL, intent(in) :: stand_lon ! standard longitude for map projection
203	REAL, intent(in) :: true_lat1 ! truelat1 for map projection
204	REAL, intent(in) :: true_lat2 ! truelat2 for map projection
205	REAL, intent(in) :: known_lat ! latitude of domain origin point (i,j)=(1,1)
206	REAL, intent(in) :: known_lon ! longigude of domain origin point (i,j)=(1,1)
207	REAL, intent(in) :: dxm ! grid size in x (meters)
208	REAL, intent(in) :: dym ! grid size in y (meters)
209	REAL, intent(in) :: rovg ! constant rho over g
210	REAL, intent(in) :: t0 ! background temperature
211	INTEGER, intent(in) :: e_we ! max grid index in south-north coordinate
212	INTEGER, intent(in) :: e_sn ! max grid index in west-east coordinate
213	INTEGER, intent(in) :: sn_maxcg ! maximum coarse grid south-north coordinate
214	INTEGER, intent(in) :: we_maxcg ! maximum coarse grid west-east coordinate
215	INTEGER, intent(in) :: map_proj ! map projection index
216	INTEGER, intent(in) :: parent_grid_ratio ! parent to nest grid ration
217	INTEGER, intent(in) :: i_parent_start ! starting i coordinate in parent domain
218	INTEGER, intent(in) :: j_parent_start ! starting j coordinate in parent domain
219	LOGICAL, intent(in) :: iprt ! print flag
220
221	!*** DECLARATIONS FOR IMPLICIT NONE
222	integer :: n, ndum, nopen, nlast, nvol, idate, imm, iss
223	integer :: nsta ! number of stations held in timeobs array
224	integer :: nstaw ! # stations within the actual time window
225	integer :: ips ! # stations to report printout
226	integer :: meas_count, imc, njend, njc, njcc, julob
227	real :: hourob, rjulob
228	real :: xhour, tback, tforwd, rjdate1, timanl1, rtimob
229	real :: rj, ri, elevation, pressure_data
230	real :: pressure_qc, height_data, height_qc, temperature_data
231	real :: temperature_qc, u_met_data, u_met_qc, v_met_data
232	real :: v_met_qc, rh_data, rh_qc, r_data, slp_data, slp_qc
233	real :: ref_pres_data, ref_pres_qc, psfc_data, psfc_qc
234	real :: precip_data, precip_qc, tbar, twdop
235	real*8 :: tempob
236	INTEGER, EXTERNAL :: nvals_le_limit ! for finding #obs with timeobs <= tforwd
237
238	! Local variables
239	TYPE (PROJ_INFO) :: obs_proj ! Structure for obs projection information.
240	character*14 date_char
241	character*40 platform,source,id,namef
242	character*2 fonc
243	real latitude,longitude
244	real lat_prt(100),lon_prt(100)
245	logical is_sound,bogus
246	LOGICAL OPENED,exist
247	integer :: ieof(5),ifon(5)
248	data ieof/0,0,0,0,0/
249	data ifon/0,0,0,0,0/
250	integer :: nmove, nvola
251	DATA NMOVE/0/,NVOLA/61/
252
253	if(ieof(inest).eq.2.and.fdob%nstat.eq.0)then
254	IF (iprt) print *,'returning from in4dob'
255	return
256	endif
257	IF (iprt) print *,'start in4dob ',inest,xtime
258	IF(nudge_opt.NE.1)RETURN
259
260	! if start time, or restart time, set obs array to missing value
261	IF(KTAU.EQ.0.OR.KTAU.EQ.KTAUR) THEN
262	DO N=1,NIOBF
263	TIMEOB(N)=99999.
264	ENDDO
265	ENDIF
266	! set number of obs=0 if at start or restart
267	IF(KTAU.EQ.KTAUR)fdob%NSTAT=0
268	NSTA=fdob%NSTAT
269	XHOUR=(XTIME-DTMIN)/60.
270	XHOUR=AMAX1(XHOUR,0.0)
271
272	10 CONTINUE
273
274	! DEFINE THE MAX LIMITS OF THE WINDOW
275	TBACK=XHOUR-fdob%WINDOW
276	TFORWD=XHOUR+fdob%WINDOW
277
278	if (iprt) write(6,*) 'TBACK = ',tback,' TFORWD = ',tforwd
279
280	IF(NSTA.NE.0) THEN
281	NDUM=0
282	t_window : DO N=1,NSTA+1
283	IF((TIMEOB(N)-TBACK).LT.0) THEN
284	TIMEOB(N)=99999.
285	ENDIF
286	IF(TIMEOB(N).LT.9.E4) EXIT t_window
287	NDUM=N
288	ENDDO t_window
289
290	! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
291	IF (iprt) print *,'ndum at 20=',ndum
292	NDUM=ABS(NDUM)
293	NMOVE=NIOBF-NDUM
294	IF(NMOVE.GT.0 .AND. NDUM.NE.0 ) THEN
295	DO N=1,NMOVE
296	VAROBS(1,N)=VAROBS(1,N+NDUM)
297	VAROBS(2,N)=VAROBS(2,N+NDUM)
298	VAROBS(3,N)=VAROBS(3,N+NDUM)
299	VAROBS(4,N)=VAROBS(4,N+NDUM)
300	VAROBS(5,N)=VAROBS(5,N+NDUM)
301	! RIO is the west-east coordinate. RJO is south-north. (ajb)
302	RJO(N)=RJO(N+NDUM)
303	RIO(N)=RIO(N+NDUM)
304	RKO(N)=RKO(N+NDUM)
305	TIMEOB(N)=TIMEOB(N+NDUM)
306	nlevs_ob(n)=nlevs_ob(n+ndum)
307	lev_in_ob(n)=lev_in_ob(n+ndum)
308	plfo(n)=plfo(n+ndum)
309	elevob(n)=elevob(n+ndum)
310	ENDDO
311	ENDIF
312	NOPEN=NMOVE+1
313	IF(NOPEN.LE.NIOBF) THEN
314	DO N=NOPEN,NIOBF
315	VAROBS(1,N)=99999.
316	VAROBS(2,N)=99999.
317	VAROBS(3,N)=99999.
318	VAROBS(4,N)=99999.
319	VAROBS(5,N)=99999.
320	RIO(N)=99999.
321	RJO(N)=99999.
322	RKO(N)=99999.
323	TIMEOB(N)=99999.
324	nlevs_ob(n)=99999.
325	lev_in_ob(n)=99999.
326	plfo(n)=99999.
327	elevob(n)=99999.
328	ENDDO
329	ENDIF
330	ENDIF
331
332	! Compute map projection info.
333	call set_projection(obs_proj, map_proj, cen_lat, cen_lon, &
334	true_lat1, true_lat2, stand_lon, &
335	known_lat, known_lon, &
336	e_we, e_sn, dxm, dym )
337
338	! FIND THE LAST OBS IN THE LIST
339	NLAST=0
340	last_ob : DO N=1,NIOBF
341	! print *,'nlast,n,timeob(n)=',nlast,n,timeob(n)
342	IF(TIMEOB(N).GT.9.E4) EXIT last_ob
343	NLAST=N
344	ENDDO last_ob
345
346	! print *,'in4dob, after 90 ',nlast,ktau,ktaur,nsta
347	! open file if at beginning or restart
348	IF(KTAU.EQ.0.OR.KTAU.EQ.KTAUR) THEN
349	fdob%RTLAST=-999.
350	INQUIRE (NVOLA+INEST-1,OPENED=OPENED)
351	IF (.NOT. OPENED) THEN
352	ifon(inest)=1
353	write(fonc(1:2),'(i2)')ifon(inest)
354	if(fonc(1:1).eq.' ')fonc(1:1)='0'
355	INQUIRE (file='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2) &
356	,EXIST=exist)
357	if(exist)then
358	IF (iprt) THEN
359	print *,'opening first fdda obs file, fonc=', &
360	fonc,' inest=',inest
361	print *,'ifon=',ifon(inest)
362	ENDIF
363	OPEN(NVOLA+INEST-1, &
364	FILE='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2), &
365	FORM='FORMATTED',STATUS='OLD')
366	else
367	! no first file to open
368	IF (iprt) print *,'there are no fdda obs files to open'
369	return
370	endif
371
372	ENDIF
373	ENDIF !end if(KTAU.EQ.0.OR.KTAU.EQ.KTAUR)
374	! print *,'at jc check1'
375
376	!**********************************************************************
377	! -------------- BIG 100 LOOP OVER N --------------
378	!**********************************************************************
379	! NOW CHECK TO SEE IF EXTRA DATA MUST BE READ IN FROM THE
380	! DATA FILE. CONTINUE READING UNTIL THE REACHING THE EOF
381	! (DATA TIME IS NEGATIVE) OR FIRST TIME PAST TFORWD. THE
382	! LAST OBS CURRENTLY AVAILABLE IS IN N=NMOVE.
383	N=NLAST
384	IF(N.EQ.0)GOTO 110
385
386	1001 continue
387
388	! ieof=2 means no more files
389	! print *,'after 1001,n,timeob(n)=',n,timeob(n)
390
391	IF(IEOF(inest).GT.1) then
392	GOTO 130
393	endif
394
395	100 CONTINUE
396	!ajb 20070116 bugfix for situation that first obs is not in the domain
397	IF(N.ne.0) THEN
398	IF(TIMEOB(N).GT.TFORWD.and.timeob(n).lt.99999.) THEN
399	GOTO 130
400	ENDIF
401	ENDIF
402
403	! OBSFILE: no more data in the obsfile
404	if(ieof(inest).eq.1 )then
405	ieof(inest)=2
406	goto 130
407	endif
408
409	!**********************************************************************
410	! -------------- 110 SUBLOOP OVER N --------------
411	!**********************************************************************
412	! THE TIME OF THE MOST RECENTLY ACQUIRED OBS IS .LE. TFORWD,
413	! SO CONTINUE READING
414	110 continue
415	IF(N.GT.NIOBF-1)GOTO 120
416	! REPLACE NVOLA WITH LUN 70, AND USE NVOLA AS A FILE COUNTER
417	NVOL=NVOLA+INEST-1
418	IF(fdob%IEODI.EQ.1)GOTO 111
419	read(nvol,101,end=111,err=111)date_char
420	101 FORMAT(1x,a14)
421
422	n=n+1
423
424	read(date_char(3:10),'(i8)')idate
425	read(date_char(11:12),'(i2)')imm
426	read(date_char(13:14),'(i2)')iss
427	! output is rjdate (jjjhh.) and timanl (time in minutes since model start)
428	call julgmt(idate,rjdate1,timanl1,julday,gmt,0)
429	rtimob=rjdate1+real(imm)/60.+real(iss)/3600.
430	timeob(n)=rtimob
431	timeob(n) = int(timeob(n)*1000)/1000.0
432
433	! CONVERT TIMEOB FROM JULIAN DATE AND GMT FORM TO FORECAST
434	! TIME IN HOURS (EX. TIMEOB=13002.4 REPRESENTS JULDAY 130
435	! AND GMT (HOUR) = 2.4)
436	JULOB=TIMEOB(N)/100.+0.000001
437	RJULOB=FLOAT(JULOB)*100.
438	tempob = (timeob(n)*1000.)
439	tempob = int(tempob)
440	tempob = tempob/1000.
441	timeob(n) = tempob
442	HOUROB=TIMEOB(N)-RJULOB
443	TIMEOB(N)=FLOAT(JULOB-JULDAY)*24.-GMT+HOUROB
444	rtimob=timeob(n)
445
446	! print *,'read in ob ',n,timeob(n),rtimob
447	IF(IDYNIN.EQ.1.AND.TIMEOB(N)*60..GT.fdob%DATEND) THEN
448	IF (iprt) THEN
449	PRINT*,' IN4DOB: FOR INEST = ',INEST,' AT XTIME = ',XTIME, &
450	' TIMEOB = ',TIMEOB(N)*60.,' AND DATEND = ',fdob%DATEND,' :'
451	PRINT*,' END-OF-DATA FLAG SET FOR OBS-NUDGING', &
452	' DYNAMIC INITIALIZATION'
453	ENDIF
454	fdob%IEODI=1
455	TIMEOB(N)=99999.
456	rtimob=timeob(n)
457	ENDIF
458	read(nvol,102)latitude,longitude
459	! save lat and long for printout
460	if(n.le.100) then
461	lat_prt(n) = latitude
462	lon_prt(n) = longitude
463	endif
464	102 FORMAT(2x,2(f7.2,3x))
465
466	! if(ifon.eq.4)print *,'ifon=4',latitude,longitude
467	! this works only for lc projection
468	! yliu: add llxy for all 3 projection
469
470	!ajb Arguments ri and rj have been switched from MM5 orientation.
471
472	CALL latlon_to_ij(obs_proj, latitude, longitude, ri, rj)
473
474	!ajb ri and rj are referenced to the non-staggered grid (not mass-pt!).
475	! (For MM5, they were referenced to the dot grid.)
476
477	ri = ri + .5 !ajb Adjust from mass-pt to non-staggered grid.
478	rj = rj + .5 !ajb Adjust from mass-pt to non-staggered grid.
479
480	rio(n)=ri
481	rjo(n)=rj
482
483	read(nvol,1021)id,namef
484	1021 FORMAT(2x,2(a40,3x))
485	read(nvol,103)platform,source,elevation,is_sound,bogus,meas_count
486	103 FORMAT( 2x,2(a16,2x),f8.0,2x,2(l4,2x),i5)
487
488	! write(6,*) '----- OBS description ----- '
489	! write(6,*) 'platform,source,elevation,is_sound,bogus,meas_count:'
490	! write(6,*) platform,source,elevation,is_sound,bogus,meas_count
491
492	! yliu
493	elevob(n)=elevation
494	! jc
495	! change platform from synop to profiler when needed
496	if(namef(2:9).eq.'PROFILER')platform(7:14)='PROFILER'
497	! yliu
498	if(namef(2:6).eq.'ACARS')platform(7:11)='ACARS'
499	if(namef(1:7).eq.'SATWNDS') platform(1:11)='SATWNDS '
500	if(namef(1:8).eq.'CLASS DA')platform(7:10)='TEMP'
501	! yliu end
502
503	rko(n)=-99.
504	!yliu 20050706
505	! if((platform(7:11).eq.'METAR').or.(platform(7:11).eq.'SPECI').or.
506	! 1 (platform(7:10).eq.'SHIP').or.(platform(7:11).eq.'SYNOP').or.
507	! 1 (platform(1:4).eq.'SAMS'))
508	! 1 rko(n)=1.0
509	if(.NOT. is_sound) rko(n)=1.0
510	!yliu 20050706 end
511
512	! plfo is inFORMATion on what platform. May use this later in adjusting weights
513	plfo(n)=99.
514	if(platform(7:11).eq.'METAR')plfo(n)=1.
515	if(platform(7:11).eq.'SPECI')plfo(n)=2.
516	if(platform(7:10).eq.'SHIP')plfo(n)=3.
517	if(platform(7:11).eq.'SYNOP')plfo(n)=4.
518	if(platform(7:10).eq.'TEMP')plfo(n)=5.
519	if(platform(7:11).eq.'PILOT')plfo(n)=6.
520	if(platform(1:7).eq.'SATWNDS')plfo(n)=7.
521	if(platform(1:4).eq.'SAMS')plfo(n)=8.
522	if(platform(7:14).eq.'PROFILER')plfo(n)=9.
523	! yliu: ACARS->SATWINDS
524	if(platform(7:11).eq.'ACARS')plfo(n)=7.
525	! yliu: end
526	if(plfo(n).eq.99.) then
527	IF (iprt) print *,'n=',n,' unknown ob of type',platform
528	endif
529
530	!======================================================================
531	!======================================================================
532	! THIS PART READS SOUNDING INFO
533	IF(is_sound)THEN
534	nlevs_ob(n)=real(meas_count)
535	lev_in_ob(n)=1.
536	do imc=1,meas_count
537	! write(6,*) '0 inest = ',inest,' n = ',n
538	! the sounding has one header, many levels. This part puts it into
539	! "individual" observations. There's no other way for nudob to deal
540	! with it.
541	if(imc.gt.1)then ! sub-loop over N
542	n=n+1
543	if(n.gt.niobf)goto 120
544	nlevs_ob(n)=real(meas_count)
545	lev_in_ob(n)=real(imc)
546	timeob(n)=rtimob
547	rio(n)=ri
548	rjo(n)=rj
549	rko(n)=-99.
550	plfo(n)=plfo(n-imc+1)
551	elevob(n)=elevation
552	endif
553
554	read(nvol,104)pressure_data,pressure_qc, &
555	height_data,height_qc, &
556	temperature_data,temperature_qc, &
557	u_met_data,u_met_qc, &
558	v_met_data,v_met_qc, &
559	rh_data,rh_qc
560	104 FORMAT( 1x,6(f11.3,1x,f11.3,1x))
561
562	! yliu: Ensemble - add disturbance to upr obs
563	! if(plfo(n).eq.5.or.plfo(n).eq.6.or.plfo(n).eq.9) then FORE07E08
564	! if(imc.eq.1) then FORE07E08
565	! call srand(n)
566	! t_rand =- (rand(2)-0.5)*6
567	! call srand(n+100000)
568	! u_rand =- (rand(2)-0.5)*6
569	! call srand(n+200000)
570	! v_rand =- (rand(2)-0.5)*6
571	! endif FORE07E08
572	! if(temperature_qc.ge.0..and.temperature_qc.lt.30000..and.
573	! & temperature_data .gt. -88880.0 )
574	! & temperature_data = temperature_data + t_rand
575	! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
576	! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
577	! make sure at least 1 of the components is .ne.0
578	! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
579	! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
580	! u_met_data = u_met_data + u_rand
581	! v_met_data = v_met_data + v_rand
582	! endif
583	! endif FORE07E08
584	! yliu: Ens test - end
585
586
587	! jc
588	! hardwire to switch -777777. qc to 0. here temporarily
589	! -777777. is a sounding level that no qc was done on.
590
591	if(temperature_qc.eq.-777777.)temperature_qc=0.
592	if(pressure_qc.eq.-777777.)pressure_qc=0.
593	if(height_qc.eq.-777777.)height_qc=0.
594	if(u_met_qc.eq.-777777.)u_met_qc=0.
595	if(v_met_qc.eq.-777777.)v_met_qc=0.
596	if(rh_qc.eq.-777777.)rh_qc=0.
597	if(temperature_data.eq.-888888.)temperature_qc=-888888.
598	if(pressure_data.eq.-888888.)pressure_qc=-888888.
599	if(height_data.eq.-888888.)height_qc=-888888.
600	if(u_met_data.eq.-888888.)u_met_qc=-888888.
601	if(v_met_data.eq.-888888.)v_met_qc=-888888.
602	if(rh_data.eq.-888888.)rh_qc=-888888.
603
604	! jc
605	! Hardwire so that only use winds in pilot obs (no winds from temp) and
606	! only use temperatures and rh in temp obs (no temps from pilot obs)
607	! Do this because temp and pilot are treated as 2 platforms, but pilot
608	! has most of the winds, and temp has most of the temps. If use both,
609	! the second will smooth the effect of the first. Usually temps come in after
610	! pilots. pilots usually don't have any temps, but temp obs do have some
611	! winds usually.
612	! plfo=5 is TEMP ob, range sounding is an exception
613	!yliu start -- comment out to test with merged PILOT and TEMP and
614	! do not use obs interpolated by little_r
615	! if(plfo(n).eq.5. .and. namef(1:8).ne.'CLASS DA')then
616	! u_met_data=-888888.
617	! v_met_data=-888888.
618	! u_met_qc=-888888.
619	! v_met_qc=-888888.
620	! endif
621	if(plfo(n).eq.5..and.(u_met_qc.eq.256..or.v_met_qc.eq.256.))then
622	u_met_data=-888888.
623	v_met_data=-888888.
624	u_met_qc=-888888.
625	v_met_qc=-888888.
626	endif
627	!yliu end
628	! plfo=6 is PILOT ob
629	if(plfo(n).eq.6.)then
630	temperature_data=-888888.
631	rh_data=-888888.
632	temperature_qc=-888888.
633	rh_qc=-888888.
634	endif
635
636	!ajb Store potential temperature for WRF
637	if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
638
639	if(pressure_qc.ge.0..and.pressure_qc.lt.30000.)then
640
641	varobs(3,n) = &
642	temperature_data(100000./pressure_data)*RCP - t0
643
644	! write(6,*) 'reading data for N = ',n,' RCP = ',rcp
645	! write(6,*) 'temperature_data = ',temperature_data
646	! write(6,*) 'pressure_data = ',pressure_data
647	! write(6,*) 'varobs(3,n) = ',varobs(3,n)
648
649	else
650	varobs(3,n)=-888888.
651	endif
652
653	else
654	varobs(3,n)=-888888.
655	endif
656
657	if(pressure_qc.ge.0..and.pressure_qc.lt.30000.)then
658	! if(pressure_qc.ge.0.)then
659	varobs(5,n)=pressure_data
660	else
661	varobs(5,n)=-888888.
662	IF (iprt) THEN
663	print ,'******* PROBLEM ***********'
664	print *,'sounding, p undefined',latitude,longitude
665	ENDIF
666	endif
667	if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
668	! don't use data above 80 mb
669	if((varobs(5,n).gt.0.).and.(varobs(5,n).le.8.))then
670	u_met_data=-888888.
671	v_met_data=-888888.
672	u_met_qc=-888888.
673	v_met_qc=-888888.
674	temperature_data=-888888.
675	temperature_qc=-888888.
676	rh_data=-888888.
677	rh_qc=-888888.
678	endif
679
680	! yliu: add special processing of NPN and Range profiler
681	! only little_r interpolated and QC-ed data is used
682	if(namef(2:9).eq."PROFILER") then
683	if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
684	(v_met_qc.ge.0..and.v_met_qc.lt.30000.))then
685	!!yliu little_r already rotated the winds
686	! call vect(longitude,u_met_data,v_met_data,xlonc,xlatc,xn)
687	varobs(1,n)=u_met_data
688	varobs(2,n)=v_met_data
689	else
690	varobs(1,n)=-888888.
691	varobs(2,n)=-888888.
692	endif
693	else
694	if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
695	(v_met_qc.ge.0..and.v_met_qc.lt.30000.))then
696	!!yliu little_r already rotated the winds
697	! call vect(longitude,u_met_data,v_met_data,xlonc,xlatc,xn)
698	varobs(1,n)=u_met_data
699	varobs(2,n)=v_met_data
700	else
701	varobs(1,n)=-888888.
702	varobs(2,n)=-888888.
703	endif
704	endif
705	r_data=-888888.
706
707	if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
708	if((pressure_qc.ge.0.).and.(temperature_qc.ge.0.).and. &
709	(pressure_qc.lt.30000.).and.(temperature_qc.lt.30000.))then
710	call rh2r(rh_data,temperature_data,pressure_data*.01, &
711	r_data,0) ! yliu, change last arg from 1 to 0
712	else
713	! print *,'rh, but no t or p to convert',temperature_qc, &
714	! pressure_qc,n
715	r_data=-888888.
716	endif
717	endif
718	varobs(4,n)=r_data
719	enddo ! end do imc=1,meas_count
720	! print *,'--- sdng n=',n,nlevs_ob(n),lev_in_ob(n),timeob(n)
721	! read in non-sounding obs
722
723	ELSEIF(.NOT.is_sound)THEN
724	nlevs_ob(n)=1.
725	lev_in_ob(n)=1.
726	read(nvol,105)slp_data,slp_qc, &
727	ref_pres_data,ref_pres_qc, &
728	height_data,height_qc, &
729	temperature_data,temperature_qc, &
730	u_met_data,u_met_qc, &
731	v_met_data,v_met_qc, &
732	rh_data,rh_qc, &
733	psfc_data,psfc_qc, &
734	precip_data,precip_qc
735	105 FORMAT( 1x,9(f11.3,1x,f11.3,1x))
736
737	! Ensemble: add disturbance to sfc obs
738	! call srand(n)
739	! t_rand =+ (rand(2)-0.5)*5
740	! call srand(n+100000)
741	! u_rand =+ (rand(2)-0.5)*5
742	! call srand(n+200000)
743	! v_rand =+ (rand(2)-0.5)*5
744	! if(temperature_qc.ge.0..and.temperature_qc.lt.30000. .and.
745	! & temperature_data .gt. -88880.0 )
746	! & temperature_data = temperature_data + t_rand
747	! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
748	! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
749	! make sure at least 1 of the components is .ne.0
750	! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
751	! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
752	! u_met_data = u_met_data + u_rand
753	! v_met_data = v_met_data + v_rand
754	! endif
755	! yliu: Ens test - end
756
757	!Lilis
758
759	! calculate psfc if slp is there
760	if((psfc_qc.lt.0.).and.(slp_qc.ge.0..and.slp_qc.lt.30000.).and. &
761	(temperature_qc.ge.0..and.temperature_qc.lt.30000.).and. &
762	(slp_data.gt.90000.))then
763	tbar=temperature_data+0.5elevation.0065
764	psfc_data=slp_dataexp(-elevation/(rovgtbar))
765	varobs(5,n)=psfc_data*1.e-3
766	psfc_qc=0.
767	endif
768
769	!c No Very rough estimate of psfc from sfc elevation if UUtah ob and elev>1000m
770	! estimate psfc from temp and elevation
771	! Do not know sfc pressure in model at this point.
772	! if((psfc_qc.lt.0.).and.(elevation.gt.1000.).and.
773	! 1 (temperature_qc.ge.0..and.temperature_qc.lt.30000.)
774	! 1 .and.(platform(7:16).eq.'SYNOP PRET'))then
775	if((psfc_qc.lt.0.).and. &
776	(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
777	tbar=temperature_data+0.5elevation.0065
778	psfc_data=100000.exp(-elevation/(rovgtbar))
779	varobs(5,n)=psfc_data*1.e-3
780	psfc_qc=0.
781	endif
782
783	if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and.(psfc_data.gt.70000. &
784	.and.psfc_data.lt.105000.))then
785	varobs(5,n)=psfc_data
786	else
787	varobs(5,n)=-888888.
788	endif
789	if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
790
791	!Lilie
792
793
794	!ajb Store potential temperature for WRF
795	if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
796
797	if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and. &
798	(psfc_data.gt.70000. .and.psfc_data.lt.105000.))then
799
800	varobs(3,n) = &
801	temperature_data(100000./psfc_data)*RCP - t0
802	else
803	varobs(3,n)=-888888.
804	endif
805	else
806	varobs(3,n)=-888888.
807	endif
808
809	! if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and.(psfc_data.gt.70000. &
810	! .and.psfc_data.lt.105000.))then
811	! varobs(5,n)=psfc_data
812	! else
813	! varobs(5,n)=-888888.
814	! endif
815	! if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
816
817	if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
818	(v_met_qc.ge.0..and.v_met_qc.lt.30000.).and. &
819	! make sure at least 1 of the components is .ne.0
820	(u_met_data.ne.0..or.v_met_data.ne.0.))then
821	!!yliu little_r already rotated the winds
822	!!yliu call vect(longitude,u_met_data,v_met_data,xlonc,xlatc,xn)
823	varobs(1,n)=u_met_data
824	varobs(2,n)=v_met_data
825	else
826	varobs(1,n)=-888888.
827	varobs(2,n)=-888888.
828	endif
829	!! calculate psfc if slp is there
830	! if((psfc_qc.lt.0.).and.(slp_qc.ge.0..and.slp_qc.lt.30000.).and. &
831	! (temperature_qc.ge.0..and.temperature_qc.lt.30000.).and. &
832	! (slp_data.gt.90000.))then
833	! tbar=temperature_data+0.5elevation.0065
834	! psfc_data=slp_dataexp(-elevation/(rovgtbar))
835	! varobs(5,n)=psfc_data*1.e-3
836	! psfc_qc=0.
837	! endif
838
839	!!c No Very rough estimate of psfc from sfc elevation if UUtah ob and elev>1000m
840	!! estimate psfc from temp and elevation
841	!! Do not know sfc pressure in model at this point.
842	!! if((psfc_qc.lt.0.).and.(elevation.gt.1000.).and.
843	!! 1 (temperature_qc.ge.0..and.temperature_qc.lt.30000.)
844	!! 1 .and.(platform(7:16).eq.'SYNOP PRET'))then
845	! if((psfc_qc.lt.0.).and. &
846	! (temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
847	! tbar=temperature_data+0.5elevation.0065
848	! psfc_data=100000.exp(-elevation/(rovgtbar))
849	! varobs(5,n)=psfc_data*1.e-3
850	! psfc_qc=0.
851	! endif
852
853	! jc
854	! if a ship ob has rh<70%, then throw out
855
856	if(plfo(n).eq.3..and.rh_qc.ge.0..and.rh_data.lt.70.)then
857	rh_qc=-888888.
858	rh_data=-888888.
859	endif
860	!
861	r_data=-888888.
862	if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
863	if((psfc_qc.ge.0..and.psfc_qc.lt.30000.) &
864	.and.(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
865	! rh_data=amin1(rh_data,96.) ! yliu: do not allow surface to be saturated
866	call rh2r(rh_data,temperature_data,psfc_data*.01, &
867	r_data,0) ! yliu, change last arg from 1 to 0
868	else
869	! print *,'rh, but no t or p',temperature_data,
870	! 1 psfc_data,n
871	r_data=-888888.
872	endif
873	endif
874	varobs(4,n)=r_data
875	ELSE
876	IF (iprt) THEN
877	print *,' ====== '
878	print *,' NO Data Found '
879	ENDIF
880	ENDIF !end if(is_sound)
881	! END OF SFC OBS INPUT SECTION
882	!======================================================================
883	!======================================================================
884	! check if ob time is too early (only applies to beginning)
885	IF(RTIMOB.LT.TBACK-fdob%WINDOW)then
886	IF (iprt) print *,'ob too early'
887	n=n-1
888	GOTO 110
889	ENDIF
890
891	! check if this ob is a duplicate
892	! this check has to be before other checks
893	njend=n-1
894	if(is_sound)njend=n-meas_count
895	do njc=1,njend
896	! Check that time, lat, lon, and platform all match exactly.
897	! Platforms 1-4 (surface obs) can match with each other. Otherwise,
898	! platforms have to match exactly.
899	if( (timeob(n).eq.timeob(njc)) .and. &
900	(rio(n).eq.rio(njc)) .and. &
901	(rjo(n).eq.rjo(njc)) .and. &
902	(plfo(njc).ne.99.) ) then
903	!yliu: if two sfc obs are departed less than 1km, consider they are redundant
904	! (abs(rio(n)-rio(njc))*dscg.gt.1000.) &
905	! .or. (abs(rjo(n)-rjo(njc))*dscg.gt.1000.) &
906	! .or. (plfo(njc).eq.99.) )goto 801
907	!yliu end
908	! If platforms different, and either > 4, jump out
909	if( ( (plfo(n).le.4.).and.(plfo(njc).le.4.) ) .or. &
910	(plfo(n).eq.plfo(njc)) ) then
911
912	! if not a sounding, and levels are the same then replace first occurrence
913	if((.not.is_sound).and.(rko(njc).eq.rko(n))) then
914	! print *,'dup single ob-replace ',n,inest,
915	! plfo(n),plfo(njc)
916	! this is the sfc ob replacement part
917	VAROBS(1,njc)=VAROBS(1,n)
918	VAROBS(2,njc)=VAROBS(2,n)
919	VAROBS(3,njc)=VAROBS(3,n)
920	VAROBS(4,njc)=VAROBS(4,n)
921	VAROBS(5,njc)=VAROBS(5,n)
922	! don't need to switch these because they're the same
923	! RIO(njc)=RIO(n)
924	! RJO(njc)=RJO(n)
925	! RKO(njc)=RKO(n)
926	! TIMEOB(njc)=TIMEOB(n)
927	! nlevs_ob(njc)=nlevs_ob(n)
928	! lev_in_ob(njc)=lev_in_ob(n)
929	! plfo(njc)=plfo(n)
930	! end sfc ob replacement part
931
932	n=n-1
933	goto 100
934	! It's harder to fix the soundings, since the number of levels may be different
935	! The easiest thing to do is to just set the first occurrence to all missing, and
936	! keep the second occurrence, or vice versa.
937	! For temp or profiler keep the second, for pilot keep the one with more levs
938	! This is for a temp or prof sounding, equal to same
939	! also if a pilot, but second one has more obs
940	elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
941	( (plfo(njc).eq.5.).or.(plfo(njc).eq.9.).or. &
942	( (plfo(njc).eq.6.).and. &
943	(nlevs_ob(n).ge.nlevs_ob(njc)) ) ) )then
944	IF (iprt) THEN
945	print *,'duplicate sounding - eliminate first occurrence', &
946	n,inest,meas_count,nlevs_ob(njc), &
947	latitude,longitude,plfo(njc)
948	ENDIF
949	if(lev_in_ob(njc).ne.1.) then
950	IF (iprt) THEN
951	print , 'problem ****** - dup sndg ', &
952	lev_in_ob(njc),nlevs_ob(njc)
953	ENDIF
954	endif
955	! n=n-meas_count
956	! set the first sounding ob to missing
957	do njcc=njc,njc+nint(nlevs_ob(njc))-1
958	VAROBS(1,njcc)=-888888.
959	VAROBS(2,njcc)=-888888.
960	VAROBS(3,njcc)=-888888.
961	VAROBS(4,njcc)=-888888.
962	VAROBS(5,njcc)=-888888.
963	plfo(njcc)=99.
964	enddo
965	goto 100
966	! if a pilot, but first one has more obs
967	elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
968	(plfo(njc).eq.6.).and. &
969	(nlevs_ob(n).lt.nlevs_ob(njc)) )then
970	IF (iprt) THEN
971	print *, &
972	'duplicate pilot sounding - eliminate second occurrence', &
973	n,inest,meas_count,nlevs_ob(njc), &
974	latitude,longitude,plfo(njc)
975	ENDIF
976	if(lev_in_ob(njc).ne.1.) then
977	IF (iprt) THEN
978	print , 'problem ****** - dup sndg ', &
979	lev_in_ob(njc),nlevs_ob(njc)
980	ENDIF
981	endif
982	n=n-meas_count
983
984	!ajb Reset timeob for discarded indices.
985	do imc = n+1, n+meas_count
986	timeob(imc) = 99999.
987	enddo
988	goto 100
989	! This is for a single-level satellite upper air ob - replace first
990	elseif( (is_sound).and. &
991	(nlevs_ob(njc).eq.1.).and. &
992	(nlevs_ob(n).eq.1.).and. &
993	(varobs(5,njc).eq.varobs(5,n)).and. &
994	(plfo(njc).eq.7.).and.(plfo(n).eq.7.) ) then
995	IF (iprt) print *, &
996	'duplicate single lev sat-wind ob - replace first',n, &
997	inest,meas_count,varobs(5,n)
998	! this is the single ua ob replacement part
999	VAROBS(1,njc)=VAROBS(1,n)
1000	VAROBS(2,njc)=VAROBS(2,n)
1001	VAROBS(3,njc)=VAROBS(3,n)
1002	VAROBS(4,njc)=VAROBS(4,n)
1003	VAROBS(5,njc)=VAROBS(5,n)
1004	! don't need to switch these because they're the same
1005	! RIO(njc)=RIO(n)
1006	! RJO(njc)=RJO(n)
1007	! RKO(njc)=RKO(n)
1008	! TIMEOB(njc)=TIMEOB(n)
1009	! nlevs_ob(njc)=nlevs_ob(n)
1010	! lev_in_ob(njc)=lev_in_ob(n)
1011	! plfo(njc)=plfo(n)
1012	! end single ua ob replacement part
1013	n=n-1
1014	goto 100
1015	else
1016	IF (iprt) THEN
1017	print *,'duplicate location, but no match otherwise',n,njc, &
1018	plfo(n),varobs(5,n),nlevs_ob(n),lev_in_ob(n), &
1019	plfo(njc),varobs(5,njc),nlevs_ob(njc),lev_in_ob(njc)
1020	ENDIF
1021	endif
1022	endif
1023	endif
1024	! end of njc do loop
1025	enddo
1026
1027	! check if ob is a sams ob that came in via UUtah - discard
1028	if( plfo(n).eq.4..and.(platform(7:16).eq.'SYNOP PRET').and. &
1029	(id(7:15).eq.'METNET= 3') )then
1030	! print *,'elim metnet=3',latitude,longitude,rtimob
1031	n=n-1
1032	goto 100
1033	endif
1034
1035	! check if ob is in the domain
1036	if( (ri.lt.2.).or.(ri.gt.real(e_we-1)).or.(rj.lt.2.).or. &
1037	(rj.gt.real(e_sn-1)) ) then
1038	! if (iprt) write(6,*) 'Obs out of coarse mesh domain'
1039	! write(6,*) 'we_maxcg-1 = ',real(we_maxcg-1)
1040	! write(6,*) 'sn_maxcg-1 = ',real(sn_maxcg-1)
1041
1042	! n=n-1
1043	! if(is_sound)n=n-meas_count+1
1044
1045	n=n-meas_count
1046	!ajb Reset timeob for discarded indices.
1047	do imc = n+1, n+meas_count
1048	timeob(imc) = 99999.
1049	enddo
1050	goto 100
1051	endif
1052
1053	! check if an upper air ob is too high
1054	! the ptop here is hardwired
1055	! this check has to come after other checks - usually the last few
1056	! upper air obs are too high
1057	! if(is_sound)then
1058	! njc=meas_count
1059	! do jcj=meas_count,1,-1
1060	! 6. is 60 mb - hardwired
1061	! if((varobs(5,n).lt.6.).and.(varobs(5,n).gt.0.))then
1062	! print *,'obs too high - eliminate,n,p=',n,varobs(5,n)
1063	! n=n-1
1064	! else
1065	! if(varobs(5,n).gt.0.)goto 100
1066	! endif
1067	! enddo
1068	! endif
1069	!
1070	IF(TIMEOB(N).LT.fdob%RTLAST) THEN
1071	IF (iprt) THEN
1072	PRINT *,'2 OBS ARE NOT IN CHRONOLOGICAL ORDER'
1073	PRINT *,'NEW YEAR?'
1074	print *,'timeob,rtlast,n=',timeob(n),fdob%rtlast,n
1075	ENDIF
1076	STOP 111
1077	ELSE
1078	fdob%RTLAST=TIMEOB(N)
1079	ENDIF
1080	GOTO 100
1081	111 CONTINUE
1082	!**********************************************************************
1083	! -------------- END BIG 100 LOOP OVER N --------------
1084	!**********************************************************************
1085	IF (iprt) write(6,5403) NVOL,XTIME
1086	IEOF(inest)=1
1087
1088	close(NVOLA+INEST-1)
1089	IF (iprt) print *,'closed fdda file for inest=',inest,nsta
1090
1091	! if(nsta.eq.1.and.timeob(1).gt.9.e4)nsta=0
1092	goto 1001
1093
1094	! THE OBSERVATION ARRAYS ARE FULL AND THE MOST RECENTLY
1095	! ACQUIRED OBS STILL HAS TIMEOB .LE. TFORWD. SO START
1096	! DECREASING THE SIZE OF THE WINDOW
1097	! get here if too many obs
1098	120 CONTINUE
1099	IF (iprt) THEN
1100	write(6,121) N,NIOBF
1101	write(6,122)
1102	ENDIF
1103	STOP 122
1104	fdob%WINDOW=fdob%WINDOW-0.1*TWINDO
1105	IF(TWINDO.LT.0)STOP 120
1106	! IF THE WINDOW BECOMES NEGATIVE, THE INCOMING DATA IS
1107	! PROBABLY GARBLED. STOP.
1108	GOTO 10
1109	!
1110	! READ CYCLE IS COMPLETED. DETERMINE THE NUMBER OF OBS IN
1111	! THE CURRENT WINDOW
1112	!
1113	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1114	! BUT FIRST, WHEN KTAU.EQ.0 (OR IN GENERAL, KTAUR), DISCARD THE
1115	! "OLD" OBS FIRST...
1116	130 CONTINUE
1117
1118	! get here if at end of file, or if obs time is beyond what we
1119	! need right now
1120	IF(KTAU.EQ.KTAUR)THEN
1121	NSTA=0
1122	keep_obs : DO N=1,NIOBF
1123
1124	! try to keep all obs, but just don't use yet
1125	! (don't want to throw away last obs read in - especially if
1126	! its a sounding, in which case it looks like many obs)
1127	IF(TIMEOB(N).GT.9.e4) EXIT keep_obs
1128	if(timeob(n).gt.tforwd) then
1129	if(iprt) write(6,951)inest,n,timeob(n),tforwd
1130	951 FORMAT('saving ob beyond window,inest,n,timeob,tforwd=', &
1131	2i5,2f13.4)
1132	endif
1133	NSTA=N
1134	ENDDO keep_obs
1135
1136	NDUM=0
1137	! make time=99999. if ob is too old
1138	! print *,'tback,nsta=',tback,nsta
1139	old_obs : DO N=1,NSTA+1
1140	IF((TIMEOB(N)-TBACK).LT.0)THEN
1141	TIMEOB(N)=99999.
1142	ENDIF
1143	! print *,'n,ndum,timeob=',n,ndum,timeob(n)
1144	IF(TIMEOB(N).LT.9.E4) EXIT old_obs
1145	NDUM=N
1146	ENDDO old_obs
1147
1148	! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
1149	IF (iprt) THEN
1150	print *,'after 190 ndum=',ndum,nsta
1151	ENDIF
1152	NDUM=ABS(NDUM)
1153	NMOVE=NIOBF-NDUM
1154	IF( NMOVE.GT.0 .AND. NDUM.NE.0) THEN
1155	DO N=1,NMOVE
1156	VAROBS(1,N)=VAROBS(1,N+NDUM)
1157	VAROBS(2,N)=VAROBS(2,N+NDUM)
1158	VAROBS(3,N)=VAROBS(3,N+NDUM)
1159	VAROBS(4,N)=VAROBS(4,N+NDUM)
1160	VAROBS(5,N)=VAROBS(5,N+NDUM)
1161	RJO(N)=RJO(N+NDUM)
1162	RIO(N)=RIO(N+NDUM)
1163	RKO(N)=RKO(N+NDUM)
1164	TIMEOB(N)=TIMEOB(N+NDUM)
1165	nlevs_ob(n)=nlevs_ob(n+ndum)
1166	lev_in_ob(n)=lev_in_ob(n+ndum)
1167	plfo(n)=plfo(n+ndum)
1168	ENDDO
1169	ENDIF
1170	! moved obs up. now fill remaining space with 99999.
1171	NOPEN=NMOVE+1
1172	IF(NOPEN.LE.NIOBF) THEN
1173	DO N=NOPEN,NIOBF
1174	VAROBS(1,N)=99999.
1175	VAROBS(2,N)=99999.
1176	VAROBS(3,N)=99999.
1177	VAROBS(4,N)=99999.
1178	VAROBS(5,N)=99999.
1179	RIO(N)=99999.
1180	RJO(N)=99999.
1181	RKO(N)=99999.
1182	TIMEOB(N)=99999.
1183	ENDDO
1184	ENDIF
1185	ENDIF
1186	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1187	NSTA=0
1188	! print *,'nsta at restart setting is ',nsta
1189	! recalculate nsta after moving things around
1190	recalc : DO N=1,NIOBF
1191	! try to save all obs - don't throw away latest read in
1192	IF(TIMEOB(N).GT.9.e4) EXIT recalc
1193	NSTA=N
1194	! nsta=n-1 ! yliu test
1195	ENDDO recalc
1196
1197	! Find the number of stations that are actually within the time window.
1198	nstaw = nvals_le_limit(nsta, timeob, tforwd)
1199
1200	! Print obs information, according to nobs_prt, but limit to 100 max.
1201	if (iprt) then
1202	ips = min(nstaw,nobs_prt,100)
1203	if(ips.gt.0) then
1204	write(6,'(/a,i4,a,i2)') 'MASS-PT LOCATIONS FOR FIRST',ips, &
1205	' OBSERVATIONS FOR NEST ',inest
1206	write(6,*) ' OBS# I J LAT LON TIME(hrs)'
1207	endif
1208	! Note: rio and rjo are referenced to non-staggered grid (not mass-point!)
1209	! Hence subtract .5 from each to get mass-point coords.
1210	do n=1,ips
1211	write(6,'(2x,i4,2f8.3,2f8.2,3x,f8.5)') &
1212	n,rio(n)-.5,rjo(n)-.5,lat_prt(n),lon_prt(n),timeob(n)
1213	enddo
1214	endif
1215
1216	IF (iprt) write(6,160) KTAU,XTIME,NSTAW
1217	IF(KTAU.EQ.KTAUR)THEN
1218	IF(nudge_opt.EQ.1)THEN
1219	TWDOP=TWINDO*60.
1220	IF (iprt) THEN
1221	write(6,1449) INEST,RINXY,RINSIG,TWDOP
1222	IF(ISWIND.EQ.1) write(6,1450) GIV
1223	IF(ISTEMP.EQ.1) write(6,1451) GIT
1224	IF(ISMOIS.EQ.1) write(6,1452) GIQ
1225	IF(ISPSTR.EQ.1) write(6,1453) GIP
1226	ENDIF
1227	ENDIF
1228	ENDIF
1229	IF(KTAU.EQ.KTAUR)THEN
1230	IF (iprt) THEN
1231	write(6,553)
1232	write(6,554)
1233	ENDIF
1234	IF(fdob%IWTSIG.NE.1)THEN
1235	IF (iprt) THEN
1236	write(6,555)
1237	write(6,556) fdob%RINFMNRINXY,fdob%RINFMXRINXY,fdob%PFREE*10.
1238	ENDIF
1239	IF(fdob%RINFMN.GT.fdob%RINFMX)STOP 556
1240	! IS MINIMUM GREATER THAN MAXIMUM?
1241	IF (iprt) write(6,557) fdob%DPSMX*10.,fdob%DCON
1242	IF(fdob%DPSMX.GT.10.)STOP 557
1243	ENDIF
1244	ENDIF
1245	! IS DPSMX IN CB?
1246
1247	IF(KTAU.EQ.KTAUR)THEN
1248	IF (iprt) write(6,601) INEST,IONF
1249	ENDIF
1250	fdob%NSTAT=NSTA
1251	fdob%NSTAW=NSTAW
1252
1253	555 FORMAT(1X,' ABOVE THE SURFACE LAYER, OBS NUDGING IS PERFORMED', &
1254	' ON PRESSURE LEVELS,')
1255	556 FORMAT(1X,' WHERE RINXY VARIES LINEARLY FROM ',E11.3,' KM AT', &
1256	' THE SURFACE TO ',E11.3,' KM AT ',F7.2,' MB AND ABOVE')
1257	557 FORMAT(1X,' IN THE SURFACE LAYER, WXY IS A FUNCTION OF ', &
1258	'DPSMX = ',F7.2,' MB WITH DCON = ',E11.3, &
1259	' - SEE SUBROUTINE NUDOB')
1260	601 FORMAT('0','FOR EFFICIENCY, THE OBS NUDGING FREQUENCY ', &
1261	'FOR MESH #',I2,' IS ',1I2,' CGM TIMESTEPS ',/)
1262	121 FORMAT('0',' WARNING: NOBS = ',I4,' IS GREATER THAN NIOBF = ', &
1263	I4,': INCREASE PARAMETER NIOBF')
1264	5403 FORMAT(1H0,'-------------EOF REACHED FOR NVOL = ',I3, &
1265	' AND XTIME = ',F10.2,'-------------------')
1266	122 FORMAT(1X,' ...OR THE CODE WILL REDUCE THE TIME WINDOW')
1267	160 FORMAT('0','****** CALL IN4DOB AT KTAU = ',I5,' AND XTIME = ', &
1268	F10.2,': NSTA = ',I7,' ******')
1269	1449 FORMAT(1H0,'*****NUDGING INDIVIDUAL OBS ON MESH #',I2, &
1270	' WITH RINXY = ', &
1271	E11.3,' KM, RINSIG = ',E11.3,' AND TWINDO (HALF-PERIOD) = ', &
1272	E11.3,' MIN')
1273	1450 FORMAT(1X,'NUDGING IND. OBS WINDS WITH GIV = ',E11.3)
1274	1451 FORMAT(1X,'NUDGING IND. OBS TEMPERATURE WITH GIT = ',E11.3)
1275	1452 FORMAT(1X,'NUDGING IND. OBS MOISTURE WITH GIQ = ',E11.3)
1276	1453 FORMAT(1X,'NUDGING IND. OBS SURFACE PRESSURE WITH GIP = ,'E11.3)
1277	553 FORMAT(1X,'BY DEFAULT: OBS NUDGING OF TEMPERATURE AND MOISTURE ', &
1278	'IS RESTRICTED TO ABOVE THE BOUNDARY LAYER')
1279	554 FORMAT(1X,'...WHILE OBS NUDGING OF WIND IS INDEPENDENT OF THE ', &
1280	'BOUNDARY LAYER')
1281
1282	RETURN
1283	END SUBROUTINE in4dob
1284
1285	SUBROUTINE julgmt(mdate,julgmtn,timanl,julday,gmt,ind)
1286	! CONVERT MDATE YYMMDDHH TO JULGMT (JULIAN DAY * 100. +GMT)
1287	! AND TO TIMANL (TIME IN MINUTES WITH RESPECT TO MODEL TIME)
1288	! IF IND=0 INPUT MDATE, OUTPUT JULGMTN AND TIMANL
1289	! IF IND=1 INPUT TIMANL, OUTPUT JULGMTN
1290	! IF IND=2 INPUT JULGMTN, OUTPUT TIMANL
1291	INTEGER, intent(in) :: MDATE
1292	REAL, intent(out) :: JULGMTN
1293	REAL, intent(out) :: TIMANL
1294	INTEGER, intent(in) :: JULDAY
1295	REAL, intent(in) :: GMT
1296	INTEGER, intent(in) :: IND
1297
1298	!*** DECLARATIONS FOR IMPLICIT NONE
1299	real :: MO(12), rjulanl, houranl, rhr
1300
1301	integer :: iyr, idate1, imo, idy, ihr, my1, my2, my3, ileap
1302	integer :: juldayn, juldanl, idymax, mm
1303
1304
1305	IF(IND.EQ.2)GOTO 150
1306	IYR=INT(MDATE/1000000.+0.001)
1307	IDATE1=MDATE-IYR*1000000
1308	IMO=INT(IDATE1/10000.+0.001)
1309	IDY=INT((IDATE1-IMO*10000.)/100.+0.001)
1310	IHR=IDATE1-IMO10000-IDY100
1311	MO(1)=31
1312	MO(2)=28
1313	! IS THE YEAR A LEAP YEAR? (IN THIS CENTURY)
1314	IYR=IYR+1900
1315	MY1=MOD(IYR,4)
1316	MY2=MOD(IYR,100)
1317	MY3=MOD(IYR,400)
1318	ILEAP=0
1319	! jc
1320	! IF(MY1.EQ.0.AND.MY2.NE.0.OR.MY3.EQ.0)THEN
1321	IF(MY1.EQ.0)THEN
1322	ILEAP=1
1323	MO(2)=29
1324	ENDIF
1325	IF(IND.EQ.1)GOTO 200
1326	MO(3)=31
1327	MO(4)=30
1328	MO(5)=31
1329	MO(6)=30
1330	MO(7)=31
1331	MO(8)=31
1332	MO(9)=30
1333	MO(10)=31
1334	MO(11)=30
1335	MO(12)=31
1336	JULDAYN=0
1337	DO 100 MM=1,IMO-1
1338	JULDAYN=JULDAYN+MO(MM)
1339	100 CONTINUE
1340
1341	IF(IHR.GE.24)THEN
1342	IDY=IDY+1
1343	IHR=IHR-24
1344	ENDIF
1345	JULGMTN=(JULDAYN+IDY)*100.+IHR
1346	! CONVERT JULGMT TO TIMANL WRT MODEL TIME IN MINUTES (XTIME)
1347	150 CONTINUE
1348	JULDANL=INT(JULGMTN/100.+0.000001)
1349	RJULANL=FLOAT(JULDANL)*100.
1350	HOURANL=JULGMTN-RJULANL
1351	TIMANL=(FLOAT(JULDANL-JULDAY)24.-GMT+HOURANL)60.
1352	RETURN
1353	200 CONTINUE
1354	RHR=GMT+TIMANL/60.+0.000001
1355	IDY=JULDAY
1356	IDYMAX=365+ILEAP
1357	300 IF(RHR.GE.24.0)THEN
1358	RHR=RHR-24.0
1359	IDY=IDY+1
1360	GOTO 300
1361	ENDIF
1362	IF(IDY.GT.IDYMAX)IDY=IDY-IDYMAX
1363	JULGMTN=FLOAT(IDY)*100.+RHR
1364	RETURN
1365	END SUBROUTINE julgmt
1366
1367	SUBROUTINE rh2r(rh,t,p,r,iice)
1368
1369	! convert rh to r
1370	! if iice=1, use saturation with respect to ice
1371	! rh is 0-100.
1372	! r is g/g
1373	! t is K
1374	! p is mb
1375	!
1376	REAL, intent(in) :: rh
1377	REAL, intent(in) :: t
1378	REAL, intent(in) :: p
1379	REAL, intent(out) :: r
1380	INTEGER, intent(in) :: iice
1381
1382	!*** DECLARATIONS FOR IMPLICIT NONE
1383	real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
1384	real esat, rsat
1385
1386	eps=0.62197
1387	e0=6.1078
1388	eslcon1=17.2693882
1389	eslcon2=35.86
1390	esicon1=21.8745584
1391	esicon2=7.66
1392	t0=260.
1393
1394	! print *,'rh2r input=',rh,t,p
1395	rh1=rh*.01
1396
1397	if(iice.eq.1.and.t.le.t0)then
1398	esat=e0exp(esicon1(t-273.16)/(t-esicon2))
1399	else
1400	esat=e0exp(eslcon1(t-273.16)/(t-eslcon2))
1401	endif
1402	rsat=eps*esat/(p-esat)
1403	! print *,'rsat,esat=',rsat,esat
1404	r=rh1*rsat
1405
1406	! print *,'rh2r rh,t,p,r=',rh1,t,p,r
1407
1408	return
1409	END SUBROUTINE rh2r
1410
1411	SUBROUTINE rh2rb(rh,t,p,r,iice)
1412
1413	! convert rh to r
1414	! if iice=1, use daturation with respect to ice
1415	! rh is 0-100.
1416	! r is g/g
1417	! t is K
1418	! p is mb
1419
1420	REAL, intent(in) :: rh
1421	REAL, intent(in) :: t
1422	REAL, intent(in) :: p
1423	REAL, intent(out) :: r
1424	INTEGER, intent(in) :: iice
1425
1426	!*** DECLARATIONS FOR IMPLICIT NONE
1427	real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
1428	real esat, rsat
1429
1430	eps=0.622
1431	e0=6.112
1432	eslcon1=17.67
1433	eslcon2=29.65
1434	esicon1=22.514
1435	esicon2=6.15e3
1436	t0=273.15
1437
1438	print *,'rh2r input=',rh,t,p
1439	rh1=rh*.01
1440
1441	if(iice.eq.1.and.t.le.t0)then
1442	esat=e0*exp(esicon1-esicon2/t)
1443	else
1444	esat=e0exp(eslcon1(t-t0)/(t-eslcon2))
1445	endif
1446	rsat=eps*esat/(p-esat)
1447	! print *,'rsat,esat=',rsat,esat
1448	r=rh1epsrsat/(eps+rsat*(1.-rh1))
1449
1450	print *,'rh2r rh,t,p,r=',rh1,t,p,r
1451
1452	return
1453	END SUBROUTINE rh2rb
1454
1455	SUBROUTINE set_projection (obs_proj, map_proj, cen_lat, cen_lon, &
1456	true_lat1, true_lat2, stand_lon, &
1457	known_lat, known_lon, &
1458	e_we, e_sn, dxm, dym )
1459
1460	USE module_llxy
1461
1462	!*************************************************************************
1463	! Purpose: Set map projection information which will be used to map the
1464	! observation (lat,lon) location to its corresponding (x,y)
1465	! location on the WRF (coarse) grid. using the selected map
1466	! projection (e.g., Lambert, Mercator, Polar Stereo, etc).
1467	!*************************************************************************
1468
1469	IMPLICIT NONE
1470
1471	TYPE(PROJ_INFO), intent(out) :: obs_proj ! structure for obs projection info.
1472	INTEGER, intent(in) :: map_proj ! map projection index
1473	REAL, intent(in) :: cen_lat ! center latitude for map projection
1474	REAL, intent(in) :: cen_lon ! center longiture for map projection
1475	REAL, intent(in) :: true_lat1 ! truelat1 for map projection
1476	REAL, intent(in) :: true_lat2 ! truelat2 for map projection
1477	REAL, intent(in) :: stand_lon ! standard longitude for map projection
1478	INTEGER, intent(in) :: e_we ! max grid index in south-north coordinate
1479	INTEGER, intent(in) :: e_sn ! max grid index in west-east coordinate
1480	REAL, intent(in) :: known_lat ! latitude of domain origin point (i,j)=(1,1)
1481	REAL, intent(in) :: known_lon ! longigude of domain origin point (i,j)=(1,1)
1482	REAL, intent(in) :: dxm ! grid size in x (meters)
1483	REAL, intent(in) :: dym ! grid size in y (meters)
1484
1485	! Set up map transformation structure
1486	CALL map_init(obs_proj)
1487
1488	! Mercator
1489	IF (map_proj == PROJ_MERC) THEN
1490	CALL map_set(PROJ_MERC, obs_proj, &
1491	truelat1 = true_lat1, &
1492	lat1 = known_lat, &
1493	lon1 = known_lon, &
1494	knowni = 1., &
1495	knownj = 1., &
1496	dx = dxm)
1497
1498	! Lambert conformal
1499	ELSE IF (map_proj == PROJ_LC) THEN
1500	CALL map_set(PROJ_LC, obs_proj, &
1501	truelat1 = true_lat1, &
1502	truelat2 = true_lat2, &
1503	stdlon = stand_lon, &
1504	lat1 = known_lat, &
1505	lon1 = known_lon, &
1506	knowni = 1., &
1507	knownj = 1., &
1508	dx = dxm)
1509
1510	! Polar stereographic
1511	ELSE IF (map_proj == PROJ_PS) THEN
1512	CALL map_set(PROJ_PS, obs_proj, &
1513	truelat1 = true_lat1, &
1514	stdlon = stand_lon, &
1515	lat1 = known_lat, &
1516	lon1 = known_lon, &
1517	knowni = 1., &
1518	knownj = 1., &
1519	dx = dxm)
1520	! Cassini (global ARW)
1521	ELSE IF (map_proj == PROJ_CASSINI) THEN
1522	CALL map_set(PROJ_CASSINI, obs_proj, &
1523	latinc = dym360.0/(2.0EARTH_RADIUS_M*PI), &
1524	loninc = dxm360.0/(2.0EARTH_RADIUS_M*PI), &
1525	lat1 = known_lat, &
1526	lon1 = known_lon, &
1527	! We still need to get POLE_LAT and POLE_LON metadata variables before
1528	! this will work for rotated poles.
1529	lat0 = 90.0, &
1530	lon0 = 0.0, &
1531	knowni = 1., &
1532	knownj = 1., &
1533	stdlon = stand_lon)
1534
1535	! Rotated latitude-longitude
1536	ELSE IF (map_proj == PROJ_ROTLL) THEN
1537	CALL map_set(PROJ_ROTLL, obs_proj, &
1538	! I have no idea how this should work for NMM nested domains
1539	ixdim = e_we-1, &
1540	jydim = e_sn-1, &
1541	phi = real(e_sn-2)*dym/2.0, &
1542	lambda = real(e_we-2)*dxm, &
1543	lat1 = cen_lat, &
1544	lon1 = cen_lon, &
1545	latinc = dym, &
1546	loninc = dxm, &
1547	stagger = HH)
1548
1549	END IF
1550
1551	! write(6,*) 'ajb init: map_proj = ',map_proj
1552	! write(6,*) 'ajb: after setting map:'
1553	! write(6,*) 'truelat1 = ',obs_proj%truelat1
1554	! write(6,*) 'truelat2 = ',obs_proj%truelat2
1555	! write(6,*) 'stdlon = ',obs_proj%stdlon
1556	! write(6,*) 'lat1 = ',obs_proj%lat1
1557	! write(6,*) 'lon1 = ',obs_proj%lon1
1558	! write(6,*) 'knowni = ',obs_proj%knowni
1559	! write(6,*) 'knownj = ',obs_proj%knownj
1560	! write(6,*) 'dx = ',obs_proj%dx
1561
1562	END SUBROUTINE set_projection
1563
1564	INTEGER FUNCTION nvals_le_limit(isize, values, limit)
1565	!------------------------------------------------------------------------------
1566	! PURPOSE: Return the number of values in a (real) non-decreasing array that
1567	! are less than or equal to the specified upper limit.
1568	! NOTE: It is important that the array is non-decreasing!
1569	!
1570	!------------------------------------------------------------------------------
1571	IMPLICIT NONE
1572
1573	INTEGER, INTENT(IN) :: isize ! Size of input array
1574	REAL, INTENT(IN) :: values(isize) ! Input array of reals
1575	REAL, INTENT(IN) :: limit ! Upper limit
1576
1577	! Local variables
1578	integer :: n
1579
1580	! Search the array from largest to smallest values.
1581	find_nvals: DO n = isize, 1, -1
1582	if(values(n).le.limit) EXIT find_nvals
1583	ENDDO find_nvals
1584	nvals_le_limit = n
1585
1586	RETURN
1587	END FUNCTION nvals_le_limit
1588
1589	#endif
1590	!-----------------------------------------------------------------------
1591	! End subroutines for in4dob
1592	!-----------------------------------------------------------------------

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

Download in other formats:

Original Format