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 @ 2860

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

correction of array bounds
  • Property svn:executable set to *
File size: 34.8 KB
Line 
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     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              zdistFirst = distance(coord_t,coordFirst)
304              zdistLast = distance(coord_t,coordLast)
305              IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN
306                 sec%listPoint(1) = POINT_SECTION(ji,jj)
307                 zdistante=zdistFirst
308              ENDIF
309              IF( zdistLast .LT. zdistmesh .AND. zdistLast .LT. zdistante2 )THEN
310                 endingPoint = POINT_SECTION(ji,jj)
311                 zdistante2=zdistLast
312              ENDIF
313           ENDDO
314        ENDDO
315
316        IF( sec%listPoint(1)%I == endingPoint%I .AND. sec%listPoint(1)%J == endingPoint%J )THEN
317           sec%listPoint(1) = POINT_SECTION(-1,-1)
318           endingPoint      = POINT_SECTION(-1,-1)
319           coordFirst       = coordSec(1)
320           coordLast        = coordSec(2)
321           sec%nb_point     = 0
322        ENDIF
323
324     ELSE
325        !If there is no intersection section/mesh
326        sec%listPoint(1) = POINT_SECTION(-1,-1)
327        endingPoint      = POINT_SECTION(-1,-1)
328        coordFirst       = coordSec(1)
329        coordLast        = coordSec(2)
330        sec%nb_point     = 0
331     ENDIF 
332
333     !debug
334     CALL write_debug(jsec,"extremities of section in the grid : ")
335     ji=sec%listPoint(1)%I ; jj=sec%listPoint(1)%J
336     IF( sec%nb_point .ne. 0 )THEN
337        ji=sec%listPoint(1)%I ; jj=sec%listPoint(1)%J
338        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) 
339        CALL write_debug(jsec,cltmp)
340        ji=endingPoint%I ; jj=endingPoint%J
341        WRITE(cltmp,'(A15,X,i4.4,X,i4.4,X,f8.3,X,f8.3)')'Last  point: ',endingPoint,glamf(ji,jj),gphif(ji,jj)
342        CALL write_debug(jsec,cltmp)
343        !
344        coord_a=pointToCoordF(sec%listPoint(1)) ; coord_b=pointToCoordF(endingPoint)
345        ll_test = .FALSE.
346        IF ( ll_date_domain .AND. ABS( coord_a%lon - coord_b%lon ).GT. 180) ll_test= .TRUE.
347        zdistante=distance2(coord_a,coord_b ,ll_test )
348        WRITE(cltmp,'(A20,f10.3)' )'distance between IJ-extremities : ',zdistante
349        CALL write_debug(jsec,cltmp)
350        !
351        CALL write_debug(jsec,"Initial extremities : ") 
352        WRITE(cltmp,'( 2(f9.3),A3,2(f9.3) )')coordFirst,'---',coordLast
353        CALL write_debug(jsec,cltmp)
354        ll_test = .FALSE.
355        IF( ll_date_domain .AND. ABS(coordFirst%lon - coordLast%lon).GT. 180)ll_test= .TRUE.
356        zdistante=distance2(coordFirst,coordLast,ll_test)
357        WRITE(cltmp,'(A30,f10.3)')' distance between initial extremities : ',zdistante
358        CALL write_debug(jsec,cltmp)
359        CALL write_debug(jsec,"                  ")
360     ELSE
361        WRITE(cltmp,'(A50)' )"no intersection between section and mesh"
362     ENDIF
363
364     !==========================================================!
365     ! IV. Find the serie of mesh's points that form the section!
366     !==========================================================!
367     CALL write_debug(jsec,"Find the serie of mesh's points that form the section")
368
369     IF( sec%nb_point .ne. 0 )THEN
370     !IF( nb_inmesh .NE. 0 )THEN
371
372        !The serie of mesh's points that form the section will 'link'
373        !sec%listPoint(1) to endingPoint: it will be stored in
374        !sec%listPoint(jseg)
375        !
376        !We take place on the fist point (sec%listPoint(1))
377        ! a.  We find the 4 adjacent points (North, South, East, West)
378        ! b.  Compute distance between current point and endingPoint
379        ! c.  Compute distance between the 4 adjacent points and endingPoint
380        ! d.  Select the points which are closer to end-point than current point
381        ! e.1 If at least one point is selected, select the point which is closest to original section among selected points
382        ! e.2 If no point is selected, select the point which is the closest to end-point
383        ! f. save next point and direction of velocity.
384        ! g. Save nextPoint and go to nextPoint
385        !
386        !We get out of this loop if:
387        !    - we are on endingPoint
388        !    - the number of points (jseg) that link sec%listPoint(1) to endingPoint is
389        !      twice greater than number of section/mesh intersection (nb_inmesh):
390        !      it could be possible if thr algorithm can't link endingPoint (bug).
391
392        !initialize distnew value (with distance between section's extremities)
393        zdistnew  = distance(coordFirst,coordLast,sec%ll_date_line) 
394        prevPoint  = POINT_SECTION(0,0)
395        jseg       = 1
396
397        DO WHILE ( (  sec%listPoint(jseg)%I .NE.  endingPoint%I    &
398                 .OR. sec%listPoint(jseg)%J .NE. endingPoint%J   ) &
399                 .AND. jseg .LT. 500 .AND. sec%listPoint(jseg)%I .GT. 0  )         
400   
401           ! a. find the 4 adjacent points (North, South, East, West)
402           !---------------------------------------------------------
403           SouthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J-1)
404           NorthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
405           WestPoint  = POINT_SECTION(sec%listPoint(jseg)%I-1,sec%listPoint(jseg)%J)
406           EstPoint   = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
407
408           !debug
409           CALL write_debug(jsec,"---------------")
410           WRITE(cltmp,100)'Current points: ',sec%listPoint(jseg), & 
411                         glamf(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J), &
412                         gphif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)
413           CALL write_debug(jsec,cltmp)
414           CALL write_debug(jsec,"E/W/N/S points")
415           WRITE(cltmp,102)EstPoint,WestPoint,NorthPoint,SouthPoint
416           CALL write_debug(jsec,cltmp)
417           itest=MIN(MAX(EstPoint%I,0),jpi+1) ; jtest=MIN(MAX(EstPoint%J,0),jpj+1)
418           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
419              WRITE(cltmp,101)'Est',glamf(itest,jtest),gphif(itest,jtest)
420              CALL write_debug(jsec,cltmp)
421           ELSE
422              CALL write_debug(jsec,"Est point out of domain")
423           ENDIF
424           !
425           itest=MIN(MAX(WestPoint%I,0),jpi+1) ; jtest=MIN(MAX(WestPoint%J,0),jpj+1)
426           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
427              WRITE(cltmp,101)'West',glamf(itest,jtest),gphif(itest,jtest)
428              CALL write_debug(jsec,cltmp)
429           ELSE
430              CALL write_debug(jsec,"West point out of domain")
431           ENDIF
432           !
433           itest=MIN(MAX(NorthPoint%I,0),jpi+1) ; jtest=MIN(MAX(NorthPoint%J,0),jpj+1)
434           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
435              WRITE(cltmp,101)'North',glamf(itest,jtest),gphif(itest,jtest)
436              CALL write_debug(jsec,cltmp)
437           ELSE
438              CALL write_debug(jsec,"North point out of domain")
439           ENDIF
440           !
441           itest=MIN(MAX(SouthPoint%I,0),jpi+1) ; jtest=MIN(MAX(SouthPoint%J,0),jpj+1)
442           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN
443              WRITE(cltmp,101)'South',glamf(itest,jtest),gphif(itest,jtest)
444              CALL write_debug(jsec,cltmp)
445           ELSE
446              CALL write_debug(jsec,"South point out of domain")
447           ENDIF
448           !
449           !
450100 FORMAT ( A15,2(i4.4," "),2(f7.3," ") )
451101 FORMAT ( A6,2(f7.3," "))
452102 FORMAT ( "E ",i4.4,' ',i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 )
453
454               
455           !we  are on end-point
456           !--------------------
457           IF(      SouthPoint%I==endingPoint%I .AND. SouthPoint%J==endingPoint%J )THEN
458               jseg = jseg+1 ; sec%listPoint(jseg) = SouthPoint
459           ELSE IF( NorthPoint%I==endingPoint%I .AND. NorthPoint%J==endingPoint%J )THEN
460               jseg = jseg+1 ; sec%listPoint(jseg) = NorthPoint
461           ELSE IF(  WestPoint%I==endingPoint%I .AND.  WestPoint%J==endingPoint%J )THEN
462               jseg = jseg+1 ; sec%listPoint(jseg) = WestPoint
463           ELSE IF(   EstPoint%I==endingPoint%I .AND.   EstPoint%J==endingPoint%J )THEN
464               jseg = jseg+1 ; sec%listPoint(jseg) = EstPoint
465
466           ELSE
467           !we  are NOT on end-point
468           !------------------------
469
470              ! b. distance between current point and endingPoint
471              !--------------------------------------------------
472              zdistante = zdistnew
473
474              ! c. compute distance between the 4 adjacent points and endingPoint
475              !------------------------------------------------------------------
476              ! BE CARREFUL! When the domain crosses the date line (ll_date_domain):
477              ! When we will compute distances between W/E/S/N points and endingPoint,
478              ! we have to check if theses 4 lines crosses the date line
479              ! (test: ABS(coordTemp%lon - coordLast%lon).GT. 180)
480              ! If it's true,we have to add 360° to coordLast%long to compute the
481              ! good distance through the date line and not through the complementary
482              ! in the mesh.
483       
484              ! c.1 compute distWest: distance between west point and endingPoint
485              !----------------------
486              zdistWest2 = 99999999.9 ;  zdistWest3 = 99999999.9 
487              IF( sec%listPoint(jseg)%I .NE. 1 )THEN
488                 !When we are on the west side of the mesh (i=1),we can't go to the west.
489                 coordTemp = pointToCoordF(WestPoint) 
490                 ll_test = .FALSE.
491                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test = .TRUE.
492                 zdistWest2  = distance2(pointToCoordF(WestPoint) ,coordLast ,ll_test)
493              ENDIF
494
495              ! c.2 compute distEst: distance between east point and endingPoint
496              !---------------------
497              zdistEst2 = 99999999.9 ;  zdistEst3 = 99999999.9
498              IF( sec%listPoint(jseg)%I .EQ. jpi )THEN
499                 !We test if the current point is on the east side of the mesh
500                 ! The method is done such as we go toward east to link
501                 ! sec%listPoint(1) to endingPoint.
502                 ! So, if the section crosses the overlapping band (ll_overlap_sec=T),
503                 ! we won't have to stop if the current point is on the EAST side of the mesh:
504                 ! we have to follow the construction of the section on the
505                 ! WEST side of the mesh
506                 IF( ll_overlap_sec  )THEN
507                    !section crosses the overlapping band
508                    !So EstPoint is on the west side of the mesh
509                    EstPoint = POINT_SECTION(3,sec%listPoint(jseg)%J)
510                    zdistEst2= distance2(pointToCoordF(EstPoint)  ,coordLast ,.FALSE.)
511                 ENDIF
512              ELSE
513                 coordTemp = pointToCoordF(EstPoint) 
514                 ll_test = .FALSE.
515                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
516                 zdistEst2 = distance2(pointToCoordF(EstPoint)  ,coordLast ,ll_test )
517              ENDIF
518
519              ! c.3 compute distSouth: distance between south point and endingPoint
520              !-----------------------
521              zdistSouth2 = 99999999.9 ; zdistSouth3 = 99999999.9
522              IF( sec%listPoint(jseg)%J .NE. 1 )THEN
523                 !When we are on the south side of the mesh (j=1),we can't go to the south.
524                 coordTemp=pointToCoordF(SouthPoint)
525                 ll_test = .FALSE.
526                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
527                    zdistSouth2 = distance2(pointToCoordF(SouthPoint),coordlast ,ll_test )
528              ENDIF
529
530              ! c.4 compute distNorth: distance between north and endingPoint
531              !-----------------------
532              zdistNorth2 = 99999999.9 ; zdistNorth3 = 99999999.9 
533              IF( sec%listPoint(jseg)%J .NE. jpj )THEN
534                 !When we are on the north side of the mesh (j=jpj),we can't go to the south.
535                 coordTemp=pointToCoordF(NorthPoint)
536                 ll_test = .FALSE.
537                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
538                 zdistNorth2 = distance2(pointToCoordF(NorthPoint),coordlast ,ll_test )
539              ENDIF
540
541              ! d. select the points which are closer to end-point than current point
542              !----------------------------------------------------------------------
543              zdistref=distance2(pointToCoordF(sec%listPoint(jseg)),coordlast ,ll_test )
544              WRITE(cltmp,'( A56,f10.3 )' )'distance between actual point and last point: zdistref = ',zdistref
545              CALL write_debug(jsec,cltmp)
546              lest   = .FALSE. ; IF( zdistEst2   .LE. zdistref ) lest  = .TRUE.
547              lwest  = .FALSE. ; IF( zdistwest2  .LE. zdistref ) lwest = .TRUE.
548              lnorth = .FALSE. ; IF( zdistnorth2 .LE. zdistref ) lnorth= .TRUE.
549              lsouth = .FALSE. ; IF( zdistsouth2 .LE. zdistref ) lsouth= .TRUE.
550
551              !debug
552              IF( .NOT. lest   )CALL write_debug(jsec,'Est   point eliminated')
553              IF( .NOT. lwest  )CALL write_debug(jsec,'West  point eliminated')
554              IF( .NOT. lnorth )CALL write_debug(jsec,'North point eliminated')
555              IF( .NOT. lsouth )CALL write_debug(jsec,'South point eliminated')
556
557              l_oldmethod = .FALSE.
558
559              IF( ( COUNT((/lest/))+COUNT((/lwest/))+COUNT((/lnorth/))+COUNT((/lsouth/)) ) .NE. 0 )THEN
560
561                 ! e.1 If at least one point is selected, select the point
562                 !     which is the closest to original section among selected points
563                 !-------------------------------------------------------------------
564
565                 zdistWest3  = 9999999.9
566                 IF( lwest )zdistWest3  = &
567                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(WestPoint) ,pointToCoordF(endingPoint) ,lkdebug )
568                 zdistEst3   = 9999999.9
569                 IF( lest )zdistEst3   =  &
570                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(EstPoint)  ,pointToCoordF(endingPoint) ,lkdebug )
571                 zdistSouth3 = 9999999.9
572                 IF( lsouth )zdistSouth3 = &
573                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(SouthPoint),pointToCoordF(endingPoint) ,lkdebug )
574                 zdistNorth3 = 9999999.9
575                 IF( lnorth )zdistNorth3 = &
576                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(NorthPoint),pointToCoordF(endingPoint) ,lkdebug )
577
578                 zdistEst3   = zdistEst3   + (1-COUNT((/lest/))  )*9999999.9
579                 zdistWest3  = zdistWest3  + (1-COUNT((/lwest/)) )*9999999.9
580                 zdistNorth3 = zdistNorth3 + (1-COUNT((/lnorth/)))*9999999.9
581                 zdistSouth3 = zdistSouth3 + (1-COUNT((/lsouth/)))*9999999.9
582
583                 zdistnew = MIN(zdistEst3,zdistWest3,zdistnorth3,zdistSouth3)
584
585              ELSE 
586
587                 ! e.2 If no point is selected, select the point which is the closest to end-point
588                 !--------------------------------------------------------------------------------
589                 l_oldmethod = .TRUE.
590
591                 !debug
592                 WRITE(cltmp,'(A30,i3.3)' )'SEARCH NEW POINT WITH OLD METHOD: ',jsec
593                 CALL write_debug(jsec,cltmp)
594                   
595                 ! on passe à l'ancienne methode: le point parmies les 4 pts (NSWE)  qui se rapproche
596                 ! le + du dernier pt
597                 !-----------------------------
598                 !be carreful! we can't go backward.
599
600                 zdistNorth = zdistNorth2    ; zdistSouth = zdistSouth2
601                 zdistEst   = zdistEst2      ; zdistWest  = zdistWest2 
602
603                 IF(     prevPoint%I .EQ. NorthPoint%I .AND. prevPoint%J .EQ. NorthPoint%J) THEN
604                     zdistnew = MIN(zdistEst,zdistWest,zdistSouth)
605                 ELSE IF(prevPoint%I .EQ. SouthPoint%I .AND. prevPoint%J .EQ. SouthPoint%J) THEN
606                    zdistnew = MIN(zdistEst,zdistWest,zdistNorth)
607                 ELSE IF(prevPoint%I .EQ. WestPoint%I  .AND. prevPoint%J .EQ. WestPoint%J ) THEN
608                    zdistnew = MIN(zdistEst,zdistNorth,zdistSouth)
609                 ELSE IF(prevPoint%I .EQ. EstPoint%I   .AND. prevPoint%J .EQ. EstPoint%J  ) THEN
610                    zdistnew = MIN(zdistWest,zdistNorth,zdistSouth)
611                 ELSE
612                    zdistnew = MIN(zdistEst,zdistWest,zdistNorth,zdistSouth)
613                 ENDIF             
614
615              ENDIF
616
617              !debug
618              WRITE(cltmp,'(A11, f8.3)')'zdistref = ',zdistref
619              CALL write_debug(jsec,cltmp)
620              WRITE(cltmp, 103         )'distance2 :',zdistEst2,zdistWest2,zdistNorth2,zdistSouth2 
621              CALL write_debug(jsec,cltmp)
622              WRITE(cltmp, 103         )'distance3 :',zdistEst3,zdistWest3,zdistNorth3,zdistSouth3
623              CALL write_debug(jsec,cltmp)
624              WRITE(cltmp,'(A11, f8.3)')"zdistnew = ",zdistnew
625              CALL write_debug(jsec,cltmp)
626
627103 FORMAT (A12,"E",f12.3," W",f12.3," N",f12.3," S",f12.3)
628
629              !f. save next point and direction of velocity.
630              !---------------------------------------------
631              !nextPoint will be the one which is the closest to endingPoint.
632              !sec%direction will be direction between current point and nextPoint:
633              !It will be used to compute velocity through the segment [CurrentPoint,nextPoint}:
634              !                 
635              !A:Current Point    NorthPoint(I,J+1)   Nextpoint=NorthPoint(I,J+1) => sec%direction=3
636              !                   |                   Nextpoint=SouthPoint(I,J-1) => sec%direction=2   
637              !                   |                   Nextpoint=WestPoint(I-1,J)  => sec%direction=0   
638              !                   |==>V(I,J+1)        Nextpoint=EastPoint(I+1,J)  => sec%direction=1
639              !            U(I,J) |       U(I+1,J)
640              !            ^      |       ^
641              !   West_____|______A_______|_____EstPoint
642              !   Point           |(I,J)        (I+1,J) 
643              !   (I-1,J)         |               
644              !                   |==>V(I,J)
645              !                   |
646              !              SoutPoint(I,J-1)
647              IF( l_oldmethod )THEN
648                 IF( zdistnew == zdistWest )      THEN
649                      sec%direction(jseg)=0 ; nextPoint = WestPoint
650                 ELSE IF( zdistnew == zdistEst )  THEN
651                      sec%direction(jseg)=1 ; nextPoint = EstPoint
652                 ELSE IF( zdistnew == zdistSouth )THEN
653                      sec%direction(jseg)=2 ; nextPoint = SouthPoint
654                 ELSE IF( zdistnew == zdistNorth )THEN
655                      sec%direction(jseg)=3 ; nextPoint= NorthPoint
656                 ENDIF
657              ELSE
658                 IF( zdistnew == zdistWest3 )      THEN
659                      sec%direction(jseg)=0 ; nextPoint = WestPoint
660                 ELSE IF( zdistnew == zdistEst3 )  THEN
661                      sec%direction(jseg)=1 ; nextPoint = EstPoint
662                 ELSE IF( zdistnew == zdistSouth3 )THEN
663                      sec%direction(jseg)=2 ; nextPoint = SouthPoint
664                 ELSE IF( zdistnew == zdistNorth3 )THEN
665                      sec%direction(jseg)=3 ; nextPoint= NorthPoint
666                 ENDIF
667              ENDIF
668
669              WRITE(cltmp,'(A11, 2(i4.4,1X) )')'nextPoint = ', nextPoint
670              CALL write_debug(jsec,cltmp)
671       
672              !f. Save nextPoint and go to nextPoint
673              !-------------------------------------
674              prevPoint = sec%listPoint(jseg)
675              jseg = jseg+1                   !increment of number of segments that form the section
676              sec%listPoint(jseg) = nextPoint !Save next point
677   
678           ENDIF ! southP/northP/WestP/EstP == endingpoint ?
679
680        ENDDO                  !End of loop on jseg
681        sec%nb_point = jseg    !Save the number of segments that form the section
682
683
684     ELSE ! nb_inmesh == 0
685        DO jseg=1,nb_point_max 
686           sec%listPoint(:)=POINT_SECTION(0,0)
687        ENDDO
688        sec%direction(:)=0.
689        sec%nb_point = 0
690     ENDIF
691
692     !debug     
693     IF( sec%nb_point .ne. 0 )THEN
694        CALL write_debug(jsec,"-------------------------------------")
695        CALL write_debug(jsec,"list of points in the grid : ")
696        DO jseg=1,sec%nb_point 
697           ji=sec%listPoint(jseg)%I ; jj=sec%listPoint(jseg)%J
698           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)
699           CALL write_debug(jsec,cltmp)
700        ENDDO
701   
702        !test if we are one end-point
703        IF( sec%listPoint(sec%nb_point)%I .NE. endingPoint%J .AND. sec%listPoint(sec%nb_point)%J .NE. endingPoint%J )THEN   
704           PRINT*,TRIM(sec%name)," NOT ARRIVED TO endingPoint FOR jsec =  ",jsec
705        ENDIF
706     ENDIF
707
708     !now compute new slopeSection with ij-coordinates of first and last point
709     IF( sec%nb_point .ne. 0 )THEN
710        IF (  sec%listPoint(sec%nb_point)%I .NE.  sec%listPoint(1)%I ) THEN
711           sec%slopeSection = ( sec%listPoint(sec%nb_point)%J - sec%listPoint(1)%J ) /  &
712                              ( sec%listPoint(sec%nb_point)%I - sec%listPoint(1)%I )     
713        ELSE
714           sec%slopeSection = 10000._wp
715        ENDIF
716     ENDIF
717 
718  END SUBROUTINE compsec
719
720END MODULE compute_sections 
Note: See TracBrowser for help on using the repository browser.