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