New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

compute_sections.f90 in branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src – NEMO

Context Navigation

source: branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90 @ 2877

Last change on this file since 2877 was 2877, checked in by cbricaud, 13 years ago
coding rules
Property svn:executable set to ``*
File size: 35.0 KB

Line
1	MODULE compute_sections
2	!!=====================================================================
3	!! * MODULE diadct *
4	!! Ocean diagnostics: Compute the transport trough a sec.
5	!!
6	!!
7	!! History: 2011: cbricaud Mercator-Ocean
8	!!===============================================================
9	!! * Modules used
10	USE declarations
11	USE sections_tools
12
13	IMPLICIT NONE
14	PRIVATE
15
16	!! * Routine accessibility
17	PUBLIC compsec
18
19	!! * Module variables
20
21	CONTAINS
22
23	SUBROUTINE compsec(jsec,sec,lkdebug)
24	!!---------------------------------------------------------------------
25	!! * ROUTINE compsec *
26	!!
27	!! ** Purpose : Compute the serie of mesh's points that represents the section
28	!! defined by its extremities.
29	!!
30	!! ** Method :
31	!! I. Found which cells of the mesh the section is crossing
32	!! II. Classification of the intersections mesh/section
33	!! -first classification west to east
34	!! -second classification south to north
35	!! III. Find extremities of section in the mesh
36	!! IV. Find the serie of mesh's points that form the section
37	!! ** Input: sec : the section to compute
38	!!
39	!! ** Output:
40	!!---------------------------------------------------------------------
41	!! * Arguments
42	INTEGER,INTENT(IN) :: jsec !number of the section
43	TYPE(SECTION), INTENT(INOUT) :: sec !information about the section
44	LOGICAL ,INTENT(IN) :: lkdebug !debug or not debug this section
45
46	!! * Local variables
47	INTEGER :: &
48	ji,jj , & ! dummy loop argument
49	jseg , & ! loop on segments taht form the section
50	nb_inmesh , & ! number of intersection between section and the mesh
51	nmesh ! number of cells in processor domain
52	INTEGER :: itest , jtest ! dummy integer
53	REAL(wp),SAVE :: zdistmesh ! Taller cell of the mesh in ocean
54	REAL(wp) :: &
55	zdistEst , zdistNorth , zdistWest , zdistSouth , &! temporary scalars
56	zdistEst2 , zdistNorth2, zdistWest2, zdistSouth2, &! temporary scalars
57	zdistEst3 , zdistNorth3, zdistWest3, zdistSouth3, &! temporary scalars
58	zdistFirst , zdistLast , zdistref , &! "
59	zdistante , zdistante2 , zdistnew , zdistnew2 , &! "
60	zdeltai , zdeltaj ! "
61	LOGICAL :: &
62	ll_overlap_sec_left = .FALSE. , ll_overlap_sec_right = .FALSE. ,&! temporary logical
63	ll_date_domain_left = .FALSE. , ll_date_domain_right = .FALSE. ,&! "
64	ll_overlap_sec = .FALSE. , ll_date_domain = .FALSE. ,&! "
65	ll_test = .FALSE. ! "
66	LOGICAL :: lest, lwest, lnorth, lsouth
67	LOGICAL :: l_oldmethod
68	CHARACTER(len=120) :: cltmp
69	TYPE(COORD_SECTION) :: &
70	coord_a , coord_b , coord_d , coord_t , &!temporary point (long/lat)
71	coordFirst, coordLast, coordTemp ! "
72	TYPE(COORD_SECTION), DIMENSION(nb_point_max) :: coordSec !intersections mesh/section
73	TYPE(POINT_SECTION) :: &
74	endingPoint, prevPoint , nextPoint, &!temporary point (I/J)
75	SouthPoint , NorthPoint, EstPoint , WestPoint ! "
76	!!---------------------------------------------------------------------
77	IF( jsec==1 )THEN
78	PRINT*,' '
79	PRINT*,'COMPUTE SECTIONS'
80	PRINT*,'----------------'
81	ENDIF
82
83	!debug
84	CALL write_debug(jsec,'compute section:')
85	CALL write_debug(jsec,'================')
86
87
88	!==================!
89	!0. Initializations!
90	!==================!
91
92	nmesh = jpi*jpj ! number of cells in processor domain
93	nb_inmesh = 0 !initialize number of intersection between section and the mesh
94	zdistEst =0. ; zdistNorth=0. ; zdistWest=0. ; zdistSouth=0. ! temporary scalars
95	zdistFirst=0. ; zdistLast =0. ! "
96	zdistante =0. ; zdistante2=0. ; zdistnew=0. ; zdistnew2=0. ! "
97	zdeltai=0. ; zdeltaj=0.
98	coord_a = COORD_SECTION( 0., 0. ) ; coord_b = COORD_SECTION( 0., 0. )
99	coord_d = COORD_SECTION( 0., 0. ) ; coord_t = COORD_SECTION( 0., 0. )
100	coordFirst = COORD_SECTION( 0., 0. ) ; coordLast = COORD_SECTION( 0., 0. )
101	coordTemp = COORD_SECTION( 0., 0. ) ; coordSec = COORD_SECTION( 0., 0. )
102	endingPoint = POINT_SECTION( -1, -1 )
103	prevPoint = POINT_SECTION( -1, -1 ) ; nextPoint = POINT_SECTION( -1, -1 )
104	SouthPoint = POINT_SECTION( -1, -1 ) ; NorthPoint = POINT_SECTION( -1, -1 )
105	EstPoint = POINT_SECTION( -1, -1 ) ; WestPoint = POINT_SECTION( -1, -1 )
106
107	IF( jsec == 1 )THEN
108	!Found the taller cell of the mesh in ocean (used in compsec)
109	zdistmesh=0.
110	DO jj=2,jpj
111	DO ji=2,jpi
112	IF( zdistmesh .LT. &
113	(e1t(ji,jj)/(cos(gphit(ji,jj))*110000)) ) THEN
114	zdistmesh = e1t(ji,jj)/(cos(gphit(ji,jj))*110000)
115	ENDIF
116	ENDDO
117	ENDDO
118
119	ENDIF
120
121	!debug
122	WRITE(cltmp,'(A10,f10.3)')'zdistmesh = ',zdistmesh
123	CALL write_debug(jsec,cltmp)
124
125	!===========================================================!
126	!I. Found which cells of the mesh the section is crossing !
127	!===========================================================!
128
129	!loop on the mesh
130	DO jj=2,jpj
131	DO ji=2,jpi
132	!-----------------------------------------------------------
133	!For each cell of the mesh, we find the intersection between
134	!the section and the cell as described below:
135	!The cell is defined by 4 points A B C D
136	!and the section by S1 and S2
137	!
138	!Section\ ________
139	! \ B\| \|
140	! \ \| \|
141	! \\|one cell\|
142	! .....\ \|
143	! \|\ \|
144	! \| \ \|
145	! A\|__\_____\|D
146	! ........\
147	!----------------\-------------------------------------------
148
149	!definition of points A,B,D (we don't need C)
150	coord_a = COORD_SECTION( glamf(ji-1,jj-1) ,gphif(ji-1,jj-1) )
151	coord_b = COORD_SECTION( glamf(ji-1,jj) ,gphif(ji-1,jj) )
152	coord_d = COORD_SECTION( glamf(ji,jj-1) ,gphif(ji,jj-1) )
153
154	!size of the cell
155	zdeltai= glamf(ji-1,jj) - glamf(ji-1,jj-1) !zdeltai=[AB]
156	zdeltaj= gphif(ji,jj-1) - gphif(ji-1,jj-1) !zdeltaj=[AD]
157
158	IF( ABS(zdeltai) .LE. 2.*zdistmesh .OR. &
159	ABS(zdeltaj) .LE. 2.*zdistmesh )THEN
160
161	!intersection section/[AB]=?
162	coordTemp = intersec(sec,coord_a,coord_b) !compute intersection
163	coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
164	coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
165
166	IF( coordTemp%lon .NE. -9999 )THEN
167	IF( nb_inmesh+1 .GT. nb_point_max )THEN
168	PRINT*,"WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh
169	ELSE
170	nb_inmesh=nb_inmesh+1
171	coordSec(nb_inmesh) = coordTemp !store the intersection's coordinates
172
173	!We need to know if the section crosses the overlapping band.
174
175	!Fist we look if there is an intersection mesh/section
176	!just on the left of the overlapping band:
177	IF( coordTemp%lon .GT. glamf(1,1)-5 .AND. coordTemp%lon .LT. glamf(1,1) ) &
178	& ll_overlap_sec_left = .TRUE.
179	!And we look if there is an intersection mesh/section
180	!just on the right of the overlapping band:
181	IF( coordTemp%lon .GT. glamf(jpi,1) .AND. coordTemp%lon .LT. glamf(1,1)+5 ) &
182	& ll_overlap_sec_right = .TRUE.
183	ENDIF
184	ENDIF
185
186	!intersection section/[AD]=?
187	coordTemp=intersec(sec,coord_a,coord_d) !compute intersection
188	coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
189	coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
190
191	IF( coordTemp%lon .NE. -9999 )THEN
192	IF( nb_inmesh+1 .GT. nb_point_max )THEN
193	PRINT*, "WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh
194	ELSE
195	nb_inmesh=nb_inmesh+1
196	coordSec(nb_inmesh)=coordTemp
197
198	!We need to know if the section crosses the overlapping band:
199	!same test as above
200	IF( coordTemp%lon .GE. glamf(1,1)-3 .AND. coordTemp%lon .LE. glamf(1,1) ) &
201	& ll_overlap_sec_left = .TRUE.
202	IF( coordTemp%lon .GE. glamf(jpi,1) .AND. coordTemp%lon .LE. glamf(jpi,1)+3) &
203	& ll_overlap_sec_right = .TRUE.
204	ENDIF
205	ENDIF
206
207	ENDIF
208
209	!We need to know if the domain crosses the date line:
210	!Fist, we search a mesh point that is just one the left of date line:
211	IF( glamf(ji-1,jj-1) .GT. 175 .AND. glamf(ji-1,jj-1) .LT. 180 ) &
212	& ll_date_domain_left = .TRUE.
213	!And we search a mesh point that is just one the right of date line:
214	IF( glamf(ji-1,jj-1) .GT. -180 .AND. glamf(ji-1,jj-1) .LT. -175 ) &
215	& ll_date_domain_right = .TRUE.
216
217	ENDDO
218	ENDDO !End of the loop on the mesh
219
220
221	!Crossing section/overlapping band (we need to know it for later):
222	!-----------------------------------------------------------------
223	!If there is one intersection mesh/section just on the left of
224	!the overlapping band (ll_overlap_sec_left = .TRUE.)
225	!AND there is one just the right of the overlapping band
226	!(ll_overlap_sec_right = .TRUE.),
227	!so the section crosses the overlapping band.
228	ll_overlap_sec = ll_overlap_sec_left .AND. ll_overlap_sec_right
229
230	!Crossing of the domain and the date line (we need to know it for later):
231	!------------------------------------------------------------------------
232	!If there is one point of the domain that is just on the left of the date line
233	!(ll_date_domain_left = .TRUE.) AND one point that is just on the right of the
234	!date line (ll_date_domain_right = .TRUE. )
235	!So the domain crosses the date line:
236	ll_date_domain = ll_date_domain_left .AND. ll_date_domain_right
237
238	!=====================================================!
239	! II. Classification of the intersections mesh/section!
240	!=====================================================!
241
242	! -first classification west to east
243	! -second classification south to north
244	!CAUTION: routine qcksrt doesn't work in the same way if the section
245	!and the domain crosse the date line (sec%ll_date_line=T and ll_date_domain=T)
246
247	IF( sec%ll_date_line .AND. ll_date_domain )THEN
248
249	!we add 360° to negative longitudes to have a good classification
250	DO jseg=1,nb_inmesh
251	IF( coordSec(jseg)%lon .LT. 0 ) coordSec(jseg)%lon=coordSec(jseg)%lon+360.
252	ENDDO
253	IF( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon ) THEN
254	CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh)
255	ELSE
256	CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh)
257	ENDIF
258	DO jseg=1,nb_inmesh
259	IF( coordSec(jseg)%lon .GT. 180 ) coordSec(jseg)%lon=coordSec(jseg)%lon-360.
260	ENDDO
261
262	ELSE
263
264	IF( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon )THEN
265	CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh)
266	ELSE
267	CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh)
268	ENDIF
269
270	ENDIF
271
272	!debug
273	WRITE(cltmp,'(A20,i3.3)')'number intersections = ',nb_inmesh ; CALL write_debug(jsec,cltmp)
274	CALL write_debug(jsec,'List of intersections between grid and section: ')
275	DO jseg=1,nb_inmesh
276	WRITE(cltmp,'(i4.4,1X,2(f8.3,1X) )')jseg,coordSec(jseg) ; CALL write_debug(jsec,cltmp)
277	ENDDO
278
279	!=====================================================!
280	! III. Find extremities of section in the mesh !
281	!=====================================================!
282	!we can find section's extremities in the mesh only if
283	!there is intersection between section and mesh (nb_inmesh .ne. 0)
284
285	IF( nb_inmesh .ne. 0 )THEN
286	coordFirst = coordSec(1)
287	coordLast = coordSec(nb_inmesh)
288	sec%nb_point = nb_inmesh
289	sec%listPoint(1) = POINT_SECTION(-1,-1)
290	zdistante = 1000.
291	zdistante2 = 1000.
292
293	!First, we find the point of the mesh that is the closest
294	!to the first intersection section/mesh (=coordFirst=coordSec(1)):
295	!this point will be called sec%listPoint(1).
296	!Then, we find the point of the mesh that is the closest
297	!to the last intersection section/mesh (coordLast=coordSec(nb_inmesh))
298	!this point will be called endingPoint.
299
300	DO jj=1,jpj
301	DO ji=1,jpi
302	coord_t=COORD_SECTION(glamf(ji,jj),gphif(ji,jj))
303	!cbr zdistFirst = distance(coord_t,coordFirst)
304	zdistFirst = distance2(coord_t,coordFirst)
305	!cbr zdistLast = distance(coord_t,coordLast)
306	zdistLast = distance2(coord_t,coordLast)
307	IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN
308	sec%listPoint(1) = POINT_SECTION(ji,jj)
309	zdistante=zdistFirst
310	ENDIF
311	IF( zdistLast .LT. zdistmesh .AND. zdistLast .LT. zdistante2 )THEN
312	endingPoint = POINT_SECTION(ji,jj)
313	zdistante2=zdistLast
314	ENDIF
315	ENDDO
316	ENDDO
317
318	IF( sec%listPoint(1)%I == endingPoint%I .AND. sec%listPoint(1)%J == endingPoint%J )THEN
319	sec%listPoint(1) = POINT_SECTION(-1,-1)
320	endingPoint = POINT_SECTION(-1,-1)
321	coordFirst = coordSec(1)
322	coordLast = coordSec(2)
323	sec%nb_point = 0
324	ENDIF
325
326	ELSE
327	!If there is no intersection section/mesh
328	sec%listPoint(1) = POINT_SECTION(-1,-1)
329	endingPoint = POINT_SECTION(-1,-1)
330	coordFirst = coordSec(1)
331	coordLast = coordSec(2)
332	sec%nb_point = 0
333	ENDIF
334
335	!debug
336	CALL write_debug(jsec,"extremities of section in the grid : ")
337	ji=sec%listPoint(1)%I ; jj=sec%listPoint(1)%J
338	IF( sec%nb_point .ne. 0 )THEN
339	ji=sec%listPoint(1)%I ; jj=sec%listPoint(1)%J
340	WRITE(cltmp,'(A15,X,i4.4,X,i4.4,X,f8.3,X,f8.3)')'First point: ',sec%listPoint(1),glamf(ji,jj),gphif(ji,jj)
341	CALL write_debug(jsec,cltmp)
342	ji=endingPoint%I ; jj=endingPoint%J
343	WRITE(cltmp,'(A15,X,i4.4,X,i4.4,X,f8.3,X,f8.3)')'Last point: ',endingPoint,glamf(ji,jj),gphif(ji,jj)
344	CALL write_debug(jsec,cltmp)
345	!
346	coord_a=pointToCoordF(sec%listPoint(1)) ; coord_b=pointToCoordF(endingPoint)
347	ll_test = .FALSE.
348	IF ( ll_date_domain .AND. ABS( coord_a%lon - coord_b%lon ).GT. 180) ll_test= .TRUE.
349	zdistante=distance2(coord_a,coord_b ,ll_test )
350	WRITE(cltmp,'(A20,f10.3)' )'distance between IJ-extremities : ',zdistante
351	CALL write_debug(jsec,cltmp)
352	!
353	CALL write_debug(jsec,"Initial extremities : ")
354	WRITE(cltmp,'( 2(f9.3),A3,2(f9.3) )')coordFirst,'---',coordLast
355	CALL write_debug(jsec,cltmp)
356	ll_test = .FALSE.
357	IF( ll_date_domain .AND. ABS(coordFirst%lon - coordLast%lon).GT. 180)ll_test= .TRUE.
358	zdistante=distance2(coordFirst,coordLast,ll_test)
359	WRITE(cltmp,'(A30,f10.3)')' distance between initial extremities : ',zdistante
360	CALL write_debug(jsec,cltmp)
361	CALL write_debug(jsec," ")
362	ELSE
363	WRITE(cltmp,'(A50)' )"no intersection between section and mesh"
364	ENDIF
365
366	!==========================================================!
367	! IV. Find the serie of mesh's points that form the section!
368	!==========================================================!
369	CALL write_debug(jsec,"Find the serie of mesh's points that form the section")
370
371	IF( sec%nb_point .ne. 0 )THEN
372	!IF( nb_inmesh .NE. 0 )THEN
373
374	!The serie of mesh's points that form the section will 'link'
375	!sec%listPoint(1) to endingPoint: it will be stored in
376	!sec%listPoint(jseg)
377	!
378	!We take place on the fist point (sec%listPoint(1))
379	! a. We find the 4 adjacent points (North, South, East, West)
380	! b. Compute distance between current point and endingPoint
381	! c. Compute distance between the 4 adjacent points and endingPoint
382	! d. Select the points which are closer to end-point than current point
383	! e.1 If at least one point is selected, select the point which is closest to original section among selected points
384	! e.2 If no point is selected, select the point which is the closest to end-point
385	! f. save next point and direction of velocity.
386	! g. Save nextPoint and go to nextPoint
387	!
388	!We get out of this loop if:
389	! - we are on endingPoint
390	! - the number of points (jseg) that link sec%listPoint(1) to endingPoint is
391	! twice greater than number of section/mesh intersection (nb_inmesh):
392	! it could be possible if thr algorithm can't link endingPoint (bug).
393
394	!initialize distnew value (with distance between section's extremities)
395	!cbr zdistnew = distance(coordFirst,coordLast,sec%ll_date_line)
396	zdistnew = distance2(coordFirst,coordLast,sec%ll_date_line)
397	prevPoint = POINT_SECTION(0,0)
398	jseg = 1
399
400	DO WHILE ( ( sec%listPoint(jseg)%I .NE. endingPoint%I &
401	.OR. sec%listPoint(jseg)%J .NE. endingPoint%J ) &
402	.AND. jseg .LT. 500 .AND. sec%listPoint(jseg)%I .GT. 0 )
403
404	! a. find the 4 adjacent points (North, South, East, West)
405	!---------------------------------------------------------
406	SouthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J-1)
407	NorthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
408	WestPoint = POINT_SECTION(sec%listPoint(jseg)%I-1,sec%listPoint(jseg)%J)
409	EstPoint = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
410
411	!debug
412	CALL write_debug(jsec,"---------------")
413	WRITE(cltmp,100)'Current points: ',sec%listPoint(jseg), &
414	glamf(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J), &
415	gphif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)
416	CALL write_debug(jsec,cltmp)
417	CALL write_debug(jsec,"E/W/N/S points")
418	WRITE(cltmp,102)EstPoint,WestPoint,NorthPoint,SouthPoint
419	CALL write_debug(jsec,cltmp)
420	itest=MIN(MAX(EstPoint%I,0),jpi+1) ; jtest=MIN(MAX(EstPoint%J,0),jpj+1)
421	IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
422	WRITE(cltmp,101)'Est',glamf(itest,jtest),gphif(itest,jtest)
423	CALL write_debug(jsec,cltmp)
424	ELSE
425	CALL write_debug(jsec,"Est point out of domain")
426	ENDIF
427	!
428	itest=MIN(MAX(WestPoint%I,0),jpi+1) ; jtest=MIN(MAX(WestPoint%J,0),jpj+1)
429	IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
430	WRITE(cltmp,101)'West',glamf(itest,jtest),gphif(itest,jtest)
431	CALL write_debug(jsec,cltmp)
432	ELSE
433	CALL write_debug(jsec,"West point out of domain")
434	ENDIF
435	!
436	itest=MIN(MAX(NorthPoint%I,0),jpi+1) ; jtest=MIN(MAX(NorthPoint%J,0),jpj+1)
437	IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
438	WRITE(cltmp,101)'North',glamf(itest,jtest),gphif(itest,jtest)
439	CALL write_debug(jsec,cltmp)
440	ELSE
441	CALL write_debug(jsec,"North point out of domain")
442	ENDIF
443	!
444	itest=MIN(MAX(SouthPoint%I,0),jpi+1) ; jtest=MIN(MAX(SouthPoint%J,0),jpj+1)
445	IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
446	WRITE(cltmp,101)'South',glamf(itest,jtest),gphif(itest,jtest)
447	CALL write_debug(jsec,cltmp)
448	ELSE
449	CALL write_debug(jsec,"South point out of domain")
450	ENDIF
451	!
452	!
453	100 FORMAT ( A15,2(i4.4," "),2(f7.3," ") )
454	101 FORMAT ( A6,2(f7.3," "))
455	102 FORMAT ( "E ",i4.4,' ',i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 )
456
457
458	!we are on end-point
459	!--------------------
460	IF( SouthPoint%I==endingPoint%I .AND. SouthPoint%J==endingPoint%J )THEN
461	jseg = jseg+1 ; sec%listPoint(jseg) = SouthPoint
462	ELSE IF( NorthPoint%I==endingPoint%I .AND. NorthPoint%J==endingPoint%J )THEN
463	jseg = jseg+1 ; sec%listPoint(jseg) = NorthPoint
464	ELSE IF( WestPoint%I==endingPoint%I .AND. WestPoint%J==endingPoint%J )THEN
465	jseg = jseg+1 ; sec%listPoint(jseg) = WestPoint
466	ELSE IF( EstPoint%I==endingPoint%I .AND. EstPoint%J==endingPoint%J )THEN
467	jseg = jseg+1 ; sec%listPoint(jseg) = EstPoint
468
469	ELSE
470	!we are NOT on end-point
471	!------------------------
472
473	! b. distance between current point and endingPoint
474	!--------------------------------------------------
475	zdistante = zdistnew
476
477	! c. compute distance between the 4 adjacent points and endingPoint
478	!------------------------------------------------------------------
479	! BE CARREFUL! When the domain crosses the date line (ll_date_domain):
480	! When we will compute distances between W/E/S/N points and endingPoint,
481	! we have to check if theses 4 lines crosses the date line
482	! (test: ABS(coordTemp%lon - coordLast%lon).GT. 180)
483	! If it's true,we have to add 360° to coordLast%long to compute the
484	! good distance through the date line and not through the complementary
485	! in the mesh.
486
487	! c.1 compute distWest: distance between west point and endingPoint
488	!----------------------
489	zdistWest2 = 99999999.9 ; zdistWest3 = 99999999.9
490	IF( sec%listPoint(jseg)%I .NE. 1 )THEN
491	!When we are on the west side of the mesh (i=1),we can't go to the west.
492	coordTemp = pointToCoordF(WestPoint)
493	ll_test = .FALSE.
494	IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test = .TRUE.
495	zdistWest2 = distance2(pointToCoordF(WestPoint) ,coordLast ,ll_test)
496	ENDIF
497
498	! c.2 compute distEst: distance between east point and endingPoint
499	!---------------------
500	zdistEst2 = 99999999.9 ; zdistEst3 = 99999999.9
501	IF( sec%listPoint(jseg)%I .EQ. jpi )THEN
502	!We test if the current point is on the east side of the mesh
503	! The method is done such as we go toward east to link
504	! sec%listPoint(1) to endingPoint.
505	! So, if the section crosses the overlapping band (ll_overlap_sec=T),
506	! we won't have to stop if the current point is on the EAST side of the mesh:
507	! we have to follow the construction of the section on the
508	! WEST side of the mesh
509	IF( ll_overlap_sec )THEN
510	!section crosses the overlapping band
511	!So EstPoint is on the west side of the mesh
512	EstPoint = POINT_SECTION(3,sec%listPoint(jseg)%J)
513	zdistEst2= distance2(pointToCoordF(EstPoint) ,coordLast ,.FALSE.)
514	ENDIF
515	ELSE
516	coordTemp = pointToCoordF(EstPoint)
517	ll_test = .FALSE.
518	IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
519	zdistEst2 = distance2(pointToCoordF(EstPoint) ,coordLast ,ll_test )
520	ENDIF
521
522	! c.3 compute distSouth: distance between south point and endingPoint
523	!-----------------------
524	zdistSouth2 = 99999999.9 ; zdistSouth3 = 99999999.9
525	IF( sec%listPoint(jseg)%J .NE. 1 )THEN
526	!When we are on the south side of the mesh (j=1),we can't go to the south.
527	coordTemp=pointToCoordF(SouthPoint)
528	ll_test = .FALSE.
529	IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
530	zdistSouth2 = distance2(pointToCoordF(SouthPoint),coordlast ,ll_test )
531	ENDIF
532
533	! c.4 compute distNorth: distance between north and endingPoint
534	!-----------------------
535	zdistNorth2 = 99999999.9 ; zdistNorth3 = 99999999.9
536	IF( sec%listPoint(jseg)%J .NE. jpj )THEN
537	!When we are on the north side of the mesh (j=jpj),we can't go to the south.
538	coordTemp=pointToCoordF(NorthPoint)
539	ll_test = .FALSE.
540	IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
541	zdistNorth2 = distance2(pointToCoordF(NorthPoint),coordlast ,ll_test )
542	ENDIF
543
544	! d. select the points which are closer to end-point than current point
545	!----------------------------------------------------------------------
546	zdistref=distance2(pointToCoordF(sec%listPoint(jseg)),coordlast ,ll_test )
547	WRITE(cltmp,'( A56,f10.3 )' )'distance between actual point and last point: zdistref = ',zdistref
548	CALL write_debug(jsec,cltmp)
549	lest = .FALSE. ; IF( zdistEst2 .LE. zdistref ) lest = .TRUE.
550	lwest = .FALSE. ; IF( zdistwest2 .LE. zdistref ) lwest = .TRUE.
551	lnorth = .FALSE. ; IF( zdistnorth2 .LE. zdistref ) lnorth= .TRUE.
552	lsouth = .FALSE. ; IF( zdistsouth2 .LE. zdistref ) lsouth= .TRUE.
553
554	!debug
555	IF( .NOT. lest )CALL write_debug(jsec,'Est point eliminated')
556	IF( .NOT. lwest )CALL write_debug(jsec,'West point eliminated')
557	IF( .NOT. lnorth )CALL write_debug(jsec,'North point eliminated')
558	IF( .NOT. lsouth )CALL write_debug(jsec,'South point eliminated')
559
560	l_oldmethod = .FALSE.
561
562	IF( ( COUNT((/lest/))+COUNT((/lwest/))+COUNT((/lnorth/))+COUNT((/lsouth/)) ) .NE. 0 )THEN
563
564	! e.1 If at least one point is selected, select the point
565	! which is the closest to original section among selected points
566	!-------------------------------------------------------------------
567
568	zdistWest3 = 9999999.9
569	IF( lwest )zdistWest3 = &
570	distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(WestPoint) ,pointToCoordF(endingPoint) ,lkdebug )
571	zdistEst3 = 9999999.9
572	IF( lest )zdistEst3 = &
573	distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(EstPoint) ,pointToCoordF(endingPoint) ,lkdebug )
574	zdistSouth3 = 9999999.9
575	IF( lsouth )zdistSouth3 = &
576	distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(SouthPoint),pointToCoordF(endingPoint) ,lkdebug )
577	zdistNorth3 = 9999999.9
578	IF( lnorth )zdistNorth3 = &
579	distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(NorthPoint),pointToCoordF(endingPoint) ,lkdebug )
580
581	zdistEst3 = zdistEst3 + (1-COUNT((/lest/)) )*9999999.9
582	zdistWest3 = zdistWest3 + (1-COUNT((/lwest/)) )*9999999.9
583	zdistNorth3 = zdistNorth3 + (1-COUNT((/lnorth/)))*9999999.9
584	zdistSouth3 = zdistSouth3 + (1-COUNT((/lsouth/)))*9999999.9
585
586	zdistnew = MIN(zdistEst3,zdistWest3,zdistnorth3,zdistSouth3)
587
588	ELSE
589
590	! e.2 If no point is selected, select the point which is the closest to end-point
591	!--------------------------------------------------------------------------------
592	l_oldmethod = .TRUE.
593
594	!debug
595	WRITE(cltmp,'(A30,i3.3)' )'SEARCH NEW POINT WITH OLD METHOD: ',jsec
596	CALL write_debug(jsec,cltmp)
597
598	! on passe à l'ancienne methode: le point parmies les 4 pts (NSWE) qui se rapproche
599	! le + du dernier pt
600	!-----------------------------
601	!be carreful! we can't go backward.
602
603	zdistNorth = zdistNorth2 ; zdistSouth = zdistSouth2
604	zdistEst = zdistEst2 ; zdistWest = zdistWest2
605
606	IF( prevPoint%I .EQ. NorthPoint%I .AND. prevPoint%J .EQ. NorthPoint%J) THEN
607	zdistnew = MIN(zdistEst,zdistWest,zdistSouth)
608	ELSE IF(prevPoint%I .EQ. SouthPoint%I .AND. prevPoint%J .EQ. SouthPoint%J) THEN
609	zdistnew = MIN(zdistEst,zdistWest,zdistNorth)
610	ELSE IF(prevPoint%I .EQ. WestPoint%I .AND. prevPoint%J .EQ. WestPoint%J ) THEN
611	zdistnew = MIN(zdistEst,zdistNorth,zdistSouth)
612	ELSE IF(prevPoint%I .EQ. EstPoint%I .AND. prevPoint%J .EQ. EstPoint%J ) THEN
613	zdistnew = MIN(zdistWest,zdistNorth,zdistSouth)
614	ELSE
615	zdistnew = MIN(zdistEst,zdistWest,zdistNorth,zdistSouth)
616	ENDIF
617
618	ENDIF
619
620	!debug
621	WRITE(cltmp,'(A11, f8.3)')'zdistref = ',zdistref
622	CALL write_debug(jsec,cltmp)
623	WRITE(cltmp, 103 )'distance2 :',zdistEst2,zdistWest2,zdistNorth2,zdistSouth2
624	CALL write_debug(jsec,cltmp)
625	WRITE(cltmp, 103 )'distance3 :',zdistEst3,zdistWest3,zdistNorth3,zdistSouth3
626	CALL write_debug(jsec,cltmp)
627	WRITE(cltmp,'(A11, f8.3)')"zdistnew = ",zdistnew
628	CALL write_debug(jsec,cltmp)
629
630	103 FORMAT (A12,"E",f12.3," W",f12.3," N",f12.3," S",f12.3)
631
632	!f. save next point and direction of velocity.
633	!---------------------------------------------
634	!nextPoint will be the one which is the closest to endingPoint.
635	!sec%direction will be direction between current point and nextPoint:
636	!It will be used to compute velocity through the segment [CurrentPoint,nextPoint}:
637	!
638	!A:Current Point NorthPoint(I,J+1) Nextpoint=NorthPoint(I,J+1) => sec%direction=3
639	! \| Nextpoint=SouthPoint(I,J-1) => sec%direction=2
640	! \| Nextpoint=WestPoint(I-1,J) => sec%direction=0
641	! \|==>V(I,J+1) Nextpoint=EastPoint(I+1,J) => sec%direction=1
642	! U(I,J) \| U(I+1,J)
643	! ^ \| ^
644	! West_____\|______A_______\|_____EstPoint
645	! Point \|(I,J) (I+1,J)
646	! (I-1,J) \|
647	! \|==>V(I,J)
648	! \|
649	! SoutPoint(I,J-1)
650	IF( l_oldmethod )THEN
651	IF( zdistnew == zdistWest ) THEN
652	sec%direction(jseg)=0 ; nextPoint = WestPoint
653	ELSE IF( zdistnew == zdistEst ) THEN
654	sec%direction(jseg)=1 ; nextPoint = EstPoint
655	ELSE IF( zdistnew == zdistSouth )THEN
656	sec%direction(jseg)=2 ; nextPoint = SouthPoint
657	ELSE IF( zdistnew == zdistNorth )THEN
658	sec%direction(jseg)=3 ; nextPoint= NorthPoint
659	ENDIF
660	ELSE
661	IF( zdistnew == zdistWest3 ) THEN
662	sec%direction(jseg)=0 ; nextPoint = WestPoint
663	ELSE IF( zdistnew == zdistEst3 ) THEN
664	sec%direction(jseg)=1 ; nextPoint = EstPoint
665	ELSE IF( zdistnew == zdistSouth3 )THEN
666	sec%direction(jseg)=2 ; nextPoint = SouthPoint
667	ELSE IF( zdistnew == zdistNorth3 )THEN
668	sec%direction(jseg)=3 ; nextPoint= NorthPoint
669	ENDIF
670	ENDIF
671
672	WRITE(cltmp,'(A11, 2(i4.4,1X) )')'nextPoint = ', nextPoint
673	CALL write_debug(jsec,cltmp)
674
675	!f. Save nextPoint and go to nextPoint
676	!-------------------------------------
677	prevPoint = sec%listPoint(jseg)
678	jseg = jseg+1 !increment of number of segments that form the section
679	sec%listPoint(jseg) = nextPoint !Save next point
680
681	ENDIF ! southP/northP/WestP/EstP == endingpoint ?
682
683	ENDDO !End of loop on jseg
684	sec%nb_point = jseg !Save the number of segments that form the section
685
686
687	ELSE ! nb_inmesh == 0
688	DO jseg=1,nb_point_max
689	sec%listPoint(:)=POINT_SECTION(0,0)
690	ENDDO
691	sec%direction(:)=0.
692	sec%nb_point = 0
693	ENDIF
694
695	!debug
696	IF( sec%nb_point .ne. 0 )THEN
697	CALL write_debug(jsec,"-------------------------------------")
698	CALL write_debug(jsec,"list of points in the grid : ")
699	DO jseg=1,sec%nb_point
700	ji=sec%listPoint(jseg)%I ; jj=sec%listPoint(jseg)%J
701	WRITE(cltmp, '(i4.4,X,i4.4,X,i4.4,X,f8.3,X,f8.3)' )jseg,ji,jj,glamf(ji,jj),gphif(ji,jj)
702	CALL write_debug(jsec,cltmp)
703	ENDDO
704
705	!test if we are one end-point
706	IF( sec%listPoint(sec%nb_point)%I .NE. endingPoint%J .AND. sec%listPoint(sec%nb_point)%J .NE. endingPoint%J )THEN
707	PRINT*,TRIM(sec%name)," NOT ARRIVED TO endingPoint FOR jsec = ",jsec
708	ENDIF
709	ENDIF
710
711	!now compute new slopeSection with ij-coordinates of first and last point
712	sec%slopeSection = 0 ! default value
713	IF( sec%nb_point .ne. 0 )THEN
714	IF ( sec%listPoint(sec%nb_point)%I .NE. sec%listPoint(1)%I ) THEN
715	sec%slopeSection = ( sec%listPoint(sec%nb_point)%J - sec%listPoint(1)%J ) / &
716	( sec%listPoint(sec%nb_point)%I - sec%listPoint(1)%I )
717	ELSE
718	sec%slopeSection = 10000._wp
719	ENDIF
720	ENDIF
721
722	END SUBROUTINE compsec
723
724	END MODULE compute_sections

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

Download in other formats: