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, msf, 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) :: msf(ims:ime) ! Map scale (only used for vars u & v)
1335 INTEGER, intent(in) :: iswind ! Nudge flag for wind
1336 INTEGER, intent(in) :: istemp ! Nudge flag for temperature
1337 INTEGER, intent(in) :: ismois ! Nudge flag for moisture
1338 REAL, intent(in) :: giv ! Coefficient for wind
1339 REAL, intent(in) :: git ! Coefficient for temperature
1340 REAL, intent(in) :: giq ! Coefficient for moisture
1341 REAL, INTENT(INOUT) :: aten( ims:ime, kms:kme)
1342 REAL, INTENT(INOUT) :: savwt( nndgv, ims:ime, kms:kme )
1343 INTEGER, INTENT(IN) :: kpblt(its:ite)
1344 INTEGER, INTENT(IN) :: nscan ! number of scans
1345 LOGICAL, INTENT(IN) :: iprt ! print flag
1346
1347 ! Local variables
1348 integer :: mm(maxdom)
1349 real :: ra(niobf)
1350 real :: rb(niobf)
1351 real :: psurf(niobf)
1352 real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
1353 real :: rscale(ims:ime) ! For converting to rho-coupled units.
1354 ! real :: tscale(ims:ime,kms:kme) ! For converting to potential temp. units.
1355 real :: reserf(100)
1356 character*40 name
1357 character*3 chr_hr
1358
1359 !*** DECLARATIONS FOR IMPLICIT NONE
1360 integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
1361 integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
1362 integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
1363 integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
1364 integer :: komin,komax,nn,nhi,nlo,nnjc
1365 integer :: i_s,i_e
1366 integer :: istq
1367 real :: gfactor,rfactor,gridx,gridy,rindx,schnes,ris
1368 real :: grfacx,grfacy
1369 real :: fdtim,tw1,tw2,tconst,timewt,timewt2,ttim,dift,pob
1370 real :: ri,rj,rx,ry,rsq,wtij,pdfac,erfivr,dk,slope,rinfac
1371 real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
1372
1373 real :: scratch
1374
1375 ! print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
1376 ! if(ivar.ne.4)return
1377 !yliu start -- for multi-scans: NSCAN=0: original
1378 ! NSCAN=1: added a scan with a larger Ri and smaller G
1379 ! if(NSCAN.ne.0 .and. NSCAN.ne.1) stop
1380 ! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
1381 if(NSCAN.ne.0) then
1382 IF (iprt) write(6,*) 'SAVWT must be resized for NSCAN=1'
1383 stop
1384 endif
1385 IPLO=0 + NSCAN*4
1386 GFACTOR=1. + NSCAN*(-1. + 0.33333)
1387 RFACTOR=1. + NSCAN*(-1. + 3.0)
1388 !yliu end
1389 ! jc
1390
1391 ! return if too close to j boundary
1392 if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
1393 ! write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
1394 ! $ ' too close to boundary.'
1395 return
1396 endif
1397 if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
1398 ! write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
1399 ! $ ' too close to boundary.'
1400 return
1401 endif
1402
1403 ! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
1404 ICUT=0
1405 IF(INEST.GT.1)ICUT=1
1406 i_s = max0(2+icut,its)
1407 i_e = min0(ide-1-icut,ite)
1408
1409 IPL=IVAR + IPLO !yliu +IPLO
1410
1411 ! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
1412
1413 INPF=(IRATIO**LEVIDN(INEST))*NPFI
1414 INFR=(IRATIO**LEVIDN(INEST))*IONF
1415
1416 GRIDX=0.0
1417 GRIDY=0.0
1418 IGRID=0
1419 IF(IVAR.GE.3)THEN
1420 GRIDX=0.5
1421 GRIDY=0.5
1422 IGRID=1
1423 ELSEIF(IVAR.eq.1) THEN
1424 GRIDY=0.5
1425 IGRID=1
1426 ELSEIF(IVAR.eq.2) THEN
1427 GRIDX=0.5
1428 IGRID=1
1429 ENDIF
1430
1431 ! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
1432 ! KILOMETERS TO GRID LENGTHS, RINDX
1433
1434 RINDX=RINXY*1000./DX * RFACTOR !yliu *RFACTOR
1435
1436 ! jc
1437 ! make horizontal radius vary per nest
1438 ! rindx=rindx/float(inest)
1439 ! yliu test1 -- En 3, 4
1440 ! rindx=rindx/float(3**(in-1)) !YLIU
1441 ! jc
1442 ! make horizontal radius vary per nest
1443 ! schnes=1/float(inest) !JC
1444 ! yliu test1 -- En 3, 4 !YLIU
1445 schnes=1/float(3**(inest-1)) !JC
1446 ! reduce the Rinf in the nested grid proportionally
1447 rindx=rindx*schnes
1448 ! rinfmn =1., rinfmx=2., pfree=50 in param.F
1449 ! yliu test: for upper-air data, use larger influence radii
1450 ! Essentially increase the slope -- the same radii
1451 ! at 500mb and above as the coarse mesh and the
1452 ! same small radii at sfc as that for sfc obs
1453 fdob%rinfmx=2. *1.0 /schnes !YLIU
1454 ! rinfmx=1.2*1.0 /schnes !YLIU
1455 ! jc
1456 RIS=RINDX*RINDX
1457 IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
1458 IF(MOD(KTAU,INFR).NE.0)GOTO 126
1459 5 CONTINUE
1460 IF (iprt) THEN
1461 IF(J.EQ.10) write(6,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
1462 ENDIF
1463 6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,', &
1464 'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2, &
1465 ' rindx=',f4.1)
1466
1467 ! SET RA AND RB
1468 IF(INEST.EQ.1) THEN
1469
1470 ! SET RA AND RB FOR THE COARSE MESH...
1471 DO N=1,NSTAT
1472 RA(N)=RIO(N)-GRIDX
1473 RB(N)=RJO(N)-GRIDY
1474 ENDDO
1475
1476 ELSE
1477
1478 ! SET RA AND RB FOR THE FINE MESH CASE...
1479 DO N=1,NSTAT
1480
1481 ! COMPUTE THE OBS LOCATION WITH RESPECT TO THIS MESH (INEST)...
1482 NML=INEST
1483 MM(LEVIDN(INEST)+1)=INEST
1484 ! WORKING TOWARD COARSER MESHES, DETERMINE THE HIERARCHY OF MOTHER
1485 ! MESHES WITH RESPECT TO EACH MOTHER MESH STARTING AT MESH "INEST"...
1486 ! THAT IS, DETERMINE ITS MOTHER, GRANDMOTHER, GREAT-GRANDMOTHER, ETC.
1487 ! OUT TO THE COARSE GRID MESH (INEST=1).
1488 ! LEVIDN HOLDS THE NEST LEVEL AND PARID HOLDS THE MOTHER MESH FOR EACH
1489 ! GRID (E.G., FOR 3 MESHES AND 2 NEST LEVELS, INEST=1 IS THE COARSE GRID
1490 ! MESH, INEST=2 HAS LEVIDN(2)=1 AND PARID(2)=1, AND INEST=3 HAS LEVIDN(3)=2
1491 ! AND PARID(3)=2...)
1492 DO NNL=LEVIDN(INEST),1,-1
1493 MM(NNL)=PARID(NML)
1494 NML=MM(NNL)
1495 ENDDO
1496
1497 ! MM(1) MUST BE THE COARSE GRID MESH (INEST=0)
1498
1499 IF(MM(1).NE.1) then
1500 IF (iprt) write(6,*) 'stop 21 in nudob: inest = ',inest
1501 STOP 21
1502 ENDIF
1503 RA(N)=RIO(N)
1504 RB(N)=RJO(N)
1505 DO NNL=1,LEVIDN(INEST)
1506 GRFACX=0.
1507 GRFACY=0.
1508 ! COMPUTE THE OBS LOCATION WITH RESPECT TO THE INNER GRID IN DOT-POINT
1509 ! SPACE (GRID=0.). WHEN WE REACH MESH IN, THEN USE "GRID" TO GO TO
1510 ! CROSS OR DOT DEPENDING ON THE VARIABLE...
1511 IF(NNL.EQ.LEVIDN(INEST)) THEN
1512 GRFACX=GRIDX
1513 GRFACY=GRIDY
1514 ENDIF
1515
1516 RA(N)=(RA(N)-FLOAT(i_parent_start(MM(NNL+1))))* &
1517 FLOAT(IRATIO)+1.-GRFACX
1518 RB(N)=(RB(N)-FLOAT(j_parent_start(MM(NNL+1))))* &
1519 FLOAT(IRATIO)+1.-GRFACY
1520
1521 ENDDO
1522
1523 ENDDO ! END LOOP OVER OBS STATIONS FOR FINE MESH CASE
1524
1525 ENDIF ! END SECTION FOR SETTING RA AND RB
1526
1527
1528 ! OUTPUT OBS PER GRID EVERY HOUR
1529 if ( mod(xtime,60.).gt.56. .and. ivar.eq.3 .and. j.eq.10) then
1530 IF (iprt) print *,'outputting obs number on grid ', &
1531 inest,' at time=',xtime
1532 write(chr_hr(1:3),'(i3)')nint(xtime/60.)
1533 if(chr_hr(1:1).eq.' ')chr_hr(1:1)='0'
1534 if(chr_hr(2:2).eq.' ')chr_hr(2:2)='0'
1535 IF (iprt) print *,'chr_hr=',chr_hr(1:3),nint(xtime/60.)
1536 open(91,file= &
1537 'obs_g'//char(inest+ichar('0'))//'_'//chr_hr(1:3), &
1538 form='FORMATted',status='unknown')
1539 write(91,911)nstat
1540 write(6,911)nstat
1541 911 FORMAT('total obs=',i8)
1542 nsthis=0
1543 nsmetar=0
1544 nsspeci=0
1545 nsship=0
1546 nssynop=0
1547 nstemp=0
1548 nspilot=0
1549 nssatwnds=0
1550 nssams=0
1551 nsprofs=0
1552 ! print *,'ide,jde=',ide,jde
1553 do jjjn=1,nstat
1554 ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
1555 FDTIM=XTIME-DTMIN
1556 ! CONVERT TWINDO AND TIMEOB FROM HOURS TO MINUTES:
1557 TW1=TWINDO/2.*60.
1558 TW2=TWINDO*60.
1559 TCONST=1./TW1
1560 TIMEWT2=0.0
1561 TTIM=TIMEOB(jjjn)*60.
1562 !***********TTIM=TARGET TIME IN MINUTES
1563 DIFT=ABS(FDTIM-TTIM)
1564 IF(DIFT.LE.TW1)TIMEWT2=1.0
1565
1566 IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
1567 IF(FDTIM.LT.TTIM)TIMEWT2=(FDTIM-(TTIM-TW2))*TCONST
1568 IF(FDTIM.GT.TTIM)TIMEWT2=((TTIM+TW2)-FDTIM)*TCONST
1569 ENDIF
1570
1571 ! print *,'timewt2=',timewt2,ttim,fdtim
1572 if (ra(jjjn).ge.1. .and. rb(jjjn).ge.1. &
1573 .and.ra(jjjn).le.real(ide) .and. rb(jjjn).le.real(jde) &
1574 .and.timewt2.gt.0.) then
1575 if(lev_in_ob(jjjn).eq.1.)nsthis=nsthis+1
1576 if(plfo(jjjn).eq.1.)nsmetar=nsmetar+1
1577 if(plfo(jjjn).eq.2.)nsspeci=nsspeci+1
1578 if(plfo(jjjn).eq.3.)nsship=nsship+1
1579 if(plfo(jjjn).eq.4.)nssynop=nssynop+1
1580 if(plfo(jjjn).eq.5..and.lev_in_ob(jjjn).eq.1.) nstemp=nstemp+1
1581 if(plfo(jjjn).eq.6..and.lev_in_ob(jjjn).eq.1.) nspilot=nspilot+1
1582 if(plfo(jjjn).eq.7.)nssatwnds=nssatwnds+1
1583 if(plfo(jjjn).eq.8.)nssams=nssams+1
1584 if(plfo(jjjn).eq.9..and.lev_in_ob(jjjn).eq.1.) nsprofs=nsprofs+1
1585 endif
1586 enddo
1587 write(91,912)nsthis
1588 write(6,912)nsthis
1589 912 FORMAT('total obs on this grid=',i8)
1590 write(91,921)nsmetar
1591 write(6,921)nsmetar
1592 921 FORMAT('total metar obs on this grid=',i8)
1593 write(91,922)nsspeci
1594 write(6,922)nsspeci
1595 922 FORMAT('total special obs on this grid=',i8)
1596 write(91,923)nsship
1597 write(6,923)nsship
1598 923 FORMAT('total ship obs on this grid=',i8)
1599 write(91,924)nssynop
1600 write(6,924)nssynop
1601 924 FORMAT('total synop obs on this grid=',i8)
1602 write(91,925)nstemp
1603 write(6,925)nstemp
1604 925 FORMAT('total temp obs on this grid=',i8)
1605 write(91,926)nspilot
1606 write(6,926)nspilot
1607 926 FORMAT('total pilot obs on this grid=',i8)
1608 write(91,927)nssatwnds
1609 write(6,927)nssatwnds
1610 927 FORMAT('total sat-wind obs on this grid=',i8)
1611 write(91,928)nssams
1612 write(6,928)nssams
1613 928 FORMAT('total sams obs on this grid=',i8)
1614 write(91,929)nsprofs
1615 write(6,929)nsprofs
1616 929 FORMAT('total profiler obs on this grid=',i8)
1617 close(91)
1618 endif ! END OUTPUT OBS PER GRID EVERY HOUR
1619
1620
1621 ! INITIALIZE WEIGHTING ARRAYS TO ZERO
1622 DO I=its,ite
1623 DO K=1,kte
1624 WT(I,K)=0.0
1625 WT2ERR(I,K)=0.0
1626 ENDDO
1627 ENDDO
1628
1629 ! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
1630 ! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
1631 !
1632 ! COMPUTE P* AT OBS LOCATION (RA,RB). DO THIS AS SEPARATE VECTOR LOOP H
1633 ! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
1634
1635 ! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
1636 ! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
1637 ! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
1638 DO N=1,NSTAT
1639 IF(IVAR.GE.3)THEN
1640 PSURF(N)=ERRF(6,N)
1641 ELSE
1642 IF(IVAR.EQ.1)THEN
1643 PSURF(N)=ERRF(7,N) ! U-points
1644 ELSE
1645 PSURF(N)=ERRF(8,N) ! V-points
1646 ENDIF
1647 ENDIF
1648 ENDDO
1649
1650 ! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
1651 ! J-STRIP
1652
1653 MAXJ=J+IFIX(RINDX*fdob%RINFMX+0.99) !ajb
1654 MINJ=J-IFIX(RINDX*fdob%RINFMX+0.99) !ajb
1655
1656 ! jc comment out this? want to use obs beyond the domain?
1657 ! MAXJ=MIN0(JL-IGRID,MAXJ) !yliu
1658 ! MINJ=MAX0(1,MINJ) !yliu
1659
1660 n=1
1661
1662 !***********************************************************************
1663 DO nnn=1,NSTAT ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
1664 !***********************************************************************
1665 ! Soundings are consecutive obs, but they need to be treated as a single
1666 ! entity. Thus change the looping to nnn, with n defined separately.
1667
1668
1669 !yliu
1670 ! note for sfc data: nsndlev=1 and njcsnd=1
1671 nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
1672
1673 ! yliu start -- set together with the other parts
1674 ! test: do the sounding levels as individual obs
1675 ! nsndlev=1
1676 ! yliu end
1677 njcsnd=nsndlev
1678 ! set pob here, to be used later
1679 pob=varobs(5,n)
1680
1681 ! print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
1682 ! print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
1683 ! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
1684 ! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP. THIS
1685 ! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
1686 ! ATION.
1687
1688 !yliu: Skip bad obs if it is sfc or single level sounding.
1689 !yliu: Before this (020918), a snd will be skipped if its first
1690 !yliu level has bad data- often true due to elevation
1691
1692 IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
1693 ! print *, " bad obs skipped"
1694
1695 ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
1696 ! print *, " skipped obs far away from this J-slice"
1697
1698 !----------------------------------------------------------------------
1699 ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
1700 !----------------------------------------------------------------------
1701
1702 ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
1703 FDTIM=XTIME-DTMIN
1704 ! TWINDO IS IN MINUTES:
1705 TW1=TWINDO/2.*60.
1706 TW2=TWINDO*60.
1707 TCONST=1./TW1
1708 TIMEWT=0.0
1709 TTIM=TIMEOB(N)*60.
1710 !***********TTIM=TARGET TIME IN MINUTES
1711 DIFT=ABS(FDTIM-TTIM)
1712 IF(DIFT.LE.TW1)TIMEWT=1.0
1713 IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
1714 IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
1715 IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
1716 ENDIF
1717
1718 ! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
1719 ! FOR THE VERTICAL WEIGHTING, WTSIG
1720
1721 ! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
1722 !ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
1723 !ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
1724
1725 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1726 rko(n) = errf(9,n) !ajb 20021210
1727 #endif
1728 KOB=nint(RKO(N))
1729 KOB=MIN0(kte,KOB)
1730 KOB=MAX0(1,KOB)
1731
1732 ! ASSIMILATE SURFACE LAYER DATA ON SIGMA
1733 IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
1734 DO K=1,kte
1735 WTSIG(K)=0.0
1736 ENDDO
1737 ! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
1738 ! WTSIG(1)=1.0
1739 ! WTSIG(2)=0.67
1740 ! WTSIG(3)=0.33
1741 ! KOMIN=3
1742 ! KOMAX=1
1743 ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
1744 ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
1745 ! fix this because kpblt at 1 and il is 0
1746 MAXI=IFIX(RA(N)+0.99+RINDX)
1747 MAXI=MIN0(ide-1,MAXI)
1748 MINI=IFIX(RA(N)-RINDX-0.99)
1749 MINI=MAX0(2,MINI)
1750 !yliu start
1751 ! use also obs outside of this domain -- surface obs
1752 ! if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
1753 ! & RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
1754 ! print *, " skipped obs far away from this domain"
1755 ! currently can use obs within this domain or ones very close to (1/3
1756 ! influence of radius in the coarse domain) this
1757 ! domain. In later case, use BC column value to approximate the model value
1758 ! at obs point -- ERRF need model field in errob.F !!
1759 if ( RA(N).GE.(0.-RINDX/3) &
1760 .and. RA(N).LE.float(ide)+RINDX/3 &
1761 .and. RB(N).GE.(0.-RINDX/3) &
1762 .and. RB(N).LE.float(jde)+RINDX/3) then
1763
1764 ! or use obs within this domain only
1765 ! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
1766 ! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
1767 ! print *, " skipped obs far outside of this domain"
1768 ! if(j.eq.3 .and. ivar.eq.3) then
1769 ! write(6,*) 'N = ',n,' exit 120 3'
1770 ! endif
1771 !yliu end
1772 !
1773 ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
1774 ! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
1775 ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
1776 RJ=FLOAT(J)
1777 RX=RJ-RB(N)
1778 ! WEIGHTS FOR THE 3-D VARIABLES
1779 ERFIVR=ERRF(IVAR,N)
1780 !
1781 !JM I will be local, because it indexes into PDOC, WT, and others
1782
1783 ! if((ivar.eq.1 .or. ivar.eq.3) .and. n.le.200) then
1784 ! write(6,'(a,i3,a,i3)')'SURF OBS NEAR: N = ',n,' nest = ',inest
1785 ! write(6,'(a,f10.3,a,f10.3,a,i3,a,i3,a,i3,a,i2)')
1786 ! $ ' RA =',RA(N),' RB =',RB(N),' J =',j,
1787 ! $ ' MINI =',MINI,' MAXI =',MAXI,' NEST =',inest
1788 ! endif
1789
1790 DO I=max0(its,MINI),min0(ite,MAXI)
1791
1792 RI=FLOAT(I)
1793 RY=RI-RA(N)
1794 RIS=RINDX*RINDX
1795 RSQ=RX*RX+RY*RY
1796 ! DPRIM=SQRT(RSQ)
1797 ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
1798 ! D=DPRIM+RINDX*DCON*ABS(PSBO(N)-PDOC(I,J))
1799 ! DSQ=D*D
1800 ! WTIJ=(RIS-DSQ)/(RIS+DSQ)
1801 wtij=(ris-rsq)/(ris+rsq)
1802 scratch = (abs(psurf(n)-.001*pbase(i,1))*fdob%DCON)
1803 pdfac=1.-AMIN1(1.0,scratch)
1804 wtij=wtij*pdfac
1805 WTIJ=AMAX1(0.0,WTIJ)
1806
1807 ! try making sfc obs weighting go thru pbl
1808 ! jc kpbl is at dot or cross only - need to interpolate?
1809 ! wtsig(1)=1.
1810 komax=max0(3,kpblt(i))
1811
1812 ! jc arbitrary check here
1813 IF (iprt) THEN
1814 if (kpblt(i).gt.25 .and. ktau.ne.0) &
1815 write(6,552)inest,i,j,kpblt(i)
1816 552 FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
1817 ENDIF
1818
1819 if(kpblt(i).gt.25) komax=3
1820 komin=1
1821 dk=float(komax)
1822
1823 do k=komin,komax
1824
1825 wtsig(k)=float(komax-k+1)/dk
1826 WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ
1827
1828 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ*WTSIG(K) &
1829 *WTSIG(K)*ERFIVR
1830
1831 ! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
1832 !
1833 ! write(6,'(a,i2,a,f8.3,a,f8.3,a,f8.3,a,f8.3,a,f8.3)')
1834 ! 'Surface obs, after: k = ',k, &
1835 ! ' WT2ERR = ',wt2err(i,k), &
1836 ! ' TIMEWT = ',timewt, &
1837 ! ' WTIJ = ',wtij, &
1838 ! ' WSIG = ',wtsig(k), &
1839 ! ' ERFIVR = ',erfivr
1840 ! endif
1841
1842 enddo
1843
1844 ENDDO
1845
1846 ! print *, " Surface "
1847
1848 endif ! end check for obs in domain
1849 ! END SURFACE-LAYER U OR V OBS NUDGING
1850
1851 ELSE
1852 ! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
1853 !
1854 ! print *,'in upper air section'
1855 ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
1856 ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
1857 ! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
1858 ! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
1859 !ajb SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
1860
1861 slope = (fdob%RINFMN-fdob%RINFMX)/(psurf(n)-fdob%PFREE)
1862
1863 RINFAC=SLOPE*POB+fdob%RINFMX-SLOPE*fdob%pfree
1864 RINFAC=AMAX1(RINFAC,fdob%RINFMN)
1865 RINFAC=AMIN1(RINFAC,fdob%RINFMX)
1866 !yliu: for multilevel upper-air data, take the maximum
1867 ! for the I loop.
1868 if(nsndlev.gt.1) RINFAC = fdob%RINFMX
1869 !yliu end
1870
1871 MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
1872 MAXI=MIN0(ide-IGRID,MAXI)
1873 MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
1874 MINI=MAX0(1,MINI)
1875
1876 ! yliu start
1877 ! use also obs outside of but close to this domain -- upr data
1878 ! if( RA(N).LT.(0.-RINFAC*RINDX)
1879 ! & .or. RA(N).GT.float(IL)+RINFAC*RINDX
1880 ! & .or. RB(N).LT.(0.-RINFAC*RINDX)
1881 ! & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then
1882 ! print *, " skipped obs far away from this I-range"
1883 ! currently can use obs within this domain or ones very close to (1/3
1884 ! influence of radius in the coarse domain) this
1885 ! domain. In later case, use BC column value to approximate the model value
1886 ! at obs point -- ERRF need model field in errob.F !!
1887
1888 !cc if (i.eq.39 .and. j.eq.34) then
1889 !cc write(6,*) 'RA(N) = ',ra(n)
1890 !cc write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
1891 !cc endif
1892 if( RA(N).GE.(0.-RINFAC*RINDX/3) &
1893 .and. RA(N).LE.float(ide)+RINFAC*RINDX/3 &
1894 .and. RB(N).GE.(0.-RINFAC*RINDX/3) &
1895 .and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
1896 ! or use obs within this domain only
1897 ! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
1898 ! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
1899 ! print *, " skipped obs far outside of this domain"
1900
1901 ! yliu end
1902 ! is this 2 needed here - kpbl not used?
1903 ! MINI=MAX0(2,MINI)
1904
1905 ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
1906 ! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
1907 ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
1908 RJ=FLOAT(J)
1909 RX=RJ-RB(N)
1910 ! WEIGHTS FOR THE 3-D VARIABLES
1911 !
1912 ERFIVR=ERRF(IVAR,N)
1913 ! jc
1914 nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
1915 ! yliu start
1916 ! test: do the sounding levels as individual obs
1917 ! nsndlev=1
1918 ! yliu end
1919 njcsnd=nsndlev
1920 !
1921 DO I=max0(its,MINI),min0(ite,MAXI)
1922 ! jc
1923 RI=FLOAT(I)
1924 RY=RI-RA(N)
1925 RIS=RINDX*RINFAC*RINDX*RINFAC
1926 RSQ=RX*RX+RY*RY
1927 ! yliu test: for upper-air data, keep D1 influence radii
1928 ! RIS=RIS /schnes /schnes
1929 WTIJ=(RIS-RSQ)/(RIS+RSQ)
1930 WTIJ=AMAX1(0.0,WTIJ)
1931 ! weight ob in vertical with +- 50 mb
1932 ! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
1933 if(nsndlev.eq.1) then
1934 rinprs=7.5
1935 else
1936 rinprs=3.0
1937 endif
1938 ! yliu end
1939 !
1940 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
1941 ! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
1942 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
1943
1944 if(nsndlev.eq.1)then
1945 !----------------------------------------------------------------------
1946 ! --- HANDLE 1-LEVEL OBSERVATIONS ---
1947 !----------------------------------------------------------------------
1948
1949 ! if(I.eq.MINI) print *, " Single snd "
1950 ! ERFIVR is the residual (difference) between the ob and the model
1951 ! at that point. We can analyze that residual up and down.
1952 ! First find komin for ob.
1953 !yliu start -- in the old code, komax and komin were reversed!
1954 do k=kte,1,-1
1955 pijk = .001*(pbase(i,k)+pp(i,k))
1956 ! print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
1957 if(pijk.ge.(pob+rinprs)) then
1958 komin=k
1959 go to 325
1960 endif
1961 enddo
1962 komin=1
1963 325 continue
1964 ! now find komax for ob
1965 do k=3,kte
1966 pijk = .001*(pbase(i,k)+pp(i,k))
1967 if(pijk.le.(pob-rinprs)) then
1968 komax=k
1969 go to 326
1970 endif
1971 enddo
1972 komax=kte ! yliu 20050706
1973 326 continue
1974
1975 ! yliu: single-level upper-air data will act either above or below the PBL top
1976 ! komax=min0(kpblt(i), komax)
1977 if(komax.gt.kpblt(i)) komin=max0(kpblt(i), komin)
1978 if(komin.lt.kpblt(i)) komax=min0(kpblt(i), komax)
1979 ! yliu end
1980 !
1981 ! print *,'1 level, komin,komax=',komin,komax
1982 ! if(i.eq.MINI) then
1983 ! print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
1984 ! ERFIVR=0
1985 ! endif
1986 do k=1,kte
1987 reserf(k)=0.0
1988 wtsig(k)=0.0
1989 enddo
1990 !yliu end
1991
1992 !cc if (i.eq.39 .and. j.eq.34) then
1993 !cc write(6,*) ' komin = ', komin,' komax = ',komax
1994 !cc endif
1995
1996 do k=komin,komax
1997 pijk = .001*(pbase(i,k)+pp(i,k))
1998 reserf(k)=erfivr
1999 wtsig(k)=1.-abs(pijk-pob)/rinprs
2000 wtsig(k)=amax1(wtsig(k),0.0)
2001 ! print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
2002 ! Now calculate WT and WT2ERR for each i,j,k point cajb
2003 WT(I,K)=WT(I,K)+TIMEWT*WTIJ*wtsig(k)
2004
2005 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ* &
2006 reserf(k)*wtsig(k)*wtsig(k)
2007 enddo
2008
2009 else
2010 !----------------------------------------------------------------------
2011 ! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
2012 !----------------------------------------------------------------------
2013 !yliu start
2014 ! if(I.eq.MINI) print *, " Multi-level snd "
2015 ! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
2016 ! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
2017 if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
2018 IF (iprt) THEN
2019 print *, "n = ",n,"nsndlev = ",nsndlev
2020 print *, "nlevs_ob,lev_in_ob", &
2021 nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
2022 print *, "in nudobs.F: sounding level messed up, stopping"
2023 ENDIF
2024 stop
2025 endif
2026 !yliu end
2027 ! This is for a multi-level observation
2028 ! The trick here is that the sounding is "one ob". You don't
2029 ! want multiple levels to each be treated like separate
2030 ! and independent observations.
2031 ! At each i,j want to interpolate sounding to the model levels at that
2032 ! particular point.
2033 komin=1
2034 komax=kte-2
2035 ! this loop goes to 1501
2036 ! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
2037 !yliu start
2038 do k=1,kte
2039 reserf(k)=0.0
2040 wtsig(k)=0.0
2041 enddo
2042 !yliu end
2043
2044 do k=komax,komin,-1
2045
2046 pijk = .001*(pbase(i,k)+pp(i,k))
2047
2048 ! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
2049 if(pijk.gt.varobs(5,n)) then
2050 go to 1501
2051 endif
2052
2053 ! if sigma level pressure is .lt. than the highest ob level, don't interpolate
2054 if(pijk.le.varobs(5,n+nsndlev-1)) then
2055 go to 1501
2056 endif
2057
2058 ! now interpolate sounding to this k
2059 ! yliu start-- recalculate WTij for each k-level
2060 !ajb SLOPE = (fdob%RINFMN-fdob%RINFMX)/(pdoc(i,j)+PTOP-fdob%PFREE)
2061 slope = (fdob%RINFMN-fdob%RINFMX)/ (.001*pbase(i,1)-fdob%PFREE)
2062 RINFAC=SLOPE*pijk+fdob%RINFMX-SLOPE*fdob%PFREE
2063 RINFAC=AMAX1(RINFAC,fdob%RINFMN)
2064 RINFAC=AMIN1(RINFAC,fdob%RINFMX)
2065 RIS=RINDX*RINFAC*RINDX*RINFAC
2066 RSQ=RX*RX+RY*RY
2067
2068 ! for upper-air data, keep D1 influence radii
2069 ! RIS=RIS /schnes /schnes
2070 WTIJ=(RIS-RSQ)/(RIS+RSQ)
2071 WTIJ=AMAX1(0.0,WTIJ)
2072 ! yliu end
2073
2074 ! this loop goes to 1503
2075 do nn=2,nsndlev
2076 ! only set pobhi if varobs(ivar) is ok
2077 pobhi=-888888.
2078
2079 if(varobs(ivar,n+nn-1).gt.-800000. &
2080 .and. varobs(5,n+nn-1).gt.-800000.) then
2081 pobhi=varobs(5,n+nn-1)
2082 nhi=n+nn-1
2083 if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.20.) then
2084 go to 1502 ! within 200mb of obs height
2085 endif
2086 endif
2087
2088 enddo
2089
2090 ! did not find any ob above within 100 mb, so jump out
2091 go to 1501
2092 1502 continue
2093
2094 nlo=nhi-1
2095 do nnjc=nhi-1,n,-1
2096 if(varobs(ivar,nnjc).gt.-800000. &
2097 .and. varobs(5,nnjc).gt.-800000.) then
2098 poblo=varobs(5,nnjc)
2099 nlo=nnjc
2100 if(poblo.gt.pijk .and. abs(poblo-pijk).lt.20.) then
2101 go to 1505 ! within 200mb of obs height
2102 endif
2103 endif
2104 enddo
2105 !yliu end --
2106
2107 ! did not find any ob below within 200 mb, so jump out
2108 go to 1501
2109 1505 continue
2110
2111 ! interpolate to model level
2112 pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
2113 reserf(k)=errf(ivar,nlo)+ &
2114 (errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
2115 wtsig(k)=1.
2116
2117 1501 continue
2118
2119 ! now calculate WT and WT2ERR for each i,j,k point cajb
2120 WT(I,K)=WT(I,K)+TIMEWT*WTIJ*wtsig(k)
2121
2122 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ* &
2123 reserf(k)*wtsig(k)*wtsig(k)
2124
2125 ! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
2126 !
2127 ! if(wt(i,k) .ne. 0.0) then
2128 ! scratch = WT2ERR(I,K)/WT(I,K)
2129 ! else
2130 ! scratch = 999.
2131 ! endif
2132 !
2133 ! write(6,'(a,i2,a,f8.3,a,f4.2,a,f7.4,a,f4.2,a,f5.3,a,f7.4)')
2134 ! $ 'Multi-level obs: k = ',k,
2135 ! $ ' WT2ERR = ',wt2err(i,k),
2136 ! $ ' WTIJ = ',wtij,
2137 ! $ ' RSF = ',reserf(k),
2138 ! $ ' WSIG = ',wtsig(k),
2139 ! $ ' WT = ',wt(i,k),
2140 ! $ ' W2EOWT = ',scratch
2141 ! endif
2142
2143
2144 ! end do k
2145 enddo ! enddo k levels
2146 ! end multi-levels
2147 endif ! end if(nsndlev.eq.1)
2148 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2149 ! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
2150 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2151 !
2152 ENDDO ! END DO MINI,MAXI LOOP
2153
2154 endif ! check for obs in domain
2155
2156 ! END OF NUDGING TO OBS ON PRESSURE LEVELS
2157
2158 ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
2159
2160 !----------------------------------------------------------------------
2161 ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
2162 !----------------------------------------------------------------------
2163
2164 ! n=n+1
2165 n=n+njcsnd
2166
2167 !yliu 1202 continue
2168 if(n.gt.nstat)then
2169 ! print *,'n,nstat=',n,nstat,ivar,j
2170 go to 1203
2171 endif
2172 ! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
2173
2174 !***********************************************************************
2175 ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
2176 !***********************************************************************
2177
2178 1203 continue
2179
2180 ! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
2181 ! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
2182 ! THE ATEN ARRAY
2183 ! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
2184 ! THEY ARE USED BELOW IN THE DENOMINATOR.
2185 DO K=kts,kte
2186 DO I=its,ite
2187 IF(WT(I,K).EQ.0)THEN
2188 WT2ERR(I,K)=0.0
2189 ENDIF
2190 IF(WT(I,K).EQ.0)THEN
2191 WT(I,K)=1.0
2192 ENDIF
2193 ENDDO
2194 ENDDO
2195
2196 126 CONTINUE
2197
2198 IF(IVAR.GE.3)GOTO 170
2199 ! this is for u,v
2200 ! 3-D DOT POINT TENDENCIES
2201
2202 ! Calculate scales for converting nudge factor from u (v)
2203 ! to rho_u (or rho_v) units.
2204
2205 call calc_rcouple_scales(mu,msf,rscale,ims,ime,its,ite)
2206
2207 DO K=1,kte
2208
2209 DO I=i_s,i_e
2210
2211 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2212 W2EOWT=WT2ERR(I,K)/WT(I,K)
2213 ELSE
2214 W2EOWT=SAVWT(IPL,I,K)
2215 ENDIF
2216
2217 ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2218 ! scratch = GIV*RSCALE(I)*W2EOWT*fdob%TFACI*ISWIND*GFACTOR
2219 ! write(6,*) 'ATEN calc: k = ',k
2220 ! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
2221 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2222 ! $ ' W2EOWT = ',w2eowt
2223 ! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
2224 ! $ ' GFACTOR = ',gfactor
2225 ! endif
2226 !
2227 ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2228 ! scratch = GIV*RSCALE(I)*W2EOWT*fdob%TFACI*ISWIND*GFACTOR
2229 ! write(6,*) 'ATEN calc: k = ',k
2230 ! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
2231 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2232 ! $ ' W2EOWT = ',w2eowt
2233 ! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
2234 ! $ ' GFACTOR = ',gfactor
2235 ! endif
2236
2237 ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
2238 *W2EOWT*fdob%TFACI &
2239 *ISWIND *GFACTOR !yliu *GFACTOR
2240
2241 ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2242 ! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
2243 ! endif
2244 ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2245 ! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
2246 ! endif
2247
2248 ENDDO
2249 ENDDO
2250
2251 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2252 DO K=1,kte
2253 DO I=its,ite
2254 SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2255 ENDDO
2256 ENDDO
2257 ENDIF
2258
2259 RETURN
2260
2261 170 CONTINUE
2262
2263 ! 3-D CROSS-POINT TENDENCIES
2264 ! this is for t (ivar=3) and q (ivsr=4)
2265 IF(3-IVAR.LT.0)THEN
2266 GITQ=GIQ
2267 ELSE
2268 GITQ=GIT
2269 ENDIF
2270 IF(3-IVAR.LT.0)THEN
2271 ISTQ=ISMOIS
2272 ELSE
2273 ISTQ=ISTEMP
2274 ENDIF
2275
2276 DO K=1,kte
2277 DO I=i_s,i_e
2278 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2279 W2EOWT=WT2ERR(I,K)/WT(I,K)
2280 ELSE
2281 W2EOWT=SAVWT(IPL,I,K)
2282 ENDIF
2283
2284 ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
2285 ! scratch = GITQ*MU(I)*W2EOWT*fdob%TFACI*ISTQ*GFACTOR
2286 ! write(6,*) 'ATEN calc: k = ',k
2287 ! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
2288 ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
2289 ! $ ' W2EOWT = ',w2eowt
2290 ! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
2291 ! $ ' GFACTOR = ',gfactor
2292 ! endif
2293 !
2294 ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2295 ! scratch = GITQ*MU(I)*W2EOWT*fdob%TFACI*ISTQ*GFACTOR
2296 ! write(6,*) 'ATEN calc: k = ',k
2297 ! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
2298 ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
2299 ! $ ' W2EOWT = ',w2eowt
2300 ! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
2301 ! $ ' GFACTOR = ',gfactor
2302 ! endif
2303
2304 ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
2305 *W2EOWT*fdob%TFACI*ISTQ *GFACTOR !yliu *GFACTOR
2306
2307 ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
2308 ! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
2309 ! endif
2310 ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2311 ! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
2312 ! endif
2313
2314 ENDDO
2315 ENDDO
2316
2317 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
2318 DO K=1,kte
2319 DO I=its,ite
2320 SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2321 ENDDO
2322 ENDDO
2323 ENDIF
2324
2325 RETURN
2326 END SUBROUTINE nudob
2327
2328 SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
2329 !-----------------------------------------------------------------------
2330 IMPLICIT NONE
2331 !-----------------------------------------------------------------------
2332
2333 INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
2334 INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
2335 REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
2336 REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
2337 REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
2338
2339 ! Local variables
2340 integer :: i
2341
2342 ! Calculate scales to be used for producing rho-coupled nudging factors.
2343 do i = its,ite
2344 rscale(i) = a(i)/msf(i)
2345 enddo
2346
2347 RETURN
2348 END SUBROUTINE calc_rcouple_scales
2349
2350 !ajb: Not used
2351 SUBROUTINE set_real_array(rscale, value, ims,ime, its,ite)
2352 !-----------------------------------------------------------------------
2353 IMPLICIT NONE
2354 !-----------------------------------------------------------------------
2355
2356 INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
2357 INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
2358 REAL, INTENT(IN) :: value ! Constant array value
2359 REAL, INTENT(OUT) :: rscale( ims:ime ) ! Output array
2360
2361 ! Local variables
2362 integer :: i
2363
2364 ! Set array to constant value
2365 do i = its,ite
2366 rscale(i) = value
2367 enddo
2368
2369 RETURN
2370 END SUBROUTINE set_real_array
2371
2372 !ajb: Not used
2373 SUBROUTINE calc_pottemp_scales(ivar, rcp, pb, p, tscale, &
2374 ims,ime, its,ite, &
2375 kms,kme, kts,kte)
2376 !-----------------------------------------------------------------------
2377 IMPLICIT NONE
2378 !-----------------------------------------------------------------------
2379
2380 INTEGER, INTENT(IN) :: ims,ime, kms,kme ! Memory dimensions
2381 INTEGER, INTENT(IN) :: its,ite, kts,kte ! Tile dimensions
2382 INTEGER, INTENT(IN) :: ivar ! Variable identifier
2383 REAL, INTENT(IN) :: rcp ! Constant (2./7.)
2384 REAL, INTENT(IN) :: pb(ims:ime, kms:kme) ! Base pressure (Pa) array
2385 REAL, INTENT(IN) :: p(ims:ime, kms:kme) ! Pressure pert. (Pa) array
2386 REAL, INTENT(OUT) :: tscale(ims:ime, kms:kme) ! Scales for pot. temp.
2387 ! Local variables
2388 integer :: i,k
2389
2390 if(ivar.eq.3) then
2391
2392 ! Calculate scales to be used for producing potential temperature nudging factors.
2393 do k = kts,kte
2394 do i = its,ite
2395 tscale(i,k) = ( 1000000. / ( pb(i,k)+p(i,k)) )**rcp
2396 enddo
2397 enddo
2398 else
2399 ! Set to 1. for moisture scaling.
2400 do k = kts,kte
2401 do i = its,ite
2402 tscale(i,k) = 1.0
2403 enddo
2404 enddo
2405 endif
2406
2407 RETURN
2408 END SUBROUTINE calc_pottemp_scales
2409 #endif
2410
2411 END MODULE module_fddaobs_rtfdda
2412