New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sections_tools.f90 in branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src – NEMO

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

Last change on this file since 2858 was 2849, checked in by cbricaud, 13 years ago

tools to compute sections pathway

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