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