grid_tools  1.13.0
filter_topo.F90
Go to the documentation of this file.
1 
4 
10 program 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 
87 contains
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 
2063 end program filter_topo
cart_to_latlon
subroutine cart_to_latlon(np, q, xs, ys)
???
Definition: filter_topo.F90:498
cell_center2
subroutine cell_center2(q1, q2, q3, q4, e2)
???
Definition: filter_topo.F90:599
read_topo_file
subroutine read_topo_file(regional)
???
Definition: filter_topo.F90:1008
utils::fill_regional_halo
subroutine fill_regional_halo(data, halo)
This routine extrapolate geolat_c and geolon_c halo points for the regional standalone grid...
Definition: utils.F90:71
mid_pt3_cart
subroutine mid_pt3_cart(p1, p2, e)
???
Definition: filter_topo.F90:431
cos_angle
real function cos_angle(p1, p2, p3)
???
Definition: filter_topo.F90:545
fill_agrid_scalar_corners
subroutine fill_agrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
???
Definition: filter_topo.F90:252
latlon2xyz
subroutine latlon2xyz(p, e)
???
Definition: filter_topo.F90:464
great_circle_dist
real function great_circle_dist(q1, q2, radius)
???
Definition: filter_topo.F90:97
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
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
utils::nested
logical nested
If true, process a global grid with a nest.
Definition: utils.F90:22
utils::handle_err
subroutine handle_err(status, string)
Prints an error message to standard output, then halts program execution with a bad status...
Definition: utils.F90:107
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_bgrid_scalar_corners
subroutine fill_bgrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
???
Definition: filter_topo.F90:302
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
utils::read_namelist
subroutine read_namelist
Read the program namelist file.
Definition: utils.F90:37
utils::grid_file
character(len=512) grid_file
Path/name of the grid mosaic file.
Definition: utils.F90:18
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
fill_agrid_xy_corners
subroutine fill_agrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
???
Definition: filter_topo.F90:353
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
utils::stretch_fac
real stretch_fac
Grid stretching factor.
Definition: utils.F90:27
utils
Module that contains general utility routines.
Definition: utils.F90:8
fill_dgrid_xy_corners
subroutine fill_dgrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
???
Definition: filter_topo.F90:385
read_grid_file
subroutine read_grid_file(regional)
???
Definition: filter_topo.F90:630
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
check
subroutine check(status)
Check results of netCDF call.
Definition: filter_topo.F90:1999
compute_filter_constants
subroutine compute_filter_constants
Compute resolution-dependent values for the filtering.
Definition: filter_topo.F90:2015
filter_topo
program filter_topo
This program does ???
Definition: filter_topo.F90:10
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