source: branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90 @ 5476

Last change on this file since 5476 was 2949, checked in by cbricaud, 10 years ago

corrections

  • 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 intermediate 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 or 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=constant
306              zX = sec%coordSec(1)%lon
307              zY = za2 * zX + zb2
308          ELSE IF ( za2 == 1.e20) THEN ! Case X=constant
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 careful! 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     !!              using "numerical recipies " routine.
390     !!---------------------------------------------------------------------
391     !! * arguments
392     INTEGER,INTENT(IN)                 :: n           ! number of points
393     REAL(wp),DIMENSION(n),INTENT(INOUT):: arr1,arr2   ! longitude and latitude
394
395     !! * Local declarations
396     INTEGER , PARAMETER                :: m = 7, nstack = 500
397     REAL(wp), PARAMETER                :: fm=7875.,fa=211.,fc=1663.,fmi=1./fm
398     INTEGER                            :: il, ir, i, iq, ist    !local integers
399     INTEGER                            :: jj                    !loop indices
400     REAL(wp)                           :: zfx, za, zb           !local real
401     INTEGER,DIMENSION(nstack)          :: istack(nstack)        !temp array
402     !---------------------------------------------------
403     ist=0
404     il=1
405     ir=n
406     zfx=0.
407     !
408 10  IF( ir-il .LT. m )THEN
409
410        DO jj = il+1,ir
411           za = arr1(jj)
412           zb = arr2(jj)
413           DO i = jj-1,1,-1
414              IF( arr1(i) .LE. za )GOTO 12
415              arr1(i+1) = arr1(i)
416              arr2(i+1) = arr2(i)
417           ENDDO
418           i=0
419 12        arr1(i+1) = za
420           arr2(i+1) = zb
421        ENDDO
422        IF( ist .EQ. 0 )RETURN
423        ir  = istack(ist)
424        il  = istack(ist-1)
425        ist = ist-2
426     ELSE
427        i  = il
428        jj = ir
429        zfx = MOD(zfx*fa+fc,fm)
430        iq = il+(ir-il+1)*(zfx*fmi)
431        za  = arr1(iq)
432        arr1(iq) = arr1(il)
433        zb = arr2(iq)
434        arr2(iq) = arr2(il)
435 20     CONTINUE
436 21     IF( jj .GT. 0 )THEN
437           IF( za .LT. arr1(jj) )THEN
438              jj = jj-1
439              GOTO 21
440           ENDIF
441        ENDIF
442        IF( jj .LE. i )THEN
443           arr1(i) = za
444           arr2(i) = zb
445           GOTO 30
446        ENDIF
447        arr1(i) = arr1(jj)
448        arr2(i) = arr2(jj)
449        i=i+1
450 22     IF( i .LE. n )THEN
451           IF( za .GT. arr1(i) )THEN
452              i = i+1
453              GOTO 22
454           ENDIF
455        ENDIF
456        IF( jj .LE.i )THEN
457           arr1(jj) = za
458           arr2(jj) = zb
459           i        = jj
460           GOTO 30
461        ENDIF
462        arr1(jj) = arr1(i)
463        arr2(jj) = arr2(i)
464        jj = jj - 1
465        GOTO 20
466 30     ist = ist + 2
467        IF( ist .GT. nstack )PAUSE 'nstack too small'
468        IF( ir-i .GE. i-1 )THEN
469           istack(ist  ) = ir
470           istack(ist-1) = i + 1
471           ir = i-1
472        ELSE
473           istack(ist  ) = i-1
474           istack(ist-1) = il
475           il = i+1
476        ENDIF
477     ENDIF
478     GOTO 10
479  END SUBROUTINE qcksrt
480
481  SUBROUTINE write_debug(ksec,cd_write)
482     !!---------------------------------------------------------------------
483     !!         ***  SUBROUTINE write_debug  ***
484     !!
485     !! ** Purpose: write debug messages
486     !!
487     !!---------------------------------------------------------------------
488     !! * arguments
489     INTEGER           :: ksec       !number of section
490     CHARACTER(len=*)  :: cd_write   !message to write
491
492     !! * Local declarations
493     INTEGER           :: iunit
494     LOGICAL           :: llok
495     CHARACTER(len=80) :: clfilename
496   
497     !!---------------------------------------------------------------------
498     IF( ksec == nsecdebug .OR. nsecdebug == -1 )THEN
499
500
501        !open / verify is debug output file is open
502
503        clfilename=TRIM(secs(ksec)%name)//"_debug"
504        iunit = num_sec_debug(ksec)
505
506        IF( iunit .EQ. 0 )THEN
507           PRINT*,"Open debug file: "//TRIM(clfilename)
508           iunit = 100 + ksec
509           num_sec_debug(ksec) = iunit 
510 
511           CALL file_open(iunit,clfilename,llok,cdform="FORMATTED",cdstatus="REPLACE",cdaction="WRITE")
512       
513           IF( .NOT. llok )THEN
514               PRINT*,"Can not open TRIM(clfilename)." ; STOP
515           ENDIF   
516        ENDIF
517     
518       WRITE(iunit,*)TRIM(cd_write)
519
520     ENDIF
521
522  END SUBROUTINE write_debug
523
524  SUBROUTINE file_open(knum,cdfile,ldok,cdform,cdstatus,cdaction)
525     !!---------------------------------------------------------------------
526     !!         ***  ROUTINE file_open ***
527     !!
528     !! ** Purpose: open a file
529     !!
530     !!---------------------------------------------------------------------
531     !! * arguments
532     INTEGER          ,INTENT(IN)  :: knum       ! file unit number
533     CHARACTER(len=*),INTENT(IN)   :: cdfile     ! file name to open
534     LOGICAL          ,INTENT(OUT) :: ldok       ! =.TRUE. if file exists and is corectly opened
535     CHARACTER(len=*),INTENT(IN)   :: cdform,   &! FORM arguments for OPEN function
536                                      cdstatus, &! STATUS arguments for OPEN function
537                                      cdaction   ! ACTION arguments for OPEN function
538                                     
539     !! * Local declarations
540     INTEGER :: iost
541     LOGICAL :: llbon=.FALSE. !check existence of file
542     !!---------------------------------------------------------------------
543     ldok = .FALSE.                        ! initialization; file not found and not opened
544
545     !check presence of file
546     IF( cdstatus .EQ. "OLD" ) THEN
547        INQUIRE( FILE=cdfile, EXIST=llbon ) ! check presence of namelist
548        IF( llbon )THEN ;  PRINT*,TRIM(cdfile)//' EXISTS'
549        ELSE            ;  PRINT*,TRIM(cdfile)//' NOT EXISTS' ; STOP
550        ENDIF
551     ENDIF
552
553     !open file
554     OPEN(UNIT=knum,FILE=TRIM(cdfile),FORM=cdform,status=cdstatus,action=cdaction,iostat=iost)
555     IF ( iost == 0 )THEN
556          ldok=.TRUE.
557     ELSE
558          PRINT*,TRIM(cdfile)//' bad opening. STOP.'; STOP
559     ENDIF
560
561  END SUBROUTINE file_open
562
563END MODULE sections_tools 
Note: See TracBrowser for help on using the repository browser.