Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

wrf_fddaobs_in.F @ 3576

Last change on this file since 3576 was 2759, checked in by aslmd, 3 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 61.5 KB

Rev	Line
[2759]	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