Context Navigation

haversine.f90 @ 6892

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_global/haversine.f90 @ 6892

Download in other formats: