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

Last change on this file since 2877 was 2877, checked in by cbricaud, 9 years ago

coding rules

  • Property svn:executable set to *
File size: 22.7 KB
Line 
1MODULE sections_tools
2   !!=====================================================================
3   !!                       ***  MODULE  sections_tools  ***
4   !!
5   !! History: 2011: cbricaud 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 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 
25CONTAINS
26
27  TYPE(COORD_SECTION) FUNCTION pointToCoordF(p)
28     !!---------------------------------------------------------------------
29     !!         ***  FUNCTION pointToCoordF  ***
30     !!           
31     !! ** Purpose: define a point with geographical coordinates
32     !!
33     !!---------------------------------------------------------------------
34     !! * arguments
35     TYPE(POINT_SECTION), INTENT(IN) ::  p
36
37     !!---------------------------------------------------------------------
38
39     pointToCoordF = COORD_SECTION(glamf(p%I,p%J),gphif(p%I,p%J))
40
41     RETURN
42
43  END FUNCTION pointToCoordF
44
45  REAL(wp) FUNCTION distance(coordA,coordB,ld_date_line)
46     !!---------------------------------------------------------------------
47     !!         ***  FUNCTION distance  ***
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     !!
53     !!---------------------------------------------------------------------
54     !! * arguments
55     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB
56     LOGICAL,INTENT(IN),OPTIONAL     :: ld_date_line
57
58     !! * Local declarations
59     INTEGER ::idateline             ! local integer
60
61     !-----------------------------------------------
62     !
63     idateline=0
64     IF( PRESENT( ld_date_line))THEN
65       IF( ld_date_line )idateline=1
66     ENDIF
67     !
68     !function output
69     distance = (sqrt((coordA%lon-coordB%lon-360*idateline )**2 + (coordA%lat-coordB%lat )**2))
70     !
71     RETURN
72     !
73  END FUNCTION distance
74
75  REAL(wp) FUNCTION distance2(coordA,coordB,ld_date_line)
76     !!---------------------------------------------------------------------
77     !!         ***  FUNCTION distance  ***
78     !!
79     !! ** Purpose:compute distance between coordA and coordB
80     !!            We add 360 to coordB%long if the line (coordA,coordB)
81     !!            crosses the date-line.
82     !!
83     !!---------------------------------------------------------------------
84     !! * arguments
85     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB
86     LOGICAL,INTENT(IN),OPTIONAL     :: ld_date_line
87
88     !! * Local declarations
89     REAL(wp) :: zrad, zrayon
90     REAL(wp) :: zpi = 3.141592653589793
91     REAL(wp) :: zlam1, zlam2, zphi1, zphi2
92
93     !-----------------------------------------------
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) )
103     !
104     RETURN
105  END FUNCTION distance2
106
107  REAL(wp) FUNCTION distance3(coordA,coordB,coordC,ll_debug) 
108     !!---------------------------------------------------------------------
109     !!         ***  FUNCTION cosinus  ***
110     !!
111     !! ** Purpose: compute distance between a point B and a line (AC)
112     !!
113     !! ** Methode: use Al-Kashi theorem
114     !!
115     !!              B                   A first point of the section
116     !!             / \                  B intermediar point on the section
117     !!            / | \                 C last point of the section
118     !!           /  |  \               
119     !!          /   |   \               angle=ACB
120     !!         /    |    \             
121     !!        /     |     \             distance3=Bd
122     !!       /      |      \
123     !!    C  -------d------- A 
124     !!
125     !!---------------------------------------------------------------------
126     !! * arguments
127     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB,coordC
128     LOGICAL                         :: ll_debug
129
130     !! * Local declarations
131     REAL(wp)                        :: ztmp
132     REAL(wp)                        :: za, zb, zc, zangle
133
134     !-----------------------------------------------
135     za = SQRT( (coordB%lon-coordC%lon)**2 + (coordB%lat-coordC%lat)**2 )
136     zb = SQRT( (coordA%lon-coordC%lon)**2 + (coordA%lat-coordC%lat)**2 )
137     zc = SQRT( (coordA%lon-coordB%lon)**2 + (coordA%lat-coordB%lat)**2 )
138     !
139     IF( za /= 0. .AND. zb /= 0. )THEN
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)) 
144        ENDIF
145     ELSE
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
152     RETURN
153  END FUNCTION distance3
154
155  REAL(wp) FUNCTION slope_coeff(coordA,coordB,ld_dateline)
156     !!---------------------------------------------------------------------
157     !!         ***  Function slope_coeff  ***
158     !!
159     !! ** Purpose: Compute slope coefficient of the line (coordA,coordB)
160     !!
161     !! ** Method:
162     !!           Usual method :
163     !!                   slope_coeff = (latB-latA)/(lonB-lonA)
164     !!           Special case: the segment [A,B] crosses the date-line (lddate=T)
165     !!                   slope_coeff = (latB-latA)/(360+lonB-lonA)
166     !!---------------------------------------------------------------------
167     !! * arguments
168     TYPE(COORD_SECTION), INTENT(IN) :: coordA, coordB
169     LOGICAL,INTENT(IN),OPTIONAL     :: ld_dateline
170
171     !! * Local declarations
172     REAL(wp)                        :: zcoeff
173     INTEGER                         :: idateline
174
175     !!---------------------------------------------------------
176     !initialization
177     zcoeff = 1.e20
178     idateline=0
179     IF( PRESENT( ld_dateline))THEN
180       IF( ld_dateline )idateline=1
181     ENDIF
182
183     !compute slope coefficient
184     IF ( coordB%lon-coordA%lon .NE. 0.) &
185         zcoeff=(coordB%lat-coordA%lat)/(360*idateline+coordB%lon-coordA%lon)
186
187     !output
188     slope_coeff = zcoeff
189     RETURN
190
191  END FUNCTION slope_coeff 
192
193  TYPE(COORD_SECTION) FUNCTION intersec(sec,coord_a,coord_b)
194     !!---------------------------------------------------------------------
195     !!         ***  Function intersec  ***
196     !!
197     !! ** Purpose:Return the coordinates of the intersection point,
198     !!            between segment [a,b] and the section.
199     !!            If no intersection the point = -9999.
200     !!
201     !! ** Method:(coord_a,coord_b) => y=za2*x+zb2  (2)
202     !!            sec              => y=za1*x+zb1  (1)
203     !!    Intersection = (X,Y) solves (1) and (2)
204     !!    Vefify that (X,Y) is in [coord_a,coord_b] AND in the section
205     !!         
206     !! ** Action: 1. compute za1, za1
207     !!            2. compute zb1, zb2
208     !!            3. compute X and Y
209     !!            4.Verify that (zX,zY) is in [coord_a,coord_b]
210     !!              and in [sec%coordSec(1),sec%coordSec(2)]   
211     !!
212     !! History: Author: 10/05 Matthieu Laborie
213     !!          Additions:
214     !!           05-2007: (C Bricaud) add special cases
215     !!                                (crossing date line)
216     !-----------------------------------------------------------------------------
217     !! * arguments
218     TYPE(SECTION), INTENT(IN)        :: sec
219     TYPE(COORD_SECTION)  , INTENT(IN):: coord_a, coord_b
220
221     !! * Local declarations
222     TYPE(COORD_SECTION) :: coordInter
223     REAL(wp)            :: za1,za2, &! slope coefficient
224                            zb1,zb2, &!
225                            zX,zY     ! coordinates of intersection
226     LOGICAL             :: ll_date_line=.FALSE. !segment [a,b] crosses the date-line ?
227
228     !----------------------------------------------------------------------------
229
230     !=================!
231     !0. INITALIZATION !
232     !=================!
233     coordInter=COORD_SECTION(-9999.,-9999.) 
234     ll_date_line=.FALSE.   
235     !we need to know if [coord_a,coord_b] crosses the date-line
236     if(coord_a%lon-coord_b%lon .GT. 180)ll_date_line=.TRUE.
237
238     !=======================!
239     !1. compute za1 and za2 !
240     !=======================!
241     za1=sec%slopeSection
242     za2=slope_coeff(coord_a,coord_b,ll_date_line)
243 
244     !=======================!
245     !2. compute zb1 and zb2 !
246     !=======================!
247
248     ! Compute coefficient b for a straight line y=a*x+b
249     !        Usual method:  knowing value of a, we compute b with coordinates of 1 point:
250     !                       b=latA-a.lonA ou bien b=latB-a.lonB
251     !        Particular case: the straight line crosses the date line; so it is in 2 parts:
252     !                         one on the left of the date-line and one the right
253     !                         then,b can could have 2 values:
254     ! As the date-line can be crossed by the section and the segment of the mesh [a,b],
255     ! we have to check if the segment [a,b] crosses the date-line (we know it already for the
256     ! section(sec%ll_date_line): ll_date_line
257     ! Then, there are 4 cases for computing b1 and b2 ( sec%ll_date_line=T/F AND ll_date_line=T/F )
258     !=========================================================================================!
259     ! CASE A :                               !CASE C:
260     ! sec%ll_date_line=T AND ll_date_line =T !sec%ll_date_line=F AND ll_date_line =T
261     !                                        !     |  -180 !+180   |
262     ! this case doesn't exist                !     |       !       |
263     !========================================! a___|_______!_______|__b
264     !CASE B:                                 !     |       !       |
265     !sec%ll_date_line=T AND ll_date_line=F   !     |       !       |
266     !     \     |                            ! section  date-line  section
267     !   a__\__b |                            ! b1: usual method
268     !       \   |                            ! b2: depend of longitude of
269     !        \  |                            !     sec%coordSec(1) et de  sec%coordSec(2)
270     !         \ |                            !=================================================!
271     !          \|                            !CASE D:
272     !     +180  \  -180                      !sec%ll_date_line=F AND ll_date_line =F
273     !           |\                           !      |       !       |
274     !           | \                          !      | +180  ! -180  |
275     !           |  \                         !  a___|___b   !   a___|___b
276     !           |a__\__b                     !      |       !       |
277     !           |    \                       !      |       !       |
278     !    date-line   section                 !   section date-line section
279     !b2: usual method                        !b1: usual method
280     !b1: depend of longitude of a and b.     !b2: usual method
281     !==========================================================================================!
282
283     IF( sec%ll_date_line )THEN
284          IF( ll_date_line )THEN
285              !CASE A: this case doesn't exist
286              !-------
287              zb1=1.e20
288              zb2=1.e20
289          ELSE 
290              !CASE B:
291              !-------
292              !compute zb2:
293              !zb2 = coord_b%lat - za2 * coord_b%lon = coord_a%lat - za2 * coord_a%lon
294              zb2=coord_b%lat - za2 * coord_b%lon 
295
296              !compute zb1:
297              IF( coord_a%lon .GT. 0 .AND. coord_b%lon .GT. 0  )THEN
298                  zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
299              ELSE 
300                  zb1 = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
301              ENDIF
302          ENDIF
303     ELSE
304          IF( ll_date_line )THEN
305              !CASE C:
306              !-------
307              !Compute zb1:
308              !zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
309              !    = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
310              zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
311
312              !Compute zb2:
313              IF( sec%coordSec(1)%lon .GT. 0 .AND. sec%coordSec(2)%lon .GT. 0 )THEN
314                  zb2 = coord_a%lat - za2 * coord_a%lon
315              ELSE 
316                  zb2 = coord_b%lat - za2 * coord_b%lon
317              ENDIF
318          ELSE
319              !CASE D:
320              !-------
321              !compute zb1:
322              !zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon 
323              !    = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
324              zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
325
326              !compute zb2:
327              !zb2 = coord_b%lat - za2 * coord_b%lon = coord_a%lat - za2 * coord_a%lon
328              zb2 = coord_b%lat - za2 * coord_b%lon 
329          ENDIF
330     ENDIF
331
332     IF( (za1 - za2) .NE. 0) THEN
333          !================================!
334          !3. Compute intersection (zX,zY) !
335          !================================!
336          IF( za1 == 1.e20 ) THEN ! Case X=cste
337              zX = sec%coordSec(1)%lon
338              zY = za2 * zX + zb2
339          ELSE IF ( za2 == 1.e20) THEN ! Case X=cste
340              zX = coord_b%lon
341              zY = za1 * zX + zb1
342          ELSE ! Case zY=A*zX+B
343              zX = (zb2 - zb1) / (za1 - za2)
344              zY = (za1 * zX ) + zb1
345          ENDIF
346           
347          !==============================================!
348          !4.Verify that (zX,zY) is in [coord_a,coord_b] !
349          !  and in [sec%coordSec(1),sec%coordSec(2)]    !
350          !==============================================!
351          !Be carreful! The test is not the same for all configurations
352          IF( sec%ll_date_line )THEN
353
354              IF( ll_date_line )THEN
355                  !CASE A: this case doesn't exist
356                  !-------
357                  coordInter=COORD_SECTION(-9999.,-9999.)
358              ELSE
359                  !CASE B:
360                  !-------   
361                  IF( zX .GE. MIN(coord_a%lon,coord_b%lon ) .AND. &
362                      zX .LE. MAX(coord_a%lon,coord_b%lon ) .AND. &
363                      zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
364                      zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
365                    ((zX .LE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
366                      zX .GE. -180) .OR.  &
367                     (zX .GE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
368                      zX .LE. 180)) .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          ELSE
375
376              IF( ll_date_line )THEN
377                  !CASE C:
378                  !-------
379                  IF( ((zX .LE. MIN(coord_a%lon,coord_b%lon ) .AND. &
380                        zX .GE. -180 )  .OR. &
381                       (zX .GE. MAX(coord_a%lon,coord_b%lon ) .AND. &
382                        zX .LE. 180 )) .AND. & 
383                        zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
384                        zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
385                        zX .GE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
386                        zX .LE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
387                        zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
388                        zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
389                        coordInter = COORD_SECTION(zX, zY)   
390                  ENDIF
391              ELSE
392                  !CASE D:
393                  !-------
394                  IF( zX .GE. MIN(coord_a%lon,coord_b%lon ) .AND. &
395                      zX .LE. MAX(coord_a%lon,coord_b%lon ) .AND. &
396                      zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
397                      zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
398                      zX .GE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
399                      zX .LE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
400                      zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
401                      zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
402                      coordInter = COORD_SECTION(zX, zY)
403                  ENDIF             
404              ENDIF
405          ENDIF   
406           
407     ENDIF
408     
409     !output
410     intersec = coordInter
411     RETURN
412   
413  END FUNCTION intersec
414
415  SUBROUTINE qcksrt(arr1,arr2,n)
416     !!---------------------------------------------------------------------
417     !!         ***  SUBROUTINE qcksrt  ***
418     !!
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
426
427     !! * Local declarations
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
434     !---------------------------------------------------
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
511  END SUBROUTINE qcksrt
512
513  SUBROUTINE write_debug(ksec,cd_write)
514     !!---------------------------------------------------------------------
515     !!         ***  SUBROUTINE write_debug  ***
516     !!
517     !! ** Purpose: write debug messages
518     !!
519     !!---------------------------------------------------------------------
520     !! * arguments
521     INTEGER           :: ksec       !number of section
522     CHARACTER(len=*)  :: cd_write   !message to write
523
524     !! * Local declarations
525     INTEGER           :: iunit
526     LOGICAL           :: llok
527     CHARACTER(len=80) :: clfilename
528   
529     !!---------------------------------------------------------------------
530     IF( ksec == nsecdebug .OR. nsecdebug == -1 )THEN
531
532
533        !open / verify is debug output file is open
534
535        clfilename=TRIM(secs(ksec)%name)//"_debug"
536        iunit = num_sec_debug(ksec)
537
538        IF( iunit .EQ. 0 )THEN
539           PRINT*,"Open debug file: "//TRIM(clfilename)
540           iunit = 100 + ksec
541           num_sec_debug(ksec) = iunit 
542 
543           CALL file_open(iunit,clfilename,llok,cdform="FORMATTED",cdstatus="REPLACE",cdaction="WRITE")
544       
545           IF( .NOT. llok )THEN
546               PRINT*,"Can not open TRIM(clfilename)." ; STOP
547           ENDIF   
548        !ELSE
549        !   PRINT*,TRIM(clfilename)//" already opened" 
550        ENDIF
551     
552       !write
553       WRITE(iunit,*)TRIM(cd_write)
554
555     ENDIF
556
557  END SUBROUTINE write_debug
558
559  SUBROUTINE file_open(knum,cdfile,ldok,cdform,cdstatus,cdaction)
560     !!---------------------------------------------------------------------
561     !!         ***  ROUTINE file_open ***
562     !!
563     !! ** Purpose: open a file
564     !!
565     !!---------------------------------------------------------------------
566     !! * arguments
567     INTEGER          ,INTENT(IN)  :: knum       ! file unit number
568     CHARACTER(len=*),INTENT(IN)   :: cdfile     ! file name to open
569     LOGICAL          ,INTENT(OUT) :: ldok       ! =.TRUE. if file exists and is corectly opened
570     CHARACTER(len=*),INTENT(IN)   :: cdform,   &! FORM arguments for OPEN function
571                                      cdstatus, &! STATUS arguments for OPEN function
572                                      cdaction   ! ACTION arguments for OPEN function
573                                     
574     !! * Local declarations
575     INTEGER :: iost
576     LOGICAL :: llbon=.FALSE. !check existence of file
577     !!---------------------------------------------------------------------
578     ldok = .FALSE.                        ! initialization; file not found and not opened
579
580     !check presence of file
581     IF( cdstatus .EQ. "OLD" ) THEN
582        INQUIRE( FILE=cdfile, EXIST=llbon ) ! check presence of namelist
583        IF( llbon )THEN ;  PRINT*,TRIM(cdfile)//' EXISTS'
584        ELSE            ;  PRINT*,TRIM(cdfile)//' NOT EXISTS' ; STOP
585        ENDIF
586     ENDIF
587
588     !open file
589     OPEN(UNIT=knum,FILE=TRIM(cdfile),FORM=cdform,status=cdstatus,action=cdaction,iostat=iost)
590     IF ( iost == 0 )THEN
591          ldok=.TRUE.
592     ELSE
593          PRINT*,TRIM(cdfile)//' bad opening. STOP.'; STOP
594     ENDIF
595
596  END SUBROUTINE file_open
597
598END MODULE sections_tools 
Note: See TracBrowser for help on using the repository browser.