Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

wrf_fddaobs_in.F @ 2756

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

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

Download in other formats:

Original Format