gsi_thinning.f90
References to this file elsewhere.
1 module gsi_thinning
2 !$$$ subprogram documenation block
3 !
4 ! abstract: This module contains code which may be used to selectively
5 ! thin satellite data.
6 !
7 ! program history log:
8 ! 2006-10-28 Jianjun Xu, developed from GSI
9 ! 2007-03-30 Zhiquan Liu, modify and add comments
10 !
11 !
12 ! Subroutines Included:
13 ! sub makegvals - set up for superob weighting
14 ! sub makegrids - set up thinning grids
15 ! sub map2tgrids - map observation to location on thinning grid
16 ! sub destroygrids - deallocate thinning grid arrays
17 ! sub destroy_sfc - deallocate full horizontal surface arrays
18 !
19 ! Variable Definitions:
20 ! def mlat - number of latitudes in thinning grid
21 ! def mlon - number of longitudes in thinning grid
22 ! def maxthin - maximum number of obs to retain in thinning grid box
23 ! def itxmax - total number of points in thinning grid
24 ! def istart_val - starting location on thinning grid for superobs (not used)
25 ! def icount - observation counter
26 ! def ibest_obs - index of "best" quality obs in thinning grid box
27 ! def glat - latitudes on thinning grid
28 ! def glon - longitudes on thinning grid
29 ! def hll - (i,j) index of point on thinning grid
30 ! def score_crit - "best" quality obs score in thinning grid box
31 !
32 !$$$ end documentation block
33 use gsi_kinds, only: r_kind,i_kind
34 use gsi_constants, only: deg2rad,rearth_equator,zero,two,pi,half,one,quarter,&
35 rad2deg
36 implicit none
37
38 real(r_kind),parameter:: r90 = 90.0_r_kind
39 real(r_kind),parameter:: r360 = 360.0_r_kind
40 real(r_kind),parameter:: r999 = 999.0_r_kind
41 real(r_kind),parameter:: r1000 = 1000.0_r_kind
42
43 ! lat/lon range inside tile
44 real(r_kind) rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid
45
46 type thinning_type
47 ! mlat: lat #, mlonx: max lon #, itxmax: grid box #
48 integer(i_kind) mlat,maxthin,itxmax,dthin,mlonx,mlony
49 integer(i_kind),dimension(0:51):: istart_val
50
51 ! mlon(mlat): lon # in each lat
52 integer(i_kind),allocatable,dimension(:):: mlon,icount,ibest_obs
53 integer(i_kind),allocatable,dimension(:,:):: isli
54
55 ! glat(mlat): lat #, glon(mlat,mlonx), hll(mlat,mlonx)
56 integer(i_kind),allocatable,dimension(:,:) :: hll
57 real(r_kind),allocatable,dimension(:) :: glat
58 real(r_kind),allocatable,dimension(:,:) :: glon,sli,sno
59 real(r_kind),allocatable,dimension(:) :: score_crit
60 end type thinning_type
61
62 type(thinning_type), allocatable :: thinning_grid(:)
63
64 contains
65
66 subroutine makegrids (n,rmesh)
67 ! compute dimention of thinning box
68 ! output (mlat,mlonx,istart_val)
69 implicit none
70
71 integer(i_kind), intent(in) :: n ! sensor index
72 real(r_kind), intent(in) :: rmesh ! thinning box size
73
74 logical odd
75 integer(i_kind) i,ii,j,k,nlat,nlon
76 integer(i_kind) icnt,mlonj
77 real(r_kind) delonx,delat,dgv,dx,dy
78 real(r_kind) twopi,halfpi,dlon_g,dlat_g,dlon_e,dlat_e
79 real(r_kind) factor,factors,delon
80 real(r_kind) rkm2dg,glatm,glatx
81
82 ! Initialize variables, set constants
83 thinning_grid(n)%dthin = 1
84 thinning_grid(n)%maxthin=thinning_grid(n)%dthin
85
86 thinning_grid(n)%istart_val=0
87 twopi = two*pi
88 halfpi = pi*half
89 rkm2dg = r360/(twopi*rearth_equator)*r1000
90
91 dx = rmesh*rkm2dg
92 dy = dx
93 thinning_grid(n)%mlat = dlat_grid/dy + half
94 thinning_grid(n)%mlonx = dlon_grid/dx + half
95 delat = dlat_grid/thinning_grid(n)%mlat
96 delonx= dlon_grid/thinning_grid(n)%mlonx
97 dgv = delat*half
98
99 thinning_grid(n)%mlat=max(2,thinning_grid(n)%mlat)
100 thinning_grid(n)%mlonx=max(2,thinning_grid(n)%mlonx)
101
102 do ii=1,thinning_grid(n)%maxthin
103 thinning_grid(n)%istart_val(ii+1)=thinning_grid(n)%istart_val(ii)
104 icnt=0
105 do j = 1,thinning_grid(n)%mlat
106 glatx = rlat_min + (j-1)*delat
107 glatx = glatx*deg2rad
108 glatm = glatx + dgv*deg2rad
109 factor = abs(cos(abs(glatm)))
110 mlonj = nint(thinning_grid(n)%mlonx*factor)
111 mlonj = max(2,mlonj)
112 do i = 1,mlonj
113 icnt=icnt+1
114 thinning_grid(n)%istart_val(ii+1)=thinning_grid(n)%istart_val(ii+1)+1
115 enddo
116 enddo
117 end do
118
119 ! making thinning box
120 ! output: mlon(mlat),glat(mlat),glon(mlonx,mlat),hll(mlonx,mlat)
121
122 allocate(thinning_grid(n)%mlon(thinning_grid(n)%mlat), &
123 thinning_grid(n)%glat(thinning_grid(n)%mlat), &
124 thinning_grid(n)%glon(thinning_grid(n)%mlonx,thinning_grid(n)%mlat), &
125 thinning_grid(n)%hll(thinning_grid(n)%mlonx,thinning_grid(n)%mlat))
126
127 ! Set up thinning grid lon & lat. The lon & lat represent the location of the
128 ! lower left corner of the thinning grid box.
129
130 thinning_grid(n)%itxmax=0
131 do j = 1,thinning_grid(n)%mlat
132 thinning_grid(n)%glat(j) = rlat_min + (j-1)*delat
133 thinning_grid(n)%glat(j) = thinning_grid(n)%glat(j)*deg2rad
134 glatm = thinning_grid(n)%glat(j) + dgv*deg2rad
135
136 factor = abs(cos(abs(glatm)))
137 mlonj = nint(thinning_grid(n)%mlonx*factor)
138 thinning_grid(n)%mlon(j) = max(2,mlonj)
139 delon = dlon_grid/thinning_grid(n)%mlon(j)
140
141 thinning_grid(n)%glat(j) = min(max(-halfpi,thinning_grid(n)%glat(j)),halfpi)
142 do i = 1,thinning_grid(n)%mlon(j)
143 thinning_grid(n)%itxmax=thinning_grid(n)%itxmax+1
144 thinning_grid(n)%hll(i,j)=thinning_grid(n)%itxmax
145 thinning_grid(n)%glon(i,j) = rlon_min + (i-1)*delon
146 thinning_grid(n)%glon(i,j) = thinning_grid(n)%glon(i,j)*deg2rad
147 thinning_grid(n)%glon(i,j) = min(max(zero,thinning_grid(n)%glon(i,j)),twopi)
148 enddo
149 !write(6,'(f10.5,i8,2i10)') glat(j)*rad2deg, mlon(j),hll(1,j),hll(mlon(j),j)
150 !write(6,'(10f8.3)') (glon(i,j)*rad2deg,i=1,mlon(j))
151
152 end do
153
154 ! Allocate and initialize arrays
155 allocate(thinning_grid(n)%icount(thinning_grid(n)%itxmax))
156 allocate(thinning_grid(n)%ibest_obs(thinning_grid(n)%itxmax))
157 allocate(thinning_grid(n)%score_crit(thinning_grid(n)%itxmax))
158
159 do j=1,thinning_grid(n)%itxmax
160 thinning_grid(n)%icount(j) = 0
161 thinning_grid(n)%ibest_obs(j) = 0
162 thinning_grid(n)%score_crit(j) = 9.99e6_r_kind
163 end do
164
165 return
166 end subroutine makegrids
167
168 subroutine map2grids(n,dlat_earth,dlon_earth,crit1,iobs,itx,ithin,itt,iobsout,iuse)
169 !$$$ subprogram documentation block
170 ! . . . .
171 ! subprogram: map2grids
172 ! prgmmr: treadon org: np23 date: 2002-10-17
173 !
174 ! abstract: This routine maps observations to the thinning grid.
175 !
176 ! program history log:
177 ! 2002-10-17 treadon
178 ! 2004-06-22 treadon - update documentation
179 ! 2004-07-23 derber - modify code to thin obs as read in
180 ! 2004-12-08 li, xu - fix bug --> set iuse=.true. when use_all=.true.
181 ! 2005-10-14 treadon - variable name change (dlat0,dlon0) --> d*_earth
182 ! 2006-03-25 kistler - define iobsout for the case use_all=.true.
183 !
184 ! input argument list:
185 ! n - sensor index
186 ! dlat_earth - earth relative observation latitude (radians)
187 ! dlon_earth - earth relative observation longitude (radians)
188 ! crit1 - quality indicator for observation (smaller = better)
189 ! ithin - number of obs to retain per thinning grid box
190 !
191 ! output argument list:
192 ! iobs - observation counter
193 ! itx - combined (i,j) index of observation on thinning grid
194 ! itt - superobs thinning counter
195 ! iobsout- location for observation to be put
196 ! iuse - .true. if observation should be used
197 !
198 implicit none
199 logical, intent(out) :: iuse
200 integer(i_kind), intent(out) :: itt,itx
201 integer(i_kind), intent(in) :: ithin,n
202 integer(i_kind), intent(inout) :: iobs,iobsout
203 real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1
204
205 integer(i_kind) :: ix,iy
206 real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy
207 real(r_kind) dist1,crit
208
209 ! Compute (i,j) indices of coarse mesh grid (grid number 1) which
210 ! contains the current observation.
211 dlat1=dlat_earth
212 dlon1=dlon_earth
213
214 call grdcrd(dlat1,1,thinning_grid(n)%glat,thinning_grid(n)%mlat,1)
215 iy=int(dlat1)
216 dy=dlat1-iy
217 iy=max(1,min(iy,thinning_grid(n)%mlat))
218
219 call grdcrd(dlon1,1,thinning_grid(n)%glon(1,iy),thinning_grid(n)%mlon(iy),1)
220 ix=int(dlon1)
221 dx=dlon1-ix
222 ix=max(1,min(ix,thinning_grid(n)%mlon(iy)))
223
224 dxx=half-min(dx,one-dx)
225 dyy=half-min(dy,one-dy)
226 dist1=dxx*dxx+dyy*dyy+half
227 itx=thinning_grid(n)%hll(ix,iy)
228 itt=thinning_grid(n)%istart_val(ithin)+itx
229 if(ithin == 0) itt=0
230
231 ! Increment obs counter on coarse mesh grid. Also accumulate observation
232 ! score and distance functions
233
234 thinning_grid(n)%icount(itx)=thinning_grid(n)%icount(itx)+1
235 ! dist1=one - quarter*(dista + distb) !dist1 is min at grid box center and
236 !ranges from 1 (at corners)to
237 !.5 (at center of box)
238 crit=crit1*dist1
239 iuse=.false.
240
241 if(thinning_grid(n)%icount(itx) == 1)then
242
243 ! Increment obs counter
244
245 iuse=.true.
246 iobs=iobs+1
247 thinning_grid(n)%score_crit(itx)= crit
248 thinning_grid(n)%ibest_obs(itx) = iobs
249 iobsout=iobs
250
251 end if
252 if(crit < thinning_grid(n)%score_crit(itx) .and. thinning_grid(n)%icount(itx) > 1)then
253 iuse=.true.
254 thinning_grid(n)%score_crit(itx)= crit
255 iobsout=thinning_grid(n)%ibest_obs(itx)
256 end if
257
258 return
259 end subroutine map2grids
260
261 subroutine map2tgrid(n,dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse)
262 ! input argument list:
263 ! dlat_earth - earth relative observation latitude (radians)
264 ! dlon_earth - earth relative observation longitude (radians)
265 ! crit1 - quality indicator for observation (smaller = better)
266 ! ithin - number of obs to retain per thinning grid box
267 !
268 ! output argument list:
269 ! iobs - observation counter
270 ! itx - combined (i,j) index of observation on thinning grid
271 ! itt - superobs thinning counter
272 ! iobsout- location for observation to be put
273 ! iuse - .true. if observation should be used
274 !
275 implicit none
276 logical,intent(out):: iuse
277 integer(i_kind),intent(in):: ithin,n
278 integer(i_kind),intent(out):: itt,itx
279 real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1
280 real(r_kind),intent(out):: dist1
281
282 integer(i_kind) ix,iy
283 real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy
284
285
286 ! Compute (i,j) indices of coarse mesh grid (grid number 1) which
287 ! contains the current observation.
288 dlat1=dlat_earth
289 dlon1=dlon_earth
290
291 call grdcrd(dlat1,1,thinning_grid(n)%glat,thinning_grid(n)%mlat,1)
292 iy=int(dlat1)
293 dy=dlat1-iy
294 iy=max(1,min(iy,thinning_grid(n)%mlat))
295
296 call grdcrd(dlon1,1,thinning_grid(n)%glon(1,iy),thinning_grid(n)%mlon(iy),1)
297 ix=int(dlon1)
298 dx=dlon1-ix
299 ix=max(1,min(ix,thinning_grid(n)%mlon(iy)))
300
301 dxx=half-min(dx,one-dx)
302 dyy=half-min(dy,one-dy)
303 dist1=dxx*dxx+dyy*dyy+half
304 itx=thinning_grid(n)%hll(ix,iy)
305 itt=thinning_grid(n)%istart_val(ithin)+itx
306 if(ithin == 0) itt=0
307 iuse=.true.
308 if(dist1*crit1 > thinning_grid(n)%score_crit(itx) .and. thinning_grid(n)%icount(itx) == 0)iuse=.false.
309
310 !write(6,'(a,3f10.3)') 'dlat_earth dlon_earth crit1 ',dlat_earth*rad2deg,dlon_earth*rad2deg,crit1
311 !write(6,'(a,2i5,3f10.3,i10,e12.5,2x,L)') 'ix iy',ix,iy,dx,dy,dist1,itx,score_crit(itx),iuse
312 return
313 end subroutine map2tgrid
314
315 subroutine grdcrd(d,nd,x,nx,flg)
316 implicit none
317 integer(i_kind),intent(in) :: nd,nx,flg
318 real(r_kind),intent(inout) :: d
319 real(r_kind),dimension(nx), intent(in) :: x
320
321 integer(i_kind) :: id,ix
322 ! Treat "normal" case in which nx>1
323 if(nx>1) then
324 if (flg.eq.1) then
325
326 ! Case in which x is in increasing order
327 if(d<=x(1)) then
328 ix=1
329 else
330 ix=isrchf(nx-1,x,d,flg)-1
331 end if
332 if(ix==nx) ix=ix-1
333
334 else if (flg.eq.(-1)) then
335
336 ! Case in which x is in decreasing order
337 if(d>=x(1)) then
338 ix=1
339 else
340 ix=isrchf(nx-1,x,d,flg)-1
341 end if
342 end if
343 d=float(ix)+(d-x(ix))/(x(ix+1)-x(ix))
344
345 ! Treat special case of nx=1
346 elseif (nx==1) then
347 d = one
348 endif
349
350 return
351 end subroutine grdcrd
352
353 subroutine checkob(n,dist1,crit1,itx,iuse)
354 !
355 ! input argument list:
356 ! dist1 - quality indicator for distance (smaller = better)
357 ! crit1 - quality indicator for observation (smaller = better)
358 ! itx - combined (i,j) index of observation on thinning grid
359 ! iuse - .true. if observation should be used
360 !
361 ! output argument list:
362 ! iuse - .true. if observation should be used
363 !
364 implicit none
365 logical,intent(inout):: iuse
366 integer(i_kind),intent(in):: n,itx
367 real(r_kind),intent(in):: dist1,crit1
368
369 ! If data (no thinning), simply return to calling routine
370 if(.not. iuse .or. thinning_grid(n)%icount(itx)==0)return
371 if(crit1*dist1 > thinning_grid(n)%score_crit(itx))iuse=.false.
372
373 return
374 end subroutine checkob
375
376 subroutine finalcheck(n,dist1,crit1,iobs,itx,iobsout,iuse,sis)
377 !
378 ! input argument list:
379 ! dist1 - quality indicator for distance (smaller = better)
380 ! crit1 - quality indicator for observation (smaller = better)
381 ! itx - combined (i,j) index of observation on thinning grid
382 ! iobs - observation counter
383 ! iuse - .true. if observation should be used
384 ! sis - sensor/instrument/satellite
385 !
386 ! output argument list:
387 ! iobs - observation counter
388 ! iobsout- location for observation to be put
389 ! iuse - .true. if observation should be used
390 !
391 implicit none
392 logical,intent(inout):: iuse
393 integer(i_kind),intent(inout):: iobs,iobsout
394 integer(i_kind),intent(in):: n,itx
395 real(r_kind),intent(in):: dist1,crit1
396 character(20),intent(in):: sis
397
398 real(r_kind) crit
399
400 if(.not. iuse)return
401
402 ! If using all data (no thinning), simply return to calling routine
403
404
405 crit=crit1*dist1
406
407 if(thinning_grid(n)%icount(itx) == 0)then
408
409 ! Increment obs counter
410
411 if(iobs < thinning_grid(n)%itxmax)then
412 iobs=iobs+1
413 iobsout=iobs
414 thinning_grid(n)%score_crit(itx)= crit
415 thinning_grid(n)%ibest_obs(itx) = iobs
416 thinning_grid(n)%icount(itx)=thinning_grid(n)%icount(itx)+1
417 else
418 iuse = .false.
419 write(6,*)' ndata > maxobs when reading data for ',sis,thinning_grid(n)%itxmax
420 end if
421
422 else if(crit < thinning_grid(n)%score_crit(itx))then
423 thinning_grid(n)%score_crit(itx)= crit
424 iobsout=thinning_grid(n)%ibest_obs(itx)
425 thinning_grid(n)%icount(itx)=thinning_grid(n)%icount(itx)+1
426 else
427 iuse = .false.
428 end if
429
430
431 return
432 end subroutine finalcheck
433
434 subroutine destroygrids(n)
435 implicit none
436 integer(i_kind), intent(in) :: n
437 deallocate(thinning_grid(n)%mlon,thinning_grid(n)%glat, &
438 thinning_grid(n)%glon,thinning_grid(n)%hll)
439 deallocate(thinning_grid(n)%icount)
440 deallocate(thinning_grid(n)%ibest_obs)
441 deallocate(thinning_grid(n)%score_crit)
442 return
443 end subroutine destroygrids
444
445 subroutine destroy_sfc(n)
446 implicit none
447 integer(i_kind), intent(in) :: n
448 deallocate(thinning_grid(n)%sli,thinning_grid(n)%sno,thinning_grid(n)%isli)
449 return
450 end subroutine destroy_sfc
451
452 function isrchf(nx1,x,y,flg)
453 ! . . . .
454 ! subprogram: isrchf
455 ! prgmmr: parrish org: np22 date: 1990-10-11
456 !
457 ! abstract: get grid coordinates from monotonically increasing or
458 ! decreasing points
459 !
460 ! program history log:
461 ! 2005-03-07 treadon - add doc block
462 !
463 ! input argument list:
464 ! nx1 - number of input points
465 ! x - grid values
466 ! y - target value
467 ! flg - marks order of values in x
468 ! (1=increasing, -1=decreasing)
469 !
470 ! output argument list:
471 ! isrchf - array index of input grid value near target value
472 !
473 implicit none
474 integer(i_kind):: isrchf
475 integer(i_kind),intent(in):: nx1
476 integer(i_kind),intent(in):: flg
477 real(r_kind),intent(in):: y
478 real(r_kind),dimension(nx1),intent(in):: x
479
480 integer(i_kind) k
481
482 if(flg.eq.1) then
483 do k=1,nx1
484 if(y<=x(k)) then
485 isrchf=k
486
487 go to 100
488 end if
489 end do
490 else
491 do k=1,nx1
492 if(y>=x(k)) then
493 isrchf=k
494 go to 100
495 end if
496 end do
497 end if
498
499 isrchf=nx1+1
500 if(nx1<=0) isrchf=0
501
502 100 continue
503 return
504 end function isrchf
505
506 end module gsi_thinning