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