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