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,rj,ri,xlatc,xlonc,map_proj,
435 call llxy (latitude,longitude,ri,rj,xlatc_cg,xlonc_cg,map_proj, &
436 true_lat1,true_lat2,dscg,xn,sn_maxcg,we_maxcg, &
437 1,1,1)
438
439 !ajb ri and rj are referenced to the non-staggered grid. (For MM5, they
440 ! were referenced to the dot grid.)
441
442 rio(n)=ri
443 rjo(n)=rj
444
445 if (iprt) THEN
446 if(n.le.10) then
447 write(6,'(/,a,i5,a,f8.2,a,f8.2,a,f8.2,a,f8.2/)') &
448 ' OBS N = ',n, &
449 ' RIO = ',rio(n), &
450 ' RJO = ',rjo(n), &
451 ' LAT = ',latitude, &
452 ' LON = ',longitude
453 endif
454 endif
455
456 read(nvol,1021)id,namef
457 1021 FORMAT(2x,2(a40,3x))
458 read(nvol,103)platform,source,elevation,is_sound,bogus,meas_count
459 103 FORMAT( 2x,2(a16,2x),f8.0,2x,2(l4,2x),i5)
460
461 ! write(6,*) '----- OBS description ----- '
462 ! write(6,*) 'platform,source,elevation,is_sound,bogus,meas_count:'
463 ! write(6,*) platform,source,elevation,is_sound,bogus,meas_count
464
465 ! yliu
466 elevob(n)=elevation
467 ! jc
468 ! change platform from synop to profiler when needed
469 if(namef(2:9).eq.'PROFILER')platform(7:14)='PROFILER'
470 ! yliu
471 if(namef(2:6).eq.'ACARS')platform(7:11)='ACARS'
472 if(namef(1:7).eq.'SATWNDS') platform(1:11)='SATWNDS '
473 if(namef(1:8).eq.'CLASS DA')platform(7:10)='TEMP'
474 ! yliu end
475
476 rko(n)=-99.
477 !yliu 20050706
478 ! if((platform(7:11).eq.'METAR').or.(platform(7:11).eq.'SPECI').or.
479 ! 1 (platform(7:10).eq.'SHIP').or.(platform(7:11).eq.'SYNOP').or.
480 ! 1 (platform(1:4).eq.'SAMS'))
481 ! 1 rko(n)=1.0
482 if(.NOT. is_sound) rko(n)=1.0
483 !yliu 20050706 end
484
485 ! plfo is inFORMATion on what platform. May use this later in adjusting weights
486 plfo(n)=99.
487 if(platform(7:11).eq.'METAR')plfo(n)=1.
488 if(platform(7:11).eq.'SPECI')plfo(n)=2.
489 if(platform(7:10).eq.'SHIP')plfo(n)=3.
490 if(platform(7:11).eq.'SYNOP')plfo(n)=4.
491 if(platform(7:10).eq.'TEMP')plfo(n)=5.
492 if(platform(7:11).eq.'PILOT')plfo(n)=6.
493 if(platform(1:7).eq.'SATWNDS')plfo(n)=7.
494 if(platform(1:4).eq.'SAMS')plfo(n)=8.
495 if(platform(7:14).eq.'PROFILER')plfo(n)=9.
496 ! yliu: ACARS->SATWINDS
497 if(platform(7:11).eq.'ACARS')plfo(n)=7.
498 ! yliu: end
499 if(plfo(n).eq.99.) then
500 IF (iprt) print *,'n=',n,' unknown ob of type',platform
501 endif
502
503 !======================================================================
504 !======================================================================
505 ! THIS PART READS SOUNDING INFO
506 IF(is_sound)THEN
507 nlevs_ob(n)=real(meas_count)
508 lev_in_ob(n)=1.
509 do imc=1,meas_count
510 ! write(6,*) '0 inest = ',inest,' n = ',n
511 ! the sounding has one header, many levels. This part puts it into
512 ! "individual" observations. There's no other way for nudob to deal
513 ! with it.
514 if(imc.gt.1)then ! sub-loop over N
515 n=n+1
516 if(n.gt.niobf)goto 120
517 nlevs_ob(n)=real(meas_count)
518 lev_in_ob(n)=real(imc)
519 timeob(n)=rtimob
520 rio(n)=ri
521 rjo(n)=rj
522 rko(n)=-99.
523 plfo(n)=plfo(n-imc+1)
524 elevob(n)=elevation
525 endif
526
527 read(nvol,104)pressure_data,pressure_qc, &
528 height_data,height_qc, &
529 temperature_data,temperature_qc, &
530 u_met_data,u_met_qc, &
531 v_met_data,v_met_qc, &
532 rh_data,rh_qc
533 104 FORMAT( 1x,6(f11.3,1x,f11.3,1x))
534
535 ! yliu: Ensemble - add disturbance to upr obs
536 ! if(plfo(n).eq.5.or.plfo(n).eq.6.or.plfo(n).eq.9) then FORE07E08
537 ! if(imc.eq.1) then FORE07E08
538 ! call srand(n)
539 ! t_rand =- (rand(2)-0.5)*6
540 ! call srand(n+100000)
541 ! u_rand =- (rand(2)-0.5)*6
542 ! call srand(n+200000)
543 ! v_rand =- (rand(2)-0.5)*6
544 ! endif FORE07E08
545 ! if(temperature_qc.ge.0..and.temperature_qc.lt.30000..and.
546 ! & temperature_data .gt. -88880.0 )
547 ! & temperature_data = temperature_data + t_rand
548 ! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
549 ! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
550 ! make sure at least 1 of the components is .ne.0
551 ! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
552 ! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
553 ! u_met_data = u_met_data + u_rand
554 ! v_met_data = v_met_data + v_rand
555 ! endif
556 ! endif FORE07E08
557 ! yliu: Ens test - end
558
559
560 ! jc
561 ! hardwire to switch -777777. qc to 0. here temporarily
562 ! -777777. is a sounding level that no qc was done on.
563
564 if(temperature_qc.eq.-777777.)temperature_qc=0.
565 if(pressure_qc.eq.-777777.)pressure_qc=0.
566 if(height_qc.eq.-777777.)height_qc=0.
567 if(u_met_qc.eq.-777777.)u_met_qc=0.
568 if(v_met_qc.eq.-777777.)v_met_qc=0.
569 if(rh_qc.eq.-777777.)rh_qc=0.
570 if(temperature_data.eq.-888888.)temperature_qc=-888888.
571 if(pressure_data.eq.-888888.)pressure_qc=-888888.
572 if(height_data.eq.-888888.)height_qc=-888888.
573 if(u_met_data.eq.-888888.)u_met_qc=-888888.
574 if(v_met_data.eq.-888888.)v_met_qc=-888888.
575 if(rh_data.eq.-888888.)rh_qc=-888888.
576
577 ! jc
578 ! Hardwire so that only use winds in pilot obs (no winds from temp) and
579 ! only use temperatures and rh in temp obs (no temps from pilot obs)
580 ! Do this because temp and pilot are treated as 2 platforms, but pilot
581 ! has most of the winds, and temp has most of the temps. If use both,
582 ! the second will smooth the effect of the first. Usually temps come in after
583 ! pilots. pilots usually don't have any temps, but temp obs do have some
584 ! winds usually.
585 ! plfo=5 is TEMP ob, range sounding is an exception
586 !yliu start -- comment out to test with merged PILOT and TEMP and
587 ! do not use obs interpolated by little_r
588 ! if(plfo(n).eq.5. .and. namef(1:8).ne.'CLASS DA')then
589 ! u_met_data=-888888.
590 ! v_met_data=-888888.
591 ! u_met_qc=-888888.
592 ! v_met_qc=-888888.
593 ! endif
594 if(plfo(n).eq.5..and.(u_met_qc.eq.256..or.v_met_qc.eq.256.))then
595 u_met_data=-888888.
596 v_met_data=-888888.
597 u_met_qc=-888888.
598 v_met_qc=-888888.
599 endif
600 !yliu end
601 ! plfo=6 is PILOT ob
602 if(plfo(n).eq.6.)then
603 temperature_data=-888888.
604 rh_data=-888888.
605 temperature_qc=-888888.
606 rh_qc=-888888.
607 endif
608
609 !ajb Store potential temperature for WRF
610 if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
611
612 if(pressure_qc.ge.0..and.pressure_qc.lt.30000.)then
613
614 varobs(3,n) = &
615 temperature_data*(100000./pressure_data)**RCP - t0
616
617 ! write(6,*) 'reading data for N = ',n,' RCP = ',rcp
618 ! write(6,*) 'temperature_data = ',temperature_data
619 ! write(6,*) 'pressure_data = ',pressure_data
620 ! write(6,*) 'varobs(3,n) = ',varobs(3,n)
621
622 else
623 varobs(3,n)=-888888.
624 endif
625
626 else
627 varobs(3,n)=-888888.
628 endif
629
630 if(pressure_qc.ge.0..and.pressure_qc.lt.30000.)then
631 ! if(pressure_qc.ge.0.)then
632 varobs(5,n)=pressure_data
633 else
634 varobs(5,n)=-888888.
635 IF (iprt) THEN
636 print *,'********** PROBLEM *************'
637 print *,'sounding, p undefined',latitude,longitude
638 ENDIF
639 endif
640 if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
641 ! don't use data above 80 mb
642 if((varobs(5,n).gt.0.).and.(varobs(5,n).le.8.))then
643 u_met_data=-888888.
644 v_met_data=-888888.
645 u_met_qc=-888888.
646 v_met_qc=-888888.
647 temperature_data=-888888.
648 temperature_qc=-888888.
649 rh_data=-888888.
650 rh_qc=-888888.
651 endif
652
653 ! yliu: add special processing of NPN and Range profiler
654 ! only little_r interpolated and QC-ed data is used
655 if(namef(2:9).eq."PROFILER") then
656 if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
657 (v_met_qc.ge.0..and.v_met_qc.lt.30000.))then
658 !!yliu little_r already rotated the winds
659 ! call vect(longitude,u_met_data,v_met_data,xlonc,xlatc,xn)
660 varobs(1,n)=u_met_data
661 varobs(2,n)=v_met_data
662 else
663 varobs(1,n)=-888888.
664 varobs(2,n)=-888888.
665 endif
666 else
667 if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
668 (v_met_qc.ge.0..and.v_met_qc.lt.30000.))then
669 !!yliu little_r already rotated the winds
670 ! call vect(longitude,u_met_data,v_met_data,xlonc,xlatc,xn)
671 varobs(1,n)=u_met_data
672 varobs(2,n)=v_met_data
673 else
674 varobs(1,n)=-888888.
675 varobs(2,n)=-888888.
676 endif
677 endif
678 r_data=-888888.
679
680 if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
681 if((pressure_qc.ge.0.).and.(temperature_qc.ge.0.).and. &
682 (pressure_qc.lt.30000.).and.(temperature_qc.lt.30000.))then
683 call rh2r(rh_data,temperature_data,pressure_data*.01, &
684 r_data,0) ! yliu, change last arg from 1 to 0
685 else
686 ! print *,'rh, but no t or p to convert',temperature_qc, &
687 ! pressure_qc,n
688 r_data=-888888.
689 endif
690 endif
691 varobs(4,n)=r_data
692 enddo ! end do imc=1,meas_count
693 ! print *,'--- sdng n=',n,nlevs_ob(n),lev_in_ob(n),timeob(n)
694 ! read in non-sounding obs
695
696 ELSEIF(.NOT.is_sound)THEN
697 nlevs_ob(n)=1.
698 lev_in_ob(n)=1.
699 read(nvol,105)slp_data,slp_qc, &
700 ref_pres_data,ref_pres_qc, &
701 height_data,height_qc, &
702 temperature_data,temperature_qc, &
703 u_met_data,u_met_qc, &
704 v_met_data,v_met_qc, &
705 rh_data,rh_qc, &
706 psfc_data,psfc_qc, &
707 precip_data,precip_qc
708 105 FORMAT( 1x,9(f11.3,1x,f11.3,1x))
709
710 ! Ensemble: add disturbance to sfc obs
711 ! call srand(n)
712 ! t_rand =+ (rand(2)-0.5)*5
713 ! call srand(n+100000)
714 ! u_rand =+ (rand(2)-0.5)*5
715 ! call srand(n+200000)
716 ! v_rand =+ (rand(2)-0.5)*5
717 ! if(temperature_qc.ge.0..and.temperature_qc.lt.30000. .and.
718 ! & temperature_data .gt. -88880.0 )
719 ! & temperature_data = temperature_data + t_rand
720 ! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
721 ! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
722 ! make sure at least 1 of the components is .ne.0
723 ! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
724 ! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
725 ! u_met_data = u_met_data + u_rand
726 ! v_met_data = v_met_data + v_rand
727 ! endif
728 ! yliu: Ens test - end
729
730 !ajb Store potential temperature for WRF
731 if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
732
733 if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and. &
734 (psfc_data.gt.70000. .and.psfc_data.lt.105000.))then
735
736 varobs(3,n) = &
737 temperature_data*(100000./psfc_data)**RCP - t0
738 else
739 varobs(3,n)=-888888.
740 endif
741 else
742 varobs(3,n)=-888888.
743 endif
744
745 if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and.(psfc_data.gt.70000. &
746 .and.psfc_data.lt.105000.))then
747 varobs(5,n)=psfc_data
748 else
749 varobs(5,n)=-888888.
750 endif
751 if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
752
753 if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
754 (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and. &
755 ! make sure at least 1 of the components is .ne.0
756 (u_met_data.ne.0..or.v_met_data.ne.0.))then
757 !!yliu little_r already rotated the winds
758 !!yliu call vect(longitude,u_met_data,v_met_data,xlonc,xlatc,xn)
759 varobs(1,n)=u_met_data
760 varobs(2,n)=v_met_data
761 else
762 varobs(1,n)=-888888.
763 varobs(2,n)=-888888.
764 endif
765 ! calculate psfc if slp is there
766 if((psfc_qc.lt.0.).and.(slp_qc.ge.0..and.slp_qc.lt.30000.).and. &
767 (temperature_qc.ge.0..and.temperature_qc.lt.30000.).and. &
768 (slp_data.gt.90000.))then
769 tbar=temperature_data+0.5*elevation*.0065
770 psfc_data=slp_data*exp(-elevation/(rovg*tbar))
771 varobs(5,n)=psfc_data*1.e-3
772 psfc_qc=0.
773 endif
774
775 !c *No* **Very rough** estimate of psfc from sfc elevation if UUtah ob and elev>1000m
776 ! estimate psfc from temp and elevation
777 ! Do not know sfc pressure in model at this point.
778 ! if((psfc_qc.lt.0.).and.(elevation.gt.1000.).and.
779 ! 1 (temperature_qc.ge.0..and.temperature_qc.lt.30000.)
780 ! 1 .and.(platform(7:16).eq.'SYNOP PRET'))then
781 if((psfc_qc.lt.0.).and. &
782 (temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
783 tbar=temperature_data+0.5*elevation*.0065
784 psfc_data=100000.*exp(-elevation/(rovg*tbar))
785 varobs(5,n)=psfc_data*1.e-3
786 psfc_qc=0.
787 endif
788 ! jc
789 ! if a ship ob has rh<70%, then throw out
790 if(plfo(n).eq.3..and.rh_qc.ge.0..and.rh_data.lt.70.)then
791 rh_qc=-888888.
792 rh_data=-888888.
793 endif
794 !
795 r_data=-888888.
796 if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
797 if((psfc_qc.ge.0..and.psfc_qc.lt.30000.) &
798 .and.(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
799 ! rh_data=amin1(rh_data,96.) ! yliu: do not allow surface to be saturated
800 call rh2r(rh_data,temperature_data,psfc_data*.01, &
801 r_data,0) ! yliu, change last arg from 1 to 0
802 else
803 ! print *,'rh, but no t or p',temperature_data,
804 ! 1 psfc_data,n
805 r_data=-888888.
806 endif
807 endif
808 varobs(4,n)=r_data
809 ELSE
810 IF (iprt) THEN
811 print *,' ====== '
812 print *,' NO Data Found '
813 ENDIF
814 ENDIF !end if(is_sound)
815 ! END OF SFC OBS INPUT SECTION
816 !======================================================================
817 !======================================================================
818 ! check if ob time is too early (only applies to beginning)
819 IF(RTIMOB.LT.TBACK-fdob%WINDOW)then
820 IF (iprt) print *,'ob too early'
821 n=n-1
822 GOTO 110
823 ENDIF
824
825 ! check if this ob is a duplicate
826 ! this check has to be before other checks
827 njend=n-1
828 if(is_sound)njend=n-meas_count
829 do njc=1,njend
830 ! Check that time, lat, lon, and platform all match exactly.
831 ! Platforms 1-4 (surface obs) can match with each other. Otherwise,
832 ! platforms have to match exactly.
833 if( (timeob(n).eq.timeob(njc)) .and. &
834 (rio(n).eq.rio(njc)) .and. &
835 (rjo(n).eq.rjo(njc)) .and. &
836 (plfo(njc).ne.99.) ) then
837 !yliu: if two sfc obs are departed less than 1km, consider they are redundant
838 ! (abs(rio(n)-rio(njc))*dscg.gt.1000.) &
839 ! .or. (abs(rjo(n)-rjo(njc))*dscg.gt.1000.) &
840 ! .or. (plfo(njc).eq.99.) )goto 801
841 !yliu end
842 ! If platforms different, and either > 4, jump out
843 if( ( (plfo(n).le.4.).and.(plfo(njc).le.4.) ) .or. &
844 (plfo(n).eq.plfo(njc)) ) then
845
846 ! if not a sounding, and levels are the same then replace first occurrence
847 if((.not.is_sound).and.(rko(njc).eq.rko(n))) then
848 ! print *,'dup single ob-replace ',n,inest,
849 ! plfo(n),plfo(njc)
850 ! this is the sfc ob replacement part
851 VAROBS(1,njc)=VAROBS(1,n)
852 VAROBS(2,njc)=VAROBS(2,n)
853 VAROBS(3,njc)=VAROBS(3,n)
854 VAROBS(4,njc)=VAROBS(4,n)
855 VAROBS(5,njc)=VAROBS(5,n)
856 ! don't need to switch these because they're the same
857 ! RIO(njc)=RIO(n)
858 ! RJO(njc)=RJO(n)
859 ! RKO(njc)=RKO(n)
860 ! TIMEOB(njc)=TIMEOB(n)
861 ! nlevs_ob(njc)=nlevs_ob(n)
862 ! lev_in_ob(njc)=lev_in_ob(n)
863 ! plfo(njc)=plfo(n)
864 ! end sfc ob replacement part
865
866 n=n-1
867 goto 100
868 ! It's harder to fix the soundings, since the number of levels may be different
869 ! The easiest thing to do is to just set the first occurrence to all missing, and
870 ! keep the second occurrence, or vice versa.
871 ! For temp or profiler keep the second, for pilot keep the one with more levs
872 ! This is for a temp or prof sounding, equal to same
873 ! also if a pilot, but second one has more obs
874 elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
875 ( (plfo(njc).eq.5.).or.(plfo(njc).eq.9.).or. &
876 ( (plfo(njc).eq.6.).and. &
877 (nlevs_ob(n).ge.nlevs_ob(njc)) ) ) )then
878 IF (iprt) THEN
879 print *,'duplicate sounding - eliminate first occurrence', &
880 n,inest,meas_count,nlevs_ob(njc), &
881 latitude,longitude,plfo(njc)
882 ENDIF
883 if(lev_in_ob(njc).ne.1.) then
884 IF (iprt) THEN
885 print *, 'problem ******* - dup sndg ', &
886 lev_in_ob(njc),nlevs_ob(njc)
887 ENDIF
888 endif
889 ! n=n-meas_count
890 ! set the first sounding ob to missing
891 do njcc=njc,njc+nint(nlevs_ob(njc))-1
892 VAROBS(1,njcc)=-888888.
893 VAROBS(2,njcc)=-888888.
894 VAROBS(3,njcc)=-888888.
895 VAROBS(4,njcc)=-888888.
896 VAROBS(5,njcc)=-888888.
897 plfo(njcc)=99.
898 enddo
899 goto 100
900 ! if a pilot, but first one has more obs
901 elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
902 (plfo(njc).eq.6.).and. &
903 (nlevs_ob(n).lt.nlevs_ob(njc)) )then
904 IF (iprt) THEN
905 print *, &
906 'duplicate pilot sounding - eliminate second occurrence', &
907 n,inest,meas_count,nlevs_ob(njc), &
908 latitude,longitude,plfo(njc)
909 ENDIF
910 if(lev_in_ob(njc).ne.1.) then
911 IF (iprt) THEN
912 print *, 'problem ******* - dup sndg ', &
913 lev_in_ob(njc),nlevs_ob(njc)
914 ENDIF
915 endif
916 n=n-meas_count
917
918 !ajb Reset timeob for discarded indices.
919 do imc = n+1, n+meas_count
920 timeob(imc) = 99999.
921 enddo
922 goto 100
923 ! This is for a single-level satellite upper air ob - replace first
924 elseif( (is_sound).and. &
925 (nlevs_ob(njc).eq.1.).and. &
926 (nlevs_ob(n).eq.1.).and. &
927 (varobs(5,njc).eq.varobs(5,n)).and. &
928 (plfo(njc).eq.7.).and.(plfo(n).eq.7.) ) then
929 IF (iprt) print *, &
930 'duplicate single lev sat-wind ob - replace first',n, &
931 inest,meas_count,varobs(5,n)
932 ! this is the single ua ob replacement part
933 VAROBS(1,njc)=VAROBS(1,n)
934 VAROBS(2,njc)=VAROBS(2,n)
935 VAROBS(3,njc)=VAROBS(3,n)
936 VAROBS(4,njc)=VAROBS(4,n)
937 VAROBS(5,njc)=VAROBS(5,n)
938 ! don't need to switch these because they're the same
939 ! RIO(njc)=RIO(n)
940 ! RJO(njc)=RJO(n)
941 ! RKO(njc)=RKO(n)
942 ! TIMEOB(njc)=TIMEOB(n)
943 ! nlevs_ob(njc)=nlevs_ob(n)
944 ! lev_in_ob(njc)=lev_in_ob(n)
945 ! plfo(njc)=plfo(n)
946 ! end single ua ob replacement part
947 n=n-1
948 goto 100
949 else
950 IF (iprt) THEN
951 print *,'duplicate location, but no match otherwise',n,njc, &
952 plfo(n),varobs(5,n),nlevs_ob(n),lev_in_ob(n), &
953 plfo(njc),varobs(5,njc),nlevs_ob(njc),lev_in_ob(njc)
954 ENDIF
955 endif
956 endif
957 endif
958 ! end of njc do loop
959 enddo
960
961 ! check if ob is a sams ob that came in via UUtah - discard
962 if( plfo(n).eq.4..and.(platform(7:16).eq.'SYNOP PRET').and. &
963 (id(7:15).eq.'METNET= 3') )then
964 ! print *,'elim metnet=3',latitude,longitude,rtimob
965 n=n-1
966 goto 100
967 endif
968
969 ! check if ob is in coarse mesh domain (061404 switched sn/we)
970 if( (ri.lt.2.).or.(ri.gt.real(we_maxcg-1)).or.(rj.lt.2.).or. &
971 (rj.gt.real(sn_maxcg-1)) ) then
972
973 ! if (iprt) write(6,*) 'Obs out of coarse mesh domain'
974 ! write(6,*) 'we_maxcg-1 = ',real(we_maxcg-1)
975 ! write(6,*) 'sn_maxcg-1 = ',real(sn_maxcg-1)
976
977 ! n=n-1
978 ! if(is_sound)n=n-meas_count+1
979
980 n=n-meas_count
981 !ajb Reset timeob for discarded indices.
982 do imc = n+1, n+meas_count
983 timeob(imc) = 99999.
984 enddo
985 goto 100
986 endif
987
988 ! check if an upper air ob is too high
989 ! the ptop here is hardwired
990 ! this check has to come after other checks - usually the last few
991 ! upper air obs are too high
992 ! if(is_sound)then
993 ! njc=meas_count
994 ! do jcj=meas_count,1,-1
995 ! 6. is 60 mb - hardwired
996 ! if((varobs(5,n).lt.6.).and.(varobs(5,n).gt.0.))then
997 ! print *,'obs too high - eliminate,n,p=',n,varobs(5,n)
998 ! n=n-1
999 ! else
1000 ! if(varobs(5,n).gt.0.)goto 100
1001 ! endif
1002 ! enddo
1003 ! endif
1004 !
1005 IF(TIMEOB(N).LT.fdob%RTLAST) THEN
1006 IF (iprt) THEN
1007 PRINT *,'2 OBS ARE NOT IN CHRONOLOGICAL ORDER'
1008 PRINT *,'NEW YEAR?'
1009 print *,'timeob,rtlast,n=',timeob(n),fdob%rtlast,n
1010 ENDIF
1011 STOP 111
1012 ELSE
1013 fdob%RTLAST=TIMEOB(N)
1014 ENDIF
1015 GOTO 100
1016 111 CONTINUE
1017 !**********************************************************************
1018 ! -------------- END BIG 100 LOOP OVER N --------------
1019 !**********************************************************************
1020 IF (iprt) write(6,5403) NVOL,XTIME
1021 IEOF(inest)=1
1022
1023 close(NVOLA+INEST-1)
1024 IF (iprt) print *,'closed fdda file for inest=',inest,nsta
1025
1026 ! if(nsta.eq.1.and.timeob(1).gt.9.e4)nsta=0
1027 goto 1001
1028
1029 ! THE OBSERVATION ARRAYS ARE FULL AND THE MOST RECENTLY
1030 ! ACQUIRED OBS STILL HAS TIMEOB .LE. TFORWD. SO START
1031 ! DECREASING THE SIZE OF THE WINDOW
1032 ! get here if too many obs
1033 120 CONTINUE
1034 IF (iprt) THEN
1035 write(6,121) N,NIOBF
1036 write(6,122)
1037 ENDIF
1038 STOP 122
1039 fdob%WINDOW=fdob%WINDOW-0.1*TWINDO
1040 IF(TWINDO.LT.0)STOP 120
1041 ! IF THE WINDOW BECOMES NEGATIVE, THE INCOMING DATA IS
1042 ! PROBABLY GARBLED. STOP.
1043 GOTO 10
1044 !
1045 ! READ CYCLE IS COMPLETED. DETERMINE THE NUMBER OF OBS IN
1046 ! THE CURRENT WINDOW
1047 !
1048 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1049 ! BUT FIRST, WHEN KTAU.EQ.0 (OR IN GENERAL, KTAUR), DISCARD THE
1050 ! "OLD" OBS FIRST...
1051 130 CONTINUE
1052
1053 ! get here if at end of file, or if obs time is beyond what we
1054 ! need right now
1055 IF(KTAU.EQ.KTAUR)THEN
1056 NSTA=0
1057 keep_obs : DO N=1,NIOBF
1058
1059 ! try to keep all obs, but just don't use yet
1060 ! (don't want to throw away last obs read in - especially if
1061 ! its a sounding, in which case it looks like many obs)
1062 IF(TIMEOB(N).GT.9.e4) EXIT keep_obs
1063 if(timeob(n).gt.tforwd) then
1064 if(iprt) write(6,951)inest,n,timeob(n),tforwd
1065 951 FORMAT('saving ob beyond window,inest,n,timeob,tforwd=', &
1066 2i5,2f13.4)
1067 endif
1068 NSTA=N
1069 ENDDO keep_obs
1070
1071 NDUM=0
1072 ! make time=99999. if ob is too old
1073 ! print *,'tback,nsta=',tback,nsta
1074 old_obs : DO N=1,NSTA+1
1075 IF((TIMEOB(N)-TBACK).LT.0)THEN
1076 TIMEOB(N)=99999.
1077 ENDIF
1078 ! print *,'n,ndum,timeob=',n,ndum,timeob(n)
1079 IF(TIMEOB(N).LT.9.E4) EXIT old_obs
1080 NDUM=N
1081 ENDDO old_obs
1082
1083 ! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
1084 IF (iprt) THEN
1085 print *,'after 190 ndum=',ndum,nsta
1086 print *,'timeob=',timeob(1),timeob(2)
1087 ENDIF
1088 NDUM=ABS(NDUM)
1089 NMOVE=NIOBF-NDUM
1090 IF( NMOVE.GT.0 .AND. NDUM.NE.0) THEN
1091 DO N=1,NMOVE
1092 VAROBS(1,N)=VAROBS(1,N+NDUM)
1093 VAROBS(2,N)=VAROBS(2,N+NDUM)
1094 VAROBS(3,N)=VAROBS(3,N+NDUM)
1095 VAROBS(4,N)=VAROBS(4,N+NDUM)
1096 VAROBS(5,N)=VAROBS(5,N+NDUM)
1097 RJO(N)=RJO(N+NDUM)
1098 RIO(N)=RIO(N+NDUM)
1099 RKO(N)=RKO(N+NDUM)
1100 TIMEOB(N)=TIMEOB(N+NDUM)
1101 nlevs_ob(n)=nlevs_ob(n+ndum)
1102 lev_in_ob(n)=lev_in_ob(n+ndum)
1103 plfo(n)=plfo(n+ndum)
1104 ENDDO
1105 ENDIF
1106 ! moved obs up. now fill remaining space with 99999.
1107 NOPEN=NMOVE+1
1108 IF(NOPEN.LE.NIOBF) THEN
1109 DO N=NOPEN,NIOBF
1110 VAROBS(1,N)=99999.
1111 VAROBS(2,N)=99999.
1112 VAROBS(3,N)=99999.
1113 VAROBS(4,N)=99999.
1114 VAROBS(5,N)=99999.
1115 RIO(N)=99999.
1116 RJO(N)=99999.
1117 RKO(N)=99999.
1118 TIMEOB(N)=99999.
1119 ENDDO
1120 ENDIF
1121 ENDIF
1122 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1123 NSTA=0
1124 ! print *,'nsta at restart setting is ',nsta
1125 ! recalculate nsta after moving things around
1126 recalc : DO N=1,NIOBF
1127 ! try to save all obs - don't throw away latest read in
1128 IF(TIMEOB(N).GT.9.e4) EXIT recalc
1129 NSTA=N
1130 ! nsta=n-1 ! yliu test
1131 ENDDO recalc
1132
1133 IF (iprt) write(6,160) KTAU,XTIME,NSTA
1134 IF(KTAU.EQ.KTAUR)THEN
1135 IF(nudge_opt.EQ.1)THEN
1136 TWDOP=TWINDO*60.
1137 IF (iprt) THEN
1138 write(6,1449) INEST,RINXY,RINSIG,TWDOP
1139 IF(ISWIND.EQ.1) write(6,1450) GIV
1140 IF(ISTEMP.EQ.1) write(6,1451) GIT
1141 IF(ISMOIS.EQ.1) write(6,1452) GIQ
1142 IF(ISPSTR.EQ.1) write(6,1453) GIP
1143 ENDIF
1144 ENDIF
1145 ENDIF
1146 IF(KTAU.EQ.KTAUR)THEN
1147 IF (iprt) THEN
1148 write(6,553)
1149 write(6,554)
1150 ENDIF
1151 IF(fdob%IWTSIG.NE.1)THEN
1152 IF (iprt) THEN
1153 write(6,555)
1154 write(6,556) fdob%RINFMN*RINXY,fdob%RINFMX*RINXY,fdob%PFREE*10.
1155 ENDIF
1156 IF(fdob%RINFMN.GT.fdob%RINFMX)STOP 556
1157 ! IS MINIMUM GREATER THAN MAXIMUM?
1158 IF (iprt) write(6,557) fdob%DPSMX*10.,fdob%DCON
1159 IF(fdob%DPSMX.GT.10.)STOP 557
1160 ENDIF
1161 ENDIF
1162 ! IS DPSMX IN CB?
1163
1164 IF(KTAU.EQ.KTAUR)THEN
1165 IF (iprt) write(6,601) INEST,IONF
1166 ENDIF
1167 fdob%NSTAT=NSTA
1168
1169 555 FORMAT(1X,' ABOVE THE SURFACE LAYER, OBS NUDGING IS PERFORMED', &
1170 ' ON PRESSURE LEVELS,')
1171 556 FORMAT(1X,' WHERE RINXY VARIES LINEARLY FROM ',E11.3,' KM AT', &
1172 ' THE SURFACE TO ',E11.3,' KM AT ',F7.2,' MB AND ABOVE')
1173 557 FORMAT(1X,' IN THE SURFACE LAYER, WXY IS A FUNCTION OF ', &
1174 'DPSMX = ',F7.2,' MB WITH DCON = ',E11.3, &
1175 ' - SEE SUBROUTINE NUDOB')
1176 601 FORMAT('0','FOR EFFICIENCY, THE OBS NUDGING FREQUENCY ', &
1177 'FOR MESH #',I2,' IS ',1I2,' CGM TIMESTEPS ',/)
1178 121 FORMAT('0',' WARNING: NOBS = ',I4,' IS GREATER THAN NIOBF = ', &
1179 I4,': INCREASE PARAMETER NIOBF')
1180 5403 FORMAT(1H0,'-------------EOF REACHED FOR NVOL = ',I3, &
1181 ' AND XTIME = ',F10.2,'-------------------')
1182 122 FORMAT(1X,' ...OR THE CODE WILL REDUCE THE TIME WINDOW')
1183 160 FORMAT('0','****** CALL IN4DOB AT KTAU = ',I5,' AND XTIME = ', &
1184 F10.2,': NSTA = ',I7,' ******')
1185 1449 FORMAT(1H0,'*****NUDGING INDIVIDUAL OBS ON MESH #',I2, &
1186 ' WITH RINXY = ', &
1187 E11.3,' KM, RINSIG = ',E11.3,' AND TWINDO (HALF-PERIOD) = ', &
1188 E11.3,' MIN')
1189 1450 FORMAT(1X,'NUDGING IND. OBS WINDS WITH GIV = ',E11.3)
1190 1451 FORMAT(1X,'NUDGING IND. OBS TEMPERATURE WITH GIT = ',E11.3)
1191 1452 FORMAT(1X,'NUDGING IND. OBS MOISTURE WITH GIQ = ',E11.3)
1192 1453 FORMAT(1X,'NUDGING IND. OBS SURFACE PRESSURE WITH GIP = ,'E11.3)
1193 553 FORMAT(1X,'BY DEFAULT: OBS NUDGING OF TEMPERATURE AND MOISTURE ', &
1194 'IS RESTRICTED TO ABOVE THE BOUNDARY LAYER')
1195 554 FORMAT(1X,'...WHILE OBS NUDGING OF WIND IS INDEPENDENT OF THE ', &
1196 'BOUNDARY LAYER')
1197
1198 RETURN
1199 END SUBROUTINE in4dob
1200
1201 SUBROUTINE julgmt(mdate,julgmtn,timanl,julday,gmt,ind)
1202 ! CONVERT MDATE YYMMDDHH TO JULGMT (JULIAN DAY * 100. +GMT)
1203 ! AND TO TIMANL (TIME IN MINUTES WITH RESPECT TO MODEL TIME)
1204 ! IF IND=0 INPUT MDATE, OUTPUT JULGMTN AND TIMANL
1205 ! IF IND=1 INPUT TIMANL, OUTPUT JULGMTN
1206 ! IF IND=2 INPUT JULGMTN, OUTPUT TIMANL
1207 INTEGER, intent(in) :: MDATE
1208 REAL, intent(out) :: JULGMTN
1209 REAL, intent(out) :: TIMANL
1210 INTEGER, intent(in) :: JULDAY
1211 REAL, intent(in) :: GMT
1212 INTEGER, intent(in) :: IND
1213
1214 !*** DECLARATIONS FOR IMPLICIT NONE
1215 real :: MO(12), rjulanl, houranl, rhr
1216
1217 integer :: iyr, idate1, imo, idy, ihr, my1, my2, my3, ileap
1218 integer :: juldayn, juldanl, idymax, mm
1219
1220
1221 IF(IND.EQ.2)GOTO 150
1222 IYR=INT(MDATE/1000000.+0.001)
1223 IDATE1=MDATE-IYR*1000000
1224 IMO=INT(IDATE1/10000.+0.001)
1225 IDY=INT((IDATE1-IMO*10000.)/100.+0.001)
1226 IHR=IDATE1-IMO*10000-IDY*100
1227 MO(1)=31
1228 MO(2)=28
1229 ! IS THE YEAR A LEAP YEAR? (IN THIS CENTURY)
1230 IYR=IYR+1900
1231 MY1=MOD(IYR,4)
1232 MY2=MOD(IYR,100)
1233 MY3=MOD(IYR,400)
1234 ILEAP=0
1235 ! jc
1236 ! IF(MY1.EQ.0.AND.MY2.NE.0.OR.MY3.EQ.0)THEN
1237 IF(MY1.EQ.0)THEN
1238 ILEAP=1
1239 MO(2)=29
1240 ENDIF
1241 IF(IND.EQ.1)GOTO 200
1242 MO(3)=31
1243 MO(4)=30
1244 MO(5)=31
1245 MO(6)=30
1246 MO(7)=31
1247 MO(8)=31
1248 MO(9)=30
1249 MO(10)=31
1250 MO(11)=30
1251 MO(12)=31
1252 JULDAYN=0
1253 DO 100 MM=1,IMO-1
1254 JULDAYN=JULDAYN+MO(MM)
1255 100 CONTINUE
1256
1257 IF(IHR.GE.24)THEN
1258 IDY=IDY+1
1259 IHR=IHR-24
1260 ENDIF
1261 JULGMTN=(JULDAYN+IDY)*100.+IHR
1262 ! CONVERT JULGMT TO TIMANL WRT MODEL TIME IN MINUTES (XTIME)
1263 150 CONTINUE
1264 JULDANL=INT(JULGMTN/100.+0.000001)
1265 RJULANL=FLOAT(JULDANL)*100.
1266 HOURANL=JULGMTN-RJULANL
1267 TIMANL=(FLOAT(JULDANL-JULDAY)*24.-GMT+HOURANL)*60.
1268 RETURN
1269 200 CONTINUE
1270 RHR=GMT+TIMANL/60.+0.000001
1271 IDY=JULDAY
1272 IDYMAX=365+ILEAP
1273 300 IF(RHR.GE.24.0)THEN
1274 RHR=RHR-24.0
1275 IDY=IDY+1
1276 GOTO 300
1277 ENDIF
1278 IF(IDY.GT.IDYMAX)IDY=IDY-IDYMAX
1279 JULGMTN=FLOAT(IDY)*100.+RHR
1280 RETURN
1281 END SUBROUTINE julgmt
1282
1283 SUBROUTINE vect(xlon,e1,e2,xlonc,xlatc,xn)
1284
1285 ! THIS ROUTINE CONVERTS INCOMING U AND V COMPS INTO MAP U AND V COMPS.
1286 ! iproj is projection (1=lamconf, 2=polarst, 3=mercator)
1287 ! xlonc is center longitude
1288 ! xn is cone factor (.716 for current lc)
1289 !
1290 REAL, intent(in) :: XLON
1291 REAL, intent(inout) :: E1
1292 REAL, intent(inout) :: E2
1293 REAL, intent(in) :: xlonc
1294 REAL, intent(in) :: xlatc
1295 REAL, intent(in) :: xn
1296
1297 !*** DECLARATIONS FOR IMPLICIT NONE
1298 real :: pi, degran, u, v, xlonr, angle
1299
1300 pi=3.1415926535
1301 DEGRAN=PI/180.
1302 !
1303 !
1304 u=e1
1305 v=e2
1306 XLONR=XLONC-XLON
1307 IF(XLONR.GT.180.) XLONR=XLONR-360.
1308 IF(XLONR.LT.-180.) XLONR=XLONR+360.
1309 ANGLE=XLONR*XN*DEGRAN
1310 IF (xlatC.LT.0.0) ANGLE=-ANGLE
1311 E1=V*SIN(ANGLE)+U*COS(ANGLE)
1312 E2=V*COS(ANGLE)-U*SIN(ANGLE)
1313 RETURN
1314 END SUBROUTINE vect
1315
1316 SUBROUTINE rh2r(rh,t,p,r,iice)
1317
1318 ! convert rh to r
1319 ! if iice=1, use saturation with respect to ice
1320 ! rh is 0-100.
1321 ! r is g/g
1322 ! t is K
1323 ! p is mb
1324 !
1325 REAL, intent(in) :: rh
1326 REAL, intent(in) :: t
1327 REAL, intent(in) :: p
1328 REAL, intent(out) :: r
1329 INTEGER, intent(in) :: iice
1330
1331 !*** DECLARATIONS FOR IMPLICIT NONE
1332 real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
1333 real esat, rsat
1334
1335 eps=0.62197
1336 e0=6.1078
1337 eslcon1=17.2693882
1338 eslcon2=35.86
1339 esicon1=21.8745584
1340 esicon2=7.66
1341 t0=260.
1342
1343 ! print *,'rh2r input=',rh,t,p
1344 rh1=rh*.01
1345
1346 if(iice.eq.1.and.t.le.t0)then
1347 esat=e0*exp(esicon1*(t-273.16)/(t-esicon2))
1348 else
1349 esat=e0*exp(eslcon1*(t-273.16)/(t-eslcon2))
1350 endif
1351 rsat=eps*esat/(p-esat)
1352 ! print *,'rsat,esat=',rsat,esat
1353 r=rh1*rsat
1354
1355 ! print *,'rh2r rh,t,p,r=',rh1,t,p,r
1356
1357 return
1358 END SUBROUTINE rh2r
1359
1360 SUBROUTINE rh2rb(rh,t,p,r,iice)
1361
1362 ! convert rh to r
1363 ! if iice=1, use daturation with respect to ice
1364 ! rh is 0-100.
1365 ! r is g/g
1366 ! t is K
1367 ! p is mb
1368
1369 REAL, intent(in) :: rh
1370 REAL, intent(in) :: t
1371 REAL, intent(in) :: p
1372 REAL, intent(out) :: r
1373 INTEGER, intent(in) :: iice
1374
1375 !*** DECLARATIONS FOR IMPLICIT NONE
1376 real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
1377 real esat, rsat
1378
1379 eps=0.622
1380 e0=6.112
1381 eslcon1=17.67
1382 eslcon2=29.65
1383 esicon1=22.514
1384 esicon2=6.15e3
1385 t0=273.15
1386
1387 print *,'rh2r input=',rh,t,p
1388 rh1=rh*.01
1389
1390 if(iice.eq.1.and.t.le.t0)then
1391 esat=e0*exp(esicon1-esicon2/t)
1392 else
1393 esat=e0*exp(eslcon1*(t-t0)/(t-eslcon2))
1394 endif
1395 rsat=eps*esat/(p-esat)
1396 ! print *,'rsat,esat=',rsat,esat
1397 r=rh1*eps*rsat/(eps+rsat*(1.-rh1))
1398
1399 print *,'rh2r rh,t,p,r=',rh1,t,p,r
1400
1401 return
1402 END SUBROUTINE rh2rb
1403
1404 SUBROUTINE llxy_lam (xlat,xlon,X,Y,xlatc, xlonc,xn,ds, &
1405 imax, jmax, true_lat1, true_lat2 )
1406
1407 !*** DECLARATIONS FOR IMPLICIT NONE
1408 real :: pi, conv, a
1409
1410 PARAMETER(pi=3.14159,CONV=180./pi,a =6370.)
1411
1412 REAL TRUE_LAT1 , PHI1 , POLE , XLATC , PHIC , XC , YC , &
1413 XN , FLP , XLON , XLONC , PSX , XLAT , R , XX , YY , &
1414 CENTRI , CENTRJ , X , DS , Y , TRUE_LAT2
1415
1416 INTEGER IMAX , JMAX
1417
1418 ! Calculate x and y given latitude and longitude for Lambert conformal projection
1419
1420 IF(ABS(true_lat1).GT.90) THEN
1421 PHI1 = 90. - 30.
1422 ELSE
1423 PHI1 = 90. - true_lat1
1424 END IF
1425 POLE = 90.0
1426 IF ( XLATC.LT.0.0 ) THEN
1427 PHI1 = -PHI1
1428 POLE = -POLE
1429 END IF
1430 PHIC = ( POLE - XLATC )/CONV
1431 PHI1 = PHI1/CONV
1432 XC = 0.0
1433 YC = -A/XN*SIN(PHI1)*(TAN(PHIC/2.0)/TAN(PHI1/2.0))**XN
1434
1435 ! CALCULATE X,Y COORDS. RELATIVE TO POLE
1436
1437 FLP = XN*( XLON - XLONC )/CONV
1438 PSX = ( POLE - XLAT )/CONV
1439 R = -A/XN*SIN(PHI1)*(TAN(PSX/2.0)/TAN(PHI1/2.0))**XN
1440 IF ( XLATC.LT.0.0 ) THEN
1441 XX = R*SIN(FLP)
1442 YY = R*COS(FLP)
1443 ELSE
1444 XX = -R*SIN(FLP)
1445 YY = R*COS(FLP)
1446 END IF
1447
1448 ! TRANSFORM (1,1) TO THE ORIGIN
1449
1450 CENTRI = FLOAT(IMAX + 1)/2.0
1451 CENTRJ = FLOAT(JMAX + 1)/2.0
1452 X = ( XX - XC )/DS + CENTRJ
1453 Y = ( YY - YC )/DS + CENTRI
1454
1455 return
1456 END SUBROUTINE llxy_lam
1457
1458 SUBROUTINE llxy(xlat,xlon,x,y,xlatc,xlonc,kproj,psi1,psi2,ds, &
1459 xn,sn_max,we_max,parent_grid_ratio, &
1460 i_parent_start,j_parent_start)
1461
1462 IMPLICIT NONE
1463
1464 ! CALCULATE X AND Y GIVEN LATITUDE AND LONGITUDE.
1465 ! PETER HOWELLS, NCAR, 1984
1466
1467 REAL, intent(in) :: xlat
1468 REAL, intent(inout) :: xlon
1469 REAL, intent(out) :: x
1470 REAL, intent(out) :: y
1471 REAL, intent(in) :: xlatc
1472 REAL, intent(in) :: xlonc
1473 INTEGER, intent(in) :: kproj
1474 REAL, intent(in) :: psi1
1475 REAL, intent(in) :: psi2
1476 REAL, intent(in) :: ds
1477 REAL, intent(in) :: xn
1478 INTEGER, intent(in) :: sn_max
1479 INTEGER, intent(in) :: we_max
1480 INTEGER, intent(in) :: parent_grid_ratio
1481 INTEGER, intent(in) :: i_parent_start
1482 INTEGER, intent(in) :: j_parent_start
1483
1484 !*** DECLARATIONS FOR IMPLICIT NONE
1485 real conv, a, phi1, pole, c2, xc, phicr, cell, yc, xlatr
1486 real phic, xx, yy, centri, centrj, ylon, flp, psx, r
1487 integer imax, jmax, imapst, jmapst
1488
1489 ! write(6,*) 'enter llxy'
1490 ! write(6,*) 'enter llxy: xlatc = ',xlatc,' xlonc = ',xlonc
1491 ! write(6,*) 'xlat = ',xlat,' xlon = ',xlon
1492 ! write(6,*) 'psi1 = ',psi1,' psi2 = ',psi2
1493 ! write(6,*) 'xn = ',xn,' kproj = ',kproj,' ds = ',ds
1494 ! write(6,*) 'sn_max = ',sn_max,' we_max = ',we_max
1495 ! write(6,*) 'parent_grid_ratio = ',parent_grid_ratio
1496 ! write(6,*) 'i_parent_start = ',i_parent_start
1497 ! write(6,*) 'j_parent_start = ',j_parent_start
1498
1499 conv = 57.29578
1500 a = 6370.0
1501 ! imax = sn_max*parent_grid_ratio+1
1502 ! jmax = we_max*parent_grid_ratio+1
1503 imax = sn_max*parent_grid_ratio !ajb for WRF
1504 jmax = we_max*parent_grid_ratio !ajb for WRF
1505 imapst= (j_parent_start-1)*parent_grid_ratio+1
1506 jmapst= (i_parent_start-1)*parent_grid_ratio+1
1507 phi1 = 90.0-psi2
1508 pole = 90.0
1509
1510 if ( xlatc.lt.0.0 ) then
1511 phi1 = -90.0-psi2
1512 pole = -pole
1513 endif
1514
1515 if (kproj.eq.3) then
1516 ! MERCATOR PROJECTION
1517 C2 = A*COS(PSI1)
1518 XC = 0.0
1519 PHICR = XLATC/CONV
1520 CELL = COS(PHICR)/(1.0+SIN(PHICR))
1521 YC = - C2*ALOG(CELL)
1522 IF (XLAT.NE.-90.) THEN
1523 XLATR = XLAT/CONV
1524 CELL = COS(XLATR)/(1.0+SIN(XLATR))
1525 YY = -C2*ALOG(CELL)
1526 IF (XLONC.LT.0.0) THEN
1527 IF (XLON.GT.0.0) XLON=XLON-360.
1528 ELSE
1529 IF (XLON.LT.0.0) XLON=360.+XLON
1530 ENDIF
1531 XX = C2*(XLON-XLONC)/CONV
1532 ENDIF
1533
1534 ELSE IF (KPROJ.EQ.1) THEN
1535 ! LAMBERT-COMFORMAL or POLAR-STEREO PROJECTION
1536 PHIC = ( POLE - XLATC )/CONV
1537 PHI1 = PHI1/CONV
1538 XC = 0.0
1539 YC = -A/XN*SIN(PHI1)*(TAN(PHIC/2.0)/TAN(PHI1/2.0))**XN
1540
1541 ! CALCULATE X,Y COORDS. RELATIVE TO POLE
1542
1543 YLON = XLON - XLONC
1544 IF(YLON.GT.180) YLON = YLON - 360.
1545 IF(YLON.LT.-180) YLON = YLON + 360.
1546 FLP = XN*YLON/CONV
1547 PSX = ( POLE - XLAT )/CONV
1548 R = -A/XN*SIN(PHI1)*(TAN(PSX/2.0)/TAN(PHI1/2.0))**XN
1549 IF ( XLATC.LT.0.0 ) THEN
1550 XX = R*SIN(FLP)
1551 YY = R*COS(FLP)
1552 ELSE
1553 XX = -R*SIN(FLP)
1554 YY = R*COS(FLP)
1555 END IF
1556 END IF
1557
1558 ! TRANSFORM (1,1) TO THE ORIGIN
1559
1560 CENTRI = FLOAT(IMAX + 1)/2.0
1561 CENTRJ = FLOAT(JMAX + 1)/2.0
1562 X = ( XX - XC )/DS + CENTRJ - jmapst + 1
1563 Y = ( YY - YC )/DS + CENTRI - imapst + 1
1564
1565 RETURN
1566 END SUBROUTINE llxy
1567
1568 SUBROUTINE llxy_try(xlat,xlon,x,y,xlatc,xlonc,kproj,psi1a,psi2,ds, &
1569 xn,sn_max,we_max,parent_grid_ratio, &
1570 i_parent_start,j_parent_start)
1571
1572 IMPLICIT NONE
1573
1574 ! CALCULATE X AND Y GIVEN LATITUDE AND LONGITUDE.
1575 ! PETER HOWELLS, NCAR, 1984
1576
1577 REAL, intent(in) :: xlat
1578 REAL, intent(inout) :: xlon
1579 REAL, intent(out) :: x
1580 REAL, intent(out) :: y
1581 REAL, intent(in) :: xlatc
1582 REAL, intent(in) :: xlonc
1583 INTEGER, intent(in) :: kproj
1584 REAL, intent(in) :: psi1a
1585 REAL, intent(in) :: psi2
1586 REAL, intent(in) :: ds
1587 REAL, intent(in) :: xn
1588 INTEGER, intent(in) :: sn_max
1589 INTEGER, intent(in) :: we_max
1590 INTEGER, intent(in) :: parent_grid_ratio
1591 INTEGER, intent(in) :: i_parent_start
1592 INTEGER, intent(in) :: j_parent_start
1593
1594 !*** DECLARATIONS FOR IMPLICIT NONE
1595 real conv, a, phi1, pole, c2, xc, phicr, cell, yc, xlatr
1596 real phic, xx, yy, centri, centrj, ylon, flp, psx, r
1597 integer imax, jmax, imapst, jmapst
1598 integer jmxc,imxc,iratio
1599 real ric0,rjc0,rs,yind,xind,rix,rjx,psi1
1600
1601 ! write(6,*) 'enter llxy'
1602 ! write(6,*) 'enter llxy: xlatc = ',xlatc,' xlonc = ',xlonc
1603 ! write(6,*) 'psi1 = ',psi1a,' psi2 = ',psi2
1604 ! write(6,*) 'xn = ',xn,' kproj = ',kproj,' ds = ',ds
1605 ! write(6,*) 'sn_max = ',sn_max,' we_max = ',we_max
1606 ! write(6,*) 'parent_grid_ratio = ',parent_grid_ratio
1607 ! write(6,*) 'i_parent_start = ',i_parent_start
1608 ! write(6,*) 'j_parent_start = ',j_parent_start
1609
1610 conv = 57.29578
1611 a = 6370.0
1612 imax = sn_max*parent_grid_ratio+1
1613 jmax = we_max*parent_grid_ratio+1
1614 imapst= (j_parent_start-1)*parent_grid_ratio+1
1615 jmapst= (i_parent_start-1)*parent_grid_ratio+1
1616 iratio = parent_grid_ratio
1617 imxc = imax
1618 jmxc = jmax
1619 rix=imapst
1620 rjx=jmapst
1621
1622 phi1 = psi1a
1623 if(xlatc .gt. 0.) then
1624 psi1=90.-phi1
1625 else
1626 psi1=-90+abs(phi1)
1627 endif
1628
1629 ric0=(imxc+1.)*0.5
1630 rjc0=(jmxc+1.)*0.5
1631
1632 ! if(kproj .eq. 1) then
1633 if(kproj .le. 2) then
1634 if(xlatc .gt. 0.) then
1635 rs=a/xn*sin(psi1/conv)*(tan((90.-xlat)/conv/2.) &
1636 /tan(psi1/2./conv))**xn
1637 yc=-a/xn*sin(psi1/conv)*(tan((90.-xlatc)/conv/2.) &
1638 /tan(psi1/2./conv))**xn
1639 else
1640 rs=a/xn*sin(psi1/conv)*(tan((-90.-xlat)/conv/2.) &
1641 /tan(psi1/2./conv))**xn
1642 yc=-a/xn*sin(psi1/conv)*(tan((-90.-xlatc)/conv/2.) &
1643 /tan(psi1/2./conv))**xn
1644 endif
1645 ! elseif(kproj .eq. 2) then
1646 ! if(xlatc .gt. 0.) then
1647 ! rs=a*sin((90.-xlat)/conv)*(1.+cos(psi1/conv))
1648 ! & /(1.+cos((90.-xlat)/conv))
1649 ! yc=-a*sin((90.-xlatc)/conv)*(1.+cos(psi1/conv))
1650 ! & /(1.+cos((90.-xlatc)/conv))
1651 ! else
1652 ! rs=a*sin((-90.-xlat)/conv)*(1.+cos(psi1/conv))
1653 ! & /(1.+cos((-90.-xlat)/conv))
1654 ! yc=-a*sin((-90.-xlatc)/conv)*(1.+cos(psi1/conv))
1655 ! & /(1.+cos((-90.-xlatc)/conv))
1656 ! endif
1657 else
1658 psi1 = 0 ! yliu added
1659 c2=a*cos(psi1/conv)
1660 yc=c2*alog((1.+sin(xlatc/conv))/cos(xlatc/conv))
1661 endif
1662
1663 if(kproj .le. 2) then
1664 y=(ric0-(yc/ds+rs/ds*cos(xn*(xlon-xlonc)/conv))-rix) &
1665 *iratio+1.0
1666 if(xlatc .gt. 0.) then
1667 x=(rjc0+rs/ds*sin(xn*(xlon-xlonc)/conv)-rjx)*iratio+1.0
1668 else
1669 x=(rjc0-rs/ds*sin(xn*(xlon-xlonc)/conv)-rjx)*iratio+1.0
1670 endif
1671 else
1672 y=c2*alog((1.+sin(xlat/conv))/cos(xlat/conv))
1673 x=c2*(xlon-xlonc)/conv
1674 y=(ric0+(y-yc)/ds-rix)*iratio+1.0
1675 x=(rjc0+x/ds-rjx)*iratio+1.0
1676 endif
1677 write(6,*) 'xj = ',x, 'yi = ',y,"xlat",xlat,"xlon",xlon
1678
1679 RETURN
1680 END SUBROUTINE llxy_try
1681
1682 #endif
1683 !-----------------------------------------------------------------------
1684 ! End subroutines for in4dob
1685 !-----------------------------------------------------------------------