Changeset 2877


Ignore:
Timestamp:
2011-09-29T18:54:19+02:00 (9 years ago)
Author:
cbricaud
Message:

coding rules

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  
    301301           DO ji=1,jpi   
    302302              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) 
    305307              IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN 
    306308                 sec%listPoint(1) = POINT_SECTION(ji,jj) 
     
    391393 
    392394        !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)  
    394397        prevPoint  = POINT_SECTION(0,0) 
    395398        jseg       = 1 
     
    707710 
    708711     !now compute new slopeSection with ij-coordinates of first and last point  
     712     sec%slopeSection = 0  ! default value 
    709713     IF( sec%nb_point .ne. 0 )THEN 
    710714        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  
    2222  !! 
    2323  !! ** Purpose :   Read a coordinate and a meshmask file in NetCDF format  
     24  !! 
    2425  !!---------------------------------------------------------------------       
    2526  PRINT*,'              ' 
     
    5455  SUBROUTINE getdim(cdfile) 
    5556  !!---------------------------------------------------------------------- 
    56   !!              ***  ROUTINE coord_mesh_read  *** 
     57  !!              ***  ROUTINE getdim  *** 
    5758  !! 
    58   !! ** Purpose :   Read a coordinate and a meshmask file in NetCDF format 
     59  !! ** Purpose :   get dimsensions of a netcdf file 
     60  !! 
    5961  !!---------------------------------------------------------------------- 
    6062  !! * Arguments 
     
    110112  !! 
    111113  !! ** Purpose :   Read a coordinate and a meshmask file in NetCDF format 
     114  !! 
    112115  !!---------------------------------------------------------------------- 
    113116  !! * Arguments 
  • branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readsections.f90

    r2849 r2877  
    2222     !!         ***  ROUTINE read_list_sections *** 
    2323     !! 
    24      !! ** Purpose 
     24     !! ** Purpose: read ascii file 'list_sections.ascii' that contains 
     25     !!             sections description 
    2526     !! 
    26      !! ** Method 
    27      !! 
    28      !! ** Input 
    29      !! 
    30      !! ** Action 
    31      !! 
    32      !! History 
    3327     !!--------------------------------------------------------------------- 
    3428     !! * arguments 
     
    4539     CHARACTER(len=9)    :: cdstrpond 
    4640     CHARACTER(len=110)  :: clname,cdsecname,cltmp 
    47      REAL,DIMENSION(nb_type_class):: zclass_value 
    48      TYPE(COORD_SECTION) :: coord_point1,coord_point2,coordTemp 
    49      TYPE(COORD_SECTION), DIMENSION(2)::coord_sec 
     41     REAL,DIMENSION(nb_type_class)     :: zclass_value 
     42     TYPE(COORD_SECTION)               :: coord_point1,coord_point2,coordTemp 
     43     TYPE(COORD_SECTION), DIMENSION(2) :: coord_sec 
    5044     !!--------------------------------------------------------------------- 
    5145     PRINT*,'              ' 
  • branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90

    r2849 r2877  
    1313 
    1414   !! * 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) 
    2424  
    2525CONTAINS 
     
    2828     !!--------------------------------------------------------------------- 
    2929     !!         ***  FUNCTION pointToCoordF  *** 
    30      !! ** Purpose 
     30     !!            
     31     !! ** Purpose: define a point with geographical coordinates 
    3132     !! 
    3233     !!--------------------------------------------------------------------- 
    3334     !! * arguments 
    3435     TYPE(POINT_SECTION), INTENT(IN) ::  p 
    35      !!--------------------------------------------------------------------- 
     36 
     37     !!--------------------------------------------------------------------- 
     38 
    3639     pointToCoordF = COORD_SECTION(glamf(p%I,p%J),gphif(p%I,p%J)) 
     40 
    3741     RETURN 
     42 
    3843  END FUNCTION pointToCoordF 
    3944 
     
    4146     !!--------------------------------------------------------------------- 
    4247     !!         ***  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     !! 
    4653     !!--------------------------------------------------------------------- 
    4754     !! * arguments 
    4855     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB 
    4956     LOGICAL,INTENT(IN),OPTIONAL     :: ld_date_line 
     57 
    5058     !! * Local declarations  
    51      INTEGER ::idateline     
     59     INTEGER ::idateline             ! local integer 
     60 
    5261     !----------------------------------------------- 
     62     ! 
    5363     idateline=0 
    5464     IF( PRESENT( ld_date_line))THEN 
    5565       IF( ld_date_line )idateline=1 
    5666     ENDIF 
     67     ! 
     68     !function output 
    5769     distance = (sqrt((coordA%lon-coordB%lon-360*idateline )**2 + (coordA%lat-coordB%lat )**2)) 
     70     ! 
    5871     RETURN 
     72     ! 
    5973  END FUNCTION distance 
    6074 
     
    6276     !!--------------------------------------------------------------------- 
    6377     !!         ***  FUNCTION distance  *** 
     78     !! 
    6479     !! ** Purpose:compute distance between coordA and coordB 
    6580     !!            We add 360 to coordB%long if the line (coordA,coordB) 
    6681     !!            crosses the date-line. 
     82     !! 
    6783     !!--------------------------------------------------------------------- 
    6884     !! * arguments 
    6985     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB 
    7086     LOGICAL,INTENT(IN),OPTIONAL     :: ld_date_line 
     87 
    7188     !! * 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 
    7593     !----------------------------------------------- 
    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) ) 
    84103     ! 
    85104     RETURN 
     
    89108     !!--------------------------------------------------------------------- 
    90109     !!         ***  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 
    92114     !! 
    93115     !!              B                   A first point of the section 
     
    105127     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB,coordC 
    106128     LOGICAL                         :: ll_debug 
     129 
    107130     !! * Local declarations 
    108131     REAL(wp)                        :: ztmp 
    109      REAL(wp)                        :: za,zb,zc,angle 
     132     REAL(wp)                        :: za, zb, zc, zangle 
     133 
    110134     !----------------------------------------------- 
    111135     za = SQRT( (coordB%lon-coordC%lon)**2 + (coordB%lat-coordC%lat)**2 ) 
     
    114138     ! 
    115139     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))  
    123144        ENDIF 
    124145     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 
    130152     RETURN 
    131153  END FUNCTION distance3 
     
    146168     TYPE(COORD_SECTION), INTENT(IN) :: coordA, coordB 
    147169     LOGICAL,INTENT(IN),OPTIONAL     :: ld_dateline 
     170 
    148171     !! * Local declarations 
    149172     REAL(wp)                        :: zcoeff 
    150173     INTEGER                         :: idateline 
     174 
    151175     !!--------------------------------------------------------- 
    152176     !initialization 
     
    194218     TYPE(SECTION), INTENT(IN)        :: sec 
    195219     TYPE(COORD_SECTION)  , INTENT(IN):: coord_a, coord_b 
     220 
    196221     !! * Local declarations 
    197222     TYPE(COORD_SECTION) :: coordInter 
     
    200225                            zX,zY     ! coordinates of intersection 
    201226     LOGICAL             :: ll_date_line=.FALSE. !segment [a,b] crosses the date-line ? 
     227 
    202228     !---------------------------------------------------------------------------- 
    203229 
     
    391417     !!         ***  SUBROUTINE qcksrt  *** 
    392418     !! 
    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 
    406426 
    407427     !! * 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 
    413434     !--------------------------------------------------- 
    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 
    488511  END SUBROUTINE qcksrt 
    489512 
    490513  SUBROUTINE write_debug(ksec,cd_write) 
    491514     !!--------------------------------------------------------------------- 
    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     !! 
    503519     !!--------------------------------------------------------------------- 
    504520     !! * arguments 
     
    545561     !!         ***  ROUTINE file_open *** 
    546562     !! 
    547      !! ** Purpose 
    548      !! 
    549      !! ** Method 
    550      !! 
    551      !! ** Input 
    552      !! 
    553      !! ** Action 
    554      !! 
    555      !! Hisotry 
     563     !! ** Purpose: open a file 
     564     !! 
    556565     !!--------------------------------------------------------------------- 
    557566     !! * arguments 
Note: See TracChangeset for help on using the changeset viewer.