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