Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grids.f90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 17 months ago
Moved oasis-mct_5.0 in oasis3-mct/branches directory.
File size: 47.7 KB

Line
1	!****
2	! ************************
3	! * OASIS MODULE *
4	! * ------------ *
5	! ************************
6	!****
7	!***********************************************************************
8	! This module belongs to the SCRIP library. It is modified to run
9	! within OASIS.
10	! Modifications:
11	! - routine does not read SCRIP grid description files, but gets
12	! arrays from the calling routine
13	! - unit conversion : only from degrees (OASIS unit) to radians
14	! - some allocated array will be freed in the end to allow multiple
15	! calls of SCRIP
16	! - map-methods and restriction-types are written in capital letters
17	! - added bin definition for reduced grid
18	!
19	! Modified by V. Gayler, M&D 20.09.2001
20	! Modified by D. Declat, CERFACS 27.06.2002
21	!***********************************************************************
22	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23	!
24	! This module reads in and initializes two grids for remapping.
25	! NOTE: grid1 must be the master grid -- the grid that determines
26	! which cells participate (e.g. land mask) and the fractional
27	! area of grid2 cells that participate in the remapping.
28	!
29	!-----------------------------------------------------------------------
30	!
31	! CVS:$Id: grids.f 2826 2010-12-10 11:14:21Z valcke $
32	!
33	! Copyright (c) 1997, 1998 the Regents of the University of
34	! California.
35	!
36	! This software and ancillary information (herein called software)
37	! called SCRIP is made available under the terms described here.
38	! The software has been approved for release with associated
39	! LA-CC Number 98-45.
40	!
41	! Unless otherwise indicated, this software has been authored
42	! by an employee or employees of the University of California,
43	! operator of the Los Alamos National Laboratory under Contract
44	! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S.
45	! Government has rights to use, reproduce, and distribute this
46	! software. The public may copy and use this software without
47	! charge, provided that this Notice and any statement of authorship
48	! are reproduced on all copies. Neither the Government nor the
49	! University makes any warranty, express or implied, or assumes
50	! any liability or responsibility for the use of this software.
51	!
52	! If software is modified to produce derivative works, such modified
53	! software should be clearly marked, so as not to confuse it with
54	! the version available from Los Alamos National Laboratory.
55	!
56	!***********************************************************************
57
58	module grids
59
60	!-----------------------------------------------------------------------
61
62	use kinds_mod ! defines data types
63	use constants ! common constants
64	use iounits ! I/O unit manager
65	USE mod_oasis_flush
66
67	implicit none
68
69	!-----------------------------------------------------------------------
70	!
71	! variables that describe each grid
72	!
73	!-----------------------------------------------------------------------
74
75	integer (kind=int_kind), save :: grid1_size, grid2_size, & ! total points on each grid
76	grid1_rank, grid2_rank, & ! rank of each grid
77	grid1_corners, grid2_corners ! number of corners
78	! for each grid cell
79
80	integer (kind=int_kind), dimension(:), allocatable, save :: grid1_dims, grid2_dims ! size of each grid dimension
81
82	character(char_len), save :: grid1_name, grid2_name ! name for each grid
83
84	character (char_len), save :: grid1_units, & ! units for grid coords (degs/radians)
85	grid2_units ! units for grid coords
86
87	real (kind=dbl_kind), parameter :: deg2rad = pi/180. ! conversion for deg to rads
88
89	!-----------------------------------------------------------------------
90	!
91	! grid coordinates and masks
92	!
93	!-----------------------------------------------------------------------
94
95	logical (kind=log_kind), dimension(:), allocatable, save :: grid1_mask, & ! flag which cells participate
96	grid2_mask ! flag which cells participate
97
98	real (kind=dbl_kind), dimension(:), allocatable, save :: grid1_center_lat, & ! lat/lon coordinates for
99	grid1_center_lon, & ! each grid center in radians
100	grid2_center_lat, &
101	grid2_center_lon, &
102	grid1_area, & ! tot area of each grid1 cell
103	grid2_area, & ! tot area of each grid2 cell
104	grid1_area_in, & ! area of grid1 cell from file
105	grid2_area_in, & ! area of grid2 cell from file
106	grid1_frac, & ! fractional area of grid cells
107	grid2_frac ! participating in remapping
108
109	real (kind=dbl_kind), dimension(:,:), allocatable, save :: grid1_corner_lat, & ! lat/lon coordinates for
110	grid1_corner_lon, & ! each grid corner in radians
111	grid2_corner_lat, &
112	grid2_corner_lon
113
114	logical (kind=log_kind), save :: luse_grid_centers, & ! use centers for bounding boxes
115	luse_grid1_area, & ! use area from grid file
116	luse_grid2_area, & ! use area from grid file
117	lstore_grid1_area, & ! store area from grid file
118	lstore_grid2_area ! store area from grid file
119
120	real (kind=dbl_kind), dimension(:,:), allocatable, save :: grid1_bound_box, & ! lat/lon bounding box for use
121	grid2_bound_box ! in restricting grid searches
122
123	integer, dimension(:), allocatable, save :: grid1_bbox_per, & ! bbox position wrt [0,2pi]
124	grid2_bbox_per
125
126	!-----------------------------------------------------------------------
127	!
128	! bins for restricting searches
129	!
130	!-----------------------------------------------------------------------
131
132	character (char_len), save :: restrict_type ! type of bins to use
133
134	integer (kind=int_kind), save :: num_srch_bins, & ! num of bins for restricted srch
135	num_srch_red ! num of bins for reduced case
136
137	integer (kind=int_kind), dimension(:,:), allocatable, save :: bin_addr1, & ! min,max adds for grid1 cells in this lat bin
138	bin_addr2 ! min,max adds for grid2 cells in this lat bin
139	integer (kind=int_kind), dimension(:,:), allocatable, save :: bin_addr1_r ! min,max adds for red grid1 cells
140
141	real(kind=dbl_kind), dimension(:,:), allocatable, save :: bin_lats, & ! min,max latitude for each search bin
142	bin_lons ! min,max longitude for each search bin
143	real(kind=dbl_kind), dimension(:,:), allocatable, save :: bin_lats_r ! min,max lat for each search bin for red grid
144
145	!***********************************************************************
146
147	contains
148
149	!***********************************************************************
150
151	subroutine grid_init(m_method, rst_type, n_srch_bins, &
152	src_size, dst_size, src_dims, dst_dims, &
153	src_rank, dst_rank, ncrn_src, ncrn_dst, &
154	src_mask, dst_mask, src_name, dst_name, &
155	src_lat, src_lon, dst_lat, dst_lon, &
156	src_corner_lat, src_corner_lon, &
157	dst_corner_lat, dst_corner_lon, &
158	lstore_src_area, src_area, &
159	lstore_dst_area, dst_area, &
160	ilogunit, ilogprt)
161
162	!-----------------------------------------------------------------------
163	!
164	! this routine gets grid info from routine scriprmp and makes any
165	! necessary changes (e.g. for 0,2pi longitude range)
166	!
167	!-----------------------------------------------------------------------
168
169	!-----------------------------------------------------------------------
170	!
171	! input variables
172	!
173	!-----------------------------------------------------------------------
174
175	integer (kind=int_kind), intent(in) :: n_srch_bins, & ! num of bins for restricted srch
176	src_size, & ! source grid size
177	dst_size, & ! target grid size
178	src_rank, & ! source grid rank
179	dst_rank, & ! target grid rank
180	src_dims(src_rank),& ! source grid dimensions
181	dst_dims(dst_rank),& ! target grid dimensions
182	ncrn_src, & ! number of corners of a source grid cell
183	ncrn_dst, & ! number of corners of a target grid cell
184	src_mask(src_size),& ! source grid mask (master mask)
185	dst_mask(dst_size) ! target grid mask
186
187	integer(kind=int_kind), intent(in), optional :: ilogunit
188	integer(kind=int_kind), intent(in), optional :: ilogprt
189
190	character(len=*), intent(in) :: m_method, & ! remapping method
191	rst_type, & ! restriction type
192	src_name, & ! source grid name
193	dst_name ! target grid name
194
195	logical, intent(in) :: lstore_src_area, & ! store source area from grid file
196	lstore_dst_area ! store dest area from grid file
197
198	real (kind=real_kind), intent (in) :: src_lat(src_size), & ! source grid latitudes
199	src_lon(src_size), & ! sourde grid longitudes
200	dst_lat(dst_size), & ! target grid latitudes
201	dst_lon(dst_size), & ! target grid longitudes
202	src_corner_lat(ncrn_src,src_size), &
203	src_corner_lon(ncrn_src,src_size), &
204	dst_corner_lat(ncrn_dst,dst_size), &
205	dst_corner_lon(ncrn_dst,dst_size), &
206	src_area(src_size),& ! true source grid areas
207	dst_area(dst_size) ! true target grid areas
208
209	!-----------------------------------------------------------------------
210	!
211	! local variables
212	!
213	!-----------------------------------------------------------------------
214
215	integer (kind=int_kind) :: n, & ! loop counter
216	next_n, nnext_n, &
217	nele, & ! element loop counter
218	i,j, & ! logical 2d addresses
219	ip1,jp1, n_add, e_add, ne_add, nx, ny
220
221	real (kind=dbl_kind) :: dlat, dlon, & ! lat/lon intervals for search bins
222	vec1_lat, vec1_lon, & ! vectors and cross products used
223	vec2_lat, vec2_lon, & ! during grid check
224	vecc_lat, vecc_lon, & ! during grid check
225	cross_product, cross_product_ctr
226
227	logical (kind=log_kind) :: ll_error
228
229	real (kind=dbl_kind), dimension(4) ::tmp_lats, tmp_lons ! temps for computing bounding boxes
230	character(len=*),parameter :: subname = 'scrip:grid_init '
231
232	logical (kind=log_kind) :: lp_debug_grid_corners = .false.
233
234	!
235	!-----------------------------------------------------------------------
236	!
237	if (present(ilogunit)) then
238	nulou = ilogunit
239	endif
240
241	if (present(ilogprt)) then
242	nlogprt = ilogprt
243	endif
244
245	IF (nlogprt .GE. 2) THEN
246	WRITE (UNIT = nulou,FMT = *)' '
247	WRITE (UNIT = nulou,FMT = *)'Entering routine grid_init'
248	WRITE (UNIT = nulou,FMT = *)' '
249	CALL OASIS_FLUSH_SCRIP(nulou)
250	ENDIF
251	!
252	!-----------------------------------------------------------------------
253	!
254	! allocate grid coordinates/masks and read data
255	!
256	!-----------------------------------------------------------------------
257
258	! write(nulou,*) subname,trim(m_method)
259	! write(nulou,*) subname,src_size,dst_size,src_rank,dst_rank
260
261	lstore_grid1_area = .false.
262	lstore_grid2_area = .false.
263
264	select case(m_method)
265	case ('CONSERV')
266	luse_grid_centers = .false.
267	lstore_grid1_area = lstore_src_area
268	lstore_grid2_area = lstore_dst_area
269	case ('BILINEAR','BILINEARNF')
270	luse_grid_centers = .true.
271	case ('BICUBIC','BICUBICNF')
272	luse_grid_centers = .true.
273	case ('DISTWGT','DISTWGTNF')
274	luse_grid_centers = .true.
275	case ('GAUSWGT','GAUSWGTNF')
276	luse_grid_centers = .true.
277	case ('LOCCUNIF', 'LOCCDIST', 'LOCCGAUS')
278	luse_grid_centers = .true.
279	lstore_grid1_area = lstore_src_area
280	lstore_grid2_area = lstore_dst_area
281	case default
282	stop 'unknown mapping method'
283	end select
284
285	allocate( grid1_mask (src_size), &
286	grid2_mask (dst_size), &
287	grid1_center_lat(src_size), &
288	grid1_center_lon(src_size), &
289	grid2_center_lat(dst_size), &
290	grid2_center_lon(dst_size), &
291	grid1_area (src_size), &
292	grid2_area (dst_size), &
293	grid1_frac (src_size), &
294	grid2_frac (dst_size), &
295	grid1_dims (src_rank), &
296	grid2_dims (dst_rank), &
297	grid1_bound_box (4 , src_size), &
298	grid2_bound_box (4 , dst_size), &
299	grid1_bbox_per (src_size), &
300	grid2_bbox_per (dst_size))
301
302	if (lstore_grid1_area) allocate (grid1_area_in(src_size))
303	if (lstore_grid2_area) allocate (grid2_area_in(dst_size))
304
305	if (.not. luse_grid_centers) then
306	allocate( grid1_corner_lat(ncrn_src, src_size), &
307	grid1_corner_lon(ncrn_src, src_size), &
308	grid2_corner_lat(ncrn_dst, dst_size), &
309	grid2_corner_lon(ncrn_dst, dst_size))
310	endif
311
312	!-----------------------------------------------------------------------
313	!
314	! copy input data to module data
315	!
316	!-----------------------------------------------------------------------
317
318	restrict_type = rst_type
319	num_srch_bins = n_srch_bins
320	grid1_size = src_size
321	grid2_size = dst_size
322	grid1_dims = src_dims
323	grid2_dims = dst_dims
324	grid1_rank = src_rank
325	grid2_rank = dst_rank
326	grid1_corners = ncrn_src
327	grid2_corners = ncrn_dst
328	grid1_name = src_name
329	grid2_name = dst_name
330	grid1_center_lat = src_lat
331	grid1_center_lon = src_lon
332	grid2_center_lat = dst_lat
333	grid2_center_lon = dst_lon
334
335	if (.not. luse_grid_centers) then
336	grid1_corner_lat = src_corner_lat
337	grid1_corner_lon = src_corner_lon
338	grid2_corner_lat = dst_corner_lat
339	grid2_corner_lon = dst_corner_lon
340	endif
341
342	if (lstore_grid1_area) then
343	grid1_area_in = src_area
344	endif
345	if (lstore_grid2_area) then
346	grid2_area_in = dst_area
347	endif
348
349	grid1_area = zero
350	grid1_frac = zero
351	grid2_area = zero
352	grid2_frac = zero
353
354	!-----------------------------------------------------------------------
355	!
356	! initialize logical mask and convert lat/lon units if required
357	!
358	!-----------------------------------------------------------------------
359	where (src_mask == 1)
360	grid1_mask = .true.
361	elsewhere
362	grid1_mask = .false.
363	endwhere
364
365	where (dst_mask == 1)
366	grid2_mask = .true.
367	elsewhere
368	grid2_mask = .false.
369	endwhere
370
371	!
372	!* -- convert unit from degrees (OASIS unit) to radians
373	!
374	grid1_center_lat = grid1_center_lat*deg2rad
375	grid1_center_lon = grid1_center_lon*deg2rad
376	grid2_center_lat = grid2_center_lat*deg2rad
377	grid2_center_lon = grid2_center_lon*deg2rad
378
379	if (.not. luse_grid_centers) then
380	grid1_corner_lat = grid1_corner_lat*deg2rad
381	grid1_corner_lon = grid1_corner_lon*deg2rad
382	grid2_corner_lat = grid2_corner_lat*deg2rad
383	grid2_corner_lon = grid2_corner_lon*deg2rad
384	endif
385
386	grid1_units='radians'
387	grid2_units='radians'
388
389	!-----------------------------------------------------------------------
390	!
391	! change corners periodicity according to center position
392	!
393	!-----------------------------------------------------------------------
394
395	if (.not. luse_grid_centers) then
396
397	do nele = 1, grid1_size
398	if (maxval(grid1_corner_lon(:,nele)) - &
399	&minval(grid1_corner_lon(:,nele)) > pi) then
400	!AP write(nulou,*) 'WARNING fixing wrong corners periodicity for grid1, nele =',nele
401	do n = 1,ncrn_src
402	if (grid1_corner_lon(n,nele) < grid1_center_lon(nele) - pi) then
403	grid1_corner_lon(n,nele) = grid1_corner_lon(n,nele) + pi2
404	else
405	if (grid1_corner_lon(n,nele) > grid1_center_lon(nele) + pi) then
406	grid1_corner_lon(n,nele) = grid1_corner_lon(n,nele) - pi2
407	end if
408	end if
409	end do
410	end if
411	end do
412
413	do nele = 1, grid2_size
414	if (maxval(grid2_corner_lon(:,nele)) - &
415	&minval(grid2_corner_lon(:,nele)) > pi) then
416	!AP write(nulou,*) 'WARNING fixing wrong corners periodicity for grid2, nele =',nele
417	do n = 1,ncrn_dst
418	if (grid2_corner_lon(n,nele) < grid2_center_lon(nele) - pi) then
419	grid2_corner_lon(n,nele) = grid2_corner_lon(n,nele) + pi2
420	else
421	if (grid2_corner_lon(n,nele) > grid2_center_lon(nele) + pi) then
422	grid2_corner_lon(n,nele) = grid2_corner_lon(n,nele) - pi2
423	end if
424	end if
425	end do
426	end if
427	end do
428
429	end if
430
431	!-----------------------------------------------------------------------
432	!
433	! convert center longitudes to 0,2pi interval preserving corners
434	!
435	!-----------------------------------------------------------------------
436
437	do n = 1, grid1_size
438	if (grid1_center_lon(n) .gt. pi2) then
439	grid1_center_lon(n) = grid1_center_lon(n) - pi2
440	if (.not. luse_grid_centers) &
441	& grid1_corner_lon(:,n) = grid1_corner_lon(:,n) - pi2
442	end if
443	end do
444
445	do n = 1, grid1_size
446	if (grid1_center_lon(n) .lt. zero) then
447	grid1_center_lon(n) = grid1_center_lon(n) + pi2
448	if (.not. luse_grid_centers) &
449	& grid1_corner_lon(:,n) = grid1_corner_lon(:,n) + pi2
450	end if
451	end do
452
453	do n = 1, grid2_size
454	if (grid2_center_lon(n) .gt. pi2) then
455	grid2_center_lon(n) = grid2_center_lon(n) - pi2
456	if (.not. luse_grid_centers) &
457	& grid2_corner_lon(:,n) = grid2_corner_lon(:,n) - pi2
458	end if
459	end do
460
461	do n = 1, grid2_size
462	if (grid2_center_lon(n) .lt. zero) then
463	grid2_center_lon(n) = grid2_center_lon(n) + pi2
464	if (.not. luse_grid_centers) &
465	& grid2_corner_lon(:,n) = grid2_corner_lon(:,n) + pi2
466	end if
467	end do
468
469	!-----------------------------------------------------------------------
470	!
471	! make sure input latitude range is within the machine values
472	! for +/- pi/2
473	!
474	!-----------------------------------------------------------------------
475
476	where (grid1_center_lat > pih) grid1_center_lat = pih
477	where (grid1_center_lat < -pih) grid1_center_lat = -pih
478	where (grid2_center_lat > pih) grid2_center_lat = pih
479	where (grid2_center_lat < -pih) grid2_center_lat = -pih
480
481	if (.not. luse_grid_centers) then
482	where (grid1_corner_lat > pih) grid1_corner_lat = pih
483	where (grid1_corner_lat < -pih) grid1_corner_lat = -pih
484	where (grid2_corner_lat > pih) grid2_corner_lat = pih
485	where (grid2_corner_lat < -pih) grid2_corner_lat = -pih
486	endif
487
488	!-----------------------------------------------------------------------
489	!
490	! check corner counterclockwise numbering and consistency
491	!
492	!-----------------------------------------------------------------------
493
494	if (lp_debug_grid_corners .and. .not. luse_grid_centers) then
495
496	ll_error = .false.
497	do nele = 1, grid1_size
498	corner_loop1: do n=1,ncrn_src
499	next_n = mod(n,ncrn_src) + 1
500	nnext_n = mod(next_n,ncrn_src) + 1
501	vec1_lat = grid1_corner_lat(next_n, nele) - grid1_corner_lat(n,nele)
502	vec1_lon = grid1_corner_lon(next_n, nele) - grid1_corner_lon(n,nele)
503	vec2_lat = grid1_corner_lat(nnext_n,nele) - grid1_corner_lat(n,nele)
504	vec2_lon = grid1_corner_lon(nnext_n,nele) - grid1_corner_lon(n,nele)
505	vecc_lat = grid1_center_lat(nele) - grid1_corner_lat(n,nele)
506	vecc_lon = grid1_center_lon(nele) - grid1_corner_lon(n,nele)
507
508	cross_product = vec1_lonvec2_lat - vec2_lonvec1_lat
509	cross_product_ctr = vec1_lonvecc_lat - vecc_lonvec1_lat
510	!***
511	!*** if cross product is less than zero, this cell
512	!*** doesn't work
513	!***
514	if (cross_product < zero) then
515	if (grid1_mask(nele)) then
516	write(nulou,'(A,I7,A)') ' ERROR : valid cell ', nele,'&
517	& from src grid '//trim(grid1_name)//' does not respect corner order '
518	ll_error = .true.
519	else
520	write(nulou,'(A,I7,A)') ' WARNING : masked cell ', nele,'&
521	& from src grid '//trim(grid1_name)//' does not respect corner order '
522	end if
523	write(nulou,'(A,2F9.2)') ' grid_center lon and lat (deg) ',&
524	& grid1_center_lon(nele)/deg2rad,&
525	& grid1_center_lat(nele)/deg2rad
526	do next_n = 1,ncrn_src
527	write(nulou,'(A,I3,A,2F9.2)') ' grid_corner ',next_n,' lon and lat (deg) ',&
528	& grid1_corner_lon(next_n,nele)/deg2rad,&
529	& grid1_corner_lat(next_n,nele)/deg2rad
530	end do
531	exit corner_loop1
532	end if
533
534	if (cross_product_ctr < zero) then
535	if (abs(grid1_center_lat(nele)) .ge. 1.48) then
536	if (grid1_mask(nele)) then
537	write(nulou,'(A,I7,A)') ' WARNING : valid polar cell ', nele,'&
538	& from src grid '//trim(grid1_name)//' could lead to'
539	else
540	write(nulou,'(A,I7,A)') ' WARNING : masked polar cell ', nele,'&
541	& from src grid '//trim(grid1_name)//' could lead to'
542	end if
543	write(nulou,'(2A)') ' a loss in precision if',&
544	& ' a polar projection is not active in remap_conserv'
545	else
546	if (grid1_mask(nele)) then
547	write(nulou,'(A,I7,A)') ' WARNING : valid cell ', nele,'&
548	& from src grid '//trim(grid1_name)//' does not contain its center'
549	else
550	write(nulou,'(A,I7,A)') ' WARNING : masked cell ', nele,'&
551	& from src grid '//trim(grid1_name)//' does not contain its center'
552	end if
553	end if
554	write(nulou,'(A,2F9.2)') ' grid_center lon and lat (deg) ',&
555	& grid1_center_lon(nele)/deg2rad,&
556	& grid1_center_lat(nele)/deg2rad
557	do next_n = 1,ncrn_src
558	write(nulou,'(A,I3,A,2F9.2)') ' grid_corner ',next_n,' lon and lat (deg) ',&
559	& grid1_corner_lon(next_n,nele)/deg2rad,&
560	& grid1_corner_lat(next_n,nele)/deg2rad
561	end do
562	exit corner_loop1
563	end if
564
565	end do corner_loop1
566
567	end do
568
569	do nele = 1, grid2_size
570	corner_loop2: do n=1,ncrn_dst
571	next_n = mod(n,ncrn_dst) + 1
572	nnext_n = mod(next_n,ncrn_dst) + 1
573	vec1_lat = grid2_corner_lat(next_n, nele) - grid2_corner_lat(n,nele)
574	vec1_lon = grid2_corner_lon(next_n, nele) - grid2_corner_lon(n,nele)
575	vec2_lat = grid2_corner_lat(nnext_n,nele) - grid2_corner_lat(n,nele)
576	vec2_lon = grid2_corner_lon(nnext_n,nele) - grid2_corner_lon(n,nele)
577	vecc_lat = grid2_center_lat(nele) - grid2_corner_lat(n,nele)
578	vecc_lon = grid2_center_lon(nele) - grid2_corner_lon(n,nele)
579
580	cross_product = vec1_lonvec2_lat - vec2_lonvec1_lat
581	cross_product_ctr = vec1_lonvecc_lat - vecc_lonvec1_lat
582	!***
583	!*** if cross product is less than zero, this cell
584	!*** doesn't work
585	!***
586	if (cross_product < zero) then
587	if (grid2_mask(nele)) then
588	write(nulou,'(A,I7,A)') ' ERROR : valid cell ', nele,&
589	& ' from dst grid '//trim(grid2_name)//' does not respect corner order '
590	ll_error = .true.
591	else
592	write(nulou,'(A,I7,A)') ' WARNING : masked cell ', nele,&
593	& ' from dst grid '//trim(grid2_name)//' does not respect corner order '
594	end if
595	write(nulou,'(A,2F9.2)') ' grid_center lon and lat (deg) ',&
596	& grid2_center_lon(nele)/deg2rad,&
597	& grid2_center_lat(nele)/deg2rad
598	do next_n = 1,ncrn_dst
599	write(nulou,'(A,I3,A,2F9.2)') ' grid_corner ',next_n,' lon and lat (deg) ',&
600	& grid2_corner_lon(next_n,nele)/deg2rad,&
601	& grid2_corner_lat(next_n,nele)/deg2rad
602	end do
603	exit corner_loop2
604	end if
605
606	if (cross_product_ctr < zero) then
607	if (abs(grid2_center_lat(nele)) .ge. 1.48) then
608	if (grid2_mask(nele)) then
609	write(nulou,'(A,I7,A)') ' WARNING : valid polar cell ', nele,'&
610	& from dst grid '//trim(grid2_name)//' could lead to'
611	else
612	write(nulou,'(A,I7,A)') ' WARNING : masked polar cell ', nele,'&
613	& from dst grid '//trim(grid2_name)//' could lead to'
614	end if
615	write(nulou,'(2A)') ' a loss in precision if',&
616	& ' a polar projection is not active in remap_conserv'
617	else
618	if (grid2_mask(nele)) then
619	write(nulou,'(A,I7,A)') ' WARNING : valid cell ', nele,'&
620	& from dst grid '//trim(grid2_name)//' does not contain its center'
621	else
622	write(nulou,'(A,I7,A)') ' WARNING : masked cell ', nele,'&
623	& from dst grid '//trim(grid2_name)//' does not contain its center'
624	end if
625	end if
626	write(nulou,'(A,2F9.2)') ' grid_center lon and lat (deg) ',&
627	& grid2_center_lon(nele)/deg2rad,&
628	& grid2_center_lat(nele)/deg2rad
629	do next_n = 1,ncrn_dst
630	write(nulou,'(A,I3,A,2F9.2)') ' grid_corner ',next_n,' lon and lat (deg) ',&
631	& grid2_corner_lon(next_n,nele)/deg2rad,&
632	& grid2_corner_lat(next_n,nele)/deg2rad
633	end do
634	exit corner_loop2
635	end if
636
637	end do corner_loop2
638
639	end do
640
641	if (ll_error) stop ('Wrong corners order')
642
643	end if
644
645	!-----------------------------------------------------------------------
646	!
647	! compute bounding boxes for restricting future grid searches
648	!
649	!-----------------------------------------------------------------------
650
651	grid1_bbox_per(:) = 0
652	grid2_bbox_per(:) = 0
653
654	if (.not. luse_grid_centers) then
655
656	grid1_bound_box(1,:) = minval(grid1_corner_lat, DIM=1)
657	grid1_bound_box(2,:) = maxval(grid1_corner_lat, DIM=1)
658	grid1_bound_box(3,:) = minval(grid1_corner_lon, DIM=1)
659	grid1_bound_box(4,:) = maxval(grid1_corner_lon, DIM=1)
660
661	grid2_bound_box(1,:) = minval(grid2_corner_lat, DIM=1)
662	grid2_bound_box(2,:) = maxval(grid2_corner_lat, DIM=1)
663	grid2_bound_box(3,:) = minval(grid2_corner_lon, DIM=1)
664	grid2_bound_box(4,:) = maxval(grid2_corner_lon, DIM=1)
665	else
666
667	nx = grid1_dims(1)
668	ny = grid1_dims(2)
669
670	do n=1,grid1_size
671
672	!*** find N,S and NE points to this grid point
673
674	j = (n - 1)/nx +1
675	i = n - (j-1)*nx
676
677	if (i < nx) then
678	ip1 = i + 1
679	elseif (ny == 1) then
680	ip1 = i
681	else
682	!*** assume cyclic
683	ip1 = 1
684	!*** but if it is not, correct
685	dlat = abs(grid1_center_lat(n) - grid1_center_lat(n-1))
686	dlon = abs(grid1_center_lon(n) - grid1_center_lon(n-1))
687	e_add = (j - 1)*nx + ip1
688	if (abs(grid1_center_lat(e_add) - &
689	grid1_center_lat(n )) > 10.0 * dlat .or. &
690	abs(grid1_center_lon(e_add) - &
691	grid1_center_lon(n )) > 10.0 * dlon) then
692	ip1 = i
693	endif
694	endif
695
696	if (j < ny) then
697	jp1 = j+1
698	elseif (ny == 1) then
699	jp1 = j
700	else
701	!*** assume cyclic
702	jp1 = 1
703	!*** but if it is not, correct
704	dlat = abs(grid1_center_lat(n) - grid1_center_lat(n-nx))
705	dlon = abs(grid1_center_lon(n) - grid1_center_lon(n-nx))
706	n_add = (jp1 - 1)*nx + i
707	if (abs(grid1_center_lat(n_add) - &
708	grid1_center_lat(n )) > 10.0 * dlat .or. &
709	abs(grid1_center_lon(n_add) - &
710	grid1_center_lon(n )) > 10.0 * dlon) then
711	jp1 = j
712	endif
713	endif
714
715	n_add = (jp1 - 1)*nx + i
716	e_add = (j - 1)*nx + ip1
717	ne_add = (jp1 - 1)*nx + ip1
718
719	!*** find N,S and NE lat/lon coords and check bounding box
720
721	tmp_lats(1) = grid1_center_lat(n)
722	tmp_lats(2) = grid1_center_lat(e_add)
723	tmp_lats(3) = grid1_center_lat(ne_add)
724	tmp_lats(4) = grid1_center_lat(n_add)
725
726	tmp_lons(1) = grid1_center_lon(n)
727	tmp_lons(2) = grid1_center_lon(e_add)
728	tmp_lons(3) = grid1_center_lon(ne_add)
729	tmp_lons(4) = grid1_center_lon(n_add)
730
731	grid1_bound_box(1,n) = minval(tmp_lats)
732	grid1_bound_box(2,n) = maxval(tmp_lats)
733	grid1_bound_box(3,n) = minval(tmp_lons)
734	grid1_bound_box(4,n) = maxval(tmp_lons)
735	if (abs(grid1_bound_box(4,n) - grid1_bound_box(3,n)) > pi - 20*epsilon(1.)) then
736	grid1_bound_box(3,n) = maxval(tmp_lons)
737	grid1_bound_box(4,n) = minval(tmp_lons)
738	grid1_bbox_per(n) = 1
739	end if
740	end do
741
742	nx = grid2_dims(1)
743	ny = grid2_dims(2)
744
745	do n=1,grid2_size
746
747	!*** find N,S and NE points to this grid point
748
749	j = (n - 1)/nx +1
750	i = n - (j-1)*nx
751
752	if (i < nx) then
753	ip1 = i + 1
754	elseif (ny == 1) then
755	ip1 = i
756	else
757	!*** assume cyclic
758	ip1 = 1
759	!*** but if it is not, correct
760	dlat = abs(grid2_center_lat(n) - grid2_center_lat(n-1))
761	dlon = abs(grid2_center_lon(n) - grid2_center_lon(n-1))
762	e_add = (j - 1)*nx + ip1
763	if (abs(grid2_center_lat(e_add) - &
764	grid2_center_lat(n )) > 10.0 * dlat .or. &
765	abs(grid2_center_lon(e_add) - &
766	grid2_center_lon(n )) > 10.0 * dlon) then
767	ip1 = i
768	endif
769	endif
770
771	if (j < ny) then
772	jp1 = j+1
773	elseif (ny == 1) then
774	jp1 = j
775	else
776	!*** assume cyclic
777	jp1 = 1
778	!*** but if it is not, correct
779	dlat = abs(grid2_center_lat(n) - grid2_center_lat(n-nx))
780	dlon = abs(grid2_center_lon(n) - grid2_center_lon(n-nx))
781	n_add = (jp1 - 1)*nx + i
782	if (abs(grid2_center_lat(n_add) - &
783	grid2_center_lat(n )) > 10.0 * dlat .or. &
784	abs(grid2_center_lon(n_add) - &
785	grid2_center_lon(n )) > 10.0 * dlon) then
786	jp1 = j
787	endif
788	endif
789
790	n_add = (jp1 - 1)*nx + i
791	e_add = (j - 1)*nx + ip1
792	ne_add = (jp1 - 1)*nx + ip1
793
794	!*** find N,S and NE lat/lon coords and check bounding box
795
796	tmp_lats(1) = grid2_center_lat(n)
797	tmp_lats(2) = grid2_center_lat(e_add)
798	tmp_lats(3) = grid2_center_lat(ne_add)
799	tmp_lats(4) = grid2_center_lat(n_add)
800
801	tmp_lons(1) = grid2_center_lon(n)
802	tmp_lons(2) = grid2_center_lon(e_add)
803	tmp_lons(3) = grid2_center_lon(ne_add)
804	tmp_lons(4) = grid2_center_lon(n_add)
805
806	grid2_bound_box(1,n) = minval(tmp_lats)
807	grid2_bound_box(2,n) = maxval(tmp_lats)
808	grid2_bound_box(3,n) = minval(tmp_lons)
809	grid2_bound_box(4,n) = maxval(tmp_lons)
810	if (abs(grid2_bound_box(4,n) - grid2_bound_box(3,n)) > pi - 20*epsilon(1.)) then
811	grid2_bound_box(3,n) = maxval(tmp_lons)
812	grid2_bound_box(4,n) = minval(tmp_lons)
813	grid2_bbox_per(n) = 1
814	end if
815	end do
816
817	endif
818
819	do n = 1,grid1_size
820	if (grid1_bbox_per(n) == 0 .and. &
821	& abs(grid1_bound_box(4,n) - grid1_bound_box(3,n)) > pi - 20*epsilon(1.)) then
822	grid1_bound_box(3,n) = zero
823	grid1_bound_box(4,n) = pi2
824	if (grid1_center_lat(n) > pi/3.) grid1_bound_box(2,n) = pih
825	if (grid1_center_lat(n) < - pi/3.) grid1_bound_box(1,n) = -pih
826	end if
827	end do
828
829	do n = 1,grid2_size
830	if (grid2_bbox_per(n) == 0 .and. &
831	& abs(grid2_bound_box(4,n) - grid2_bound_box(3,n)) > pi - 20*epsilon(1.)) then
832	grid2_bound_box(3,n) = zero
833	grid2_bound_box(4,n) = pi2
834	if (grid2_center_lat(n) > pi/3.) grid2_bound_box(2,n) = pih
835	if (grid2_center_lat(n) < - pi/3.) grid2_bound_box(1,n) = -pih
836	end if
837	end do
838
839	where (grid1_bound_box(3,:) < 0) grid1_bbox_per(:) = -1
840	where (grid1_bound_box(4,:) > pi2) grid1_bbox_per(:) = 1
841	where (grid2_bound_box(3,:) < 0) grid2_bbox_per(:) = -1
842	where (grid2_bound_box(4,:) > pi2) grid2_bbox_per(:) = 1
843
844	!***
845	!*** try to check for cells that overlap poles or do not encompass centers
846	!***
847
848	do n = 1, grid1_size
849	if (grid1_center_lat(n) > grid1_bound_box(2,n)) then
850	if (grid1_center_lat(n) > pi/3.) then
851	grid1_bound_box(2,n) = pih
852	else
853	grid1_bound_box(2,n) = grid1_center_lat(n)
854	end if
855	end if
856	end do
857
858	do n = 1, grid1_size
859	if (grid1_center_lat(n) < grid1_bound_box(1,n)) then
860	if (grid1_center_lat(n) < - pi/3.) then
861	grid1_bound_box(1,n) = -pih
862	else
863	grid1_bound_box(1,n) = grid1_center_lat(n)
864	end if
865	end if
866
867	end do
868
869	do n = 1, grid2_size
870	if (grid2_center_lat(n) > grid2_bound_box(2,n)) then
871	if (grid2_center_lat(n) > pi/3.) then
872	grid2_bound_box(2,n) = pih
873	else
874	grid2_bound_box(2,n) = grid2_center_lat(n)
875	end if
876	end if
877	end do
878
879	do n = 1, grid2_size
880	if (grid2_center_lat(n) < grid2_bound_box(1,n)) then
881	IF (grid2_center_lat(n) < - pi/3.) THEN
882	grid2_bound_box(1,n) = -pih
883	else
884	grid2_bound_box(1,n) = grid2_center_lat(n)
885	end if
886	end if
887	end do
888
889	!-----------------------------------------------------------------------
890	!
891	! convert corner longitudes to 0,2pi interval (to be compliant with
892	! old grids.f90 arithmetics. This section could be dropped.
893	!
894	!-----------------------------------------------------------------------
895
896	if (.not. luse_grid_centers) then
897	where (grid1_corner_lon .gt. pi2) grid1_corner_lon = grid1_corner_lon - pi2
898	where (grid1_corner_lon .lt. zero) grid1_corner_lon = grid1_corner_lon + pi2
899	where (grid2_corner_lon .gt. pi2) grid2_corner_lon = grid2_corner_lon - pi2
900	where (grid2_corner_lon .lt. zero) grid2_corner_lon = grid2_corner_lon + pi2
901	endif
902
903	!-----------------------------------------------------------------------
904	!
905	! set up and assign address ranges to search bins in order to
906	! further restrict later searches
907	!
908	!-----------------------------------------------------------------------
909
910	select case (restrict_type)
911
912	case ('LATITUDE')
913	IF (nlogprt .GE. 2) &
914	write(nulou,*) 'Using latitude bins to restrict search.'
915
916	allocate(bin_addr1(2,num_srch_bins))
917	allocate(bin_addr2(2,num_srch_bins))
918	allocate(bin_lats (2,num_srch_bins))
919	allocate(bin_lons (2,num_srch_bins))
920
921	dlat = pi/num_srch_bins
922
923	do n=1,num_srch_bins
924	bin_lats(1,n) = (n-1)*dlat - pih
925	bin_lats(2,n) = n*dlat - pih
926	bin_lons(1,n) = zero
927	bin_lons(2,n) = pi2
928	bin_addr1(1,n) = grid1_size + 1
929	bin_addr1(2,n) = 0
930	bin_addr2(1,n) = grid2_size + 1
931	bin_addr2(2,n) = 0
932	end do
933
934	do nele=1,grid1_size
935	do n=1,num_srch_bins
936	if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and. &
937	grid1_bound_box(2,nele) >= bin_lats(1,n)) then
938	bin_addr1(1,n) = min(nele,bin_addr1(1,n))
939	bin_addr1(2,n) = max(nele,bin_addr1(2,n))
940	endif
941	end do
942	end do
943
944	do nele=1,grid2_size
945	do n=1,num_srch_bins
946	if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and. &
947	grid2_bound_box(2,nele) >= bin_lats(1,n)) then
948	bin_addr2(1,n) = min(nele,bin_addr2(1,n))
949	bin_addr2(2,n) = max(nele,bin_addr2(2,n))
950	endif
951	end do
952	end do
953
954	case ('LATLON')
955	IF (nlogprt .GE. 2) &
956	write(nulou,*) 'Using lat/lon boxes to restrict search.'
957
958	dlat = pi /num_srch_bins
959	dlon = pi2/num_srch_bins
960
961	allocate(bin_addr1(2,num_srch_bins*num_srch_bins))
962	allocate(bin_addr2(2,num_srch_bins*num_srch_bins))
963	allocate(bin_lats (2,num_srch_bins*num_srch_bins))
964	allocate(bin_lons (2,num_srch_bins*num_srch_bins))
965
966	n = 0
967	do j=1,num_srch_bins
968	do i=1,num_srch_bins
969	n = n + 1
970
971	bin_lats(1,n) = (j-1)*dlat - pih
972	bin_lats(2,n) = j*dlat - pih
973	bin_lons(1,n) = (i-1)*dlon
974	bin_lons(2,n) = i*dlon
975	bin_addr1(1,n) = grid1_size + 1
976	bin_addr1(2,n) = 0
977	bin_addr2(1,n) = grid2_size + 1
978	bin_addr2(2,n) = 0
979	end do
980	end do
981
982	num_srch_bins = num_srch_bins**2
983
984	do nele=1,grid1_size
985	do n=1,num_srch_bins
986	if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and. &
987	grid1_bound_box(2,nele) >= bin_lats(1,n) .and. &
988	grid1_bound_box(3,nele) <= bin_lons(2,n) .and. &
989	grid1_bound_box(4,nele) >= bin_lons(1,n)) then
990	bin_addr1(1,n) = min(nele,bin_addr1(1,n))
991	bin_addr1(2,n) = max(nele,bin_addr1(2,n))
992	endif
993	end do
994	end do
995
996	do nele=1,grid2_size
997	do n=1,num_srch_bins
998	if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and. &
999	grid2_bound_box(2,nele) >= bin_lats(1,n) .and. &
1000	grid2_bound_box(3,nele) <= bin_lons(2,n) .and. &
1001	grid2_bound_box(4,nele) >= bin_lons(1,n)) then
1002	bin_addr2(1,n) = min(nele,bin_addr2(1,n))
1003	bin_addr2(2,n) = max(nele,bin_addr2(2,n))
1004	endif
1005	end do
1006	end do
1007
1008	case ('REDUCED')
1009	!EM write(nulou,*) '\| Using reduced bins to restrict search. Reduced grids with'
1010	!EM write(nulou,) '\| a maximum of 500NBRBINS latitude circles can be treated'
1011
1012	allocate(bin_addr2(2,num_srch_bins))
1013	allocate(bin_lats (2,num_srch_bins))
1014	allocate(bin_lons (2,num_srch_bins))
1015
1016	dlat = pi/num_srch_bins
1017
1018	do n=1,num_srch_bins
1019	bin_lats(1,n) = (n-1)*dlat - pih
1020	bin_lats(2,n) = n*dlat - pih
1021	bin_lons(1,n) = zero
1022	bin_lons(2,n) = pi2
1023	bin_addr2(1,n) = grid2_size + 1
1024	bin_addr2(2,n) = 0
1025	end do
1026
1027	do nele=1,grid2_size
1028	do n=1,num_srch_bins
1029	if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and. &
1030	grid2_bound_box(2,nele) >= bin_lats(1,n)) then
1031	bin_addr2(1,n) = min(nele,bin_addr2(1,n))
1032	bin_addr2(2,n) = max(nele,bin_addr2(2,n))
1033	endif
1034	end do
1035	end DO
1036
1037	! Count number of search bins
1038	num_srch_red = 1
1039	do nele=1,grid1_size-1
1040	IF (grid1_center_lat(nele+1) /= grid1_center_lat(nele)) &
1041	num_srch_red = num_srch_red+1
1042	end DO
1043
1044	IF (nlogprt .GE. 2) &
1045	write(nulou,*) 'Number of automatic search bins : ', num_srch_red
1046
1047	allocate(bin_addr1_r(4,num_srch_red))
1048	allocate(bin_lats_r (2,num_srch_red))
1049
1050
1051	bin_addr1_r(1,1) = 0
1052	bin_lats_r(1,1) = grid1_center_lat(1)
1053	num_srch_red = 1
1054	do nele=1,grid1_size-1
1055	if (grid1_center_lat(nele+1) == grid1_center_lat(nele)) THEN
1056	bin_addr1_r(2,num_srch_red) = nele+1
1057	!EM why not simply : else ?
1058	else IF (grid1_center_lat(nele+1) /= grid1_center_lat(nele)) THEN
1059	bin_addr1_r(1,num_srch_red+1) = nele+1
1060	bin_lats_r(2,num_srch_red) = grid1_center_lat(nele+1)
1061	bin_lats_r(1,num_srch_red+1) = bin_lats_r(2,num_srch_red)
1062	num_srch_red = num_srch_red+1
1063	ENDIF
1064	end DO
1065
1066	DO nele = 1,num_srch_red-1
1067	bin_addr1_r(3,nele)=bin_addr1_r(1,nele+1)
1068	bin_addr1_r(4,nele)=bin_addr1_r(2,nele+1)
1069	enddo
1070	bin_addr1_r(3,num_srch_red) = bin_addr1_r(1,num_srch_red-1)
1071	bin_addr1_r(4,num_srch_red) = bin_addr1_r(2,num_srch_red-1)
1072
1073	case default
1074	stop 'unknown search restriction method'
1075	end select
1076	!
1077	IF (nlogprt .GE. 2) THEN
1078	WRITE (UNIT = nulou,FMT = *)' '
1079	WRITE (UNIT = nulou,FMT = *)'Leaving routine grid_init'
1080	WRITE (UNIT = nulou,FMT = *)' '
1081	CALL OASIS_FLUSH_SCRIP(nulou)
1082	ENDIF
1083	!
1084	!-----------------------------------------------------------------------
1085
1086	end subroutine grid_init
1087
1088	!***********************************************************************
1089
1090	subroutine free_grids
1091
1092	!-----------------------------------------------------------------------
1093	!
1094	! subroutine to deallocate allocated arrays
1095	!
1096	!-----------------------------------------------------------------------
1097	!
1098	IF (nlogprt .GE. 2) THEN
1099	WRITE (UNIT = nulou,FMT = *)' '
1100	WRITE (UNIT = nulou,FMT = *)'Entering routine free_grid'
1101	WRITE (UNIT = nulou,FMT = *)' '
1102	CALL OASIS_FLUSH_SCRIP(nulou)
1103	ENDIF
1104	!
1105	deallocate(grid1_mask, grid2_mask, &
1106	grid1_center_lat, grid1_center_lon, &
1107	grid2_center_lat, grid2_center_lon, &
1108	grid1_area, grid2_area, &
1109	grid1_frac, grid2_frac, &
1110	grid1_dims, grid2_dims)
1111
1112	if (lstore_grid1_area) deallocate (grid1_area_in)
1113	if (lstore_grid2_area) deallocate (grid2_area_in)
1114
1115	IF (restrict_TYPE == 'REDUCED') then
1116	deallocate( grid1_bound_box, grid2_bound_box, &
1117	grid1_bbox_per, grid2_bbox_per, &
1118	bin_addr1_r, bin_addr2, &
1119	bin_lons, bin_lats, &
1120	bin_lats_r)
1121	else
1122	deallocate( grid1_bound_box, grid2_bound_box, &
1123	grid1_bbox_per, grid2_bbox_per, &
1124	bin_addr1, bin_addr2, &
1125	bin_lats, bin_lons)
1126	endif
1127
1128	if (.not. luse_grid_centers) then
1129	deallocate(grid1_corner_lat, grid1_corner_lon, &
1130	grid2_corner_lat, grid2_corner_lon)
1131	ENDIF
1132	!
1133	IF (nlogprt .GE. 2) THEN
1134	WRITE (UNIT = nulou,FMT = *)' '
1135	WRITE (UNIT = nulou,FMT = *)'Leaving routine free_grid'
1136	WRITE (UNIT = nulou,FMT = *)' '
1137	CALL OASIS_FLUSH_SCRIP(nulou)
1138	ENDIF
1139	!-----------------------------------------------------------------------
1140
1141	end subroutine free_grids
1142
1143	!***********************************************************************
1144
1145	end module grids
1146
1147	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1148

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/lib/scrip/src/grids.f90 @ 6331

Download in other formats: