- Timestamp:
- 2011-09-26T10:54:50+02:00 (13 years ago)
- Location:
- branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90
r2849 r2858 58 58 zdistante , zdistante2 , zdistnew , zdistnew2 , &! " 59 59 zdeltai , zdeltaj ! " 60 ! REAL(wp),DIMENSION(:,:),ALLOCATABLE,SAVE :: zmask61 60 LOGICAL :: & 62 61 ll_overlap_sec_left = .FALSE. , ll_overlap_sec_right = .FALSE. ,&! temporary logical … … 75 74 SouthPoint , NorthPoint, EstPoint , WestPoint ! " 76 75 !!--------------------------------------------------------------------- 77 IF 76 IF( jsec==1 )THEN 78 77 PRINT*,' ' 79 78 PRINT*,'COMPUTE SECTIONS' … … 104 103 SouthPoint = POINT_SECTION( -1, -1 ) ; NorthPoint = POINT_SECTION( -1, -1 ) 105 104 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. & 117 112 (e1t(ji,jj)/(cos(gphit(ji,jj))*110000)) ) THEN 118 113 zdistmesh = e1t(ji,jj)/(cos(gphit(ji,jj))*110000) 119 120 121 114 ENDIF 115 ENDDO 116 ENDDO 122 117 123 !Compute array zmask used later to avoid array overflow124 !ALLOCATE(zmask(jpi,jpj))125 !zmask(:,:)=tmask(:,:,1)126 !DO jj=3,nlcj-2127 ! DO ji=3,nlci-2128 ! 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 ! ENDDO136 !ENDDO137 138 118 ENDIF 139 !PRINT*,"min max zmask ",MINVAL(zmask),MAXVAL(zmask)140 119 141 120 !debug … … 148 127 149 128 !loop on the mesh 150 DO jj=2, nlcj151 DO ji=2, nlci129 DO jj=2,jpj 130 DO ji=2,jpi 152 131 !----------------------------------------------------------- 153 132 !For each cell of the mesh, we find the intersection between … … 176 155 zdeltaj= gphif(ji,jj-1) - gphif(ji-1,jj-1) !zdeltaj=[AD] 177 156 178 IF (ABS(zdeltai) .LE. 2.*zdistmesh .OR. &179 ABS(zdeltaj) .LE. 2.*zdistmesh )THEN157 IF( ABS(zdeltai) .LE. 2.*zdistmesh .OR. & 158 ABS(zdeltaj) .LE. 2.*zdistmesh )THEN 180 159 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) ) & 198 200 & 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) & 223 202 & ll_overlap_sec_right = .TRUE. 224 225 226 227 ENDIF203 ENDIF 204 ENDIF 205 206 ENDIF 228 207 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 mesh238 ENDDO239 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 240 219 241 220 !Crossing section/overlapping band (we need to know it for later): … … 266 245 267 246 IF( sec%ll_date_line .AND. ll_date_domain )THEN 247 268 248 !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 280 261 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 286 269 ENDIF 287 270 … … 300 283 301 284 IF( nb_inmesh .ne. 0 )THEN 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 DO jj=1,nlej317 DO ji=1,nlei318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 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,jpi 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 339 322 340 323 ELSE 341 342 343 344 345 346 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 347 330 ENDIF 348 331 … … 380 363 IF( nb_inmesh .NE. 0 )THEN 381 364 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 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 ) 410 393 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 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) 432 415 433 416 100 FORMAT ( A15,2(i4.4," "),2(f7.3," ") ) 434 417 101 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 )418 102 FORMAT ( "E ",i4.4,' ',i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 ) 436 419 437 420 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 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. 466 449 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 IF( sec%listPoint(jseg)%I .EQ. nlci )THEN482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 IF( sec%listPoint(jseg)%J .NE. nlcj )THEN517 !When we are on the north side of the mesh (j=nlcj),we can't go to the south.518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 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) 577 560 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 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) 609 592 610 593 103 FORMAT (A12,"E",f12.3," W",f12.3," N",f12.3," S",f12.3) 611 594 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) 653 637 654 655 656 657 658 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 659 643 660 661 662 663 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 664 648 665 649 666 650 ELSE ! nb_inmesh == 0 667 668 669 670 671 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 672 656 ENDIF 673 657 … … 676 660 CALL write_debug(jsec,"list of points in the grid : ") 677 661 DO jseg=1,sec%nb_point 678 679 680 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) 681 665 ENDDO 682 666 -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90
r2849 r2858 18 18 INTEGER, PUBLIC, PARAMETER :: numdctout=2 ! Unit for output file 19 19 20 INTEGER, PUBLIC :: jpidta,jpjdta,jpizoom,jpjzoom !dimensions des fichiers lus 21 INTEGER, PUBLIC :: jpi,jpj,jpiglo,jpjglo !dimensions du domaine 22 INTEGER, PUBLIC :: nlei,nlej,nlci,nlcj 20 INTEGER, PUBLIC :: jpi,jpj ! domain dimensions 23 21 INTEGER, PUBLIC :: nb_sec ! Number of section read in input file 24 22 INTEGER, PUBLIC :: nsecdebug = 0 ! Number of the section to debug 25 23 26 INTEGER, PUBLIC, DIMENSION(:) , ALLOCATABLE :: mig27 INTEGER, PUBLIC, DIMENSION(:) , ALLOCATABLE :: mjg28 24 REAL(wp), PUBLIC ,DIMENSION(:,:) , ALLOCATABLE :: glamf,gphif,glamt,gphit,e1t 29 ! REAL(wp), PUBLIC ,DIMENSION(:,:,:), ALLOCATABLE :: tmask30 25 INTEGER, PUBLIC ,DIMENSION(nb_sec_max) :: num_sec_debug 31 26 … … 54 49 zlay ! level classes (99 if you don't want) 55 50 REAL(wp) :: slopeSection ! coeff directeur de la section 56 INTEGER :: nb_point ! nb de points de la section51 INTEGER :: nb_point ! section's number of point 57 52 TYPE(POINT_SECTION),DIMENSION(nb_point_max) :: listPoint ! list of point in the section 58 53 END TYPE SECTION -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90
r2849 r2858 7 7 !! 8 8 !! 9 !! History: 2011: cbricaud Mercator-Ocean9 !! History: 09/2011: Clement Bricaud ( Mercator-Ocean ) 10 10 !! 11 11 !!============================================================================== … … 34 34 PRINT*,'CREATION OF SECTIONS FOR NEMO diadct.F90 ROUTINE' 35 35 PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 36 37 !----------------------!38 !0. Read arguments !39 !----------------------!40 PRINT*,' '41 PRINT*,'READ ARGUMENTS'42 PRINT*,'--------------'43 44 !check number of arguments and display usage message if wrong45 narg=iargc()46 PRINT*,'narg= ',narg47 IF ( narg /= 4 ) THEN48 PRINT *,' Usage : generate_sections jpidta jpjdta jpizoom jpjzoom '49 STOP50 ENDIF51 52 ! read on line arguments53 CALL getarg(1,cdum) ; READ(cdum,*) jpidta54 CALL getarg(2,cdum) ; READ(cdum,*) jpjdta55 CALL getarg(3,cdum) ; READ(cdum,*) jpizoom56 CALL getarg(4,cdum) ; READ(cdum,*) jpjzoom57 58 PRINT*,'jpidta jpjdta =',jpidta,jpjdta59 PRINT*,'jpizoom jpjzoom=',jpizoom,jpjzoom60 61 !------------------!62 !0. INITIALISATION !63 !------------------!64 PRINT*,' '65 PRINT*,'DOMAIN SIZE'66 PRINT*,'--------------'67 68 !Domain size69 jpiglo = jpidta-jpizoom+1 ; jpjglo = jpjdta-jpjzoom+170 jpi = jpiglo ; jpj = jpjglo71 nlci = jpiglo ; nlcj = jpjglo72 nlei = jpiglo ; nlej = jpjglo73 74 PRINT*,'jpiglo jpjglo = ',jpiglo,jpjglo75 PRINT*,'jpi jpj = ',jpi ,jpj76 PRINT*,'nlci nlcj = ',nlci,nlcj77 36 78 37 !-------------------! … … 108 67 CALL read_coord_mesh 109 68 69 PRINT*,'domain sizes: ' 70 PRINT*,'jpi jpj = ',jpi ,jpj 110 71 PRINT*,'domain boundaries: ' 111 72 PRINT*,' 1 1 ',glamt(1,1),gphit(1,1) … … 135 96 ENDDO 136 97 137 !---------------------- !138 !5. ecriture du fichier!139 !---------------------- !98 !--------------------------------! 99 !5.Write section_ijglobal.diadct ! 100 !--------------------------------! 140 101 CALL write_sections 141 102 -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90
r2849 r2858 14 14 15 15 PUBLIC read_coord_mesh 16 PRIVATE read_ncdf17 16 18 17 CONTAINS … … 24 23 !! ** Purpose : Read a coordinate and a meshmask file in NetCDF format 25 24 !!--------------------------------------------------------------------- 26 !! * Local declarations27 INTEGER ::ji,jj28 29 !!----------------------------------------------------------------------30 25 PRINT*,' ' 31 26 PRINT*,'READ COORDINATES AND MESHMASK' 32 27 PRINT*,'-----------------------------' 33 28 29 ! Get coordinates dimensions 30 CALL getdim(cdfile="coordinates.nc") 34 31 35 !Allocate coordinate and meshmaskarray with domain size32 !Allocate coordinates array with domain size 36 33 ALLOCATE(glamt(jpi,jpj)) ; ALLOCATE(gphit(jpi,jpj)) 37 34 ALLOCATE(glamf(jpi,jpj)) ; ALLOCATE(gphif(jpi,jpj)) 38 ALLOCATE(e1t(jpi,jpj) ) !; ALLOCATE(tmask(jpi,jpj,1))35 ALLOCATE(e1t(jpi,jpj) ) 39 36 40 37 !Read glamt 41 CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamt",ksize=(/jpi,jpj,1,1/), & 42 kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=glamt) 38 CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamt",ksize=(/jpi,jpj,1,1/),ptab=glamt) 43 39 44 40 !Read gphit 45 CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphit",ksize=(/jpi,jpj,1,1/), & 46 kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=gphit) 41 CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphit",ksize=(/jpi,jpj,1,1/),ptab=gphit) 47 42 48 43 !Read glamf 49 CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamf",ksize=(/jpi,jpj,1,1/), & 50 kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=glamf) 44 CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamf",ksize=(/jpi,jpj,1,1/),ptab=glamf) 51 45 52 46 !Read gphif 53 CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphif",ksize=(/jpi,jpj,1,1/), & 54 kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=gphif) 47 CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphif",ksize=(/jpi,jpj,1,1/),ptab=gphif) 55 48 56 49 !Read e1t 57 CALL read_ncdf(cdfile="coordinates.nc",cdvar="e1t",ksize=(/jpi,jpj,1,1/), & 58 kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=e1t) 50 CALL read_ncdf(cdfile="coordinates.nc",cdvar="e1t",ksize=(/jpi,jpj,1,1/),ptab=e1t) 59 51 60 ! tmask(:,:,:)=161 62 !compute mig and mjg arrays63 ALLOCATE(mig(jpi))64 ALLOCATE(mjg(jpj))65 ! local domain indices ==> data domain indices66 DO ji = 1, jpi ; mig(ji) = ji + jpizoom - 1 ; ENDDO67 DO jj = 1, jpj ; mjg(jj) = jj + jpjzoom - 1 ; ENDDO68 69 70 52 END SUBROUTINE read_coord_mesh 71 53 72 SUBROUTINE read_ncdf(cdfile,cdvar,ksize,kstart,kcount,ptab) 54 SUBROUTINE getdim(cdfile) 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE coord_mesh_read *** 57 !! 58 !! ** Purpose : Read a coordinate and a meshmask file in NetCDF format 59 !!---------------------------------------------------------------------- 60 !! * Arguments 61 CHARACTER(*),INTENT(IN):: cdfile 62 63 !! * Local declarations 64 INTEGER :: ncid ! file unit 65 INTEGER :: idims ! number of dimensions 66 INTEGER :: istatus, id_var ! dummy variable 67 CHARACTER(len=30) :: clname ! dimension name 68 INTEGER, ALLOCATABLE,DIMENSION(:) :: ndim ! dimensions value 69 !!---------------------------------------------------------------------- 70 71 !Open file 72 istatus=NF90_OPEN(TRIM(cdfile),nf90_nowrite,ncid) 73 74 IF( istatus/= NF90_NOERR )THEN 75 PRINT*,TRIM(cdfile),' not found.stop ' ; STOP 76 ELSE 77 78 ! read number of dimensions 79 istatus=NF90_INQUIRE(ncid,ndimensions=idims) 80 81 ALLOCATE( ndim(idims) ) 82 83 ! read each dimension 84 PRINT*,' File dimensions: ' 85 DO id_var=1,idims 86 istatus=NF90_Inquire_Dimension(ncid,id_var,clname,ndim(id_var)) 87 PRINT*,' ',id_var,clname,ndim(id_var) 88 ENDDO 89 90 !close 91 istatus=NF90_CLOSE( ncid ) 92 IF( istatus/=NF90_NOERR )THEN 93 PRINT*,'Problem for closing ',TRIM(cdfile);STOP 94 ELSE 95 PRINT*,' close ',TRIM(cdfile),' OK' 96 ENDIF 97 98 ENDIF 99 100 !domain dimensions 101 jpi = ndim(1) 102 jpj = ndim(2) 103 104 DEALLOCATE( ndim ) 105 END SUBROUTINE getdim 106 107 SUBROUTINE read_ncdf(cdfile,cdvar,ksize,ptab,kstart,kcount) 73 108 !!---------------------------------------------------------------------- 74 109 !! *** ROUTINE coord_mesh_read *** … … 81 116 INTEGER,DIMENSION(4),INTENT(IN) :: ksize 82 117 INTEGER,DIMENSION(4),INTENT(IN),OPTIONAL :: kstart,kcount 83 REAL(wp),DIMENSION(ksize(1),ksize(2),ksize(3),ksize(4)),INTENT(OUT) ,OPTIONAL:: ptab118 REAL(wp),DIMENSION(ksize(1),ksize(2),ksize(3),ksize(4)),INTENT(OUT):: ptab 84 119 85 120 !! * Local declarations … … 97 132 istatus=NF90_OPEN(TRIM(cdfile),nf90_nowrite,ncid) 98 133 IF( istatus/= NF90_NOERR )THEN 99 PRINT*,TRIM(cdfile),' not found.stop ' ; stop 100 ELSE 101 PRINT*,' File dimensions: ' 102 DO id_var=1,4 103 istatus=NF90_Inquire_Dimension(ncid,id_var,clname,len) 104 PRINT*,' ',id_var,clname,len 105 ENDDO 134 PRINT*,TRIM(cdfile),' not found.stop ' ; STOP 106 135 ENDIF 107 136 … … 111 140 istatus=NF90_INQ_VARID (ncid,TRIM(cdvar),id_var) 112 141 IF( istatus/=NF90_NOERR )THEN 113 PRINT*,TRIM(cdvar),' not found in ',TRIM(cdfile),'.stop'; stop142 PRINT*,TRIM(cdvar),' not found in ',TRIM(cdfile),'.stop';STOP 114 143 ENDIF 115 144 … … 126 155 127 156 IF( istatus/=NF90_NOERR )THEN 128 PRINT*,'Problem for reading ',TRIM(cdvar),' in ',TRIM(cdfile); stop157 PRINT*,'Problem for reading ',TRIM(cdvar),' in ',TRIM(cdfile); STOP 129 158 ELSE 130 159 PRINT*,' read ',TRIM(cdvar),' OK' -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90
r2849 r2858 52 52 llok=.FALSE. 53 53 clname='section_ijglobal.diadct' 54 !CALL file_open(numdctout,clname,llok)55 54 CALL file_open(numdctout,clname,llok,cdform="UNFORMATTED",cdstatus="REPLACE",cdaction="WRITE") 56 55 … … 77 76 IF( secs(jsec)%nb_point .NE. 0 )THEN 78 77 DO jseg=1,secs(jsec)%nb_point 79 point=POINT_SECTION( mig(secs(jsec)%listPoint(jseg)%I) , mjg(secs(jsec)%listPoint(jseg)%J) ) 80 i1 = mig(secs(jsec)%listPoint(jseg)%I) ; i2 = mjg(secs(jsec)%listPoint(jseg)%J) 81 !WRITE(numdctout)point 78 i1 = secs(jsec)%listPoint(jseg)%I ; i2 = secs(jsec)%listPoint(jseg)%J 82 79 WRITE(numdctout)i1,i2 83 80 ENDDO … … 85 82 ENDIF 86 83 87 88 !---------------------89 !WRITE(500,*)jsec90 !WRITE(500,*)secs(jsec)%name91 !WRITE(500,*)secs(jsec)%llstrpond92 !WRITE(500,*)secs(jsec)%ll_ice_section93 !WRITE(500,*)secs(jsec)%ll_date_line94 !WRITE(500,*)secs(jsec)%coordSec95 !WRITE(500,*)secs(jsec)%nb_class96 !WRITE(500,*)secs(jsec)%zsigi97 !WRITE(500,*)secs(jsec)%zsigp98 !WRITE(500,*)secs(jsec)%zsal99 !WRITE(500,*)secs(jsec)%ztem100 !WRITE(500,*)secs(jsec)%zlay101 !WRITE(500,*)secs(jsec)%slopeSection102 !WRITE(500,*)secs(jsec)%nb_point103 !IF( secs(jsec)%nb_point .NE. 0 )THEN104 ! DO jseg=1,secs(jsec)%nb_point105 ! point=POINT_SECTION( mig(secs(jsec)%listPoint(jseg)%I) , mjg(secs(jsec)%listPoint(jseg)%J) )106 ! WRITE(500,*)point107 ! ENDDO108 ! WRITE(500,*)secs(jsec)%direction(1:secs(jsec)%nb_point)109 !ENDIF110 111 !---------------------112 84 ENDDO !end of loop on sections 113 85
Note: See TracChangeset
for help on using the changeset viewer.