Changeset 2877
- Timestamp:
- 2011-09-29T18:54:19+02:00 (13 years ago)
- Location:
- branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90
r2860 r2877 301 301 DO ji=1,jpi 302 302 coord_t=COORD_SECTION(glamf(ji,jj),gphif(ji,jj)) 303 zdistFirst = distance(coord_t,coordFirst) 304 zdistLast = distance(coord_t,coordLast) 303 !cbr zdistFirst = distance(coord_t,coordFirst) 304 zdistFirst = distance2(coord_t,coordFirst) 305 !cbr zdistLast = distance(coord_t,coordLast) 306 zdistLast = distance2(coord_t,coordLast) 305 307 IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN 306 308 sec%listPoint(1) = POINT_SECTION(ji,jj) … … 391 393 392 394 !initialize distnew value (with distance between section's extremities) 393 zdistnew = distance(coordFirst,coordLast,sec%ll_date_line) 395 !cbr zdistnew = distance(coordFirst,coordLast,sec%ll_date_line) 396 zdistnew = distance2(coordFirst,coordLast,sec%ll_date_line) 394 397 prevPoint = POINT_SECTION(0,0) 395 398 jseg = 1 … … 707 710 708 711 !now compute new slopeSection with ij-coordinates of first and last point 712 sec%slopeSection = 0 ! default value 709 713 IF( sec%nb_point .ne. 0 )THEN 710 714 IF ( sec%listPoint(sec%nb_point)%I .NE. sec%listPoint(1)%I ) THEN -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90
r2858 r2877 22 22 !! 23 23 !! ** Purpose : Read a coordinate and a meshmask file in NetCDF format 24 !! 24 25 !!--------------------------------------------------------------------- 25 26 PRINT*,' ' … … 54 55 SUBROUTINE getdim(cdfile) 55 56 !!---------------------------------------------------------------------- 56 !! *** ROUTINE coord_mesh_read***57 !! *** ROUTINE getdim *** 57 58 !! 58 !! ** Purpose : Read a coordinate and a meshmask file in NetCDF format 59 !! ** Purpose : get dimsensions of a netcdf file 60 !! 59 61 !!---------------------------------------------------------------------- 60 62 !! * Arguments … … 110 112 !! 111 113 !! ** Purpose : Read a coordinate and a meshmask file in NetCDF format 114 !! 112 115 !!---------------------------------------------------------------------- 113 116 !! * Arguments -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readsections.f90
r2849 r2877 22 22 !! *** ROUTINE read_list_sections *** 23 23 !! 24 !! ** Purpose 24 !! ** Purpose: read ascii file 'list_sections.ascii' that contains 25 !! sections description 25 26 !! 26 !! ** Method27 !!28 !! ** Input29 !!30 !! ** Action31 !!32 !! History33 27 !!--------------------------------------------------------------------- 34 28 !! * arguments … … 45 39 CHARACTER(len=9) :: cdstrpond 46 40 CHARACTER(len=110) :: clname,cdsecname,cltmp 47 REAL,DIMENSION(nb_type_class) :: zclass_value48 TYPE(COORD_SECTION) :: coord_point1,coord_point2,coordTemp49 TYPE(COORD_SECTION), DIMENSION(2) ::coord_sec41 REAL,DIMENSION(nb_type_class) :: zclass_value 42 TYPE(COORD_SECTION) :: coord_point1,coord_point2,coordTemp 43 TYPE(COORD_SECTION), DIMENSION(2) :: coord_sec 50 44 !!--------------------------------------------------------------------- 51 45 PRINT*,' ' -
branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90
r2849 r2877 13 13 14 14 !! * Routine accessibility 15 PUBLIC pointToCoordF 16 PUBLIC distance 17 PUBLIC distance2 18 PUBLIC distance3 19 PUBLIC intersec 20 PUBLIC slope_coeff 21 PUBLIC qcksrt 22 PUBLIC write_debug 23 PUBLIC file_open 15 PUBLIC pointToCoordF ! define a point with geographical coordinates 16 PUBLIC distance ! compute distance between 2 points 17 PUBLIC distance2 ! compute distance between 2 points 18 PUBLIC distance3 ! compute distance between a point and a line 19 PUBLIC intersec ! intersection between 2 lines 20 PUBLIC slope_coeff ! slope coefficient of a line 21 PUBLIC qcksrt ! organize a list of point 22 PUBLIC write_debug ! write debug messages 23 PUBLIC file_open ! open a file ( ascii, bin) 24 24 25 25 CONTAINS … … 28 28 !!--------------------------------------------------------------------- 29 29 !! *** FUNCTION pointToCoordF *** 30 !! ** Purpose 30 !! 31 !! ** Purpose: define a point with geographical coordinates 31 32 !! 32 33 !!--------------------------------------------------------------------- 33 34 !! * arguments 34 35 TYPE(POINT_SECTION), INTENT(IN) :: p 35 !!--------------------------------------------------------------------- 36 37 !!--------------------------------------------------------------------- 38 36 39 pointToCoordF = COORD_SECTION(glamf(p%I,p%J),gphif(p%I,p%J)) 40 37 41 RETURN 42 38 43 END FUNCTION pointToCoordF 39 44 … … 41 46 !!--------------------------------------------------------------------- 42 47 !! *** FUNCTION distance *** 43 !! ** Purpose:compute distance between coordA and coordB 44 !! We add 360] to coordB%long if the line (coordA,coordB) 45 !! crosses the date-line. 48 !! 49 !! ** Purpose: Compute distance between coordA and coordB 50 !! We add 360° to coordB%long if the line (coordA,coordB) 51 !! crosses the date-line. 52 !! 46 53 !!--------------------------------------------------------------------- 47 54 !! * arguments 48 55 TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB 49 56 LOGICAL,INTENT(IN),OPTIONAL :: ld_date_line 57 50 58 !! * Local declarations 51 INTEGER ::idateline 59 INTEGER ::idateline ! local integer 60 52 61 !----------------------------------------------- 62 ! 53 63 idateline=0 54 64 IF( PRESENT( ld_date_line))THEN 55 65 IF( ld_date_line )idateline=1 56 66 ENDIF 67 ! 68 !function output 57 69 distance = (sqrt((coordA%lon-coordB%lon-360*idateline )**2 + (coordA%lat-coordB%lat )**2)) 70 ! 58 71 RETURN 72 ! 59 73 END FUNCTION distance 60 74 … … 62 76 !!--------------------------------------------------------------------- 63 77 !! *** FUNCTION distance *** 78 !! 64 79 !! ** Purpose:compute distance between coordA and coordB 65 80 !! We add 360 to coordB%long if the line (coordA,coordB) 66 81 !! crosses the date-line. 82 !! 67 83 !!--------------------------------------------------------------------- 68 84 !! * arguments 69 85 TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB 70 86 LOGICAL,INTENT(IN),OPTIONAL :: ld_date_line 87 71 88 !! * Local declarations 72 REAL(wp) :: rad, zrayon 73 REAL(wp) :: rpi = 3.141592653589793 74 REAL(wp) :: lam1,lam2,phi1,phi2 89 REAL(wp) :: zrad, zrayon 90 REAL(wp) :: zpi = 3.141592653589793 91 REAL(wp) :: zlam1, zlam2, zphi1, zphi2 92 75 93 !----------------------------------------------- 76 zrayon=6367000.0 77 rad=rpi/180. 78 lam1=coordA%lon*rad 79 lam2=coordB%lon*rad 80 phi1=coordA%lat*rad 81 phi2=coordB%lat*rad 82 ! 83 distance2=0.001*2.*zrayon *ASIN(SQRT(SIN((phi2-phi1)/2.0)**2.0+COS(phi1)*COS(phi2)*SIN((lam2-lam1)/2.0)**2.0) ) 94 zrayon = 6367000.0 95 zrad = zpi/180. 96 zlam1 = coordA%lon*zrad 97 zlam2 = coordB%lon*zrad 98 zphi1 = coordA%lat*zrad 99 zphi2 = coordB%lat*zrad 100 ! 101 !function output 102 distance2=0.001*2.*zrayon *ASIN(SQRT(SIN((zphi2-zphi1)/2.0)**2.0+COS(zphi1)*COS(zphi2)*SIN((zlam2-zlam1)/2.0)**2.0) ) 84 103 ! 85 104 RETURN … … 89 108 !!--------------------------------------------------------------------- 90 109 !! *** FUNCTION cosinus *** 91 !! théorème d'Al-Kashi 110 !! 111 !! ** Purpose: compute distance between a point B and a line (AC) 112 !! 113 !! ** Methode: use Al-Kashi theorem 92 114 !! 93 115 !! B A first point of the section … … 105 127 TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB,coordC 106 128 LOGICAL :: ll_debug 129 107 130 !! * Local declarations 108 131 REAL(wp) :: ztmp 109 REAL(wp) :: za,zb,zc,angle 132 REAL(wp) :: za, zb, zc, zangle 133 110 134 !----------------------------------------------- 111 135 za = SQRT( (coordB%lon-coordC%lon)**2 + (coordB%lat-coordC%lat)**2 ) … … 114 138 ! 115 139 IF( za /= 0. .AND. zb /= 0. )THEN 116 ztmp=( za**2 + zb**2 - zc**2 ) / ( 2.0*za*zb ) 117 ztmp=MIN(ztmp,1.00_wp) 118 !angle=ABS(acos( ( za**2 + zb**2 - zc**2 ) / ( 2.0*za*zb ) ) ) 119 IF( ztmp==1.00_wp )THEN 120 angle=0.0 121 ELSE 122 angle=ABS(acos(ztmp)) 140 ztmp = ( za**2 + zb**2 - zc**2 ) / ( 2.0*za*zb ) 141 ztmp = MIN(ztmp,1.00_wp) 142 IF( ztmp==1.00_wp )THEN ; zangle = 0.0 143 ELSE ; zangle = ABS(acos(ztmp)) 123 144 ENDIF 124 145 ELSE 125 PRINT*,'za zb =',za,zb ; STOP 126 ENDIF 127 ! 128 distance3=za*ABS(sin(angle)) 129 ! 146 PRINT*,'You should not have this situation: za zb =',za,zb ; STOP 147 ENDIF 148 ! 149 distance3 = za*ABS(sin(zangle)) 150 ! 151 !function output 130 152 RETURN 131 153 END FUNCTION distance3 … … 146 168 TYPE(COORD_SECTION), INTENT(IN) :: coordA, coordB 147 169 LOGICAL,INTENT(IN),OPTIONAL :: ld_dateline 170 148 171 !! * Local declarations 149 172 REAL(wp) :: zcoeff 150 173 INTEGER :: idateline 174 151 175 !!--------------------------------------------------------- 152 176 !initialization … … 194 218 TYPE(SECTION), INTENT(IN) :: sec 195 219 TYPE(COORD_SECTION) , INTENT(IN):: coord_a, coord_b 220 196 221 !! * Local declarations 197 222 TYPE(COORD_SECTION) :: coordInter … … 200 225 zX,zY ! coordinates of intersection 201 226 LOGICAL :: ll_date_line=.FALSE. !segment [a,b] crosses the date-line ? 227 202 228 !---------------------------------------------------------------------------- 203 229 … … 391 417 !! *** SUBROUTINE qcksrt *** 392 418 !! 393 !! ** Purpose 394 !! 395 !! ** Method 396 !! 397 !! ** Input 398 !! 399 !! ** Action 400 !! 401 !! History 402 !!--------------------------------------------------------------------- 403 !! * arguments 404 INTEGER,INTENT(IN) :: n 405 REAL(wp),DIMENSION(n),INTENT(INOUT):: arr1,arr2 419 !! ** Purpose : organize a list point by latitude , and after by longitude 420 !! 421 !! History 422 !!--------------------------------------------------------------------- 423 !! * arguments 424 INTEGER,INTENT(IN) :: n ! number of points 425 REAL(wp),DIMENSION(n),INTENT(INOUT):: arr1,arr2 ! longitude and latitude 406 426 407 427 !! * Local declarations 408 REAL(wp) :: fx,a,b 409 INTEGER,PARAMETER :: m=7,nstack=500 410 REAL(wp),PARAMETER :: fm=7875.,fa=211.,fc=1663.,fmi=1./fm 411 INTEGER :: l,ir,j,i,iq,jstack 412 INTEGER,DIMENSION(nstack) :: istack(nstack) 428 INTEGER , PARAMETER :: m = 7, nstack = 500 429 REAL(wp), PARAMETER :: fm=7875.,fa=211.,fc=1663.,fmi=1./fm 430 INTEGER :: il, ir, i, iq, ist !local integers 431 INTEGER :: jj !loop indice 432 REAL(wp) :: zfx, za, zb !local real 433 INTEGER,DIMENSION(nstack) :: istack(nstack) !temp array 413 434 !--------------------------------------------------- 414 jstack=0 415 l=1 416 ir=n 417 fx=0. 418 10 if(ir-l .lt.m)then 419 do j=l+1,ir 420 a=arr1(j) 421 b=arr2(j) 422 do i=j-1,1,-1 423 if (arr1(i).le.a) goto 12 424 arr1(i+1)=arr1(i) 425 arr2(i+1)=arr2(i) 426 enddo 427 i=0 428 12 arr1(i+1)=a 429 arr2(i+1)=b 430 enddo 431 if(jstack.eq.0)return 432 ir=istack(jstack) 433 l=istack(jstack-1) 434 jstack=jstack-2 435 else 436 i=l 437 j=ir 438 fx=MOD(fx*fa+fc,fm) 439 iq=l+(ir-l+1)*(fx*fmi) 440 a=arr1(iq) 441 arr1(iq)=arr1(l) 442 b=arr2(iq) 443 arr2(iq)=arr2(l) 444 20 continue 445 21 if (j.gt.0) then 446 if(a.lt.arr1(j))then 447 j=j-1 448 goto 21 449 endif 450 endif 451 if (j.le.i)then 452 arr1(i)=a 453 arr2(i)=b 454 goto 30 455 endif 456 arr1(i)=arr1(j) 457 arr2(i)=arr2(j) 458 i=i+1 459 22 if (i.le.n) then 460 if (a.gt.arr1(i)) then 461 i=i+1 462 goto 22 463 endif 464 endif 465 if (j.le.i)then 466 arr1(j)=a 467 arr2(j)=b 468 i=j 469 goto 30 470 endif 471 arr1(j)=arr1(i) 472 arr2(j)=arr2(i) 473 j=j-1 474 goto 20 475 30 jstack=jstack+2 476 if(jstack.gt.nstack)pause 'nstack trop petit' 477 if (ir-i.ge.i-1) then 478 istack(jstack)=ir 479 istack(jstack-1)=i+1 480 ir=i-1 481 else 482 istack(jstack)=i-1 483 istack(jstack-1)=l 484 l=i+1 485 endif 486 endif 487 goto 10 435 ist=0 436 il=1 437 ir=n 438 zfx=0. 439 ! 440 10 IF( ir-il .LT. m )THEN 441 442 DO jj = il+1,ir 443 za = arr1(jj) 444 zb = arr2(jj) 445 DO i = jj-1,1,-1 446 IF( arr1(i) .LE. za )GOTO 12 447 arr1(i+1) = arr1(i) 448 arr2(i+1) = arr2(i) 449 ENDDO 450 i=0 451 12 arr1(i+1) = za 452 arr2(i+1) = zb 453 ENDDO 454 IF( ist .EQ. 0 )RETURN 455 ir = istack(ist) 456 il = istack(ist-1) 457 ist = ist-2 458 ELSE 459 i = il 460 jj = ir 461 zfx = MOD(zfx*fa+fc,fm) 462 iq = il+(ir-il+1)*(zfx*fmi) 463 za = arr1(iq) 464 arr1(iq) = arr1(il) 465 zb = arr2(iq) 466 arr2(iq) = arr2(il) 467 20 CONTINUE 468 21 IF( jj .GT. 0 )THEN 469 IF( za .LT. arr1(jj) )THEN 470 jj = jj-1 471 GOTO 21 472 ENDIF 473 ENDIF 474 IF( jj .LE. i )THEN 475 arr1(i) = za 476 arr2(i) = zb 477 GOTO 30 478 ENDIF 479 arr1(i) = arr1(jj) 480 arr2(i) = arr2(jj) 481 i=i+1 482 22 IF( i .LE. n )THEN 483 IF( za .GT. arr1(i) )THEN 484 i = i+1 485 GOTO 22 486 ENDIF 487 ENDIF 488 IF( jj .LE.i )THEN 489 arr1(jj) = za 490 arr2(jj) = zb 491 i = jj 492 GOTO 30 493 ENDIF 494 arr1(jj) = arr1(i) 495 arr2(jj) = arr2(i) 496 jj = jj - 1 497 GOTO 20 498 30 ist = ist + 2 499 IF( ist .GT. nstack )PAUSE 'nstack too small' 500 IF( ir-i .GE. i-1 )THEN 501 istack(ist ) = ir 502 istack(ist-1) = i + 1 503 ir = i-1 504 ELSE 505 istack(ist ) = i-1 506 istack(ist-1) = il 507 il = i+1 508 ENDIF 509 ENDIF 510 GOTO 10 488 511 END SUBROUTINE qcksrt 489 512 490 513 SUBROUTINE write_debug(ksec,cd_write) 491 514 !!--------------------------------------------------------------------- 492 !! *** SUBROUTINE qcksrt *** 493 !! 494 !! ** Purpose 495 !! 496 !! ** Method 497 !! 498 !! ** Input 499 !! 500 !! ** Action 501 !! 502 !! History 515 !! *** SUBROUTINE write_debug *** 516 !! 517 !! ** Purpose: write debug messages 518 !! 503 519 !!--------------------------------------------------------------------- 504 520 !! * arguments … … 545 561 !! *** ROUTINE file_open *** 546 562 !! 547 !! ** Purpose 548 !! 549 !! ** Method 550 !! 551 !! ** Input 552 !! 553 !! ** Action 554 !! 555 !! Hisotry 563 !! ** Purpose: open a file 564 !! 556 565 !!--------------------------------------------------------------------- 557 566 !! * arguments
Note: See TracChangeset
for help on using the changeset viewer.