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 trunk/NEMOGCM/TOOLS/SECTIONS_DIADCT/src – NEMO

Context Navigation

source: trunk/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:executable set to ``*
File size: 33.7 KB

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

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

Download in other formats: