source: branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90 @ 2878

Last change on this file since 2878 was 2878, checked in by cbricaud, 9 years ago
  • Property svn:executable set to *
File size: 21.7 KB
Line 
1MODULE sections_tools
2   !!=====================================================================
3   !!                       ***  MODULE  sections_tools  ***
4   !!
5   !! History: 2011: Clement Bricaud, Mercator-Ocean
6   !!
7   !!=====================================================================
8   !! * Modules used
9   USE declarations
10 
11   IMPLICIT NONE
12   PRIVATE
13
14   !! * Routine accessibility
15   PUBLIC pointToCoordF         ! define a point with geographical coordinates
16   PUBLIC distance2             ! compute distance between 2 points
17   PUBLIC distance3             ! compute distance between a point and a line
18   PUBLIC intersec              ! intersection between 2 lines
19   PUBLIC slope_coeff           ! slope coefficient of a line
20   PUBLIC qcksrt                ! organize a list of point
21   PUBLIC write_debug           ! write debug messages
22   PUBLIC file_open             ! open a file ( ascii, bin)
23 
24CONTAINS
25
26  TYPE(COORD_SECTION) FUNCTION pointToCoordF(p)
27     !!---------------------------------------------------------------------
28     !!         ***  FUNCTION pointToCoordF  ***
29     !!           
30     !! ** Purpose: define a point with geographical coordinates
31     !!
32     !!---------------------------------------------------------------------
33     !! * arguments
34     TYPE(POINT_SECTION), INTENT(IN) ::  p
35
36     !!---------------------------------------------------------------------
37
38     pointToCoordF = COORD_SECTION(glamf(p%I,p%J),gphif(p%I,p%J))
39
40     RETURN
41
42  END FUNCTION pointToCoordF
43
44  REAL(wp) FUNCTION distance2(coordA,coordB,ld_date_line)
45     !!---------------------------------------------------------------------
46     !!         ***  FUNCTION distance  ***
47     !!
48     !! ** Purpose:compute distance between coordA and coordB
49     !!            We add 360 to coordB%long if the line (coordA,coordB)
50     !!            crosses the date-line.
51     !!
52     !!---------------------------------------------------------------------
53     !! * arguments
54     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB
55     LOGICAL,INTENT(IN),OPTIONAL     :: ld_date_line
56
57     !! * Local declarations
58     REAL(wp) :: zrad, zrayon
59     REAL(wp) :: zpi = 3.141592653589793
60     REAL(wp) :: zlam1, zlam2, zphi1, zphi2
61
62     !-----------------------------------------------
63     zrayon = 6367000.0
64     zrad   = zpi/180.
65     zlam1  = coordA%lon*zrad
66     zlam2  = coordB%lon*zrad
67     zphi1  = coordA%lat*zrad
68     zphi2  = coordB%lat*zrad
69     !
70     !function output
71     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) )
72     !
73     RETURN
74  END FUNCTION distance2
75
76  REAL(wp) FUNCTION distance3(coordA,coordB,coordC,ll_debug) 
77     !!---------------------------------------------------------------------
78     !!         ***  FUNCTION cosinus  ***
79     !!
80     !! ** Purpose: compute distance between a point B and a line (AC)
81     !!
82     !! ** Methode: use Al-Kashi theorem
83     !!
84     !!              B                   A first point of the section
85     !!             / \                  B intermediar point on the section
86     !!            / | \                 C last point of the section
87     !!           /  |  \               
88     !!          /   |   \               angle=ACB
89     !!         /    |    \             
90     !!        /     |     \             distance3=Bd
91     !!       /      |      \
92     !!    C  -------d------- A 
93     !!
94     !!---------------------------------------------------------------------
95     !! * arguments
96     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB,coordC
97     LOGICAL                         :: ll_debug
98
99     !! * Local declarations
100     REAL(wp)                        :: ztmp
101     REAL(wp)                        :: za, zb, zc, zangle
102
103     !-----------------------------------------------
104     za = SQRT( (coordB%lon-coordC%lon)**2 + (coordB%lat-coordC%lat)**2 )
105     zb = SQRT( (coordA%lon-coordC%lon)**2 + (coordA%lat-coordC%lat)**2 )
106     zc = SQRT( (coordA%lon-coordB%lon)**2 + (coordA%lat-coordB%lat)**2 )
107     !
108     IF( za /= 0. .AND. zb /= 0. )THEN
109        ztmp = ( za**2 + zb**2  - zc**2 ) / ( 2.0*za*zb )
110        ztmp = MIN(ztmp,1.00_wp)
111        IF( ztmp==1.00_wp )THEN ; zangle = 0.0 
112        ELSE                    ; zangle = ABS(acos(ztmp)) 
113        ENDIF
114     ELSE
115        PRINT*,'You should not have this situation: za zb =',za,zb ; STOP
116     ENDIF
117     !
118     distance3 = za*ABS(sin(zangle)) 
119     !
120     !function output
121     RETURN
122  END FUNCTION distance3
123
124  REAL(wp) FUNCTION slope_coeff(coordA,coordB,ld_dateline)
125     !!---------------------------------------------------------------------
126     !!         ***  Function slope_coeff  ***
127     !!
128     !! ** Purpose: Compute slope coefficient of the line (coordA,coordB)
129     !!
130     !! ** Method:
131     !!           Usual method :
132     !!                   slope_coeff = (latB-latA)/(lonB-lonA)
133     !!           Special case: the segment [A,B] crosses the date-line (lddate=T)
134     !!                   slope_coeff = (latB-latA)/(360+lonB-lonA)
135     !!---------------------------------------------------------------------
136     !! * arguments
137     TYPE(COORD_SECTION), INTENT(IN) :: coordA, coordB
138     LOGICAL,INTENT(IN),OPTIONAL     :: ld_dateline
139
140     !! * Local declarations
141     REAL(wp)                        :: zcoeff
142     INTEGER                         :: idateline
143
144     !!---------------------------------------------------------
145     !initialization
146     zcoeff = 1.e20
147     idateline=0
148     IF( PRESENT( ld_dateline))THEN
149       IF( ld_dateline )idateline=1
150     ENDIF
151
152     !compute slope coefficient
153     IF ( coordB%lon-coordA%lon .NE. 0.) &
154         zcoeff=(coordB%lat-coordA%lat)/(360*idateline+coordB%lon-coordA%lon)
155
156     !output
157     slope_coeff = zcoeff
158     RETURN
159
160  END FUNCTION slope_coeff 
161
162  TYPE(COORD_SECTION) FUNCTION intersec(sec,coord_a,coord_b)
163     !!---------------------------------------------------------------------
164     !!         ***  Function intersec  ***
165     !!
166     !! ** Purpose:Return the coordinates of the intersection point,
167     !!            between segment [a,b] and the section.
168     !!            If no intersection the point = -9999.
169     !!
170     !! ** Method:(coord_a,coord_b) => y=za2*x+zb2  (2)
171     !!            sec              => y=za1*x+zb1  (1)
172     !!    Intersection = (X,Y) solves (1) and (2)
173     !!    Vefify that (X,Y) is in [coord_a,coord_b] AND in the section
174     !!         
175     !! ** Action: 1. compute za1, za1
176     !!            2. compute zb1, zb2
177     !!            3. compute X and Y
178     !!            4.Verify that (zX,zY) is in [coord_a,coord_b]
179     !!              and in [sec%coordSec(1),sec%coordSec(2)]   
180     !!
181     !! History: Author: 10/05 Matthieu Laborie
182     !!          Additions:
183     !!           05-2007: (C Bricaud) add special cases
184     !!                                (crossing date line)
185     !-----------------------------------------------------------------------------
186     !! * arguments
187     TYPE(SECTION), INTENT(IN)        :: sec
188     TYPE(COORD_SECTION)  , INTENT(IN):: coord_a, coord_b
189
190     !! * Local declarations
191     TYPE(COORD_SECTION) :: coordInter
192     REAL(wp)            :: za1,za2, &! slope coefficient
193                            zb1,zb2, &!
194                            zX,zY     ! coordinates of intersection
195     LOGICAL             :: ll_date_line=.FALSE. !segment [a,b] crosses the date-line ?
196
197     !----------------------------------------------------------------------------
198
199     !=================!
200     !0. INITALIZATION !
201     !=================!
202     coordInter=COORD_SECTION(-9999.,-9999.) 
203     ll_date_line=.FALSE.   
204     !we need to know if [coord_a,coord_b] crosses the date-line
205     if(coord_a%lon-coord_b%lon .GT. 180)ll_date_line=.TRUE.
206
207     !=======================!
208     !1. compute za1 and za2 !
209     !=======================!
210     za1=sec%slopeSection
211     za2=slope_coeff(coord_a,coord_b,ll_date_line)
212 
213     !=======================!
214     !2. compute zb1 and zb2 !
215     !=======================!
216
217     ! Compute coefficient b for a straight line y=a*x+b
218     !        Usual method:  knowing value of a, we compute b with coordinates of 1 point:
219     !                       b=latA-a.lonA ou bien b=latB-a.lonB
220     !        Particular case: the straight line crosses the date line; so it is in 2 parts:
221     !                         one on the left of the date-line and one the right
222     !                         then,b can could have 2 values:
223     ! As the date-line can be crossed by the section and the segment of the mesh [a,b],
224     ! we have to check if the segment [a,b] crosses the date-line (we know it already for the
225     ! section(sec%ll_date_line): ll_date_line
226     ! Then, there are 4 cases for computing b1 and b2 ( sec%ll_date_line=T/F AND ll_date_line=T/F )
227     !=========================================================================================!
228     ! CASE A :                               !CASE C:
229     ! sec%ll_date_line=T AND ll_date_line =T !sec%ll_date_line=F AND ll_date_line =T
230     !                                        !     |  -180 !+180   |
231     ! this case doesn't exist                !     |       !       |
232     !========================================! a___|_______!_______|__b
233     !CASE B:                                 !     |       !       |
234     !sec%ll_date_line=T AND ll_date_line=F   !     |       !       |
235     !     \     |                            ! section  date-line  section
236     !   a__\__b |                            ! b1: usual method
237     !       \   |                            ! b2: depend of longitude of
238     !        \  |                            !     sec%coordSec(1) et de  sec%coordSec(2)
239     !         \ |                            !=================================================!
240     !          \|                            !CASE D:
241     !     +180  \  -180                      !sec%ll_date_line=F AND ll_date_line =F
242     !           |\                           !      |       !       |
243     !           | \                          !      | +180  ! -180  |
244     !           |  \                         !  a___|___b   !   a___|___b
245     !           |a__\__b                     !      |       !       |
246     !           |    \                       !      |       !       |
247     !    date-line   section                 !   section date-line section
248     !b2: usual method                        !b1: usual method
249     !b1: depend of longitude of a and b.     !b2: usual method
250     !==========================================================================================!
251
252     IF( sec%ll_date_line )THEN
253          IF( ll_date_line )THEN
254              !CASE A: this case doesn't exist
255              !-------
256              zb1=1.e20
257              zb2=1.e20
258          ELSE 
259              !CASE B:
260              !-------
261              !compute zb2:
262              !zb2 = coord_b%lat - za2 * coord_b%lon = coord_a%lat - za2 * coord_a%lon
263              zb2=coord_b%lat - za2 * coord_b%lon 
264
265              !compute zb1:
266              IF( coord_a%lon .GT. 0 .AND. coord_b%lon .GT. 0  )THEN
267                  zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
268              ELSE 
269                  zb1 = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
270              ENDIF
271          ENDIF
272     ELSE
273          IF( ll_date_line )THEN
274              !CASE C:
275              !-------
276              !Compute zb1:
277              !zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
278              !    = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
279              zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
280
281              !Compute zb2:
282              IF( sec%coordSec(1)%lon .GT. 0 .AND. sec%coordSec(2)%lon .GT. 0 )THEN
283                  zb2 = coord_a%lat - za2 * coord_a%lon
284              ELSE 
285                  zb2 = coord_b%lat - za2 * coord_b%lon
286              ENDIF
287          ELSE
288              !CASE D:
289              !-------
290              !compute zb1:
291              !zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon 
292              !    = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
293              zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
294
295              !compute zb2:
296              !zb2 = coord_b%lat - za2 * coord_b%lon = coord_a%lat - za2 * coord_a%lon
297              zb2 = coord_b%lat - za2 * coord_b%lon 
298          ENDIF
299     ENDIF
300
301     IF( (za1 - za2) .NE. 0) THEN
302          !================================!
303          !3. Compute intersection (zX,zY) !
304          !================================!
305          IF( za1 == 1.e20 ) THEN ! Case X=cste
306              zX = sec%coordSec(1)%lon
307              zY = za2 * zX + zb2
308          ELSE IF ( za2 == 1.e20) THEN ! Case X=cste
309              zX = coord_b%lon
310              zY = za1 * zX + zb1
311          ELSE ! Case zY=A*zX+B
312              zX = (zb2 - zb1) / (za1 - za2)
313              zY = (za1 * zX ) + zb1
314          ENDIF
315           
316          !==============================================!
317          !4.Verify that (zX,zY) is in [coord_a,coord_b] !
318          !  and in [sec%coordSec(1),sec%coordSec(2)]    !
319          !==============================================!
320          !Be carreful! The test is not the same for all configurations
321          IF( sec%ll_date_line )THEN
322
323              IF( ll_date_line )THEN
324                  !CASE A: this case doesn't exist
325                  !-------
326                  coordInter=COORD_SECTION(-9999.,-9999.)
327              ELSE
328                  !CASE B:
329                  !-------   
330                  IF( zX .GE. MIN(coord_a%lon,coord_b%lon ) .AND. &
331                      zX .LE. MAX(coord_a%lon,coord_b%lon ) .AND. &
332                      zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
333                      zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
334                    ((zX .LE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
335                      zX .GE. -180) .OR.  &
336                     (zX .GE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
337                      zX .LE. 180)) .AND. &
338                      zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
339                      zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
340                     coordInter = COORD_SECTION(zX, zY)
341                  ENDIF
342              ENDIF
343          ELSE
344
345              IF( ll_date_line )THEN
346                  !CASE C:
347                  !-------
348                  IF( ((zX .LE. MIN(coord_a%lon,coord_b%lon ) .AND. &
349                        zX .GE. -180 )  .OR. &
350                       (zX .GE. MAX(coord_a%lon,coord_b%lon ) .AND. &
351                        zX .LE. 180 )) .AND. & 
352                        zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
353                        zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
354                        zX .GE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
355                        zX .LE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
356                        zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
357                        zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
358                        coordInter = COORD_SECTION(zX, zY)   
359                  ENDIF
360              ELSE
361                  !CASE D:
362                  !-------
363                  IF( zX .GE. MIN(coord_a%lon,coord_b%lon ) .AND. &
364                      zX .LE. MAX(coord_a%lon,coord_b%lon ) .AND. &
365                      zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
366                      zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
367                      zX .GE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
368                      zX .LE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
369                      zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
370                      zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
371                      coordInter = COORD_SECTION(zX, zY)
372                  ENDIF             
373              ENDIF
374          ENDIF   
375           
376     ENDIF
377     
378     !output
379     intersec = coordInter
380     RETURN
381   
382  END FUNCTION intersec
383
384  SUBROUTINE qcksrt(arr1,arr2,n)
385     !!---------------------------------------------------------------------
386     !!         ***  SUBROUTINE qcksrt  ***
387     !!
388     !! ** Purpose : organize a list point by latitude , and after by longitude
389     !!
390     !!     History
391     !!---------------------------------------------------------------------
392     !! * arguments
393     INTEGER,INTENT(IN)                 :: n           ! number of points
394     REAL(wp),DIMENSION(n),INTENT(INOUT):: arr1,arr2   ! longitude and latitude
395
396     !! * Local declarations
397     INTEGER , PARAMETER                :: m = 7, nstack = 500
398     REAL(wp), PARAMETER                :: fm=7875.,fa=211.,fc=1663.,fmi=1./fm
399     INTEGER                            :: il, ir, i, iq, ist    !local integers
400     INTEGER                            :: jj                    !loop indice
401     REAL(wp)                           :: zfx, za, zb           !local real
402     INTEGER,DIMENSION(nstack)          :: istack(nstack)        !temp array
403     !---------------------------------------------------
404     ist=0
405     il=1
406     ir=n
407     zfx=0.
408     !
409 10  IF( ir-il .LT. m )THEN
410
411        DO jj = il+1,ir
412           za = arr1(jj)
413           zb = arr2(jj)
414           DO i = jj-1,1,-1
415              IF( arr1(i) .LE. za )GOTO 12
416              arr1(i+1) = arr1(i)
417              arr2(i+1) = arr2(i)
418           ENDDO
419           i=0
420 12        arr1(i+1) = za
421           arr2(i+1) = zb
422        ENDDO
423        IF( ist .EQ. 0 )RETURN
424        ir  = istack(ist)
425        il  = istack(ist-1)
426        ist = ist-2
427     ELSE
428        i  = il
429        jj = ir
430        zfx = MOD(zfx*fa+fc,fm)
431        iq = il+(ir-il+1)*(zfx*fmi)
432        za  = arr1(iq)
433        arr1(iq) = arr1(il)
434        zb = arr2(iq)
435        arr2(iq) = arr2(il)
436 20     CONTINUE
437 21     IF( jj .GT. 0 )THEN
438           IF( za .LT. arr1(jj) )THEN
439              jj = jj-1
440              GOTO 21
441           ENDIF
442        ENDIF
443        IF( jj .LE. i )THEN
444           arr1(i) = za
445           arr2(i) = zb
446           GOTO 30
447        ENDIF
448        arr1(i) = arr1(jj)
449        arr2(i) = arr2(jj)
450        i=i+1
451 22     IF( i .LE. n )THEN
452           IF( za .GT. arr1(i) )THEN
453              i = i+1
454              GOTO 22
455           ENDIF
456        ENDIF
457        IF( jj .LE.i )THEN
458           arr1(jj) = za
459           arr2(jj) = zb
460           i        = jj
461           GOTO 30
462        ENDIF
463        arr1(jj) = arr1(i)
464        arr2(jj) = arr2(i)
465        jj = jj - 1
466        GOTO 20
467 30     ist = ist + 2
468        IF( ist .GT. nstack )PAUSE 'nstack too small'
469        IF( ir-i .GE. i-1 )THEN
470           istack(ist  ) = ir
471           istack(ist-1) = i + 1
472           ir = i-1
473        ELSE
474           istack(ist  ) = i-1
475           istack(ist-1) = il
476           il = i+1
477        ENDIF
478     ENDIF
479     GOTO 10
480  END SUBROUTINE qcksrt
481
482  SUBROUTINE write_debug(ksec,cd_write)
483     !!---------------------------------------------------------------------
484     !!         ***  SUBROUTINE write_debug  ***
485     !!
486     !! ** Purpose: write debug messages
487     !!
488     !!---------------------------------------------------------------------
489     !! * arguments
490     INTEGER           :: ksec       !number of section
491     CHARACTER(len=*)  :: cd_write   !message to write
492
493     !! * Local declarations
494     INTEGER           :: iunit
495     LOGICAL           :: llok
496     CHARACTER(len=80) :: clfilename
497   
498     !!---------------------------------------------------------------------
499     IF( ksec == nsecdebug .OR. nsecdebug == -1 )THEN
500
501
502        !open / verify is debug output file is open
503
504        clfilename=TRIM(secs(ksec)%name)//"_debug"
505        iunit = num_sec_debug(ksec)
506
507        IF( iunit .EQ. 0 )THEN
508           PRINT*,"Open debug file: "//TRIM(clfilename)
509           iunit = 100 + ksec
510           num_sec_debug(ksec) = iunit 
511 
512           CALL file_open(iunit,clfilename,llok,cdform="FORMATTED",cdstatus="REPLACE",cdaction="WRITE")
513       
514           IF( .NOT. llok )THEN
515               PRINT*,"Can not open TRIM(clfilename)." ; STOP
516           ENDIF   
517        !ELSE
518        !   PRINT*,TRIM(clfilename)//" already opened" 
519        ENDIF
520     
521       !write
522       WRITE(iunit,*)TRIM(cd_write)
523
524     ENDIF
525
526  END SUBROUTINE write_debug
527
528  SUBROUTINE file_open(knum,cdfile,ldok,cdform,cdstatus,cdaction)
529     !!---------------------------------------------------------------------
530     !!         ***  ROUTINE file_open ***
531     !!
532     !! ** Purpose: open a file
533     !!
534     !!---------------------------------------------------------------------
535     !! * arguments
536     INTEGER          ,INTENT(IN)  :: knum       ! file unit number
537     CHARACTER(len=*),INTENT(IN)   :: cdfile     ! file name to open
538     LOGICAL          ,INTENT(OUT) :: ldok       ! =.TRUE. if file exists and is corectly opened
539     CHARACTER(len=*),INTENT(IN)   :: cdform,   &! FORM arguments for OPEN function
540                                      cdstatus, &! STATUS arguments for OPEN function
541                                      cdaction   ! ACTION arguments for OPEN function
542                                     
543     !! * Local declarations
544     INTEGER :: iost
545     LOGICAL :: llbon=.FALSE. !check existence of file
546     !!---------------------------------------------------------------------
547     ldok = .FALSE.                        ! initialization; file not found and not opened
548
549     !check presence of file
550     IF( cdstatus .EQ. "OLD" ) THEN
551        INQUIRE( FILE=cdfile, EXIST=llbon ) ! check presence of namelist
552        IF( llbon )THEN ;  PRINT*,TRIM(cdfile)//' EXISTS'
553        ELSE            ;  PRINT*,TRIM(cdfile)//' NOT EXISTS' ; STOP
554        ENDIF
555     ENDIF
556
557     !open file
558     OPEN(UNIT=knum,FILE=TRIM(cdfile),FORM=cdform,status=cdstatus,action=cdaction,iostat=iost)
559     IF ( iost == 0 )THEN
560          ldok=.TRUE.
561     ELSE
562          PRINT*,TRIM(cdfile)//' bad opening. STOP.'; STOP
563     ENDIF
564
565  END SUBROUTINE file_open
566
567END MODULE sections_tools 
Note: See TracBrowser for help on using the repository browser.