grid_tools 1.14.0
Loading...
Searching...
No Matches
filter_topo.F90
Go to the documentation of this file.
1
4
10program filter_topo
11
12 use omp_lib
13 use utils
14
15 implicit none
16
17#include <netcdf.inc>
18
19#ifdef NO_QUAD_PRECISION
20 ! 64-bit precision (kind=8)
21 integer, parameter:: f_p = selected_real_kind(15)
22#else
23 ! Higher precision (kind=16) for grid geometrical factors:
24 integer, parameter:: f_p = selected_real_kind(20)
25#endif
26
27 integer, parameter :: xdir=1
28 integer, parameter :: ydir=2
29 real, parameter :: pi = 3.14159265358979323846d0
30 real, parameter :: radius = 6371200.0
31 real, parameter :: big_number=1.d8
32 real, parameter :: tiny_number=1.d-8
33
34 real:: cd4 ! Dimensionless coeff for del-4 difussion (with FCT)
35 real:: peak_fac ! overshoot factor for the mountain peak
36 real:: max_slope ! max allowable terrain slope: 1 --> 45 deg
37
38 integer :: n_del2_weak, tid, nthreads
39
40 integer :: ntiles = 0
41
42 real da_min
43
44 real, allocatable :: oro(:,:,:), mask(:,:,:)
45 real, allocatable :: dx(:,:,:), dy(:,:,:)
46 real, allocatable :: dxa(:,:,:), dya(:,:,:)
47 real, allocatable :: dxc(:,:,:), dyc(:,:,:)
48 real, allocatable :: area(:,:,:)
49 real, allocatable :: sin_sg(:,:,:,:)
50
51 integer :: is,ie,js,je,isd,ied,jsd,jed
52 integer,parameter :: ng = 3
53 integer :: nx, ny, npx, npy, nx_nest, ny_nest, npx_nest, npy_nest, is_nest, ie_nest, js_nest, je_nest, isd_nest, ied_nest, jsd_nest, jed_nest
54
55!$OMP PARALLEL PRIVATE(TID)
56 tid = omp_get_thread_num()
57 if (tid == 0) then
58 nthreads = omp_get_num_threads()
59 print*
60 print*,'- BEGIN EXECUTION WITH NUMBER OF THREADS = ',nthreads
61 print*
62 endif
63!$OMP END PARALLEL
64
65 !--- read namelist
66 call read_namelist()
67
68 !--- compute filter constants according to grid resolution.
69 call compute_filter_constants
70
71 !--- read the target grid.
72 call read_grid_file(regional)
73
74 !--- read the topography data
75 call read_topo_file(regional)
76
77 !--- filter the data
78 call fv3_zs_filter(is,ie,js,je,isd,ied,jsd,jed,npx,npy,npx,ntiles,grid_type, &
79 stretch_fac, nested, area, dxa, dya, dx, dy, dxc, dyc, sin_sg, oro, regional )
80
81 !--- write out the data
82 call write_topo_file(is,ie,js,je,ntiles,oro(is:ie,js:je,:),regional )
83
84 print*
85 print*,'- NORMAL TERMINATION.'
86
87contains
88
96 real function great_circle_dist( q1, q2, radius )
97 real, intent(IN) :: q1(2), q2(2)
98 real, intent(IN), optional :: radius
99
100 real (f_p):: p1(2), p2(2)
101 real (f_p):: beta
102 integer n
103
104 do n=1,2
105 p1(n) = q1(n)
106 p2(n) = q2(n)
107 enddo
108
109 beta = asin( sqrt( sin((p1(2)-p2(2))/2.)**2 + cos(p1(2))*cos(p2(2))* &
110 sin((p1(1)-p2(1))/2.)**2 ) ) * 2.
111
112 if ( present(radius) ) then
113 great_circle_dist = radius * beta
114 else
115 great_circle_dist = beta ! Returns the angle
116 endif
117
118 end function great_circle_dist
119
137 real function spherical_angle(p1, p2, p3)
138
139 real p1(3), p2(3), p3(3)
140
141 real (f_p):: e1(3), e2(3), e3(3)
142 real (f_p):: px, py, pz
143 real (f_p):: qx, qy, qz
144 real (f_p):: angle, ddd
145 integer n
146
147 do n=1,3
148 e1(n) = p1(n)
149 e2(n) = p2(n)
150 e3(n) = p3(n)
151 enddo
152
153 !-------------------------------------------------------------------
154 ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
155 !-------------------------------------------------------------------
156 ! Vector P:
157 px = e1(2)*e2(3) - e1(3)*e2(2)
158 py = e1(3)*e2(1) - e1(1)*e2(3)
159 pz = e1(1)*e2(2) - e1(2)*e2(1)
160 ! Vector Q:
161 qx = e1(2)*e3(3) - e1(3)*e3(2)
162 qy = e1(3)*e3(1) - e1(1)*e3(3)
163 qz = e1(1)*e3(2) - e1(2)*e3(1)
164
165 ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz)
166
167 if ( ddd <= 0.0d0 ) then
168 angle = 0.d0
169 else
170 ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd)
171 if ( abs(ddd)>1.d0) then
172 angle = 2.d0*atan(1.0) ! 0.5*pi
173 !FIX (lmh) to correctly handle co-linear points (angle near pi or 0)
174 if (ddd < 0.d0) then
175 angle = 4.d0*atan(1.0d0) !should be pi
176 else
177 angle = 0.d0
178 end if
179 else
180 angle = acos( ddd )
181 endif
182 endif
183
184 spherical_angle = angle
185
186 end function spherical_angle
187
188
198 real function get_area(p1, p4, p2, p3, radius)
199 !-----------------------------------------------
200 real, intent(in), dimension(2):: p1, p2, p3, p4
201 real, intent(in), optional:: radius
202 !-----------------------------------------------
203 real e1(3), e2(3), e3(3)
204 real ang1, ang2, ang3, ang4
205
206 ! S-W: 1
207 call latlon2xyz(p1, e1) ! p1
208 call latlon2xyz(p2, e2) ! p2
209 call latlon2xyz(p4, e3) ! p4
210 ang1 = spherical_angle(e1, e2, e3)
211 !----
212 ! S-E: 2
213 !----
214 call latlon2xyz(p2, e1)
215 call latlon2xyz(p3, e2)
216 call latlon2xyz(p1, e3)
217 ang2 = spherical_angle(e1, e2, e3)
218 !----
219 ! N-E: 3
220 !----
221 call latlon2xyz(p3, e1)
222 call latlon2xyz(p4, e2)
223 call latlon2xyz(p2, e3)
224 ang3 = spherical_angle(e1, e2, e3)
225 !----
226 ! N-W: 4
227 !----
228 call latlon2xyz(p4, e1)
229 call latlon2xyz(p3, e2)
230 call latlon2xyz(p1, e3)
231 ang4 = spherical_angle(e1, e2, e3)
232
233 if ( present(radius) ) then
234 get_area = (ang1 + ang2 + ang3 + ang4 - 2.*pi) * radius**2
235 else
236 get_area = ang1 + ang2 + ang3 + ang4 - 2.*pi
237 endif
238
239 end function get_area
240
251 subroutine fill_agrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
252 integer, intent(in) :: ng, npx, npy, isd, jsd
253 integer, intent(in) :: fill
254 real, DIMENSION(isd:,jsd:,:), intent(INOUT):: q
255
256 integer :: i, j
257
258 select case (fill)
259 case (xdir)
260 do j=1,ng
261 do i=1,ng
262 q(1-i ,1-j ,:) = q(1-j ,i ,:) !SW Corner
263 q(1-i ,npy-1+j,:) = q(1-j ,npy-1-i+1,:) !NW Corner
264 q(npx-1+i,1-j ,:) = q(npx-1+j,i ,:) !SE Corner
265 q(npx-1+i,npy-1+j,:) = q(npx-1+j,npy-1-i+1,:) !NE Corner
266 enddo
267 enddo
268 case (ydir)
269 do j=1,ng
270 do i=1,ng
271 q(1-j ,1-i ,:) = q(i ,1-j ,:) !SW Corner
272 q(1-j ,npy-1+i,:) = q(i ,npy-1+j,:) !NW Corner
273 q(npx-1+j,1-i ,:) = q(npx-1-i+1,1-j ,:) !SE Corner
274 q(npx-1+j,npy-1+i,:) = q(npx-1-i+1,npy-1+j,:) !NE Corner
275 enddo
276 enddo
277 case default
278 do j=1,ng
279 do i=1,ng
280 q(1-j ,1-i ,:) = q(i ,1-j ,:) !SW Corner
281 q(1-j ,npy-1+i,:) = q(i ,npy-1+j,:) !NW Corner
282 q(npx-1+j,1-i ,:) = q(npx-1-i+1,1-j ,:) !SE Corner
283 q(npx-1+j,npy-1+i,:) = q(npx-1-i+1,npy-1+j,:) !NE Corner
284 enddo
285 enddo
286 end select
287
288
289 end subroutine fill_agrid_scalar_corners
290
301 subroutine fill_bgrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
302 integer, intent(in) :: ng, npx, npy, isd, jsd
303 integer, intent(in) :: fill
304 real, DIMENSION(isd:,jsd:,:), intent(INOUT):: q
305
306 integer :: i, j
307
308 select case (fill)
309 case (xdir)
310 do j=1,ng
311 do i=1,ng
312 q(1-i ,1-j ,:) = q(1-j ,i+1 ,:) !SW Corner
313 q(1-i ,npy+j,:) = q(1-j ,npy-i ,:) !NW Corner
314 q(npx+i,1-j ,:) = q(npx+j,i+1 ,:) !SE Corner
315 q(npx+i,npy+j,:) = q(npx+j,npy-i ,:) !NE Corner
316 enddo
317 enddo
318 case (ydir)
319 do j=1,ng
320 do i=1,ng
321 q(1-j ,1-i ,:) = q(i+1 ,1-j ,:) !SW Corner
322 q(1-j ,npy+i,:) = q(i+1 ,npy+j ,:) !NW Corner
323 q(npx+j,1-i ,:) = q(npx-i,1-j ,:) !SE Corner
324 q(npx+j,npy+i,:) = q(npx-i,npy+j ,:) !NE Corner
325 enddo
326 enddo
327 case default
328 do j=1,ng
329 do i=1,ng
330 q(1-i ,1-j ,:) = q(1-j ,i+1 ,:) !SW Corner
331 q(1-i ,npy+j,:) = q(1-j ,npy-i ,:) !NW Corner
332 q(npx+i,1-j ,:) = q(npx+j,i+1 ,:) !SE Corner
333 q(npx+i,npy+j,:) = q(npx+j,npy-i ,:) !NE Corner
334 enddo
335 enddo
336 end select
337
338
339
340 end subroutine fill_bgrid_scalar_corners
341
352 subroutine fill_agrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
353 integer, intent(in) :: ng, npx, npy, isd, jsd
354 real, DIMENSION(isd:,jsd:,:), intent(INOUT):: x
355 real, DIMENSION(isd:,jsd:,:), intent(INOUT):: y
356 integer :: i,j
357
358 do j=1,ng
359 do i=1,ng
360 x(1-i ,1-j ,:) = y(1-j ,i ,:) !SW Corner
361 x(1-i ,npy-1+j,:) = y(1-j ,npy-1-i+1,:) !NW Corner
362 x(npx-1+i,1-j ,:) = y(npx-1+j,i ,:) !SE Corner
363 x(npx-1+i,npy-1+j,:) = y(npx-1+j,npy-1-i+1,:) !NE Corner
364
365 y(1-j ,1-i ,:) = x(i ,1-j ,:) !SW Corner
366 y(1-j ,npy-1+i,:) = x(i ,npy-1+j,:) !NW Corner
367 y(npx-1+j,1-i ,:) = x(npx-1-i+1,1-j ,:) !SE Corner
368 y(npx-1+j,npy-1+i,:) = x(npx-1-i+1,npy-1+j,:) !NE Corner
369 enddo
370 enddo
371
372 end subroutine fill_agrid_xy_corners
373
384 subroutine fill_dgrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
385 integer, intent(in) :: ng, npx, npy, isd, jsd
386 real, DIMENSION(isd:,jsd:,:), intent(INOUT):: x
387 real, DIMENSION(isd:,jsd:,:), intent(INOUT):: y
388 integer :: i,j
389
390 do j=1,ng
391 do i=1,ng
392 x(1-i ,1-j , :) = y(1-j ,i , :) !SW Corner
393 x(1-i ,npy+j , :) = y(1-j ,npy-i, :) !NW Corner
394 x(npx-1+i,1-j , :) = y(npx+j,i , :) !SE Corner
395 x(npx-1+i,npy+j , :) = y(npx+j,npy-i, :) !NE Corner
396 y(1-i ,1-j , :) = x(j ,1-i , :) !SW Corner
397 y(1-i ,npy-1+j, :) = x(j ,npy+i, :) !NW Corner
398 y(npx+i ,1-j , :) = x(npx-j,1-i , :) !SE Corner
399 y(npx+i ,npy-1+j, :) = x(npx-j,npy+i, :) !NE Corner
400 enddo
401 enddo
402
403 end subroutine fill_dgrid_xy_corners
404
411 subroutine mid_pt_sphere(p1, p2, pm)
412 real, intent(IN) :: p1(2), p2(2)
413 real, intent(OUT) :: pm(2)
414 !------------------------------------------
415 real :: e1(3), e2(3), e3(3)
416
417 call latlon2xyz(p1, e1)
418 call latlon2xyz(p2, e2)
419 call mid_pt3_cart(e1, e2, e3)
420 call cart_to_latlon(1, e3, pm(1), pm(2))
421
422 end subroutine mid_pt_sphere
423
430 subroutine mid_pt3_cart(p1, p2, e)
431 real, intent(IN) :: p1(3), p2(3)
432 real, intent(OUT) :: e(3)
433 !
434 real (f_p):: q1(3), q2(3)
435 real (f_p):: dd, e1, e2, e3
436 integer k
437
438 do k=1,3
439 q1(k) = p1(k)
440 q2(k) = p2(k)
441 enddo
442
443 e1 = q1(1) + q2(1)
444 e2 = q1(2) + q2(2)
445 e3 = q1(3) + q2(3)
446
447 dd = sqrt( e1**2 + e2**2 + e3**2 )
448 e1 = e1 / dd
449 e2 = e2 / dd
450 e3 = e3 / dd
451
452 e(1) = e1
453 e(2) = e2
454 e(3) = e3
455
456 end subroutine mid_pt3_cart
457
463 subroutine latlon2xyz(p, e)
464 !
465 ! Routine to map (lon, lat) to (x,y,z)
466 !
467 real, intent(in) :: p(2)
468 real, intent(out):: e(3)
469
470 integer n
471 real (f_p):: q(2)
472 real (f_p):: e1, e2, e3
473
474 do n=1,2
475 q(n) = p(n)
476 enddo
477
478 e1 = cos(q(2)) * cos(q(1))
479 e2 = cos(q(2)) * sin(q(1))
480 e3 = sin(q(2))
481 !-----------------------------------
482 ! Truncate to the desired precision:
483 !-----------------------------------
484 e(1) = e1
485 e(2) = e2
486 e(3) = e3
487
488 end subroutine latlon2xyz
489
497 subroutine cart_to_latlon(np, q, xs, ys)
498 ! vector version of cart_to_latlon1
499 integer, intent(in):: np
500 real, intent(inout):: q(3,np)
501 real, intent(inout):: xs(np), ys(np)
502 ! local
503 real, parameter:: esl=1.d-10
504 real (f_p):: p(3)
505 real (f_p):: dist, lat, lon
506 integer i,k
507
508 do i=1,np
509 do k=1,3
510 p(k) = q(k,i)
511 enddo
512 dist = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
513 do k=1,3
514 p(k) = p(k) / dist
515 enddo
516
517 if ( (abs(p(1))+abs(p(2))) < esl ) then
518 lon = real(0.,kind=f_p)
519 else
520 lon = atan2( p(2), p(1) ) ! range [-pi,pi]
521 endif
522
523 if ( lon < 0.) lon = real(2.,kind=f_p)*pi + lon
524 ! RIGHT_HAND system:
525 lat = asin(p(3))
526
527 xs(i) = lon
528 ys(i) = lat
529 ! q Normalized:
530 do k=1,3
531 q(k,i) = p(k)
532 enddo
533 enddo
534
535 end subroutine cart_to_latlon
536
544 real function cos_angle(p1, p2, p3)
545 ! As spherical_angle, but returns the cos(angle)
546 ! p3
547 ! ^
548 ! |
549 ! |
550 ! p1 ---> p2
551 !
552 real, intent(in):: p1(3), p2(3), p3(3)
553
554 real (f_p):: e1(3), e2(3), e3(3)
555 real (f_p):: px, py, pz
556 real (f_p):: qx, qy, qz
557 real (f_p):: angle, ddd
558 integer n
559
560 do n=1,3
561 e1(n) = p1(n)
562 e2(n) = p2(n)
563 e3(n) = p3(n)
564 enddo
565
566 !-------------------------------------------------------------------
567 ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
568 !-------------------------------------------------------------------
569 ! Vector P:= e1 X e2
570 px = e1(2)*e2(3) - e1(3)*e2(2)
571 py = e1(3)*e2(1) - e1(1)*e2(3)
572 pz = e1(1)*e2(2) - e1(2)*e2(1)
573
574 ! Vector Q: e1 X e3
575 qx = e1(2)*e3(3) - e1(3)*e3(2)
576 qy = e1(3)*e3(1) - e1(1)*e3(3)
577 qz = e1(1)*e3(2) - e1(2)*e3(1)
578
579 ! ddd = sqrt[ (P*P) (Q*Q) ]
580 ddd = sqrt( (px**2+py**2+pz**2)*(qx**2+qy**2+qz**2) )
581 if ( ddd > 0.d0 ) then
582 angle = (px*qx+py*qy+pz*qz) / ddd
583 else
584 angle = 1.d0
585 endif
586 cos_angle = angle
587
588 end function cos_angle
589
598 subroutine cell_center2(q1, q2, q3, q4, e2)
599 real, intent(in ) :: q1(2), q2(2), q3(2), q4(2)
600 real, intent(out) :: e2(2)
601 ! Local
602 real p1(3), p2(3), p3(3), p4(3)
603 real ec(3)
604 real dd
605 integer k
606
607 call latlon2xyz(q1, p1)
608 call latlon2xyz(q2, p2)
609 call latlon2xyz(q3, p3)
610 call latlon2xyz(q4, p4)
611
612 do k=1,3
613 ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
614 enddo
615 dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
616
617 do k=1,3
618 ec(k) = ec(k) / dd
619 enddo
620
621 call cart_to_latlon(1, ec, e2(1), e2(2))
622
623 end subroutine cell_center2
624
629 subroutine read_grid_file(regional)
630
631 logical, intent(in) :: regional ! Is this a regional run?
632 integer :: fsize=65536
633 integer :: status, ncid, id_dim, id_var, ncid2, t
634 integer :: ni, nj, i, j, tw, te, ip
635 real :: g1(2), g2(2), g3(2), g4(2), g5(2)
636 real :: p1(3), p3(3)
637 real :: p_lL(2), p_uL(2), p_lR(2), p_uR(2)
638 character(len=512) :: tile_file
639 real, allocatable, dimension(:,:) :: tmpvar, geolon_c_nest, geolat_c_nest
640 real, allocatable, dimension(:,:,:) :: geolon_c, geolat_c
641 real, allocatable, dimension(:,:,:) :: geolon_t, geolat_t, cos_sg, grid3
642 integer :: start(4), nread(4)
643
644 print*, "Read the grid from file "//trim(grid_file)
645
646 status=nf__open(trim(grid_file),nf_nowrite,fsize,ncid)
647 call handle_err(status, 'Open file '//trim(grid_file) )
648
649 status=nf_inq_dimid(ncid, 'ntiles', id_dim)
650 call handle_err(status, 'inquire dimension ntiles from file '//trim(grid_file) )
651 status=nf_inq_dimlen(ncid,id_dim,ntiles)
652 call handle_err(status, 'inquire dimension ntiles length from file '//trim(grid_file) )
653
654 if( ntiles == 6) then
655 print*, " read_grid_file: This is a global grid."
656 elseif( ntiles >= 6 )then
657 print*, " read_grid_file: This is a nested grid."
658 print*, " filter only the globe."
659 ntiles=6
660 elseif( ntiles == 1 )then
661 print*, " read_grid_file: This is a standalone regional grid."
662 endif
663
664 !--- loop through ntiles and make sure the grid size match between all the tiles.
665
666 start(:) = 1
667 nread(:) = 1
668
669 do t = 1, ntiles
670 start(2) = t; nread(1) = 255
671 status = nf_inq_varid(ncid, 'gridfiles', id_var)
672 call handle_err(status, 'inquire varid of gridfiles from file '//trim(grid_file) )
673
674 !--- Obtain the grid file name from the mosaic file.
675 status = nf_get_vara_text(ncid, id_var, start, nread, tile_file )
676 call handle_err(status, 'get value of gridfiles from file '//trim(grid_file) )
677
678 status=nf__open(trim(tile_file),nf_nowrite,fsize,ncid2)
679 call handle_err(status, 'Open file '//trim(tile_file) )
680
681 status=nf_inq_dimid(ncid2, 'nx', id_dim)
682 call handle_err(status, 'inquire dimension nx from file '//trim(grid_file) )
683 status=nf_inq_dimlen(ncid2,id_dim,ni)
684 call handle_err(status, 'inquire dimension nx length from file '//trim(grid_file) )
685 status=nf_inq_dimid(ncid2, 'ny', id_dim)
686 call handle_err(status, 'inquire dimension ny from file '//trim(grid_file) )
687 status=nf_inq_dimlen(ncid2,id_dim,nj)
688 call handle_err(status, 'inquire dimension ny length '//'from file '//trim(grid_file) )
689 if( t == 1 ) then
690 ! ni and nj must be even
691 if(mod(ni,2) .NE. 0 .or. mod(nj,2) .NE. 0) &
692 call handle_err(-1, "read_grid_file: ni and nj must be even")
693
694 nx = ni/2
695 ny = nj/2
696 npx = nx + 1
697 npy = ny + 1
698 is = 1 ; ie = nx
699 js = 1 ; je = ny
700 isd=is-ng; ied=ie+ng
701 jsd=js-ng; jed=je+ng
702
703 allocate(tmpvar(ni+1,nj+1))
704 allocate(geolon_c(isd:ied+1,jsd:jed+1,6))
705 allocate(geolat_c(isd:ied+1,jsd:jed+1,6))
706 else if ( t == 7 ) then ! nested grid
707 if(mod(ni,2) .NE. 0 .or. mod(nj,2) .NE. 0) &
708 call handle_err(-1, "read_grid_file: ni and nj must be even")
709
710 nx_nest = ni/2
711 ny_nest = nj/2
712 npx_nest = nx_nest + 1
713 npy_nest = ny_nest + 1
714 is_nest = 1 ; ie_nest = nx
715 js_nest = 1 ; je_nest = ny
716 isd_nest=is_nest-ng; ied_nest=ie_nest+ng
717 jsd_nest=js_nest-ng; jed_nest=je_nest+ng
718 deallocate(tmpvar)
719 allocate(tmpvar(ni+1,nj+1))
720 allocate(geolon_c_nest(isd:ied+1,jsd:jed+1))
721 allocate(geolat_c_nest(isd:ied+1,jsd:jed+1))
722 else
723 !-- make sure ni and nj match between tiles
724 if(ni .ne. nx*2 .OR. nj .ne. ny*2) &
725 call handle_err(-1, "mismatch of grid size between tiles")
726 endif
727
728 status=nf_inq_varid(ncid2, 'x', id_var)
729 call handle_err(status, 'inquire varid of x from file '//trim(grid_file) )
730 status=nf_get_var_double(ncid2, id_var, tmpvar)
731 call handle_err(status, 'inquire data of x from file '//trim(grid_file) )
732 if(t==7) then
733 geolon_c_nest(1:npx,1:npy) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
734 else
735 geolon_c(1:npx,1:npy,t) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
736 endif
737
738 status=nf_inq_varid(ncid2, 'y', id_var)
739 call handle_err(status, 'inquire varid of y from file '//trim(grid_file) )
740 status=nf_get_var_double(ncid2, id_var, tmpvar)
741 call handle_err(status, 'inquire data of y from file '//trim(grid_file) )
742 if(t==7) then
743 geolat_c_nest(1:npx,1:npy) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
744 else
745 geolat_c(1:npx,1:npy,t) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
746 endif
747 status = nf_close(ncid2)
748 call handle_err(status, "close file "//trim(tile_file))
749 enddo
750
751 deallocate(tmpvar)
752
753 status = nf_close(ncid)
754 call handle_err(status, "close file "//trim(grid_file))
755
756 is = 1 ; ie = nx
757 js = 1 ; je = ny
758 isd=is-ng; ied=ie+ng
759 jsd=js-ng; jed=je+ng
760
761 if( .not. regional ) then
762 call fill_cubic_grid_halo(geolon_c, geolon_c, ng, 1, 1, 1, 1)
763 call fill_cubic_grid_halo(geolat_c, geolat_c, ng, 1, 1, 1, 1)
764 if(.not. nested) call fill_bgrid_scalar_corners(geolon_c, ng, npx, npy, isd, jsd, xdir)
765 if(.not. nested) call fill_bgrid_scalar_corners(geolat_c, ng, npx, npy, isd, jsd, ydir)
766 else
767 call fill_regional_halo(geolon_c, ng)
768 call fill_regional_halo(geolat_c, ng)
769 endif
770
771 !--- compute grid cell center
772 allocate(geolon_t(isd:ied,jsd:jed,ntiles), geolat_t(isd:ied,jsd:jed,ntiles))
773
774 geolon_t(:,:,:) = -1.e25
775 geolat_t(:,:,:) = -1.e25
776
777 do t = 1, ntiles
778!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,G2,G3,G4,G5)
779 do j=js,je ; do i=is,ie
780 g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
781 g2(1) = geolon_c(i+1,j,t); g2(2) = geolat_c(i+1,j,t)
782 g3(1) = geolon_c(i,j+1,t); g3(2) = geolat_c(i,j+1,t)
783 g4(1) = geolon_c(i+1,j+1,t); g4(2) = geolat_c(i+1,j+1,t)
784 call cell_center2(g1, g2, g3, g4, g5 )
785 geolon_t(i,j,t) = g5(1)
786 geolat_t(i,j,t) = g5(2)
787 enddo ; enddo
788!$OMP END PARALLEL DO
789 enddo
790
791 if( .not. regional ) then
792 call fill_cubic_grid_halo(geolon_t, geolon_t, ng, 0, 0, 1, 1)
793 call fill_cubic_grid_halo(geolat_t, geolat_t, ng, 0, 0, 1, 1)
794 if (.not. nested) call fill_agrid_scalar_corners(geolon_t, ng, npx, npy, isd, jsd, xdir)
795 if (.not. nested) call fill_agrid_scalar_corners(geolat_t, ng, npx, npy, isd, jsd, ydir)
796 endif
797
798 !--- compute dx, dy
799 allocate(dx(isd:ied,jsd:jed+1,ntiles))
800 allocate(dy(isd:ied+1,jsd:jed,ntiles))
801 do t = 1, ntiles
802!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,G2)
803 do j = js, je+1 ; do i = is, ie
804 g1(1) = geolon_c(i ,j,t)
805 g1(2) = geolat_c(i ,j,t)
806 g2(1) = geolon_c(i+1,j,t)
807 g2(2) = geolat_c(i+1,j,t)
808 dx(i,j,t) = great_circle_dist( g2, g1, radius )
809 enddo ; enddo
810!$OMP END PARALLEL DO
811 enddo
812 do t = 1, ntiles
813!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,G2)
814 do j = js, je
815 do i = is, ie+1
816 g1(1) = geolon_c(i,j, t)
817 g1(2) = geolat_c(i,j, t)
818 g2(1) = geolon_c(i,j+1,t)
819 g2(2) = geolat_c(i,j+1,t)
820 dy(i,j,t) = great_circle_dist( g2, g1, radius )
821 enddo
822 enddo
823!$OMP END PARALLEL DO
824 enddo
825
826 if( .not. regional ) then
827 !--- make sure it is consitent between tiles. The following maybe not necessary.
828 do t = 1, ntiles
829 if(mod(t,2) ==0) then ! tile 2 4 6
830 tw = t - 1
831 te = t + 2
832 if(te > ntiles) te = te - ntiles
833 dy(is, js:je,t) = dy(ie+1,js:je,tw) ! west boundary
834 dy(ie+1, js:je, t) = dx(ie:is:-1,js, te) ! east boundary
835 else
836 tw = t - 2
837 if( tw <= 0) tw = tw + ntiles
838 te = t + 1
839 dy(is, js:je, t) = dx(ie:is:-1, je+1, tw) ! west boundary
840 dy(ie+1, js:je,t) = dy(1,js:je,te) ! east boundary
841 endif
842 enddo
843
844 call fill_cubic_grid_halo(dx, dy, ng, 0, 1, 1, 1)
845 call fill_cubic_grid_halo(dy, dx, ng, 1, 0, 1, 1)
846
847 if (.not. nested) call fill_dgrid_xy_corners(dx, dy, ng, npx, npy, isd, jsd)
848 endif
849
850 !--- compute dxa and dya -----
851 allocate(dxa(isd:ied,jsd:jed,ntiles))
852 allocate(dya(isd:ied,jsd:jed,ntiles))
853 do t = 1, ntiles
854!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,G2,G3,G4)
855 do j=js,je ; do i=is,ie
856 g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
857 g2(1) = geolon_c(i,j+1,t); g2(2) = geolat_c(i,j+1,t)
858 call mid_pt_sphere(g1, g2, g3)
859 g1(1) = geolon_c(i+1,j,t); g1(2) = geolat_c(i+1,j,t)
860 g2(1) = geolon_c(i+1,j+1,t); g2(2) = geolat_c(i+1,j+1,t)
861 call mid_pt_sphere(g1, g2, g4)
862 dxa(i,j,t) = great_circle_dist( g4, g3, radius )
863 g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
864 g2(1) = geolon_c(i+1,j,t); g2(2) = geolat_c(i+1,j,t)
865 call mid_pt_sphere(g1, g2, g3)
866 g1(1) = geolon_c(i,j+1,t); g1(2) = geolat_c(i,j+1,t)
867 g2(1) = geolon_c(i+1,j+1,t); g2(2) = geolat_c(i+1,j+1,t)
868 call mid_pt_sphere(g1, g2, g4)
869 dya(i,j,t) = great_circle_dist( g4, g3, radius )
870 enddo; enddo
871!$OMP END PARALLEL DO
872 enddo
873
874 if( .not.regional ) then
875 call fill_cubic_grid_halo(dxa, dya, ng, 0, 0, 1, 1)
876 call fill_cubic_grid_halo(dya, dxa, ng, 0, 0, 1, 1)
877
878 if (.not. nested) call fill_agrid_xy_corners(dxa, dya, ng, npx, npy, isd, jsd)
879 endif
880
881 !--- compute dxc and dyc
882 allocate(dxc(isd:ied+1,jsd:jed,ntiles))
883 allocate(dyc(isd:ied,jsd:jed+1,ntiles))
884 do t = 1, ntiles
885
886!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,G2)
887 do j=jsd,jed
888 do i=isd+1,ied
889 g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
890 g2(1) = geolon_c(i-1,j,t); g2(2) = geolat_c(i-1,j,t)
891 dxc(i,j,t) = great_circle_dist(g1, g2, radius)
892 enddo
893 dxc(isd,j,t) = dxc(isd+1,j,t)
894 dxc(ied+1,j,t) = dxc(ied,j,t)
895 enddo
896!$OMP END PARALLEL DO
897
898!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,G2)
899 do j=jsd+1,jed
900 do i=isd,ied
901 g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
902 g2(1) = geolon_c(i,j-1,t); g2(2) = geolat_c(i,j-1,t)
903 dyc(i,j,t) = great_circle_dist(g1, g2, radius)
904 enddo
905 enddo
906!$OMP END PARALLEL DO
907
908!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I)
909 do i=isd,ied
910 dyc(i,jsd,t) = dyc(i,jsd+1,t)
911 dyc(i,jed+1,t) = dyc(i,jed,t)
912 end do
913!$OMP END PARALLEL DO
914 enddo
915
916 !--- compute area
917 allocate(area(isd:ied,jsd:jed,ntiles))
918 do t = 1, ntiles
919!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,p_lL,p_uL,p_lr,p_uR)
920 do j=js,je
921 do i=is,ie
922 p_ll(1) = geolon_c(i ,j ,t) ; p_ll(2) = geolat_c(i ,j ,t)
923 p_ul(1) = geolon_c(i ,j+1,t) ; p_ul(2) = geolat_c(i ,j+1,t)
924 p_lr(1) = geolon_c(i+1,j ,t) ; p_lr(2) = geolat_c(i+1,j ,t)
925 p_ur(1) = geolon_c(i+1,j+1,t) ; p_ur(2) = geolat_c(i+1,j+1,t)
926
927 ! Spherical Excess Formula
928 area(i,j,t) = get_area(p_ll, p_ul, p_lr, p_ur, radius)
929 enddo
930 enddo
931!$OMP END PARALLEL DO
932 enddo
933
934 if( .not.regional ) then
935 call fill_cubic_grid_halo(area, area, ng, 0, 0, 1, 1)
936 endif
937
938 da_min = minval(area(is:ie,js:je,:))
939
940 !--- compute sin_sg
941 allocate(sin_sg(4,isd:ied,jsd:jed,ntiles))
942 allocate(cos_sg(4,isd:ied,jsd:jed))
943 allocate(grid3(3, npx, npy))
944 cos_sg(:,:,:) = big_number
945 sin_sg(:,:,:,:) = tiny_number
946
947 ! 9---4---8
948 ! | |
949 ! 1 5 3
950 ! | |
951 ! 6---2---7
952 do t = 1, ntiles
953
954!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1)
955 do j=js,je+1
956 do i = is,ie+1
957 g1(1) = geolon_c(i,j,t)
958 g1(2) = geolat_c(i,j,t)
959 call latlon2xyz(g1, grid3(:,i,j))
960 enddo
961 enddo
962!$OMP END PARALLEL DO
963
964!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,G1,P1,P3)
965 do j=js,je
966 do i=is,ie
967 g1(1) = geolon_t(i,j,t); g1(2) = geolat_t(i,j,t)
968 call latlon2xyz(g1, p3) ! righ-hand system consistent with grid3
969 call mid_pt3_cart(grid3(1,i,j), grid3(1,i,j+1), p1)
970 cos_sg(1,i,j) = cos_angle( p1, p3, grid3(1,i,j+1) )
971 call mid_pt3_cart(grid3(1,i,j), grid3(1,i+1,j), p1)
972 cos_sg(2,i,j) = cos_angle( p1, grid3(1,i+1,j), p3 )
973 call mid_pt3_cart(grid3(1,i+1,j), grid3(1,i+1,j+1), p1)
974 cos_sg(3,i,j) = cos_angle( p1, p3, grid3(1,i+1,j) )
975 call mid_pt3_cart(grid3(1,i,j+1), grid3(1,i+1,j+1), p1)
976 cos_sg(4,i,j) = cos_angle( p1, grid3(1,i,j+1), p3 )
977 enddo
978 enddo
979!$OMP END PARALLEL DO
980
981 do ip=1,4
982!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J)
983 do j=js,je
984 do i=is,ie
985 sin_sg(ip,i,j,t) = min(1.0, sqrt( max(0., 1.-cos_sg(ip,i,j)**2) ) )
986 enddo
987 enddo
988!$OMP END PARALLEL DO
989 enddo
990 enddo
991
992 if( .not.regional ) then
993 do ip=1,4
994 call fill_cubic_grid_halo(sin_sg(ip,:,:,:), sin_sg(ip,:,:,:), ng, 0, 0, 1, 1)
995 enddo
996 endif
997
998 deallocate(cos_sg, grid3, geolon_c, geolat_c, geolon_t, geolat_t)
999
1000
1001 end subroutine read_grid_file
1002
1007 subroutine read_topo_file(regional)
1008
1009 logical,intent(in) :: regional ! Is this a run with a regional domain?
1010 integer :: fsize=65536
1011 integer :: status, ncid, id_var, ndim, dimsiz, nt
1012 character(len=256) :: tile_file
1013 character(len=32) :: text
1014 integer :: len, t, dims(2)
1015 real :: tmp(is:ie,js:je)
1016
1017 allocate(oro(isd:ied,jsd:jed,ntiles))
1018 allocate(mask(isd:ied,jsd:jed,ntiles))
1019 oro = -big_number
1020 mask = 0
1021
1022 !--- make sure topo_file suffix is not ".nc"
1023 len = len_trim(topo_file)
1024 if( index(topo_file, '.nc', back=.true.) == len-2) then
1025 call handle_err(-1, "remove .nc from namelist topo_file="//trim(topo_file) )
1026 endif
1027
1028 !--- loop through each tile file to get the orography
1029 do nt = 1, ntiles
1030
1031 if( regional ) then
1032 t = nt + 6 ! The single regional tile must be #7 for now.
1033 else
1034 t = nt
1035 endif
1036
1037 write(text, '(i1.1)' ) t
1038 tile_file = trim(topo_file)//'.tile'//trim(text)//'.nc'
1039 status=nf__open(trim(tile_file),nf_nowrite,fsize,ncid)
1040 call handle_err(status, 'Open file '//trim(tile_file) )
1041
1042 status=nf_inq_varid(ncid, topo_field, id_var)
1043 call handle_err(status, 'inquire varid of '//trim(topo_field)//' from file '//trim(tile_file) )
1044
1045 status = nf_inq_varndims(ncid, id_var, ndim)
1046 call handle_err(status, 'inquire ndims of '//trim(topo_field)//' from file '//trim(tile_file) )
1047
1048 if(ndim .NE. 2) call handle_err(-1, 'ndims of '//trim(topo_field)//' from file '// &
1049 trim(tile_file)//' should be 2')
1050
1051 ! get data dimension and should match grid file size
1052 status = nf_inq_vardimid(ncid, id_var,dims);
1053 call handle_err(status, 'inquire dimid of '//trim(topo_field)//' from file '//trim(tile_file) )
1054
1055 status = nf_inq_dimlen(ncid, dims(1), dimsiz)
1056 call handle_err(status, 'inquire first dimension length of '//trim(topo_field)//' from file '//trim(tile_file) )
1057 if(dimsiz .NE. nx) call handle_err(-1, "mismatch of lon dimension size between "// &
1058 trim(grid_file)//' and '//trim(tile_file) )
1059
1060 status = nf_inq_dimlen(ncid, dims(2), dimsiz)
1061 call handle_err(status, 'inquire second dimension length of '//trim(topo_field)//' from file '//trim(tile_file) )
1062
1063 if(dimsiz .NE. ny) call handle_err(-1, "mismatch of lat dimension size between "// &
1064 trim(grid_file)//' and '//trim(tile_file) )
1065
1066 status = nf_get_var_double(ncid, id_var, oro(is:ie,js:je,nt))
1067 call handle_err(status, 'get the value of '//trim(topo_field)//' from file '//trim(tile_file) )
1068
1069 status=nf_inq_varid(ncid, mask_field, id_var)
1070 call handle_err(status, 'inquire varid of '//trim(mask_field)//' from file '//trim(tile_file) )
1071
1072 status = nf_get_var_double(ncid, id_var, tmp)
1073 call handle_err(status, 'get the value of '//trim(mask_field)//' from file '//trim(tile_file) )
1074
1075 mask(is:ie,js:je,nt) = tmp
1076
1077 status = nf_close(ncid)
1078 call handle_err(status, "close file "//trim(tile_file))
1079 enddo
1080
1081 if( .not.regional ) then
1082 !--- update halo
1083 call fill_cubic_grid_halo(oro, oro, ng, 0, 0, 1, 1)
1084 call fill_cubic_grid_halo(mask, mask, ng, 0, 0, 1, 1)
1085 endif
1086
1087 if( regional ) then
1088 call fill_regional_halo(oro, ng)
1089 oro(:,:,:) = max(oro(:,:,:),0.)
1090 call fill_regional_halo(mask, ng)
1091 mask(:,:,:) = min(max(mask(:,:,:),0.),1.)
1092 endif
1093
1094
1095 end subroutine read_topo_file
1096
1107 subroutine write_topo_file(is,ie,js,je,ntiles,q,regional)
1108 integer, intent(in) :: is,ie,js,je,ntiles
1109 real, intent(in) :: q(is:ie,js:je,ntiles)
1110 logical, intent(in) :: regional
1111
1112 integer :: fsize=65536
1113 integer :: nt, t, status, ncid, id_var
1114 character(len=256) :: tile_file
1115 character(len=3) :: text
1116 !--- loop through each tile file to update topo_field
1117
1118 do nt = 1, ntiles
1119
1120 if( regional ) then
1121 t = nt + 6
1122 else
1123 t = nt
1124 endif
1125
1126 write(text, '(i1.1)' ) t
1127 tile_file = trim(topo_file)//'.tile'//trim(text)//'.nc'
1128 status=nf__open(trim(tile_file),nf_write,fsize,ncid)
1129 call handle_err(status, 'write_topo_file: Open file '//trim(tile_file) )
1130
1131 status=nf_inq_varid(ncid, topo_field, id_var)
1132 call handle_err(status, 'write_topo_file:inquire varid of '//trim(topo_field)//' from file '//trim(tile_file) )
1133
1134 status = nf_put_var_double(ncid, id_var, q(:,:,nt))
1135 call handle_err(status, 'write_topo_file: put the value of '//trim(topo_field)//' from file '//trim(tile_file) )
1136
1137 status = nf_close(ncid)
1138 call handle_err(status, "write_topo_file: close file "//trim(tile_file))
1139 enddo
1140
1141
1142 end subroutine write_topo_file
1143
1155 subroutine fill_cubic_grid_halo(data, data2, halo, ioff, joff, sign1, sign2)
1156 integer, intent(in) :: halo
1157 real, dimension(1-halo:,1-halo:,:), intent(inout) :: data, data2
1158 integer, intent(in) :: ioff, joff, sign1, sign2
1159 integer :: lw, le, ls, ln
1160 integer :: i, tile
1161
1162 ntiles = size(data,3)
1163
1164 do tile = 1, ntiles
1165 if(mod(tile,2) == 0) then ! tile 2, 4, 6
1166 lw = tile - 1; le = tile + 2; ls = tile - 2; ln = tile + 1
1167 if(le > 6 ) le = le - 6
1168 if(ls < 1 ) ls = ls + 6
1169 if(ln > 6 ) ln = ln - 6
1170 data(1-halo:0, 1:ny+joff, tile) = data(nx-halo+1:nx, 1:ny+joff, lw) ! west
1171 do i = 1, halo
1172 data(nx+i+ioff, 1:ny+joff, tile) = sign1*data2(nx+joff:1:-1, i+ioff, le) ! east
1173 end do
1174 do i = 1, halo
1175 data(1:nx+ioff, 1-i, tile) = sign2*data2(nx-i+1, ny+ioff:1:-1, ls) ! south
1176 end do
1177 data(1:nx+ioff, ny+1+joff:ny+halo+joff, tile) = data(1:nx+ioff, 1+joff:halo+joff, ln) ! north
1178 else ! tile 1, 3, 5
1179 lw = tile - 2; le = tile + 1; ls = tile - 1; ln = tile + 2
1180 if(lw < 1 ) lw = lw + 6
1181 if(ls < 1 ) ls = ls + 6
1182 if(ln > 6 ) ln = ln - 6
1183 do i = 1, halo
1184 data(1-i, 1:ny+joff, tile) = sign1*data2(nx+joff:1:-1, ny-i+1, lw) ! west
1185 end do
1186 data(nx+1+ioff:nx+halo+ioff, 1:ny+joff, tile) = data(1+ioff:halo+ioff, 1:ny+joff, le) ! east
1187 data(1:nx+ioff, 1-halo:0, tile) = data(1:nx+ioff, ny-halo+1:ny, ls) ! south
1188 do i = 1, halo
1189 data(1:nx+ioff, ny+i+joff, tile) = sign2*data2(i+joff, ny+ioff:1:-1, ln) ! north
1190 end do
1191 end if
1192 enddo
1193
1194 end subroutine fill_cubic_grid_halo
1195
1224 subroutine fv3_zs_filter (is, ie, js, je, isd, ied, jsd, jed, npx, npy, npx_global, ntiles, &
1225 grid_type, stretch_fac, nested, area, dxa, dya, dx, dy, dxc, dyc, &
1226 sin_sg, phis, regional )
1227 integer, intent(in) :: is, ie, js, je, ntiles
1228 integer, intent(in) :: isd, ied, jsd, jed, npx, npy, npx_global, grid_type
1229 real, intent(in), dimension(isd:ied,jsd:jed, ntiles)::area, dxa, dya
1230 real, intent(in), dimension(isd:ied, jsd:jed+1, ntiles):: dx, dyc
1231 real, intent(in), dimension(isd:ied+1,jsd:jed, ntiles):: dy, dxc
1232
1233 real, intent(IN):: sin_sg(4,isd:ied,jsd:jed,ntiles)
1234 real, intent(IN):: stretch_fac
1235 logical, intent(IN) :: nested, regional
1236 real, intent(inout):: phis(isd:ied,jsd:jed,ntiles)
1237 real:: cd2
1238 integer mdim, n_del2, n_del4
1239
1240 mdim = nint( real(npx_global) * min(10., stretch_fac) )
1241
1242 ! Del-2: high resolution only
1243 if ( npx_global<=97 ) then
1244 n_del2 = 0
1245 elseif ( npx_global<=193 ) then
1246 n_del2 = 1
1247 else
1248 n_del2 = 2
1249 endif
1250 cd2 = 0.16*da_min
1251 ! Applying strong 2-delta-filter:
1252 if ( n_del2 > 0 ) &
1253 call two_delta_filter(is,ie,js,je,isd,ied,jsd,jed, npx, npy, ntiles, phis, area, &
1254 dx, dy, dxa, dya, dxc, dyc, sin_sg, cd2, zero_ocean, &
1255 .true.,0, grid_type, mask, nested, n_del2, regional)
1256
1257 ! MFCT Del-4:
1258 if ( mdim<=193 ) then
1259 n_del4 = 1
1260 elseif ( mdim<=1537 ) then
1261 n_del4 = 2
1262 else
1263 n_del4 = 3
1264 endif
1265 call del4_cubed_sphere(is,ie,js,je,isd,ied,jsd,jed,npx, npy, ntiles, &
1266 phis, area, dx, dy, dxc, dyc, sin_sg, n_del4, zero_ocean, mask, nested, regional)
1267 ! Applying weak 2-delta-filter:
1268 cd2 = 0.12*da_min
1269 call two_delta_filter(is,ie,js,je,isd,ied,jsd,jed,npx, npy, ntiles, &
1270 phis, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd2, zero_ocean, &
1271 .true., 1, grid_type, mask, nested, n_del2_weak, regional)
1272
1273 end subroutine fv3_zs_filter
1274
1307 subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, &
1308 q, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd, zero_ocean, &
1309 check_slope, filter_type, grid_type, mask, nested, ntmax, regional)
1310 integer, intent(in) :: is, ie, js, je
1311 integer, intent(in) :: isd, ied, jsd, jed
1312 integer, intent(in) :: npx, npy, grid_type
1313 integer, intent(in) :: ntmax, ntiles
1314 integer, intent(in) :: filter_type ! 0: strong, 1: weak
1315 real, intent(in) :: cd
1316 ! INPUT arrays
1317 real, intent(in)::area(isd:ied, jsd:jed, ntiles)
1318 real, intent(in):: dx(isd:ied, jsd:jed+1, ntiles)
1319 real, intent(in):: dy(isd:ied+1,jsd:jed, ntiles)
1320 real, intent(in):: dxa(isd:ied, jsd:jed, ntiles)
1321 real, intent(in):: dya(isd:ied, jsd:jed, ntiles)
1322 real, intent(in):: dxc(isd:ied+1,jsd:jed, ntiles)
1323 real, intent(in):: dyc(isd:ied, jsd:jed+1, ntiles)
1324 real, intent(in):: sin_sg(4,isd:ied,jsd:jed, ntiles)
1325 real, intent(in):: mask(isd:ied, jsd:jed, ntiles) ! 0==water, 1==land
1326 logical, intent(in):: zero_ocean, check_slope
1327 logical, intent(in):: nested, regional
1328 ! OUTPUT arrays
1329 real, intent(inout):: q(isd:ied, jsd:jed,ntiles)
1330 ! Local:
1331 real, parameter:: p1 = 7./12.
1332 real, parameter:: p2 = -1./12.
1333 real, parameter:: c1 = -2./14.
1334 real, parameter:: c2 = 11./14.
1335 real, parameter:: c3 = 5./14.
1336
1337 real:: ddx(is:ie+1,js:je), ddy(is:ie,js:je+1)
1338 logical:: extm(is-1:ie+1)
1339 logical:: ext2(is:ie,js-1:je+1)
1340 real:: a1(is-1:ie+2)
1341 real:: a2(is:ie,js-1:je+2)
1342 real:: a3(is:ie,js:je,ntiles)
1343 real:: smax, m_slope, fac
1344 integer:: i,j, nt, t
1345 integer:: is1, ie2, js1, je2
1346
1347 if ( .not. nested .and. grid_type<3 ) then
1348 is1 = max(3,is-1); ie2 = min(npx-2,ie+2)
1349 js1 = max(3,js-1); je2 = min(npy-2,je+2)
1350 else
1351 is1 = is-1; ie2 = ie+2
1352 js1 = js-1; je2 = je+2
1353 end if
1354
1355 if ( check_slope ) then
1356 m_slope = max_slope
1357 else
1358 m_slope = 10.
1359 endif
1360
1361
1362 do nt=1, ntmax
1363 if( .not.regional ) then
1364 call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1365 endif
1366
1367 ! Check slope
1368 if ( nt==1 .and. check_slope ) then
1369 do t = 1, ntiles
1370 do j=js,je
1371 do i=is,ie+1
1372 ddx(i,j) = (q(i,j,t) - q(i-1,j,t))/dxc(i,j,t)
1373 ddx(i,j) = abs(ddx(i,j))
1374 enddo
1375 enddo
1376 do j=js,je+1
1377 do i=is,ie
1378 ddy(i,j) = (q(i,j,t) - q(i,j-1,t))/dyc(i,j,t)
1379 ddy(i,j) = abs(ddy(i,j))
1380 enddo
1381 enddo
1382 do j=js,je
1383 do i=is,ie
1384 a3(i,j,t) = max( ddx(i,j), ddx(i+1,j), ddy(i,j), ddy(i,j+1) )
1385 enddo
1386 enddo
1387 enddo
1388 smax = maxval(a3(is:ie,js:je,:))
1389 write(*,*) 'Before filter: Max_slope=', smax
1390 endif
1391
1392
1393 ! First step: average the corners:
1394 if ( .not. nested .and. nt==1 ) then
1395 do t = 1, ntiles
1396 q(1,1,t) = (q(1,1,t)*area(1,1,t)+q(0,1,t)*area(0,1,t)+q(1,0,t)*area(1,0,t)) &
1397 / ( area(1,1,t)+ area(0,1,t)+ area(1,0,t) )
1398 q(0,1,t) = q(1,1,t)
1399 q(1,0,t) = q(1,1,t)
1400
1401 q(ie, 1,t) = (q(ie,1,t)*area(ie,1,t)+q(npx,1,t)*area(npx,1,t)+q(ie,0,t)*area(ie,0,t)) &
1402 / ( area(ie,1,t)+ area(npx,1,t)+ area(ie,0,t))
1403 q(npx,1,t) = q(ie,1,t)
1404 q(ie, 0,t) = q(ie,1,t)
1405
1406 q(1, je,t) = (q(1,je,t)*area(1,je,t)+q(0,je,t)*area(0,je,t)+q(1,npy,t)*area(1,npy,t)) &
1407 / ( area(1,je,t)+ area(0,je,t)+ area(1,npy,t))
1408 q(0, je,t) = q(1,je,t)
1409 q(1,npy,t) = q(1,je,t)
1410
1411 q(ie, je,t) = (q(ie,je,t)*area(ie,je,t)+q(npx,je,t)*area(npx,je,t)+q(ie,npy,t)*area(ie,npy,t)) &
1412 / ( area(ie,je,t)+ area(npx,je,t)+ area(ie,npy,t))
1413 q(npx,je,t) = q(ie,je,t)
1414 q(ie,npy,t) = q(ie,je,t)
1415 enddo
1416 if( .not.regional ) then
1417 call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1418 endif
1419 endif
1420
1421 do t = 1, ntiles
1422 ! x-diffusive flux:
1423 do j=js,je
1424
1425 do i=is1, ie2
1426 a1(i) = p1*(q(i-1,j,t)+q(i,j,t)) + p2*(q(i-2,j,t)+q(i+1,j,t))
1427 enddo
1428
1429 if ( (.not. (nested .or. regional)) .and. grid_type<3 ) then
1430 a1(0) = c1*q(-2,j,t) + c2*q(-1,j,t) + c3*q(0,j,t)
1431 a1(1) = 0.5*(((2.*dxa(0,j,t)+dxa(-1,j,t))*q(0,j,t)-dxa(0,j,t)*q(-1,j,t))/(dxa(-1,j,t)+dxa(0,j,t)) &
1432 + ((2.*dxa(1,j,t)+dxa( 2,j,t))*q(1,j,t)-dxa(1,j,t)*q( 2,j,t))/(dxa(1, j,t)+dxa(2,j,t)))
1433 a1(2) = c3*q(1,j,t) + c2*q(2,j,t) +c1*q(3,j,t)
1434
1435 a1(npx-1) = c1*q(npx-3,j,t) + c2*q(npx-2,j,t) + c3*q(npx-1,j,t)
1436 a1(npx) = 0.5*(((2.*dxa(npx-1,j,t)+dxa(npx-2,j,t))*q(npx-1,j,t)-dxa(npx-1,j,t)*q(npx-2,j,t)) &
1437 /(dxa(npx-2,j,t)+dxa(npx-1,j,t)) &
1438 + ((2.*dxa(npx, j,t)+dxa(npx+1,j,t))*q(npx, j,t)-dxa(npx, j,t)*q(npx+1,j,t))/ &
1439 (dxa(npx, j,t)+dxa(npx+1,j,t)))
1440 a1(npx+1) = c3*q(npx,j,t) + c2*q(npx+1,j,t) + c1*q(npx+2,j,t)
1441 endif
1442
1443 if ( regional .and. grid_type<3 ) then
1444 a1(0) = c1*q(-1,j,t) + c2*q(-1,j,t) + c3*q(0,j,t)
1445 a1(2) = c3*q(1,j,t) + c2*q(2,j,t) + c1*q(3,j,t)
1446 a1(1) = 0.5*(a1(0) + a1(2))
1447
1448 a1(npx-1) = c1*q(npx-3,j,t) + c2*q(npx-2,j,t) + c3*q(npx-1,j,t)
1449 a1(npx+1) = c3*q(npx,j,t) + c2*q(npx+1,j,t) + c1*q(npx+2,j,t)
1450 a1(npx) = 0.5*(a1(npx-1)+a1(npx+1))
1451 endif
1452
1453 if ( filter_type == 0 ) then
1454 do i=is-1, ie+1
1455 if( abs(3.*(a1(i)+a1(i+1)-2.*q(i,j,t))) > abs(a1(i)-a1(i+1)) ) then
1456 extm(i) = .true.
1457 else
1458 extm(i) = .false.
1459 endif
1460 enddo
1461 else
1462 do i=is-1, ie+1
1463 if ( (a1(i)-q(i,j,t))*(a1(i+1)-q(i,j,t)) > 0. ) then
1464 extm(i) = .true.
1465 else
1466 extm(i) = .false.
1467 endif
1468 enddo
1469 endif
1470
1471 do i=is,ie+1
1472 ddx(i,j) = (q(i-1,j,t)-q(i,j,t))/dxc(i,j,t)
1473 if ( extm(i-1).and.extm(i) ) then
1474 ddx(i,j) = 0.5*(sin_sg(3,i-1,j,t)+sin_sg(1,i,j,t))*dy(i,j,t)*ddx(i,j)
1475 elseif ( abs(ddx(i,j)) > m_slope ) then
1476 fac = min(1., max( 0.1, ( abs(ddx(i,j))-m_slope )/m_slope) )
1477 ddx(i,j) = fac*0.5*(sin_sg(3,i-1,j,t)+sin_sg(1,i,j,t))*dy(i,j,t)*ddx(i,j)
1478 else
1479 ddx(i,j) = 0.
1480 endif
1481 enddo
1482 enddo ! do j=js,je
1483
1484 ! y-diffusive flux:
1485 do j=js1,je2
1486 do i=is,ie
1487 a2(i,j) = p1*(q(i,j-1,t)+q(i,j,t)) + p2*(q(i,j-2,t)+q(i,j+1,t))
1488 enddo
1489 enddo
1490 if ( (.not. (nested .or. regional)) .and. grid_type<3 ) then
1491 do i=is,ie
1492 a2(i,0) = c1*q(i,-2,t) + c2*q(i,-1,t) + c3*q(i,0,t)
1493 a2(i,1) = 0.5*(((2.*dya(i,0,t)+dya(i,-1,t))*q(i,0,t)-dya(i,0,t)*q(i,-1,t))/(dya(i,-1,t)+dya(i,0,t)) &
1494 + ((2.*dya(i,1,t)+dya(i, 2,t))*q(i,1,t)-dya(i,1,t)*q(i, 2,t))/(dya(i, 1,t)+dya(i,2,t)))
1495 a2(i,2) = c3*q(i,1,t) + c2*q(i,2,t) + c1*q(i,3,t)
1496 enddo
1497
1498 do i=is,ie
1499 a2(i,npy-1) = c1*q(i,npy-3,t) + c2*q(i,npy-2,t) + c3*q(i,npy-1,t)
1500 a2(i,npy) = 0.5*(((2.*dya(i,npy-1,t)+dya(i,npy-2,t))*q(i,npy-1,t)-dya(i,npy-1,t)*q(i,npy-2,t))/ &
1501 (dya(i,npy-2,t)+dya(i,npy-1,t)) &
1502 + ((2.*dya(i,npy,t)+dya(i,npy+1,t))*q(i,npy,t)-dya(i,npy,t)*q(i,npy+1,t))/&
1503 (dya(i,npy,t)+dya(i,npy+1,t)))
1504 a2(i,npy+1) = c3*q(i,npy,t) + c2*q(i,npy+1,t) + c1*q(i,npy+2,t)
1505 enddo
1506 endif
1507
1508 if ( regional .and. grid_type<3 ) then
1509 do i=is,ie
1510 a2(i,0) = c1*q(i,-2,t) + c2*q(i,-1,t) + c3*q(i,0,t)
1511 a2(i,2) = c3*q(i,1,t) + c2*q(i,2,t) + c1*q(i,3,t)
1512 a2(i,1) = 0.5*(a2(i,0) + a2(i,2))
1513 enddo
1514
1515 do i=is,ie
1516 a2(i,npy-1) = c1*q(i,npy-3,t) + c2*q(i,npy-2,t) + c3*q(i,npy-1,t)
1517 a2(i,npy+1) = c3*q(i,npy,t) + c2*q(i,npy+1,t) + c1*q(i,npy+2,t)
1518 a2(i,npy) = 0.5*(a2(i,npy-1)+a2(i,npy+1))
1519 enddo
1520 endif
1521
1522 if ( filter_type == 0 ) then
1523 do j=js-1,je+1
1524 do i=is,ie
1525 if( abs(3.*(a2(i,j)+a2(i,j+1)-2.*q(i,j,t))) > abs(a2(i,j)-a2(i,j+1)) ) then
1526 ext2(i,j) = .true.
1527 else
1528 ext2(i,j) = .false.
1529 endif
1530 enddo
1531 enddo
1532 else
1533 do j=js-1,je+1
1534 do i=is,ie
1535 if ( (a2(i,j)-q(i,j,t))*(a2(i,j+1)-q(i,j,t)) > 0. ) then
1536 ext2(i,j) = .true.
1537 else
1538 ext2(i,j) = .false.
1539 endif
1540 enddo
1541 enddo
1542 endif
1543
1544 do j=js,je+1
1545 do i=is,ie
1546 ddy(i,j) = (q(i,j-1,t)-q(i,j,t))/dyc(i,j,t)
1547 if ( ext2(i,j-1) .and. ext2(i,j) ) then
1548 ddy(i,j) = 0.5*(sin_sg(4,i,j-1,t)+sin_sg(2,i,j,t))*dx(i,j,t)*ddy(i,j)
1549 elseif ( abs(ddy(i,j))>m_slope ) then
1550 fac = min(1., max(0.1,(abs(ddy(i,j))-m_slope)/m_slope))
1551 ddy(i,j) = fac*0.5*(sin_sg(4,i,j-1,t)+sin_sg(2,i,j,t))*dx(i,j,t)*ddy(i,j)
1552 else
1553 ddy(i,j) = 0.
1554 endif
1555 enddo
1556 enddo
1557
1558 if ( zero_ocean ) then
1559 ! Limit diffusive flux over water cells:
1560 do j=js,je
1561 do i=is,ie+1
1562 ddx(i,j) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * ddx(i,j)
1563 enddo
1564 enddo
1565 do j=js,je+1
1566 do i=is,ie
1567 ddy(i,j) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * ddy(i,j)
1568 enddo
1569 enddo
1570 endif
1571
1572 do j=js,je
1573 do i=is,ie
1574 q(i,j,t) = q(i,j,t) + cd/area(i,j,t)*(ddx(i,j)-ddx(i+1,j)+ddy(i,j)-ddy(i,j+1))
1575 enddo
1576 enddo
1577 enddo ! do t = 1, ntiles
1578 enddo ! nt=1, ntmax
1579
1580! Check slope
1581 if ( check_slope ) then
1582 if( .not.regional ) then
1583 call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1584 endif
1585 do t = 1, ntiles
1586 do j=js,je
1587 do i=is,ie+1
1588 ddx(i,j) = (q(i,j,t) - q(i-1,j,t))/dxc(i,j,t)
1589 ddx(i,j) = abs(ddx(i,j))
1590 enddo
1591 enddo
1592 do j=js,je+1
1593 do i=is,ie
1594 ddy(i,j) = (q(i,j,t) - q(i,j-1,t))/dyc(i,j,t)
1595 ddy(i,j) = abs(ddy(i,j))
1596 enddo
1597 enddo
1598 do j=js,je
1599 do i=is,ie
1600 a3(i,j,t) = max( ddx(i,j), ddx(i+1,j), ddy(i,j), ddy(i,j+1) )
1601 enddo
1602 enddo
1603 enddo
1604 smax = maxval(a3(is:ie,js:je,:))
1605 write(*,*) 'After filter: Max_slope=', smax
1606 endif
1607
1608 end subroutine two_delta_filter
1609
1637 subroutine del2_cubed_sphere(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles,&
1638 q, area, dx, dy, dxc, dyc, sin_sg, nmax, cd, zero_ocean, mask, nested, regional)
1639 integer, intent(in) :: is, ie, js, je
1640 integer, intent(in) :: isd, ied, jsd, jed
1641 integer, intent(in):: npx, npy, ntiles
1642 integer, intent(in):: nmax
1643 real, intent(in):: cd
1644 logical, intent(in):: zero_ocean
1645 ! INPUT arrays
1646 real, intent(in)::area(isd:ied, jsd:jed, ntiles)
1647 real, intent(in):: dx(isd:ied, jsd:jed+1, ntiles)
1648 real, intent(in):: dy(isd:ied+1,jsd:jed, ntiles)
1649 real, intent(in):: dxc(isd:ied+1,jsd:jed, ntiles)
1650 real, intent(in):: dyc(isd:ied, jsd:jed+1, ntiles)
1651 real, intent(IN):: sin_sg(4,isd:ied,jsd:jed, ntiles)
1652 real, intent(in):: mask(isd:ied, jsd:jed, ntiles) ! 0==water, 1==land
1653 logical, intent(IN) :: nested, regional
1654
1655 ! OUTPUT arrays
1656 real, intent(inout):: q(is-ng:ie+ng, js-ng:je+ng, ntiles)
1657 ! Local:
1658 real ddx(is:ie+1,js:je), ddy(is:ie,js:je+1)
1659 integer i,j,n,t
1660
1661 if( .not.regional ) then
1662 call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1663 endif
1664
1665 do t = 1, ntiles
1666 ! First step: average the corners:
1667 if ( .not. nested) then
1668 q(1,1,t) = (q(1,1,t)*area(1,1,t)+q(0,1,t)*area(0,1,t)+q(1,0,t)*area(1,0,t)) &
1669 / ( area(1,1,t)+ area(0,1,t)+ area(1,0,t) )
1670 q(0,1,t) = q(1,1,t)
1671 q(1,0,t) = q(1,1,t)
1672 endif
1673 if ( .not. nested) then
1674 q(ie, 1,t) = (q(ie,1,t)*area(ie,1,t)+q(npx,1,t)*area(npx,1,t)+q(ie,0,t)*area(ie,0,t)) &
1675 / ( area(ie,1,t)+ area(npx,1,t)+ area(ie,0,t))
1676 q(npx,1,t) = q(ie,1,t)
1677 q(ie, 0,t) = q(ie,1,t)
1678 endif
1679 if ( .not. nested ) then
1680 q(ie, je,t) = (q(ie,je,t)*area(ie,je,t)+q(npx,je,t)*area(npx,je,t)+q(ie,npy,t)*area(ie,npy,t)) &
1681 / ( area(ie,je,t)+ area(npx,je,t)+ area(ie,npy,t))
1682 q(npx,je,t) = q(ie,je,t)
1683 q(ie,npy,t) = q(ie,je,t)
1684 endif
1685 if ( .not. nested) then
1686 q(1, je,t) = (q(1,je,t)*area(1,je,t)+q(0,je,t)*area(0,je,t)+q(1,npy,t)*area(1,npy,t)) &
1687 / ( area(1,je,t)+ area(0,je,t)+ area(1,npy,t))
1688 q(0, je,t) = q(1,je,t)
1689 q(1,npy,t) = q(1,je,t)
1690 endif
1691 enddo
1692
1693 do n=1,nmax
1694 if( .not.regional ) then
1695 if( n>1 ) call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1696 endif
1697 do t = 1, ntiles
1698 do j=js,je
1699 do i=is,ie+1
1700 ddx(i,j) = 0.5*(sin_sg(3,i-1,j,t)+sin_sg(1,i,j,t))*dy(i,j,t)*(q(i-1,j,t)-q(i,j,t))/dxc(i,j,t)
1701 enddo
1702 enddo
1703 do j=js,je+1
1704 do i=is,ie
1705 ddy(i,j) = dx(i,j,t)*(q(i,j-1,t)-q(i,j,t))/dyc(i,j,t) &
1706 *0.5*(sin_sg(4,i,j-1,t)+sin_sg(2,i,j,t))
1707 enddo
1708 enddo
1709
1710 if ( zero_ocean ) then
1711 ! Limit diffusive flux over ater cells:
1712 do j=js,je
1713 do i=is,ie+1
1714 ddx(i,j) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * ddx(i,j)
1715 enddo
1716 enddo
1717 do j=js,je+1
1718 do i=is,ie
1719 ddy(i,j) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * ddy(i,j)
1720 enddo
1721 enddo
1722 endif
1723
1724 do j=js,je
1725 do i=is,ie
1726 q(i,j,t) = q(i,j,t) + cd/area(i,j,t)*(ddx(i,j)-ddx(i+1,j)+ddy(i,j)-ddy(i,j+1))
1727 enddo
1728 enddo
1729 enddo
1730 enddo
1731
1732 end subroutine del2_cubed_sphere
1733
1760 subroutine del4_cubed_sphere(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, &
1761 q, area, dx, dy, dxc, dyc, sin_sg, nmax, zero_ocean, mask, nested, regional)
1762 integer, intent(in) :: is, ie, js, je
1763 integer, intent(in) :: isd, ied, jsd, jed
1764 integer, intent(in) :: npx, npy, nmax, ntiles
1765 logical, intent(in) :: zero_ocean
1766 real, intent(in):: mask(isd:ied, jsd:jed, ntiles) ! 0==water, 1==land
1767 real, intent(in)::area(isd:ied, jsd:jed, ntiles)
1768 real, intent(in):: dx(isd:ied, jsd:jed+1, ntiles)
1769 real, intent(in):: dy(isd:ied+1,jsd:jed, ntiles)
1770 real, intent(in):: dxc(isd:ied+1,jsd:jed, ntiles)
1771 real, intent(in):: dyc(isd:ied, jsd:jed+1, ntiles)
1772 real, intent(IN):: sin_sg(4,isd:ied,jsd:jed, ntiles)
1773 real, intent(inout):: q(isd:ied, jsd:jed, ntiles)
1774 logical, intent(IN) :: nested, regional
1775 ! Local:
1776 ! diffusivity
1777 real :: diff(is-1:ie+1,js-1:je+1, ntiles)
1778 ! diffusive fluxes:
1779 real :: fx1(is:ie+1,js:je), fy1(is:ie,js:je+1)
1780 real :: fx2(is:ie+1,js:je,ntiles), fy2(is:ie,js:je+1,ntiles)
1781 real :: fx4(is:ie+1,js:je,ntiles), fy4(is:ie,js:je+1,ntiles)
1782 real, dimension(isd:ied,jsd:jed,ntiles):: d2, win, wou
1783 real, dimension(is:ie,js:je, ntiles) :: qlow, qmin, qmax
1784 real, parameter:: esl = 1.e-20
1785 integer i,j, n, t
1786
1787 ! On a nested grid the haloes are not filled. Set to zero.
1788 d2 = 0.
1789 win = 0.
1790 wou = 0.
1791
1792 do t = 1, ntiles
1793 do j=js-1,je+1 ; do i=is-1,ie+1
1794 diff(i,j,t) = cd4*area(i,j,t) ! area dependency is needed for stretched grid
1795 enddo; enddo
1796
1797 do j=js,je ; do i=is,ie
1798 qmax(i,j,t) = q(i,j,t) * peak_fac
1799 qmin(i,j,t) = q(i,j,t) / peak_fac
1800 enddo; enddo
1801 enddo
1802
1803 do n=1,nmax
1804 if( .not.regional ) then
1805 call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1806 endif
1807
1808 ! First step: average the corners:
1809 if ( .not. nested .and. n==1 ) then
1810 do t = 1, ntiles
1811 q(1,1,t) = (q(1,1,t)*area(1,1,t)+q(0,1,t)*area(0,1,t)+q(1,0,t)*area(1,0,t)) &
1812 / ( area(1,1,t)+ area(0,1,t)+ area(1,0,t) )
1813 q(0,1,t) = q(1,1,t)
1814 q(1,0,t) = q(1,1,t)
1815 q(0,0,t) = q(1,1,t)
1816
1817 q(ie, 1,t) = (q(ie,1,t)*area(ie,1,t)+q(npx,1,t)*area(npx,1,t)+q(ie,0,t)*area(ie,0,t)) &
1818 / ( area(ie,1,t)+ area(npx,1,t)+ area(ie,0,t))
1819 q(npx,1,t) = q(ie,1,t)
1820 q(ie, 0,t) = q(ie,1,t)
1821 q(npx,0,t) = q(ie,1,t)
1822
1823 q(1, je,t) = (q(1,je,t)*area(1,je,t)+q(0,je,t)*area(0,je,t)+q(1,npy,t)*area(1,npy,t)) &
1824 / ( area(1,je,t)+ area(0,je,t)+ area(1,npy,t))
1825 q(0, je,t) = q(1,je,t)
1826 q(1,npy,t) = q(1,je,t)
1827 q(0,npy,t) = q(1,je,t)
1828
1829 q(ie, je,t) = (q(ie,je,t)*area(ie,je,t)+q(npx,je,t)*area(npx,je,t)+q(ie,npy,t)*area(ie,npy,t)) &
1830 / ( area(ie,je,t)+ area(npx,je,t)+ area(ie,npy,t))
1831 q(npx, je,t) = q(ie,je,t)
1832 q(ie, npy,t) = q(ie,je,t)
1833 q(npx,npy,t) = q(ie,je,t)
1834 enddo
1835 if( .not.regional ) then
1836 call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1837 endif
1838 endif
1839
1840 do t = 1, ntiles
1841
1842 !--------------
1843 ! Compute del-2
1844 !--------------
1845 ! call copy_corners(q, npx, npy, 1)
1846 do j=js,je
1847 do i=is,ie+1
1848 fx2(i,j,t) = 0.25*(diff(i-1,j,t)+diff(i,j,t))*dy(i,j,t)*(q(i-1,j,t)-q(i,j,t))/dxc(i,j,t) &
1849 *(sin_sg(1,i,j,t)+sin_sg(3,i-1,j,t))
1850 enddo
1851 enddo
1852
1853 ! call copy_corners(q, npx, npy, 2)
1854 do j=js,je+1
1855 do i=is,ie
1856 fy2(i,j,t) = 0.25*(diff(i,j-1,t)+diff(i,j,t))*dx(i,j,t)*(q(i,j-1,t)-q(i,j,t))/dyc(i,j,t) &
1857 *(sin_sg(2,i,j,t)+sin_sg(4,i,j-1,t))
1858 enddo
1859 enddo
1860
1861 do j=js,je
1862 do i=is,ie
1863 d2(i,j,t) = (fx2(i,j,t)-fx2(i+1,j,t)+fy2(i,j,t)-fy2(i,j+1,t)) / area(i,j,t)
1864 enddo
1865 enddo
1866
1867 ! qlow == low order monotonic solution
1868 if ( zero_ocean ) then
1869 ! Limit diffusive flux over water cells:
1870 do j=js,je
1871 do i=is,ie+1
1872 fx1(i,j) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * fx2(i,j,t)
1873 enddo
1874 enddo
1875 do j=js,je+1
1876 do i=is,ie
1877 fy1(i,j) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * fy2(i,j,t)
1878 enddo
1879 enddo
1880 do j=js,je
1881 do i=is,ie
1882 qlow(i,j,t) = q(i,j,t) + (fx1(i,j)-fx1(i+1,j)+fy1(i,j)-fy1(i,j+1)) / area(i,j,t)
1883 d2(i,j,t) = diff(i,j,t) * d2(i,j,t)
1884 enddo
1885 enddo
1886 else
1887 do j=js,je
1888 do i=is,ie
1889 qlow(i,j,t) = q(i,j,t) + d2(i,j,t)
1890 d2(i,j,t) = diff(i,j,t) * d2(i,j,t)
1891 enddo
1892 enddo
1893 endif
1894 enddo
1895 if( .not.regional ) then
1896 call fill_cubic_grid_halo(d2, d2, ng, 0, 0, 1, 1)
1897 endif
1898
1899 !---------------------
1900 ! Compute del4 fluxes:
1901 !---------------------
1902 ! call copy_corners(d2, npx, npy, 1)
1903 do t = 1, ntiles
1904 do j=js,je
1905 do i=is,ie+1
1906 fx4(i,j,t) = 0.5*(sin_sg(3,i-1,j,t)+sin_sg(1,i,j,t))*dy(i,j,t)*(d2(i,j,t)-d2(i-1,j,t))/dxc(i,j,t)-fx2(i,j,t)
1907 enddo
1908 enddo
1909
1910 ! call copy_corners(d2, npx, npy, 2)
1911 do j=js,je+1
1912 do i=is,ie
1913 fy4(i,j,t) = dx(i,j,t)*(d2(i,j,t)-d2(i,j-1,t))/dyc(i,j,t) &
1914 *0.5*(sin_sg(2,i,j,t)+sin_sg(4,i,j-1,t))-fy2(i,j,t)
1915 enddo
1916 enddo
1917
1918 do j=js,je
1919 do i=is,ie
1920 qmin(i,j,t) = min(qmin(i,j,t), q(i-1,j-1,t), q(i,j-1,t), q(i+1,j-1,t), &
1921 q(i-1,j ,t), q(i,j ,t), q(i+1,j ,t), &
1922 q(i-1,j+1,t), q(i,j+1,t), q(i+1,j+1,t) )
1923 qmax(i,j,t) = max(qmax(i,j,t), q(i-1,j-1,t), q(i,j-1,t), q(i+1,j-1,t), &
1924 q(i-1,j ,t), q(i,j ,t), q(i+1,j ,t), &
1925 q(i-1,j+1,t), q(i,j+1,t), q(i+1,j+1,t) )
1926 enddo
1927 enddo
1928
1929 !----------------
1930 ! Flux limitting:
1931 !----------------
1932 do j=js,je
1933 do i=is,ie
1934 win(i,j,t) = max(0.,fx4(i, j,t)) - min(0.,fx4(i+1,j,t)) + &
1935 max(0.,fy4(i, j,t)) - min(0.,fy4(i,j+1,t)) + esl
1936 wou(i,j,t) = max(0.,fx4(i+1,j,t)) - min(0.,fx4(i, j,t)) + &
1937 max(0.,fy4(i,j+1,t)) - min(0.,fy4(i, j,t)) + esl
1938 win(i,j,t) = max(0., qmax(i,j,t) - qlow(i,j,t)) / win(i,j,t)*area(i,j,t)
1939 wou(i,j,t) = max(0., qlow(i,j,t) - qmin(i,j,t)) / wou(i,j,t)*area(i,j,t)
1940 enddo
1941 enddo
1942 enddo
1943 if( .not.regional ) then
1944 call fill_cubic_grid_halo(win, win, ng, 0, 0, 1, 1)
1945 call fill_cubic_grid_halo(wou, wou, ng, 0, 0, 1, 1)
1946 endif
1947 do t = 1, ntiles
1948 do j=js,je
1949 do i=is,ie+1
1950 if ( fx4(i,j,t) > 0. ) then
1951 fx4(i,j,t) = min(1., wou(i-1,j,t), win(i,j,t)) * fx4(i,j,t)
1952 else
1953 fx4(i,j,t) = min(1., win(i-1,j,t), wou(i,j,t)) * fx4(i,j,t)
1954 endif
1955 enddo
1956 enddo
1957 do j=js,je+1
1958 do i=is,ie
1959 if ( fy4(i,j,t) > 0. ) then
1960 fy4(i,j,t) = min(1., wou(i,j-1,t), win(i,j,t)) * fy4(i,j,t)
1961 else
1962 fy4(i,j,t) = min(1., win(i,j-1,t), wou(i,j,t)) * fy4(i,j,t)
1963 endif
1964 enddo
1965 enddo
1966
1967
1968 if ( zero_ocean ) then
1969 ! Limit diffusive flux over ocean cells:
1970 do j=js,je
1971 do i=is,ie+1
1972 fx4(i,j,t) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * fx4(i,j,t)
1973 enddo
1974 enddo
1975 do j=js,je+1
1976 do i=is,ie
1977 fy4(i,j,t) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * fy4(i,j,t)
1978 enddo
1979 enddo
1980 endif
1981
1982 ! Update:
1983 do j=js,je
1984 do i=is,ie
1985 q(i,j,t) = qlow(i,j,t) + (fx4(i,j,t)-fx4(i+1,j,t)+fy4(i,j,t)-fy4(i,j+1,t))/area(i,j,t)
1986 enddo
1987 enddo
1988 enddo
1989 enddo ! end n-loop
1990
1991
1992 end subroutine del4_cubed_sphere
1993
1998 subroutine check(status)
1999 use netcdf
2000 integer,intent(in) :: status
2001!
2002 if(status /= nf90_noerr) then
2003 write(0,*) ' check netcdf status = ', status
2004 write(0,'("error ", a)') trim(nf90_strerror(status))
2005 write(0,*) "Stopped"
2006 stop 4
2007 endif
2008 end subroutine check
2009
2014 subroutine compute_filter_constants
2015
2016 ! set the given values for various cube resolutions (c48, c96, c192, c384, c768, c1152, c3072)
2017
2018 integer,parameter :: nres=8
2019 integer :: index1,index2,n
2020
2021 real :: factor
2022
2023 real,dimension(1:nres) :: cube_res=(/48.,96.,128.,192.,384.,768.,1152.,3072./)
2024
2025 real,dimension(1:nres) :: n_del2_weak_vals=(/4.,8.,8.,12.,12.,16.,20.,24./)
2026 real,dimension(1:nres) :: cd4_vals =(/0.12,0.12,0.13,0.15,0.15,0.15,0.15,0.15/)
2027 real,dimension(1:nres) :: max_slope_vals =(/0.12,0.12,0.12,0.12,0.12,0.12,0.16,0.30/)
2028 real,dimension(1:nres) :: peak_fac_vals =(/1.1,1.1,1.1,1.05,1.0,1.0,1.0,1.0/)
2029
2030 if(res<cube_res(1))then
2031 index1 = 1
2032 index2 = 1
2033 factor = 0.
2034 elseif(res>cube_res(nres))then
2035 index1 = nres
2036 index2 = nres
2037 factor = 0.
2038 else
2039 do n=2,nres
2040 if(res<=cube_res(n))then
2041 index2 = n
2042 index1 = n-1
2043 factor = (res-cube_res(n-1))/(cube_res(n)-cube_res(n-1))
2044 exit
2045 endif
2046 enddo
2047 endif
2048
2049 n_del2_weak = nint(n_del2_weak_vals(index1)+factor*(n_del2_weak_vals(index2)-n_del2_weak_vals(index1)))
2050 cd4 = cd4_vals(index1)+factor*(cd4_vals(index2)-cd4_vals(index1))
2051 max_slope = max_slope_vals(index1)+factor*(max_slope_vals(index2)-max_slope_vals(index1))
2052 peak_fac = peak_fac_vals(index1)+factor*(peak_fac_vals(index2)-peak_fac_vals(index1))
2053
2054 print*,''
2055 print*,'- FILTER COEFFICIENTS FOR RESOLUTION ', res
2056 print*,'- CD4: ', cd4
2057 print*,'- N_DEL2_WEAK: ', n_del2_weak
2058 print*,'- MAX_SLOPE: ', max_slope
2059 print*,'- PEAK_FAC: ', peak_fac
2060
2061 end subroutine compute_filter_constants
2062
2063end program filter_topo
mid_pt3_cart
subroutine mid_pt3_cart(p1, p2, e)
???
Definition filter_topo.F90:431
fill_bgrid_scalar_corners
subroutine fill_bgrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
???
Definition filter_topo.F90:302
del2_cubed_sphere
subroutine del2_cubed_sphere(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, q, area, dx, dy, dxc, dyc, sin_sg, nmax, cd, zero_ocean, mask, nested, regional)
???
Definition filter_topo.F90:1639
latlon2xyz
subroutine latlon2xyz(p, e)
???
Definition filter_topo.F90:464
cart_to_latlon
subroutine cart_to_latlon(np, q, xs, ys)
???
Definition filter_topo.F90:498
two_delta_filter
subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, q, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd, zero_ocean, check_slope, filter_type, grid_type, mask, nested, ntmax, regional)
???
Definition filter_topo.F90:1310
fill_agrid_xy_corners
subroutine fill_agrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
???
Definition filter_topo.F90:353
compute_filter_constants
subroutine compute_filter_constants
Compute resolution-dependent values for the filtering.
Definition filter_topo.F90:2015
fv3_zs_filter
subroutine fv3_zs_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, npx_global, ntiles, grid_type, stretch_fac, nested, area, dxa, dya, dx, dy, dxc, dyc, sin_sg, phis, regional)
???
Definition filter_topo.F90:1227
check
subroutine check(status)
Check results of netCDF call.
Definition filter_topo.F90:1999
read_grid_file
subroutine read_grid_file(regional)
???
Definition filter_topo.F90:630
del4_cubed_sphere
subroutine del4_cubed_sphere(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, q, area, dx, dy, dxc, dyc, sin_sg, nmax, zero_ocean, mask, nested, regional)
???
Definition filter_topo.F90:1762
filter_topo
program filter_topo
This program does ???
Definition filter_topo.F90:10
spherical_angle
real function spherical_angle(p1, p2, p3)
???
Definition filter_topo.F90:138
get_area
real function get_area(p1, p4, p2, p3, radius)
???
Definition filter_topo.F90:199
fill_cubic_grid_halo
subroutine fill_cubic_grid_halo(data, data2, halo, ioff, joff, sign1, sign2)
This routine fill the halo points for the cubic grid.
Definition filter_topo.F90:1156
fill_dgrid_xy_corners
subroutine fill_dgrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
???
Definition filter_topo.F90:385
read_topo_file
subroutine read_topo_file(regional)
???
Definition filter_topo.F90:1008
great_circle_dist
real function great_circle_dist(q1, q2, radius)
???
Definition filter_topo.F90:97
write_topo_file
subroutine write_topo_file(is, ie, js, je, ntiles, q, regional)
Replace the topo_field.
Definition filter_topo.F90:1108
mid_pt_sphere
subroutine mid_pt_sphere(p1, p2, pm)
???
Definition filter_topo.F90:412
cos_angle
real function cos_angle(p1, p2, p3)
???
Definition filter_topo.F90:545
cell_center2
subroutine cell_center2(q1, q2, q3, q4, e2)
???
Definition filter_topo.F90:599
fill_agrid_scalar_corners
subroutine fill_agrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
???
Definition filter_topo.F90:252
utils
Module that contains general utility routines.
Definition utils.F90:8
utils::regional
logical regional
If true, process a stand-alone regional grid.
Definition utils.F90:23
utils::grid_type
integer grid_type
Grid type.
Definition utils.F90:25
utils::stretch_fac
real stretch_fac
Grid stretching factor.
Definition utils.F90:27
utils::read_namelist
subroutine read_namelist
Read the program namelist file.
Definition utils.F90:37
utils::nested
logical nested
If true, process a global grid with a nest.
Definition utils.F90:22