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

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

Last change on this file since 2849 was 2849, checked in by cbricaud, 13 years ago

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