Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_global/haversine.f90 @ 5143

Last change on this file since 5143 was 3965, checked in by jan.polcher, 8 years ago

Merge with trunk at revision3959.
This includes all the developments made for CMIP6 and passage to XIOS2.
All conflicts are resolved and the code compiles.

But it still does not link because of an "undefined reference to `_intel_fast_memmove'"

File size: 41.4 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE : This module provides a few basic operations for convex polygons on the sphere and in longitude
3	! latitude coordinates. These operations are based on the haversine equations for great circle arcs.
4	!
5	! Definition of polygons : the basic assumption is that they describe grid boxes and are thus
6	! provided with more points than strictly necessary. Each polygon starts at index 1 on a vertex and
7	! alternates mid-points of segments and vertices. The mid-point of segments are kept as they are
8	! useful elements for the operations on the grid boxes.
9	!
10	! The module provides the following subroutine and functions :
11	! haversine_reglatlontoploy : Lists the polygons and all their attributes (centre point,
12	! neighbouring grids, indices in i and j on the original grid) for
13	! regular lon/lat grid.
14	! haversine_regxytoploy : As above but for grid boxes on the sphere projected onto a regular
15	! X/Y grid.
16	! haversine_singlepointpoly : Computes the polygone around a given point with a known area.
17	! haversine_polyheadings : Compute the heading out of the polygon for each vertex and mid-point of
18	! segment.
19	! haversine_polysort : Sort the polygons so that all points are in clockwise order.
20	! haversine_polyseglen : Compute the length of all segments of the polygon using great circles.
21	! haversine_polyseglen : Compute the length of all segments of the polygon on a regular lat lon grid.
22	! haversine_clockwise : Get the indices which sort a polygon in a clockwise order starting from an
23	! initial vertex.
24	! haversine_polyarea : Computes the area covered by a polygon on the sphere.
25	! haversine_laloarea : Compute the area for a regular lat lon grid.
26	! haversine_xyarea : Compute the area for the special case where the grid box size in x and y are already
27	! given by the projection.
28	! haversine_heading : Initial heading between start point and end point along a great circle.
29	! haversine_distance : Compute the distance between 2 points along the great circle.
30	! haversine_diffangle : Computes the angle between two headings.
31	! haversine_radialdis : Compute the coordinates found in a given heading and distance.
32	! haversine_dtor : Degrees to radians
33	! haversine_rtod : Radians to degrees
34	!
35	! CONTACT : jan.polcher@lmd.jussieu.fr
36	!
37	! LICENCE : IPSL (2016)
38	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
39	!
40	!>\BRIEF
41	!!
42	!! RECENT CHANGE(S): None
43	!!
44	!! REFERENCE(S) : None
45	!!
46	!_ ================================================================================================================================
47	MODULE haversine
48
49	USE defprec
50	USE constantes_var
51	USE module_llxy
52
53	IMPLICIT NONE
54
55	PRIVATE
56	PUBLIC :: haversine_reglatlontoploy, haversine_regxytoploy, haversine_singlepointploy, &
57	& haversine_polyheadings, haversine_polysort, &
58	& haversine_polyseglen, haversine_laloseglen, haversine_clockwise, haversine_polyarea, &
59	& haversine_laloarea, haversine_xyarea, haversine_diffangle, haversine_heading
60
61	CONTAINS
62	!! =============================================================================================================================
63	!! SUBROUTINE: haversine_reglatlontoploy
64	!!
65	!>\BRIEF Lists the polygons and all their attributes (centre point,
66	! neighbouring grids, indices in i and j on the original grid) for
67	! regular lon/lat grid.
68	!!
69	!! DESCRIPTION:
70	!!
71	!! \n
72	!_ ==============================================================================================================================
73	SUBROUTINE haversine_reglatlontoploy(iim, jjm, lon, lat, nbpt, index_loc, global, &
74	& nbseg, lonpoly, latpoly, center, neighb_loc, iorig, jorig)
75	!
76	! This subroutine constructs a series of polygons out of the grid boxes which are listed
77	! in the index_loc array. The polygons are ordered according to the indexing space and independently
78	! of the geographical coordinates.
79	!
80	! 0 interface
81	!
82	IMPLICIT NONE
83	!
84	! 0.1 input !
85	! Size of cartesian grid
86	INTEGER(i_std), INTENT(in) :: iim, jjm
87	! Longitudes on cartesian grid
88	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: lon
89	! Latitudes on cartesian grid
90	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: lat
91	! Size of the gathered domain
92	INTEGER(i_std), INTENT(in) :: nbpt
93	! Index of land point on 2D map (in local position)
94	INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: index_loc
95	! Is it a global grid ?
96	LOGICAL, INTENT(in) :: global
97	! Number of segments
98	INTEGER(i_std), INTENT(in) :: nbseg
99	!
100	! 0.2 Ouput
101	!
102	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: lonpoly
103	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: latpoly
104	REAL(r_std), DIMENSION(nbpt,2), INTENT(out) :: center
105	INTEGER(i_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: neighb_loc
106	INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: iorig, jorig
107	!
108	!
109	! 0.3 Local variables
110	!
111	INTEGER(i_std), DIMENSION(iim,jjm) :: correspondance
112	INTEGER(i_std) :: i, ip, jp
113	INTEGER(i_std) :: ipm1, ipp1, jpm1, jpp1
114	REAL(r_std) :: dlonm1, dlonp1, dlatm1, dlatp1
115	!
116	!
117	correspondance(:,:) = -1
118	DO i = 1, nbpt
119	!
120	! 1 find numbers of the latitude and longitude of each point
121	!
122	! index of latitude
123	jp = INT( (index_loc(i)-1) /iim ) + 1
124	! index of longitude
125	ip = index_loc(i) - ( jp-1 ) * iim
126	!
127	!correspondance(ip,jp) = kindex(i)
128	!
129	correspondance(ip,jp) = i
130	iorig(i)=ip
131	jorig(i)=jp
132	!
133	ENDDO
134	!
135	!
136	! Go through all the points and build the polygone which defines the grid area.
137	! This polygone will include the mid-point of each segment so that
138	! we can later compute the direction of the normal to the segment.
139	!
140	neighb_loc(:,:) = -1
141	!
142	DO i = 1, nbpt
143	!
144	ip = iorig(i)
145	jp = jorig(i)
146	!
147	ipm1 = ip-1
148	ipp1 = ip+1
149	jpm1 = jp-1
150	jpp1 = jp+1
151	!
152	! Compute the longitude increments
153	!
154	IF ( ipp1 <= iim ) THEN
155	dlonp1 = (lon(ipp1,jp)-lon(ip,jp))/2.0
156	ELSE IF ( ipm1 > 0 ) THEN
157	dlonp1 = (lon(ip,jp)-lon(ipm1,jp))/2.0
158	IF ( global ) ipp1=1
159	ELSE
160	dlonp1 = undef_sechiba
161	ENDIF
162	IF ( ipm1 > 0 ) THEN
163	dlonm1 = (lon(ip,jp)-lon(ipm1,jp))/2.0
164	ELSE IF ( ipp1 <= iim ) THEN
165	dlonm1 = (lon(ipp1,jp)-lon(ip,jp))/2.0
166	IF ( global ) ipm1=iim
167	ELSE
168	dlonm1 = undef_sechiba
169	ENDIF
170	!
171	! Verify that we have at least one valid longitude increment. Else we do not have enough
172	! points in the grid to estimate the position of the vertices in longitude.
173	!
174	IF ( dlonp1 >= undef_sechiba-1 ) dlonp1 = dlonm1
175	IF ( dlonm1 >= undef_sechiba-1 ) dlonm1 = dlonp1
176	IF ( dlonp1 >= undef_sechiba-1 .AND. dlonm1 >= undef_sechiba-1 ) THEN
177	CALL ipslerr(3, "haversine_reglatlontoploy", "There are not enogh point in longitude", &
178	& "to estimate the bounds of the polygone of the grid box.",&
179	& "Please choose a larger grid.")
180	ENDIF
181	!
182	! Compute the latitude increments
183	!
184	IF ( jpp1 <= jjm ) THEN
185	dlatp1 = (lat(ip,jpp1)-lat(ip,jp))/2.0
186	ELSE IF ( jpm1 > 0 ) THEN
187	dlatp1 = (lat(ip,jp)-lat(ip,jpm1))/2.0
188	ELSE
189	dlatp1 = undef_sechiba
190	ENDIF
191	IF ( jpm1 > 0 ) THEN
192	dlatm1 = (lat(ip,jp)-lat(ip,jpm1))/2.0
193	ELSE IF ( jpp1 <= jjm ) THEN
194	dlatm1 = (lat(ip,jpp1)-lat(ip,jp))/2.0
195	ELSE
196	dlatm1 = undef_sechiba
197	ENDIF
198	!
199	! Verify that we have at least one valid latitude increment. Else we do not have enough
200	! points in the grid to estimate the position of the vertices in latitude.
201	!
202	IF ( dlatp1 >= undef_sechiba-1 ) dlatp1 = dlatm1
203	IF ( dlatm1 >= undef_sechiba-1 ) dlatm1 = dlatp1
204	IF ( dlatp1 >= undef_sechiba-1 .AND. dlatm1 >= undef_sechiba-1 ) THEN
205	CALL ipslerr(3, "haversine_reglatlontoploy", "There are not enogh point in latitude", &
206	& "to estimate the bounds of the polygone of the grid box.",&
207	& "Please choose a larger grid.")
208	ENDIF
209	!
210	! The longitude of all elements of the polygone
211	!
212	lonpoly(i,1) = lon(ip,jp)-dlonm1
213	lonpoly(i,2) = lon(ip,jp)
214	lonpoly(i,3) = lon(ip,jp)+dlonp1
215	lonpoly(i,4) = lon(ip,jp)+dlonp1
216	lonpoly(i,5) = lon(ip,jp)+dlonp1
217	lonpoly(i,6) = lon(ip,jp)
218	lonpoly(i,7) = lon(ip,jp)-dlonm1
219	lonpoly(i,8) = lon(ip,jp)-dlonm1
220	!
221	! The longitude of all elements of the polygone
222	!
223	latpoly(i,1) = lat(ip,jp)-dlatp1
224	latpoly(i,2) = lat(ip,jp)-dlatp1
225	latpoly(i,3) = lat(ip,jp)-dlatp1
226	latpoly(i,4) = lat(ip,jp)
227	latpoly(i,5) = lat(ip,jp)+dlatm1
228	latpoly(i,6) = lat(ip,jp)+dlatm1
229	latpoly(i,7) = lat(ip,jp)+dlatm1
230	latpoly(i,8) = lat(ip,jp)
231	!
232	! Keep the center of the gridbox
233	!
234	center(i,1) = lon(ip,jp)
235	center(i,2) = lat(ip,jp)
236	!
237	! Get the neighbours when they exist in the list of land points
238	! There are no neighbours over the North or South poles.
239	!
240	IF ( ipm1 > 0 .AND. jpm1 > 0 ) THEN
241	neighb_loc(i,1) = correspondance(ipm1,jpm1)
242	ELSE
243	neighb_loc(i,1) = -2
244	ENDIF
245	IF ( jpm1 > 0 ) THEN
246	neighb_loc(i,2) = correspondance(ip,jpm1)
247	ELSE
248	neighb_loc(i,2) = -2
249	ENDIF
250	IF ( ipp1 <= iim .AND. jpm1 > 0 ) THEN
251	neighb_loc(i,3) = correspondance(ipp1,jpm1)
252	ELSE
253	neighb_loc(i,3) = -2
254	ENDIF
255	IF ( ipp1 <= iim ) THEN
256	neighb_loc(i,4) = correspondance(ipp1,jp)
257	ELSE
258	neighb_loc(i,4) = -2
259	ENDIF
260	IF ( ipp1 <= iim .AND. jpp1 <= jjm ) THEN
261	neighb_loc(i,5) = correspondance(ipp1,jpp1)
262	ELSE
263	neighb_loc(i,5) = -2
264	ENDIF
265	IF ( jpp1 <= jjm ) THEN
266	neighb_loc(i,6) = correspondance(ip,jpp1)
267	ELSE
268	neighb_loc(i,6) = -2
269	ENDIF
270	IF ( ipm1 > 0 .AND. jpp1 <= jjm ) THEN
271	neighb_loc(i,7) = correspondance(ipm1,jpp1)
272	ELSE
273	neighb_loc(i,7) = -2
274	ENDIF
275	IF ( ipm1 > 0 ) THEN
276	neighb_loc(i,8) = correspondance(ipm1,jp)
277	ELSE
278	neighb_loc(i,8) = -2
279	ENDIF
280	!
281	ENDDO
282	!
283	END SUBROUTINE haversine_reglatlontoploy
284	!!
285	!! =============================================================================================================================
286	!! SUBROUTINE: haversine_regxytoploy
287	!!
288	!>\BRIEF Same as haversine_reglatlontoploy but for grid boxes on the sphere projected onto a regular
289	! X/Y grid.
290	!!
291	!! DESCRIPTION: Keep in mind that in these projections the straight line assumed between 2 points is not always
292	!! the great circle. Thus the distance computed by the haversine formula might deviate a little.
293	!!
294	!! \n
295	!_ ==============================================================================================================================
296	SUBROUTINE haversine_regxytoploy(iim, jjm, lon, lat, nbpt, index_loc, proj, &
297	& nbseg, lonpoly, latpoly, center, neighb_loc, iorig, jorig)
298	!
299	! This subroutine constructs a series of polygons out of the grid boxes which are listed
300	! in the index array. This version will go directly from the indexing space to the coordinate as we know that
301	! we are dealing with a projection of the sphere to the plane where the regular grid is created.
302	! The polygons are ordered according to the indexing space and independently
303	! of the geographical coordinates.
304	!
305	! 0 interface
306	!
307	IMPLICIT NONE
308	!
309	! 0.1 input !
310	! Size of cartesian grid
311	INTEGER(i_std), INTENT(in) :: iim, jjm
312	! Longitudes on cartesian grid
313	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: lon
314	! Latitudes on cartesian grid
315	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: lat
316	! Size of the gathered domain
317	INTEGER(i_std), INTENT(in) :: nbpt
318	! Index of land point on 2D map (in local position)
319	INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: index_loc
320	! Projection ID
321	type (proj_info), DIMENSION(1) :: proj
322	! Number of segments
323	INTEGER(i_std), INTENT(in) :: nbseg
324	!
325	! 0.2 Ouput
326	!
327	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: lonpoly
328	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: latpoly
329	REAL(r_std), DIMENSION(nbpt,2), INTENT(out) :: center
330	INTEGER(i_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: neighb_loc
331	INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: iorig, jorig
332	!
333	!
334	! 0.3 Local variables
335	!
336	INTEGER(i_std), DIMENSION(iim,jjm) :: correspondance
337	INTEGER(i_std) :: i, ip, jp
338	INTEGER(i_std) :: ipm1, ipp1, jpm1, jpp1
339	REAL(r_std) :: dlonm1, dlonp1, dlatm1, dlatp1
340	!
341	!
342	!
343	correspondance(:,:) = -1
344	DO i = 1, nbpt
345	!
346	! 1 find numbers of the latitude and longitude of each point
347	!
348	! index of latitude
349	jp = INT( (index_loc(i)-1) /iim ) + 1
350	! index of longitude
351	ip = index_loc(i) - ( jp-1 ) * iim
352	!
353	!correspondance(ip,jp) = kindex(i)
354	!
355	correspondance(ip,jp) = i
356	iorig(i)=ip
357	jorig(i)=jp
358	!
359	ENDDO
360	!
361	! Go through all the points and build the polygone which defines the grid area.
362	! This polygone will include the mid-point of each segment so that
363	! we can later compute the direction of the normal to the segment.
364	!
365	neighb_loc(:,:) = -1
366	!
367	DO i = 1, nbpt
368	! index of latitude
369	jp = INT( (index_loc(i)-1) /iim ) + 1
370	! index of longitude
371	ip = index_loc(i) - ( jp-1 ) * iim
372	!
373	ipm1 = ip-1
374	ipp1 = ip+1
375	jpm1 = jp-1
376	jpp1 = jp+1
377	!
378	! Get the longitude and latitude throug the projection
379	! The range of possible values for projection depends on the module
380	! which defines these projections. For module_llxy the range is 0-203.
381	!
382	IF ( proj(1)%code > 0 .AND. proj(1)%code < 203 ) THEN
383	CALL ij_to_latlon(proj(1), ip-0.5, jp-0.5, latpoly(i,1), lonpoly(i,1))
384	CALL ij_to_latlon(proj(1), ip+0.0, jp-0.5, latpoly(i,2), lonpoly(i,2))
385	CALL ij_to_latlon(proj(1), ip+0.5, jp-0.5, latpoly(i,3), lonpoly(i,3))
386	CALL ij_to_latlon(proj(1), ip+0.5, jp+0.0, latpoly(i,4), lonpoly(i,4))
387	CALL ij_to_latlon(proj(1), ip+0.5, jp+0.5, latpoly(i,5), lonpoly(i,5))
388	CALL ij_to_latlon(proj(1), ip+0.0, jp+0.5, latpoly(i,6), lonpoly(i,6))
389	CALL ij_to_latlon(proj(1), ip-0.5, jp+0.5, latpoly(i,7), lonpoly(i,7))
390	CALL ij_to_latlon(proj(1), ip-0.5, jp+0.0, latpoly(i,8), lonpoly(i,8))
391	ELSE
392	CALL ipslerr(3, "haversine_regxytoploy", "Unknown projection code", &
393	& "Check proj(1)%code","")
394	ENDIF
395	!
396	! Keep the center of the gridbox
397	!
398	center(i,1) = lon(ip,jp)
399	center(i,2) = lat(ip,jp)
400	!
401	! Get the neighbours when they exist in the list of land points
402	! There are no neighbours over the North or South poles.
403	!
404	IF ( ipm1 > 0 .AND. jpm1 > 0 ) THEN
405	neighb_loc(i,1) = correspondance(ipm1,jpm1)
406	ELSE
407	neighb_loc(i,1) = -2
408	ENDIF
409	IF ( jpm1 > 0 ) THEN
410	neighb_loc(i,2) = correspondance(ip,jpm1)
411	ELSE
412	neighb_loc(i,2) = -2
413	ENDIF
414	IF ( ipp1 <= iim .AND. jpm1 > 0 ) THEN
415	neighb_loc(i,3) = correspondance(ipp1,jpm1)
416	ELSE
417	neighb_loc(i,3) = -2
418	ENDIF
419	IF ( ipp1 <= iim ) THEN
420	neighb_loc(i,4) = correspondance(ipp1,jp)
421	ELSE
422	neighb_loc(i,4) = -2
423	ENDIF
424	IF ( ipp1 <= iim .AND. jpp1 <= jjm ) THEN
425	neighb_loc(i,5) = correspondance(ipp1,jpp1)
426	ELSE
427	neighb_loc(i,5) = -2
428	ENDIF
429	IF ( jpp1 <= jjm ) THEN
430	neighb_loc(i,6) = correspondance(ip,jpp1)
431	ELSE
432	neighb_loc(i,6) = -2
433	ENDIF
434	IF ( ipm1 > 0 .AND. jpp1 <= jjm ) THEN
435	neighb_loc(i,7) = correspondance(ipm1,jpp1)
436	ELSE
437	neighb_loc(i,7) = -2
438	ENDIF
439	IF ( ipm1 > 0 ) THEN
440	neighb_loc(i,8) = correspondance(ipm1,jp)
441	ELSE
442	neighb_loc(i,8) = -2
443	ENDIF
444	!
445	ENDDO
446	!
447	END SUBROUTINE haversine_regxytoploy
448	!! =============================================================================================================================
449	!! SUBROUTINE: haversine_singlepointploy
450	!!
451	!>\BRIEF Lists the polygons and all their attributes (centre point,
452	! neighbouring grids, indices in i and j on the original grid) for
453	! a single point. A regular lon/lat grid is assumed.
454	!!
455	!! DESCRIPTION:
456	!!
457	!! \n
458	!_ ==============================================================================================================================
459	SUBROUTINE haversine_singlepointploy(iim, jjm, lon, lat, nbpt, index_loc, global, &
460	& nbseg, lonpoly, latpoly, center, neighb_loc, iorig, jorig)
461	!
462	! This subroutine constructs a series of polygons out of the grid boxe which is provided.
463	! The polygon is ordered according to the indexing space and independently
464	! of the geographical coordinates.
465	! It uses the same interface as haversine_reglatlontoploy but can only used with iim=jjm=1
466	!
467	! 0 interface
468	!
469	IMPLICIT NONE
470	!
471	! 0.1 input !
472	! Size of cartesian grid
473	INTEGER(i_std), INTENT(in) :: iim, jjm
474	! Longitudes on cartesian grid
475	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: lon
476	! Latitudes on cartesian grid
477	REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: lat
478	! Size of the gathered domain
479	INTEGER(i_std), INTENT(in) :: nbpt
480	! Index of land point on 2D map (in local position)
481	INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: index_loc
482	! Is it a global grid ?
483	LOGICAL, INTENT(in) :: global
484	! Number of segments
485	INTEGER(i_std), INTENT(in) :: nbseg
486	!
487	! 0.2 Ouput
488	!
489	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: lonpoly
490	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: latpoly
491	REAL(r_std), DIMENSION(nbpt,2), INTENT(out) :: center
492	INTEGER(i_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: neighb_loc
493	INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: iorig, jorig
494	!
495	!
496	! 0.3 Local variables
497	!
498	REAL(r_std) :: area, rlon, rlat
499	REAL(r_std), DIMENSION(2) :: coord
500	!
501	IF ( iim .NE. 1 .AND. jjm .NE. 1 ) THEN
502	CALL ipslerr(3, "haversine_singlepointploy", "Can only be used if iim=jjm=1", &
503	& "Please ensure this routine is called in the","right conditions")
504	ENDIF
505	!
506	iorig(1)=1
507	jorig(1)=1
508	!
509	area = 111111.0*111111.0
510	rlon= SQRT(area)/2.0
511	rlat=rlon
512	WRITE(,) "Area : ", area/1.0e6
513	!
514	! Set all the variables defining the polygone of the specified area
515	!
516	neighb_loc(:,:) = -1
517	!
518	! Northern point
519	coord = haversine_radialdis(lon(1,1), lat(1,1), 0.0, rlon)
520	lonpoly(1,2) = coord(1)
521	latpoly(1,2) = coord(2)
522	! Eastern point
523	coord = haversine_radialdis(lon(1,1), lat(1,1), 90.0, rlon)
524	lonpoly(1,4) = coord(1)
525	latpoly(1,4) = coord(2)
526	! Souther point
527	coord = haversine_radialdis(lon(1,1), lat(1,1), 180.0, rlon)
528	lonpoly(1,6) = coord(1)
529	latpoly(1,6) = coord(2)
530	! Souther point
531	coord = haversine_radialdis(lon(1,1), lat(1,1), 270.0, rlon)
532	lonpoly(1,8) = coord(1)
533	latpoly(1,8) = coord(2)
534	!
535	! Doing the corners
536	!
537	! North West
538	lonpoly(1,1) = lonpoly(1,8)
539	latpoly(1,1) = latpoly(1,2)
540	! North East
541	lonpoly(1,3) = lonpoly(1,4)
542	latpoly(1,3) = latpoly(1,2)
543	! South East
544	lonpoly(1,5) = lonpoly(1,4)
545	latpoly(1,5) = latpoly(1,6)
546	! South West
547	lonpoly(1,7) = lonpoly(1,8)
548	latpoly(1,7) = latpoly(1,6)
549	!
550	center(1,1) = lon(1,1)
551	center(1,2) = lat(1,1)
552	!
553	END SUBROUTINE haversine_singlepointploy
554	!!
555	!! =============================================================================================================================
556	!! SUBROUTINE: haversine_polyheadings
557	!!
558	!>\BRIEF Compute the heading out of the polygon for each vertex and mid-point of
559	! segment.
560	!!
561	!! DESCRIPTION: This heading is computed by using the great circle between the centre of the polygon and
562	!! the point on the boundary considered. The direction is the one facing outwards from the polygon.
563	!!
564	!! \n
565	!_ ==============================================================================================================================
566	!!
567	!!
568	SUBROUTINE haversine_polyheadings(nbpt, nbseg, lonpoly, latpoly, center, headings)
569	!
570	! 0.1 Input variables
571	! Size of the gathered domain
572	INTEGER(i_std), INTENT(in) :: nbpt
573	! Number of segments
574	INTEGER(i_std), INTENT(in) :: nbseg
575	!
576	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: lonpoly
577	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: latpoly
578	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: center
579	!
580	! 0.2 Output variables
581	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(out) :: headings
582	!
583	! 0.3 Local variables
584	!
585	INTEGER(i_std) :: i, ns
586	!
587	! We compute for each vertex of our polygon (actual vertices and mid-points) the direction to the
588	! centre. The we add 180. to get the opposite direction.
589	!
590	DO i=1,nbpt
591	DO ns=1,nbseg*2
592	headings(i,ns) = MOD(haversine_heading(lonpoly(i,ns), latpoly(i,ns), center(i,1), center(i,2))+180.0, 360.0)
593	ENDDO
594	ENDDO
595	!
596	END SUBROUTINE haversine_polyheadings
597	!!
598	!! =============================================================================================================================
599	!! SUBROUTINE: haversine_polysort
600	!!
601	!>\BRIEF Sort the polygons so that all points are in clockwise order.
602	!!
603	!! DESCRIPTION:
604	!!
605	!! \n
606	!_ ==============================================================================================================================
607	!!
608	SUBROUTINE haversine_polysort(nbpt, nbseg, lonpoly, latpoly, headings, neighb)
609	!
610	! This subroutine is foreseen for polygones which start with a vertex and then alternate
611	! with mid-points of the segments. The heading at a vertex is the direction (along the
612	! great circle) between the center of the polygone and this vertex. At a segment mid-point
613	! the direction is the normal facing outward.
614	!
615	! 0.1 Input Variables
616	! Size of the gathered domain
617	INTEGER(i_std), INTENT(in) :: nbpt
618	! Number of segments
619	INTEGER(i_std), INTENT(in) :: nbseg
620	!
621	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(inout) :: lonpoly
622	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(inout) :: latpoly
623	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(inout) :: headings
624	INTEGER(i_std), DIMENSION(nbpt,nbseg*2), INTENT(inout) :: neighb
625	!
626	! 0.2 Local variables
627	!
628	INTEGER(i_std) :: i, ic, startindex
629	INTEGER(i_std), DIMENSION(nbseg*2) :: reindex
630	REAL(r_std) :: starthead
631	REAL(r_std), DIMENSION(nbseg*2) :: lon_loc
632	REAL(r_std), DIMENSION(nbseg*2) :: lat_loc
633	REAL(r_std), DIMENSION(nbseg*2) :: head_loc
634	REAL(r_std), DIMENSION(nbseg*2) :: negb_loc
635	!
636	DO i=1,nbpt
637	head_loc(:) = headings(i,:)
638	lon_loc(:) = lonpoly(i,:)
639	lat_loc(:) = latpoly(i,:)
640	negb_loc(:) = neighb(i,:)
641	!
642	! The first vertice of our polygone needs to be the heading closest to
643	! North (0 degree). No difference is done between vertices and mid-points
644	! of segments.
645	!
646	starthead=0.0
647	!
648	CALL haversine_clockwise(nbseg, head_loc, starthead, reindex)
649	!
650	!
651	headings(i,:) = head_loc(reindex(:))
652	lonpoly(i,:) = lon_loc(reindex(:))
653	latpoly(i,:) = lat_loc(reindex(:))
654	neighb(i,:) = negb_loc(reindex(:))
655	!
656	ENDDO
657	!
658	END SUBROUTINE haversine_polysort
659	!!
660	!! =============================================================================================================================
661	!! SUBROUTINE: haversine_polyseglen
662	!!
663	!>\BRIEF Compute the length of all segments of the polygon using the great circle.
664	!!
665	!! DESCRIPTION:
666	!!
667	!! \n
668	!_ ==============================================================================================================================
669	!!
670	SUBROUTINE haversine_polyseglen(nbpt, nbseg, lonpoly, latpoly, seglength)
671	!
672	! Computes the segment length for each of the polygones. These are
673	! polygones with middle points given for each segment. This we need
674	! to take only every other point.
675	!
676	!
677	! 0.1 Input Variables
678	! Size of the gathered domain
679	INTEGER(i_std), INTENT(in) :: nbpt
680	! Number of segments
681	INTEGER(i_std), INTENT(in) :: nbseg
682	!
683	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: lonpoly
684	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: latpoly
685	REAL(r_std), DIMENSION(nbpt,nbseg), INTENT(out) :: seglength
686	!
687	! 0.2 Local variables
688	!
689	INTEGER(i_std) :: i, iv, istart, iend, iseg, ioff
690	REAL(r_std) :: slpm1, slpp1
691	!
692	DO i=1,nbpt
693	iseg = 0
694	!
695	! Find if the first element is a vertex or mid-point of segment.
696	! Get the 2 headings out of the start point. If these headings are larger than 135 degrees
697	! then probably this point is a segment mid-point.
698	!
699	slpm1 = haversine_heading(lonpoly(i,1), latpoly(i,1), lonpoly(i,nbseg2), latpoly(i,nbseg2))
700	slpp1 = haversine_heading(lonpoly(i,1), latpoly(i,1), lonpoly(i,2), latpoly(i,2))
701	!
702	IF ( ABS(MOD(slpp1-slpm1, 360.0)) > 135.0 ) THEN
703	! The polygon starts with a segment mid-point
704	ioff = -1
705	ELSE
706	ioff = 0
707	ENDIF
708	!
709	DO iv=1,nbseg*2,2
710	!
711	istart=MODULO((iv-1)+ioff,nbseg*2)+1
712	iend=MODULO((iv-1)+ioff+2,nbseg*2)+1
713	iseg = iseg + 1
714	!
715	seglength(i,iseg) = haversine_distance(lonpoly(i,istart), latpoly(i,istart), &
716	& lonpoly(i,iend), latpoly(i,iend))
717	ENDDO
718	ENDDO
719	!
720	END SUBROUTINE haversine_polyseglen
721	!!
722	!! =============================================================================================================================
723	!! SUBROUTINE: haversine_laloseglen
724	!!
725	!>\BRIEF Compute the length of all segments of the polygon when on a regular Lat Lon grid.
726	!!
727	!! DESCRIPTION:
728	!!
729	!! \n
730	!_ ==============================================================================================================================
731	!!
732	SUBROUTINE haversine_laloseglen(nbpt, nbseg, lonpoly, latpoly, seglength)
733	!
734	! Computes the segment length for each of the polygones. These are
735	! polygones with middle points given for each segment. This we need
736	! to take only every other point.
737	!
738	!
739	! 0.1 Input Variables
740	! Size of the gathered domain
741	INTEGER(i_std), INTENT(in) :: nbpt
742	! Number of segments
743	INTEGER(i_std), INTENT(in) :: nbseg
744	!
745	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: lonpoly
746	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: latpoly
747	REAL(r_std), DIMENSION(nbpt,nbseg), INTENT(out) :: seglength
748	!
749	! 0.2 Local variables
750	!
751	INTEGER(i_std) :: i, iv, istart, iend, ioff, iseg
752	REAL(r_std) :: coslat, slpm1, slpp1
753	!
754	DO i=1,nbpt
755	iseg = 0
756	!
757	! Find if the first element is a vertex or mid-point of segment.
758	! Get the 2 headings out of the start point. If these headings are larger than 135 degrees
759	! then probably this point is a segment mid-point.
760	!
761	slpm1 = haversine_heading(lonpoly(i,1), latpoly(i,1), lonpoly(i,nbseg2), latpoly(i,nbseg2))
762	slpp1 = haversine_heading(lonpoly(i,1), latpoly(i,1), lonpoly(i,2), latpoly(i,2))
763	!
764	IF ( ABS(MOD(slpp1-slpm1, 360.0)) > 135.0 ) THEN
765	! The polygon starts with a segment mid-point
766	ioff = -1
767	ELSE
768	ioff = 0
769	ENDIF
770	!
771	DO iv=1,nbseg*2,2
772	!
773	istart=MODULO((iv-1)+ioff,nbseg*2)+1
774	iend=MODULO((iv-1)+ioff+2,nbseg*2)+1
775	iseg = iseg + 1
776	!
777	!
778	IF ( ABS(lonpoly(i,istart)-lonpoly(i,iend)) < EPSILON(lonpoly) ) THEN
779	!
780	! Distance along a meridian
781	!
782	seglength(i,iseg) = ABS(latpoly(i,istart) - latpoly(i,iend)) * &
783	pi/180. * R_Earth
784	!
785	ELSE IF ( ABS(latpoly(i,istart)-latpoly(i,iend)) < EPSILON(latpoly) ) THEN
786	!
787	! Distance along a circle of constant latitude
788	!
789	coslat = MAX(COS(latpoly(i,istart) * pi/180.), mincos)
790	seglength(i,iseg) = ABS( lonpoly(i,istart) - lonpoly(i,iend) ) * &
791	pi/180. * R_Earth * coslat
792	!
793	ELSE
794	CALL ipslerr(3, "haversine_laloseglen", "The polygon here does not originate from a regular", &
795	& "latitude longitude grid.","")
796	ENDIF
797	!
798	ENDDO
799	ENDDO
800	!
801	END SUBROUTINE haversine_laloseglen
802	!!
803	!! =============================================================================================================================
804	!! SUBROUTINE: haversine_clockwise
805	!!
806	!>\BRIEF Get the indices which sort a polygon in a clockwise order starting from an
807	! initial vertex given by start.
808	!!
809	!! DESCRIPTION:
810	!!
811	!! \n
812	!_ ==============================================================================================================================
813	!!
814	SUBROUTINE haversine_clockwise(nbseg, heading, start, sortindex)
815	!
816	! Find the order of the polygone vertices which start at "start" and
817	! follow in a clockwise direction.
818	!
819	! 0.1 Input Variables
820	! Number of segments
821	INTEGER(i_std), INTENT(in) :: nbseg
822	REAL(r_std), DIMENSION(nbseg*2), INTENT(in) :: heading
823	REAL(r_std), INTENT(in) :: start
824	!
825	! 0.2 Output variable
826	!
827	INTEGER(i_std), DIMENSION(nbseg*2), INTENT(out) :: sortindex
828	!
829	! 0.3 Local variables
830	!
831	INTEGER(i_std) :: is, js, imin(1)
832	REAL(r_std) :: delang
833	REAL(r_std) :: undef = 9999999999.99999
834	REAL(r_std), DIMENSION(nbseg*2) :: workhead
835	!
836	delang = 360.0/(nbseg*2)
837	!
838	workhead(:) = heading(:)
839	!
840	DO is=1,nbseg*2
841	!
842	! Compute the difference of heading to the next target angle
843	!
844	DO js=1,nbseg*2
845	IF ( workhead(js) < undef ) THEN
846	workhead(js) = MOD(heading(js)-(start+(is-1)*delang)+360.0, 360.0)
847	! Transfer to -180:180 interval
848	IF (workhead(js) > 180.0) workhead(js)=workhead(js)-360.0
849	ENDIF
850	ENDDO
851	!
852	! Locate the vertex closest to that target angle
853	!
854	imin=MINLOC(ABS(workhead))
855	sortindex(is) = imin(1)
856	!
857	! Mask this vertex so that it is skipped in the next iteration
858	!
859	workhead(imin(1)) = undef
860	!
861	ENDDO
862	!
863	END SUBROUTINE haversine_clockwise
864	!!
865	!! =============================================================================================================================
866	!! SUBROUTINE: haversine_polyarea
867	!!
868	!>\BRIEF Computes the area covered by a polygon on the sphere.
869	!!
870	!! DESCRIPTION: Computes the area of each polygon based on Girard's theorem. It
871	!! states that the area of a polygon of great circles is R**2 times
872	!! the sum of the angles between the polygons minus (N-2)*pi where N
873	!! x is number of corners.
874	!!
875	!! \n
876	!_ ==============================================================================================================================
877	!!
878	SUBROUTINE haversine_polyarea(nbpt, nbseg, lonpoly, latpoly, area)
879	!
880	!
881	! 0.1 Input Variables
882	! Size of the gathered domain
883	INTEGER(i_std), INTENT(in) :: nbpt
884	! Number of segments
885	INTEGER(i_std), INTENT(in) :: nbseg
886	!
887	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: lonpoly
888	REAL(r_std), DIMENSION(nbpt,nbseg*2), INTENT(in) :: latpoly
889	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: area
890	!
891	! 0.2 Local variables
892	!
893	INTEGER(i_std) :: i, ia, ib1, ib2, iseg
894	REAL(r_std) :: beta1, beta2
895	REAL(r_std), DIMENSION(nbseg) :: angles
896	!
897	DO i=1,nbpt
898	iseg=0
899	DO ia=1,nbseg*2,2
900	!! Index of next vertex, i.e. ia+1
901	ib1=MOD(ia+1, nbseg*2)+1
902	!! Index of previous vertex, i.e. ia-2
903	ib2=MOD(ia+nbseg2-3, nbseg2)+1
904	iseg=iseg+1
905	!
906	beta1=haversine_dtor(haversine_heading(lonpoly(i,ia), latpoly(i,ia), lonpoly(i,ib1), latpoly(i,ib1)))
907	beta2=haversine_dtor(haversine_heading(lonpoly(i,ia), latpoly(i,ia), lonpoly(i,ib2), latpoly(i,ib2)))
908	!
909	angles(iseg)=acos(cos(-beta1)cos(-beta2) + sin(-beta1)sin(-beta2))
910	!
911	ENDDO
912	area(i) = (sum(angles) - (nbseg-2)pi)R_Earth**2
913	ENDDO
914	!
915	END SUBROUTINE haversine_polyarea
916	!!
917	!! =============================================================================================================================
918	!! SUBROUTINE: haversine_laloarea
919	!!
920	!>\BRIEF Computes the area of a regular latitude longitude box for which we already have the
921	!! the segment length.
922	!!
923	!! DESCRIPTION: Just verify that we have 4 segments.
924	!!
925	!! \n
926	!_ ==============================================================================================================================
927	!!
928	SUBROUTINE haversine_laloarea(nbpt, nbseg, seglen, area)
929	!
930	!
931	! 0.1 Input Variables
932	! Size of the gathered domain
933	INTEGER(i_std), INTENT(in) :: nbpt
934	! Number of segments
935	INTEGER(i_std), INTENT(in) :: nbseg
936	!
937	REAL(r_std), DIMENSION(nbpt,nbseg), INTENT(in) :: seglen
938	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: area
939	!
940	! 0.2 Local variables
941	!
942	INTEGER(i_std) :: i
943	!
944	IF ( nbseg .NE. 4 ) THEN
945	CALL ipslerr(3, "haversine_laloarea", "We need to have 4 segments in the polygons", &
946	& "to use this subroutine","")
947	ENDIF
948	!
949	DO i=1,nbpt
950	area(i) = (seglen(i,1)+seglen(i,3))/2.0*(seglen(i,2)+seglen(i,4))/2.0
951	ENDDO
952	!
953	END SUBROUTINE haversine_laloarea
954	!!
955	!! =============================================================================================================================
956	!! SUBROUTINE: haversine_laloarea
957	!!
958	!>\BRIEF Computes the area in the special case where the projection has given us the size in X and Y of each
959	!! grid box.
960	!!
961	!! DESCRIPTION:
962	!!
963	!! \n
964	!_ ==============================================================================================================================
965	!!
966	SUBROUTINE haversine_xyarea(nbpt, nbseg, ilandtoij, jlandtoij, dx, dy, area)
967	!
968	!
969	! 0.1 Input Variables
970	! Size of the gathered domain
971	INTEGER(i_std), INTENT(in) :: nbpt
972	! Number of segments
973	INTEGER(i_std), INTENT(in) :: nbseg
974	!
975	INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: ilandtoij, jlandtoij
976	REAL(r_std), DIMENSION(:,:), INTENT(in) :: dx, dy
977	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: area
978	!
979	! 0.2 Local variables
980	!
981	INTEGER(i_std) :: il
982	!
983	DO il=1,nbpt
984	area(il) = dx(ilandtoij(il),jlandtoij(il))*dy(ilandtoij(il),jlandtoij(il))
985	ENDDO
986	!
987	END SUBROUTINE haversine_xyarea
988	!!
989	!! =============================================================================================================================
990	!! FUNCTIONS: haversine_heading, haversine_distance, haversine_dtor, haversine_rtod
991	!!
992	!>\BRIEF Various functions to help with the calculations in this module.
993	!!
994	!! DESCRIPTION:
995	!!
996	!! \n
997	!_ ==============================================================================================================================
998	REAL(r_std) FUNCTION haversine_heading(lon_start, lat_start, lon_end, lat_end)
999	!
1000	! The heading is an angle in degrees between 0 and 360.
1001	!
1002	IMPLICIT NONE
1003	REAL(r_std), INTENT(IN) :: lon_start, lat_start, lon_end, lat_end
1004	!
1005	REAL(r_std) :: dlon, lat1, lat2, y, x
1006	!
1007	dlon = haversine_dtor(lon_end-lon_start)
1008	lat1 = haversine_dtor(lat_start)
1009	lat2 = haversine_dtor(lat_end)
1010	y = sin(dlon) * cos(lat2)
1011	x = cos(lat1)* sin(lat2) - sin(lat1) * cos(lat2)* cos(dLon)
1012	haversine_heading = MOD(haversine_rtod(atan2(y,x))+360.0, 360.0)
1013	!
1014	END FUNCTION haversine_heading
1015	!!
1016	!!
1017	!!
1018	REAL(r_std) FUNCTION haversine_diffangle(x,y)
1019	!
1020	! This function compues the smallest angle between 2 headings.
1021	! A positive angle will be in mathematical positive direction (anti-clockwise)
1022	! and a negative one is measured clockwise.
1023	!
1024	IMPLICIT NONE
1025	REAL(r_std), INTENT(IN) :: x, y
1026	!
1027	REAL(r_std) :: rx, ry
1028	!
1029	rx=haversine_dtor(x)
1030	ry=haversine_dtor(y)
1031	haversine_diffangle = haversine_rtod(atan2(sin(ry-rx), cos(ry-rx)))
1032	!
1033	END FUNCTION haversine_diffangle
1034	!!
1035	!!
1036	!!
1037	REAL(r_std) FUNCTION haversine_distance(lon_start, lat_start, lon_end, lat_end)
1038	IMPLICIT NONE
1039	REAL(r_std), INTENT(IN) :: lon_start, lat_start, lon_end, lat_end
1040	!
1041	REAL(r_std) :: dlon, dlat, lat1, lat2, a, c
1042	!
1043	dlat = haversine_dtor(lat_end-lat_start)
1044	dlon = haversine_dtor(lon_end-lon_start)
1045	lat1 = haversine_dtor(lat_start)
1046	lat2 = haversine_dtor(lat_end)
1047	a = sin(dlat/2) * sin(dlat/2) + sin(dlon/2) * sin(dlon/2) * cos(lat1) * cos(lat2)
1048	c = 2 * atan2(sqrt(a), sqrt(1-a))
1049	haversine_distance = c * R_Earth
1050	!
1051	END FUNCTION haversine_distance
1052	!!
1053	!!
1054	!!
1055	FUNCTION haversine_radialdis(lon_start, lat_start, head, dis)
1056	IMPLICIT NONE
1057	REAL(r_std), INTENT(IN) :: lon_start, lat_start, head, dis
1058	REAL(r_std), DIMENSION(2) :: haversine_radialdis
1059	!
1060	REAL(r_std) :: lonout, latout, lat1, lon1, lat, lon, dlon, tc
1061	!
1062	lon1 = haversine_dtor(lon_start)
1063	lat1 = haversine_dtor(lat_start)
1064	tc = haversine_dtor(head)
1065	!
1066	lat =asin(sin(lat1)cos(dis/R_Earth)+cos(lat1)sin(dis/R_Earth)*cos(tc))
1067	dlon = atan2(sin(tc)sin(dis/R_Earth)cos(lat1),cos(dis/R_Earth)-sin(lat1)*sin(lat))
1068	lon = MOD(lon1+dlon+pi, 2*pi)-pi
1069	!
1070	latout = haversine_rtod(lat)
1071	lonout = haversine_rtod(lon)
1072	!
1073	haversine_radialdis(1) = lonout
1074	haversine_radialdis(2) = latout
1075	!
1076	END FUNCTION haversine_radialdis
1077	! --------------------------------------------------------------------
1078	! FUNCTION haversine_dtor():
1079	! This function takes a REAL argument in degree and converts it to
1080	! the equivalent radian.
1081	! --------------------------------------------------------------------
1082	REAL(r_std) FUNCTION haversine_dtor(Degree)
1083	IMPLICIT NONE
1084	REAL(r_std), INTENT(IN) :: Degree
1085	haversine_dtor = Degree * pi/180.0
1086	END FUNCTION haversine_dtor
1087	! --------------------------------------------------------------------
1088	! FUNCTION haversine_rtod():
1089	! This function takes a REAL argument in radian and converts it to
1090	! the equivalent degree.
1091	! --------------------------------------------------------------------
1092	REAL(r_std) FUNCTION haversine_rtod(Radian)
1093	IMPLICIT NONE
1094	REAL(r_std), INTENT(IN) :: Radian
1095
1096	haversine_rtod = Radian * 180.0/pi
1097	END FUNCTION haversine_rtod
1098	!!
1099	!!
1100	!!
1101	END MODULE haversine

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

Download in other formats: