grid_tools  1.13.0
filter_topo.F90
Go to the documentation of this file.
1 
4 
10 program filter_topo
11 
12  use utils
13 
14  implicit none
15 
16 #include <netcdf.inc>
17 
18 #ifdef NO_QUAD_PRECISION
19  ! 64-bit precision (kind=8)
20  integer, parameter:: f_p = selected_real_kind(15)
21 #else
22  ! Higher precision (kind=16) for grid geometrical factors:
23  integer, parameter:: f_p = selected_real_kind(20)
24 #endif
25 
26  integer, parameter :: xdir=1
27  integer, parameter :: ydir=2
28  real, parameter :: pi = 3.14159265358979323846d0
29  real, parameter :: radius = 6371200.0
30  real, parameter :: big_number=1.d8
31  real, parameter :: tiny_number=1.d-8
32 
33  real:: cd4 ! Dimensionless coeff for del-4 difussion (with FCT)
34  real:: peak_fac ! overshoot factor for the mountain peak
35  real:: max_slope ! max allowable terrain slope: 1 --> 45 deg
36 
37  integer :: n_del2_weak
38 
39  integer :: ntiles = 0
40 
41  real da_min
42 
43  real, allocatable :: oro(:,:,:), mask(:,:,:)
44  real, allocatable :: dx(:,:,:), dy(:,:,:)
45  real, allocatable :: dxa(:,:,:), dya(:,:,:)
46  real, allocatable :: dxc(:,:,:), dyc(:,:,:)
47  real, allocatable :: area(:,:,:)
48  real, allocatable :: sin_sg(:,:,:,:)
49 
50  integer :: is,ie,js,je,isd,ied,jsd,jed
51  integer,parameter :: ng = 3
52  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
53  !--- read namelist
54  call read_namelist()
55 
56  !--- compute filter constants according to grid resolution.
57  call compute_filter_constants
58 
59  !--- read the target grid.
60  call read_grid_file(regional)
61 
62  !--- read the topography data
63  call read_topo_file(regional)
64 
65  !--- filter the data
66  call fv3_zs_filter(is,ie,js,je,isd,ied,jsd,jed,npx,npy,npx,ntiles,grid_type, &
67  stretch_fac, nested, area, dxa, dya, dx, dy, dxc, dyc, sin_sg, oro, regional )
68 
69  !--- write out the data
70  call write_topo_file(is,ie,js,je,ntiles,oro(is:ie,js:je,:),regional )
71 
72 contains
73 
81  real function great_circle_dist( q1, q2, radius )
82  real, intent(IN) :: q1(2), q2(2)
83  real, intent(IN), optional :: radius
84 
85  real (f_p):: p1(2), p2(2)
86  real (f_p):: beta
87  integer n
88 
89  do n=1,2
90  p1(n) = q1(n)
91  p2(n) = q2(n)
92  enddo
93 
94  beta = asin( sqrt( sin((p1(2)-p2(2))/2.)**2 + cos(p1(2))*cos(p2(2))* &
95  sin((p1(1)-p2(1))/2.)**2 ) ) * 2.
96 
97  if ( present(radius) ) then
98  great_circle_dist = radius * beta
99  else
100  great_circle_dist = beta ! Returns the angle
101  endif
102 
103  end function great_circle_dist
104 
122  real function spherical_angle(p1, p2, p3)
123 
124  real p1(3), p2(3), p3(3)
125 
126  real (f_p):: e1(3), e2(3), e3(3)
127  real (f_p):: px, py, pz
128  real (f_p):: qx, qy, qz
129  real (f_p):: angle, ddd
130  integer n
131 
132  do n=1,3
133  e1(n) = p1(n)
134  e2(n) = p2(n)
135  e3(n) = p3(n)
136  enddo
137 
138  !-------------------------------------------------------------------
139  ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
140  !-------------------------------------------------------------------
141  ! Vector P:
142  px = e1(2)*e2(3) - e1(3)*e2(2)
143  py = e1(3)*e2(1) - e1(1)*e2(3)
144  pz = e1(1)*e2(2) - e1(2)*e2(1)
145  ! Vector Q:
146  qx = e1(2)*e3(3) - e1(3)*e3(2)
147  qy = e1(3)*e3(1) - e1(1)*e3(3)
148  qz = e1(1)*e3(2) - e1(2)*e3(1)
149 
150  ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz)
151 
152  if ( ddd <= 0.0d0 ) then
153  angle = 0.d0
154  else
155  ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd)
156  if ( abs(ddd)>1.d0) then
157  angle = 2.d0*atan(1.0) ! 0.5*pi
158  !FIX (lmh) to correctly handle co-linear points (angle near pi or 0)
159  if (ddd < 0.d0) then
160  angle = 4.d0*atan(1.0d0) !should be pi
161  else
162  angle = 0.d0
163  end if
164  else
165  angle = acos( ddd )
166  endif
167  endif
168 
169  spherical_angle = angle
170 
171  end function spherical_angle
172 
173 
183  real function get_area(p1, p4, p2, p3, radius)
184  !-----------------------------------------------
185  real, intent(in), dimension(2):: p1, p2, p3, p4
186  real, intent(in), optional:: radius
187  !-----------------------------------------------
188  real e1(3), e2(3), e3(3)
189  real ang1, ang2, ang3, ang4
190 
191  ! S-W: 1
192  call latlon2xyz(p1, e1) ! p1
193  call latlon2xyz(p2, e2) ! p2
194  call latlon2xyz(p4, e3) ! p4
195  ang1 = spherical_angle(e1, e2, e3)
196  !----
197  ! S-E: 2
198  !----
199  call latlon2xyz(p2, e1)
200  call latlon2xyz(p3, e2)
201  call latlon2xyz(p1, e3)
202  ang2 = spherical_angle(e1, e2, e3)
203  !----
204  ! N-E: 3
205  !----
206  call latlon2xyz(p3, e1)
207  call latlon2xyz(p4, e2)
208  call latlon2xyz(p2, e3)
209  ang3 = spherical_angle(e1, e2, e3)
210  !----
211  ! N-W: 4
212  !----
213  call latlon2xyz(p4, e1)
214  call latlon2xyz(p3, e2)
215  call latlon2xyz(p1, e3)
216  ang4 = spherical_angle(e1, e2, e3)
217 
218  if ( present(radius) ) then
219  get_area = (ang1 + ang2 + ang3 + ang4 - 2.*pi) * radius**2
220  else
221  get_area = ang1 + ang2 + ang3 + ang4 - 2.*pi
222  endif
223 
224  end function get_area
225 
236  subroutine fill_agrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
237  integer, intent(in) :: ng, npx, npy, isd, jsd
238  integer, intent(in) :: fill
239  real, DIMENSION(isd:,jsd:,:), intent(INOUT):: q
240 
241  integer :: i, j
242 
243  select case (fill)
244  case (xdir)
245  do j=1,ng
246  do i=1,ng
247  q(1-i ,1-j ,:) = q(1-j ,i ,:) !SW Corner
248  q(1-i ,npy-1+j,:) = q(1-j ,npy-1-i+1,:) !NW Corner
249  q(npx-1+i,1-j ,:) = q(npx-1+j,i ,:) !SE Corner
250  q(npx-1+i,npy-1+j,:) = q(npx-1+j,npy-1-i+1,:) !NE Corner
251  enddo
252  enddo
253  case (ydir)
254  do j=1,ng
255  do i=1,ng
256  q(1-j ,1-i ,:) = q(i ,1-j ,:) !SW Corner
257  q(1-j ,npy-1+i,:) = q(i ,npy-1+j,:) !NW Corner
258  q(npx-1+j,1-i ,:) = q(npx-1-i+1,1-j ,:) !SE Corner
259  q(npx-1+j,npy-1+i,:) = q(npx-1-i+1,npy-1+j,:) !NE Corner
260  enddo
261  enddo
262  case default
263  do j=1,ng
264  do i=1,ng
265  q(1-j ,1-i ,:) = q(i ,1-j ,:) !SW Corner
266  q(1-j ,npy-1+i,:) = q(i ,npy-1+j,:) !NW Corner
267  q(npx-1+j,1-i ,:) = q(npx-1-i+1,1-j ,:) !SE Corner
268  q(npx-1+j,npy-1+i,:) = q(npx-1-i+1,npy-1+j,:) !NE Corner
269  enddo
270  enddo
271  end select
272 
273 
274  end subroutine fill_agrid_scalar_corners
275 
286  subroutine fill_bgrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
287  integer, intent(in) :: ng, npx, npy, isd, jsd
288  integer, intent(in) :: fill
289  real, DIMENSION(isd:,jsd:,:), intent(INOUT):: q
290 
291  integer :: i, j
292 
293  select case (fill)
294  case (xdir)
295  do j=1,ng
296  do i=1,ng
297  q(1-i ,1-j ,:) = q(1-j ,i+1 ,:) !SW Corner
298  q(1-i ,npy+j,:) = q(1-j ,npy-i ,:) !NW Corner
299  q(npx+i,1-j ,:) = q(npx+j,i+1 ,:) !SE Corner
300  q(npx+i,npy+j,:) = q(npx+j,npy-i ,:) !NE Corner
301  enddo
302  enddo
303  case (ydir)
304  do j=1,ng
305  do i=1,ng
306  q(1-j ,1-i ,:) = q(i+1 ,1-j ,:) !SW Corner
307  q(1-j ,npy+i,:) = q(i+1 ,npy+j ,:) !NW Corner
308  q(npx+j,1-i ,:) = q(npx-i,1-j ,:) !SE Corner
309  q(npx+j,npy+i,:) = q(npx-i,npy+j ,:) !NE Corner
310  enddo
311  enddo
312  case default
313  do j=1,ng
314  do i=1,ng
315  q(1-i ,1-j ,:) = q(1-j ,i+1 ,:) !SW Corner
316  q(1-i ,npy+j,:) = q(1-j ,npy-i ,:) !NW Corner
317  q(npx+i,1-j ,:) = q(npx+j,i+1 ,:) !SE Corner
318  q(npx+i,npy+j,:) = q(npx+j,npy-i ,:) !NE Corner
319  enddo
320  enddo
321  end select
322 
323 
324 
325  end subroutine fill_bgrid_scalar_corners
326 
337  subroutine fill_agrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
338  integer, intent(in) :: ng, npx, npy, isd, jsd
339  real, DIMENSION(isd:,jsd:,:), intent(INOUT):: x
340  real, DIMENSION(isd:,jsd:,:), intent(INOUT):: y
341  integer :: i,j
342 
343  do j=1,ng
344  do i=1,ng
345  x(1-i ,1-j ,:) = y(1-j ,i ,:) !SW Corner
346  x(1-i ,npy-1+j,:) = y(1-j ,npy-1-i+1,:) !NW Corner
347  x(npx-1+i,1-j ,:) = y(npx-1+j,i ,:) !SE Corner
348  x(npx-1+i,npy-1+j,:) = y(npx-1+j,npy-1-i+1,:) !NE Corner
349 
350  y(1-j ,1-i ,:) = x(i ,1-j ,:) !SW Corner
351  y(1-j ,npy-1+i,:) = x(i ,npy-1+j,:) !NW Corner
352  y(npx-1+j,1-i ,:) = x(npx-1-i+1,1-j ,:) !SE Corner
353  y(npx-1+j,npy-1+i,:) = x(npx-1-i+1,npy-1+j,:) !NE Corner
354  enddo
355  enddo
356 
357  end subroutine fill_agrid_xy_corners
358 
369  subroutine fill_dgrid_xy_corners(x, y, ng, npx, npy, isd, jsd)
370  integer, intent(in) :: ng, npx, npy, isd, jsd
371  real, DIMENSION(isd:,jsd:,:), intent(INOUT):: x
372  real, DIMENSION(isd:,jsd:,:), intent(INOUT):: y
373  integer :: i,j
374 
375  do j=1,ng
376  do i=1,ng
377  x(1-i ,1-j , :) = y(1-j ,i , :) !SW Corner
378  x(1-i ,npy+j , :) = y(1-j ,npy-i, :) !NW Corner
379  x(npx-1+i,1-j , :) = y(npx+j,i , :) !SE Corner
380  x(npx-1+i,npy+j , :) = y(npx+j,npy-i, :) !NE Corner
381  y(1-i ,1-j , :) = x(j ,1-i , :) !SW Corner
382  y(1-i ,npy-1+j, :) = x(j ,npy+i, :) !NW Corner
383  y(npx+i ,1-j , :) = x(npx-j,1-i , :) !SE Corner
384  y(npx+i ,npy-1+j, :) = x(npx-j,npy+i, :) !NE Corner
385  enddo
386  enddo
387 
388  end subroutine fill_dgrid_xy_corners
389 
396  subroutine mid_pt_sphere(p1, p2, pm)
397  real, intent(IN) :: p1(2), p2(2)
398  real, intent(OUT) :: pm(2)
399  !------------------------------------------
400  real :: e1(3), e2(3), e3(3)
401 
402  call latlon2xyz(p1, e1)
403  call latlon2xyz(p2, e2)
404  call mid_pt3_cart(e1, e2, e3)
405  call cart_to_latlon(1, e3, pm(1), pm(2))
406 
407  end subroutine mid_pt_sphere
408 
415  subroutine mid_pt3_cart(p1, p2, e)
416  real, intent(IN) :: p1(3), p2(3)
417  real, intent(OUT) :: e(3)
418  !
419  real (f_p):: q1(3), q2(3)
420  real (f_p):: dd, e1, e2, e3
421  integer k
422 
423  do k=1,3
424  q1(k) = p1(k)
425  q2(k) = p2(k)
426  enddo
427 
428  e1 = q1(1) + q2(1)
429  e2 = q1(2) + q2(2)
430  e3 = q1(3) + q2(3)
431 
432  dd = sqrt( e1**2 + e2**2 + e3**2 )
433  e1 = e1 / dd
434  e2 = e2 / dd
435  e3 = e3 / dd
436 
437  e(1) = e1
438  e(2) = e2
439  e(3) = e3
440 
441  end subroutine mid_pt3_cart
442 
448  subroutine latlon2xyz(p, e)
449  !
450  ! Routine to map (lon, lat) to (x,y,z)
451  !
452  real, intent(in) :: p(2)
453  real, intent(out):: e(3)
454 
455  integer n
456  real (f_p):: q(2)
457  real (f_p):: e1, e2, e3
458 
459  do n=1,2
460  q(n) = p(n)
461  enddo
462 
463  e1 = cos(q(2)) * cos(q(1))
464  e2 = cos(q(2)) * sin(q(1))
465  e3 = sin(q(2))
466  !-----------------------------------
467  ! Truncate to the desired precision:
468  !-----------------------------------
469  e(1) = e1
470  e(2) = e2
471  e(3) = e3
472 
473  end subroutine latlon2xyz
474 
482  subroutine cart_to_latlon(np, q, xs, ys)
483  ! vector version of cart_to_latlon1
484  integer, intent(in):: np
485  real, intent(inout):: q(3,np)
486  real, intent(inout):: xs(np), ys(np)
487  ! local
488  real, parameter:: esl=1.d-10
489  real (f_p):: p(3)
490  real (f_p):: dist, lat, lon
491  integer i,k
492 
493  do i=1,np
494  do k=1,3
495  p(k) = q(k,i)
496  enddo
497  dist = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
498  do k=1,3
499  p(k) = p(k) / dist
500  enddo
501 
502  if ( (abs(p(1))+abs(p(2))) < esl ) then
503  lon = real(0.,kind=f_p)
504  else
505  lon = atan2( p(2), p(1) ) ! range [-pi,pi]
506  endif
507 
508  if ( lon < 0.) lon = real(2.,kind=f_p)*pi + lon
509  ! RIGHT_HAND system:
510  lat = asin(p(3))
511 
512  xs(i) = lon
513  ys(i) = lat
514  ! q Normalized:
515  do k=1,3
516  q(k,i) = p(k)
517  enddo
518  enddo
519 
520  end subroutine cart_to_latlon
521 
529  real function cos_angle(p1, p2, p3)
530  ! As spherical_angle, but returns the cos(angle)
531  ! p3
532  ! ^
533  ! |
534  ! |
535  ! p1 ---> p2
536  !
537  real, intent(in):: p1(3), p2(3), p3(3)
538 
539  real (f_p):: e1(3), e2(3), e3(3)
540  real (f_p):: px, py, pz
541  real (f_p):: qx, qy, qz
542  real (f_p):: angle, ddd
543  integer n
544 
545  do n=1,3
546  e1(n) = p1(n)
547  e2(n) = p2(n)
548  e3(n) = p3(n)
549  enddo
550 
551  !-------------------------------------------------------------------
552  ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
553  !-------------------------------------------------------------------
554  ! Vector P:= e1 X e2
555  px = e1(2)*e2(3) - e1(3)*e2(2)
556  py = e1(3)*e2(1) - e1(1)*e2(3)
557  pz = e1(1)*e2(2) - e1(2)*e2(1)
558 
559  ! Vector Q: e1 X e3
560  qx = e1(2)*e3(3) - e1(3)*e3(2)
561  qy = e1(3)*e3(1) - e1(1)*e3(3)
562  qz = e1(1)*e3(2) - e1(2)*e3(1)
563 
564  ! ddd = sqrt[ (P*P) (Q*Q) ]
565  ddd = sqrt( (px**2+py**2+pz**2)*(qx**2+qy**2+qz**2) )
566  if ( ddd > 0.d0 ) then
567  angle = (px*qx+py*qy+pz*qz) / ddd
568  else
569  angle = 1.d0
570  endif
571  cos_angle = angle
572 
573  end function cos_angle
574 
583  subroutine cell_center2(q1, q2, q3, q4, e2)
584  real, intent(in ) :: q1(2), q2(2), q3(2), q4(2)
585  real, intent(out) :: e2(2)
586  ! Local
587  real p1(3), p2(3), p3(3), p4(3)
588  real ec(3)
589  real dd
590  integer k
591 
592  call latlon2xyz(q1, p1)
593  call latlon2xyz(q2, p2)
594  call latlon2xyz(q3, p3)
595  call latlon2xyz(q4, p4)
596 
597  do k=1,3
598  ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
599  enddo
600  dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
601 
602  do k=1,3
603  ec(k) = ec(k) / dd
604  enddo
605 
606  call cart_to_latlon(1, ec, e2(1), e2(2))
607 
608  end subroutine cell_center2
609 
614  subroutine read_grid_file(regional)
615 
616  logical, intent(in) :: regional ! Is this a regional run?
617  integer :: fsize=65536
618  integer :: status, ncid, id_dim, id_var, ncid2, t
619  integer :: ni, nj, i, j, tw, te, ip
620  real :: g1(2), g2(2), g3(2), g4(2), g5(2)
621  real :: p1(3), p3(3)
622  real :: p_lL(2), p_uL(2), p_lR(2), p_uR(2)
623  character(len=512) :: tile_file
624  real, allocatable, dimension(:,:) :: tmpvar, geolon_c_nest, geolat_c_nest
625  real, allocatable, dimension(:,:,:) :: geolon_c, geolat_c
626  real, allocatable, dimension(:,:,:) :: geolon_t, geolat_t, cos_sg, grid3
627  integer :: start(4), nread(4)
628 
629  print*, "Read the grid from file "//trim(grid_file)
630 
631  status=nf__open(trim(grid_file),nf_nowrite,fsize,ncid)
632  call handle_err(status, 'Open file '//trim(grid_file) )
633 
634  status=nf_inq_dimid(ncid, 'ntiles', id_dim)
635  call handle_err(status, 'inquire dimension ntiles from file '//trim(grid_file) )
636  status=nf_inq_dimlen(ncid,id_dim,ntiles)
637  call handle_err(status, 'inquire dimension ntiles length from file '//trim(grid_file) )
638 
639  if( ntiles == 6) then
640  print*, " read_grid_file: This is a global grid."
641  elseif( ntiles >= 6 )then
642  print*, " read_grid_file: This is a nested grid."
643  print*, " filter only the globe."
644  ntiles=6
645  elseif( ntiles == 1 )then
646  print*, " read_grid_file: This is a standalone regional grid."
647  endif
648 
649  !--- loop through ntiles and make sure the grid size match between all the tiles.
650 
651  start(:) = 1
652  nread(:) = 1
653 
654  do t = 1, ntiles
655  start(2) = t; nread(1) = 255
656  status = nf_inq_varid(ncid, 'gridfiles', id_var)
657  call handle_err(status, 'inquire varid of gridfiles from file '//trim(grid_file) )
658 
659  !--- Obtain the grid file name from the mosaic file.
660  status = nf_get_vara_text(ncid, id_var, start, nread, tile_file )
661  call handle_err(status, 'get value of gridfiles from file '//trim(grid_file) )
662 
663  status=nf__open(trim(tile_file),nf_nowrite,fsize,ncid2)
664  call handle_err(status, 'Open file '//trim(tile_file) )
665 
666  status=nf_inq_dimid(ncid2, 'nx', id_dim)
667  call handle_err(status, 'inquire dimension nx from file '//trim(grid_file) )
668  status=nf_inq_dimlen(ncid2,id_dim,ni)
669  call handle_err(status, 'inquire dimension nx length from file '//trim(grid_file) )
670  status=nf_inq_dimid(ncid2, 'ny', id_dim)
671  call handle_err(status, 'inquire dimension ny from file '//trim(grid_file) )
672  status=nf_inq_dimlen(ncid2,id_dim,nj)
673  call handle_err(status, 'inquire dimension ny length '//'from file '//trim(grid_file) )
674  if( t == 1 ) then
675  ! ni and nj must be even
676  if(mod(ni,2) .NE. 0 .or. mod(nj,2) .NE. 0) &
677  call handle_err(-1, "read_grid_file: ni and nj must be even")
678 
679  nx = ni/2
680  ny = nj/2
681  npx = nx + 1
682  npy = ny + 1
683  is = 1 ; ie = nx
684  js = 1 ; je = ny
685  isd=is-ng; ied=ie+ng
686  jsd=js-ng; jed=je+ng
687 
688  allocate(tmpvar(ni+1,nj+1))
689  allocate(geolon_c(isd:ied+1,jsd:jed+1,6))
690  allocate(geolat_c(isd:ied+1,jsd:jed+1,6))
691  else if ( t == 7 ) then ! nested grid
692  if(mod(ni,2) .NE. 0 .or. mod(nj,2) .NE. 0) &
693  call handle_err(-1, "read_grid_file: ni and nj must be even")
694 
695  nx_nest = ni/2
696  ny_nest = nj/2
697  npx_nest = nx_nest + 1
698  npy_nest = ny_nest + 1
699  is_nest = 1 ; ie_nest = nx
700  js_nest = 1 ; je_nest = ny
701  isd_nest=is_nest-ng; ied_nest=ie_nest+ng
702  jsd_nest=js_nest-ng; jed_nest=je_nest+ng
703  deallocate(tmpvar)
704  allocate(tmpvar(ni+1,nj+1))
705  allocate(geolon_c_nest(isd:ied+1,jsd:jed+1))
706  allocate(geolat_c_nest(isd:ied+1,jsd:jed+1))
707  else
708  !-- make sure ni and nj match between tiles
709  if(ni .ne. nx*2 .OR. nj .ne. ny*2) &
710  call handle_err(-1, "mismatch of grid size between tiles")
711  endif
712 
713  status=nf_inq_varid(ncid2, 'x', id_var)
714  call handle_err(status, 'inquire varid of x from file '//trim(grid_file) )
715  status=nf_get_var_double(ncid2, id_var, tmpvar)
716  call handle_err(status, 'inquire data of x from file '//trim(grid_file) )
717  if(t==7) then
718  geolon_c_nest(1:npx,1:npy) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
719  else
720  geolon_c(1:npx,1:npy,t) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
721  endif
722 
723  status=nf_inq_varid(ncid2, 'y', id_var)
724  call handle_err(status, 'inquire varid of y from file '//trim(grid_file) )
725  status=nf_get_var_double(ncid2, id_var, tmpvar)
726  call handle_err(status, 'inquire data of y from file '//trim(grid_file) )
727  if(t==7) then
728  geolat_c_nest(1:npx,1:npy) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
729  else
730  geolat_c(1:npx,1:npy,t) = tmpvar(1:ni+1:2,1:nj+1:2)*pi/180.
731  endif
732  status = nf_close(ncid2)
733  call handle_err(status, "close file "//trim(tile_file))
734  enddo
735 
736  deallocate(tmpvar)
737 
738  status = nf_close(ncid)
739  call handle_err(status, "close file "//trim(grid_file))
740 
741  is = 1 ; ie = nx
742  js = 1 ; je = ny
743  isd=is-ng; ied=ie+ng
744  jsd=js-ng; jed=je+ng
745 
746  if( .not. regional ) then
747  call fill_cubic_grid_halo(geolon_c, geolon_c, ng, 1, 1, 1, 1)
748  call fill_cubic_grid_halo(geolat_c, geolat_c, ng, 1, 1, 1, 1)
749  if(.not. nested) call fill_bgrid_scalar_corners(geolon_c, ng, npx, npy, isd, jsd, xdir)
750  if(.not. nested) call fill_bgrid_scalar_corners(geolat_c, ng, npx, npy, isd, jsd, ydir)
751  else
752  call fill_regional_halo(geolon_c, ng)
753  call fill_regional_halo(geolat_c, ng)
754  endif
755 
756  !--- compute grid cell center
757  allocate(geolon_t(isd:ied,jsd:jed,ntiles), geolat_t(isd:ied,jsd:jed,ntiles))
758 
759  geolon_t(:,:,:) = -1.e25
760  geolat_t(:,:,:) = -1.e25
761 
762  do t = 1, ntiles
763  do j=js,je ; do i=is,ie
764  g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
765  g2(1) = geolon_c(i+1,j,t); g2(2) = geolat_c(i+1,j,t)
766  g3(1) = geolon_c(i,j+1,t); g3(2) = geolat_c(i,j+1,t)
767  g4(1) = geolon_c(i+1,j+1,t); g4(2) = geolat_c(i+1,j+1,t)
768  call cell_center2(g1, g2, g3, g4, g5 )
769  geolon_t(i,j,t) = g5(1)
770  geolat_t(i,j,t) = g5(2)
771  enddo ; enddo
772  enddo
773 
774 
775  if( .not. regional ) then
776  call fill_cubic_grid_halo(geolon_t, geolon_t, ng, 0, 0, 1, 1)
777  call fill_cubic_grid_halo(geolat_t, geolat_t, ng, 0, 0, 1, 1)
778  if (.not. nested) call fill_agrid_scalar_corners(geolon_t, ng, npx, npy, isd, jsd, xdir)
779  if (.not. nested) call fill_agrid_scalar_corners(geolat_t, ng, npx, npy, isd, jsd, ydir)
780  endif
781 
782  !--- compute dx, dy
783  allocate(dx(isd:ied,jsd:jed+1,ntiles))
784  allocate(dy(isd:ied+1,jsd:jed,ntiles))
785  do t = 1, ntiles
786  do j = js, je+1 ; do i = is, ie
787  g1(1) = geolon_c(i ,j,t)
788  g1(2) = geolat_c(i ,j,t)
789  g2(1) = geolon_c(i+1,j,t)
790  g2(2) = geolat_c(i+1,j,t)
791  dx(i,j,t) = great_circle_dist( g2, g1, radius )
792  enddo ; enddo
793  enddo
794  do t = 1, ntiles
795  do j = js, je
796  do i = is, ie+1
797  g1(1) = geolon_c(i,j, t)
798  g1(2) = geolat_c(i,j, t)
799  g2(1) = geolon_c(i,j+1,t)
800  g2(2) = geolat_c(i,j+1,t)
801  dy(i,j,t) = great_circle_dist( g2, g1, radius )
802  enddo
803  enddo
804  enddo
805 
806  if( .not. regional ) then
807  !--- make sure it is consitent between tiles. The following maybe not necessary.
808  do t = 1, ntiles
809  if(mod(t,2) ==0) then ! tile 2 4 6
810  tw = t - 1
811  te = t + 2
812  if(te > ntiles) te = te - ntiles
813  dy(is, js:je,t) = dy(ie+1,js:je,tw) ! west boundary
814  dy(ie+1, js:je, t) = dx(ie:is:-1,js, te) ! east boundary
815  else
816  tw = t - 2
817  if( tw <= 0) tw = tw + ntiles
818  te = t + 1
819  dy(is, js:je, t) = dx(ie:is:-1, je+1, tw) ! west boundary
820  dy(ie+1, js:je,t) = dy(1,js:je,te) ! east boundary
821  endif
822  enddo
823 
824  call fill_cubic_grid_halo(dx, dy, ng, 0, 1, 1, 1)
825  call fill_cubic_grid_halo(dy, dx, ng, 1, 0, 1, 1)
826 
827  if (.not. nested) call fill_dgrid_xy_corners(dx, dy, ng, npx, npy, isd, jsd)
828  endif
829 
830  !--- compute dxa and dya -----
831  allocate(dxa(isd:ied,jsd:jed,ntiles))
832  allocate(dya(isd:ied,jsd:jed,ntiles))
833  do t = 1, ntiles
834  do j=js,je ; do i=is,ie
835  g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
836  g2(1) = geolon_c(i,j+1,t); g2(2) = geolat_c(i,j+1,t)
837  call mid_pt_sphere(g1, g2, g3)
838  g1(1) = geolon_c(i+1,j,t); g1(2) = geolat_c(i+1,j,t)
839  g2(1) = geolon_c(i+1,j+1,t); g2(2) = geolat_c(i+1,j+1,t)
840  call mid_pt_sphere(g1, g2, g4)
841  dxa(i,j,t) = great_circle_dist( g4, g3, radius )
842  g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
843  g2(1) = geolon_c(i+1,j,t); g2(2) = geolat_c(i+1,j,t)
844  call mid_pt_sphere(g1, g2, g3)
845  g1(1) = geolon_c(i,j+1,t); g1(2) = geolat_c(i,j+1,t)
846  g2(1) = geolon_c(i+1,j+1,t); g2(2) = geolat_c(i+1,j+1,t)
847  call mid_pt_sphere(g1, g2, g4)
848  dya(i,j,t) = great_circle_dist( g4, g3, radius )
849  enddo; enddo
850  enddo
851 
852  if( .not.regional ) then
853  call fill_cubic_grid_halo(dxa, dya, ng, 0, 0, 1, 1)
854  call fill_cubic_grid_halo(dya, dxa, ng, 0, 0, 1, 1)
855 
856  if (.not. nested) call fill_agrid_xy_corners(dxa, dya, ng, npx, npy, isd, jsd)
857  endif
858 
859  !--- compute dxc and dyc
860  allocate(dxc(isd:ied+1,jsd:jed,ntiles))
861  allocate(dyc(isd:ied,jsd:jed+1,ntiles))
862  do t = 1, ntiles
863  do j=jsd,jed
864  do i=isd+1,ied
865  g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
866  g2(1) = geolon_c(i-1,j,t); g2(2) = geolat_c(i-1,j,t)
867  dxc(i,j,t) = great_circle_dist(g1, g2, radius)
868  enddo
869  dxc(isd,j,t) = dxc(isd+1,j,t)
870  dxc(ied+1,j,t) = dxc(ied,j,t)
871  enddo
872 
873  do j=jsd+1,jed
874  do i=isd,ied
875  g1(1) = geolon_c(i,j,t); g1(2) = geolat_c(i,j,t)
876  g2(1) = geolon_c(i,j-1,t); g2(2) = geolat_c(i,j-1,t)
877  dyc(i,j,t) = great_circle_dist(g1, g2, radius)
878  enddo
879  enddo
880  do i=isd,ied
881  dyc(i,jsd,t) = dyc(i,jsd+1,t)
882  dyc(i,jed+1,t) = dyc(i,jed,t)
883  end do
884  enddo
885 
886  !--- compute area
887  allocate(area(isd:ied,jsd:jed,ntiles))
888  do t = 1, ntiles
889  do j=js,je
890  do i=is,ie
891  p_ll(1) = geolon_c(i ,j ,t) ; p_ll(2) = geolat_c(i ,j ,t)
892  p_ul(1) = geolon_c(i ,j+1,t) ; p_ul(2) = geolat_c(i ,j+1,t)
893  p_lr(1) = geolon_c(i+1,j ,t) ; p_lr(2) = geolat_c(i+1,j ,t)
894  p_ur(1) = geolon_c(i+1,j+1,t) ; p_ur(2) = geolat_c(i+1,j+1,t)
895 
896  ! Spherical Excess Formula
897  area(i,j,t) = get_area(p_ll, p_ul, p_lr, p_ur, radius)
898  enddo
899  enddo
900  enddo
901 
902  if( .not.regional ) then
903  call fill_cubic_grid_halo(area, area, ng, 0, 0, 1, 1)
904  endif
905 
906  da_min = minval(area(is:ie,js:je,:))
907 
908  !--- compute sin_sg
909  allocate(sin_sg(4,isd:ied,jsd:jed,ntiles))
910  allocate(cos_sg(4,isd:ied,jsd:jed))
911  allocate(grid3(3, npx, npy))
912  cos_sg(:,:,:) = big_number
913  sin_sg(:,:,:,:) = tiny_number
914 
915  ! 9---4---8
916  ! | |
917  ! 1 5 3
918  ! | |
919  ! 6---2---7
920  do t = 1, ntiles
921  do j=js,je+1
922  do i = is,ie+1
923  g1(1) = geolon_c(i,j,t)
924  g1(2) = geolat_c(i,j,t)
925  call latlon2xyz(g1, grid3(:,i,j))
926  enddo
927  enddo
928  do j=js,je
929  do i=is,ie
930  g1(1) = geolon_t(i,j,t); g1(2) = geolat_t(i,j,t)
931  call latlon2xyz(g1, p3) ! righ-hand system consistent with grid3
932  call mid_pt3_cart(grid3(1,i,j), grid3(1,i,j+1), p1)
933  cos_sg(1,i,j) = cos_angle( p1, p3, grid3(1,i,j+1) )
934  call mid_pt3_cart(grid3(1,i,j), grid3(1,i+1,j), p1)
935  cos_sg(2,i,j) = cos_angle( p1, grid3(1,i+1,j), p3 )
936  call mid_pt3_cart(grid3(1,i+1,j), grid3(1,i+1,j+1), p1)
937  cos_sg(3,i,j) = cos_angle( p1, p3, grid3(1,i+1,j) )
938  call mid_pt3_cart(grid3(1,i,j+1), grid3(1,i+1,j+1), p1)
939  cos_sg(4,i,j) = cos_angle( p1, grid3(1,i,j+1), p3 )
940  enddo
941  enddo
942 
943  do ip=1,4
944  do j=js,je
945  do i=is,ie
946  sin_sg(ip,i,j,t) = min(1.0, sqrt( max(0., 1.-cos_sg(ip,i,j)**2) ) )
947  enddo
948  enddo
949  enddo
950  enddo
951 
952  if( .not.regional ) then
953  do ip=1,4
954  call fill_cubic_grid_halo(sin_sg(ip,:,:,:), sin_sg(ip,:,:,:), ng, 0, 0, 1, 1)
955  enddo
956  endif
957 
958  deallocate(cos_sg, grid3, geolon_c, geolat_c, geolon_t, geolat_t)
959 
960 
961  end subroutine read_grid_file
962 
967  subroutine read_topo_file(regional)
968 
969  logical,intent(in) :: regional ! Is this a run with a regional domain?
970  integer :: fsize=65536
971  integer :: status, ncid, id_var, ndim, dimsiz, nt
972  character(len=256) :: tile_file
973  character(len=32) :: text
974  integer :: len, t, dims(2)
975  real :: tmp(is:ie,js:je)
976 
977  allocate(oro(isd:ied,jsd:jed,ntiles))
978  allocate(mask(isd:ied,jsd:jed,ntiles))
979  oro = -big_number
980  mask = 0
981 
982  !--- make sure topo_file suffix is not ".nc"
983  len = len_trim(topo_file)
984  if( index(topo_file, '.nc', back=.true.) == len-2) then
985  call handle_err(-1, "remove .nc from namelist topo_file="//trim(topo_file) )
986  endif
987 
988  !--- loop through each tile file to get the orography
989  do nt = 1, ntiles
990 
991  if( regional ) then
992  t = nt + 6 ! The single regional tile must be #7 for now.
993  else
994  t = nt
995  endif
996 
997  write(text, '(i1.1)' ) t
998  tile_file = trim(topo_file)//'.tile'//trim(text)//'.nc'
999  status=nf__open(trim(tile_file),nf_nowrite,fsize,ncid)
1000  call handle_err(status, 'Open file '//trim(tile_file) )
1001 
1002  status=nf_inq_varid(ncid, topo_field, id_var)
1003  call handle_err(status, 'inquire varid of '//trim(topo_field)//' from file '//trim(tile_file) )
1004 
1005  status = nf_inq_varndims(ncid, id_var, ndim)
1006  call handle_err(status, 'inquire ndims of '//trim(topo_field)//' from file '//trim(tile_file) )
1007 
1008  if(ndim .NE. 2) call handle_err(-1, 'ndims of '//trim(topo_field)//' from file '// &
1009  trim(tile_file)//' should be 2')
1010 
1011  ! get data dimension and should match grid file size
1012  status = nf_inq_vardimid(ncid, id_var,dims);
1013  call handle_err(status, 'inquire dimid of '//trim(topo_field)//' from file '//trim(tile_file) )
1014 
1015  status = nf_inq_dimlen(ncid, dims(1), dimsiz)
1016  call handle_err(status, 'inquire first dimension length of '//trim(topo_field)//' from file '//trim(tile_file) )
1017  if(dimsiz .NE. nx) call handle_err(-1, "mismatch of lon dimension size between "// &
1018  trim(grid_file)//' and '//trim(tile_file) )
1019 
1020  status = nf_inq_dimlen(ncid, dims(2), dimsiz)
1021  call handle_err(status, 'inquire second dimension length of '//trim(topo_field)//' from file '//trim(tile_file) )
1022 
1023  if(dimsiz .NE. ny) call handle_err(-1, "mismatch of lat dimension size between "// &
1024  trim(grid_file)//' and '//trim(tile_file) )
1025 
1026  status = nf_get_var_double(ncid, id_var, oro(is:ie,js:je,nt))
1027  call handle_err(status, 'get the value of '//trim(topo_field)//' from file '//trim(tile_file) )
1028 
1029  status=nf_inq_varid(ncid, mask_field, id_var)
1030  call handle_err(status, 'inquire varid of '//trim(mask_field)//' from file '//trim(tile_file) )
1031 
1032  status = nf_get_var_double(ncid, id_var, tmp)
1033  call handle_err(status, 'get the value of '//trim(mask_field)//' from file '//trim(tile_file) )
1034 
1035  mask(is:ie,js:je,nt) = tmp
1036 
1037  status = nf_close(ncid)
1038  call handle_err(status, "close file "//trim(tile_file))
1039  enddo
1040 
1041  if( .not.regional ) then
1042  !--- update halo
1043  call fill_cubic_grid_halo(oro, oro, ng, 0, 0, 1, 1)
1044  call fill_cubic_grid_halo(mask, mask, ng, 0, 0, 1, 1)
1045  endif
1046 
1047  if( regional ) then
1048  call fill_regional_halo(oro, ng)
1049  oro(:,:,:) = max(oro(:,:,:),0.)
1050  call fill_regional_halo(mask, ng)
1051  mask(:,:,:) = min(max(mask(:,:,:),0.),1.)
1052  endif
1053 
1054 
1055  end subroutine read_topo_file
1056 
1067  subroutine write_topo_file(is,ie,js,je,ntiles,q,regional)
1068  integer, intent(in) :: is,ie,js,je,ntiles
1069  real, intent(in) :: q(is:ie,js:je,ntiles)
1070  logical, intent(in) :: regional
1071 
1072  integer :: fsize=65536
1073  integer :: nt, t, status, ncid, id_var
1074  character(len=256) :: tile_file
1075  character(len=3) :: text
1076  !--- loop through each tile file to update topo_field
1077 
1078  do nt = 1, ntiles
1079 
1080  if( regional ) then
1081  t = nt + 6
1082  else
1083  t = nt
1084  endif
1085 
1086  write(text, '(i1.1)' ) t
1087  tile_file = trim(topo_file)//'.tile'//trim(text)//'.nc'
1088  status=nf__open(trim(tile_file),nf_write,fsize,ncid)
1089  call handle_err(status, 'write_topo_file: Open file '//trim(tile_file) )
1090 
1091  status=nf_inq_varid(ncid, topo_field, id_var)
1092  call handle_err(status, 'write_topo_file:inquire varid of '//trim(topo_field)//' from file '//trim(tile_file) )
1093 
1094  status = nf_put_var_double(ncid, id_var, q(:,:,nt))
1095  call handle_err(status, 'write_topo_file: put the value of '//trim(topo_field)//' from file '//trim(tile_file) )
1096 
1097  status = nf_close(ncid)
1098  call handle_err(status, "write_topo_file: close file "//trim(tile_file))
1099  enddo
1100 
1101 
1102  end subroutine write_topo_file
1103 
1115  subroutine fill_cubic_grid_halo(data, data2, halo, ioff, joff, sign1, sign2)
1116  integer, intent(in) :: halo
1117  real, dimension(1-halo:,1-halo:,:), intent(inout) :: data, data2
1118  integer, intent(in) :: ioff, joff, sign1, sign2
1119  integer :: lw, le, ls, ln
1120  integer :: i, tile
1121 
1122  ntiles = size(data,3)
1123 
1124  do tile = 1, ntiles
1125  if(mod(tile,2) == 0) then ! tile 2, 4, 6
1126  lw = tile - 1; le = tile + 2; ls = tile - 2; ln = tile + 1
1127  if(le > 6 ) le = le - 6
1128  if(ls < 1 ) ls = ls + 6
1129  if(ln > 6 ) ln = ln - 6
1130  data(1-halo:0, 1:ny+joff, tile) = data(nx-halo+1:nx, 1:ny+joff, lw) ! west
1131  do i = 1, halo
1132  data(nx+i+ioff, 1:ny+joff, tile) = sign1*data2(nx+joff:1:-1, i+ioff, le) ! east
1133  end do
1134  do i = 1, halo
1135  data(1:nx+ioff, 1-i, tile) = sign2*data2(nx-i+1, ny+ioff:1:-1, ls) ! south
1136  end do
1137  data(1:nx+ioff, ny+1+joff:ny+halo+joff, tile) = data(1:nx+ioff, 1+joff:halo+joff, ln) ! north
1138  else ! tile 1, 3, 5
1139  lw = tile - 2; le = tile + 1; ls = tile - 1; ln = tile + 2
1140  if(lw < 1 ) lw = lw + 6
1141  if(ls < 1 ) ls = ls + 6
1142  if(ln > 6 ) ln = ln - 6
1143  do i = 1, halo
1144  data(1-i, 1:ny+joff, tile) = sign1*data2(nx+joff:1:-1, ny-i+1, lw) ! west
1145  end do
1146  data(nx+1+ioff:nx+halo+ioff, 1:ny+joff, tile) = data(1+ioff:halo+ioff, 1:ny+joff, le) ! east
1147  data(1:nx+ioff, 1-halo:0, tile) = data(1:nx+ioff, ny-halo+1:ny, ls) ! south
1148  do i = 1, halo
1149  data(1:nx+ioff, ny+i+joff, tile) = sign2*data2(i+joff, ny+ioff:1:-1, ln) ! north
1150  end do
1151  end if
1152  enddo
1153 
1154  end subroutine fill_cubic_grid_halo
1155 
1184  subroutine fv3_zs_filter (is, ie, js, je, isd, ied, jsd, jed, npx, npy, npx_global, ntiles, &
1185  grid_type, stretch_fac, nested, area, dxa, dya, dx, dy, dxc, dyc, &
1186  sin_sg, phis, regional )
1187  integer, intent(in) :: is, ie, js, je, ntiles
1188  integer, intent(in) :: isd, ied, jsd, jed, npx, npy, npx_global, grid_type
1189  real, intent(in), dimension(isd:ied,jsd:jed, ntiles)::area, dxa, dya
1190  real, intent(in), dimension(isd:ied, jsd:jed+1, ntiles):: dx, dyc
1191  real, intent(in), dimension(isd:ied+1,jsd:jed, ntiles):: dy, dxc
1192 
1193  real, intent(IN):: sin_sg(4,isd:ied,jsd:jed,ntiles)
1194  real, intent(IN):: stretch_fac
1195  logical, intent(IN) :: nested, regional
1196  real, intent(inout):: phis(isd:ied,jsd:jed,ntiles)
1197  real:: cd2
1198  integer mdim, n_del2, n_del4
1199 
1200  mdim = nint( real(npx_global) * min(10., stretch_fac) )
1201 
1202  ! Del-2: high resolution only
1203  if ( npx_global<=97 ) then
1204  n_del2 = 0
1205  elseif ( npx_global<=193 ) then
1206  n_del2 = 1
1207  else
1208  n_del2 = 2
1209  endif
1210  cd2 = 0.16*da_min
1211  ! Applying strong 2-delta-filter:
1212  if ( n_del2 > 0 ) &
1213  call two_delta_filter(is,ie,js,je,isd,ied,jsd,jed, npx, npy, ntiles, phis, area, &
1214  dx, dy, dxa, dya, dxc, dyc, sin_sg, cd2, zero_ocean, &
1215  .true.,0, grid_type, mask, nested, n_del2, regional)
1216 
1217  ! MFCT Del-4:
1218  if ( mdim<=193 ) then
1219  n_del4 = 1
1220  elseif ( mdim<=1537 ) then
1221  n_del4 = 2
1222  else
1223  n_del4 = 3
1224  endif
1225  call del4_cubed_sphere(is,ie,js,je,isd,ied,jsd,jed,npx, npy, ntiles, &
1226  phis, area, dx, dy, dxc, dyc, sin_sg, n_del4, zero_ocean, mask, nested, regional)
1227  ! Applying weak 2-delta-filter:
1228  cd2 = 0.12*da_min
1229  call two_delta_filter(is,ie,js,je,isd,ied,jsd,jed,npx, npy, ntiles, &
1230  phis, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd2, zero_ocean, &
1231  .true., 1, grid_type, mask, nested, n_del2_weak, regional)
1232 
1233  end subroutine fv3_zs_filter
1234 
1267  subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, &
1268  q, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd, zero_ocean, &
1269  check_slope, filter_type, grid_type, mask, nested, ntmax, regional)
1270  integer, intent(in) :: is, ie, js, je
1271  integer, intent(in) :: isd, ied, jsd, jed
1272  integer, intent(in) :: npx, npy, grid_type
1273  integer, intent(in) :: ntmax, ntiles
1274  integer, intent(in) :: filter_type ! 0: strong, 1: weak
1275  real, intent(in) :: cd
1276  ! INPUT arrays
1277  real, intent(in)::area(isd:ied, jsd:jed, ntiles)
1278  real, intent(in):: dx(isd:ied, jsd:jed+1, ntiles)
1279  real, intent(in):: dy(isd:ied+1,jsd:jed, ntiles)
1280  real, intent(in):: dxa(isd:ied, jsd:jed, ntiles)
1281  real, intent(in):: dya(isd:ied, jsd:jed, ntiles)
1282  real, intent(in):: dxc(isd:ied+1,jsd:jed, ntiles)
1283  real, intent(in):: dyc(isd:ied, jsd:jed+1, ntiles)
1284  real, intent(in):: sin_sg(4,isd:ied,jsd:jed, ntiles)
1285  real, intent(in):: mask(isd:ied, jsd:jed, ntiles) ! 0==water, 1==land
1286  logical, intent(in):: zero_ocean, check_slope
1287  logical, intent(in):: nested, regional
1288  ! OUTPUT arrays
1289  real, intent(inout):: q(isd:ied, jsd:jed,ntiles)
1290  ! Local:
1291  real, parameter:: p1 = 7./12.
1292  real, parameter:: p2 = -1./12.
1293  real, parameter:: c1 = -2./14.
1294  real, parameter:: c2 = 11./14.
1295  real, parameter:: c3 = 5./14.
1296 
1297  real:: ddx(is:ie+1,js:je), ddy(is:ie,js:je+1)
1298  logical:: extm(is-1:ie+1)
1299  logical:: ext2(is:ie,js-1:je+1)
1300  real:: a1(is-1:ie+2)
1301  real:: a2(is:ie,js-1:je+2)
1302  real:: a3(is:ie,js:je,ntiles)
1303  real:: smax, m_slope, fac
1304  integer:: i,j, nt, t
1305  integer:: is1, ie2, js1, je2
1306 
1307  if ( .not. nested .and. grid_type<3 ) then
1308  is1 = max(3,is-1); ie2 = min(npx-2,ie+2)
1309  js1 = max(3,js-1); je2 = min(npy-2,je+2)
1310  else
1311  is1 = is-1; ie2 = ie+2
1312  js1 = js-1; je2 = je+2
1313  end if
1314 
1315  if ( check_slope ) then
1316  m_slope = max_slope
1317  else
1318  m_slope = 10.
1319  endif
1320 
1321 
1322  do nt=1, ntmax
1323  if( .not.regional ) then
1324  call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1325  endif
1326 
1327  ! Check slope
1328  if ( nt==1 .and. check_slope ) then
1329  do t = 1, ntiles
1330  do j=js,je
1331  do i=is,ie+1
1332  ddx(i,j) = (q(i,j,t) - q(i-1,j,t))/dxc(i,j,t)
1333  ddx(i,j) = abs(ddx(i,j))
1334  enddo
1335  enddo
1336  do j=js,je+1
1337  do i=is,ie
1338  ddy(i,j) = (q(i,j,t) - q(i,j-1,t))/dyc(i,j,t)
1339  ddy(i,j) = abs(ddy(i,j))
1340  enddo
1341  enddo
1342  do j=js,je
1343  do i=is,ie
1344  a3(i,j,t) = max( ddx(i,j), ddx(i+1,j), ddy(i,j), ddy(i,j+1) )
1345  enddo
1346  enddo
1347  enddo
1348  smax = maxval(a3(is:ie,js:je,:))
1349  write(*,*) 'Before filter: Max_slope=', smax
1350  endif
1351 
1352 
1353  ! First step: average the corners:
1354  if ( .not. nested .and. nt==1 ) then
1355  do t = 1, ntiles
1356  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)) &
1357  / ( area(1,1,t)+ area(0,1,t)+ area(1,0,t) )
1358  q(0,1,t) = q(1,1,t)
1359  q(1,0,t) = q(1,1,t)
1360 
1361  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)) &
1362  / ( area(ie,1,t)+ area(npx,1,t)+ area(ie,0,t))
1363  q(npx,1,t) = q(ie,1,t)
1364  q(ie, 0,t) = q(ie,1,t)
1365 
1366  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)) &
1367  / ( area(1,je,t)+ area(0,je,t)+ area(1,npy,t))
1368  q(0, je,t) = q(1,je,t)
1369  q(1,npy,t) = q(1,je,t)
1370 
1371  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)) &
1372  / ( area(ie,je,t)+ area(npx,je,t)+ area(ie,npy,t))
1373  q(npx,je,t) = q(ie,je,t)
1374  q(ie,npy,t) = q(ie,je,t)
1375  enddo
1376  if( .not.regional ) then
1377  call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1378  endif
1379  endif
1380 
1381  do t = 1, ntiles
1382  ! x-diffusive flux:
1383  do j=js,je
1384 
1385  do i=is1, ie2
1386  a1(i) = p1*(q(i-1,j,t)+q(i,j,t)) + p2*(q(i-2,j,t)+q(i+1,j,t))
1387  enddo
1388 
1389  if ( (.not. (nested .or. regional)) .and. grid_type<3 ) then
1390  a1(0) = c1*q(-2,j,t) + c2*q(-1,j,t) + c3*q(0,j,t)
1391  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)) &
1392  + ((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)))
1393  a1(2) = c3*q(1,j,t) + c2*q(2,j,t) +c1*q(3,j,t)
1394 
1395  a1(npx-1) = c1*q(npx-3,j,t) + c2*q(npx-2,j,t) + c3*q(npx-1,j,t)
1396  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)) &
1397  /(dxa(npx-2,j,t)+dxa(npx-1,j,t)) &
1398  + ((2.*dxa(npx, j,t)+dxa(npx+1,j,t))*q(npx, j,t)-dxa(npx, j,t)*q(npx+1,j,t))/ &
1399  (dxa(npx, j,t)+dxa(npx+1,j,t)))
1400  a1(npx+1) = c3*q(npx,j,t) + c2*q(npx+1,j,t) + c1*q(npx+2,j,t)
1401  endif
1402 
1403  if ( regional .and. grid_type<3 ) then
1404  a1(0) = c1*q(-1,j,t) + c2*q(-1,j,t) + c3*q(0,j,t)
1405  a1(2) = c3*q(1,j,t) + c2*q(2,j,t) + c1*q(3,j,t)
1406  a1(1) = 0.5*(a1(0) + a1(2))
1407 
1408  a1(npx-1) = c1*q(npx-3,j,t) + c2*q(npx-2,j,t) + c3*q(npx-1,j,t)
1409  a1(npx+1) = c3*q(npx,j,t) + c2*q(npx+1,j,t) + c1*q(npx+2,j,t)
1410  a1(npx) = 0.5*(a1(npx-1)+a1(npx+1))
1411  endif
1412 
1413  if ( filter_type == 0 ) then
1414  do i=is-1, ie+1
1415  if( abs(3.*(a1(i)+a1(i+1)-2.*q(i,j,t))) > abs(a1(i)-a1(i+1)) ) then
1416  extm(i) = .true.
1417  else
1418  extm(i) = .false.
1419  endif
1420  enddo
1421  else
1422  do i=is-1, ie+1
1423  if ( (a1(i)-q(i,j,t))*(a1(i+1)-q(i,j,t)) > 0. ) then
1424  extm(i) = .true.
1425  else
1426  extm(i) = .false.
1427  endif
1428  enddo
1429  endif
1430 
1431  do i=is,ie+1
1432  ddx(i,j) = (q(i-1,j,t)-q(i,j,t))/dxc(i,j,t)
1433  if ( extm(i-1).and.extm(i) ) then
1434  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)
1435  elseif ( abs(ddx(i,j)) > m_slope ) then
1436  fac = min(1., max( 0.1, ( abs(ddx(i,j))-m_slope )/m_slope) )
1437  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)
1438  else
1439  ddx(i,j) = 0.
1440  endif
1441  enddo
1442  enddo ! do j=js,je
1443 
1444  ! y-diffusive flux:
1445  do j=js1,je2
1446  do i=is,ie
1447  a2(i,j) = p1*(q(i,j-1,t)+q(i,j,t)) + p2*(q(i,j-2,t)+q(i,j+1,t))
1448  enddo
1449  enddo
1450  if ( (.not. (nested .or. regional)) .and. grid_type<3 ) then
1451  do i=is,ie
1452  a2(i,0) = c1*q(i,-2,t) + c2*q(i,-1,t) + c3*q(i,0,t)
1453  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)) &
1454  + ((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)))
1455  a2(i,2) = c3*q(i,1,t) + c2*q(i,2,t) + c1*q(i,3,t)
1456  enddo
1457 
1458  do i=is,ie
1459  a2(i,npy-1) = c1*q(i,npy-3,t) + c2*q(i,npy-2,t) + c3*q(i,npy-1,t)
1460  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))/ &
1461  (dya(i,npy-2,t)+dya(i,npy-1,t)) &
1462  + ((2.*dya(i,npy,t)+dya(i,npy+1,t))*q(i,npy,t)-dya(i,npy,t)*q(i,npy+1,t))/&
1463  (dya(i,npy,t)+dya(i,npy+1,t)))
1464  a2(i,npy+1) = c3*q(i,npy,t) + c2*q(i,npy+1,t) + c1*q(i,npy+2,t)
1465  enddo
1466  endif
1467 
1468  if ( regional .and. grid_type<3 ) then
1469  do i=is,ie
1470  a2(i,0) = c1*q(i,-2,t) + c2*q(i,-1,t) + c3*q(i,0,t)
1471  a2(i,2) = c3*q(i,1,t) + c2*q(i,2,t) + c1*q(i,3,t)
1472  a2(i,1) = 0.5*(a2(i,0) + a2(i,2))
1473  enddo
1474 
1475  do i=is,ie
1476  a2(i,npy-1) = c1*q(i,npy-3,t) + c2*q(i,npy-2,t) + c3*q(i,npy-1,t)
1477  a2(i,npy+1) = c3*q(i,npy,t) + c2*q(i,npy+1,t) + c1*q(i,npy+2,t)
1478  a2(i,npy) = 0.5*(a2(i,npy-1)+a2(i,npy+1))
1479  enddo
1480  endif
1481 
1482  if ( filter_type == 0 ) then
1483  do j=js-1,je+1
1484  do i=is,ie
1485  if( abs(3.*(a2(i,j)+a2(i,j+1)-2.*q(i,j,t))) > abs(a2(i,j)-a2(i,j+1)) ) then
1486  ext2(i,j) = .true.
1487  else
1488  ext2(i,j) = .false.
1489  endif
1490  enddo
1491  enddo
1492  else
1493  do j=js-1,je+1
1494  do i=is,ie
1495  if ( (a2(i,j)-q(i,j,t))*(a2(i,j+1)-q(i,j,t)) > 0. ) then
1496  ext2(i,j) = .true.
1497  else
1498  ext2(i,j) = .false.
1499  endif
1500  enddo
1501  enddo
1502  endif
1503 
1504  do j=js,je+1
1505  do i=is,ie
1506  ddy(i,j) = (q(i,j-1,t)-q(i,j,t))/dyc(i,j,t)
1507  if ( ext2(i,j-1) .and. ext2(i,j) ) then
1508  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)
1509  elseif ( abs(ddy(i,j))>m_slope ) then
1510  fac = min(1., max(0.1,(abs(ddy(i,j))-m_slope)/m_slope))
1511  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)
1512  else
1513  ddy(i,j) = 0.
1514  endif
1515  enddo
1516  enddo
1517 
1518  if ( zero_ocean ) then
1519  ! Limit diffusive flux over water cells:
1520  do j=js,je
1521  do i=is,ie+1
1522  ddx(i,j) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * ddx(i,j)
1523  enddo
1524  enddo
1525  do j=js,je+1
1526  do i=is,ie
1527  ddy(i,j) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * ddy(i,j)
1528  enddo
1529  enddo
1530  endif
1531 
1532  do j=js,je
1533  do i=is,ie
1534  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))
1535  enddo
1536  enddo
1537  enddo ! do t = 1, ntiles
1538  enddo ! nt=1, ntmax
1539 
1540 ! Check slope
1541  if ( check_slope ) then
1542  if( .not.regional ) then
1543  call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1544  endif
1545  do t = 1, ntiles
1546  do j=js,je
1547  do i=is,ie+1
1548  ddx(i,j) = (q(i,j,t) - q(i-1,j,t))/dxc(i,j,t)
1549  ddx(i,j) = abs(ddx(i,j))
1550  enddo
1551  enddo
1552  do j=js,je+1
1553  do i=is,ie
1554  ddy(i,j) = (q(i,j,t) - q(i,j-1,t))/dyc(i,j,t)
1555  ddy(i,j) = abs(ddy(i,j))
1556  enddo
1557  enddo
1558  do j=js,je
1559  do i=is,ie
1560  a3(i,j,t) = max( ddx(i,j), ddx(i+1,j), ddy(i,j), ddy(i,j+1) )
1561  enddo
1562  enddo
1563  enddo
1564  smax = maxval(a3(is:ie,js:je,:))
1565  write(*,*) 'After filter: Max_slope=', smax
1566  endif
1567 
1568  end subroutine two_delta_filter
1569 
1597  subroutine del2_cubed_sphere(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles,&
1598  q, area, dx, dy, dxc, dyc, sin_sg, nmax, cd, zero_ocean, mask, nested, regional)
1599  integer, intent(in) :: is, ie, js, je
1600  integer, intent(in) :: isd, ied, jsd, jed
1601  integer, intent(in):: npx, npy, ntiles
1602  integer, intent(in):: nmax
1603  real, intent(in):: cd
1604  logical, intent(in):: zero_ocean
1605  ! INPUT arrays
1606  real, intent(in)::area(isd:ied, jsd:jed, ntiles)
1607  real, intent(in):: dx(isd:ied, jsd:jed+1, ntiles)
1608  real, intent(in):: dy(isd:ied+1,jsd:jed, ntiles)
1609  real, intent(in):: dxc(isd:ied+1,jsd:jed, ntiles)
1610  real, intent(in):: dyc(isd:ied, jsd:jed+1, ntiles)
1611  real, intent(IN):: sin_sg(4,isd:ied,jsd:jed, ntiles)
1612  real, intent(in):: mask(isd:ied, jsd:jed, ntiles) ! 0==water, 1==land
1613  logical, intent(IN) :: nested, regional
1614 
1615  ! OUTPUT arrays
1616  real, intent(inout):: q(is-ng:ie+ng, js-ng:je+ng, ntiles)
1617  ! Local:
1618  real ddx(is:ie+1,js:je), ddy(is:ie,js:je+1)
1619  integer i,j,n,t
1620 
1621  if( .not.regional ) then
1622  call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1623  endif
1624 
1625  do t = 1, ntiles
1626  ! First step: average the corners:
1627  if ( .not. nested) then
1628  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)) &
1629  / ( area(1,1,t)+ area(0,1,t)+ area(1,0,t) )
1630  q(0,1,t) = q(1,1,t)
1631  q(1,0,t) = q(1,1,t)
1632  endif
1633  if ( .not. nested) then
1634  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)) &
1635  / ( area(ie,1,t)+ area(npx,1,t)+ area(ie,0,t))
1636  q(npx,1,t) = q(ie,1,t)
1637  q(ie, 0,t) = q(ie,1,t)
1638  endif
1639  if ( .not. nested ) then
1640  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)) &
1641  / ( area(ie,je,t)+ area(npx,je,t)+ area(ie,npy,t))
1642  q(npx,je,t) = q(ie,je,t)
1643  q(ie,npy,t) = q(ie,je,t)
1644  endif
1645  if ( .not. nested) then
1646  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)) &
1647  / ( area(1,je,t)+ area(0,je,t)+ area(1,npy,t))
1648  q(0, je,t) = q(1,je,t)
1649  q(1,npy,t) = q(1,je,t)
1650  endif
1651  enddo
1652 
1653  do n=1,nmax
1654  if( .not.regional ) then
1655  if( n>1 ) call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1656  endif
1657  do t = 1, ntiles
1658  do j=js,je
1659  do i=is,ie+1
1660  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)
1661  enddo
1662  enddo
1663  do j=js,je+1
1664  do i=is,ie
1665  ddy(i,j) = dx(i,j,t)*(q(i,j-1,t)-q(i,j,t))/dyc(i,j,t) &
1666  *0.5*(sin_sg(4,i,j-1,t)+sin_sg(2,i,j,t))
1667  enddo
1668  enddo
1669 
1670  if ( zero_ocean ) then
1671  ! Limit diffusive flux over ater cells:
1672  do j=js,je
1673  do i=is,ie+1
1674  ddx(i,j) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * ddx(i,j)
1675  enddo
1676  enddo
1677  do j=js,je+1
1678  do i=is,ie
1679  ddy(i,j) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * ddy(i,j)
1680  enddo
1681  enddo
1682  endif
1683 
1684  do j=js,je
1685  do i=is,ie
1686  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))
1687  enddo
1688  enddo
1689  enddo
1690  enddo
1691 
1692  end subroutine del2_cubed_sphere
1693 
1720  subroutine del4_cubed_sphere(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles, &
1721  q, area, dx, dy, dxc, dyc, sin_sg, nmax, zero_ocean, mask, nested, regional)
1722  integer, intent(in) :: is, ie, js, je
1723  integer, intent(in) :: isd, ied, jsd, jed
1724  integer, intent(in) :: npx, npy, nmax, ntiles
1725  logical, intent(in) :: zero_ocean
1726  real, intent(in):: mask(isd:ied, jsd:jed, ntiles) ! 0==water, 1==land
1727  real, intent(in)::area(isd:ied, jsd:jed, ntiles)
1728  real, intent(in):: dx(isd:ied, jsd:jed+1, ntiles)
1729  real, intent(in):: dy(isd:ied+1,jsd:jed, ntiles)
1730  real, intent(in):: dxc(isd:ied+1,jsd:jed, ntiles)
1731  real, intent(in):: dyc(isd:ied, jsd:jed+1, ntiles)
1732  real, intent(IN):: sin_sg(4,isd:ied,jsd:jed, ntiles)
1733  real, intent(inout):: q(isd:ied, jsd:jed, ntiles)
1734  logical, intent(IN) :: nested, regional
1735  ! Local:
1736  ! diffusivity
1737  real :: diff(is-1:ie+1,js-1:je+1, ntiles)
1738  ! diffusive fluxes:
1739  real :: fx1(is:ie+1,js:je), fy1(is:ie,js:je+1)
1740  real :: fx2(is:ie+1,js:je,ntiles), fy2(is:ie,js:je+1,ntiles)
1741  real :: fx4(is:ie+1,js:je,ntiles), fy4(is:ie,js:je+1,ntiles)
1742  real, dimension(isd:ied,jsd:jed,ntiles):: d2, win, wou
1743  real, dimension(is:ie,js:je, ntiles) :: qlow, qmin, qmax
1744  real, parameter:: esl = 1.e-20
1745  integer i,j, n, t
1746 
1747  ! On a nested grid the haloes are not filled. Set to zero.
1748  d2 = 0.
1749  win = 0.
1750  wou = 0.
1751 
1752  do t = 1, ntiles
1753  do j=js-1,je+1 ; do i=is-1,ie+1
1754  diff(i,j,t) = cd4*area(i,j,t) ! area dependency is needed for stretched grid
1755  enddo; enddo
1756 
1757  do j=js,je ; do i=is,ie
1758  qmax(i,j,t) = q(i,j,t) * peak_fac
1759  qmin(i,j,t) = q(i,j,t) / peak_fac
1760  enddo; enddo
1761  enddo
1762 
1763  do n=1,nmax
1764  if( .not.regional ) then
1765  call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1766  endif
1767 
1768  ! First step: average the corners:
1769  if ( .not. nested .and. n==1 ) then
1770  do t = 1, ntiles
1771  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)) &
1772  / ( area(1,1,t)+ area(0,1,t)+ area(1,0,t) )
1773  q(0,1,t) = q(1,1,t)
1774  q(1,0,t) = q(1,1,t)
1775  q(0,0,t) = q(1,1,t)
1776 
1777  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)) &
1778  / ( area(ie,1,t)+ area(npx,1,t)+ area(ie,0,t))
1779  q(npx,1,t) = q(ie,1,t)
1780  q(ie, 0,t) = q(ie,1,t)
1781  q(npx,0,t) = q(ie,1,t)
1782 
1783  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)) &
1784  / ( area(1,je,t)+ area(0,je,t)+ area(1,npy,t))
1785  q(0, je,t) = q(1,je,t)
1786  q(1,npy,t) = q(1,je,t)
1787  q(0,npy,t) = q(1,je,t)
1788 
1789  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)) &
1790  / ( area(ie,je,t)+ area(npx,je,t)+ area(ie,npy,t))
1791  q(npx, je,t) = q(ie,je,t)
1792  q(ie, npy,t) = q(ie,je,t)
1793  q(npx,npy,t) = q(ie,je,t)
1794  enddo
1795  if( .not.regional ) then
1796  call fill_cubic_grid_halo(q, q, ng, 0, 0, 1, 1)
1797  endif
1798  endif
1799 
1800  do t = 1, ntiles
1801 
1802  !--------------
1803  ! Compute del-2
1804  !--------------
1805  ! call copy_corners(q, npx, npy, 1)
1806  do j=js,je
1807  do i=is,ie+1
1808  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) &
1809  *(sin_sg(1,i,j,t)+sin_sg(3,i-1,j,t))
1810  enddo
1811  enddo
1812 
1813  ! call copy_corners(q, npx, npy, 2)
1814  do j=js,je+1
1815  do i=is,ie
1816  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) &
1817  *(sin_sg(2,i,j,t)+sin_sg(4,i,j-1,t))
1818  enddo
1819  enddo
1820 
1821  do j=js,je
1822  do i=is,ie
1823  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)
1824  enddo
1825  enddo
1826 
1827  ! qlow == low order monotonic solution
1828  if ( zero_ocean ) then
1829  ! Limit diffusive flux over water cells:
1830  do j=js,je
1831  do i=is,ie+1
1832  fx1(i,j) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * fx2(i,j,t)
1833  enddo
1834  enddo
1835  do j=js,je+1
1836  do i=is,ie
1837  fy1(i,j) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * fy2(i,j,t)
1838  enddo
1839  enddo
1840  do j=js,je
1841  do i=is,ie
1842  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)
1843  d2(i,j,t) = diff(i,j,t) * d2(i,j,t)
1844  enddo
1845  enddo
1846  else
1847  do j=js,je
1848  do i=is,ie
1849  qlow(i,j,t) = q(i,j,t) + d2(i,j,t)
1850  d2(i,j,t) = diff(i,j,t) * d2(i,j,t)
1851  enddo
1852  enddo
1853  endif
1854  enddo
1855  if( .not.regional ) then
1856  call fill_cubic_grid_halo(d2, d2, ng, 0, 0, 1, 1)
1857  endif
1858 
1859  !---------------------
1860  ! Compute del4 fluxes:
1861  !---------------------
1862  ! call copy_corners(d2, npx, npy, 1)
1863  do t = 1, ntiles
1864  do j=js,je
1865  do i=is,ie+1
1866  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)
1867  enddo
1868  enddo
1869 
1870  ! call copy_corners(d2, npx, npy, 2)
1871  do j=js,je+1
1872  do i=is,ie
1873  fy4(i,j,t) = dx(i,j,t)*(d2(i,j,t)-d2(i,j-1,t))/dyc(i,j,t) &
1874  *0.5*(sin_sg(2,i,j,t)+sin_sg(4,i,j-1,t))-fy2(i,j,t)
1875  enddo
1876  enddo
1877 
1878  do j=js,je
1879  do i=is,ie
1880  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), &
1881  q(i-1,j ,t), q(i,j ,t), q(i+1,j ,t), &
1882  q(i-1,j+1,t), q(i,j+1,t), q(i+1,j+1,t) )
1883  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), &
1884  q(i-1,j ,t), q(i,j ,t), q(i+1,j ,t), &
1885  q(i-1,j+1,t), q(i,j+1,t), q(i+1,j+1,t) )
1886  enddo
1887  enddo
1888 
1889  !----------------
1890  ! Flux limitting:
1891  !----------------
1892  do j=js,je
1893  do i=is,ie
1894  win(i,j,t) = max(0.,fx4(i, j,t)) - min(0.,fx4(i+1,j,t)) + &
1895  max(0.,fy4(i, j,t)) - min(0.,fy4(i,j+1,t)) + esl
1896  wou(i,j,t) = max(0.,fx4(i+1,j,t)) - min(0.,fx4(i, j,t)) + &
1897  max(0.,fy4(i,j+1,t)) - min(0.,fy4(i, j,t)) + esl
1898  win(i,j,t) = max(0., qmax(i,j,t) - qlow(i,j,t)) / win(i,j,t)*area(i,j,t)
1899  wou(i,j,t) = max(0., qlow(i,j,t) - qmin(i,j,t)) / wou(i,j,t)*area(i,j,t)
1900  enddo
1901  enddo
1902  enddo
1903  if( .not.regional ) then
1904  call fill_cubic_grid_halo(win, win, ng, 0, 0, 1, 1)
1905  call fill_cubic_grid_halo(wou, wou, ng, 0, 0, 1, 1)
1906  endif
1907  do t = 1, ntiles
1908  do j=js,je
1909  do i=is,ie+1
1910  if ( fx4(i,j,t) > 0. ) then
1911  fx4(i,j,t) = min(1., wou(i-1,j,t), win(i,j,t)) * fx4(i,j,t)
1912  else
1913  fx4(i,j,t) = min(1., win(i-1,j,t), wou(i,j,t)) * fx4(i,j,t)
1914  endif
1915  enddo
1916  enddo
1917  do j=js,je+1
1918  do i=is,ie
1919  if ( fy4(i,j,t) > 0. ) then
1920  fy4(i,j,t) = min(1., wou(i,j-1,t), win(i,j,t)) * fy4(i,j,t)
1921  else
1922  fy4(i,j,t) = min(1., win(i,j-1,t), wou(i,j,t)) * fy4(i,j,t)
1923  endif
1924  enddo
1925  enddo
1926 
1927 
1928  if ( zero_ocean ) then
1929  ! Limit diffusive flux over ocean cells:
1930  do j=js,je
1931  do i=is,ie+1
1932  fx4(i,j,t) = max(0., min(mask(i-1,j,t), mask(i,j,t))) * fx4(i,j,t)
1933  enddo
1934  enddo
1935  do j=js,je+1
1936  do i=is,ie
1937  fy4(i,j,t) = max(0., min(mask(i,j-1,t), mask(i,j,t))) * fy4(i,j,t)
1938  enddo
1939  enddo
1940  endif
1941 
1942  ! Update:
1943  do j=js,je
1944  do i=is,ie
1945  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)
1946  enddo
1947  enddo
1948  enddo
1949  enddo ! end n-loop
1950 
1951 
1952  end subroutine del4_cubed_sphere
1953 
1958  subroutine check(status)
1959  use netcdf
1960  integer,intent(in) :: status
1961 !
1962  if(status /= nf90_noerr) then
1963  write(0,*) ' check netcdf status = ', status
1964  write(0,'("error ", a)') trim(nf90_strerror(status))
1965  write(0,*) "Stopped"
1966  stop 4
1967  endif
1968  end subroutine check
1969 
1974  subroutine compute_filter_constants
1975 
1976  ! set the given values for various cube resolutions (c48, c96, c192, c384, c768, c1152, c3072)
1977 
1978  integer,parameter :: nres=8
1979  integer :: index1,index2,n
1980 
1981  real :: factor
1982 
1983  real,dimension(1:nres) :: cube_res=(/48.,96.,128.,192.,384.,768.,1152.,3072./)
1984 
1985  real,dimension(1:nres) :: n_del2_weak_vals=(/4.,8.,8.,12.,12.,16.,20.,24./)
1986  real,dimension(1:nres) :: cd4_vals =(/0.12,0.12,0.13,0.15,0.15,0.15,0.15,0.15/)
1987  real,dimension(1:nres) :: max_slope_vals =(/0.12,0.12,0.12,0.12,0.12,0.12,0.16,0.30/)
1988  real,dimension(1:nres) :: peak_fac_vals =(/1.1,1.1,1.1,1.05,1.0,1.0,1.0,1.0/)
1989 
1990  if(res<cube_res(1))then
1991  index1 = 1
1992  index2 = 1
1993  factor = 0.
1994  elseif(res>cube_res(nres))then
1995  index1 = nres
1996  index2 = nres
1997  factor = 0.
1998  else
1999  do n=2,nres
2000  if(res<=cube_res(n))then
2001  index2 = n
2002  index1 = n-1
2003  factor = (res-cube_res(n-1))/(cube_res(n)-cube_res(n-1))
2004  exit
2005  endif
2006  enddo
2007  endif
2008 
2009  n_del2_weak = nint(n_del2_weak_vals(index1)+factor*(n_del2_weak_vals(index2)-n_del2_weak_vals(index1)))
2010  cd4 = cd4_vals(index1)+factor*(cd4_vals(index2)-cd4_vals(index1))
2011  max_slope = max_slope_vals(index1)+factor*(max_slope_vals(index2)-max_slope_vals(index1))
2012  peak_fac = peak_fac_vals(index1)+factor*(peak_fac_vals(index2)-peak_fac_vals(index1))
2013 
2014  print*,''
2015  print*,'- FILTER COEFFICIENTS FOR RESOLUTION ', res
2016  print*,'- CD4: ', cd4
2017  print*,'- N_DEL2_WEAK: ', n_del2_weak
2018  print*,'- MAX_SLOPE: ', max_slope
2019  print*,'- PEAK_FAC: ', peak_fac
2020 
2021  end subroutine compute_filter_constants
2022 
2023 end program filter_topo
cart_to_latlon
subroutine cart_to_latlon(np, q, xs, ys)
???
Definition: filter_topo.F90:483
cell_center2
subroutine cell_center2(q1, q2, q3, q4, e2)
???
Definition: filter_topo.F90:584
read_topo_file
subroutine read_topo_file(regional)
???
Definition: filter_topo.F90:968
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:416
cos_angle
real function cos_angle(p1, p2, p3)
???
Definition: filter_topo.F90:530
fill_agrid_scalar_corners
subroutine fill_agrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
???
Definition: filter_topo.F90:237
latlon2xyz
subroutine latlon2xyz(p, e)
???
Definition: filter_topo.F90:449
great_circle_dist
real function great_circle_dist(q1, q2, radius)
???
Definition: filter_topo.F90:82
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:1722
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:1599
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:123
get_area
real function get_area(p1, p4, p2, p3, radius)
???
Definition: filter_topo.F90:184
fill_bgrid_scalar_corners
subroutine fill_bgrid_scalar_corners(q, ng, npx, npy, isd, jsd, fill)
???
Definition: filter_topo.F90:287
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:1187
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:338
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:1116
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:370
read_grid_file
subroutine read_grid_file(regional)
???
Definition: filter_topo.F90:615
write_topo_file
subroutine write_topo_file(is, ie, js, je, ntiles, q, regional)
Replace the topo_field.
Definition: filter_topo.F90:1068
mid_pt_sphere
subroutine mid_pt_sphere(p1, p2, pm)
???
Definition: filter_topo.F90:397
check
subroutine check(status)
Check results of netCDF call.
Definition: filter_topo.F90:1959
compute_filter_constants
subroutine compute_filter_constants
Compute resolution-dependent values for the filtering.
Definition: filter_topo.F90:1975
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:1270