module_fddaobs_rtfdda.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_fddaobs_rtfdda
4
5 ! This obs-nudging FDDA module (RTFDDA) is developed by the
6 ! NCAR/RAL/NSAP (National Security Application Programs), under the
7 ! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is
8 ! acknowledged for releasing this capability for WRF community
9 ! research applications.
10 !
11 ! The NCAR/RAL RTFDDA module was adapted, and significantly modified
12 ! from the obs-nudging module in the standard MM5V3.1 which was originally
13 ! developed by PSU (Stauffer and Seaman, 1994).
14 !
15 ! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module
16 ! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
17 ! Nov. 2006
18 !
19 ! References:
20 !
21 ! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
22 ! implementation of obs-nudging-based FDDA into WRF for supporting
23 ! ATEC test operations. 2005 WRF user workshop. Paper 10.7.
24 !
25 ! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update
26 ! on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE
27 ! and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7.
28
29 !
30 ! Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data
31 ! assimilation. J. Appl. Meteor., 33, 416-434.
32 !
33 ! http://www.rap.ucar.edu/projects/armyrange/references.html
34 !
35
36 CONTAINS
37
38 !------------------------------------------------------------------------------
39 SUBROUTINE fddaobs_init(obs_nudge_opt, maxdom, inest, parid, &
40 dx_coarse, restart, obs_twindo, itimestep, &
41 e_sn, s_sn_cg, e_sn_cg, s_we_cg, e_we_cg, &
42 #if ( EM_CORE == 1 )
43 fdob, &
44 #endif
45 ids,ide, jds,jde, kds,kde, &
46 ims,ime, jms,jme, kms,kme, &
47 its,ite, jts,jte, kts,kte)
48 !-----------------------------------------------------------------------
49 ! This routine does initialization for real time fdda obs-nudging.
50 !
51 !-----------------------------------------------------------------------
52 USE module_domain
53 !-----------------------------------------------------------------------
54 IMPLICIT NONE
55 !-----------------------------------------------------------------------
56
57 !=======================================================================
58 ! Definitions
59 !-----------
60 INTEGER, intent(in) :: maxdom
61 INTEGER, intent(in) :: obs_nudge_opt(maxdom)
62 INTEGER, intent(in) :: ids,ide, jds,jde, kds,kde, &
63 ims,ime, jms,jme, kms,kme, &
64 its,ite, jts,jte, kts,kte
65 INTEGER, intent(in) :: inest
66 INTEGER, intent(in) :: parid(maxdom)
67 REAL ,intent(in) :: dx_coarse ! coarse-domain grid cell-size (km)
68 LOGICAL, intent(in) :: restart
69 REAL, intent(inout) :: obs_twindo
70 INTEGER, intent(in) :: itimestep
71 INTEGER, intent(in) :: e_sn ! ending south-north grid index
72 INTEGER, intent(in) :: s_sn_cg ! starting south-north coarse-grid index
73 INTEGER, intent(in) :: e_sn_cg ! ending south-north coarse-grid index
74 INTEGER, intent(in) :: s_we_cg ! starting west-east coarse-grid index
75 INTEGER, intent(in) :: e_we_cg ! ending west-east coarse-grid index
76 #if ( EM_CORE == 1 )
77 TYPE(fdob_type), intent(inout) :: fdob
78 #endif
79
80 ! Local variables
81 logical :: nudge_flag ! nudging flag for this nest
82 integer :: ktau ! current timestep
83 integer :: nest ! loop counter
84 integer :: idom ! domain id
85 integer :: parent ! parent domain
86
87 #if ( EM_CORE == 1 )
88 ! This routine should only be called once. This is a check to make
89 ! certain that initialization only happens once.
90 if (fdob%domain_init .ne. 1) then
91 ! Obs-nudging will be initialized on this call
92 fdob%domain_init = 1
93 else
94 ! Obs-nudging has already been initialized, so return
95 return
96 endif
97
98 ! Set flag for nudging on pressure (not sigma) surfaces
99 fdob%iwtsig = 0
100
101 ! Set ending nudging date (used with dynamic ramp-down) to zero.
102 fdob%datend = 0.
103
104 ! Convert twindo from minutes to hours.
105 obs_twindo = obs_twindo / 60.
106
107 ! Initialize flags.
108
109 fdob%domain_tot=0
110 do nest=1,maxdom
111 fdob%domain_tot = fdob%domain_tot + obs_nudge_opt(nest)
112 end do
113
114 ! Set parameters.
115
116 fdob%pfree = 50.0
117 fdob%rinfmn = 1.0
118 fdob%rinfmx = 2.0
119 fdob%dpsmx = 7.5
120 fdob%dcon = 1.0/fdob%dpsmx
121 fdob%xn = 0.7155668 ! cone factor
122
123 fdob%ds_cg = dx_coarse / 1000. ! coarse gridsize (km)
124 fdob%sn_maxcg = e_sn_cg - s_sn_cg + 1 ! coarse domain grid dimension in N-S
125 fdob%we_maxcg = e_we_cg - s_we_cg + 1 ! coarse domain grid dimension in W-E
126 fdob%sn_end = e_sn - 1 ! ending S-N grid coordinate
127
128 ! Calculate the nest levels, levidn. Note that each nest
129 ! must know the nest levels levidn(maxdom) of each domain.
130 do nest=1,maxdom
131
132 ! Initialize nest level for each domain.
133 if (nest .eq. 1) then
134 fdob%levidn(nest) = 0 ! Mother domain has nest level 0
135 else
136 fdob%levidn(nest) = 1 ! All other domains start with 1
137 endif
138 idom = nest
139 100 parent = parid(idom) ! Go up the parent tree
140 if (parent .gt. 1) then ! If not up to mother domain
141 fdob%levidn(nest) = fdob%levidn(nest) + 1
142 idom = parid(parent)
143 goto 100
144 endif
145 enddo
146
147 ! Check to see if the nudging flag has been set. If not,
148 ! simply RETURN.
149 nudge_flag = (obs_nudge_opt(inest) .eq. 1)
150 if (.not. nudge_flag) return
151
152 ktau = itimestep
153 if(restart) then
154 fdob%ktaur = ktau
155 else
156 fdob%ktaur = 0
157 endif
158
159 RETURN
160 #endif
161 END SUBROUTINE fddaobs_init
162
163 #if ( EM_CORE == 1 )
164 !-----------------------------------------------------------------------
165 SUBROUTINE errob(inest, ub, vb, tb, t0, qvb, pbase, pp, rovcp, &
166 uratx, vratx, tratx, nndgv, &
167 nerrf, niobf, maxdom, levidn, parid, nstat, &
168 iswind, &
169 istemp, ismois, ispstr, rio, rjo, rko, varobs, &
170 errf, i_parent_start, j_parent_start, &
171 ktau, iratio, npfi, iprt, &
172 ids,ide, jds,jde, kds,kde, &
173 ims,ime, jms,jme, kms,kme, &
174 its,ite, jts,jte, kts,kte )
175
176 !-----------------------------------------------------------------------
177 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
178 USE module_dm, ONLY : get_full_obs_vector
179 #endif
180
181 !-----------------------------------------------------------------------
182 IMPLICIT NONE
183 !-----------------------------------------------------------------------
184 !
185 ! PURPOSE: THIS SUBROUTINE CALCULATES THE DIFFERENCE BETWEEN THE
186 ! OBSERVED VALUES AND THE FORECAST VALUES AT THE OBSERVATION
187 ! POINTS. THE OBSERVED VALUES CLOSEST TO THE CURRENT
188 ! FORECAST TIME (XTIME) WERE DETERMINED IN SUBROUTINE
189 ! IN4DOB AND STORED IN ARRAY VAROBS. THE DIFFERENCES
190 ! CALCULATED BY SUBROUTINE ERROB WILL BE STORED IN ARRAY
191 ! ERRF FOR THE NSTA OBSERVATION LOCATIONS. MISSING
192 ! OBSERVATIONS ARE DENOTED BY THE DUMMY VALUE 99999.
193 !
194 ! HISTORY: Original author: MM5 version???
195 ! 02/04/2004 - Creation of WRF version. Al Bourgeois
196 ! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
197 !------------------------------------------------------------------------------
198
199 ! THE STORAGE ORDER IN VAROBS AND ERRF IS AS FOLLOWS:
200 ! IVAR VARIABLE TYPE(TAU-1)
201 ! ---- --------------------
202 ! 1 U error
203 ! 2 V error
204 ! 3 T error
205 ! 4 Q error
206 ! 5 Surface press error at T points (not used)
207 ! 6 Model surface press at T-points
208 ! 7 Model surface press at U-points
209 ! 8 Model surface press at V-points
210 ! 9 RKO at U-points
211
212 !-----------------------------------------------------------------------
213 !
214 ! Description of input arguments.
215 !
216 !-----------------------------------------------------------------------
217
218 INTEGER, INTENT(IN) :: inest ! Domain index.
219 INTEGER, INTENT(IN) :: nndgv ! Number of nudge variables.
220 INTEGER, INTENT(IN) :: nerrf ! Number of error fields.
221 INTEGER, INTENT(IN) :: niobf ! Number of observations.
222 INTEGER, INTENT(IN) :: maxdom ! Maximum number of domains.
223 INTEGER, INTENT(IN) :: levidn(maxdom) ! Level of nest.
224 INTEGER, INTENT(IN) :: parid(maxdom) ! Id of parent grid.
225 INTEGER, INTENT(IN) :: i_parent_start(maxdom) ! Start i index in parent domain.
226 INTEGER, INTENT(IN) :: j_parent_start(maxdom) ! Start j index in parent domain.
227 INTEGER, INTENT(IN) :: ktau
228 INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
229 INTEGER, INTENT(IN) :: npfi ! Coarse-grid diagnostics freq.
230 LOGICAL, INTENT(IN) :: iprt ! Print flag
231 INTEGER, INTENT(IN) :: nstat
232 INTEGER, intent(in) :: iswind
233 INTEGER, intent(in) :: istemp
234 INTEGER, intent(in) :: ismois
235 INTEGER, intent(in) :: ispstr
236 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
237 INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
238 INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
239
240 REAL, INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
241 REAL, INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
242 REAL, INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
243 REAL, INTENT(IN) :: t0
244 REAL, INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
245 REAL, INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme )
246 REAL, INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme ) ! Press. perturbation (Pa)
247 REAL, INTENT(IN) :: rovcp
248 REAL, INTENT(IN) :: uratx( ims:ime, jms:jme ) ! U to U10 ratio on mass points.
249 REAL, INTENT(IN) :: vratx( ims:ime, jms:jme ) ! V to V10 ratio on mass points.
250 REAL, INTENT(IN) :: tratx( ims:ime, jms:jme ) ! T to TH2 ratio on mass points.
251 REAL, INTENT(IN) :: rio(niobf) ! West-east coordinate.
252 REAL, INTENT(IN) :: rjo(niobf) ! South-north coordinate.
253 REAL, INTENT(INOUT) :: rko(niobf)
254 REAL, INTENT(INOUT) :: varobs(nndgv, niobf)
255 REAL, INTENT(INOUT) :: errf(nerrf, niobf)
256
257 ! Local variables
258 INTEGER :: iobmg(niobf) ! Obs i-coord on mass grid
259 INTEGER :: jobmg(niobf) ! Obs j-coord on mass grid
260 INTEGER :: ia(niobf)
261 INTEGER :: ib(niobf)
262 INTEGER :: ic(niobf)
263 REAL :: pbbo(kds:kde) ! column base pressure (cb) at obs loc.
264 REAL :: ppbo(kds:kde) ! column pressure perturbation (cb) at obs loc.
265
266 REAL :: ra(niobf)
267 REAL :: rb(niobf)
268 REAL :: rc(niobf)
269 REAL :: dxobmg(niobf) ! Interp. fraction (x dir) referenced to mass-grid
270 REAL :: dyobmg(niobf) ! Interp. fraction (y dir) referenced to mass-grid
271 INTEGER MM(MAXDOM)
272 INTEGER NNL
273 real :: uratio( ims:ime, jms:jme ) ! U to U10 ratio on momentum points.
274 real :: vratio( ims:ime, jms:jme ) ! V to V10 ratio on momentum points.
275 real :: pug1,pug2,pvg1,pvg2
276
277 ! Define staggers for U, V, and T grids, referenced from non-staggered grid.
278 real, parameter :: gridx_t = 0.5 ! Mass-point x stagger
279 real, parameter :: gridy_t = 0.5 ! Mass-point y stagger
280 real, parameter :: gridx_u = 0.0 ! U-point x stagger
281 real, parameter :: gridy_u = 0.5 ! U-point y stagger
282 real, parameter :: gridx_v = 0.5 ! V-point x stagger
283 real, parameter :: gridy_v = 0.0 ! V-point y stagger
284
285 real :: dummy = 99999.
286
287 real :: pbhi, pphi
288 real :: press,ttemp !ajb scratch variables
289 ! real model_temp,pot_temp !ajb scratch variables
290
291 !*** DECLARATIONS FOR IMPLICIT NONE
292 integer nsta,ivar,n,ityp
293 integer iob,job,kob,iob_ms,job_ms
294 integer k,kbot,nml,nlb,nle
295 integer iobm,jobm,iobp,jobp,kobp,inpf,i,j
296 integer i_start,i_end,j_start,j_end ! loop ranges for uratio,vratio calc.
297 integer k_start,k_end
298
299 real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
300 real pob
301 real grfacx,grfacy,uratiob,vratiob,tratiob,tratxob,fnpf
302 real stagx ! For x correction to mass-point stagger
303 real stagy ! For y correction to mass-point stagger
304
305 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
306 LOGICAL MP_LOCAL_DUMMASK(NIOBF) ! Mask for work to be done on this processor
307 LOGICAL MP_LOCAL_UOBMASK(NIOBF) ! Dot-point mask
308 LOGICAL MP_LOCAL_VOBMASK(NIOBF) ! Dot-point mask
309 LOGICAL MP_LOCAL_COBMASK(NIOBF) ! Cross-point mask
310 #endif
311 ! LOGICAL, EXTERNAL :: TILE_MASK
312
313 NSTA=NSTAT
314
315 ! FIRST, DETERMINE THE GRID TYPE CORRECTION FOR U-momentum, V-momentum,
316 ! AND MASS POINTS, AND WHEN INEST=2, CONVERT THE STORED COARSE MESH INDICES
317 ! TO THE FINE MESH INDEX EQUIVALENTS
318
319 ! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
320
321 if (iprt) then
322 write(6,'(a,i5,a,i2,a,i5,a)') '++++++CALL ERROB AT KTAU = ', &
323 KTAU,' AND INEST = ',INEST,': NSTA = ',NSTA,' ++++++'
324 endif
325
326 ERRF = 0.0 ! Zero out errf array
327
328 ! Set up loop bounds for this grid's boundary conditions
329 i_start = max( its-1,ids )
330 i_end = min( ite+1,ide-1 )
331 j_start = max( jts-1,jds )
332 j_end = min( jte+1,jde-1 )
333 k_start = kts
334 k_end = min( kte, kde-1 )
335
336 DO ITYP=1,3 ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP
337
338 ! Set grid stagger
339 IF(ITYP.EQ.1) THEN ! U-POINT CASE
340 GRIDX = gridx_u
341 GRIDY = gridy_u
342 ELSE IF(ITYP.EQ.2) THEN ! V-POINT CASE
343 GRIDX = gridx_v
344 GRIDY = gridy_v
345 ELSE ! MASS-POINT CASE
346 GRIDX = gridx_t
347 GRIDY = gridy_t
348 ENDIF
349
350 ! Compute URATIO and VRATIO fields on momentum (u,v) points.
351 IF(ityp.eq.1)THEN
352 call upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, uratx, uratio)
353 ELSE IF (ityp.eq.2) THEN
354 call vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, vratx, vratio)
355 ENDIF
356
357 IF(INEST.EQ.1) THEN ! COARSE MESH CASE...
358 DO N=1,NSTA
359 RA(N)=RIO(N)-GRIDX
360 RB(N)=RJO(N)-GRIDY
361 IA(N)=RA(N)
362 IB(N)=RB(N)
363 IOB=MAX0(1,IA(N))
364 IOB=MIN0(IOB,ide-1)
365 JOB=MAX0(1,IB(N))
366 JOB=MIN0(JOB,jde-1)
367 DXOB=RA(N)-FLOAT(IA(N))
368 DYOB=RB(N)-FLOAT(IB(N))
369
370 ! Save mass-point arrays for computing rko for all var types
371 if(ityp.eq.1) then
372 iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
373 jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
374 dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
375 dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
376 endif
377 iob_ms = iobmg(n)
378 job_ms = jobmg(n)
379 dxob_ms = dxobmg(n)
380 dyob_ms = dyobmg(n)
381
382
383 !if(n.eq.1 .and. iprt) then
384 ! write(6,*) 'ERROB - COARSE MESH:'
385 ! write(6,'(a,i1,a,i1,4(a,f5.2),2(a,i3),2(a,f6.3))') 'OBS= ',n, &
386 ! ' ityp= ',ityp, &
387 ! ' ra= ',ra(n),' rb= ',rb(n), &
388 ! ' rio= ',rio(n),' rjo= ',rjo(n), &
389 ! ' iob= ',iob,' job= ',job, &
390 ! ' dxob= ',dxob,' dyob= ',dyob
391 ! write(6,'(a,i3,a,i3,a,f5.2,a,f5.2)') &
392 ! ' iob_ms= ',iob_ms,' job_ms= ',job_ms, &
393 ! ' dxob_ms= ',dxob_ms,' dyob_ms= ',dyob_ms
394 !endif
395
396 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
397 ! Set mask for obs to be handled by this processor
398 MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
399
400 IF ( MP_LOCAL_DUMMASK(N) ) THEN
401 #endif
402
403 ! Interpolate pressure to obs location column and convert from Pa to cb.
404
405 do k = kds, kde
406 pbbo(k) = .001*( &
407 (1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) + &
408 DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
409 DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
410 DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
411 ppbo(k) = .001*( &
412 (1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) + &
413 DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
414 DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
415 DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
416
417 ! write(6,'(a,i2,2(a,f9.3)') ' k= ',k,' pbbo= ',pbbo(k),' ppbo= ',ppbo(k)
418 enddo
419
420 !ajb 20040119: Note based on bugfix for dot/cross points split across processors,
421 !ajb which was actually a serial code fix: The ityp=2 (v-points) and
422 !ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
423 !ajb case rko! This is necessary for bit-for-bit reproducability
424 !ajb with the parallel run. (coarse mesh)
425
426
427 if(abs(rko(n)+99).lt.1.)then
428 pob = varobs(5,n)
429
430 if(pob .gt.-800000.)then
431 do k=k_end-1,1,-1
432 kbot = k
433 if(pob .le. pbbo(k)+ppbo(k)) then
434 goto 199
435 endif
436 enddo
437 199 continue
438
439 pphi = ppbo(kbot+1)
440 pbhi = pbbo(kbot+1)
441
442 rko(n) = real(kbot+1)- &
443 ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
444
445 rko(n)=max(rko(n),1.0)
446 endif
447 endif
448
449 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
450 ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
451 #endif
452
453 RC(N)=RKO(N)
454
455 ENDDO ! END COARSE MESH LOOP OVER NSTA
456
457 ELSE ! FINE MESH CASE
458
459 ! CONVERT (I,J,K) OF OBSERVATIONS TO THE EQUIVALENT FINE MESH VALUES.
460 DO N=1,NSTA
461
462 ! COMPUTE THE OBS LOCATION WITH RESPECT TO THIS MESH (INEST)...
463 NML=INEST
464 MM(LEVIDN(INEST)+1)=INEST
465 ! WORKING TOWARD COARSER MESHES, DETERMINE THE HIERARCHY OF MOTHER
466 ! MESHES WITH RESPECT TO EACH MOTHER MESH STARTING AT MESH "IN"...
467 ! THAT IS, DETERMINE ITS MOTHER, GRANDMOTHER, GREAT-GRANDMOTHER, ETC.
468 ! OUT TO THE COARSE GRID MESH (INEST=1).
469 ! LEVIDN HOLDS THE NEST LEVEL AND PARID HOLDS THE MOTHER MESH FOR EACH
470 ! GRID (E.G., FOR 3 MESHES AND 2 NEST LEVELS, IN=1 IS THE COARSE GRID
471 ! MESH, IN=2 HAS LEVIDN(2)=1 AND PARID(2)=1, AND IN=3 HAS LEVIDN(3)=2
472 ! AND PARID(3)=2...)
473 DO NNL=LEVIDN(INEST),1,-1
474 MM(NNL)=PARID(NML)
475 NML=MM(NNL)
476 ENDDO
477
478 ! NOTE: MM(1) MUST BE THE COARSE GRID MESH (INEST=0)
479 IF(MM(1).NE.1) then
480 if(iprt) write(6,*) 'stopping in errob: inest = ',inest
481 STOP 21
482 ENDIF
483
484 RA(N)=RIO(N)
485 RB(N)=RJO(N)
486 DO NNL=1,LEVIDN(INEST)
487 GRFACX=0.
488 GRFACY=0.
489 ! COMPUTE THE OBS LOCATION WITH RESPECT TO THE INNER GRID IN NON-
490 ! STAGGERED SPACE (GRID=0.). WHEN WE REACH MESH INEST, THEN
491 ! APPLY THE APPRPRIATE STAGGER, DEPENDING ON THE VARIABLE...
492 IF(NNL.EQ.LEVIDN(INEST)) THEN
493 GRFACX=GRIDX
494 GRFACY=GRIDY
495 ENDIF
496
497 RA(N)=(RA(N)-FLOAT(i_parent_start(MM(NNL+1))))* &
498 FLOAT(IRATIO)+1.-GRFACX
499 RB(N)=(RB(N)-FLOAT(j_parent_start(MM(NNL+1))))* &
500 FLOAT(IRATIO)+1.-GRFACY
501
502 IA(N)=RA(N)
503 IB(N)=RB(N)
504 IOB=MAX0(1,IA(N))
505 IOB=MIN0(IOB,ide-1)
506 JOB=MAX0(1,IB(N))
507 JOB=MIN0(JOB,jde-1)
508 DXOB=RA(N)-FLOAT(IA(N))
509 DYOB=RB(N)-FLOAT(IB(N))
510
511 ! Save mass-point arrays for computing rko for all var types
512 if(ityp.eq.1) then
513 stagx = grfacx - gridx_t !Correct x stagger to mass-point
514 stagy = grfacy - gridy_t !Correct y stagger to mass-point
515 iobmg(n) = MIN0(MAX0(1,int(RA(n)+stagx)),ide-1)
516 jobmg(n) = MIN0(MAX0(1,int(RB(n)+stagy)),jde-1)
517 dxobmg(n) = RA(N)+stagx-FLOAT(int(RA(N)+stagx))
518 dyobmg(n) = RB(N)+stagy-FLOAT(int(RB(N)+stagy))
519 endif
520 iob_ms = iobmg(n)
521 job_ms = jobmg(n)
522 dxob_ms = dxobmg(n)
523 dyob_ms = dyobmg(n)
524
525 !if(n.eq.1) then
526 ! write(6,*) 'ERROB - FINE MESH:'
527 ! write(6,*) 'RA = ',ra(n),' RB = ',rb(n)
528 ! write(6,'(a,i1,a,i1,4(a,f5.2),2(a,i3),2(a,f6.3))') 'OBS= ',n, &
529 ! ' ityp= ',ityp, &
530 ! ' ra= ',ra(n),' rb= ',rb(n), &
531 ! ' rio= ',rio(n),' rjo= ',rjo(n), &
532 ! ' iob= ',iob,' job= ',job, &
533 ! ' dxob= ',dxob,' dyob= ',dyob
534 ! write(6,'(a,i3,a,i3,a,f5.2,a,f5.2)') &
535 ! ' iob_ms= ',iob_ms,' job_ms= ',job_ms, &
536 ! ' dxob_ms= ',dxob_ms,' dyob_ms= ',dyob_ms
537 !endif
538
539 ENDDO ! end do nnl=1,levidn(inest)
540
541 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
542 ! Set mask for obs to be handled by this processor
543 MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
544
545 IF ( MP_LOCAL_DUMMASK(N) ) THEN
546 #endif
547
548 ! Interpolate pressure to obs location column and convert from Pa to cb.
549
550 do k = kds, kde
551 pbbo(k) = .001*( &
552 (1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) + &
553 DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
554 DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
555 DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
556 ppbo(k) = .001*( &
557 (1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) + &
558 DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
559 DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
560 DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
561
562 ! write(6,'(a,i2,2(a,f9.3)') ' k= ',k,' pbbo= ',pbbo(k),' ppbo= ',ppbo(k)
563 enddo
564
565 !ajb 20040119: Note based on bugfix for dot/cross points split across processors,
566 !ajb which was actually a serial code fix: The ityp=2 (v-points) and
567 !ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
568 !ajb case) rko! This is necessary for bit-for-bit reproducability
569 !ajb with parallel run. (fine mesh)
570
571 if(abs(rko(n)+99).lt.1.)then
572 pob = varobs(5,n)
573
574 if(pob .gt.-800000.)then
575 do k=k_end-1,1,-1
576 kbot = k
577 if(pob .le. pbbo(k)+ppbo(k)) then
578 goto 198
579 endif
580 enddo
581 198 continue
582
583 pphi = ppbo(kbot+1)
584 pbhi = pbbo(kbot+1)
585
586 rko(n) = real(kbot+1)- &
587 ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
588 rko(n)=max(rko(n),1.0)
589 endif
590 endif
591
592 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
593 ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
594 #endif
595
596 RC(N)=RKO(N)
597
598 ENDDO ! END FINE MESH LOOP OVER NSTA
599
600 ENDIF ! end if(inest.eq.1)
601
602 ! INITIALIZE THE ARRAY OF DIFFERENCES BETWEEN THE OBSERVATIONS
603 ! AND THE LOCAL FORECAST VALUES TO ZERO. FOR THE FINE MESH
604 ! ONLY, SET THE DIFFERENCE TO A LARGE DUMMY VALUE IF THE
605 ! OBSERVATION IS OUTSIDE THE FINE MESH DOMAIN.
606
607 ! SET DIFFERENCE VALUE TO A DUMMY VALUE FOR OBS POINTS OUTSIDE
608 ! CURRENT DOMAIN
609 IF(ITYP.EQ.1) THEN
610 NLB=1
611 NLE=1
612 ELSE IF(ITYP.EQ.2) THEN
613 NLB=2
614 NLE=2
615 ELSE
616 NLB=3
617 NLE=5
618 ENDIF
619 DO IVAR=NLB,NLE
620 DO N=1,NSTA
621 IF((RA(N)-1.).LT.0)THEN
622 ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
623 ENDIF
624 IF((RB(N)-1.).LT.0)THEN
625 ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
626 ENDIF
627 IF((FLOAT(ide)-2.0*GRIDX-RA(N)-1.E-10).LT.0)THEN
628 ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
629 ENDIF
630 IF((FLOAT(jde)-2.0*GRIDY-RB(N)-1.E-10).LT.0)THEN
631 ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
632 ENDIF
633 if(rc(n).lt.1.)errf(ivar,n)=errf(ivar,n)+dummy
634 ENDDO
635 ENDDO
636
637 ! NOW FIND THE EXACT OFFSET OF EACH OBSERVATION FROM THE
638 ! GRID POINT TOWARD THE LOWER LEFT
639 DO N=1,NSTA
640 IA(N)=RA(N)
641 IB(N)=RB(N)
642 IC(N)=RC(N)
643 ENDDO
644 DO N=1,NSTA
645 RA(N)=RA(N)-FLOAT(IA(N))
646 RB(N)=RB(N)-FLOAT(IB(N))
647 RC(N)=RC(N)-FLOAT(IC(N))
648 ENDDO
649 ! PERFORM A TRILINEAR EIGHT-POINT (3-D) INTERPOLATION
650 ! TO FIND THE FORECAST VALUE AT THE EXACT OBSERVATION
651 ! POINTS FOR U, V, T, AND Q.
652
653 ! Compute local masks for dot and cross points.
654 if(ityp.eq.1) then
655 DO N=1,NSTA
656 IOB=MAX0(1,IA(N))
657 IOB=MIN0(IOB,ide-1)
658 JOB=MAX0(1,IB(N))
659 JOB=MIN0(JOB,jde-1)
660 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
661 ! Set mask for U-momemtum points to be handled by this processor
662 MP_LOCAL_UOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
663 #endif
664 ENDDO
665 endif
666 if(ityp.eq.2) then
667 DO N=1,NSTA
668 IOB=MAX0(1,IA(N))
669 IOB=MIN0(IOB,ide-1)
670 JOB=MAX0(1,IB(N))
671 JOB=MIN0(JOB,jde-1)
672 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
673 ! Set mask for V-momentum points to be handled by this processor
674 MP_LOCAL_VOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
675 #endif
676 ENDDO
677 endif
678 if(ityp.eq.3) then
679 DO N=1,NSTA
680 IOB=MAX0(1,IA(N))
681 IOB=MIN0(IOB,ide-1)
682 JOB=MAX0(1,IB(N))
683 JOB=MIN0(JOB,jde-1)
684 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
685 ! Set mask for cross (mass) points to be handled by this processor
686 MP_LOCAL_COBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
687 #endif
688 ENDDO
689 endif
690
691 !**********************************************************
692 ! PROCESS U VARIABLE (IVAR=1)
693 !**********************************************************
694 IF(ITYP.EQ.1) THEN
695 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
696 DO N=1,NSTA
697 IF(MP_LOCAL_UOBMASK(N)) THEN
698 ERRF(9,N)=rko(n) !RKO is needed by neighboring processors !ajb
699 ENDIF
700 ENDDO
701 #endif
702 IF(ISWIND.EQ.1) THEN
703 DO N=1,NSTA
704 IOB=MAX0(2,IA(N))
705 IOB=MIN0(IOB,ide-1)
706 IOBM=MAX0(1,IOB-1)
707 IOBP=MIN0(ide-1,IOB+1)
708 JOB=MAX0(2,IB(N))
709 JOB=MIN0(JOB,jde-1)
710 JOBM=MAX0(1,JOB-1)
711 JOBP=MIN0(jde-1,JOB+1)
712 KOB=MIN0(K_END,IC(N))
713
714 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
715 IF(MP_LOCAL_UOBMASK(N))THEN ! Do if obs on this processor
716 #endif
717 KOBP=MIN0(KOB+1,k_end)
718 DXOB=RA(N)
719 DYOB=RB(N)
720 DZOB=RC(N)
721
722 ! Compute surface pressure values at surrounding U and V points
723 PUG1 = .5*( pbase(IOBM,1,JOB) + pbase(IOB,1,JOB) )
724 PUG2 = .5*( pbase(IOB,1,JOB) + pbase(IOBP,1,JOB) )
725
726 ! This is to correct obs to model sigma level using reverse similarity theory
727 if(rko(n).eq.1.0)then
728 uratiob=((1.-DXOB)*((1.-DYOB)*uratio(IOB,JOB)+ &
729 DYOB*uratio(IOBP,JOB) &
730 )+DXOB*((1.-DYOB)*uratio(IOB,JOBP)+ &
731 DYOB*uratio(IOBP,JOBP)))
732 else
733 uratiob=1.
734 endif
735 !YLIU Some PBL scheme do not define the vratio/uratio
736 if(abs(uratiob).lt.1.0e-3) then
737 uratiob=1.
738 endif
739
740 ! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
741 ! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
742 ! OUTSIDE THE CURRENT DOMAIN
743
744 ! U COMPONENT WIND ERROR
745 ERRF(1,N)=ERRF(1,N)+uratiob*VAROBS(1,N)-((1.-DZOB)* &
746 ((1.-DyOB)*((1.- &
747 DxOB)*UB(IOB,KOB,JOB)+DxOB*UB(IOB+1,KOB,JOB) &
748 )+DyOB*((1.-DxOB)*UB(IOB,KOB,JOB+1)+DxOB* &
749 UB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
750 *UB(IOB,KOBP,JOB)+DxOB*UB(IOB+1,KOBP,JOB))+ &
751 DyOB*((1.-DxOB)*UB(IOB,KOBP,JOB+1)+DxOB* &
752 UB(IOB+1,KOBP,JOB+1))))
753
754 ! if(n.le.10) then
755 ! write(6,*)
756 ! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF1 at ',iob,job,kob, &
757 ! ' N = ',n,' inest = ',inest
758 ! write(6,*) 'VAROBS(1,N) = ',varobs(1,n)
759 ! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
760 ! write(6,*) 'UB(IOB,KOB,JOB) = ',UB(IOB,KOB,JOB)
761 ! write(6,*) 'UB(IOB+1,KOB,JOB) = ',UB(IOB+1,KOB,JOB)
762 ! write(6,*) 'UB(IOB,KOB,JOB+1) = ',UB(IOB,KOB,JOB+1)
763 ! write(6,*) 'UB(IOB+1,KOB,JOB+1) = ',UB(IOB+1,KOB,JOB+1)
764 ! write(6,*) 'UB(IOB,KOBP,JOB) = ',UB(IOB,KOBP,JOB)
765 ! write(6,*) 'UB(IOB+1,KOBP,JOB) = ',UB(IOB+1,KOBP,JOB)
766 ! write(6,*) 'UB(IOB,KOBP,JOB+1) = ',UB(IOB,KOBP,JOB+1)
767 ! write(6,*) 'UB(IOB+1,KOBP,JOB+1) = ',UB(IOB+1,KOBP,JOB+1)
768 ! write(6,*) 'uratiob = ',uratiob
769 ! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
770 ! write(6,*) 'ERRF(1,N) = ',errf(1,n)
771 ! write(6,*)
772 ! endif
773
774
775 ! Store model surface pressure (not the error!) at U point.
776 ERRF(7,N)=.001*( (1.-DXOB)*PUG1 + DXOB*PUG2 )
777
778 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
779 ENDIF ! end IF( MP_LOCAL_UOBMASK(N) )
780 #endif
781 ENDDO ! END U-POINT LOOP OVER OBS
782
783 ENDIF ! end if(iswind.eq.1)
784
785 ENDIF ! ITYP=1: PROCESS U
786
787 !**********************************************************
788 ! PROCESS V VARIABLE (IVAR=2)
789 !**********************************************************
790 IF(ITYP.EQ.2) THEN
791
792 IF(ISWIND.EQ.1) THEN
793 DO N=1,NSTA
794 IOB=MAX0(2,IA(N))
795 IOB=MIN0(IOB,ide-1)
796 IOBM=MAX0(1,IOB-1)
797 IOBP=MIN0(ide-1,IOB+1)
798 JOB=MAX0(2,IB(N))
799 JOB=MIN0(JOB,jde-1)
800 JOBM=MAX0(1,JOB-1)
801 JOBP=MIN0(jde-1,JOB+1)
802 KOB=MIN0(K_END,IC(N))
803
804 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
805 IF(MP_LOCAL_VOBMASK(N))THEN ! Do if obs on this processor
806 #endif
807 KOBP=MIN0(KOB+1,k_end)
808 DXOB=RA(N)
809 DYOB=RB(N)
810 DZOB=RC(N)
811
812 ! Compute surface pressure values at surrounding U and V points
813 PVG1 = .5*( pbase(IOB,1,JOBM) + pbase(IOB,1,JOB) )
814 PVG2 = .5*( pbase(IOB,1,JOB) + pbase(IOB,1,JOBP) )
815
816 ! This is to correct obs to model sigma level using reverse similarity theory
817 if(rko(n).eq.1.0)then
818 vratiob=((1.-DXOB)*((1.-DYOB)*vratio(IOB,JOB)+ &
819 DYOB*vratio(IOBP,JOB) &
820 )+DXOB*((1.-DYOB)*vratio(IOB,JOBP)+ &
821 DYOB*vratio(IOBP,JOBP)))
822 else
823 vratiob=1.
824 endif
825 !YLIU Some PBL scheme do not define the vratio/uratio
826 if(abs(vratiob).lt.1.0e-3) then
827 vratiob=1.
828 endif
829
830 ! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
831 ! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
832 ! OUTSIDE THE CURRENT DOMAIN
833
834 ! V COMPONENT WIND ERROR
835 ERRF(2,N)=ERRF(2,N)+vratiob*VAROBS(2,N)-((1.-DZOB)* &
836 ((1.-DyOB)*((1.- &
837 DxOB)*VB(IOB,KOB,JOB)+DxOB*VB(IOB+1,KOB,JOB) &
838 )+DyOB*((1.-DxOB)*VB(IOB,KOB,JOB+1)+DxOB* &
839 VB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
840 *VB(IOB,KOBP,JOB)+DxOB*VB(IOB+1,KOBP,JOB))+ &
841 DyOB*((1.-DxOB)*VB(IOB,KOBP,JOB+1)+DxOB* &
842 VB(IOB+1,KOBP,JOB+1))))
843
844 ! Store model surface pressure (not the error!) at V point.
845 ERRF(8,N)=.001*( (1.-DYOB)*PVG1 + DYOB*PVG2 )
846
847 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
848 ENDIF ! end IF( MP_LOCAL_VOBMASK(N) )
849 #endif
850 ENDDO ! END V-POINT LOOP OVER OBS
851
852 ENDIF ! end if(iswind.eq.1)
853
854 ENDIF ! ITYP=1: PROCESS V
855
856 !**********************************************************
857 ! PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
858 !**********************************************************
859 IF(ITYP.EQ.3) THEN
860
861 IF(ISTEMP.EQ.1 .OR. ISMOIS.EQ.1) THEN
862 DO N=1,NSTA
863 IOB=MAX0(1,IA(N))
864 IOB=MIN0(IOB,ide-1)
865 JOB=MAX0(1,IB(N))
866 JOB=MIN0(JOB,jde-1)
867 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
868 IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
869 #endif
870 KOB=MIN0(k_end,IC(N))
871 KOBP=MIN0(KOB+1,K_END)
872 DXOB=RA(N)
873 DYOB=RB(N)
874 DZOB=RC(N)
875
876 ! This is to correct obs to model sigma level using reverse similarity theory
877 if(rko(n).eq.1.0)then
878 tratxob=((1.-DXOB)*((1.-DYOB)*tratx(IOB,JOB)+ &
879 DYOB*tratx(IOB+1,JOB) &
880 )+DXOB*((1.-DYOB)*tratx(IOB,JOB+1)+ &
881 DYOB*tratx(IOB+1,JOB+1)))
882 else
883 tratxob=1.
884 endif
885
886 !yliu
887 if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
888
889 !ajb testing only
890 if(iprt .and. n.eq.81) then
891 write(6,*) 'POTENTIAL TEMP FOR N=81:'
892 write(6,*)
893 write(6,*) ' K THETA TEMPERATURE', &
894 ' PBASE'
895 write(6,*)
896 do k=k_end,1,-1
897 press = pbase(iob,k,job)+pp(iob,k,job)
898 ttemp = exp ( alog(300.+TB(IOB,k,JOB)) - &
899 .2857143*alog(100000./press) )
900 write(6,*) k,300.+TB(IOB,k,JOB),ttemp,pbase(iob,k,job)
901 enddo
902 endif
903 !ajb end testing only
904
905 ! TEMPERATURE ERROR
906 ! if(n.le.10) then
907 ! write(6,*) 'before: errf(3,n) = ',errf(3,n)
908 ! endif
909 ERRF(3,N)=ERRF(3,N)+tratxob*VAROBS(3,N)-((1.-DZOB)* &
910 ((1.-DyOB)*((1.- &
911 DxOB)*(TB(IOB,KOB,JOB))+DxOB* &
912 (TB(IOB+1,KOB,JOB)))+DyOB*((1.-DxOB)* &
913 (TB(IOB,KOB,JOB+1))+DxOB* &
914 (TB(IOB+1,KOB,JOB+1))))+DZOB*((1.- &
915 DyOB)*((1.-DxOB)*(TB(IOB,KOBP,JOB))+DxOB* &
916 (TB(IOB+1,KOBP,JOB)))+DyOB*((1.-DxOB)* &
917 (TB(IOB,KOBP,JOB+1))+DxOB* &
918 (TB(IOB+1,KOBP,JOB+1)))))
919
920 ! if(n.le.10) then
921 ! write(6,*)
922 ! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF3 at ',iob,job,kob, &
923 ! ' N = ',n,' inest = ',inest
924 ! write(6,*) 'VAROBS(3,N) = ',varobs(3,n)
925 ! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
926 ! write(6,*) 'TB(IOB,KOB,JOB) = ',TB(IOB,KOB,JOB)
927 ! write(6,*) 'TB(IOB+1,KOB,JOB) = ',TB(IOB+1,KOB,JOB)
928 ! write(6,*) 'TB(IOB,KOB,JOB+1) = ',TB(IOB,KOB,JOB+1)
929 ! write(6,*) 'TB(IOB+1,KOB,JOB+1) = ',TB(IOB+1,KOB,JOB+1)
930 ! write(6,*) 'TB(IOB,KOBP,JOB) = ',TB(IOB,KOBP,JOB)
931 ! write(6,*) 'TB(IOB+1,KOBP,JOB) = ',TB(IOB+1,KOBP,JOB)
932 ! write(6,*) 'TB(IOB,KOBP,JOB+1) = ',TB(IOB,KOBP,JOB+1)
933 ! write(6,*) 'TB(IOB+1,KOBP,JOB+1) = ',TB(IOB+1,KOBP,JOB+1)
934 ! write(6,*) 'tratxob = ',tratxob
935 ! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
936 ! write(6,*) 'ERRF(3,N) = ',errf(3,n)
937 ! write(6,*)
938 ! endif
939
940
941 ! MOISTURE ERROR
942 ERRF(4,N)=ERRF(4,N)+VAROBS(4,N)-((1.-DZOB)*((1.-DyOB)*((1.- &
943 DxOB)*QVB(IOB,KOB,JOB)+DxOB* &
944 QVB(IOB+1,KOB,JOB))+DyOB*((1.-DxOB)* &
945 QVB(IOB,KOB,JOB+1)+DxOB* &
946 QVB(IOB+1,KOB,JOB+1)))+DZOB*((1.- &
947 DyOB)*((1.-DxOB)*QVB(IOB,KOBP,JOB)+DxOB &
948 *QVB(IOB+1,KOBP,JOB))+DyOB*((1.-DxOB &
949 )*QVB(IOB,KOBP,JOB+1)+DxOB* &
950 QVB(IOB+1,KOBP,JOB+1))))
951
952 ! Store model surface pressure (not the error!) at T-point
953 ERRF(6,N)= .001* &
954 ((1.-DyOB)*((1.-DxOB)*pbase(IOB,1,JOB)+DxOB* &
955 pbase(IOB+1,1,JOB))+DyOB*((1.-DxOB)* &
956 pbase(IOB,1,JOB+1)+DxOB*pbase(IOB+1,1,JOB+1) ))
957
958 if(iprt .and. n.eq.81) then
959 write(6,*) 'ERRF(6,81) calculation:'
960 write(6,*) 'iob,job = ',iob,job
961 write(6,*) 'pbase(iob,1,job) = ',pbase(iob,1,job)
962 write(6,*) 'pbase(iob+1,1,job) = ',pbase(iob+1,1,job)
963 write(6,*) 'pbase(iob,1,job+1) = ',pbase(iob,1,job+1)
964 write(6,*) 'pbase(iob+1,1,job+1) = ',pbase(iob+1,1,job+1)
965 write(6,*) 'ERRF(6,81) = ',errf(6,n)
966 ! call flush(6)
967 endif
968
969 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
970 ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
971 #endif
972 ENDDO ! END T and QV LOOP OVER OBS
973
974 ENDIF !end if(istemp.eq.1 .or. ismois.eq.1)
975
976 !**********************************************************
977 ! PROCESS SURFACE PRESSURE CROSS-POINT FIELD, IVAR=5,
978 ! USING BILINEAR FOUR-POINT 2-D INTERPOLATION
979 !**********************************************************
980 IF(ISPSTR.EQ.1) THEN
981 DO N=1,NSTA
982 IOB=MAX0(1,IA(N))
983 IOB=MIN0(IOB,ide-1)
984 JOB=MAX0(1,IB(N))
985 JOB=MIN0(JOB,jde-1)
986 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
987 IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
988 #endif
989 DXOB=RA(N)
990 DYOB=RB(N)
991 !ajb fix this (put in correct pressure calc for IOB,JOB here)
992 ERRF(5,N)=ERRF(5,N)+VAROBS(5,N)-((1.-DyOB)*((1.-DxOB)* &
993 pbase(IOB,1,JOB)+DxOB*pbase(IOB+1,1,JOB))+DyOB* &
994 ((1.-DxOB)*pbase(IOB,1,JOB+1)+DxOB* &
995 pbase(IOB+1,1,JOB+1)))
996
997 if(n.eq.81) then
998 write(6,*) 'ERRF(5,81) calculation:'
999 write(6,*) 'iob,job = ',iob,job
1000 write(6,*) 'pbase(iob,1,job) = ',pbase(iob,1,job)
1001 write(6,*) 'pbase(iob+1,1,job) = ',pbase(iob+1,1,job)
1002 write(6,*) 'pbase(iob,1,job+1) = ',pbase(iob,1,job+1)
1003 write(6,*) 'pbase(iob+1,1,job+1) = ',pbase(iob+1,1,job+1)
1004 write(6,*) 'errf(5,81) = ',errf(5,n)
1005 ! call flush(6)
1006 endif
1007
1008 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1009 ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
1010 #endif
1011
1012 ENDDO
1013
1014 ENDIF ! end if(ispstr.eq.1)
1015
1016 ENDIF ! end if(ityp.eq.3)
1017
1018 ENDDO ! END BIG LOOP
1019
1020 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1021 ! Gather the errf values calculated by the processors with the obs.
1022 CALL get_full_obs_vector(nsta, nerrf, niobf, mp_local_uobmask, &
1023 mp_local_vobmask, mp_local_cobmask, errf)
1024 #endif
1025
1026 ! DIFFERENCE BETWEEN OBS AND FCST IS COMPLETED
1027 IF(INEST.EQ.1)THEN
1028 INPF=NPFI
1029 ELSE
1030 FNPF=IRATIO**LEVIDN(INEST)
1031 INPF=FNPF*NPFI
1032 ENDIF
1033 ! Gross error check for temperature. Set all vars bad.
1034 do n=1,nsta
1035 if((abs(errf(3,n)).gt.20.).and. &
1036 (errf(3,n).gt.-800000.))then
1037
1038 errf(1,n)=-888888.
1039 errf(2,n)=-888888.
1040 errf(3,n)=-888888.
1041 errf(4,n)=-888888.
1042 varobs(1,n)=-888888.
1043 varobs(2,n)=-888888.
1044 varobs(3,n)=-888888.
1045 varobs(4,n)=-888888.
1046 endif
1047 enddo
1048
1049 ! For printout
1050 ! IF(MOD(KTAU,INPF).NE.0) THEN
1051 ! RETURN
1052 ! ENDIF
1053
1054 RETURN
1055 END SUBROUTINE errob
1056
1057 SUBROUTINE upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, &
1058 arrin, arrout)
1059 !------------------------------------------------------------------------------
1060 ! PURPOSE: This subroutine interpolates a real 2D array defined over mass
1061 ! coordinate points, to U (momentum) points.
1062 !
1063 !------------------------------------------------------------------------------
1064 IMPLICIT NONE
1065
1066 INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
1067 INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
1068 INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
1069 INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
1070 INTEGER, INTENT(IN) :: ids ! Starting i index for entire model domain
1071 INTEGER, INTENT(IN) :: ide ! Ending i index for entire model domain
1072 INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
1073 INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
1074 INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
1075 INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
1076 REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
1077 REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on U points
1078
1079 ! Local variables
1080 integer :: i, j
1081
1082 ! Do domain interior first
1083 do j = j_start, j_end
1084 do i = max(2,i_start), i_end
1085 arrout(i,j) = 0.5*(arrin(i,j)+arrin(i-1,j))
1086 enddo
1087 enddo
1088
1089 ! Do west-east boundaries
1090 if(i_start .eq. ids) then
1091 do j = j_start, j_end
1092 arrout(i_start,j) = arrin(i_start,j)
1093 enddo
1094 endif
1095 if(i_end .eq. ide-1) then
1096 do j = j_start, j_end
1097 arrout(i_end+1,j) = arrin(i_end,j)
1098 enddo
1099 endif
1100
1101 RETURN
1102 END SUBROUTINE upoint
1103
1104 SUBROUTINE vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, &
1105 arrin, arrout)
1106 !------------------------------------------------------------------------------
1107 ! PURPOSE: This subroutine interpolates a real 2D array defined over mass
1108 ! coordinate points, to V (momentum) points.
1109 !
1110 !------------------------------------------------------------------------------
1111 IMPLICIT NONE
1112
1113 INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
1114 INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
1115 INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
1116 INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
1117 INTEGER, INTENT(IN) :: jds ! Starting j index for entire model domain
1118 INTEGER, INTENT(IN) :: jde ! Ending j index for entire model domain
1119 INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
1120 INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
1121 INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
1122 INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
1123 REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
1124 REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on V points
1125
1126 ! Local variables
1127 integer :: i, j
1128
1129 ! Do domain interior first
1130 do j = max(2,j_start), j_end
1131 do i = i_start, i_end
1132 arrout(i,j) = 0.5*(arrin(i,j)+arrin(i,j-1))
1133 enddo
1134 enddo
1135
1136 ! Do south-north boundaries
1137 if(j_start .eq. jds) then
1138 do i = i_start, i_end
1139 arrout(i,j_start) = arrin(i,j_start)
1140 enddo
1141 endif
1142 if(j_end .eq. jde-1) then
1143 do i = i_start, i_end
1144 arrout(i,j_end+1) = arrin(i,j_end)
1145 enddo
1146 endif
1147
1148 RETURN
1149 END SUBROUTINE vpoint
1150
1151 LOGICAL FUNCTION TILE_MASK(iloc, jloc, its, ite, jts, jte)
1152 !------------------------------------------------------------------------------
1153 ! PURPOSE: Check to see if an i, j grid coordinate is in the tile index range.
1154 !
1155 ! Returns: TRUE if the grid coordinate (ILOC,JLOC) is in the tile defined by
1156 ! tile-range indices (its,jts) and (ite,jte)
1157 ! FALSE otherwise.
1158 !
1159 !------------------------------------------------------------------------------
1160 IMPLICIT NONE
1161
1162 INTEGER, INTENT(IN) :: iloc
1163 INTEGER, INTENT(IN) :: jloc
1164 INTEGER, INTENT(IN) :: its
1165 INTEGER, INTENT(IN) :: ite
1166 INTEGER, INTENT(IN) :: jts
1167 INTEGER, INTENT(IN) :: jte
1168
1169 ! Local variables
1170 LOGICAL :: retval
1171
1172 TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND. &
1173 jloc .LE. jte .AND. jloc .GE. jts )
1174
1175 RETURN
1176 END FUNCTION TILE_MASK
1177
1178 !-----------------------------------------------------------------------
1179 SUBROUTINE nudob(j, ivar, aten, inest, ifrest, ktau, ktaur, &
1180 xtime, mu, msf, nndgv, nerrf, niobf, maxdom, &
1181 npfi, ionf, rinxy, twindo, levidn, &
1182 parid, nstat, i_parent_start, j_parent_start, &
1183 fdob, lev_in_ob, plfo, nlevs_ob, &
1184 iratio, dx, dtmin, rio, rjo, rko, &
1185 timeob, varobs, errf, pbase, ptop, pp, &
1186 iswind, istemp, ismois, giv, git, giq, &
1187 savwt, kpblt, nscan, &
1188 iprt, &
1189 ids,ide, jds,jde, kds,kde, & ! domain dims
1190 ims,ime, jms,jme, kms,kme, & ! memory dims
1191 its,ite, jts,jte, kts,kte ) ! tile dims
1192
1193 !-----------------------------------------------------------------------
1194 USE module_model_constants
1195 USE module_domain
1196 !-----------------------------------------------------------------------
1197 IMPLICIT NONE
1198 !-----------------------------------------------------------------------
1199 !
1200 ! PURPOSE: THIS SUBROUTINE GENERATES NUDGING TENDENCIES FOR THE J-TH
1201 ! VERTICAL SLICE (I-K PLANE) FOR FOUR-DIMENSIONAL DATA
1202 ! ASSIMILATION FROM INDIVIDUAL OBSERVATIONS. THE NUDGING
1203 ! TENDENCIES ARE FOUND FROM A ONE-PASS CALCULATION OF
1204 ! WEIGHTING FACTORS SIMILAR TO THE BENJAMIN-SEAMAN OBJECTIVE
1205 ! ANALYSIS. THIS SUBROUTINE IS DESIGNED FOR RAPID EXECUTION
1206 ! AND MINIMAL STORAGE REQUIREMENTS. ALGORITHMS SHOULD BE
1207 ! VECTORIZED WHEREVER POSSIBLE.
1208 !
1209 ! HISTORY: Original author: MM5 version???
1210 ! 02/04/2004 - Creation of WRF version. Al Bourgeois
1211 ! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
1212 !------------------------------------------------------------------------------
1213 !
1214 ! NOTE: This routine was originally designed for MM5, which uses
1215 ! a nonstandard (I,J) coordinate system. For WRF, I is the
1216 ! east-west running coordinate, and J is the south-north
1217 ! running coordinate. So a "J-slab" here is west-east in
1218 ! extent, not south-north as for MM5. -ajb 06/10/2004
1219 !
1220 ! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
1221 ! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
1222 !
1223 ! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
1224 ! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
1225 ! TYPES OF FACTORS:
1226 ! 1) TIME WEIGHTING - ONLY OBSERVATIONS WITHIN A SELECTED
1227 ! TIME WINDOW (TWINDO) CENTERED AT THE CURRENT FORECAST
1228 ! TIME (XTIME) ARE USED. OBSERVATIONS CLOSEST TO
1229 ! XTIME ARE TIME-WEIGHTED MOST HEAVILY (TIMEWT)
1230 ! 2) VERTICAL WEIGHTING - NON-ZERO WEIGHTS (WTSIG) ARE
1231 ! CALCULATED WITHIN A VERTICAL REGION OF INFLUENCE
1232 ! (RINSIG).
1233 ! 3) HORIZONTAL WEIGHTING - NON-ZERO WEIGHTS (WTIJ) ARE
1234 ! CALCULATED WITHIN A RADIUS OF INFLUENCE (RINXY). THE
1235 ! VALUE OF RIN IS DEFINED IN KILOMETERS, AND CONVERTED
1236 ! TO GRID LENGTHS FOR THE APPROPRIATE MESH SIZE.
1237 !
1238 ! THE FIVE FORECAST VARIABLES ARE PROCESSED BY CHANGING THE
1239 ! VALUE OF IVAR AS FOLLOWS:
1240 ! IVAR VARIABLE(TAU-1)
1241 ! ---- ---------------
1242 ! 1 U
1243 ! 2 V
1244 ! 3 T
1245 ! 4 QV
1246 ! 5 PSB(CROSS) REMOVED IN V3
1247 ! (6) PSB(DOT)
1248 !
1249 !-----------------------------------------------------------------------
1250 !
1251 ! Description of input arguments.
1252 !
1253 !-----------------------------------------------------------------------
1254
1255 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
1256 INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
1257 INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
1258 INTEGER, INTENT(IN) :: j ! south-north running coordinate.
1259 INTEGER, INTENT(IN) :: ivar
1260 INTEGER, INTENT(IN) :: inest ! domain index
1261 LOGICAL, INTENT(IN) :: ifrest
1262 INTEGER, INTENT(IN) :: ktau
1263 INTEGER, INTENT(IN) :: ktaur
1264 REAL, INTENT(IN) :: xtime ! forecast time in minutes
1265 INTEGER, INTENT(IN) :: nndgv ! number of nudge variables
1266 INTEGER, INTENT(IN) :: nerrf ! number of error fields
1267 INTEGER, INTENT(IN) :: niobf ! number of observations
1268 INTEGER, INTENT(IN) :: maxdom ! maximum number of domains
1269 INTEGER, INTENT(IN) :: npfi
1270 INTEGER, INTENT(IN) :: ionf
1271 REAL, INTENT(IN) :: rinxy
1272 REAL, INTENT(IN) :: twindo
1273 INTEGER, INTENT(IN) :: levidn(maxdom) ! level of nest
1274 INTEGER, INTENT(IN) :: parid(maxdom) ! parent domain id
1275 INTEGER, INTENT(IN) :: nstat ! number of obs stations
1276 INTEGER, INTENT(IN) :: i_parent_start(maxdom) ! Start i index in parent domain.
1277 INTEGER, INTENT(IN) :: j_parent_start(maxdom) ! Start j index in parent domain.
1278 TYPE(fdob_type), intent(inout) :: fdob
1279 REAL, INTENT(IN) :: lev_in_ob(niobf) ! Level in sounding-type obs.
1280 REAL, intent(IN) :: plfo(niobf)
1281 REAL, INTENT(IN) :: nlevs_ob(niobf) ! Number of levels in sounding.
1282 INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
1283 REAL, INTENT(IN) :: dx ! This domain grid cell-size (m)
1284 REAL, INTENT(IN) :: dtmin
1285 REAL, INTENT(IN) :: rio(niobf)
1286 REAL, INTENT(IN) :: rjo(niobf)
1287 REAL, INTENT(INOUT) :: rko(niobf)
1288 REAL, INTENT(IN) :: timeob(niobf)
1289 REAL, INTENT(IN) :: varobs(nndgv,niobf)
1290 REAL, INTENT(IN) :: errf(nerrf, niobf)
1291 REAL, INTENT(IN) :: pbase( ims:ime, kms:kme ) ! Base pressure.
1292 REAL, INTENT(IN) :: ptop
1293 REAL, INTENT(IN) :: pp( ims:ime, kms:kme ) ! Pressure perturbation (Pa)
1294 REAL, INTENT(IN) :: mu(ims:ime) ! Air mass on u, v, or mass-grid
1295 REAL, INTENT(IN) :: msf(ims:ime) ! Map scale (only used for vars u & v)
1296 INTEGER, intent(in) :: iswind ! Nudge flag for wind
1297 INTEGER, intent(in) :: istemp ! Nudge flag for temperature
1298 INTEGER, intent(in) :: ismois ! Nudge flag for moisture
1299 REAL, intent(in) :: giv ! Coefficient for wind
1300 REAL, intent(in) :: git ! Coefficient for temperature
1301 REAL, intent(in) :: giq ! Coefficient for moisture
1302 REAL, INTENT(INOUT) :: aten( ims:ime, kms:kme)
1303 REAL, INTENT(INOUT) :: savwt( nndgv, ims:ime, kms:kme )
1304 INTEGER, INTENT(IN) :: kpblt(its:ite)
1305 INTEGER, INTENT(IN) :: nscan ! number of scans
1306 LOGICAL, INTENT(IN) :: iprt ! print flag
1307
1308 ! Local variables
1309 integer :: mm(maxdom)
1310 real :: ra(niobf)
1311 real :: rb(niobf)
1312 real :: psurf(niobf)
1313 real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
1314 real :: rscale(ims:ime) ! For converting to rho-coupled units.
1315 ! real :: tscale(ims:ime,kms:kme) ! For converting to potential temp. units.
1316 real :: reserf(100)
1317 character*40 name
1318 character*3 chr_hr
1319
1320 !*** DECLARATIONS FOR IMPLICIT NONE
1321 integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
1322 integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
1323 integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
1324 integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
1325 integer :: komin,komax,nn,nhi,nlo,nnjc
1326 integer :: i_s,i_e
1327 integer :: istq
1328 real :: gfactor,rfactor,gridx,gridy,rindx,schnes,ris
1329 real :: grfacx,grfacy
1330 real :: fdtim,tw1,tw2,tconst,timewt,timewt2,ttim,dift,pob
1331 real :: ri,rj,rx,ry,rsq,wtij,pdfac,erfivr,dk,slope,rinfac
1332 real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
1333
1334 real :: scratch
1335
1336 ! print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
1337 ! if(ivar.ne.4)return
1338 !yliu start -- for multi-scans: NSCAN=0: original
1339 ! NSCAN=1: added a scan with a larger Ri and smaller G
1340 ! if(NSCAN.ne.0 .and. NSCAN.ne.1) stop
1341 ! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
1342 if(NSCAN.ne.0) then
1343 IF (iprt) write(6,*) 'SAVWT must be resized for NSCAN=1'
1344 stop
1345 endif
1346 IPLO=0 + NSCAN*4
1347 GFACTOR=1. + NSCAN*(-1. + 0.33333)
1348 RFACTOR=1. + NSCAN*(-1. + 3.0)
1349 !yliu end
1350 ! jc
1351
1352 ! return if too close to j boundary
1353 if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
1354 ! write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
1355 ! $ ' too close to boundary.'
1356 return
1357 endif
1358 if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
1359 ! write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
1360 ! $ ' too close to boundary.'
1361 return
1362 endif
1363
1364 ! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
1365 ICUT=0
1366 IF(INEST.GT.1)ICUT=1
1367 i_s = max0(2+icut,its)
1368 i_e = min0(ide-1-icut,ite)
1369
1370 IPL=IVAR + IPLO !yliu +IPLO
1371
1372 ! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
1373
1374 IF(INEST.EQ.1)THEN
1375 INPF=NPFI
1376 INFR=IONF
1377 ELSE
1378 IF(IRATIO.NE.3) THEN
1379 IF (iprt) THEN
1380 write(6,*) 'iratio = ',iratio
1381 write(6,*) 'stop 1 in nudob: iratio = ',iratio
1382 ENDIF
1383 STOP 1
1384 ENDIF
1385 INPF=(3**LEVIDN(INEST))*NPFI
1386 INFR=(3**LEVIDN(INEST))*IONF
1387 ENDIF
1388 GRIDX=0.0
1389 GRIDY=0.0
1390 IGRID=0
1391 IF(IVAR.GE.3)THEN
1392 GRIDX=0.5
1393 GRIDY=0.5
1394 IGRID=1
1395 ELSEIF(IVAR.eq.1) THEN
1396 GRIDY=0.5
1397 IGRID=1
1398 ELSEIF(IVAR.eq.2) THEN
1399 GRIDX=0.5
1400 IGRID=1
1401 ENDIF
1402
1403 ! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
1404 ! KILOMETERS TO GRID LENGTHS, RINDX
1405
1406 RINDX=RINXY*1000./DX * RFACTOR !yliu *RFACTOR
1407
1408 ! jc
1409 ! make horizontal radius vary per nest
1410 ! rindx=rindx/float(inest)
1411 ! yliu test1 -- En 3, 4
1412 ! rindx=rindx/float(3**(in-1)) !YLIU
1413 ! jc
1414 ! make horizontal radius vary per nest
1415 ! schnes=1/float(inest) !JC
1416 ! yliu test1 -- En 3, 4 !YLIU
1417 schnes=1/float(3**(inest-1)) !JC
1418 ! reduce the Rinf in the nested grid proportionally
1419 rindx=rindx*schnes
1420 ! rinfmn =1., rinfmx=2., pfree=50 in param.F
1421 ! yliu test: for upper-air data, use larger influence radii
1422 ! Essentially increase the slope -- the same radii
1423 ! at 500mb and above as the coarse mesh and the
1424 ! same small radii at sfc as that for sfc obs
1425 fdob%rinfmx=2. *1.0 /schnes !YLIU
1426 ! rinfmx=1.2*1.0 /schnes !YLIU
1427 ! jc
1428 RIS=RINDX*RINDX
1429 IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
1430 IF(MOD(KTAU,INFR).NE.0)GOTO 126
1431 5 CONTINUE
1432 IF (iprt) THEN
1433 IF(J.EQ.10) write(6,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
1434 ENDIF
1435 6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,', &
1436 'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2, &
1437 ' rindx=',f4.1)
1438
1439 ! SET RA AND RB
1440 IF(INEST.EQ.1) THEN
1441
1442 ! SET RA AND RB FOR THE COARSE MESH...
1443 DO N=1,NSTAT
1444 RA(N)=RIO(N)-GRIDX
1445 RB(N)=RJO(N)-GRIDY
1446 ENDDO
1447
1448 ELSE
1449
1450 ! SET RA AND RB FOR THE FINE MESH CASE...
1451 DO N=1,NSTAT
1452
1453 ! COMPUTE THE OBS LOCATION WITH RESPECT TO THIS MESH (INEST)...
1454 NML=INEST
1455 MM(LEVIDN(INEST)+1)=INEST
1456 ! WORKING TOWARD COARSER MESHES, DETERMINE THE HIERARCHY OF MOTHER
1457 ! MESHES WITH RESPECT TO EACH MOTHER MESH STARTING AT MESH "INEST"...
1458 ! THAT IS, DETERMINE ITS MOTHER, GRANDMOTHER, GREAT-GRANDMOTHER, ETC.
1459 ! OUT TO THE COARSE GRID MESH (INEST=1).
1460 ! LEVIDN HOLDS THE NEST LEVEL AND PARID HOLDS THE MOTHER MESH FOR EACH
1461 ! GRID (E.G., FOR 3 MESHES AND 2 NEST LEVELS, INEST=1 IS THE COARSE GRID
1462 ! MESH, INEST=2 HAS LEVIDN(2)=1 AND PARID(2)=1, AND INEST=3 HAS LEVIDN(3)=2
1463 ! AND PARID(3)=2...)
1464 DO NNL=LEVIDN(INEST),1,-1
1465 MM(NNL)=PARID(NML)
1466 NML=MM(NNL)
1467 ENDDO
1468
1469 ! MM(1) MUST BE THE COARSE GRID MESH (INEST=0)
1470
1471 IF(MM(1).NE.1) then
1472 IF (iprt) write(6,*) 'stop 21 in nudob: inest = ',inest
1473 STOP 21
1474 ENDIF
1475 RA(N)=RIO(N)
1476 RB(N)=RJO(N)
1477 DO NNL=1,LEVIDN(INEST)
1478 GRFACX=0.
1479 GRFACY=0.
1480 ! COMPUTE THE OBS LOCATION WITH RESPECT TO THE INNER GRID IN DOT-POINT
1481 ! SPACE (GRID=0.). WHEN WE REACH MESH IN, THEN USE "GRID" TO GO TO
1482 ! CROSS OR DOT DEPENDING ON THE VARIABLE...
1483 IF(NNL.EQ.LEVIDN(INEST)) THEN
1484 GRFACX=GRIDX
1485 GRFACY=GRIDY
1486 ENDIF
1487
1488 RA(N)=(RA(N)-FLOAT(i_parent_start(MM(NNL+1))))* &
1489 FLOAT(IRATIO)+1.-GRFACX
1490 RB(N)=(RB(N)-FLOAT(j_parent_start(MM(NNL+1))))* &
1491 FLOAT(IRATIO)+1.-GRFACY
1492
1493 ENDDO
1494
1495 ENDDO ! END LOOP OVER OBS STATIONS FOR FINE MESH CASE
1496
1497 ENDIF ! END SECTION FOR SETTING RA AND RB
1498
1499
1500 ! OUTPUT OBS PER GRID EVERY HOUR
1501 if ( mod(xtime,60.).gt.56. .and. ivar.eq.3 .and. j.eq.10) then
1502 IF (iprt) print *,'outputting obs number on grid ', &
1503 inest,' at time=',xtime
1504 write(chr_hr(1:3),'(i3)')nint(xtime/60.)
1505 if(chr_hr(1:1).eq.' ')chr_hr(1:1)='0'
1506 if(chr_hr(2:2).eq.' ')chr_hr(2:2)='0'
1507 IF (iprt) print *,'chr_hr=',chr_hr(1:3),nint(xtime/60.)
1508 open(91,file= &
1509 'obs_g'//char(inest+ichar('0'))//'_'//chr_hr(1:3), &
1510 form='FORMATted',status='unknown')
1511 write(91,911)nstat
1512 write(6,911)nstat
1513 911 FORMAT('total obs=',i8)
1514 nsthis=0
1515 nsmetar=0
1516 nsspeci=0
1517 nsship=0
1518 nssynop=0
1519 nstemp=0
1520 nspilot=0
1521 nssatwnds=0
1522 nssams=0
1523 nsprofs=0
1524 ! print *,'ide,jde=',ide,jde
1525 do jjjn=1,nstat
1526 ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
1527 FDTIM=XTIME-DTMIN
1528 ! CONVERT TWINDO AND TIMEOB FROM HOURS TO MINUTES:
1529 TW1=TWINDO/2.*60.
1530 TW2=TWINDO*60.
1531 TCONST=1./TW1
1532 TIMEWT2=0.0
1533 TTIM=TIMEOB(jjjn)*60.
1534 !***********TTIM=TARGET TIME IN MINUTES
1535 DIFT=ABS(FDTIM-TTIM)
1536 IF(DIFT.LE.TW1)TIMEWT2=1.0
1537
1538 IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
1539 IF(FDTIM.LT.TTIM)TIMEWT2=(FDTIM-(TTIM-TW2))*TCONST
1540 IF(FDTIM.GT.TTIM)TIMEWT2=((TTIM+TW2)-FDTIM)*TCONST
1541 ENDIF
1542
1543 ! print *,'timewt2=',timewt2,ttim,fdtim
1544 if (ra(jjjn).ge.1. .and. rb(jjjn).ge.1. &
1545 .and.ra(jjjn).le.real(ide) .and. rb(jjjn).le.real(jde) &
1546 .and.timewt2.gt.0.) then
1547 if(lev_in_ob(jjjn).eq.1.)nsthis=nsthis+1
1548 if(plfo(jjjn).eq.1.)nsmetar=nsmetar+1
1549 if(plfo(jjjn).eq.2.)nsspeci=nsspeci+1
1550 if(plfo(jjjn).eq.3.)nsship=nsship+1
1551 if(plfo(jjjn).eq.4.)nssynop=nssynop+1
1552 if(plfo(jjjn).eq.5..and.lev_in_ob(jjjn).eq.1.) nstemp=nstemp+1
1553 if(plfo(jjjn).eq.6..and.lev_in_ob(jjjn).eq.1.) nspilot=nspilot+1
1554 if(plfo(jjjn).eq.7.)nssatwnds=nssatwnds+1
1555 if(plfo(jjjn).eq.8.)nssams=nssams+1
1556 if(plfo(jjjn).eq.9..and.lev_in_ob(jjjn).eq.1.) nsprofs=nsprofs+1
1557 endif
1558 enddo
1559 write(91,912)nsthis
1560 write(6,912)nsthis
1561 912 FORMAT('total obs on this grid=',i8)
1562 write(91,921)nsmetar
1563 write(6,921)nsmetar
1564 921 FORMAT('total metar obs on this grid=',i8)
1565 write(91,922)nsspeci
1566 write(6,922)nsspeci
1567 922 FORMAT('total special obs on this grid=',i8)
1568 write(91,923)nsship
1569 write(6,923)nsship
1570 923 FORMAT('total ship obs on this grid=',i8)
1571 write(91,924)nssynop
1572 write(6,924)nssynop
1573 924 FORMAT('total synop obs on this grid=',i8)
1574 write(91,925)nstemp
1575 write(6,925)nstemp
1576 925 FORMAT('total temp obs on this grid=',i8)
1577 write(91,926)nspilot
1578 write(6,926)nspilot
1579 926 FORMAT('total pilot obs on this grid=',i8)
1580 write(91,927)nssatwnds
1581 write(6,927)nssatwnds
1582 927 FORMAT('total sat-wind obs on this grid=',i8)
1583 write(91,928)nssams
1584 write(6,928)nssams
1585 928 FORMAT('total sams obs on this grid=',i8)
1586 write(91,929)nsprofs
1587 write(6,929)nsprofs
1588 929 FORMAT('total profiler obs on this grid=',i8)
1589 close(91)
1590 endif ! END OUTPUT OBS PER GRID EVERY HOUR
1591
1592
1593 ! INITIALIZE WEIGHTING ARRAYS TO ZERO
1594 DO I=its,ite
1595 DO K=1,kte
1596 WT(I,K)=0.0
1597 WT2ERR(I,K)=0.0
1598 ENDDO
1599 ENDDO
1600
1601 ! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
1602 ! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
1603 !
1604 ! COMPUTE P* AT OBS LOCATION (RA,RB). DO THIS AS SEPARATE VECTOR LOOP H
1605 ! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
1606
1607 ! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
1608 ! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
1609 ! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
1610 DO N=1,NSTAT
1611 IF(IVAR.GE.3)THEN
1612 PSURF(N)=ERRF(6,N)
1613 ELSE
1614 IF(IVAR.EQ.1)THEN
1615 PSURF(N)=ERRF(7,N) ! U-points
1616 ELSE
1617 PSURF(N)=ERRF(8,N) ! V-points
1618 ENDIF
1619 ENDIF
1620 ENDDO
1621
1622 ! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
1623 ! J-STRIP
1624
1625 MAXJ=J+IFIX(RINDX*fdob%RINFMX+0.99) !ajb
1626 MINJ=J-IFIX(RINDX*fdob%RINFMX+0.99) !ajb
1627
1628 ! jc comment out this? want to use obs beyond the domain?
1629 ! MAXJ=MIN0(JL-IGRID,MAXJ) !yliu
1630 ! MINJ=MAX0(1,MINJ) !yliu
1631
1632 n=1
1633
1634 !***********************************************************************
1635 DO nnn=1,NSTAT ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
1636 !***********************************************************************
1637 ! Soundings are consecutive obs, but they need to be treated as a single
1638 ! entity. Thus change the looping to nnn, with n defined separately.
1639
1640
1641 !yliu
1642 ! note for sfc data: nsndlev=1 and njcsnd=1
1643 nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
1644
1645 ! yliu start -- set together with the other parts
1646 ! test: do the sounding levels as individual obs
1647 ! nsndlev=1
1648 ! yliu end
1649 njcsnd=nsndlev
1650 ! set pob here, to be used later
1651 pob=varobs(5,n)
1652 ! print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
1653 ! print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
1654 ! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
1655 ! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP. THIS
1656 ! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
1657 ! ATION.
1658
1659 !yliu: Skip bad obs if it is sfc or single level sounding.
1660 !yliu: Before this (020918), a snd will be skipped if its first
1661 !yliu level has bad data- often true due to elevation
1662
1663 IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
1664 ! print *, " bad obs skipped"
1665
1666 ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
1667 ! print *, " skipped obs far away from this J-slice"
1668
1669 !----------------------------------------------------------------------
1670 ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
1671 !----------------------------------------------------------------------
1672
1673 ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
1674 FDTIM=XTIME-DTMIN
1675 ! TWINDO IS IN MINUTES:
1676 TW1=TWINDO/2.*60.
1677 TW2=TWINDO*60.
1678 TCONST=1./TW1
1679 TIMEWT=0.0
1680 TTIM=TIMEOB(N)*60.
1681 !***********TTIM=TARGET TIME IN MINUTES
1682 DIFT=ABS(FDTIM-TTIM)
1683 IF(DIFT.LE.TW1)TIMEWT=1.0
1684 IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
1685 IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
1686 IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
1687 ENDIF
1688
1689 ! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
1690 ! FOR THE VERTICAL WEIGHTING, WTSIG
1691
1692 ! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
1693 !ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
1694 !ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
1695
1696 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1697 rko(n) = errf(9,n) !ajb 20021210
1698 #endif
1699 KOB=nint(RKO(N))
1700 KOB=MIN0(kte,KOB)
1701 KOB=MAX0(1,KOB)
1702
1703 ! ASSIMILATE SURFACE LAYER DATA ON SIGMA
1704 IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
1705 DO K=1,kte
1706 WTSIG(K)=0.0
1707 ENDDO
1708 ! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
1709 ! WTSIG(1)=1.0
1710 ! WTSIG(2)=0.67
1711 ! WTSIG(3)=0.33
1712 ! KOMIN=3
1713 ! KOMAX=1
1714 ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
1715 ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
1716 ! fix this because kpblt at 1 and il is 0
1717 MAXI=IFIX(RA(N)+0.99+RINDX)
1718 MAXI=MIN0(ide-1,MAXI)
1719 MINI=IFIX(RA(N)-RINDX-0.99)
1720 MINI=MAX0(2,MINI)
1721 !yliu start
1722 ! use also obs outside of this domain -- surface obs
1723 ! if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
1724 ! & RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
1725 ! print *, " skipped obs far away from this domain"
1726 ! currently can use obs within this domain or ones very close to (1/3
1727 ! influence of radius in the coarse domain) this
1728 ! domain. In later case, use BC column value to approximate the model value
1729 ! at obs point -- ERRF need model field in errob.F !!
1730 if ( RA(N).GE.(0.-RINDX/3) &
1731 .and. RA(N).LE.float(ide)+RINDX/3 &
1732 .and. RB(N).GE.(0.-RINDX/3) &
1733 .and. RB(N).LE.float(jde)+RINDX/3) then
1734
1735 ! or use obs within this domain only
1736 ! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
1737 ! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
1738 ! print *, " skipped obs far outside of this domain"
1739 ! if(j.eq.3 .and. ivar.eq.3) then
1740 ! write(6,*) 'N = ',n,' exit 120 3'
1741 ! endif
1742 !yliu end
1743 !
1744 ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
1745 ! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
1746 ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
1747 RJ=FLOAT(J)
1748 RX=RJ-RB(N)
1749 ! WEIGHTS FOR THE 3-D VARIABLES
1750 ERFIVR=ERRF(IVAR,N)
1751 !
1752 !JM I will be local, because it indexes into PDOC, WT, and others
1753
1754 ! if((ivar.eq.1 .or. ivar.eq.3) .and. n.le.200) then
1755 ! write(6,'(a,i3,a,i3)')'SURF OBS NEAR: N = ',n,' nest = ',inest
1756 ! write(6,'(a,f10.3,a,f10.3,a,i3,a,i3,a,i3,a,i2)')
1757 ! $ ' RA =',RA(N),' RB =',RB(N),' J =',j,
1758 ! $ ' MINI =',MINI,' MAXI =',MAXI,' NEST =',inest
1759 ! endif
1760
1761 DO I=max0(its,MINI),min0(ite,MAXI)
1762
1763 RI=FLOAT(I)
1764 RY=RI-RA(N)
1765 RIS=RINDX*RINDX
1766 RSQ=RX*RX+RY*RY
1767 ! DPRIM=SQRT(RSQ)
1768 ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
1769 ! D=DPRIM+RINDX*DCON*ABS(PSBO(N)-PDOC(I,J))
1770 ! DSQ=D*D
1771 ! WTIJ=(RIS-DSQ)/(RIS+DSQ)
1772 wtij=(ris-rsq)/(ris+rsq)
1773 scratch = (abs(psurf(n)-.001*pbase(i,1))*fdob%DCON)
1774 pdfac=1.-AMIN1(1.0,scratch)
1775 wtij=wtij*pdfac
1776 WTIJ=AMAX1(0.0,WTIJ)
1777
1778 ! try making sfc obs weighting go thru pbl
1779 ! jc kpbl is at dot or cross only - need to interpolate?
1780 ! wtsig(1)=1.
1781 komax=max0(3,kpblt(i))
1782
1783 ! jc arbitrary check here
1784 IF (iprt) THEN
1785 if (kpblt(i).gt.25 .and. ktau.ne.0) &
1786 write(6,552)inest,i,j,kpblt(i)
1787 552 FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
1788 ENDIF
1789
1790 if(kpblt(i).gt.25) komax=3
1791 komin=1
1792 dk=float(komax)
1793
1794 do k=komin,komax
1795
1796 wtsig(k)=float(komax-k+1)/dk
1797 WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ
1798
1799 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ*WTSIG(K) &
1800 *WTSIG(K)*ERFIVR
1801
1802 ! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
1803 !
1804 ! write(6,'(a,i2,a,f8.3,a,f8.3,a,f8.3,a,f8.3,a,f8.3)')
1805 ! 'Surface obs, after: k = ',k, &
1806 ! ' WT2ERR = ',wt2err(i,k), &
1807 ! ' TIMEWT = ',timewt, &
1808 ! ' WTIJ = ',wtij, &
1809 ! ' WSIG = ',wtsig(k), &
1810 ! ' ERFIVR = ',erfivr
1811 ! endif
1812
1813 enddo
1814
1815 ENDDO
1816
1817 ! print *, " Surface "
1818
1819 endif ! end check for obs in domain
1820 ! END SURFACE-LAYER U OR V OBS NUDGING
1821
1822 ELSE
1823 ! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
1824 !
1825 ! print *,'in upper air section'
1826 ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
1827 ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
1828 ! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
1829 ! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
1830 !ajb SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
1831
1832 slope = (fdob%RINFMN-fdob%RINFMX)/(psurf(n)-fdob%PFREE)
1833
1834 RINFAC=SLOPE*POB+fdob%RINFMX-SLOPE*fdob%pfree
1835 RINFAC=AMAX1(RINFAC,fdob%RINFMN)
1836 RINFAC=AMIN1(RINFAC,fdob%RINFMX)
1837 !yliu: for multilevel upper-air data, take the maximum
1838 ! for the I loop.
1839 if(nsndlev.gt.1) RINFAC = fdob%RINFMX
1840 !yliu end
1841
1842 MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
1843 MAXI=MIN0(ide-IGRID,MAXI)
1844 MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
1845 MINI=MAX0(1,MINI)
1846
1847 ! yliu start
1848 ! use also obs outside of but close to this domain -- upr data
1849 ! if( RA(N).LT.(0.-RINFAC*RINDX)
1850 ! & .or. RA(N).GT.float(IL)+RINFAC*RINDX
1851 ! & .or. RB(N).LT.(0.-RINFAC*RINDX)
1852 ! & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then
1853 ! print *, " skipped obs far away from this I-range"
1854 ! currently can use obs within this domain or ones very close to (1/3
1855 ! influence of radius in the coarse domain) this
1856 ! domain. In later case, use BC column value to approximate the model value
1857 ! at obs point -- ERRF need model field in errob.F !!
1858
1859 !cc if (i.eq.39 .and. j.eq.34) then
1860 !cc write(6,*) 'RA(N) = ',ra(n)
1861 !cc write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
1862 !cc endif
1863 if( RA(N).GE.(0.-RINFAC*RINDX/3) &
1864 .and. RA(N).LE.float(ide)+RINFAC*RINDX/3 &
1865 .and. RB(N).GE.(0.-RINFAC*RINDX/3) &
1866 .and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
1867 ! or use obs within this domain only
1868 ! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
1869 ! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
1870 ! print *, " skipped obs far outside of this domain"
1871
1872 ! yliu end
1873 ! is this 2 needed here - kpbl not used?
1874 ! MINI=MAX0(2,MINI)
1875
1876 ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
1877 ! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
1878 ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
1879 RJ=FLOAT(J)
1880 RX=RJ-RB(N)
1881 ! WEIGHTS FOR THE 3-D VARIABLES
1882 !
1883 ERFIVR=ERRF(IVAR,N)
1884 ! jc
1885 nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
1886 ! yliu start
1887 ! test: do the sounding levels as individual obs
1888 ! nsndlev=1
1889 ! yliu end
1890 njcsnd=nsndlev
1891 !
1892 DO I=max0(its,MINI),min0(ite,MAXI)
1893 ! jc
1894 RI=FLOAT(I)
1895 RY=RI-RA(N)
1896 RIS=RINDX*RINFAC*RINDX*RINFAC
1897 RSQ=RX*RX+RY*RY
1898 ! yliu test: for upper-air data, keep D1 influence radii
1899 ! RIS=RIS /schnes /schnes
1900 WTIJ=(RIS-RSQ)/(RIS+RSQ)
1901 WTIJ=AMAX1(0.0,WTIJ)
1902 ! weight ob in vertical with +- 50 mb
1903 ! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
1904 if(nsndlev.eq.1) then
1905 rinprs=7.5
1906 else
1907 rinprs=3.0
1908 endif
1909 ! yliu end
1910 !
1911 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
1912 ! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
1913 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
1914
1915 if(nsndlev.eq.1)then
1916 !----------------------------------------------------------------------
1917 ! --- HANDLE 1-LEVEL OBSERVATIONS ---
1918 !----------------------------------------------------------------------
1919
1920 ! if(I.eq.MINI) print *, " Single snd "
1921 ! ERFIVR is the residual (difference) between the ob and the model
1922 ! at that point. We can analyze that residual up and down.
1923 ! First find komin for ob.
1924 !yliu start -- in the old code, komax and komin were reversed!
1925 do k=kte,1,-1
1926 pijk = .001*(pbase(i,k)+pp(i,k))
1927 ! print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
1928 if(pijk.ge.(pob+rinprs)) then
1929 komin=k
1930 go to 325
1931 endif
1932 enddo
1933 komin=1
1934 325 continue
1935 ! now find komax for ob
1936 do k=3,kte
1937 pijk = .001*(pbase(i,k)+pp(i,k))
1938 if(pijk.le.(pob-rinprs)) then
1939 komax=k
1940 go to 326
1941 endif
1942 enddo
1943 komax=kte ! yliu 20050706
1944 326 continue
1945
1946 ! yliu: single-level upper-air data will act either above or below the PBL top
1947 ! komax=min0(kpblt(i), komax)
1948 if(komax.gt.kpblt(i)) komin=max0(kpblt(i), komin)
1949 if(komin.lt.kpblt(i)) komax=min0(kpblt(i), komax)
1950 ! yliu end
1951 !
1952 ! print *,'1 level, komin,komax=',komin,komax
1953 ! if(i.eq.MINI) then
1954 ! print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
1955 ! ERFIVR=0
1956 ! endif
1957 do k=1,kte
1958 reserf(k)=0.0
1959 wtsig(k)=0.0
1960 enddo
1961 !yliu end
1962
1963 !cc if (i.eq.39 .and. j.eq.34) then
1964 !cc write(6,*) ' komin = ', komin,' komax = ',komax
1965 !cc endif
1966
1967 do k=komin,komax
1968 pijk = .001*(pbase(i,k)+pp(i,k))
1969 reserf(k)=erfivr
1970 wtsig(k)=1.-abs(pijk-pob)/rinprs
1971 wtsig(k)=amax1(wtsig(k),0.0)
1972 ! print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
1973 ! Now calculate WT and WT2ERR for each i,j,k point cajb
1974 WT(I,K)=WT(I,K)+TIMEWT*WTIJ*wtsig(k)
1975
1976 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ* &
1977 reserf(k)*wtsig(k)*wtsig(k)
1978 enddo
1979
1980 else
1981 !----------------------------------------------------------------------
1982 ! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
1983 !----------------------------------------------------------------------
1984 !yliu start
1985 ! if(I.eq.MINI) print *, " Multi-level snd "
1986 ! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
1987 ! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
1988 if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
1989 IF (iprt) THEN
1990 print *, "n = ",n,"nsndlev = ",nsndlev
1991 print *, "nlevs_ob,lev_in_ob", &
1992 nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
1993 print *, "in nudobs.F: sounding level messed up, stopping"
1994 ENDIF
1995 stop
1996 endif
1997 !yliu end
1998 ! This is for a multi-level observation
1999 ! The trick here is that the sounding is "one ob". You don't
2000 ! want multiple levels to each be treated like separate
2001 ! and independent observations.
2002 ! At each i,j want to interpolate sounding to the model levels at that
2003 ! particular point.
2004 komin=1
2005 komax=kte-2
2006 ! this loop goes to 1501
2007 ! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
2008 !yliu start
2009 do k=1,kte
2010 reserf(k)=0.0
2011 wtsig(k)=0.0
2012 enddo
2013 !yliu end
2014
2015 do k=komax,komin,-1
2016
2017 pijk = .001*(pbase(i,k)+pp(i,k))
2018
2019 ! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
2020 if(pijk.gt.varobs(5,n)) then
2021 go to 1501
2022 endif
2023
2024 ! if sigma level pressure is .lt. than the highest ob level, don't interpolate
2025 if(pijk.le.varobs(5,n+nsndlev-1)) then
2026 go to 1501
2027 endif
2028
2029 ! now interpolate sounding to this k
2030 ! yliu start-- recalculate WTij for each k-level
2031 !ajb SLOPE = (fdob%RINFMN-fdob%RINFMX)/(pdoc(i,j)+PTOP-fdob%PFREE)
2032 slope = (fdob%RINFMN-fdob%RINFMX)/ (.001*pbase(i,1)-fdob%PFREE)
2033 RINFAC=SLOPE*pijk+fdob%RINFMX-SLOPE*fdob%PFREE
2034 RINFAC=AMAX1(RINFAC,fdob%RINFMN)
2035 RINFAC=AMIN1(RINFAC,fdob%RINFMX)
2036 RIS=RINDX*RINFAC*RINDX*RINFAC
2037 RSQ=RX*RX+RY*RY
2038
2039 ! for upper-air data, keep D1 influence radii
2040 ! RIS=RIS /schnes /schnes
2041 WTIJ=(RIS-RSQ)/(RIS+RSQ)
2042 WTIJ=AMAX1(0.0,WTIJ)
2043 ! yliu end
2044
2045 ! this loop goes to 1503
2046 do nn=2,nsndlev
2047 ! only set pobhi if varobs(ivar) is ok
2048 pobhi=-888888.
2049
2050 if(varobs(ivar,n+nn-1).gt.-800000. &
2051 .and. varobs(5,n+nn-1).gt.-800000.) then
2052 pobhi=varobs(5,n+nn-1)
2053 nhi=n+nn-1
2054 if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.20.) then
2055 go to 1502 ! within 200mb of obs height
2056 endif
2057 endif
2058
2059 enddo
2060
2061 ! did not find any ob above within 100 mb, so jump out
2062 go to 1501
2063 1502 continue
2064
2065 nlo=nhi-1
2066 do nnjc=nhi-1,n,-1
2067 if(varobs(ivar,nnjc).gt.-800000. &
2068 .and. varobs(5,nnjc).gt.-800000.) then
2069 poblo=varobs(5,nnjc)
2070 nlo=nnjc
2071 if(poblo.gt.pijk .and. abs(poblo-pijk).lt.20.) then
2072 go to 1505 ! within 200mb of obs height
2073 endif
2074 endif
2075 enddo
2076 !yliu end --
2077
2078 ! did not find any ob below within 200 mb, so jump out
2079 go to 1501
2080 1505 continue
2081
2082 ! interpolate to model level
2083 pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
2084 reserf(k)=errf(ivar,nlo)+ &
2085 (errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
2086 wtsig(k)=1.
2087
2088 1501 continue
2089
2090 ! now calculate WT and WT2ERR for each i,j,k point cajb
2091 WT(I,K)=WT(I,K)+TIMEWT*WTIJ*wtsig(k)
2092
2093 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ* &
2094 reserf(k)*wtsig(k)*wtsig(k)
2095
2096 ! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
2097 !
2098 ! if(wt(i,k) .ne. 0.0) then
2099 ! scratch = WT2ERR(I,K)/WT(I,K)
2100 ! else
2101 ! scratch = 999.
2102 ! endif
2103 !
2104 ! write(6,'(a,i2,a,f8.3,a,f4.2,a,f7.4,a,f4.2,a,f5.3,a,f7.4)')
2105 ! $ 'Multi-level obs: k = ',k,
2106 ! $ ' WT2ERR = ',wt2err(i,k),
2107 ! $ ' WTIJ = ',wtij,
2108 ! $ ' RSF = ',reserf(k),
2109 ! $ ' WSIG = ',wtsig(k),
2110 ! $ ' WT = ',wt(i,k),
2111 ! $ ' W2EOWT = ',scratch
2112 ! endif
2113
2114
2115 ! end do k
2116 enddo ! enddo k levels
2117 ! end multi-levels
2118 endif ! end if(nsndlev.eq.1)
2119 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2120 ! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
2121 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2122 !
2123 ENDDO ! END DO MINI,MAXI LOOP
2124
2125 endif ! check for obs in domain
2126
2127 ! END OF NUDGING TO OBS ON PRESSURE LEVELS
2128
2129 ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
2130
2131 !----------------------------------------------------------------------
2132 ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
2133 !----------------------------------------------------------------------
2134
2135 ! n=n+1
2136 n=n+njcsnd
2137
2138 !yliu 1202 continue
2139 if(n.gt.nstat)then
2140 ! print *,'n,nstat=',n,nstat,ivar,j
2141 go to 1203
2142 endif
2143 ! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
2144
2145 !***********************************************************************
2146 ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
2147 !***********************************************************************
2148
2149 1203 continue
2150
2151 ! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
2152 ! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
2153 ! THE ATEN ARRAY
2154 ! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
2155 ! THEY ARE USED BELOW IN THE DENOMINATOR.
2156 DO K=kts,kte
2157 DO I=its,ite
2158 IF(WT(I,K).EQ.0)THEN
2159 WT2ERR(I,K)=0.0
2160 ENDIF
2161 IF(WT(I,K).EQ.0)THEN
2162 WT(I,K)=1.0
2163 ENDIF
2164 ENDDO
2165 ENDDO
2166
2167 126 CONTINUE
2168
2169 IF(IVAR.GE.3)GOTO 170
2170 ! this is for u,v
2171 ! 3-D DOT POINT TENDENCIES
2172
2173 ! Calculate scales for converting nudge factor from u (v)
2174 ! to rho_u (or rho_v) units.
2175
2176 call calc_rcouple_scales(mu,msf,rscale,ims,ime,its,ite)
2177
2178 DO K=1,kte
2179
2180 DO I=i_s,i_e
2181
2182 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2183 W2EOWT=WT2ERR(I,K)/WT(I,K)
2184 ELSE
2185 W2EOWT=SAVWT(IPL,I,K)
2186 ENDIF
2187
2188 ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2189 ! scratch = GIV*RSCALE(I)*W2EOWT*fdob%TFACI*ISWIND*GFACTOR
2190 ! write(6,*) 'ATEN calc: k = ',k
2191 ! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
2192 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2193 ! $ ' W2EOWT = ',w2eowt
2194 ! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
2195 ! $ ' GFACTOR = ',gfactor
2196 ! endif
2197 !
2198 ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2199 ! scratch = GIV*RSCALE(I)*W2EOWT*fdob%TFACI*ISWIND*GFACTOR
2200 ! write(6,*) 'ATEN calc: k = ',k
2201 ! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
2202 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2203 ! $ ' W2EOWT = ',w2eowt
2204 ! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
2205 ! $ ' GFACTOR = ',gfactor
2206 ! endif
2207
2208 ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
2209 *W2EOWT*fdob%TFACI &
2210 *ISWIND *GFACTOR !yliu *GFACTOR
2211
2212 ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2213 ! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
2214 ! endif
2215 ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2216 ! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
2217 ! endif
2218
2219 ENDDO
2220 ENDDO
2221
2222 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2223 DO K=1,kte
2224 DO I=its,ite
2225 SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2226 ENDDO
2227 ENDDO
2228 ENDIF
2229
2230 RETURN
2231
2232 170 CONTINUE
2233
2234 ! 3-D CROSS-POINT TENDENCIES
2235 ! this is for t (ivar=3) and q (ivsr=4)
2236 IF(3-IVAR.LT.0)THEN
2237 GITQ=GIQ
2238 ELSE
2239 GITQ=GIT
2240 ENDIF
2241 IF(3-IVAR.LT.0)THEN
2242 ISTQ=ISMOIS
2243 ELSE
2244 ISTQ=ISTEMP
2245 ENDIF
2246
2247 DO K=1,kte
2248 DO I=i_s,i_e
2249 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2250 W2EOWT=WT2ERR(I,K)/WT(I,K)
2251 ELSE
2252 W2EOWT=SAVWT(IPL,I,K)
2253 ENDIF
2254
2255 ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
2256 ! scratch = GITQ*MU(I)*W2EOWT*fdob%TFACI*ISTQ*GFACTOR
2257 ! write(6,*) 'ATEN calc: k = ',k
2258 ! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
2259 ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
2260 ! $ ' W2EOWT = ',w2eowt
2261 ! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
2262 ! $ ' GFACTOR = ',gfactor
2263 ! endif
2264 !
2265 ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2266 ! scratch = GITQ*MU(I)*W2EOWT*fdob%TFACI*ISTQ*GFACTOR
2267 ! write(6,*) 'ATEN calc: k = ',k
2268 ! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
2269 ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
2270 ! $ ' W2EOWT = ',w2eowt
2271 ! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
2272 ! $ ' GFACTOR = ',gfactor
2273 ! endif
2274
2275 ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
2276 *W2EOWT*fdob%TFACI*ISTQ *GFACTOR !yliu *GFACTOR
2277
2278 ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
2279 ! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
2280 ! endif
2281 ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2282 ! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
2283 ! endif
2284
2285 ENDDO
2286 ENDDO
2287
2288 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
2289 DO K=1,kte
2290 DO I=its,ite
2291 SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2292 ENDDO
2293 ENDDO
2294 ENDIF
2295
2296 RETURN
2297 END SUBROUTINE nudob
2298
2299 SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
2300 !-----------------------------------------------------------------------
2301 IMPLICIT NONE
2302 !-----------------------------------------------------------------------
2303
2304 INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
2305 INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
2306 REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
2307 REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
2308 REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
2309
2310 ! Local variables
2311 integer :: i
2312
2313 ! Calculate scales to be used for producing rho-coupled nudging factors.
2314 do i = its,ite
2315 rscale(i) = a(i)/msf(i)
2316 enddo
2317
2318 RETURN
2319 END SUBROUTINE calc_rcouple_scales
2320
2321 !ajb: Not used
2322 SUBROUTINE set_real_array(rscale, value, ims,ime, its,ite)
2323 !-----------------------------------------------------------------------
2324 IMPLICIT NONE
2325 !-----------------------------------------------------------------------
2326
2327 INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
2328 INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
2329 REAL, INTENT(IN) :: value ! Constant array value
2330 REAL, INTENT(OUT) :: rscale( ims:ime ) ! Output array
2331
2332 ! Local variables
2333 integer :: i
2334
2335 ! Set array to constant value
2336 do i = its,ite
2337 rscale(i) = value
2338 enddo
2339
2340 RETURN
2341 END SUBROUTINE set_real_array
2342
2343 !ajb: Not used
2344 SUBROUTINE calc_pottemp_scales(ivar, rcp, pb, p, tscale, &
2345 ims,ime, its,ite, &
2346 kms,kme, kts,kte)
2347 !-----------------------------------------------------------------------
2348 IMPLICIT NONE
2349 !-----------------------------------------------------------------------
2350
2351 INTEGER, INTENT(IN) :: ims,ime, kms,kme ! Memory dimensions
2352 INTEGER, INTENT(IN) :: its,ite, kts,kte ! Tile dimensions
2353 INTEGER, INTENT(IN) :: ivar ! Variable identifier
2354 REAL, INTENT(IN) :: rcp ! Constant (2./7.)
2355 REAL, INTENT(IN) :: pb(ims:ime, kms:kme) ! Base pressure (Pa) array
2356 REAL, INTENT(IN) :: p(ims:ime, kms:kme) ! Pressure pert. (Pa) array
2357 REAL, INTENT(OUT) :: tscale(ims:ime, kms:kme) ! Scales for pot. temp.
2358 ! Local variables
2359 integer :: i,k
2360
2361 if(ivar.eq.3) then
2362
2363 ! Calculate scales to be used for producing potential temperature nudging factors.
2364 do k = kts,kte
2365 do i = its,ite
2366 tscale(i,k) = ( 1000000. / ( pb(i,k)+p(i,k)) )**rcp
2367 enddo
2368 enddo
2369 else
2370 ! Set to 1. for moisture scaling.
2371 do k = kts,kte
2372 do i = its,ite
2373 tscale(i,k) = 1.0
2374 enddo
2375 enddo
2376 endif
2377
2378 RETURN
2379 END SUBROUTINE calc_pottemp_scales
2380 #endif
2381
2382 END MODULE module_fddaobs_rtfdda
2383