FV3DYCORE  Version 1.1.0
fv_surf_map.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the FV3 dynamical core.
5 !*
6 !* The FV3 dynamical core is free software: you can redistribute it
7 !* and/or modify it under the terms of the
8 !* GNU Lesser General Public License as published by the
9 !* Free Software Foundation, either version 3 of the License, or
10 !* (at your option) any later version.
11 !*
12 !* The FV3 dynamical core is distributed in the hope that it will be
13 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
14 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15 !* See the GNU General Public License for more details.
16 !*
17 !* You should have received a copy of the GNU Lesser General Public
18 !* License along with the FV3 dynamical core.
19 !* If not, see <http://www.gnu.org/licenses/>.
20 !***********************************************************************
21 
23 
24 ! <table>
25 ! <tr>
26 ! <th>Module Name</th>
27 ! <th>Functions Included</th>
28 ! </tr>
29 ! <tr>
30 ! <td>constants_mod</td>
31 ! <td>grav, radius, pi=>pi_8</td>
32 ! </tr>
33 ! <tr>
34 ! <td>fms_mod</td>
35 ! <td>file_exist, check_nml_error,open_namelist_file, close_file,
36 ! stdlog, mpp_pe, mpp_root_pe, FATAL, error_mesg</td>
37 ! </tr>
38 ! <tr>
39 ! <td>fv_arrays_mod</td>
40 ! <td>fv_grid_bounds_type, R_GRID</td>
41 ! </tr>
42 ! <tr>
43 ! <td>fv_grid_utils_mod</td>
44 ! <td>great_circle_dist, latlon2xyz, v_prod, normalize_vect,
45 ! g_sum, global_mx, vect_cross</td>
46 ! </tr>
47 ! <tr>
48 ! <td>fv_mp_mod</td>
49 ! <td>ng,mp_stop, mp_reduce_min, mp_reduce_max, is_master</td>
50 ! </tr>
51 ! <tr>
52 ! <td>fv_timing_mod</td>
53 ! <td>timing_on, timing_off</td>
54 ! </tr>
55 ! <tr>
56 ! <td>mpp_mod</td>
57 ! <td>get_unit, input_nml_file, mpp_error,
58 ! mpp_pe, mpp_chksum, stdout</td>
59 ! </tr>
60 ! <tr>
61 ! <td>mpp_domains_mod</td>
62 ! <td>mpp_update_domains, domain2d</td>
63 ! </tr>
64 ! </table>
65 
66  use fms_mod, only: file_exist, check_nml_error, &
67  open_namelist_file, close_file, stdlog, &
68  mpp_pe, mpp_root_pe, fatal, error_mesg
69  use mpp_mod, only: get_unit, input_nml_file, mpp_error
70  use mpp_domains_mod, only: mpp_update_domains, domain2d
71  use constants_mod, only: grav, radius, pi=>pi_8
72 
75  use fv_mp_mod, only: ng
76  use fv_mp_mod, only: mp_stop, mp_reduce_min, mp_reduce_max, is_master
79 
80  implicit none
81 
82  private
83 !-----------------------------------------------------------------------
84 ! NAMELIST
85 ! Name, resolution, and format of XXmin USGS datafile
86 ! 1min ---------> 1.85 km
87 ! nlon = 10800 * 2
88 ! nlat = 5400 * 2
89 ! 2min ---------> 3.7 km
90 ! nlon = 10800
91 ! nlat = 5400
92 ! 5min
93 ! nlon = 4320
94 ! nlat = 2160
95 ! surf_format: netcdf (default)
96 ! binary
97 ! New NASA SRTM30 data: SRTM30.nc
98 ! nlon = 43200
99 ! nlat = 21600
100  logical:: zs_filter = .true.
101  logical:: zero_ocean = .true. ! if true, no diffusive flux into water/ocean area
102  integer :: nlon = 21600
103  integer :: nlat = 10800
104  real:: cd4 = 0.15
105  real:: cd2 = -1.
106  real:: peak_fac = 1.05
107  real:: max_slope = 0.15
109  integer:: n_del2_weak = 12
110  integer:: n_del2_strong = -1
111  integer:: n_del4 = -1
112 
113 
114  character(len=128):: surf_file = "INPUT/topo1min.nc"
115  character(len=6) :: surf_format = 'netcdf'
116  logical :: namelist_read = .false.
117 
118  real(kind=R_GRID) da_min
119  real cos_grid
120  character(len=3) :: grid_string = ''
121 
122  namelist /surf_map_nml/ surf_file,surf_format,nlon,nlat, zero_ocean, zs_filter, &
124 !
125  real, allocatable:: zs_g(:,:), sgh_g(:,:), oro_g(:,:)
126 
127  public sgh_g, oro_g, zs_g
128  public surfdrv
130 
131  contains
132 
133  subroutine surfdrv(npx, npy, grid, agrid, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, phis, &
134  stretch_fac, nested, npx_global, domain,grid_number, bd, regional)
136  implicit none
137 #include <netcdf.inc>
138  integer, intent(in):: npx, npy
139 
140  ! INPUT arrays
141  type(fv_grid_bounds_type), intent(IN) :: bd
142  real(kind=R_GRID), intent(in)::area(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
143  real, intent(in):: dx(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng+1)
144  real, intent(in):: dy(bd%is-ng:bd%ie+ng+1, bd%js-ng:bd%je+ng)
145  real, intent(in), dimension(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)::dxa, dya
146  real, intent(in)::dxc(bd%is-ng:bd%ie+ng+1, bd%js-ng:bd%je+ng)
147  real, intent(in)::dyc(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng+1)
148 
149  real(kind=R_GRID), intent(in):: grid(bd%is-ng:bd%ie+ng+1, bd%js-ng:bd%je+ng+1,2)
150  real(kind=R_GRID), intent(in):: agrid(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng,2)
151  real, intent(IN):: sin_sg(bd%isd:bd%ied,bd%jsd:bd%jed,9)
152  real(kind=R_GRID), intent(IN):: stretch_fac
153  logical, intent(IN) :: nested, regional
154  integer, intent(IN) :: npx_global
155  type(domain2d), intent(INOUT) :: domain
156  integer, intent(IN) :: grid_number
157 
158  ! OUTPUT arrays
159  real, intent(out):: phis(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
160 ! Local:
161  real, allocatable :: z2(:,:)
162 ! Position of edges of the box containing the original data point:
163  integer londim
164  integer latdim
165  character(len=80) :: topoflnm
166  real(kind=4), allocatable :: ft(:,:), zs(:,:)
167  real, allocatable :: lon1(:), lat1(:)
168  real dx1, dx2, dy1, dy2, lats, latn, r2d
169  real(kind=R_GRID) da_max
170  real zmean, z2mean, delg, rgrav
171 ! real z_sp, f_sp, z_np, f_np
172  integer i, j, n, mdim
173  integer igh, jt
174  integer ncid, lonid, latid, ftopoid, htopoid
175  integer jstart, jend, start(4), nread(4)
176  integer status
177 
178  integer :: is, ie, js, je
179  integer :: isd, ied, jsd, jed
180  real phis_coarse(bd%isd:bd%ied, bd%jsd:bd%jed)
181  real wt
182 
183  is = bd%is
184  ie = bd%ie
185  js = bd%js
186  je = bd%je
187  isd = bd%isd
188  ied = bd%ied
189  jsd = bd%jsd
190  jed = bd%jed
191  if (nested) then
192  !Divide all by grav
193  rgrav = 1./grav
194  do j=jsd,jed
195  do i=isd,ied
196  phis(i,j) = phis(i,j)*rgrav
197  enddo
198  enddo
199  !Save interpolated coarse-grid data for blending
200  do j=jsd,jed
201  do i=isd,ied
202  phis_coarse(i,j) = phis(i,j)
203  enddo
204  enddo
205  endif
206 
207  do j=js,je
208  do i=is,ie
209  phis(i,j) = 0.0
210  enddo
211  enddo
212 
213  call read_namelist
214 
215  if (grid_number > 1) write(grid_string, '(A, I1)') ' g', grid_number
216 
217 
218 !
219 ! surface file must be in NetCDF format
220 !
221  if ( file_exist(surf_file) ) then
222  if (surf_format == "netcdf") then
223 
224  status = nf_open(surf_file, nf_nowrite, ncid)
225  if (status .ne. nf_noerr) call handle_err(status)
226 
227  status = nf_inq_dimid(ncid, 'lon', lonid)
228  if (status .ne. nf_noerr) call handle_err(status)
229  status = nf_inq_dimlen(ncid, lonid, londim)
230  if (status .ne. nf_noerr) call handle_err(status)
231  nlon = londim
232 
233  status = nf_inq_dimid(ncid, 'lat', latid)
234  if (status .ne. nf_noerr) call handle_err(status)
235  status = nf_inq_dimlen(ncid, latid, latdim)
236  if (status .ne. nf_noerr) call handle_err(status)
237  nlat = latdim
238 
239  if ( is_master() ) then
240  if ( nlon==43200 ) then
241  write(*,*) 'Opening NASA datset file:', surf_file, surf_format, nlon, nlat
242  else
243  write(*,*) 'Opening USGS datset file:', surf_file, surf_format, nlon, nlat
244  endif
245  endif
246 
247  else
248  call error_mesg ( 'surfdrv','Raw IEEE data format no longer supported !!!', fatal )
249  endif
250  else
251  call error_mesg ( 'surfdrv','surface file '//trim(surf_file)//' not found !', fatal )
252  endif
253 
254  allocate ( lat1(nlat+1) )
255  allocate ( lon1(nlon+1) )
256 
257  r2d = 180./pi
258 
259  cos_grid = cos( 2.*pi/real(nlat) ) ! two-data_grid distance
260 
261  dx1 = 2.*pi/real(nlon)
262  dy1 = pi/real(nlat)
263 
264  do i=1,nlon+1
265  lon1(i) = dx1 * real(i-1) ! between 0 2pi
266  enddo
267 
268  lat1(1) = - 0.5*pi
269  lat1(nlat+1) = 0.5*pi
270  do j=2,nlat
271  lat1(j) = -0.5*pi + dy1*(j-1)
272  enddo
273 
274 !-------------------------------------
275 ! Compute raw phis and oro
276 !-------------------------------------
277  call timing_on('map_to_cubed')
278 
279  if (surf_format == "netcdf") then
280 
281 ! Find latitude strips reading data
282  lats = pi/2.
283  latn = -pi/2.
284  do j=js,je
285  do i=is,ie
286  lats = min( lats, grid(i,j,2), grid(i+1,j,2), grid(i,j+1,2), grid(i+1,j+1,2), agrid(i,j,2) )
287  latn = max( latn, grid(i,j,2), grid(i+1,j,2), grid(i,j+1,2), grid(i+1,j+1,2), agrid(i,j,2) )
288  enddo
289  enddo
290 
291 ! Enlarge the search zone:
292 ! To account for the curvature of the coordinates:
293 
294  !I have had trouble running c90 with 600 pes unless the search region is expanded
295  ! due to failures in finding latlon points in the source data.
296  !This sets a larger search region if the number of cells on a PE is too small.
297  !(Alternately you can just cold start the topography using a smaller number of PEs)
298  if (min(je-js+1,ie-is+1) < 15) then
299  delg = max( 0.4*(latn-lats), pi/real(npx_global-1), 2.*pi/real(nlat) )
300  else
301  delg = max( 0.2*(latn-lats), pi/real(npx_global-1), 2.*pi/real(nlat) )
302  endif
303  lats = max( -0.5*pi, lats - delg )
304  latn = min( 0.5*pi, latn + delg )
305 
306  jstart = 1
307  do j=2,nlat
308  if ( lats < lat1(j) ) then
309  jstart = j-1
310  exit
311  endif
312  enddo
313  jstart = max(jstart-1, 1)
314 
315  jend = nlat
316  do j=2,nlat
317  if ( latn < lat1(j+1) ) then
318  jend = j+1
319  exit
320  endif
321  enddo
322  jend = min(jend+1, nlat)
323 
324  jt = jend - jstart + 1
325  igh = nlon/8 + nlon/(2*(npx_global-1))
326 
327  if (is_master()) write(*,*) 'Terrain dataset =', nlon, 'jt=', jt
328  if (is_master()) write(*,*) 'igh (terrain ghosting)=', igh
329 
330  status = nf_inq_varid(ncid, 'ftopo', ftopoid)
331  if (status .ne. nf_noerr) call handle_err(status)
332  nread = 1; start = 1
333  nread(1) = nlon
334  start(2) = jstart; nread(2) = jend - jstart + 1
335 
336  allocate ( ft(-igh:nlon+igh,jt) )
337  status = nf_get_vara_real(ncid, ftopoid, start, nread, ft(1:nlon,1:jt))
338  if (status .ne. nf_noerr) call handle_err(status)
339 
340  do j=1,jt
341  do i=-igh,0
342  ft(i,j) = ft(i+nlon,j)
343  enddo
344  do i=nlon+1,nlon+igh
345  ft(i,j) = ft(i-nlon,j)
346  enddo
347  enddo
348 
349  status = nf_inq_varid(ncid, 'htopo', htopoid)
350  if (status .ne. nf_noerr) call handle_err(status)
351  allocate ( zs(-igh:nlon+igh,jt) )
352  status = nf_get_vara_real(ncid, htopoid, start, nread, zs(1:nlon,1:jt))
353  if (status .ne. nf_noerr) call handle_err(status)
354  status = nf_close(ncid)
355  if (status .ne. nf_noerr) call handle_err(status)
356 ! Ghost Data
357  do j=1,jt
358  do i=-igh,0
359  zs(i,j) = zs(i+nlon,j)
360  enddo
361  do i=nlon+1,nlon+igh
362  zs(i,j) = zs(i-nlon,j)
363  enddo
364  enddo
365 
366  endif
367 
368 ! special SP treatment:
369 ! if ( jstart == 1 ) then
370 ! call zonal_mean(nlon, zs(1,1), z_sp)
371 ! call zonal_mean(nlon, ft(1,1), f_sp)
372 ! endif
373 
374  allocate ( oro_g(isd:ied, jsd:jed) )
375  allocate ( sgh_g(isd:ied, jsd:jed) )
376  call timing_on('map_to_cubed')
377  call map_to_cubed_raw(igh, nlon, jt, lat1(jstart:jend+1), lon1, zs, ft, grid, agrid, &
378  phis, oro_g, sgh_g, npx, npy, jstart, jend, stretch_fac, nested, npx_global, bd, regional)
379  if (is_master()) write(*,*) 'map_to_cubed_raw: master PE done'
380  call timing_off('map_to_cubed')
381 
382  deallocate ( zs )
383  deallocate ( ft )
384  deallocate ( lon1 )
385  deallocate ( lat1 )
386 
387  allocate ( zs_g(is:ie, js:je) )
388  allocate ( z2(is:ie,js:je) )
389  do j=js,je
390  do i=is,ie
391  zs_g(i,j) = phis(i,j)
392  z2(i,j) = phis(i,j)**2
393  enddo
394  enddo
395 !--------
396 ! Filter:
397 !--------
398  call global_mx(real(phis,kind=R_GRID), ng, da_min, da_max, bd)
399  zmean = g_sum(domain, zs_g(is:ie,js:je), is, ie, js, je, ng, area, 1)
400  z2mean = g_sum(domain, z2(is:ie,js:je) , is, ie, js, je, ng, area, 1)
401  if ( is_master() ) then
402  write(*,*) 'Before filter ZS', trim(grid_string), ' min=', da_min, ' Max=', da_max,' Mean=',zmean
403  write(*,*) '*** Mean variance', trim(grid_string), ' *** =', z2mean
404  endif
405 
406  !On a nested grid blend coarse-grid and nested-grid
407  ! orography near boundary.
408  ! This works only on the height of topography; assume
409  ! land fraction and sub-grid variance unchanged
410 
411  ! Here, we blend in the four cells nearest to the boundary.
412  ! In the halo we set the value to that interpolated from the coarse
413  ! grid. (Previously this was erroneously not being done, which was causing
414  ! the halo to be filled with unfiltered nested-grid terrain, creating
415  ! ugly edge artifacts.)
416  if (nested) then
417  if (is_master()) write(*,*) 'Blending nested and coarse grid topography'
418  do j=jsd,jed
419  do i=isd,ied
420  wt = max(0.,min(1.,real(5 - min(i,j,npx-i,npy-j,5))/5. ))
421  phis(i,j) = (1.-wt)*phis(i,j) + wt*phis_coarse(i,j)
422 
423  enddo
424  enddo
425  endif
426 
427  call global_mx(real(oro_g,kind=R_GRID), ng, da_min, da_max, bd)
428  if ( is_master() ) write(*,*) 'ORO', trim(grid_string), ' min=', da_min, ' Max=', da_max
429 
430  call global_mx(area, ng, da_min, da_max, bd)
431  call timing_on('Terrain_filter')
432 ! Del-2: high resolution only
433  if ( zs_filter ) then
434  if(is_master()) then
435  write(*,*) 'Applying terrain filters. zero_ocean is', zero_ocean
436  endif
437  call fv3_zs_filter (bd, isd, ied, jsd, jed, npx, npy, npx_global, &
438  stretch_fac, nested, domain, area, dxa, dya, dx, dy, dxc, dyc, grid, &
439  agrid, sin_sg, phis, oro_g, regional)
440  call mpp_update_domains(phis, domain)
441  endif ! end terrain filter
442  call timing_off('Terrain_filter')
443 
444  do j=js,je
445  do i=is,ie
446  z2(i,j) = phis(i,j)**2
447  end do
448  end do
449 
450  call global_mx(real(phis,kind=R_GRID), ng, da_min, da_max, bd)
451  zmean = g_sum(domain, phis(is:ie,js:je), is, ie, js, je, ng, area, 1)
452  z2mean = g_sum(domain, z2, is, ie, js, je, ng, area, 1)
453  deallocate ( z2 )
454 
455  if ( is_master() ) then
456  write(*,*) 'After filter Phis', trim(grid_string), ' min=', da_min, ' Max=', da_max, 'Mean=', zmean
457  write(*,*) '*** Mean variance', trim(grid_string), ' *** =', z2mean
458  endif
459 
460  !FOR NESTING: Unless we fill the outermost halo with topography
461  ! interpolated from the coarse grid, the values of phi there
462  ! will be wrong because they are just z instead of g*z; they
463  ! have not had gravity properly multiplied in so we can get phi
464 
465  ! For now we compute phis and sgh_g on the full data domain; for
466  ! nested grids this allows us to do the smoothing near the boundary
467  ! without having to fill the boundary halo from the coarse grid
468 
469  !ALSO for nesting: note that we are smoothing the terrain using
470  ! the nested-grid's outer halo filled with the terrain computed
471  ! directly from the input file computed here, and then
472  ! replacing it with interpolated topography in fv_restart, so
473  ! as to be consistent when doing the boundary condition
474  ! interpolation. We would ideally replace the nested-grid
475  ! halo topography BEFORE smoothing, which could be more
476  ! consistent, but this would require moving calls to
477  ! nested_grid_BC in this routine.
478 
479  do j=jsd,jed
480  do i=isd,ied
481  phis(i,j) = grav * phis(i,j)
482  if ( sgh_g(i,j) <= 0. ) then
483  sgh_g(i,j) = 0.
484  else
485  sgh_g(i,j) = sqrt(sgh_g(i,j))
486  endif
487  end do
488  end do
489 
490  call global_mx(real(sgh_g,kind=R_GRID), ng, da_min, da_max, bd)
491  if ( is_master() ) write(*,*) 'Before filter SGH', trim(grid_string), ' min=', da_min, ' Max=', da_max
492 
493 
494 !-----------------------------------------------
495 ! Filter the standard deviation of mean terrain:
496 !-----------------------------------------------
497  call global_mx(area, ng, da_min, da_max, bd)
498 
499  if(zs_filter) call del4_cubed_sphere(npx, npy, sgh_g, area, dx, dy, dxc, dyc, sin_sg, 1, zero_ocean, oro_g, nested, domain, bd, regional)
500 
501  call global_mx(real(sgh_g,kind=R_GRID), ng, da_min, da_max, bd)
502  if ( is_master() ) write(*,*) 'After filter SGH', trim(grid_string), ' min=', da_min, ' Max=', da_max
503  do j=js,je
504  do i=is,ie
505  sgh_g(i,j) = max(0., sgh_g(i,j))
506  enddo
507  enddo
508 
509  end subroutine surfdrv
510 
511  subroutine fv3_zs_filter (bd, isd, ied, jsd, jed, npx, npy, npx_global, &
512  stretch_fac, nested, domain, area, dxa, dya, dx, dy, dxc, dyc, grid, &
513  agrid, sin_sg, phis, oro ,regional)
514  integer, intent(in):: isd, ied, jsd, jed, npx, npy, npx_global
515  type(fv_grid_bounds_type), intent(IN) :: bd
516  real(kind=R_GRID), intent(in), dimension(isd:ied,jsd:jed)::area
517  real, intent(in), dimension(isd:ied,jsd:jed)::dxa, dya
518  real, intent(in), dimension(isd:ied, jsd:jed+1):: dx, dyc
519  real, intent(in), dimension(isd:ied+1,jsd:jed):: dy, dxc
520 
521  real(kind=R_GRID), intent(in):: grid(isd:ied+1, jsd:jed+1,2)
522  real(kind=R_GRID), intent(in):: agrid(isd:ied, jsd:jed, 2)
523  real, intent(IN):: sin_sg(isd:ied,jsd:jed,9)
524  real(kind=R_GRID), intent(IN):: stretch_fac
525  logical, intent(IN) :: nested, regional
526  real, intent(inout):: phis(isd:ied,jsd,jed)
527  real, intent(inout):: oro(isd:ied,jsd,jed)
528  type(domain2d), intent(INOUT) :: domain
529  integer mdim
530  real(kind=R_GRID) da_max
531 
532  if (is_master()) print*, ' Calling FV3_zs_filter...'
533 
534  if (.not. namelist_read) call read_namelist !when calling from external_ic
535  call global_mx(area, ng, da_min, da_max, bd)
536 
537  mdim = nint( real(npx_global) * min(10., stretch_fac) )
538 
539 ! Del-2: high resolution only
540 ! call del2_cubed_sphere(npx, npy, phis, area, dx, dy, dxc, dyc, sin_sg, n_del2, cd2, zero_ocean, oro, nested, domain, bd)
541  if (n_del2_strong < 0) then
542  if ( npx_global<=97) then
543  n_del2_strong = 0
544  elseif ( npx_global<=193 ) then
545  n_del2_strong = 1
546  else
547  n_del2_strong = 2
548  endif
549  endif
550  if (cd2 < 0.) cd2 = 0.16*da_min
551 ! Applying strong 2-delta-filter:
552  if ( n_del2_strong > 0 ) &
553  call two_delta_filter(npx, npy, phis, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd2, zero_ocean, &
554  .true., 0, oro, nested, domain, bd, n_del2_strong, regional)
555 
556 ! MFCT Del-4:
557  if (n_del4 < 0) then
558  if ( mdim<=193 ) then
559  n_del4 = 1
560  elseif ( mdim<=1537 ) then
561  n_del4 = 2
562  else
563  n_del4 = 3
564  endif
565  endif
566  call del4_cubed_sphere(npx, npy, phis, area, dx, dy, dxc, dyc, sin_sg, n_del4, zero_ocean, oro, nested, domain, bd, regional)
567 ! Applying weak 2-delta-filter:
568  cd2 = 0.12*da_min
569  call two_delta_filter(npx, npy, phis, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd2, zero_ocean, &
570  .true., 1, oro, nested, domain, bd, n_del2_weak, regional)
571 
572 
573  end subroutine fv3_zs_filter
574 
575 
576  subroutine two_delta_filter(npx, npy, q, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd, zero_ocean, &
577  check_slope, filter_type, oro, nested, domain, bd, ntmax, regional)
578  type(fv_grid_bounds_type), intent(IN) :: bd
579  integer, intent(in):: npx, npy
580  integer, intent(in):: ntmax
581  integer, intent(in):: filter_type
582  real, intent(in):: cd
583 ! INPUT arrays
584  real(kind=R_GRID), intent(in)::area(bd%isd:bd%ied, bd%jsd:bd%jed)
585  real, intent(in):: dx(bd%isd:bd%ied, bd%jsd:bd%jed+1)
586  real, intent(in):: dy(bd%isd:bd%ied+1,bd%jsd:bd%jed)
587  real, intent(in):: dxa(bd%isd:bd%ied, bd%jsd:bd%jed)
588  real, intent(in):: dya(bd%isd:bd%ied, bd%jsd:bd%jed)
589  real, intent(in):: dxc(bd%isd:bd%ied+1,bd%jsd:bd%jed)
590  real, intent(in):: dyc(bd%isd:bd%ied, bd%jsd:bd%jed+1)
591  real, intent(in):: sin_sg(bd%isd:bd%ied,bd%jsd:bd%jed,9)
592  real, intent(in):: oro(bd%isd:bd%ied, bd%jsd:bd%jed)
593  logical, intent(in):: zero_ocean, check_slope
594  logical, intent(in):: nested, regional
595  type(domain2d), intent(inout) :: domain
596 ! OUTPUT arrays
597  real, intent(inout):: q(bd%isd:bd%ied, bd%jsd:bd%jed)
598 ! Local:
599  real, parameter:: p1 = 7./12.
600  real, parameter:: p2 = -1./12.
601  real, parameter:: c1 = -2./14.
602  real, parameter:: c2 = 11./14.
603  real, parameter:: c3 = 5./14.
604  real:: ddx(bd%is:bd%ie+1,bd%js:bd%je), ddy(bd%is:bd%ie,bd%js:bd%je+1)
605  logical:: extm(bd%is-1:bd%ie+1)
606  logical:: ext2(bd%is:bd%ie,bd%js-1:bd%je+1)
607  real:: a1(bd%is-1:bd%ie+2)
608  real:: a2(bd%is:bd%ie,bd%js-1:bd%je+2)
609  real(kind=R_GRID):: a3(bd%is:bd%ie,bd%js:bd%je)
610  real(kind=R_GRID):: smax, smin
611  real:: m_slope, fac
612  integer:: i,j, nt
613  integer:: is, ie, js, je
614  integer:: isd, ied, jsd, jed
615  integer:: is1, ie2, js1, je2
616 
617  is = bd%is
618  ie = bd%ie
619  js = bd%js
620  je = bd%je
621  isd = bd%isd
622  ied = bd%ied
623  jsd = bd%jsd
624  jed = bd%jed
625 
626  if ( nested ) then
627  is1 = is-1; ie2 = ie+2
628  js1 = js-1; je2 = je+2
629  else
630  is1 = max(3,is-1); ie2 = min(npx-2,ie+2)
631  js1 = max(3,js-1); je2 = min(npy-2,je+2)
632  end if
633 
634  if ( check_slope ) then
635  m_slope = max_slope
636  else
637  m_slope = 10.
638  endif
639 
640 
641  do 777 nt=1, ntmax
642  call mpp_update_domains(q, domain)
643 
644 ! Check slope
645  if ( nt==1 .and. check_slope ) then
646  do j=js,je
647  do i=is,ie+1
648  ddx(i,j) = (q(i,j) - q(i-1,j))/dxc(i,j)
649  ddx(i,j) = abs(ddx(i,j))
650  enddo
651  enddo
652  do j=js,je+1
653  do i=is,ie
654  ddy(i,j) = (q(i,j) - q(i,j-1))/dyc(i,j)
655  ddy(i,j) = abs(ddy(i,j))
656  enddo
657  enddo
658  do j=js,je
659  do i=is,ie
660  a3(i,j) = max( ddx(i,j), ddx(i+1,j), ddy(i,j), ddy(i,j+1) )
661  enddo
662  enddo
663  call global_mx(a3, 0, smin, smax, bd)
664  if ( is_master() ) write(*,*) 'Before filter: Max_slope=', smax
665  endif
666 
667 ! First step: average the corners:
668  if ( .not. (nested .or. regional) .and. nt==1 ) then
669  if ( is==1 .and. js==1 ) then
670  q(1,1) = (q(1,1)*area(1,1)+q(0,1)*area(0,1)+q(1,0)*area(1,0)) &
671  / ( area(1,1)+ area(0,1)+ area(1,0) )
672  q(0,1) = q(1,1)
673  q(1,0) = q(1,1)
674  endif
675  if ( (ie+1)==npx .and. js==1 ) then
676  q(ie, 1) = (q(ie,1)*area(ie,1)+q(npx,1)*area(npx,1)+q(ie,0)*area(ie,0)) &
677  / ( area(ie,1)+ area(npx,1)+ area(ie,0))
678  q(npx,1) = q(ie,1)
679  q(ie, 0) = q(ie,1)
680  endif
681  if ( is==1 .and. (je+1)==npy ) then
682  q(1, je) = (q(1,je)*area(1,je)+q(0,je)*area(0,je)+q(1,npy)*area(1,npy)) &
683  / ( area(1,je)+ area(0,je)+ area(1,npy))
684  q(0, je) = q(1,je)
685  q(1,npy) = q(1,je)
686  endif
687  if ( (ie+1)==npx .and. (je+1)==npy ) then
688  q(ie, je) = (q(ie,je)*area(ie,je)+q(npx,je)*area(npx,je)+q(ie,npy)*area(ie,npy)) &
689  / ( area(ie,je)+ area(npx,je)+ area(ie,npy))
690  q(npx,je) = q(ie,je)
691  q(ie,npy) = q(ie,je)
692  endif
693  call mpp_update_domains(q, domain)
694  endif
695 
696 ! x-diffusive flux:
697  do 333 j=js,je
698 
699  do i=is1, ie2
700  a1(i) = p1*(q(i-1,j)+q(i,j)) + p2*(q(i-2,j)+q(i+1,j))
701  enddo
702 
703  if ( .not. (nested .or. regional) ) then
704  if ( is==1 ) then
705  a1(0) = c1*q(-2,j) + c2*q(-1,j) + c3*q(0,j)
706  a1(1) = 0.5*(((2.*dxa(0,j)+dxa(-1,j))*q(0,j)-dxa(0,j)*q(-1,j))/(dxa(-1,j)+dxa(0,j)) &
707  + ((2.*dxa(1,j)+dxa( 2,j))*q(1,j)-dxa(1,j)*q( 2,j))/(dxa(1, j)+dxa(2,j)))
708  a1(2) = c3*q(1,j) + c2*q(2,j) +c1*q(3,j)
709  endif
710  if ( (ie+1)==npx ) then
711  a1(npx-1) = c1*q(npx-3,j) + c2*q(npx-2,j) + c3*q(npx-1,j)
712  a1(npx) = 0.5*(((2.*dxa(npx-1,j)+dxa(npx-2,j))*q(npx-1,j)-dxa(npx-1,j)*q(npx-2,j))/(dxa(npx-2,j)+dxa(npx-1,j)) &
713  + ((2.*dxa(npx, j)+dxa(npx+1,j))*q(npx, j)-dxa(npx, j)*q(npx+1,j))/(dxa(npx, j)+dxa(npx+1,j)))
714  a1(npx+1) = c3*q(npx,j) + c2*q(npx+1,j) + c1*q(npx+2,j)
715  endif
716  endif
717 
718  if ( filter_type == 0 ) then
719  do i=is-1, ie+1
720  if( abs(3.*(a1(i)+a1(i+1)-2.*q(i,j))) > abs(a1(i)-a1(i+1)) ) then
721  extm(i) = .true.
722  else
723  extm(i) = .false.
724  endif
725  enddo
726  else
727  do i=is-1, ie+1
728  if ( (a1(i)-q(i,j))*(a1(i+1)-q(i,j)) > 0. ) then
729  extm(i) = .true.
730  else
731  extm(i) = .false.
732  endif
733  enddo
734  endif
735 
736  do i=is,ie+1
737  ddx(i,j) = (q(i-1,j)-q(i,j))/dxc(i,j)
738  if ( extm(i-1).and.extm(i) ) then
739  ddx(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*ddx(i,j)
740  elseif ( abs(ddx(i,j)) > m_slope ) then
741  fac = min(1., max(0.1,(abs(ddx(i,j))-m_slope)/m_slope ) )
742  ddx(i,j) = fac*0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*ddx(i,j)
743  else
744  ddx(i,j) = 0.
745  endif
746  enddo
747 333 continue
748 
749 ! y-diffusive flux:
750  do j=js1,je2
751  do i=is,ie
752  a2(i,j) = p1*(q(i,j-1)+q(i,j)) + p2*(q(i,j-2)+q(i,j+1))
753  enddo
754  enddo
755  if ( .not. (nested .or. regional) ) then
756  if( js==1 ) then
757  do i=is,ie
758  a2(i,0) = c1*q(i,-2) + c2*q(i,-1) + c3*q(i,0)
759  a2(i,1) = 0.5*(((2.*dya(i,0)+dya(i,-1))*q(i,0)-dya(i,0)*q(i,-1))/(dya(i,-1)+dya(i,0)) &
760  + ((2.*dya(i,1)+dya(i, 2))*q(i,1)-dya(i,1)*q(i, 2))/(dya(i, 1)+dya(i,2)))
761  a2(i,2) = c3*q(i,1) + c2*q(i,2) + c1*q(i,3)
762  enddo
763  endif
764  if( (je+1)==npy ) then
765  do i=is,ie
766  a2(i,npy-1) = c1*q(i,npy-3) + c2*q(i,npy-2) + c3*q(i,npy-1)
767  a2(i,npy) = 0.5*(((2.*dya(i,npy-1)+dya(i,npy-2))*q(i,npy-1)-dya(i,npy-1)*q(i,npy-2))/(dya(i,npy-2)+dya(i,npy-1)) &
768  + ((2.*dya(i,npy)+dya(i,npy+1))*q(i,npy)-dya(i,npy)*q(i,npy+1))/(dya(i,npy)+dya(i,npy+1)))
769  a2(i,npy+1) = c3*q(i,npy) + c2*q(i,npy+1) + c1*q(i,npy+2)
770  enddo
771  endif
772  endif
773 
774  if ( filter_type == 0 ) then
775  do j=js-1,je+1
776  do i=is,ie
777  if( abs(3.*(a2(i,j)+a2(i,j+1)-2.*q(i,j))) > abs(a2(i,j)-a2(i,j+1)) ) then
778  ext2(i,j) = .true.
779  else
780  ext2(i,j) = .false.
781  endif
782  enddo
783  enddo
784  else
785  do j=js-1,je+1
786  do i=is,ie
787  if ( (a2(i,j)-q(i,j))*(a2(i,j+1)-q(i,j)) > 0. ) then
788  ext2(i,j) = .true.
789  else
790  ext2(i,j) = .false.
791  endif
792  enddo
793  enddo
794  endif
795 
796  do j=js,je+1
797  do i=is,ie
798  ddy(i,j) = (q(i,j-1)-q(i,j))/dyc(i,j)
799  if ( ext2(i,j-1) .and. ext2(i,j) ) then
800  ddy(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*ddy(i,j)
801  elseif ( abs(ddy(i,j))>m_slope ) then
802  fac = min(1., max(0.1,(abs(ddy(i,j))-m_slope)/m_slope))
803  ddy(i,j) = fac*0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*ddy(i,j)
804  else
805  ddy(i,j) = 0.
806  endif
807  enddo
808  enddo
809 
810  if ( zero_ocean ) then
811 ! Limit diffusive flux over water cells:
812  do j=js,je
813  do i=is,ie+1
814  ddx(i,j) = max(0., min(oro(i-1,j), oro(i,j))) * ddx(i,j)
815  enddo
816  enddo
817  do j=js,je+1
818  do i=is,ie
819  ddy(i,j) = max(0., min(oro(i,j-1), oro(i,j))) * ddy(i,j)
820  enddo
821  enddo
822  endif
823 
824  do j=js,je
825  do i=is,ie
826  q(i,j) = q(i,j) + cd/area(i,j)*(ddx(i,j)-ddx(i+1,j)+ddy(i,j)-ddy(i,j+1))
827  enddo
828  enddo
829 777 continue
830 
831 ! Check slope
832  if ( check_slope ) then
833  call mpp_update_domains(q, domain)
834  do j=js,je
835  do i=is,ie+1
836  ddx(i,j) = (q(i,j) - q(i-1,j))/dxc(i,j)
837  ddx(i,j) = abs(ddx(i,j))
838  enddo
839  enddo
840  do j=js,je+1
841  do i=is,ie
842  ddy(i,j) = (q(i,j) - q(i,j-1))/dyc(i,j)
843  ddy(i,j) = abs(ddy(i,j))
844  enddo
845  enddo
846  do j=js,je
847  do i=is,ie
848  a3(i,j) = max( ddx(i,j), ddx(i+1,j), ddy(i,j), ddy(i,j+1) )
849  enddo
850  enddo
851  call global_mx(a3, 0, smin, smax, bd)
852  if ( is_master() ) write(*,*) 'After filter: Max_slope=', smax
853  endif
854 
855  end subroutine two_delta_filter
856 
857 
858 
859  subroutine del2_cubed_sphere(npx, npy, q, area, dx, dy, dxc, dyc, sin_sg, nmax, cd, zero_ocean, oro, nested, domain, bd, regional)
860  type(fv_grid_bounds_type), intent(IN) :: bd
861  integer, intent(in):: npx, npy
862  integer, intent(in):: nmax
863  real(kind=R_GRID), intent(in):: cd
864  logical, intent(in):: zero_ocean
865  ! INPUT arrays
866  real(kind=R_GRID), intent(in)::area(bd%isd:bd%ied, bd%jsd:bd%jed)
867  real, intent(in):: dx(bd%isd:bd%ied, bd%jsd:bd%jed+1)
868  real, intent(in):: dy(bd%isd:bd%ied+1,bd%jsd:bd%jed)
869  real, intent(in):: dxc(bd%isd:bd%ied+1,bd%jsd:bd%jed)
870  real, intent(in):: dyc(bd%isd:bd%ied, bd%jsd:bd%jed+1)
871  real, intent(IN):: sin_sg(bd%isd:bd%ied,bd%jsd:bd%jed,9)
872  real, intent(in):: oro(bd%isd:bd%ied, bd%jsd:bd%jed)
873  logical, intent(IN) :: nested, regional
874  type(domain2d), intent(INOUT) :: domain
875  ! OUTPUT arrays
876  real, intent(inout):: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
877 ! Local:
878  real ddx(bd%is:bd%ie+1,bd%js:bd%je), ddy(bd%is:bd%ie,bd%js:bd%je+1)
879  integer i,j,n
880 
881  integer :: is, ie, js, je
882  integer :: isd, ied, jsd, jed
883 
884  is = bd%is
885  ie = bd%ie
886  js = bd%js
887  je = bd%je
888  isd = bd%isd
889  ied = bd%ied
890  jsd = bd%jsd
891  jed = bd%jed
892 
893 
894  call mpp_update_domains(q,domain,whalo=ng,ehalo=ng,shalo=ng,nhalo=ng)
895 
896 ! First step: average the corners:
897  if ( is==1 .and. js==1 .and. .not. (nested .or. regional)) then
898  q(1,1) = (q(1,1)*area(1,1)+q(0,1)*area(0,1)+q(1,0)*area(1,0)) &
899  / ( area(1,1)+ area(0,1)+ area(1,0) )
900  q(0,1) = q(1,1)
901  q(1,0) = q(1,1)
902  endif
903  if ( (ie+1)==npx .and. js==1 .and. .not. (nested .or. regional)) then
904  q(ie, 1) = (q(ie,1)*area(ie,1)+q(npx,1)*area(npx,1)+q(ie,0)*area(ie,0)) &
905  / ( area(ie,1)+ area(npx,1)+ area(ie,0))
906  q(npx,1) = q(ie,1)
907  q(ie, 0) = q(ie,1)
908  endif
909  if ( (ie+1)==npx .and. (je+1)==npy .and. .not. (nested .or. regional) ) then
910  q(ie, je) = (q(ie,je)*area(ie,je)+q(npx,je)*area(npx,je)+q(ie,npy)*area(ie,npy)) &
911  / ( area(ie,je)+ area(npx,je)+ area(ie,npy))
912  q(npx,je) = q(ie,je)
913  q(ie,npy) = q(ie,je)
914  endif
915  if ( is==1 .and. (je+1)==npy .and. .not. (nested .or. regional)) then
916  q(1, je) = (q(1,je)*area(1,je)+q(0,je)*area(0,je)+q(1,npy)*area(1,npy)) &
917  / ( area(1,je)+ area(0,je)+ area(1,npy))
918  q(0, je) = q(1,je)
919  q(1,npy) = q(1,je)
920  endif
921 
922 
923  do n=1,nmax
924  if( n>1 ) call mpp_update_domains(q,domain,whalo=ng,ehalo=ng,shalo=ng,nhalo=ng)
925  do j=js,je
926  do i=is,ie+1
927  ddx(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*(q(i-1,j)-q(i,j))/dxc(i,j)
928  enddo
929  enddo
930  do j=js,je+1
931  do i=is,ie
932  ddy(i,j) = dx(i,j)*(q(i,j-1)-q(i,j))/dyc(i,j) &
933  *0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))
934  enddo
935  enddo
936 
937  if ( zero_ocean ) then
938 ! Limit diffusive flux over ater cells:
939  do j=js,je
940  do i=is,ie+1
941  ddx(i,j) = max(0., min(oro(i-1,j), oro(i,j))) * ddx(i,j)
942  enddo
943  enddo
944  do j=js,je+1
945  do i=is,ie
946  ddy(i,j) = max(0., min(oro(i,j-1), oro(i,j))) * ddy(i,j)
947  enddo
948  enddo
949  endif
950 
951  do j=js,je
952  do i=is,ie
953  q(i,j) = q(i,j) + cd/area(i,j)*(ddx(i,j)-ddx(i+1,j)+ddy(i,j)-ddy(i,j+1))
954  enddo
955  enddo
956  enddo
957 
958  end subroutine del2_cubed_sphere
959 
960 
961  subroutine del4_cubed_sphere(npx, npy, q, area, dx, dy, dxc, dyc, sin_sg, nmax, zero_ocean, oro, nested, domain, bd, regional)
962  type(fv_grid_bounds_type), intent(IN) :: bd
963  integer, intent(in):: npx, npy, nmax
964  logical, intent(in):: zero_ocean
965  real, intent(in):: oro(bd%isd:bd%ied, bd%jsd:bd%jed)
966  real(kind=R_GRID), intent(in)::area(bd%isd:bd%ied, bd%jsd:bd%jed)
967  real, intent(in):: dx(bd%isd:bd%ied, bd%jsd:bd%jed+1)
968  real, intent(in):: dy(bd%isd:bd%ied+1,bd%jsd:bd%jed)
969  real, intent(in):: dxc(bd%isd:bd%ied+1,bd%jsd:bd%jed)
970  real, intent(in):: dyc(bd%isd:bd%ied, bd%jsd:bd%jed+1)
971  real, intent(IN):: sin_sg(bd%isd:bd%ied,bd%jsd:bd%jed,9)
972  real, intent(inout):: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
973  logical, intent(IN) :: nested, regional
974  type(domain2d), intent(INOUT) :: domain
975 ! diffusivity
976  real :: diff(bd%is-3:bd%ie+2,bd%js-3:bd%je+2)
977 ! diffusive fluxes:
978  real :: fx1(bd%is:bd%ie+1,bd%js:bd%je), fy1(bd%is:bd%ie,bd%js:bd%je+1)
979  real :: fx2(bd%is:bd%ie+1,bd%js:bd%je), fy2(bd%is:bd%ie,bd%js:bd%je+1)
980  real :: fx4(bd%is:bd%ie+1,bd%js:bd%je), fy4(bd%is:bd%ie,bd%js:bd%je+1)
981  real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: d2, win, wou
982  real, dimension(bd%is:bd%ie,bd%js:bd%je):: qlow, qmin, qmax, q0
983  real, parameter:: esl = 1.e-20
984  integer i,j, n
985 
986  integer :: is, ie, js, je
987  integer :: isd, ied, jsd, jed
988 
989  is = bd%is
990  ie = bd%ie
991  js = bd%js
992  je = bd%je
993  isd = bd%isd
994  ied = bd%ied
995  jsd = bd%jsd
996  jed = bd%jed
997 
998  !On a nested grid the haloes are not filled. Set to zero.
999  d2 = 0.
1000  win = 0.
1001  wou = 0.
1002 
1003  do j=js-1,je+1
1004  do i=is-1,ie+1
1005  diff(i,j) = cd4*area(i,j) ! area dependency is needed for stretched grid
1006  enddo
1007  enddo
1008 
1009  do j=js,je
1010  do i=is,ie
1011  q0(i,j) = q(i,j)
1012  enddo
1013  enddo
1014 
1015  do n=1,nmax
1016  call mpp_update_domains(q,domain)
1017 
1018 ! First step: average the corners:
1019  if ( is==1 .and. js==1 .and. .not. (nested .or. regional)) then
1020  q(1,1) = (q(1,1)*area(1,1)+q(0,1)*area(0,1)+q(1,0)*area(1,0)) &
1021  / ( area(1,1)+ area(0,1)+ area(1,0) )
1022  q(0,1) = q(1,1)
1023  q(1,0) = q(1,1)
1024  q(0,0) = q(1,1)
1025  endif
1026  if ( (ie+1)==npx .and. js==1 .and. .not. (nested .or. regional)) then
1027  q(ie, 1) = (q(ie,1)*area(ie,1)+q(npx,1)*area(npx,1)+q(ie,0)*area(ie,0)) &
1028  / ( area(ie,1)+ area(npx,1)+ area(ie,0))
1029  q(npx,1) = q(ie,1)
1030  q(ie, 0) = q(ie,1)
1031  q(npx,0) = q(ie,1)
1032  endif
1033  if ( (ie+1)==npx .and. (je+1)==npy .and. .not. (nested .or. regional)) then
1034  q(ie, je) = (q(ie,je)*area(ie,je)+q(npx,je)*area(npx,je)+q(ie,npy)*area(ie,npy)) &
1035  / ( area(ie,je)+ area(npx,je)+ area(ie,npy))
1036  q(npx, je) = q(ie,je)
1037  q(ie, npy) = q(ie,je)
1038  q(npx,npy) = q(ie,je)
1039  endif
1040  if ( is==1 .and. (je+1)==npy .and. .not. (nested .or. regional)) then
1041  q(1, je) = (q(1,je)*area(1,je)+q(0,je)*area(0,je)+q(1,npy)*area(1,npy)) &
1042  / ( area(1,je)+ area(0,je)+ area(1,npy))
1043  q(0, je) = q(1,je)
1044  q(1,npy) = q(1,je)
1045  q(0,npy) = q(1,je)
1046  endif
1047 
1048  do j=js,je
1049  do i=is,ie
1050  qmin(i,j) = min(q0(i,j), q(i-1,j-1), q(i,j-1), q(i+1,j-1), &
1051  q(i-1,j ), q(i,j ), q(i+1,j ), &
1052  q(i-1,j+1), q(i,j+1), q(i+1,j+1) )
1053  qmax(i,j) = max(peak_fac*q0(i,j), q(i-1,j-1), q(i,j-1), q(i+1,j-1), &
1054  q(i-1,j ), q(i,j ), q(i+1,j ), &
1055  q(i-1,j+1), q(i,j+1), q(i+1,j+1) )
1056  enddo
1057  enddo
1058 
1059 !--------------
1060 ! Compute del-2
1061 !--------------
1062 ! call copy_corners(q, npx, npy, 1)
1063  do j=js,je
1064  do i=is,ie+1
1065  fx2(i,j) = 0.25*(diff(i-1,j)+diff(i,j))*dy(i,j)*(q(i-1,j)-q(i,j))/dxc(i,j) &
1066  *(sin_sg(i,j,1)+sin_sg(i-1,j,3))
1067  enddo
1068  enddo
1069 
1070 ! call copy_corners(q, npx, npy, 2)
1071  do j=js,je+1
1072  do i=is,ie
1073  fy2(i,j) = 0.25*(diff(i,j-1)+diff(i,j))*dx(i,j)*(q(i,j-1)-q(i,j))/dyc(i,j) &
1074  *(sin_sg(i,j,2)+sin_sg(i,j-1,4))
1075  enddo
1076  enddo
1077 
1078  do j=js,je
1079  do i=is,ie
1080  d2(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1)) / area(i,j)
1081  enddo
1082  enddo
1083 
1084 ! qlow == low order monotonic solution
1085  if ( zero_ocean ) then
1086 ! Limit diffusive flux over ater cells:
1087  do j=js,je
1088  do i=is,ie+1
1089  fx1(i,j) = max(0., min(oro(i-1,j), oro(i,j))) * fx2(i,j)
1090  enddo
1091  enddo
1092  do j=js,je+1
1093  do i=is,ie
1094  fy1(i,j) = max(0., min(oro(i,j-1), oro(i,j))) * fy2(i,j)
1095  enddo
1096  enddo
1097  do j=js,je
1098  do i=is,ie
1099  qlow(i,j) = q(i,j) + (fx1(i,j)-fx1(i+1,j)+fy1(i,j)-fy1(i,j+1)) / area(i,j)
1100  d2(i,j) = diff(i,j) * d2(i,j)
1101  enddo
1102  enddo
1103  else
1104  do j=js,je
1105  do i=is,ie
1106  qlow(i,j) = q(i,j) + d2(i,j)
1107  d2(i,j) = diff(i,j) * d2(i,j)
1108  enddo
1109  enddo
1110  endif
1111 
1112  call mpp_update_domains(d2,domain)
1113 
1114 !---------------------
1115 ! Compute del4 fluxes:
1116 !---------------------
1117 ! call copy_corners(d2, npx, npy, 1)
1118  do j=js,je
1119  do i=is,ie+1
1120  fx4(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*(d2(i,j)-d2(i-1,j))/dxc(i,j)-fx2(i,j)
1121  enddo
1122  enddo
1123 
1124 ! call copy_corners(d2, npx, npy, 2)
1125  do j=js,je+1
1126  do i=is,ie
1127  fy4(i,j) = dx(i,j)*(d2(i,j)-d2(i,j-1))/dyc(i,j) &
1128  *0.5*(sin_sg(i,j,2)+sin_sg(i,j-1,4))-fy2(i,j)
1129  enddo
1130  enddo
1131 
1132 !----------------
1133 ! Flux limitting:
1134 !----------------
1135  do j=js,je
1136  do i=is,ie
1137  win(i,j) = max(0.,fx4(i, j)) - min(0.,fx4(i+1,j)) + &
1138  max(0.,fy4(i, j)) - min(0.,fy4(i,j+1)) + esl
1139  wou(i,j) = max(0.,fx4(i+1,j)) - min(0.,fx4(i, j)) + &
1140  max(0.,fy4(i,j+1)) - min(0.,fy4(i, j)) + esl
1141  win(i,j) = max(0., qmax(i,j) - qlow(i,j)) / win(i,j)*area(i,j)
1142  wou(i,j) = max(0., qlow(i,j) - qmin(i,j)) / wou(i,j)*area(i,j)
1143  enddo
1144  enddo
1145 
1146  call mpp_update_domains(win,domain, complete=.false.)
1147  call mpp_update_domains(wou,domain, complete=.true.)
1148 
1149  do j=js,je
1150  do i=is,ie+1
1151  if ( fx4(i,j) > 0. ) then
1152  fx4(i,j) = min(1., wou(i-1,j), win(i,j)) * fx4(i,j)
1153  else
1154  fx4(i,j) = min(1., win(i-1,j), wou(i,j)) * fx4(i,j)
1155  endif
1156  enddo
1157  enddo
1158  do j=js,je+1
1159  do i=is,ie
1160  if ( fy4(i,j) > 0. ) then
1161  fy4(i,j) = min(1., wou(i,j-1), win(i,j)) * fy4(i,j)
1162  else
1163  fy4(i,j) = min(1., win(i,j-1), wou(i,j)) * fy4(i,j)
1164  endif
1165  enddo
1166  enddo
1167 
1168 
1169  if ( zero_ocean ) then
1170 ! Limit diffusive flux over ocean cells:
1171  do j=js,je
1172  do i=is,ie+1
1173  fx4(i,j) = max(0., min(oro(i-1,j), oro(i,j))) * fx4(i,j)
1174  enddo
1175  enddo
1176  do j=js,je+1
1177  do i=is,ie
1178  fy4(i,j) = max(0., min(oro(i,j-1), oro(i,j))) * fy4(i,j)
1179  enddo
1180  enddo
1181  endif
1182 
1183 ! Update:
1184  do j=js,je
1185  do i=is,ie
1186  q(i,j) = qlow(i,j) + (fx4(i,j)-fx4(i+1,j)+fy4(i,j)-fy4(i,j+1))/area(i,j)
1187  enddo
1188  enddo
1189 
1190  enddo ! end n-loop
1191 
1192  end subroutine del4_cubed_sphere
1193 
1194 
1195  subroutine map_to_cubed_raw(igh, im, jt, lat1, lon1, zs, ft, grid, agrid, &
1196  q2, f2, h2, npx, npy, jstart, jend, stretch_fac, &
1197  nested, npx_global, bd, regional)
1199 ! Input
1200  type(fv_grid_bounds_type), intent(IN) :: bd
1201  integer, intent(in):: igh, im, jt
1202  integer, intent(in):: npx, npy, npx_global
1203  real, intent(in):: lat1(jt+1)
1204  real, intent(in):: lon1(im+1)
1205  real(kind=4), intent(in), dimension(-igh:im+igh,jt):: zs, ft
1206  real(kind=R_GRID), intent(in):: grid(bd%isd:bd%ied+1, bd%jsd:bd%jed+1,2)
1207  real(kind=R_GRID), intent(in):: agrid(bd%isd:bd%ied, bd%jsd:bd%jed, 2)
1208  integer, intent(in):: jstart, jend
1209  real(kind=R_GRID), intent(IN) :: stretch_fac
1210  logical, intent(IN) :: nested, regional
1211 ! Output
1212  real, intent(out):: q2(bd%isd:bd%ied,bd%jsd:bd%jed)
1213  real, intent(out):: f2(bd%isd:bd%ied,bd%jsd:bd%jed)
1214  real, intent(out):: h2(bd%isd:bd%ied,bd%jsd:bd%jed)
1215 ! Local
1216  real :: lon_g(-igh:im+igh)
1217  real lat_g(jt), cos_g(jt)
1218  real(kind=R_GRID) e2(2)
1219  real(kind=R_GRID) grid3(3, bd%isd:bd%ied+1, bd%jsd:bd%jed+1)
1220  real(kind=R_GRID), dimension(3):: p1, p2, p3, p4, pc, pp
1221  real(kind=R_GRID), dimension(3):: vp_12, vp_23, vp_34, vp_14
1222  integer i,j, np, k
1223  integer ii, jj, i1, i2, j1, j2, min_pts
1224  real th1, aa, asum, qsum, fsum, hsum, lon_w, lon_e, lat_s, lat_n, r2d
1225  real qsp, fsp, hsp
1226  real qnp, fnp, hnp
1227  real delg, th0, tmp1, prod1, prod2, prod3, prod4
1228  integer ig_lon, jp
1229  integer:: lat_crit
1230  real pi5, pi2
1231 
1232  integer :: is, ie, js, je
1233  integer :: isd, ied, jsd, jed
1234 
1235  is = bd%is
1236  ie = bd%ie
1237  js = bd%js
1238  je = bd%je
1239  isd = bd%isd
1240  ied = bd%ied
1241  jsd = bd%jsd
1242  jed = bd%jed
1243 
1244  pi2 = pi + pi
1245  pi5 = 0.5 * pi
1246  r2d = 180./pi
1247 
1248  do i=1,im
1249  lon_g(i) = 0.5*(lon1(i)+lon1(i+1))
1250  enddo
1251  do i=-igh,0
1252  lon_g(i) = lon_g(i+im)
1253  enddo
1254  do i=im+1,im+igh
1255  lon_g(i) = lon_g(i-im)
1256  enddo
1257 
1258  do j=1,jt
1259  lat_g(j) = 0.5*(lat1(j)+lat1(j+1))
1260  cos_g(j) = cos( lat_g(j) )
1261  enddo
1262 
1263  do j=jsd,jed+1
1264  do i=isd,ied+1
1265  call latlon2xyz(real(grid(i,j,1:2),kind=R_GRID), grid3(1,i,j))
1266  enddo
1267  enddo
1268 
1269  if(is_master()) write(*,*) 'surf_map: Search started ....'
1270 
1271 ! stretch_fac * pi5/(npx-1) / (pi/nlat)
1272  lat_crit = nint( stretch_fac*real(nlat)/real(npx_global-1) )
1273  lat_crit = min( jt, max( 4, lat_crit ) )
1274 
1275  if ( jstart==1 ) then
1276  write(*,*) mpp_pe(), 'lat_crit=', r2d*lat_g(lat_crit)
1277  elseif ( jend==nlat ) then
1278 ! write(*,*) mpp_pe(), 'lat_crit=', r2d*lat_g(jt-lat_crit+1)
1279  endif
1280 
1281 !----
1282 ! SP:
1283 !----
1284  iF ( jstart == 1 ) then
1285  asum = 0.
1286  qsum = 0.
1287  fsum = 0.
1288  hsum = 0.
1289  do j=1,lat_crit
1290  aa = cos_g(j)
1291  do i=1,im
1292  asum = asum + aa
1293  qsum = qsum + zs(i,j)*aa
1294  fsum = fsum + ft(i,j)*aa
1295  enddo
1296  enddo
1297  qsp = qsum / asum
1298  fsp = fsum / asum
1299  hsum = 0.
1300  np = 0
1301  do j=1,lat_crit
1302  do i=1,im
1303  np = np + 1
1304  hsum = hsum + (qsp-zs(i,j))**2
1305  enddo
1306  enddo
1307  hsp = hsum / real(np)
1308 ! write(*,*) 'SP GID, zs_sp, f_sp, hsp=', mpp_pe(), qsp, fsp, hsp
1309  endif
1310 !----
1311 ! NP:
1312 !----
1313  iF ( jend == nlat ) then
1314  asum = 0.
1315  qsum = 0.
1316  fsum = 0.
1317  hsum = 0.
1318  do jp=jend-lat_crit+1, jend
1319  j = jp - jstart + 1
1320  aa = cos_g(j)
1321  do i=1,im
1322  asum = asum + aa
1323  qsum = qsum + zs(i,j)*aa
1324  fsum = fsum + ft(i,j)*aa
1325  enddo
1326  enddo
1327  qnp = qsum / asum
1328  fnp = fsum / asum
1329  hsum = 0.
1330  np = 0
1331  do jp=jend-lat_crit+1, jend
1332  j = jp - jstart + 1
1333  do i=1,im
1334  np = np + 1
1335  hsum = hsum + (qnp-zs(i,j))**2
1336  enddo
1337  enddo
1338  hnp = hsum / real(np)
1339 ! write(*,*) 'NP GID, zs_np, f_np, hnp=', mpp_pe(), qnp, fnp, hnp
1340  endif
1341 
1342  min_pts = 999999
1343  do j=jsd,jed
1344  do i=isd,ied
1345 
1346  !Do not go into corners
1347  if (((i < is .and. j < js) .or. &
1348  (i < is .and. j > je) .or. &
1349  (i > ie .and. j < js) .or. &
1350  (i > ie .and. j > je)) .and. .not. (nested .or. regional)) then
1351  q2(i,j) = 1.e25
1352  f2(i,j) = 1.e25
1353  h2(i,j) = 1.e25
1354  goto 4444
1355  end if
1356 
1357  if ( agrid(i,j,2) < -pi5+stretch_fac*pi5/real(npx_global-1) ) then
1358 ! SP:
1359  q2(i,j) = qsp
1360  f2(i,j) = fsp
1361  h2(i,j) = hsp
1362  goto 4444
1363  elseif ( agrid(i,j,2) > pi5-stretch_fac*pi5/real(npx_global-1) ) then
1364 ! NP:
1365  q2(i,j) = qnp
1366  f2(i,j) = fnp
1367  h2(i,j) = hnp
1368  goto 4444
1369  endif
1370 
1371  lat_s = min( grid(i,j,2), grid(i+1,j,2), grid(i,j+1,2), grid(i+1,j+1,2), agrid(i,j,2) )
1372  lat_n = max( grid(i,j,2), grid(i+1,j,2), grid(i,j+1,2), grid(i+1,j+1,2), agrid(i,j,2) )
1373 ! Enlarge the search zone:
1374  delg = max( 0.2*(lat_n-lat_s), pi/real(npx_global-1), pi2/real(nlat) )
1375  lat_s = max(-pi5, lat_s - delg)
1376  lat_n = min( pi5, lat_n + delg)
1377 
1378  j1 = nint( (pi5+lat_s)/(pi/real(nlat)) ) - 1
1379  if ( lat_s*r2d < (-90.+90./real(npx_global-1)) ) j1 = 1
1380  j1 = max(jstart, j1)
1381 
1382  j2 = nint( (pi5+lat_n)/(pi/real(nlat)) ) + 1
1383  if ( lat_n*r2d > (90.-90./real(npx_global-1)) ) j2 = nlat
1384  j2 = min(jend, j2)
1385 
1386  j1 = j1 - jstart + 1
1387  j2 = j2 - jstart + 1
1388 
1389  lon_w = min( grid(i,j,1), grid(i+1,j,1), grid(i,j+1,1), grid(i+1,j+1,1) )
1390  lon_e = max( grid(i,j,1), grid(i+1,j,1), grid(i,j+1,1), grid(i+1,j+1,1) )
1391 
1392  if ( (lon_e-lon_w) > pi ) then
1393  i1 = -nint( (pi2-lon_e)/pi2 * real(im) ) - 1
1394  i2 = nint( lon_w/pi2 * real(im) ) + 1
1395  else
1396  i1 = nint( lon_w/pi2 * real(im) ) - 1
1397  i2 = nint( lon_e/pi2 * real(im) ) + 1
1398  endif
1399 
1400 ! Enlarge the search zone:
1401  ig_lon = max(1, (i2-i1)/8)
1402  i1 = max( -igh, i1 - ig_lon)
1403  i2 = min(im+igh, i2 + ig_lon)
1404 
1405  np = 0
1406  qsum = 0.
1407  fsum = 0.
1408  hsum = 0.
1409  asum = 0.
1410 !
1411 ! 4----------3
1412 ! / /
1413 ! / pp /
1414 ! / /
1415 ! 1----------2
1416 !
1417  do k=1,3
1418  p1(k) = grid3(k,i, j)
1419  p2(k) = grid3(k,i+1,j)
1420  p3(k) = grid3(k,i+1,j+1)
1421  p4(k) = grid3(k,i,j+1)
1422  pc(k) = p1(k) + p2(k) + p3(k) + p4(k)
1423  enddo
1424  call normalize_vect( pc )
1425 
1426  th0 = min( v_prod(p1,p3), v_prod(p2, p4) )
1427  th1 = min( cos_grid, cos(0.25*acos(max(v_prod(p1,p3), v_prod(p2, p4)))))
1428 
1429  call vect_cross(vp_12, p1, p2)
1430  call vect_cross(vp_23, p2, p3)
1431  call vect_cross(vp_34, p3, p4)
1432  call vect_cross(vp_14, p1, p4)
1433 
1434  prod1 = v_prod(p3, vp_12)
1435  prod2 = v_prod(p1, vp_23)
1436  prod3 = v_prod(p1, vp_34)
1437  prod4 = v_prod(p2, vp_14)
1438 
1439  do jj=j1,j2
1440  aa = cos_g(jj)
1441  e2(2) = lat_g(jj)
1442  do ii=i1,i2
1443  e2(1) = lon_g(ii)
1444  call latlon2xyz(e2, pp)
1445 ! Check two extrems:
1446  tmp1 = v_prod(pp, pc)
1447 ! Points that are close to center:
1448  if ( tmp1 > th1 ) goto 1111 ! inside
1449 ! check to exclude points too far away:
1450  if ( tmp1 < th0 ) goto 2222 ! outside
1451 ! Check if inside the polygon
1452  if ( v_prod(pp, vp_12)*prod1 < 0. ) goto 2222
1453  if ( v_prod(pp, vp_23)*prod2 < 0. ) goto 2222
1454  if ( v_prod(pp, vp_34)*prod3 < 0. ) goto 2222
1455  if ( v_prod(pp, vp_14)*prod4 < 0. ) goto 2222
1456 1111 np = np + 1
1457  qsum = qsum + zs(ii,jj)*aa
1458  fsum = fsum + ft(ii,jj)*aa
1459  hsum = hsum + zs(ii,jj)**2
1460  asum = asum + aa
1461 2222 continue
1462  enddo
1463  enddo
1464 
1465  if ( np > 0 ) then
1466  q2(i,j) = qsum / asum
1467  f2(i,j) = fsum / asum
1468  h2(i,j) = hsum / real(np) - q2(i,j)**2
1469  min_pts = min(min_pts, np)
1470  else
1471  write(*,*) 'min and max lat_g is ', r2d*minval(lat_g), r2d*maxval(lat_g), mpp_pe()
1472  write(*,*) 'Surf_map failed for GID=', mpp_pe(), i,j, '(lon,lat)=', r2d*agrid(i,j,1),r2d*agrid(i,j,2)
1473  write(*,*) '[jstart, jend]', jstart, jend
1474  call mpp_error(fatal,'Surf_map failed')
1475  endif
1476 4444 continue
1477  enddo
1478  enddo
1479 
1480  if(is_master()) write(*,*) 'surf_map: minimum pts per cell (master PE)=', min_pts
1481  if ( min_pts <3 ) then
1482  if(is_master()) write(*,*) 'Warning: too few points used in creating the cell mean terrain !!!'
1483  endif
1484 
1485  end subroutine map_to_cubed_raw
1486 
1487 
1488 #ifdef JUNK
1489  logical function inside_p4(p1, p2, p3, p4, pp, th0)
1491  real, intent(in):: p1(3), p2(3), p3(3), p4(3)
1492  real, intent(in):: pp(3)
1493  real, intent(in):: th0
1494 ! A * B = |A| |B| cos(angle)
1495 ! Local:
1496  real(kind=R_GRID) v1(3), v2(3), vp(3)
1497  real(kind=R_GRID) tmp1
1498  integer k
1499 
1500  inside_p4 = .false.
1501 
1502 ! 1st check: to exclude points too far away:
1503  do k=1,3
1504  vp(k) = p1(k) + p2(k) + p3(k) + p4(k)
1505  enddo
1506  call normalize_vect( vp )
1507 
1508  tmp1 = v_prod(pp, vp)
1509  if ( tmp1 < th0 ) then ! outside
1510  return
1511  endif
1512 
1513 ! 1st segment: 1---2
1514  call vect_cross(vp, p1, p2)
1515  if ( v_prod(pp, vp)*v_prod(p3, vp) < 0. ) return
1516 ! 2nd segment: 2---3
1517  call vect_cross(vp, p2, p3)
1518  if ( v_prod(pp, vp)*v_prod(p1, vp) < 0. ) return
1519 ! 3rd segment: 3---4
1520  call vect_cross(vp, p3, p4)
1521  if ( v_prod(pp, vp)*v_prod(p1, vp) < 0. ) return
1522 ! 4th segment: 1---4
1523  call vect_cross(vp, p1, p4)
1524  if ( v_prod(pp, vp)*v_prod(p2, vp) < 0. ) return
1525 
1526  inside_p4 = .true.
1527 
1528  end function inside_p4
1529 #endif
1530 
1533  subroutine handle_err(status)
1534 #include <netcdf.inc>
1535  integer status
1536 
1537  if (status .ne. nf_noerr) then
1538  print *, nf_strerror(status)
1539  stop 'Stopped'
1540  endif
1541 
1542  end subroutine handle_err
1543 
1544  subroutine remove_ice_sheets (lon, lat, lfrac, bd )
1545 !---------------------------------
1546 ! Bruce Wyman's fix for Antarctic
1547 !---------------------------------
1548  type(fv_grid_bounds_type), intent(IN) :: bd
1549  real(kind=R_GRID), intent(in) :: lon(bd%isd:bd%ied,bd%jsd:bd%jed), lat(bd%isd:bd%ied,bd%jsd:bd%jed)
1550  real, intent(inout) :: lfrac(bd%isd:bd%ied,bd%jsd:bd%jed)
1551 
1552  integer :: is, ie, js, je
1553  integer :: isd, ied, jsd, jed
1554 
1555 ! lon = longitude in radians
1556 ! lat = latitude in radians
1557 ! lfrac = land-sea mask (land=1, sea=0)
1558 
1559  integer :: i, j
1560  real :: dtr, phs, phn
1561 
1562  is = bd%is
1563  ie = bd%ie
1564  js = bd%js
1565  je = bd%je
1566  isd = bd%isd
1567  ied = bd%ied
1568  jsd = bd%jsd
1569  jed = bd%jed
1570 
1571  dtr = acos(0.)/90.
1572  phs = -83.9999*dtr
1573 ! phn = -78.9999*dtr
1574  phn = -76.4*dtr
1575 
1576  do j = jsd, jed
1577  do i = isd, ied
1578  if ( lat(i,j) < phn ) then
1579  ! replace all below this latitude
1580  if ( lat(i,j) < phs ) then
1581  lfrac(i,j) = 1.0
1582  cycle
1583  endif
1584  ! replace between 270 and 360 deg
1585  if ( sin(lon(i,j)) < 0. .and. cos(lon(i,j)) > 0.) then
1586  lfrac(i,j) = 1.0
1587  cycle
1588  endif
1589  endif
1590  enddo
1591  enddo
1592  end subroutine remove_ice_sheets
1593 
1596 subroutine read_namelist
1597  integer :: unit, ierr, io
1598 ! real :: dtr, ght
1599 
1600 ! read namelist
1601 
1602  if (namelist_read) return
1603 #ifdef INTERNAL_FILE_NML
1604  read (input_nml_file, nml=surf_map_nml, iostat=io)
1605  ierr = check_nml_error(io,'surf_map_nml')
1606 #else
1607  unit = open_namelist_file( )
1608  ierr=1
1609  do while (ierr /= 0)
1610  read (unit, nml=surf_map_nml, iostat=io, end=10)
1611  ierr = check_nml_error(io,'surf_map_nml')
1612  enddo
1613  10 call close_file (unit)
1614 #endif
1615 
1616 ! write version and namelist to log file
1617 
1618  if (mpp_pe() == mpp_root_pe()) then
1619  unit = stdlog()
1620  write (unit, nml=surf_map_nml)
1621  endif
1622 
1623  namelist_read = .true.
1624 
1625 end subroutine read_namelist
1626 
1628 subroutine zonal_mean(im, p, zmean)
1629  integer, intent(in):: im
1630  real(kind=4), intent(inout):: p(im)
1631  real, intent(out):: zmean
1632  integer i
1633 
1634  zmean = 0.
1635  do i=1,im
1636  zmean = zmean + p(i)
1637  enddo
1638  zmean = zmean / real(im)
1639  do i=1,im
1640  p(i) = zmean
1641  enddo
1642 end subroutine zonal_mean
1643 
1644 
1645  end module fv_surf_map_mod
subroutine, public del2_cubed_sphere(npx, npy, q, area, dx, dy, dxc, dyc, sin_sg, nmax, cd, zero_ocean, oro, nested, domain, bd, regional)
subroutine zonal_mean(im, p, zmean)
The sugroutine &#39;zonal_mean&#39; replaces &#39;p&#39; with its zonal mean.
real max_slope
max allowable terrain slope: 1 –> 45 deg 0.15 for C768 or lower; 0.25 C1536; 0.3 for C3072 ...
character(len=3) grid_string
subroutine remove_ice_sheets(lon, lat, lfrac, bd)
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
subroutine handle_err(status)
The subroutine &#39;handle_err&#39; returns an error when it cannot find or correctly read in an external fil...
subroutine, public normalize_vect(e)
The subroutine &#39;normalize_vect&#39; makes &#39;e&#39; a unit vector.
real, dimension(:,:), allocatable, public zs_g
subroutine two_delta_filter(npx, npy, q, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, cd, zero_ocean, check_slope, filter_type, oro, nested, domain, bd, ntmax, regional)
subroutine, public fv3_zs_filter(bd, isd, ied, jsd, jed, npx, npy, npx_global, stretch_fac, nested, domain, area, dxa, dya, dx, dy, dxc, dyc, grid, agrid, sin_sg, phis, oro, regional)
subroutine, public del4_cubed_sphere(npx, npy, q, area, dx, dy, dxc, dyc, sin_sg, nmax, zero_ocean, oro, nested, domain, bd, regional)
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
The function &#39;g_sum&#39; is the fast version of &#39;globalsum&#39;.
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
subroutine read_namelist
The subroutine &#39;read_namelis&#39; reads the namelist file, writes the namelist to log file...
subroutine, public global_mx(q, n_g, qmin, qmax, bd)
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
real function, public great_circle_dist(q1, q2, radius)
character(len=6) surf_format
real(kind=r_grid) function, public v_prod(v1, v2)
integer n_del2_strong
real, dimension(:,:), allocatable, public sgh_g
subroutine, public vect_cross(e, p1, p2)
The subroutine &#39;vect_cross performs cross products of 3D vectors: e = P1 X P2.
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
integer, parameter, public ng
Definition: fv_mp_mod.F90:2716
real cd4
Dimensionless coeff for del-4 diffusion (with FCT)
subroutine, public surfdrv(npx, npy, grid, agrid, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, phis, stretch_fac, nested, npx_global, domain, grid_number, bd, regional)
logical function inside_p4(p1, p2, p3, p4, pp, th0)
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
logical namelist_read
real, dimension(:,:), allocatable, public oro_g
subroutine map_to_cubed_raw(igh, im, jt, lat1, lon1, zs, ft, grid, agrid, q2, f2, h2, npx, npy, jstart, jend, stretch_fac, nested, npx_global, bd, regional)
real(kind=r_grid) da_min
character(len=128) surf_file
real cd2
Dimensionless coeff for del-2 diffusion (-1 gives resolution-determined value)
real peak_fac
overshoot factor for the mountain peak
subroutine, public latlon2xyz(p, e, id)
The subroutine &#39;latlon2xyz&#39; maps (lon, lat) to (x,y,z)