Line data Source code
1 : module fill_module
2 : !-----------------------------------------------------------------------
3 : ! $Id$
4 : !BOP
5 : !
6 : ! !MODULE: fill_module --- utilities for filling in "bad" data
7 :
8 : use shr_kind_mod, only: r8 => shr_kind_r8
9 : use cam_logfile, only: iulog
10 : use cam_abortutils, only: endrun
11 :
12 :
13 : #ifdef NO_R16
14 : integer,parameter :: r16= selected_real_kind(12) ! 8 byte real
15 : #else
16 : integer,parameter :: r16= selected_real_kind(24) ! 16 byte real
17 : #endif
18 :
19 : !
20 : ! !PUBLIC MEMBER FUNCTIONS:
21 : public filew, fillxy, fillz, filns, pfix, fill_readnl
22 :
23 : !
24 : ! !DESCRIPTION:
25 : !
26 : ! This module provides the basic utilities to fill in regions
27 : ! with bad "data", for example slightly negative values in fields
28 : ! which must be positive, like mixing ratios. Generally this
29 : ! means borrowing positive values from neighboring cells.
30 : !
31 : ! !REVISION HISTORY:
32 : ! 99.03.01 Lin Creation
33 : ! 01.02.14 Lin Routines coalesced into this module
34 : ! 01.03.26 Sawyer Added ProTeX documentation
35 : ! 05.05.25 Sawyer Merged CAM and GEOS5 versions
36 : !
37 : !EOP
38 : !-----------------------------------------------------------------------
39 :
40 : private
41 :
42 : real(r8), parameter :: D0_0 = 0.0_r8
43 : real(r8), parameter :: D0_5 = 0.5_r8
44 : real(r8), parameter :: D1_0 = 1.0_r8
45 : real(r8), parameter :: D1_5 = 1.5_r8
46 :
47 : character(len=8) :: print_filew_warn = 'off'
48 :
49 : contains
50 :
51 1536 : subroutine fill_readnl(nlfile)
52 : use namelist_utils, only: find_group_name
53 : use units, only: getunit, freeunit
54 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_character, masterproc
55 : ! File containing namelist input.
56 : character(len=*), intent(in) :: nlfile
57 :
58 : ! Local variables
59 : integer :: unitn, ierr
60 : character(len=*), parameter :: sub = 'fill_readnl'
61 :
62 : namelist /fill_nl/ print_filew_warn
63 :
64 1536 : if (masterproc) then
65 2 : unitn = getunit()
66 2 : open( unitn, file=trim(nlfile), status='old' )
67 2 : call find_group_name(unitn, 'fill_nl', status=ierr)
68 2 : if (ierr == 0) then
69 0 : read(unitn, fill_nl, iostat=ierr)
70 0 : if (ierr /= 0) then
71 0 : call endrun(sub // ':: ERROR reading namelist fill_nl')
72 : end if
73 : end if
74 2 : close(unitn)
75 2 : call freeunit(unitn)
76 : end if
77 :
78 1536 : call mpi_bcast(print_filew_warn, len(print_filew_warn), mpi_character, mstrid, mpicom, ierr)
79 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: print_filew_warn")
80 :
81 1536 : end subroutine fill_readnl
82 :
83 : !-----------------------------------------------------------------------
84 : !BOP
85 : ! !IROUTINE: filew --- Fill from east and west neighbors; essentially
86 : ! performing local flux adjustment
87 : !
88 : ! !INTERFACE:
89 41287680 : subroutine filew(q, im, jm, jfirst, jlast, acap, ipx, tiny, cosp2)
90 :
91 : ! !USES:
92 :
93 : implicit none
94 :
95 : ! !INPUT PARAMETERS:
96 : integer im ! Longitudes
97 : integer jm ! Total latitudes
98 : integer jfirst ! Starting latitude
99 : integer jlast ! Finishing latitude
100 :
101 : real(r8) tiny ! A small number to pump up value
102 : real(r8) acap ! 1/(polar cap area)
103 : real(r8) cosp2 ! cosine(lat) at j=2
104 :
105 : ! !INPUT/OUTPUT PARAMETERS:
106 : real(r8) q(im,jfirst:jlast) ! Field to adjust
107 :
108 : ! !OUTPUT PARAMETERS:
109 : integer ipx ! Flag: 0 if Q not change, 1 if changed
110 :
111 : ! !DESCRIPTION:
112 : ! Check for "bad" data and fill from east and west neighbors
113 : !
114 : ! !REVISION HISTORY:
115 : ! 01.99.10 Lin Creation
116 : ! 01.07.30 Lin Improvement
117 : !
118 : !EOP
119 : !-----------------------------------------------------------------------
120 : !BOC
121 : ! !LOCAL VARIABLES:
122 : real(r8) d0, d1, d2
123 82575360 : real(r8) qtmp(jfirst:jlast,im)
124 : real(r8) tinyl ! local tiny mixing ratio
125 : real(r8) qmin
126 :
127 : integer i, j, jm1, ip2
128 : integer j1, j2
129 : integer imin, jmin
130 :
131 41287680 : j1 = max( jfirst, 2 )
132 41287680 : j2 = min( jlast, jm-1 )
133 41287680 : jm1 = jm-1
134 41287680 : ipx = 0
135 :
136 : ! Copy & swap direction for vectorization.
137 163860480 : do j=j1,j2
138 35464826880 : do i=1,im
139 35423539200 : qtmp(j,i) = q(i,j)
140 : enddo
141 : enddo
142 :
143 11849564160 : do i=2,im-1
144 46905384960 : do j=j1,j2
145 46864097280 : if(qtmp(j,i) < D0_0) then
146 7416321 : tinyl = max(D0_0,qtmp(j,i-1),qtmp(j,i+1))*tiny
147 7416321 : ipx = 1
148 : ! west
149 7416321 : d0 = max(D0_0,qtmp(j,i-1))
150 7416321 : d1 = min(-qtmp(j,i),d0)
151 7416321 : qtmp(j,i-1) = qtmp(j,i-1) - d1
152 7416321 : qtmp(j,i) = qtmp(j,i) + d1
153 : ! east
154 7416321 : d0 = max(D0_0,qtmp(j,i+1))
155 7416321 : d2 = min(-qtmp(j,i),d0)
156 7416321 : qtmp(j,i+1) = qtmp(j,i+1) - d2
157 7416321 : qtmp(j,i) = qtmp(j,i) + d2 + tinyl
158 : endif
159 : enddo
160 : enddo
161 :
162 163860480 : i=1
163 163860480 : do j=j1,j2
164 163860480 : if(qtmp(j,i) < D0_0) then
165 19397 : ipx = 1
166 19397 : tinyl = max(D0_0,qtmp(j,im),qtmp(j,i+1))*tiny
167 : ! west
168 19397 : d0 = max(D0_0,qtmp(j,im))
169 19397 : d1 = min(-qtmp(j,i),d0)
170 19397 : qtmp(j,im) = qtmp(j,im) - d1
171 19397 : qtmp(j,i) = qtmp(j,i) + d1
172 : ! east
173 19397 : d0 = max(D0_0,qtmp(j,i+1))
174 19397 : d2 = min(-qtmp(j,i),d0)
175 19397 : qtmp(j,i+1) = qtmp(j,i+1) - d2
176 19397 : qtmp(j,i) = qtmp(j,i) + d2 + tinyl
177 : endif
178 : enddo
179 :
180 41287680 : i=im
181 163860480 : do j=j1,j2
182 163860480 : if(qtmp(j,i) < D0_0) then
183 22897 : ipx = 1
184 22897 : tinyl = max(D0_0,qtmp(j,i-1),qtmp(j,1))*tiny
185 : ! west
186 22897 : d0 = max(D0_0,qtmp(j,i-1))
187 22897 : d1 = min(-qtmp(j,i),d0)
188 22897 : qtmp(j,i-1) = qtmp(j,i-1) - d1
189 22897 : qtmp(j,i) = qtmp(j,i) + d1
190 : ! east
191 22897 : d0 = max(D0_0,qtmp(j,1))
192 22897 : d2 = min(-qtmp(j,i),d0)
193 22897 : qtmp(j,1) = qtmp(j,1) - d2
194 22897 : qtmp(j,i) = qtmp(j,i) + d2 + tinyl
195 : endif
196 : enddo
197 :
198 41287680 : if(ipx .ne. 0) then
199 :
200 : !-----------
201 : ! Final pass
202 : !-----------
203 472701312 : do i=1,im-1
204 1875950502 : do j=j1,j2
205 1874309178 : if (qtmp(j,i) < D0_0 ) then
206 : ! Take mass from east (essentially adjusting fx(i+1,j))
207 12613412 : qtmp(j,i+1) = qtmp(j,i+1) + qtmp(j,i)
208 12613412 : qtmp(j,i) = D0_0
209 : endif
210 : enddo
211 : enddo
212 :
213 472701312 : do i=im,2,-1
214 1875950502 : do j=j1,j2
215 1874309178 : if (qtmp(j,i) < D0_0 ) then
216 : ! Take mass from west (essentially adjusting fx(i,j))
217 4826186 : qtmp(j,i-1) = qtmp(j,i-1) + qtmp(j,i)
218 4826186 : qtmp(j,i) = D0_0
219 : endif
220 : enddo
221 : enddo
222 :
223 6530694 : do j=j1,j2
224 :
225 4889370 : qmin = D0_0
226 1413027930 : do i=1, im
227 1413027930 : if (qtmp(j,i) < qmin) then
228 1005 : qmin = qtmp(j,i)
229 1005 : imin = i
230 1005 : jmin = j
231 : endif
232 : enddo
233 :
234 4889370 : if (( qmin < D0_0 ) .and. print_filew_warn == "full") then
235 0 : write(iulog,*) ' filew failed, worst i, j, qtmp, q = ', imin, jmin, qtmp(jmin,imin), q(imin,jmin)
236 : end if
237 :
238 1414669254 : do i=1,im
239 1413027930 : q(i,j) = qtmp(j,i)
240 : enddo
241 : enddo
242 :
243 : endif
244 :
245 : ! Check Poles.
246 :
247 41287680 : if ( jfirst == 1 ) then
248 645120 : if(q(1,1) < D0_0) then
249 0 : call pfix(q(1,2),q(1,1),im,ipx,acap,cosp2)
250 : else
251 : ! Check j=2
252 645120 : ip2 = 0
253 186439680 : do i=1,im
254 186439680 : if(q(i,2).lt.D0_0) then
255 : ip2 = 1
256 : go to 322
257 : endif
258 : enddo
259 : 322 continue
260 645120 : if(ip2.ne.0) call pfix(q(1,2),q(1,1),im,ipx,acap,cosp2)
261 : endif
262 : endif
263 :
264 41287680 : if ( jlast == jm ) then
265 645120 : if(q(1,jm) < D0_0) then
266 10 : call pfix(q(1,jm1),q(1,jm),im,ipx,acap,cosp2)
267 : else
268 : ! Check j=jm1
269 645110 : ip2 = 0
270 186436790 : do i=1,im
271 186436790 : if(q(i,jm1) < D0_0) then
272 : ip2 = 1
273 : go to 323
274 : endif
275 : enddo
276 : 323 continue
277 645110 : if(ip2.ne.0) call pfix(q(1,jm1),q(1,jm),im,ipx,acap,cosp2)
278 : endif
279 : endif
280 :
281 : !EOC
282 41287680 : end subroutine filew
283 : !-----------------------------------------------------------------------
284 :
285 :
286 : !-----------------------------------------------------------------------
287 : !BOP
288 : ! !IROUTINE: fillxy --- Fill from east, west, north and south neighbors
289 : !
290 : ! !INTERFACE:
291 41287680 : subroutine fillxy(q, im, jm, jfirst, jlast, acap, cosp, acosp)
292 :
293 : ! !USES:
294 :
295 : implicit none
296 :
297 : integer im ! Longitudes
298 : integer jm ! Total latitudes
299 : integer jfirst ! Starting latitude
300 : integer jlast ! Finishing latitude
301 :
302 : real(r8) acap ! ???
303 : real(r8) cosp(jm) ! ???
304 : real(r8) acosp(jm) ! ???
305 : !
306 : ! !INPUT/OUTPUT PARAMETERS:
307 : real(r8) q(im,jfirst:jlast) ! Field to adjust
308 :
309 : ! !DESCRIPTION:
310 : ! Check for "bad" data and fill from east and west neighbors
311 : !
312 : ! !BUGS:
313 : ! Currently this routine only performs the east-west fill algorithm.
314 : ! This is because the N-S fill is very hard to do in a reproducible
315 : ! fashion when the problem is decomposed by latitudes.
316 : !
317 : ! !REVISION HISTORY:
318 : ! 99.03.01 Lin Creation
319 : !
320 : !EOP
321 : !-----------------------------------------------------------------------
322 : !BOC
323 : !
324 : ! !LOCAL VARIABLES:
325 : integer ipx, ipy, j1, j2
326 : real(r8) tiny
327 : parameter( tiny = 1.e-20_r8 )
328 :
329 41287680 : call filew(q,im,jm,jfirst,jlast,acap,ipx,tiny,cosp(2))
330 :
331 : ! WS 99.08.03 : S.-J. can you clean up the j1, j2 stuff here?
332 41287680 : if(ipx.ne.0) then
333 :
334 41287680 : j1 = max( 2, jfirst )
335 41287680 : j2 = min( jm-1, jlast )
336 : !
337 : ! WS 99.08.03 : see comments in "BUGS" above
338 : !!! call filns(q,im,jm,j1,j2,cosp,acosp,ipy,tiny)
339 :
340 : ! if(ipy .ne. 0) then
341 : ! do fill zonally
342 : ! xfx is problematic
343 : ! call xfix(q,IM,JM,tiny,qt)
344 : ! endif
345 :
346 : endif
347 :
348 : !EOC
349 41287680 : end subroutine fillxy
350 : !-----------------------------------------------------------------------
351 :
352 :
353 : !-----------------------------------------------------------------------
354 : !BOP
355 : ! !IROUTINE: fillz --- Fill from neighbors below and above
356 : !
357 : ! !INTERFACE:
358 0 : subroutine fillz(im, i1, i2, km, nq, q, dp)
359 :
360 : ! !USES:
361 :
362 : implicit none
363 :
364 : ! !INPUT PARAMETERS:
365 : integer, intent(in) :: im ! No. of longitudes
366 : integer, intent(in) :: km ! No. of levels
367 : integer, intent(in) :: i1 ! Starting longitude
368 : integer, intent(in) :: i2 ! Finishing longitude
369 : integer, intent(in) :: nq ! Total number of tracers
370 : real(r8), intent(in) :: dp(im,km) ! pressure thickness
371 :
372 : ! !INPUT/OUTPUT PARAMETERS:
373 : real(r8), intent(inout) :: q(im,km,nq) ! tracer mixing ratio
374 :
375 : ! !DESCRIPTION:
376 : ! Check for "bad" data and fill from east and west neighbors
377 : !
378 : ! !BUGS:
379 : ! Currently this routine only performs the east-west fill algorithm.
380 : ! This is because the N-S fill is very hard to do in a reproducible
381 : ! fashion when the problem is decomposed by latitudes.
382 : !
383 : ! !REVISION HISTORY:
384 : ! 00.04.01 Lin Creation
385 : !
386 : !EOP
387 : !-----------------------------------------------------------------------
388 : !BOC
389 : !
390 : ! !LOCAL VARIABLES:
391 : integer i, k, ic
392 : real(r8) qup, qly, dup
393 :
394 0 : do ic=1,nq
395 : ! Top layer
396 0 : do i=i1,i2
397 0 : if( q(i,1,ic) < D0_0) then
398 0 : q(i,2,ic) = q(i,2,ic) + q(i,1,ic)*dp(i,1)/dp(i,2)
399 0 : q(i,1,ic) = D0_0
400 : endif
401 : enddo
402 :
403 : ! Interior
404 0 : do k=2,km-1
405 0 : do i=i1,i2
406 0 : if( q(i,k,ic) < D0_0 ) then
407 : ! Borrow from above
408 0 : qup = q(i,k-1,ic)*dp(i,k-1)
409 0 : qly = -q(i,k ,ic)*dp(i,k )
410 0 : dup = min( D0_5*qly, qup ) !borrow no more than 50%
411 0 : q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1)
412 : ! Borrow from below: q(i,k,ic) is still negative at this stage
413 0 : q(i,k+1,ic) = q(i,k+1,ic) + (dup-qly)/dp(i,k+1)
414 0 : q(i,k ,ic) = D0_0
415 : endif
416 : enddo
417 : enddo
418 :
419 : ! Bottom layer
420 0 : k = km
421 0 : do i=i1,i2
422 0 : if( q(i,k,ic) < D0_0) then
423 : ! Borrow from above
424 0 : qup = q(i,k-1,ic)*dp(i,k-1)
425 0 : qly = -q(i,k ,ic)*dp(i,k )
426 0 : dup = min( qly, qup )
427 0 : q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1)
428 0 : q(i,k,ic) = D0_0
429 : endif
430 : enddo
431 : enddo
432 : !EOC
433 0 : end subroutine fillz
434 : !-----------------------------------------------------------------------
435 :
436 :
437 : !-----------------------------------------------------------------------
438 : !BOP
439 : ! !IROUTINE: filns --- Fill from north and south neighbors
440 : !
441 : ! !INTERFACE:
442 0 : subroutine filns(q,im,jm,j1,j2,cosp,acosp,ipy,tiny)
443 :
444 : ! !USES:
445 :
446 : implicit none
447 :
448 : ! !INPUT PARAMETERS:
449 : integer im ! Longitudes
450 : integer jm ! Total latitudes
451 : integer j1 ! Starting latitude
452 : integer j2 ! Finishing latitude
453 :
454 :
455 : real(r8) tiny ! A small number to pump up value
456 : real(r8) cosp(*) ! ???
457 : real(r8) acosp(*) ! ???
458 :
459 : ! !INPUT/OUTPUT PARAMETERS:
460 : real(r8) q(im,*) ! Field to adjust
461 :
462 : ! !OUTPUT PARAMETERS:
463 : integer ipy ! Flag: 0 if no fill-in, 1 if fill-in
464 :
465 : ! !DESCRIPTION:
466 : ! Check for "bad" data and fill from north and south neighbors
467 : !
468 : ! !BUGS:
469 : ! Currently this routine can only be used performs when the
470 : ! problem is *not* distributed in latitude (i.e. j1=1, j2=jm).
471 : ! This is because the N-S fill is very hard to do in a reproducible
472 : ! fashion when the problem is decomposed by latitudes.
473 : !
474 : ! !REVISION HISTORY:
475 : ! 99.03.01 Lin Creation
476 : ! 05.06.30 Sawyer Removed SAVE attribute for cap1 (recalculated)
477 : !
478 : !EOP
479 : !-----------------------------------------------------------------------
480 : !BOC
481 : !
482 : ! !LOCAL VARIABLES:
483 : integer i, j
484 : ! This definition of PI as opposed to 4._r16*atan(1._r16) does not
485 : ! appear to generate non-zero differences in GEOS5 checkpoint files
486 : real(R16),parameter :: pi = 3.1415926535897932384626433832795028841971_R16
487 : real(r8) :: dp, cap1, dq, dn, ds, d0, d1, d2
488 :
489 0 : dp = pi/real(jm-1,r16)
490 0 : cap1 = im*(D1_0-cos((j1-D1_5)*dp))/dp
491 :
492 0 : ipy = 0
493 0 : do j=j1+1,j2-1
494 0 : do i=1,im
495 0 : if(q(i,j).lt.D0_0) then
496 0 : ipy = 1
497 0 : dq = - q(i,j)*cosp(j)
498 : ! North
499 0 : dn = q(i,j+1)*cosp(j+1)
500 0 : d0 = max(D0_0,dn)
501 0 : d1 = min(dq,d0)
502 0 : q(i,j+1) = (dn - d1)*acosp(j+1)
503 0 : dq = dq - d1
504 : ! South
505 0 : ds = q(i,j-1)*cosp(j-1)
506 0 : d0 = max(D0_0,ds)
507 0 : d2 = min(dq,d0)
508 0 : q(i,j-1) = (ds - d2)*acosp(j-1)
509 0 : q(i,j) = (d2 - dq)*acosp(j) + tiny
510 : endif
511 : enddo
512 : enddo
513 :
514 0 : do i=1,im
515 0 : if(q(i,j1).lt.D0_0) then
516 0 : ipy = 1
517 0 : dq = - q(i,j1)*cosp(j1)
518 : ! North
519 0 : dn = q(i,j1+1)*cosp(j1+1)
520 0 : d0 = max(D0_0,dn)
521 0 : d1 = min(dq,d0)
522 0 : q(i,j1+1) = (dn - d1)*acosp(j1+1)
523 0 : q(i,j1) = (d1 - dq)*acosp(j1) + tiny
524 : endif
525 : enddo
526 :
527 0 : j = j2
528 0 : do i=1,im
529 0 : if(q(i,j).lt.D0_0) then
530 0 : ipy = 1
531 0 : dq = - q(i,j)*cosp(j)
532 : ! South
533 0 : ds = q(i,j-1)*cosp(j-1)
534 0 : d0 = max(D0_0,ds)
535 0 : d2 = min(dq,d0)
536 0 : q(i,j-1) = (ds - d2)*acosp(j-1)
537 0 : q(i,j) = (d2 - dq)*acosp(j) + tiny
538 : endif
539 : enddo
540 :
541 : ! Check Poles.
542 0 : if(q(1,1).lt.D0_0) then
543 0 : dq = q(1,1)*cap1/real(im,r8)*acosp(j1)
544 0 : do i=1,im
545 0 : q(i,1) = tiny
546 0 : q(i,j1) = q(i,j1) + dq
547 0 : q(i,j1) = max(tiny, q(i,j1) + dq )
548 : enddo
549 : endif
550 :
551 0 : if(q(1,jm).lt.D0_0) then
552 0 : dq = q(1,jm)*cap1/real(im,r8)*acosp(j2)
553 0 : do i=1,im
554 0 : q(i,jm) = tiny
555 0 : q(i,j2) = max(tiny, q(i,j2) + dq )
556 : enddo
557 : endif
558 : !EOC
559 0 : end subroutine filns
560 : !-----------------------------------------------------------------------
561 :
562 : !-----------------------------------------------------------------------
563 : !BOP
564 : ! !IROUTINE: pfix --- fix an individual latitude-level
565 : !
566 : ! !INTERFACE:
567 10 : subroutine pfix(q, qp, im, ipx, acap, cosp2)
568 :
569 : ! !USES:
570 : implicit none
571 :
572 : ! !INPUT PARAMETERS:
573 : integer im ! Longitudes
574 : real(r8) acap ! ???
575 : real(r8) cosp2 ! ???
576 :
577 : ! !INPUT/OUTPUT PARAMETERS:
578 : real(r8) q(im) ! Latitude-level field to adjust
579 : real(r8) qp(im) ! Second latitude-level field to adjust (usually pole)
580 :
581 : ! !OUTPUT PARAMETERS:
582 : integer ipx ! Flag: 0 if Q not change, 1 if changed
583 :
584 :
585 : ! !DESCRIPTION:
586 : ! Fill one latitude-level from east and west neighbors
587 : !
588 : ! !REVISION HISTORY:
589 : ! 99.03.01 Lin Creation
590 : !
591 : !EOP
592 : !-----------------------------------------------------------------------
593 : !BOC
594 : !
595 : ! !LOCAL VARIABLES:
596 : integer i
597 : real(r8) summ, sump, pmean
598 :
599 10 : summ = D0_0
600 10 : sump = D0_0
601 2890 : do i=1,im
602 2880 : summ = summ + q(i)
603 2890 : sump = sump + qp(i)
604 : enddo
605 :
606 10 : sump = sump/im
607 10 : pmean = (sump*acap + summ*cosp2) / (acap + cosp2*im)
608 :
609 2890 : do i=1,im
610 2880 : q(i) = pmean
611 2890 : qp(i) = pmean
612 : enddo
613 :
614 10 : if( qp(1) < D0_0 ) then
615 6 : ipx = 1
616 : endif
617 :
618 : !EOC
619 10 : end subroutine pfix
620 : !-----------------------------------------------------------------------
621 :
622 : end module fill_module
|
