Line data Source code
1 : ! Separate dynamics and physics grids
2 :
3 : module dp_mapping
4 : use dimensions_mod, only: np, npsq, fv_nphys
5 : use shr_kind_mod, only: r8=>shr_kind_r8, shr_kind_cl
6 : use coordinate_systems_mod, only: spherical_polar_t
7 : use shr_const_mod, only: pi => shr_const_pi
8 : use fvm_control_volume_mod, only: fvm_struct
9 :
10 : implicit none
11 : private
12 : save
13 :
14 : public :: dp_init
15 : public :: dp_reorder
16 : public :: dp_write
17 : public :: dp_allocate
18 : public :: dp_deallocate
19 :
20 : ! Total number of physics points per spectral element
21 : ! no physgrid: nphys_pts = npsq (physics on GLL grid)
22 : ! physgrid: nphys_pts = nphys2 (physics on CSLAM grid)
23 : ! Value is set when se_fv_nphys namelist variable is read
24 : integer, public :: nphys_pts = npsq
25 :
26 : ! NOTE: dp_gid() is in space filling curve rank order
27 : ! all other global arrays are in block id (global id) order
28 : !
29 : ! dp_gid() is used to re-order data collected on root via mpi_gatherv
30 : ! into block id ordering
31 : !
32 : ! j=dp_gid(i) i = element space filling curve rank
33 : ! j = element global id = block id = history file ordering
34 : !
35 : integer, allocatable,dimension(:) :: dp_gid ! NE=240, integer*4 = 1.3MB
36 : integer, public,allocatable,dimension(:) :: dp_owner
37 :
38 : real (r8),public,allocatable,dimension(:,:,:) :: weights_all_fvm2phys
39 : integer ,public,allocatable,dimension(:,:,:) :: weights_eul_index_all_fvm2phys,weights_lgr_index_all_fvm2phys
40 : real (r8),public,allocatable,dimension(:,:,:) :: weights_all_phys2fvm
41 : integer ,public,allocatable,dimension(:,:,:) :: weights_eul_index_all_phys2fvm,weights_lgr_index_all_phys2fvm
42 : integer ,public,allocatable,dimension(:) :: jall_fvm2phys,jall_phys2fvm
43 : integer ,public :: num_weights_fvm2phys,num_weights_phys2fvm
44 :
45 :
46 : contains
47 1536 : subroutine dp_init(elem,fvm)
48 : use dimensions_mod, only: nelemd, nc, irecons_tracer
49 : use element_mod, only: element_t
50 : use spmd_utils, only: masterproc
51 : use cam_logfile, only: iulog
52 :
53 : type(element_t) , dimension(nelemd), intent(in) :: elem
54 : type (fvm_struct), dimension(nelemd), intent(in) :: fvm
55 :
56 1536 : num_weights_phys2fvm = 0
57 1536 : num_weights_fvm2phys = 0
58 1536 : if (fv_nphys>0) then
59 1536 : num_weights_phys2fvm = (nc+fv_nphys)**2
60 1536 : num_weights_fvm2phys = (nc+fv_nphys)**2
61 :
62 6144 : allocate(weights_all_fvm2phys(num_weights_fvm2phys,irecons_tracer,nelemd))
63 6144 : allocate(weights_eul_index_all_fvm2phys(num_weights_fvm2phys,2,nelemd))
64 4608 : allocate(weights_lgr_index_all_fvm2phys(num_weights_fvm2phys,2,nelemd))
65 :
66 4608 : allocate(weights_all_phys2fvm(num_weights_phys2fvm,irecons_tracer,nelemd))
67 4608 : allocate(weights_eul_index_all_phys2fvm(num_weights_phys2fvm,2,nelemd))
68 4608 : allocate(weights_lgr_index_all_phys2fvm(num_weights_phys2fvm,2,nelemd))
69 4608 : allocate(jall_fvm2phys(nelemd))
70 3072 : allocate(jall_phys2fvm(nelemd))
71 :
72 : call fvm2phys_init(elem,fvm,nc,fv_nphys,irecons_tracer,&
73 : weights_all_fvm2phys,weights_eul_index_all_fvm2phys,weights_lgr_index_all_fvm2phys,&
74 : weights_all_phys2fvm,weights_eul_index_all_phys2fvm,weights_lgr_index_all_phys2fvm,&
75 1536 : jall_fvm2phys,jall_phys2fvm)
76 :
77 1536 : if (masterproc) then
78 2 : write(iulog, *) 'dp_init: Initialized phys2fvm/fvm2phys mapping vars'
79 : end if
80 :
81 : end if
82 1536 : end subroutine dp_init
83 :
84 0 : subroutine dp_reorder(before, after)
85 : use cam_abortutils, only: endrun
86 : use dimensions_mod, only: nelem
87 : use cam_logfile, only: iulog
88 : use spmd_utils, only: masterproc
89 : use shr_sys_mod, only: shr_sys_flush
90 :
91 : real(r8), dimension(fv_nphys*fv_nphys,*), intent(in) :: before
92 : real(r8), dimension(fv_nphys*fv_nphys,*), intent(out) :: after
93 : integer :: ie
94 :
95 : ! begin
96 0 : do ie = 1,nelem
97 0 : if (dp_gid(ie) < 0) then
98 0 : if (masterproc) then
99 0 : write(iulog,*) 'ie =',ie,', dp_gid(ie) =',dp_gid(ie)
100 0 : call shr_sys_flush(iulog)
101 : end if
102 0 : call endrun('Bad element remap in dp_reorder')
103 : end if
104 0 : after(:,dp_gid(ie)) = before(:,ie)
105 : end do
106 0 : end subroutine dp_reorder
107 :
108 : !!!
109 :
110 0 : subroutine dp_allocate(elem)
111 : use dimensions_mod, only: nelem, nelemd
112 : use element_mod, only: element_t
113 : use spmd_utils, only: masterproc, masterprocid, npes
114 : use spmd_utils, only: mpicom, mpi_integer
115 :
116 : implicit none
117 : type(element_t),dimension(nelemd),intent(in) :: elem
118 :
119 : integer :: i,j,ierror
120 0 : integer,dimension(nelemd) :: lgid
121 0 : integer,dimension(:),allocatable :: displs,recvcount
122 :
123 : ! begin
124 :
125 0 : allocate(displs(npes))
126 0 : allocate(dp_gid(nelem))
127 0 : allocate(recvcount(npes))
128 : call mpi_gather(nelemd, 1, mpi_integer, recvcount, 1, mpi_integer, &
129 0 : masterprocid, mpicom, ierror)
130 0 : lgid(:) = elem(:)%globalid
131 0 : if (masterproc) then
132 0 : displs(1) = 0
133 0 : do i = 2,npes
134 0 : displs(i) = displs(i-1)+recvcount(i-1)
135 : end do
136 : end if
137 : call mpi_gatherv(lgid, nelemd, mpi_integer, dp_gid, recvcount, displs, &
138 0 : mpi_integer, masterprocid, mpicom, ierror)
139 0 : if (masterproc) then
140 0 : allocate(dp_owner(nelem))
141 0 : dp_owner(:) = -1
142 0 : do i = 1,npes
143 0 : do j = displs(i)+1,displs(i)+recvcount(i)
144 0 : dp_owner(dp_gid(j)) = i-1
145 : end do
146 : end do
147 : end if
148 0 : deallocate(displs)
149 0 : deallocate(recvcount)
150 : ! minimize global memory use
151 0 : call mpi_barrier(mpicom,ierror)
152 0 : if (.not.masterproc) then
153 0 : allocate(dp_owner(nelem))
154 : end if
155 0 : call mpi_bcast(dp_gid,nelem,mpi_integer,masterprocid,mpicom,ierror)
156 0 : call mpi_bcast(dp_owner,nelem,mpi_integer,masterprocid,mpicom,ierror)
157 0 : end subroutine dp_allocate
158 :
159 : !!!
160 :
161 0 : subroutine dp_deallocate()
162 0 : deallocate(dp_gid)
163 0 : deallocate(dp_owner)
164 0 : end subroutine dp_deallocate
165 :
166 : !!!
167 :
168 0 : subroutine dp_write(elem, fvm, grid_format, filename_in)
169 : use cam_abortutils, only: endrun
170 : use dimensions_mod, only: nelem, nelemd
171 : use element_mod, only: element_t
172 : use netcdf, only: nf90_create, nf90_close, nf90_enddef
173 : use netcdf, only: nf90_def_dim, nf90_def_var, nf90_put_var
174 : use netcdf, only: nf90_double, nf90_int, nf90_put_att
175 : use netcdf, only: nf90_noerr, nf90_strerror, nf90_clobber
176 : use spmd_utils, only: masterproc, masterprocid, mpicom, npes
177 : use spmd_utils, only: mpi_integer, mpi_real8
178 : use cam_logfile, only: iulog
179 : use shr_sys_mod, only: shr_sys_flush
180 : use dimensions_mod, only: ne
181 : use coordinate_systems_mod, only: cart2spherical
182 :
183 : ! Inputs
184 : type(element_t), intent(in) :: elem(:)
185 : type (fvm_struct), intent(in) :: fvm(:)
186 : character(len=*), intent(in) :: grid_format
187 : character(len=*), intent(in) :: filename_in
188 :
189 : real(r8), parameter :: rad2deg = 180._r8/pi
190 :
191 : ! Local variables
192 : integer :: i, ie, ierror, j, status, ivtx
193 : integer :: grid_corners_id, grid_rank_id, grid_size_id
194 : integer :: ncid
195 : integer :: grid_dims_id, grid_area_id, grid_center_lat_id
196 : integer :: grid_center_lon_id, grid_corner_lat_id
197 : integer :: grid_corner_lon_id, grid_imask_id
198 : integer :: gridsize
199 : integer :: IOrootID
200 : logical :: IOroot
201 : character(len=SHR_KIND_CL) :: errormsg
202 : character(len=SHR_KIND_CL) :: filename
203 0 : integer,allocatable,dimension(:) :: displs,recvcount
204 :
205 0 : real(r8), dimension(fv_nphys, fv_nphys, nelemd, 4, 2) :: corners
206 0 : real(r8), dimension(fv_nphys, fv_nphys, nelemd) :: lwork
207 0 : real(r8), allocatable, dimension(:) :: recvbuf
208 0 : real(r8), allocatable, dimension(:,:) :: gwork
209 : real(r8) :: x, y
210 : type (spherical_polar_t) :: sphere
211 :
212 : ! begin
213 :
214 : !! Check to see if we are doing grid output
215 0 : if (trim(grid_format) == "no") then
216 0 : if (masterproc) then
217 0 : write(iulog, *) 'dp_write: Not writing phys_grid file.'
218 : end if
219 0 : return
220 0 : else if (trim(grid_format) /= 'SCRIP') then
221 0 : if (masterproc) then
222 0 : write(errormsg, *) 'dp_write: ERROR, bad value for se_write_grid, ',&
223 0 : trim(grid_format)
224 0 : call endrun(errormsg)
225 : end if
226 : end if
227 :
228 : ! Create the NetCDF file
229 0 : if (len_trim(filename_in) == 0) then
230 0 : write(filename, '(3(a,i0),3a)') "ne", ne, "np", np, ".pg", fv_nphys, &
231 0 : "_", trim(grid_format), ".nc"
232 : else
233 0 : filename = trim(filename_in)
234 : end if
235 0 : status = nf90_create(trim(filename), nf90_clobber, ncid)
236 0 : if (status /= nf90_noerr) then
237 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
238 : end if
239 :
240 : ! PIO_put_var puts from its root node, find that (so we do our work there)
241 0 : IOrootID = masterprocid
242 0 : IOroot = masterproc
243 :
244 : ! Allocate workspace and calculate PE displacement information
245 0 : if (IOroot) then
246 0 : allocate(displs(npes))
247 0 : allocate(recvcount(npes))
248 : else
249 0 : allocate(displs(0))
250 0 : allocate(recvcount(0))
251 : end if
252 0 : gridsize = nelem * fv_nphys*fv_nphys
253 0 : if(masterproc) then
254 0 : write(iulog, *) 'Writing physics SCRIP grid file: ', trim(filename)
255 0 : write(iulog, '(a,i7,a,i3)') 'nelem = ', nelem, ', fv_nphys = ', fv_nphys
256 0 : call shr_sys_flush(iulog)
257 : end if
258 : call mpi_gather(nelemd*fv_nphys*fv_nphys, 1, mpi_integer, recvcount, 1, &
259 0 : mpi_integer, IOrootID, mpicom, ierror)
260 :
261 0 : if (IOroot) then
262 0 : displs(1) = 0
263 0 : do i = 2, npes
264 0 : displs(i) = displs(i-1)+recvcount(i-1)
265 : end do
266 0 : allocate(recvbuf(gridsize))
267 : else
268 0 : allocate(recvbuf(0))
269 : end if
270 0 : allocate(gwork(4, gridsize))
271 :
272 0 : if (IOroot) then
273 : ! Define the horizontal grid dimensions for SCRIP output
274 0 : status = nf90_def_dim(ncid, "grid_corners", 4, grid_corners_id)
275 0 : if (status /= nf90_noerr) then
276 0 : write(iulog, *) 'dp_write: Error defining dimension, grid_corners'
277 0 : call shr_sys_flush(iulog)
278 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
279 : end if
280 0 : status = nf90_def_dim(ncid, "grid_rank", 1, grid_rank_id)
281 0 : if (status /= nf90_noerr) then
282 0 : write(iulog, *) 'dp_write: Error defining dimension, grid_rank'
283 0 : call shr_sys_flush(iulog)
284 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
285 : end if
286 0 : status = nf90_def_dim(ncid, "grid_size", gridsize, grid_size_id)
287 0 : if (status /= nf90_noerr) then
288 0 : write(iulog, *) 'dp_write: Error defining dimension, grid_size'
289 0 : call shr_sys_flush(iulog)
290 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
291 : end if
292 :
293 : ! Define the coordinate variables
294 : status = nf90_def_var(ncid, "grid_dims", nf90_int, (/grid_rank_id/), &
295 0 : grid_dims_id)
296 0 : if (status /= nf90_noerr) then
297 0 : write(iulog, *) 'dp_write: Error defining variable grid_dims'
298 0 : call shr_sys_flush(iulog)
299 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
300 : end if
301 :
302 : status = nf90_def_var(ncid, "grid_area", nf90_double, &
303 0 : (/grid_size_id/), grid_area_id)
304 0 : if (status /= nf90_noerr) then
305 0 : write(iulog, *) 'dp_write: Error defining variable grid_area'
306 0 : call shr_sys_flush(iulog)
307 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
308 : end if
309 :
310 0 : status = nf90_put_att(ncid, grid_area_id, "units", "radians^2")
311 0 : if (status /= nf90_noerr) then
312 0 : write(iulog, *) 'dp_write: Error defining attributes for grid_area'
313 0 : call shr_sys_flush(iulog)
314 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
315 : end if
316 :
317 0 : status = nf90_put_att(ncid, grid_area_id, "long_name", "area weights")
318 0 : if (status /= nf90_noerr) then
319 0 : write(iulog, *) 'dp_write: Error defining attributes for grid_area'
320 0 : call shr_sys_flush(iulog)
321 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
322 : end if
323 :
324 : status = nf90_def_var(ncid, "grid_center_lat", nf90_double, &
325 0 : (/grid_size_id/), grid_center_lat_id)
326 0 : if (status /= nf90_noerr) then
327 0 : write(iulog, *) 'dp_write: Error defining variable grid_center_lat'
328 0 : call shr_sys_flush(iulog)
329 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
330 : end if
331 :
332 0 : status = nf90_put_att(ncid, grid_center_lat_id, "units", "degrees")
333 0 : if (status /= nf90_noerr) then
334 0 : write(iulog, *) 'dp_write: Error defining attributes for grid_center_lat'
335 0 : call shr_sys_flush(iulog)
336 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
337 : end if
338 :
339 : status = nf90_def_var(ncid, "grid_center_lon", nf90_double, &
340 0 : (/grid_size_id/), grid_center_lon_id)
341 0 : if (status /= nf90_noerr) then
342 0 : write(iulog, *) 'dp_write: Error defining variable grid_center_lon'
343 0 : call shr_sys_flush(iulog)
344 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
345 : end if
346 :
347 0 : status = nf90_put_att(ncid, grid_center_lon_id, "units", "degrees")
348 0 : if (status /= nf90_noerr) then
349 0 : write(iulog, *) 'dp_write: Error defining attributes for grid_center_lon'
350 0 : call shr_sys_flush(iulog)
351 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
352 : end if
353 :
354 : status = nf90_def_var(ncid, "grid_corner_lat", nf90_double, &
355 0 : (/grid_corners_id, grid_size_id/), grid_corner_lat_id)
356 0 : if (status /= nf90_noerr) then
357 0 : write(iulog, *) 'dp_write: Error defining grid_corner_lat'
358 0 : call shr_sys_flush(iulog)
359 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
360 : end if
361 :
362 0 : status = nf90_put_att(ncid, grid_corner_lat_id, "units", "degrees")
363 0 : if (status /= nf90_noerr) then
364 0 : write(iulog, *) 'dp_write: Error defining attributes for grid_corner_lat'
365 0 : call shr_sys_flush(iulog)
366 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
367 : end if
368 :
369 : status = nf90_def_var(ncid, "grid_corner_lon", nf90_double, &
370 0 : (/grid_corners_id, grid_size_id/), grid_corner_lon_id)
371 0 : if (status /= nf90_noerr) then
372 0 : write(iulog, *) 'dp_write: Error defining variable grid_corner_lon'
373 0 : call shr_sys_flush(iulog)
374 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
375 : end if
376 :
377 0 : status = nf90_put_att(ncid, grid_corner_lon_id, "units", "degrees")
378 0 : if (status /= nf90_noerr) then
379 0 : write(iulog, *) 'dp_write: Error defining attributes for grid_corner_lon'
380 0 : call shr_sys_flush(iulog)
381 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
382 : end if
383 :
384 : status = nf90_def_var(ncid, "grid_imask", nf90_double, &
385 0 : (/grid_size_id/), grid_imask_id)
386 0 : if (status /= nf90_noerr) then
387 0 : write(iulog, *) 'dp_write: Error defining variable grid_imask'
388 0 : call shr_sys_flush(iulog)
389 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
390 : end if
391 :
392 : ! End of NetCDF definitions
393 0 : status = nf90_enddef(ncid)
394 0 : if (status /= nf90_noerr) then
395 0 : write(iulog, *) 'dp_write: Error calling enddef'
396 0 : call shr_sys_flush(iulog)
397 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
398 : end if
399 : end if ! IOroot
400 :
401 : if (IOroot) then
402 0 : status = nf90_put_var(ncid, grid_dims_id, (/ gridsize /))
403 0 : if (status /= nf90_noerr) then
404 0 : write(iulog, *) 'dp_write: Error writing variable grid_dims'
405 0 : call shr_sys_flush(iulog)
406 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
407 : end if
408 : end if
409 :
410 0 : do ie=1,nelemd
411 0 : lwork(:,:,ie) = fvm(ie)%area_sphere_physgrid(:,:)
412 : end do
413 : call mpi_gatherv(lwork, size(lwork), mpi_real8, recvbuf, recvcount, &
414 0 : displs, mpi_real8, IOrootID, mpicom, ierror)
415 0 : call dp_allocate(elem)
416 0 : if (IOroot) then
417 0 : call dp_reorder(recvbuf, gwork(1,:))
418 0 : status = nf90_put_var(ncid, grid_area_id, gwork(1,:))
419 0 : if (status /= nf90_noerr) then
420 0 : write(iulog, *) 'dp_write: Error writing variable grid_area'
421 0 : call shr_sys_flush(iulog)
422 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
423 : end if
424 : end if
425 0 : do ie=1,nelemd
426 0 : lwork(:,:,ie) = rad2deg*fvm(ie)%center_cart_physgrid(:,:)%lat
427 : end do
428 : call mpi_gatherv(lwork, size(lwork), mpi_real8, recvbuf, recvcount, &
429 0 : displs, mpi_real8, IOrootID, mpicom, ierror)
430 0 : if (IOroot) then
431 0 : call dp_reorder(recvbuf, gwork(1,:))
432 0 : status = nf90_put_var(ncid, grid_center_lat_id, gwork(1,:))
433 0 : if (status /= nf90_noerr) then
434 0 : write(iulog, *) 'dp_write: Error writing variable grid_center_lat'
435 0 : call shr_sys_flush(iulog)
436 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
437 : end if
438 : end if
439 :
440 0 : do ie=1,nelemd
441 0 : lwork(:,:,ie) = rad2deg*fvm(ie)%center_cart_physgrid(:,:)%lon
442 : end do
443 : call mpi_gatherv(lwork, size(lwork), mpi_real8, recvbuf, recvcount, &
444 0 : displs, mpi_real8, IOrootID, mpicom, ierror)
445 0 : if (IOroot) then
446 0 : call dp_reorder(recvbuf, gwork(1,:))
447 0 : status = nf90_put_var(ncid, grid_center_lon_id, gwork(1,:))
448 0 : if (status /= nf90_noerr) then
449 0 : write(iulog, *) 'dp_write: Error writing variable grid_center_lon'
450 0 : call shr_sys_flush(iulog)
451 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
452 : end if
453 : end if
454 : ! compute physgrid grid corners
455 0 : do ie=1,nelemd
456 0 : do j=1,fv_nphys
457 0 : do i=1,fv_nphys
458 0 : do ivtx=1,4
459 0 : x = fvm(ie)%vtx_cart_physgrid(ivtx,1,i,j)
460 0 : y = fvm(ie)%vtx_cart_physgrid(ivtx,2,i,j)
461 0 : sphere = cart2spherical(x,y,elem(ie)%FaceNum)
462 0 : corners(i,j,ie,ivtx,1) = rad2deg * sphere%lat
463 0 : corners(i,j,ie,ivtx,2) = rad2deg * sphere%lon
464 : end do
465 : end do
466 : end do
467 : end do
468 : ! Collect all information for the grid corner latitude (counter-clockwise)
469 0 : do ivtx=1,4
470 : call mpi_gatherv(corners(:,:,:,ivtx,1), size(corners(:,:,:,ivtx,1)), mpi_real8, recvbuf, recvcount, &
471 0 : displs, mpi_real8, IOrootID, mpicom, ierror)
472 0 : if (IOroot) then
473 0 : call dp_reorder(recvbuf, gwork(ivtx,:))
474 : end if
475 : end do
476 0 : if (IOroot) then
477 0 : status = nf90_put_var(ncid, grid_corner_lat_id, gwork)
478 0 : if (status /= nf90_noerr) then
479 0 : write(iulog, *) 'dp_write: Error writing variable grid_corner_lat'
480 0 : call shr_sys_flush(iulog)
481 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
482 : end if
483 : end if
484 : ! Collect all information for the grid corner longitudes (counter-clockwise)
485 0 : do ivtx=1,4
486 : call mpi_gatherv(corners(:,:,:,ivtx,2), size(corners(:,:,:,ivtx,2)), mpi_real8, recvbuf, recvcount, &
487 0 : displs, mpi_real8, IOrootID, mpicom, ierror)
488 0 : if (IOroot) then
489 0 : call dp_reorder(recvbuf, gwork(ivtx,:))
490 : end if
491 : end do
492 0 : call dp_deallocate()
493 0 : if (IOroot) then
494 0 : status = nf90_put_var(ncid, grid_corner_lon_id, gwork)
495 0 : if (status /= nf90_noerr) then
496 0 : write(iulog, *) 'dp_write: Error writing variable grid_corner_lon'
497 0 : call shr_sys_flush(iulog)
498 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
499 : end if
500 : end if
501 :
502 : if (IOroot) then
503 0 : gwork(1,:) = 1._r8
504 0 : status = nf90_put_var(ncid, grid_imask_id, gwork(1,:))
505 0 : if (status /= nf90_noerr) then
506 0 : write(iulog, *) 'dp_write: Error writing variable grid_imask'
507 0 : call shr_sys_flush(iulog)
508 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
509 : end if
510 : end if
511 :
512 : ! call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
513 : ! Close the file
514 0 : call mpi_barrier(mpicom, ierror)
515 0 : if (IOroot) then
516 0 : status = nf90_close(ncid)
517 0 : if (status /= nf90_noerr) then
518 0 : call endrun("dp_write: "//trim(nf90_strerror(status)))
519 : end if
520 : end if
521 :
522 0 : call mpi_barrier(mpicom, ierror)
523 0 : if(masterproc) then
524 0 : write(iulog, *) 'Finished writing physics grid file: ', trim(filename)
525 0 : call shr_sys_flush(iulog)
526 : end if
527 :
528 0 : end subroutine dp_write
529 : !!!
530 :
531 1536 : subroutine fvm2phys_init(elem,fvm,fvm_nc,phys_nc,irecons,&
532 1536 : weights_all_fvm2phys,weights_eul_index_all_fvm2phys,weights_lgr_index_all_fvm2phys,&
533 1536 : weights_all_phys2fvm,weights_eul_index_all_phys2fvm,weights_lgr_index_all_phys2fvm,&
534 1536 : jall_fvm2phys,jall_phys2fvm)
535 : use dimensions_mod , only: ngpc,nelemd
536 : use fvm_overlap_mod , only: compute_weights_cell
537 : use element_mod , only: element_t
538 : type(element_t) , dimension(nelemd), intent(in) :: elem
539 : type (fvm_struct), dimension(nelemd), intent(in) :: fvm
540 : integer , intent(in) :: fvm_nc, phys_nc, irecons
541 : real (kind=r8) :: dalpha,dbeta
542 3072 : real (kind=r8), dimension(0:phys_nc+2):: xgno_phys,ygno_phys
543 3072 : real (kind=r8), dimension(0:fvm_nc+2) :: xgno_fvm,ygno_fvm
544 :
545 : real (kind=r8), dimension(ngpc):: gauss_weights, abscissae !dimension(ngauss)
546 :
547 : integer :: i,j,h
548 : integer, parameter :: nvertex = 4
549 : real (kind=r8), dimension(nvertex) :: xcell,ycell
550 :
551 : real (kind=r8) , dimension(num_weights_fvm2phys,irecons,nelemd),intent(out) :: weights_all_fvm2phys
552 : integer, dimension(num_weights_fvm2phys,2,nelemd),intent(out) :: weights_eul_index_all_fvm2phys
553 : integer, dimension(num_weights_fvm2phys,2,nelemd),intent(out) :: weights_lgr_index_all_fvm2phys
554 :
555 : real (kind=r8) , dimension(num_weights_phys2fvm,irecons,nelemd),intent(out) :: weights_all_phys2fvm
556 : integer, dimension(num_weights_phys2fvm,2,nelemd),intent(out) :: weights_eul_index_all_phys2fvm
557 : integer, dimension(num_weights_phys2fvm,2,nelemd),intent(out) :: weights_lgr_index_all_phys2fvm
558 :
559 : integer , dimension(nelemd) ,intent(out) :: jall_fvm2phys,jall_phys2fvm
560 :
561 : integer, parameter :: jmax_segments_cell = 50
562 3072 : real (kind=r8) , dimension(jmax_segments_cell,irecons) :: weights_cell
563 : integer , dimension(jmax_segments_cell,2) :: weights_eul_index_cell
564 : integer :: jcollect_cell,ie
565 3072 : real(kind=r8), dimension(phys_nc,phys_nc) :: phys_area, factor
566 3072 : real(kind=r8), dimension(fvm_nc,fvm_nc) :: fvm_area
567 :
568 1536 : xgno_phys(0) = -1D20; xgno_phys(phys_nc+2) = 1D20
569 1536 : xgno_fvm(0) = -1D20; xgno_fvm(fvm_nc+2) = 1D20
570 12336 : do ie=1,nelemd
571 10800 : dalpha = abs(elem(ie)%corners(1)%x-elem(ie)%corners(2)%x)/phys_nc !in alpha
572 10800 : dbeta = abs(elem(ie)%corners(1)%y-elem(ie)%corners(4)%y)/phys_nc !in beta
573 54000 : do i=1,phys_nc+1
574 43200 : xgno_phys(i) = tan(elem(ie)%corners(1)%x+(i-1)*dalpha)
575 54000 : ygno_phys(i) = tan(elem(ie)%corners(1)%y+(i-1)*dbeta )
576 : end do
577 :
578 10800 : dalpha = abs(elem(ie)%corners(1)%x-elem(ie)%corners(2)%x)/fvm_nc !in alpha
579 10800 : dbeta = abs(elem(ie)%corners(1)%y-elem(ie)%corners(4)%y)/fvm_nc !in beta
580 54000 : do i=1,fvm_nc+1
581 43200 : xgno_fvm(i) = tan(elem(ie)%corners(1)%x+(i-1)*dalpha)
582 54000 : ygno_fvm(i) = tan(elem(ie)%corners(1)%y+(i-1)*dbeta )
583 : end do
584 :
585 : !
586 : ! compute area using line-integrals
587 : !
588 : ! do j=1,phys_nc
589 : ! do i=1,phys_nc
590 : ! da_phys(i,j) = (I_00(xgno_phys(i+1),ygno_phys(j+1)) - I_00(xgno_phys(i ),ygno_phys(j+1)) + &
591 : ! I_00(xgno_phys(i ),ygno_phys(j )) - I_00(xgno_phys(i+1),ygno_phys(j )))
592 : ! end do
593 : ! end do
594 : !
595 : ! do j=1,fvm_nc
596 : ! do i=1,fvm_nc
597 : ! da_fvm(i,j) = (I_00(xgno_fvm(i+1),ygno_fvm(j+1)) - I_00(xgno_fvm(i ),ygno_fvm(j+1)) + &
598 : ! I_00(xgno_fvm(i ),ygno_fvm(j )) - I_00(xgno_fvm(i+1),ygno_fvm(j )))
599 : ! end do
600 : ! end do
601 :
602 10800 : gauss_weights = 0.0D0; abscissae=0.0D0!not used since line-segments are parallel to coordinate
603 :
604 10800 : jall_fvm2phys(ie)=1
605 43200 : do j=1,phys_nc
606 140400 : do i=1,phys_nc
607 97200 : xcell(1) = xgno_phys(i) ; ycell(1) = ygno_phys(j)
608 97200 : xcell(2) = xgno_phys(i) ; ycell(2) = ygno_phys(j+1)
609 97200 : xcell(3) = xgno_phys(i+1); ycell(3) = ygno_phys(j+1)
610 97200 : xcell(4) = xgno_phys(i+1); ycell(4) = ygno_phys(j)
611 :
612 : call compute_weights_cell(nvertex,.true.,&
613 : xcell,ycell,i,j,irecons,xgno_fvm,ygno_fvm,0,fvm_nc+2,&
614 : 1,fvm_nc+1,1,fvm_nc+1,&
615 : ngpc,gauss_weights,abscissae,&
616 97200 : weights_cell,weights_eul_index_cell,jcollect_cell,jmax_segments_cell)
617 :
618 129600 : if (jcollect_cell>0) then
619 0 : weights_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,:,ie) = &
620 1263600 : weights_cell(1:jcollect_cell,:)!/fvm(ie)%area_sphere_physgrid(i,j)!da_phys(i,j)
621 :
622 0 : weights_eul_index_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,:,ie) = &
623 486000 : weights_eul_index_cell(1:jcollect_cell,:)
624 194400 : weights_lgr_index_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,1,ie) = i
625 194400 : weights_lgr_index_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,2,ie) = j
626 97200 : jall_fvm2phys(ie) = jall_fvm2phys(ie)+jcollect_cell
627 : endif
628 : end do
629 : enddo
630 10800 : jall_fvm2phys(ie)=jall_fvm2phys(ie)-1
631 : !
632 : ! make sure sum of area overlap weights exactly match fvm%%area_sphere_physgrid
633 : !
634 140400 : phys_area = 0.0_r8
635 108000 : do h=1,jall_fvm2phys(ie)
636 97200 : i = weights_lgr_index_all_fvm2phys(h,1,ie); j = weights_lgr_index_all_fvm2phys(h,2,ie)
637 108000 : phys_area(i,j) = phys_area(i,j) +weights_all_fvm2phys(h,1,ie)
638 : end do
639 140400 : factor(:,:) = fvm(ie)%area_sphere_physgrid(:,:)/phys_area(:,:)
640 108000 : do h=1,jall_fvm2phys(ie)
641 97200 : i = weights_lgr_index_all_fvm2phys(h,1,ie); j = weights_lgr_index_all_fvm2phys(h,2,ie)
642 108000 : weights_all_fvm2phys(h,1,ie) = weights_all_fvm2phys(h,1,ie)*factor(i,j)
643 : end do
644 :
645 10800 : jall_phys2fvm(ie)=1
646 43200 : do j=1,fvm_nc
647 140400 : do i=1,fvm_nc
648 97200 : xcell(1) = xgno_fvm(i) ; ycell(1) = ygno_fvm(j)
649 97200 : xcell(2) = xgno_fvm(i) ; ycell(2) = ygno_fvm(j+1)
650 97200 : xcell(3) = xgno_fvm(i+1); ycell(3) = ygno_fvm(j+1)
651 97200 : xcell(4) = xgno_fvm(i+1); ycell(4) = ygno_fvm(j)
652 :
653 : call compute_weights_cell(nvertex,.true.,&
654 : xcell,ycell,i,j,irecons,xgno_phys,ygno_phys,0,phys_nc+2,&
655 : 1,phys_nc+1,1,phys_nc+1,&
656 : ngpc,gauss_weights,abscissae,&
657 97200 : weights_cell,weights_eul_index_cell,jcollect_cell,jmax_segments_cell)
658 :
659 129600 : if (jcollect_cell>0) then
660 0 : weights_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,:,ie) &
661 1263600 : = weights_cell(1:jcollect_cell,:)!/fvm(ie)%area_sphere(i,j)!da_fvm(i,j)
662 :
663 0 : weights_eul_index_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,:,ie) = &
664 486000 : weights_eul_index_cell(1:jcollect_cell,:)
665 194400 : weights_lgr_index_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,1,ie) = i
666 194400 : weights_lgr_index_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,2,ie) = j
667 97200 : jall_phys2fvm(ie) = jall_phys2fvm(ie)+jcollect_cell
668 : endif
669 : end do
670 : enddo
671 10800 : jall_phys2fvm(ie)=jall_phys2fvm(ie)-1
672 : !
673 : ! make sure sum of area overlap weights exactly matches fvm%%area_sphere_physgrid
674 : !
675 140400 : fvm_area = 0.0_r8
676 108000 : do h=1,jall_phys2fvm(ie)
677 97200 : i = weights_lgr_index_all_phys2fvm(h,1,ie); j = weights_lgr_index_all_phys2fvm(h,2,ie)
678 108000 : fvm_area(i,j) = fvm_area(i,j) +weights_all_phys2fvm(h,1,ie)
679 : end do
680 140400 : fvm_area(:,:) = fvm(ie)%area_sphere(:,:)/fvm_area(:,:)
681 109536 : do h=1,jall_phys2fvm(ie)
682 97200 : i = weights_lgr_index_all_phys2fvm(h,1,ie); j = weights_lgr_index_all_phys2fvm(h,2,ie)
683 108000 : weights_all_phys2fvm(h,1,ie) = weights_all_phys2fvm(h,1,ie)*fvm_area(i,j)
684 : end do
685 : end do
686 1536 : end subroutine fvm2phys_init
687 : end module dp_mapping
|
