Ignore:
Timestamp:
2011-09-26T10:54:50+02:00 (9 years ago)
Author:
cbricaud
Message:

cleanning, minor modifications

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90

    r2849 r2858  
    5858        zdistante  , zdistante2 , zdistnew  , zdistnew2  ,  &!         " 
    5959        zdeltai    , zdeltaj                                !         " 
    60 !     REAL(wp),DIMENSION(:,:),ALLOCATABLE,SAVE :: zmask 
    6160     LOGICAL :: &  
    6261        ll_overlap_sec_left = .FALSE. , ll_overlap_sec_right = .FALSE. ,&! temporary logical 
     
    7574        SouthPoint , NorthPoint, EstPoint , WestPoint  !        " 
    7675     !!--------------------------------------------------------------------- 
    77      IF ( jsec==1 )THEN 
     76     IF( jsec==1 )THEN 
    7877          PRINT*,'                ' 
    7978          PRINT*,'COMPUTE SECTIONS' 
     
    104103     SouthPoint  = POINT_SECTION( -1, -1 ) ; NorthPoint  = POINT_SECTION( -1, -1 ) 
    105104     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 ) 
     105 
     106     IF( jsec == 1 )THEN 
     107        !Found the taller cell of the mesh in ocean (used in compsec) 
     108        zdistmesh=0. 
     109        DO jj=2,jpj 
     110           DO ji=2,jpi 
     111              IF( zdistmesh .LT. & 
    117112                 (e1t(ji,jj)/(cos(gphit(ji,jj))*110000))  ) THEN 
    118113                 zdistmesh = e1t(ji,jj)/(cos(gphit(ji,jj))*110000) 
    119                ENDIF 
    120             ENDDO 
    121          ENDDO 
     114              ENDIF 
     115           ENDDO 
     116        ENDDO 
    122117     
    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  
    138118     ENDIF 
    139      !PRINT*,"min max zmask ",MINVAL(zmask),MAXVAL(zmask) 
    140119 
    141120     !debug 
     
    148127 
    149128     !loop on the mesh  
    150      DO jj=2,nlcj 
    151         DO ji=2,nlci 
     129     DO jj=2,jpj 
     130        DO ji=2,jpi 
    152131           !----------------------------------------------------------- 
    153132           !For each cell of the mesh, we find the intersection between  
     
    176155           zdeltaj= gphif(ji,jj-1) - gphif(ji-1,jj-1) !zdeltaj=[AD] 
    177156            
    178            IF (ABS(zdeltai) .LE. 2.*zdistmesh .OR.  &  
    179                ABS(zdeltaj) .LE. 2.*zdistmesh) THEN            
     157           IF( ABS(zdeltai) .LE. 2.*zdistmesh .OR.  &  
     158               ABS(zdeltaj) .LE. 2.*zdistmesh )THEN            
    180159   
    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)) &  
     160              !intersection section/[AB]=? 
     161              coordTemp     = intersec(sec,coord_a,coord_b)     !compute intersection  
     162              coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj)) 
     163              coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj)) 
     164  
     165              IF( coordTemp%lon .NE. -9999 )THEN 
     166                 IF( nb_inmesh+1 .GT. nb_point_max )THEN 
     167                    PRINT*,"WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh 
     168                 ELSE     
     169                    nb_inmesh=nb_inmesh+1               
     170                    coordSec(nb_inmesh) = coordTemp    !store the intersection's coordinates 
     171 
     172                    !We need to know if the section crosses the overlapping band. 
     173 
     174                    !Fist we look if there is an intersection mesh/section  
     175                    !just on the left of the overlapping band:  
     176                    IF( coordTemp%lon .GT. glamf(1,1)-5  .AND. coordTemp%lon .LT. glamf(1,1) ) &  
     177                       & ll_overlap_sec_left  = .TRUE. 
     178                    !And we look if there is an intersection mesh/section  
     179                    !just on the right of the overlapping band: 
     180                    IF( coordTemp%lon .GT. glamf(jpi,1) .AND. coordTemp%lon .LT. glamf(1,1)+5 ) & 
     181                       & ll_overlap_sec_right = .TRUE. 
     182                 ENDIF 
     183              ENDIF 
     184               
     185              !intersection section/[AD]=? 
     186              coordTemp=intersec(sec,coord_a,coord_d)       !compute intersection 
     187              coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj)) 
     188              coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj)) 
     189 
     190              IF( coordTemp%lon .NE. -9999 )THEN 
     191                 IF( nb_inmesh+1 .GT. nb_point_max )THEN 
     192                    PRINT*, "WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh 
     193                 ELSE   
     194                    nb_inmesh=nb_inmesh+1 
     195                    coordSec(nb_inmesh)=coordTemp 
     196                    
     197                    !We need to know if the section crosses the overlapping band: 
     198                    !same test as above  
     199                    IF( coordTemp%lon .GE. glamf(1,1)-3  .AND. coordTemp%lon .LE. glamf(1,1) ) & 
    198200                         & 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) &  
     201                    IF( coordTemp%lon .GE. glamf(jpi,1) .AND. coordTemp%lon .LE. glamf(jpi,1)+3) &  
    223202                         & ll_overlap_sec_right = .TRUE.  
    224                    ENDIF 
    225                ENDIF 
    226   
    227           ENDIF 
     203                 ENDIF 
     204              ENDIF 
     205  
     206           ENDIF 
    228207          
    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 
     208           !We need to know if the domain crosses the date line: 
     209           !Fist, we search a mesh point that is just one the left of date line: 
     210           IF( glamf(ji-1,jj-1) .GT.  175 .AND. glamf(ji-1,jj-1) .LT.  180 ) & 
     211              & ll_date_domain_left = .TRUE. 
     212           !And we search a mesh point that is just one the right of date line: 
     213           IF( glamf(ji-1,jj-1) .GT. -180 .AND. glamf(ji-1,jj-1) .LT. -175 ) &   
     214              & ll_date_domain_right = .TRUE. 
     215  
     216        ENDDO 
     217     ENDDO !End of the loop on the mesh 
     218 
    240219  
    241220     !Crossing section/overlapping band (we need to know it for later): 
     
    266245 
    267246     IF( sec%ll_date_line .AND. ll_date_domain )THEN 
     247 
    268248        !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 
     249        DO jseg=1,nb_inmesh 
     250           IF( coordSec(jseg)%lon .LT. 0 ) coordSec(jseg)%lon=coordSec(jseg)%lon+360. 
     251        ENDDO 
     252        IF( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon ) THEN 
     253           CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh) 
     254        ELSE 
     255           CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh) 
     256        ENDIF 
     257        DO jseg=1,nb_inmesh 
     258           IF( coordSec(jseg)%lon .GT. 180 ) coordSec(jseg)%lon=coordSec(jseg)%lon-360. 
     259        ENDDO 
     260 
    280261     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 
     262 
     263        IF( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon )THEN 
     264           CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh) 
     265        ELSE 
     266           CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh) 
     267        ENDIF 
     268 
    286269     ENDIF 
    287270 
     
    300283 
    301284     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,nle 
    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 
     285        coordFirst       = coordSec(1) 
     286        coordLast        = coordSec(nb_inmesh)  
     287        sec%nb_point     = nb_inmesh 
     288        sec%listPoint(1) = POINT_SECTION(-1,-1) 
     289        zdistante        = 1000. 
     290        zdistante2       = 1000. 
     291 
     292        !First, we find the point of the mesh that is the closest  
     293        !to the first intersection section/mesh (=coordFirst=coordSec(1)): 
     294        !this point will be called sec%listPoint(1). 
     295        !Then, we find the point of the mesh that is the closest 
     296        !to the last intersection section/mesh (coordLast=coordSec(nb_inmesh)) 
     297        !this point will be called endingPoint. 
     298 
     299        DO jj=1,jpj 
     300           DO ji=1,jp 
     301              coord_t=COORD_SECTION(glamf(ji,jj),gphif(ji,jj)) 
     302              zdistFirst = distance(coord_t,coordFirst) 
     303              zdistLast = distance(coord_t,coordLast) 
     304              IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN 
     305                 sec%listPoint(1) = POINT_SECTION(ji,jj) 
     306                 zdistante=zdistFirst 
     307              ENDIF 
     308              IF( zdistLast .LT. zdistmesh .AND. zdistLast .LT. zdistante2 )THEN 
     309                 endingPoint = POINT_SECTION(ji,jj) 
     310                 zdistante2=zdistLast 
     311              ENDIF 
     312           ENDDO 
     313        ENDDO 
     314 
     315        IF( sec%listPoint(1)%I == endingPoint%I .AND. sec%listPoint(1)%J == endingPoint%J )THEN 
     316           sec%listPoint(1) = POINT_SECTION(-1,-1) 
     317           endingPoint      = POINT_SECTION(-1,-1) 
     318           coordFirst       = coordSec(1) 
     319           coordLast        = coordSec(2) 
     320           sec%nb_point     = 0 
     321        ENDIF 
    339322 
    340323     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 
     324        !If there is no intersection section/mesh 
     325        sec%listPoint(1) = POINT_SECTION(-1,-1) 
     326        endingPoint      = POINT_SECTION(-1,-1) 
     327        coordFirst       = coordSec(1) 
     328        coordLast        = coordSec(2) 
     329        sec%nb_point     = 0 
    347330     ENDIF  
    348331 
     
    380363     IF( nb_inmesh .NE. 0 )THEN 
    381364 
    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  )          
     365        !The serie of mesh's points that form the section will 'link'  
     366        !sec%listPoint(1) to endingPoint: it will be stored in  
     367        !sec%listPoint(jseg) 
     368        ! 
     369        !We take place on the fist point (sec%listPoint(1))  
     370        ! a.  We find the 4 adjacent points (North, South, East, West) 
     371        ! b.  Compute distance between current point and endingPoint 
     372        ! c.  Compute distance between the 4 adjacent points and endingPoint 
     373        ! d.  Select the points which are closer to end-point than current point 
     374        ! e.1 If at least one point is selected, select the point which is closest to original section among selected points  
     375        ! e.2 If no point is selected, select the point which is the closest to end-point  
     376        ! f. save next point and direction of velocity. 
     377        ! g. Save nextPoint and go to nextPoint 
     378        ! 
     379        !We get out of this loop if: 
     380        !    - we are on endingPoint 
     381        !    - the number of points (jseg) that link sec%listPoint(1) to endingPoint is  
     382        !      twice greater than number of section/mesh intersection (nb_inmesh): 
     383        !      it could be possible if thr algorithm can't link endingPoint (bug). 
     384 
     385        !initialize distnew value (with distance between section's extremities)  
     386        zdistnew  = distance(coordFirst,coordLast,sec%ll_date_line)  
     387        prevPoint  = POINT_SECTION(0,0) 
     388        jseg       = 1 
     389 
     390        DO WHILE ( (  sec%listPoint(jseg)%I .NE.  endingPoint%I    & 
     391                 .OR. sec%listPoint(jseg)%J .NE. endingPoint%J   ) & 
     392                 .AND. jseg .LT. 500 .AND. sec%listPoint(jseg)%I .GT. 0  )          
    410393    
    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) 
     394           ! a. find the 4 adjacent points (North, South, East, West) 
     395           !--------------------------------------------------------- 
     396           SouthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J-1) 
     397           NorthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     398           WestPoint  = POINT_SECTION(sec%listPoint(jseg)%I-1,sec%listPoint(jseg)%J) 
     399           EstPoint   = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     400 
     401           !debug 
     402           CALL write_debug(jsec,"---------------") 
     403           WRITE(cltmp,100)'Current points: ',sec%listPoint(jseg), &  
     404                         glamf(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J), & 
     405                         gphif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) 
     406           CALL write_debug(jsec,cltmp) 
     407           CALL write_debug(jsec,"E/W/N/S points") 
     408           WRITE(cltmp,101)glamf(  EstPoint%I,EstPoint%J)  ,gphif(  EstPoint%I,  EstPoint%J), & 
     409                           glamf( WestPoint%I,WestPoint%J) ,gphif( WestPoint%I, WestPoint%J), & 
     410                           glamf(NorthPoint%I,NorthPoint%J),gphif(NorthPoint%I,NorthPoint%J) ,& 
     411                           glamf(SouthPoint%I,SouthPoint%J),gphif(SouthPoint%I,SouthPoint%J) 
     412           CALL write_debug(jsec,cltmp) 
     413           WRITE(cltmp,102)EstPoint,WestPoint,NorthPoint,SouthPoint 
     414           CALL write_debug(jsec,cltmp) 
    432415 
    433416100 FORMAT ( A15,2(i4.4," "),2(f7.3," ") ) 
    434417101 FORMAT ( "E ",2(f7.3," "),"W ",2(f7.3," "),"N ",2(f7.3," "),"S ",2(f7.3," ")) 
    435 102 FORMAT ( "E ",i4.4,' 'i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 ) 
     418102 FORMAT ( "E ",i4.4,' ',i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 ) 
    436419 
    437420                
    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. 
     421           !we  are on end-point 
     422           !-------------------- 
     423           IF(      SouthPoint%I==endingPoint%I .AND. SouthPoint%J==endingPoint%J )THEN  
     424               jseg = jseg+1 ; sec%listPoint(jseg) = SouthPoint 
     425           ELSE IF( NorthPoint%I==endingPoint%I .AND. NorthPoint%J==endingPoint%J )THEN 
     426               jseg = jseg+1 ; sec%listPoint(jseg) = NorthPoint 
     427           ELSE IF(  WestPoint%I==endingPoint%I .AND.  WestPoint%J==endingPoint%J )THEN 
     428               jseg = jseg+1 ; sec%listPoint(jseg) = WestPoint 
     429           ELSE IF(   EstPoint%I==endingPoint%I .AND.   EstPoint%J==endingPoint%J )THEN 
     430               jseg = jseg+1 ; sec%listPoint(jseg) = EstPoint 
     431 
     432           ELSE 
     433           !we  are NOT on end-point 
     434           !------------------------ 
     435 
     436              ! b. distance between current point and endingPoint 
     437              !-------------------------------------------------- 
     438              zdistante = zdistnew 
     439 
     440              ! c. compute distance between the 4 adjacent points and endingPoint 
     441              !------------------------------------------------------------------ 
     442              ! BE CARREFUL! When the domain crosses the date line (ll_date_domain): 
     443              ! When we will compute distances between W/E/S/N points and endingPoint, 
     444              ! we have to check if theses 4 lines crosses the date line 
     445              ! (test: ABS(coordTemp%lon - coordLast%lon).GT. 180) 
     446              ! If it's true,we have to add 360° to coordLast%long to compute the  
     447              ! good distance through the date line and not through the complementary 
     448              ! in the mesh. 
    466449        
    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) 
     450              ! c.1 compute distWest: distance between west point and endingPoint 
     451              !---------------------- 
     452              zdistWest2 = 99999999.9 ;  zdistWest3 = 99999999.9  
     453              IF( sec%listPoint(jseg)%I .NE. 1 )THEN 
     454                 !When we are on the west side of the mesh (i=1),we can't go to the west. 
     455                 coordTemp = pointToCoordF(WestPoint)  
     456                 ll_test = .FALSE. 
     457                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test = .TRUE. 
     458                 zdistWest2  = distance2(pointToCoordF(WestPoint) ,coordLast ,ll_test) 
     459              ENDIF 
     460 
     461              ! c.2 compute distEst: distance between east point and endingPoint 
     462              !--------------------- 
     463              zdistEst2 = 99999999.9 ;  zdistEst3 = 99999999.9 
     464              IF( sec%listPoint(jseg)%I .EQ. jpi )THEN 
     465                 !We test if the current point is on the east side of the mesh 
     466                 ! The method is done such as we go toward east to link 
     467                 ! sec%listPoint(1) to endingPoint. 
     468                 ! So, if the section crosses the overlapping band (ll_overlap_sec=T), 
     469                 ! we won't have to stop if the current point is on the EAST side of the mesh: 
     470                 ! we have to follow the construction of the section on the 
     471                 ! WEST side of the mesh 
     472                 IF( ll_overlap_sec  )THEN 
     473                    !section crosses the overlapping band  
     474                    !So EstPoint is on the west side of the mesh 
     475                    EstPoint = POINT_SECTION(3,sec%listPoint(jseg)%J) 
     476                    zdistEst2= distance2(pointToCoordF(EstPoint)  ,coordLast ,.FALSE.) 
     477                 ENDIF  
     478              ELSE 
     479                 coordTemp = pointToCoordF(EstPoint)  
     480                 ll_test = .FALSE. 
     481                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE. 
     482                 zdistEst2 = distance2(pointToCoordF(EstPoint)  ,coordLast ,ll_test ) 
     483              ENDIF 
     484 
     485              ! c.3 compute distSouth: distance between south point and endingPoint 
     486              !----------------------- 
     487              zdistSouth2 = 99999999.9 ; zdistSouth3 = 99999999.9 
     488              IF( sec%listPoint(jseg)%J .NE. 1 )THEN 
     489                 !When we are on the south side of the mesh (j=1),we can't go to the south.  
     490                 coordTemp=pointToCoordF(SouthPoint) 
     491                 ll_test = .FALSE. 
     492                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE. 
     493                    zdistSouth2 = distance2(pointToCoordF(SouthPoint),coordlast ,ll_test ) 
     494              ENDIF 
     495 
     496              ! c.4 compute distNorth: distance between north and endingPoint 
     497              !----------------------- 
     498              zdistNorth2 = 99999999.9 ; zdistNorth3 = 99999999.9  
     499              IF( sec%listPoint(jseg)%J .NE. jpj )THEN 
     500                 !When we are on the north side of the mesh (j=jpj),we can't go to the south.  
     501                 coordTemp=pointToCoordF(NorthPoint) 
     502                 ll_test = .FALSE. 
     503                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE. 
     504                 zdistNorth2 = distance2(pointToCoordF(NorthPoint),coordlast ,ll_test ) 
     505              ENDIF 
     506 
     507              ! d. select the points which are closer to end-point than current point 
     508              !---------------------------------------------------------------------- 
     509              zdistref=distance2(pointToCoordF(sec%listPoint(jseg)),coordlast ,ll_test ) 
     510              WRITE(cltmp,'( A56,f10.3 )' )'distance between actual point and last point: zdistref = ',zdistref 
     511              CALL write_debug(jsec,cltmp) 
     512              lest   = .FALSE. ; IF( zdistEst2   .LE. zdistref ) lest  = .TRUE. 
     513              lwest  = .FALSE. ; IF( zdistwest2  .LE. zdistref ) lwest = .TRUE. 
     514              lnorth = .FALSE. ; IF( zdistnorth2 .LE. zdistref ) lnorth= .TRUE. 
     515              lsouth = .FALSE. ; IF( zdistsouth2 .LE. zdistref ) lsouth= .TRUE. 
     516 
     517              !debug 
     518              IF( .NOT. lest   )CALL write_debug(jsec,'Est   point eliminated') 
     519              IF( .NOT. lwest  )CALL write_debug(jsec,'West  point eliminated') 
     520              IF( .NOT. lnorth )CALL write_debug(jsec,'North point eliminated') 
     521              IF( .NOT. lsouth )CALL write_debug(jsec,'South point eliminated') 
     522 
     523              l_oldmethod = .FALSE. 
     524 
     525              IF( ( COUNT((/lest/))+COUNT((/lwest/))+COUNT((/lnorth/))+COUNT((/lsouth/)) ) .NE. 0 )THEN 
     526 
     527                 ! e.1 If at least one point is selected, select the point  
     528                 !     which is the closest to original section among selected points  
     529                 !------------------------------------------------------------------- 
     530 
     531                 zdistWest3  = 9999999.9 
     532                 IF( lwest )zdistWest3  = & 
     533                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(WestPoint) ,pointToCoordF(endingPoint) ,lkdebug ) 
     534                 zdistEst3   = 9999999.9 
     535                 IF( lest )zdistEst3   =  & 
     536                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(EstPoint)  ,pointToCoordF(endingPoint) ,lkdebug ) 
     537                 zdistSouth3 = 9999999.9 
     538                 IF( lsouth )zdistSouth3 = & 
     539                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(SouthPoint),pointToCoordF(endingPoint) ,lkdebug ) 
     540                 zdistNorth3 = 9999999.9 
     541                 IF( lnorth )zdistNorth3 = & 
     542                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(NorthPoint),pointToCoordF(endingPoint) ,lkdebug ) 
     543 
     544                 zdistEst3   = zdistEst3   + (1-COUNT((/lest/))  )*9999999.9 
     545                 zdistWest3  = zdistWest3  + (1-COUNT((/lwest/)) )*9999999.9 
     546                 zdistNorth3 = zdistNorth3 + (1-COUNT((/lnorth/)))*9999999.9 
     547                 zdistSouth3 = zdistSouth3 + (1-COUNT((/lsouth/)))*9999999.9 
     548 
     549                 zdistnew = MIN(zdistEst3,zdistWest3,zdistnorth3,zdistSouth3) 
     550 
     551              ELSE  
     552 
     553                 ! e.2 If no point is selected, select the point which is the closest to end-point  
     554                 !-------------------------------------------------------------------------------- 
     555                 l_oldmethod = .TRUE. 
     556 
     557                 !debug 
     558                 WRITE(cltmp,'(A30,i3.3)' )'SEARCH NEW POINT WITH OLD METHOD: ',jsec 
     559                 CALL write_debug(jsec,cltmp) 
    577560                    
    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) 
     561                 ! on passe à l'ancienne methode: le point parmies les 4 pts (NSWE)  qui se rapproche 
     562                 ! le + du dernier pt 
     563                 !-----------------------------  
     564                 !be carreful! we can't go backward. 
     565 
     566                 zdistNorth = zdistNorth2    ; zdistSouth = zdistSouth2 
     567                 zdistEst   = zdistEst2      ; zdistWest  = zdistWest2  
     568 
     569                 IF(     prevPoint%I .EQ. NorthPoint%I .AND. prevPoint%J .EQ. NorthPoint%J) THEN 
     570                     zdistnew = MIN(zdistEst,zdistWest,zdistSouth) 
     571                 ELSE IF(prevPoint%I .EQ. SouthPoint%I .AND. prevPoint%J .EQ. SouthPoint%J) THEN 
     572                    zdistnew = MIN(zdistEst,zdistWest,zdistNorth) 
     573                 ELSE IF(prevPoint%I .EQ. WestPoint%I  .AND. prevPoint%J .EQ. WestPoint%J ) THEN 
     574                    zdistnew = MIN(zdistEst,zdistNorth,zdistSouth) 
     575                 ELSE IF(prevPoint%I .EQ. EstPoint%I   .AND. prevPoint%J .EQ. EstPoint%J  ) THEN 
     576                    zdistnew = MIN(zdistWest,zdistNorth,zdistSouth) 
     577                 ELSE  
     578                    zdistnew = MIN(zdistEst,zdistWest,zdistNorth,zdistSouth) 
     579                 ENDIF               
     580 
     581              ENDIF 
     582 
     583              !debug 
     584              WRITE(cltmp,'(A11, f8.3)')'zdistref = ',zdistref 
     585              CALL write_debug(jsec,cltmp) 
     586              WRITE(cltmp, 103         )'distance2 :',zdistEst2,zdistWest2,zdistNorth2,zdistSouth2  
     587              CALL write_debug(jsec,cltmp) 
     588              WRITE(cltmp, 103         )'distance3 :',zdistEst3,zdistWest3,zdistNorth3,zdistSouth3 
     589              CALL write_debug(jsec,cltmp) 
     590              WRITE(cltmp,'(A11, f8.3)')"zdistnew = ",zdistnew 
     591              CALL write_debug(jsec,cltmp) 
    609592 
    610593103 FORMAT (A12,"E",f12.3," W",f12.3," N",f12.3," S",f12.3) 
    611594 
    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) 
     595              !f. save next point and direction of velocity. 
     596              !---------------------------------------------  
     597              !nextPoint will be the one which is the closest to endingPoint. 
     598              !sec%direction will be direction between current point and nextPoint: 
     599              !It will be used to compute velocity through the segment [CurrentPoint,nextPoint}: 
     600              !                  
     601              !A:Current Point    NorthPoint(I,J+1)   Nextpoint=NorthPoint(I,J+1) => sec%direction=3 
     602              !                   |                   Nextpoint=SouthPoint(I,J-1) => sec%direction=2    
     603              !                   |                   Nextpoint=WestPoint(I-1,J)  => sec%direction=0    
     604              !                   |==>V(I,J+1)        Nextpoint=EastPoint(I+1,J)  => sec%direction=1 
     605              !            U(I,J) |       U(I+1,J) 
     606              !            ^      |       ^ 
     607              !   West_____|______A_______|_____EstPoint  
     608              !   Point           |(I,J)        (I+1,J)   
     609              !   (I-1,J)         |                
     610              !                   |==>V(I,J) 
     611              !                   | 
     612              !              SoutPoint(I,J-1) 
     613              IF( l_oldmethod )THEN 
     614                 IF( zdistnew == zdistWest )      THEN 
     615                      sec%direction(jseg)=0 ; nextPoint = WestPoint 
     616                 ELSE IF( zdistnew == zdistEst )  THEN 
     617                      sec%direction(jseg)=1 ; nextPoint = EstPoint 
     618                 ELSE IF( zdistnew == zdistSouth )THEN 
     619                      sec%direction(jseg)=2 ; nextPoint = SouthPoint 
     620                 ELSE IF( zdistnew == zdistNorth )THEN 
     621                      sec%direction(jseg)=3 ; nextPoint= NorthPoint 
     622                 ENDIF 
     623              ELSE  
     624                 IF( zdistnew == zdistWest3 )      THEN 
     625                      sec%direction(jseg)=0 ; nextPoint = WestPoint 
     626                 ELSE IF( zdistnew == zdistEst3 )  THEN 
     627                      sec%direction(jseg)=1 ; nextPoint = EstPoint 
     628                 ELSE IF( zdistnew == zdistSouth3 )THEN 
     629                      sec%direction(jseg)=2 ; nextPoint = SouthPoint 
     630                 ELSE IF( zdistnew == zdistNorth3 )THEN 
     631                      sec%direction(jseg)=3 ; nextPoint= NorthPoint 
     632                 ENDIF  
     633              ENDIF 
     634 
     635              WRITE(cltmp,'(A11, 2(i4.4,1X) )')'nextPoint = ', nextPoint 
     636              CALL write_debug(jsec,cltmp) 
    653637         
    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 
     638              !f. Save nextPoint and go to nextPoint 
     639              !------------------------------------- 
     640              prevPoint = sec%listPoint(jseg) 
     641              jseg = jseg+1                   !increment of number of segments that form the section 
     642              sec%listPoint(jseg) = nextPoint !Save next point 
    659643    
    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 
     644           ENDIF ! southP/northP/WestP/EstP == endingpoint ? 
     645 
     646        ENDDO                  !End of loop on jseg 
     647        sec%nb_point = jseg    !Save the number of segments that form the section 
    664648 
    665649 
    666650     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 
     651        DO jseg=1,nb_point_max  
     652           sec%listPoint(:)=POINT_SECTION(0,0) 
     653        ENDDO 
     654        sec%direction(:)=0. 
     655        sec%nb_point = 0 
    672656     ENDIF 
    673657 
     
    676660     CALL write_debug(jsec,"list of points in the grid : ") 
    677661     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) 
     662        ji=sec%listPoint(jseg)%I ; jj=sec%listPoint(jseg)%J 
     663        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) 
     664        CALL write_debug(jsec,cltmp) 
    681665     ENDDO 
    682666     
Note: See TracChangeset for help on using the changeset viewer.