Merge branch 'devel' of mmka.chem.univ.gda.pl:unres into devel
[unres.git] / source / wham / src-M / readrtns.F
1       subroutine read_general_data(*)
2       implicit none
3       include "DIMENSIONS"
4       include "DIMENSIONS.ZSCOPT"
5       include "DIMENSIONS.FREE"
6       include "COMMON.TORSION"
7       include "COMMON.INTERACT"
8       include "COMMON.IOUNITS"
9       include "COMMON.TIME1"
10       include "COMMON.PROT"
11       include "COMMON.PROTFILES"
12       include "COMMON.CHAIN"
13       include "COMMON.NAMES"
14       include "COMMON.FFIELD"
15       include "COMMON.ENEPS"
16       include "COMMON.WEIGHTS"
17       include "COMMON.FREE"
18       include "COMMON.CONTROL"
19       include "COMMON.ENERGIES"
20       character*800 controlcard
21       integer i,j,k,ii,n_ene_found
22       integer ind,itype1,itype2,itypf,itypsc,itypp
23       integer ilen
24       external ilen
25       character*16 ucase
26       character*16 key
27       external ucase
28
29       call card_concat(controlcard,.true.)
30       call readi(controlcard,"N_ENE",n_ene,max_ene)
31       if (n_ene.gt.max_ene) then
32         write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
33      &    max_ene
34         return1
35       endif
36       call readi(controlcard,"NPARMSET",nparmset,1)
37       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
38       call readi(controlcard,"IPARMPRINT",iparmprint,1)
39       write (iout,*) "PARMPRINT",iparmprint
40       if (nparmset.gt.max_parm) then
41         write (iout,*) "Error: parameter out of range: NPARMSET",
42      &    nparmset, Max_Parm
43         return1
44       endif
45       call readi(controlcard,"MAXIT",maxit,5000)
46       call reada(controlcard,"FIMIN",fimin,1.0d-3)
47       call readi(controlcard,"ENSEMBLES",ensembles,0)
48       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
49       write (iout,*) "Number of energy parameter sets",nparmset
50       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
51       write (iout,*) "MaxSlice",MaxSlice
52       call readi(controlcard,"NSLICE",nslice,1)
53       call flush(iout)
54       if (nslice.gt.MaxSlice) then
55         write (iout,*) "Error: parameter out of range: NSLICE",nslice,
56      &    MaxSlice
57         return1
58       endif
59       write (iout,*) "Frequency of storing conformations",
60      & (isampl(i),i=1,nparmset)
61       write (iout,*) "Maxit",maxit," Fimin",fimin
62       call readi(controlcard,"NQ",nQ,1)
63       if (nQ.gt.MaxQ) then
64         write (iout,*) "Error: parameter out of range: NQ",nq,
65      &    maxq
66         return1
67       endif
68       indpdb=0
69       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
70       call reada(controlcard,"DELTA",delta,1.0d-2)
71       call readi(controlcard,"EINICHECK",einicheck,2)
72       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
73       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
74       call readi(controlcard,"RESCALE",rescale_mode,1)
75       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
76       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
77       call readi(controlcard,'SYM',symetr,1)
78       write (iout,*) "DISTCHAINMAX",distchainmax
79       write (iout,*) "delta",delta
80       write (iout,*) "einicheck",einicheck
81       write (iout,*) "rescale_mode",rescale_mode
82       call flush(iout)
83       bxfile=index(controlcard,"BXFILE").gt.0
84       cxfile=index(controlcard,"CXFILE").gt.0
85       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
86      & bxfile=.true.
87       histfile=index(controlcard,"HISTFILE").gt.0
88       histout=index(controlcard,"HISTOUT").gt.0
89       entfile=index(controlcard,"ENTFILE").gt.0
90       zscfile=index(controlcard,"ZSCFILE").gt.0
91       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
92       write (iout,*) "with_dihed_constr ",with_dihed_constr
93       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
94       return
95       end
96 c------------------------------------------------------------------------------
97       subroutine read_efree(*)
98 C
99 C Read molecular data
100 C
101       implicit none
102       include 'DIMENSIONS'
103       include 'DIMENSIONS.ZSCOPT'
104       include 'DIMENSIONS.COMPAR'
105       include 'DIMENSIONS.FREE'
106       include 'COMMON.IOUNITS'
107       include 'COMMON.TIME1'
108       include 'COMMON.SBRIDGE'
109       include 'COMMON.CONTROL'
110       include 'COMMON.CHAIN'
111       include 'COMMON.HEADER'
112       include 'COMMON.GEO'
113       include 'COMMON.FREE'
114       character*320 controlcard,ucase
115       integer iparm,ib,i,j,npars
116       integer ilen
117       external ilen
118      
119       if (hamil_rep) then
120         npars=1
121       else
122         npars=nParmSet
123       endif
124
125       do iparm=1,npars
126
127       call card_concat(controlcard,.true.)
128       call readi(controlcard,'NT',nT_h(iparm),1)
129       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
130       call flush(iout)
131       if (nT_h(iparm).gt.MaxT_h) then
132         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),
133      &    MaxT_h
134         return1
135       endif
136       replica(iparm)=index(controlcard,"REPLICA").gt.0
137       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
138       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
139       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
140      &  replica(iparm)," umbrella ",umbrella(iparm),
141      &  " read_iset",read_iset(iparm)
142       call flush(iout)
143       do ib=1,nT_h(iparm)
144         call card_concat(controlcard,.true.)
145         call readi(controlcard,'NR',nR(ib,iparm),1)
146         if (umbrella(iparm)) then
147           nRR(ib,iparm)=1
148         else
149           nRR(ib,iparm)=nR(ib,iparm)
150         endif
151         if (nR(ib,iparm).gt.MaxR) then
152           write (iout,*)  "Error: parameter out of range: NR",
153      &      nR(ib,iparm),MaxR
154         return1
155         endif
156         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
157         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
158         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
159      &    0.0d0)
160         do i=1,nR(ib,iparm)
161           call card_concat(controlcard,.true.)
162           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
163      &      100.0d0)
164           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
165      &      0.0d0)
166         enddo
167       enddo
168       do ib=1,nT_h(iparm)
169         write (iout,*) "ib",ib," beta_h",
170      &    1.0d0/(0.001987*beta_h(ib,iparm))
171         write (iout,*) "nR",nR(ib,iparm)
172         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
173         do i=1,nR(ib,iparm)
174           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
175      &      "q0",(q0(j,i,ib,iparm),j=1,nQ)
176         enddo
177         call flush(iout)
178       enddo
179
180       enddo
181
182       if (hamil_rep) then
183
184        do iparm=2,nParmSet
185           nT_h(iparm)=nT_h(1)
186          do ib=1,nT_h(iparm)
187            nR(ib,iparm)=nR(ib,1)
188            if (umbrella(iparm)) then
189              nRR(ib,iparm)=1
190            else
191              nRR(ib,iparm)=nR(ib,1)
192            endif
193            beta_h(ib,iparm)=beta_h(ib,1)
194            do i=1,nR(ib,iparm)
195              f(i,ib,iparm)=f(i,ib,1)
196              do j=1,nQ
197                KH(j,i,ib,iparm)=KH(j,i,ib,1) 
198                Q0(j,i,ib,iparm)=Q0(j,i,ib,1) 
199              enddo
200            enddo
201            replica(iparm)=replica(1)
202            umbrella(iparm)=umbrella(1)
203            read_iset(iparm)=read_iset(1)
204          enddo
205        enddo
206         
207       endif
208
209       return
210       end
211 c-----------------------------------------------------------------------------
212       subroutine read_protein_data(*)
213       implicit none
214       include "DIMENSIONS"
215       include "DIMENSIONS.ZSCOPT"
216       include "DIMENSIONS.FREE"
217 #ifdef MPI
218       include "mpif.h"
219       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
220       include "COMMON.MPI"
221 #endif
222       include "COMMON.CHAIN"
223       include "COMMON.IOUNITS"
224       include "COMMON.PROT"
225       include "COMMON.PROTFILES"
226       include "COMMON.NAMES"
227       include "COMMON.FREE"
228       include "COMMON.OBCINKA"
229       character*64 nazwa
230       character*16000 controlcard
231       integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
232       external ilen,iroof
233       if (hamil_rep) then
234         npars=1
235       else
236         npars=nparmset
237       endif
238
239       do iparm=1,npars
240
241 C Read names of files with conformation data.
242       if (replica(iparm)) then
243         nthr = 1
244       else
245         nthr = nT_h(iparm)
246       endif
247       do ib=1,nthr
248       do ii=1,nRR(ib,iparm)
249       write (iout,*) "Parameter set",iparm," temperature",ib,
250      & " window",ii
251       call flush(iout)
252       call card_concat(controlcard,.true.) 
253       write (iout,*) controlcard(:ilen(controlcard))
254       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
255       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
256       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
257       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
258       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
259      & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
260       call reada(controlcard,"TIME_START",
261      &  time_start_collect(ii,ib,iparm),0.0d0)
262       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
263      &  1.0d10)
264       write (iout,*) "rec_start",rec_start(ii,ib,iparm),
265      & " rec_end",rec_end(ii,ib,iparm)
266       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
267      & " time_end",time_end_collect(ii,ib,iparm)
268       call flush(iout)
269       if (replica(iparm)) then
270         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
271         write (iout,*) "Number of trajectories",totraj(ii,iparm)
272         call flush(iout)
273       endif
274       if (nfile_bin(ii,ib,iparm).lt.2 
275      &    .and. nfile_asc(ii,ib,iparm).eq.0
276      &    .and. nfile_cx(ii,ib,iparm).eq.0) then
277         write (iout,*) "Error - no action specified!"
278         return1
279       endif
280       if (nfile_bin(ii,ib,iparm).gt.0) then
281         call card_concat(controlcard,.false.)
282         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
283      &   maxfile_prot,nfile_bin(ii,ib,iparm))
284 #ifdef DEBUG
285         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
286         write(iout,*) (protfiles(i,1,ii,ib,iparm),
287      &    i=1,nfile_bin(ii,ib,iparm))
288 #endif
289       endif
290       if (nfile_asc(ii,ib,iparm).gt.0) then
291         call card_concat(controlcard,.false.)
292         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
293      &   maxfile_prot,nfile_asc(ii,ib,iparm))
294 #ifdef DEBUG
295         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
296         write(iout,*) (protfiles(i,2,ii,ib,iparm),
297      &    i=1,nfile_asc(ii,ib,iparm))
298 #endif
299       else if (nfile_cx(ii,ib,iparm).gt.0) then
300         call card_concat(controlcard,.false.)
301         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
302      &   maxfile_prot,nfile_cx(ii,ib,iparm))
303 #ifdef DEBUG
304         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
305         write(iout,*) (protfiles(i,2,ii,ib,iparm),
306      &    i=1,nfile_cx(ii,ib,iparm))
307 #endif
308       endif
309       call flush(iout)
310       enddo
311       enddo
312
313       enddo
314
315       return
316       end
317 c-------------------------------------------------------------------------------
318       subroutine opentmp(islice,iunit,bprotfile_temp)
319       implicit none
320       include "DIMENSIONS"
321       include "DIMENSIONS.ZSCOPT"
322       include "DIMENSIONS.FREE"
323 #ifdef MPI
324       include "mpif.h"
325       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
326       include "COMMON.MPI"
327 #endif
328       include "COMMON.IOUNITS"
329       include "COMMON.PROTFILES"
330       include "COMMON.PROT"
331       include "COMMON.FREE"
332       character*64 bprotfile_temp
333       character*3 liczba,liczba2
334       character*2 liczba1
335       integer iunit,islice
336       integer ilen,iroof
337       external ilen,iroof
338       logical lerr
339
340       write (liczba1,'(bz,i2.2)') islice
341       write (liczba,'(bz,i3.3)') me
342 #ifdef MPI
343 c      write (iout,*) "separate_parset ",separate_parset,
344 c     &  " myparm",myparm
345       if (separate_parset) then
346       write (liczba2,'(bz,i3.3)') myparm
347       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
348      &  prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
349       open (iunit,file=bprotfile_temp,status="unknown",
350      &    form="unformatted",access="direct",recl=lenrec)
351       else
352       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
353      &  prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
354       open (iunit,file=bprotfile_temp,status="unknown",
355      &    form="unformatted",access="direct",recl=lenrec)
356       endif
357 #else
358       bprotfile_temp = scratchdir(:ilen(scratchdir))//
359      &  "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
360       open (iunit,file=bprotfile_temp,status="unknown",
361      &    form="unformatted",access="direct",recl=lenrec)
362 #endif      
363 c      write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
364 c     &  bprotfile_temp
365 c      call flush(iout)
366       return
367       end
368 c-------------------------------------------------------------------------------
369       subroutine read_database(*)
370       implicit none
371       include "DIMENSIONS"
372       include "DIMENSIONS.ZSCOPT"
373       include "DIMENSIONS.FREE"
374 #ifdef MPI
375       include "mpif.h"
376       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
377       include "COMMON.MPI"
378 #endif
379       include "COMMON.CHAIN"
380       include "COMMON.IOUNITS"
381       include "COMMON.PROTFILES"
382       include "COMMON.NAMES"
383       include "COMMON.VAR"
384       include "COMMON.GEO"
385       include "COMMON.ENEPS"
386       include "COMMON.PROT"
387       include "COMMON.INTERACT"
388       include "COMMON.FREE"
389       include "COMMON.SBRIDGE"
390       include "COMMON.OBCINKA"
391       real*4 csingle(3,maxres2)
392       character*64 nazwa,bprotfile_temp
393       character*3 liczba
394       character*2 liczba1
395       integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
396      &  ll(maxslice),mm(maxslice),if
397       integer nrec,nlines,iscor,iunit,islice
398       double precision energ
399       integer ilen,iroof
400       external ilen,iroof
401       double precision rmsdev,energia(0:max_ene),efree,eini,temp
402       double precision prop(maxQ)
403       integer ntot_all(maxslice,0:maxprocs-1)
404       integer iparm,ib,iib,ir,nprop,nthr,npars
405       double precision etot,time
406       integer ixdrf,iret 
407       logical lerr,linit
408
409       lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
410       lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
411       lenrec=lenrec2+8
412       write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
413      &  " lenrec2",lenrec2
414
415       do i=1,nQ
416         prop(i)=0.0d0
417       enddo
418       do islice=1,nslice
419         ll(islice)=0
420         mm(islice)=0
421       enddo
422       write (iout,*) "nparmset",nparmset
423       if (hamil_rep) then
424         npars=1
425       else
426         npars=nparmset
427       endif
428       do iparm=1,npars
429
430       if (replica(iparm)) then
431         nthr = 1
432       else
433         nthr = nT_h(iparm)
434       endif
435
436       do ib=1,nthr
437       do iR=1,nRR(ib,iparm)
438
439       write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
440       do islice=1,nslice
441         jj(islice)=0
442         kk(islice)=0
443       enddo
444
445       IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
446 c Read conformations from binary DA files (one per batch) and write them to 
447 c a binary DA scratchfile.
448         write (liczba,'(bz,i3.3)') me
449         do if=1,nfile_bin(iR,ib,iparm)
450           nazwa=protfiles(if,1,iR,ib,iparm)
451      &     (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
452           open (ientin,file=nazwa,status="old",form="unformatted",
453      &     access="direct",recl=lenrec2,err=1111)
454           ii=0
455           do islice=1,nslice
456             call opentmp(islice,ientout,bprotfile_temp)
457             call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
458      &        mm(islice),iR,ib,iparm)
459             close(ientout)
460           enddo
461           close(ientin)
462         enddo
463       ENDIF ! NFILE_BIN>0
464 c
465       IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
466 c Read conformations from multiple ASCII int files and write them to a binary
467 c DA scratchfile.
468         do if=1,nfile_asc(iR,ib,iparm)
469           nazwa=protfiles(if,2,iR,ib,iparm)
470      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
471           open(unit=ientin,file=nazwa,status='old',err=1111)
472           write(iout,*) "reading ",nazwa(:ilen(nazwa))
473           ii=0
474           call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
475         enddo ! if
476       ENDIF
477       IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
478 c Read conformations from cx files and write them to a binary
479 c DA scratchfile.
480         do if=1,nfile_cx(iR,ib,iparm)
481           nazwa=protfiles(if,2,iR,ib,iparm)
482      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
483           write(iout,*) "reading ",nazwa(:ilen(nazwa))
484           ii=0
485           print *,"Calling cxread"
486           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
487      &       *1111)
488           close(ientout)
489           write (iout,*) "exit cxread"
490           call flush(iout)
491         enddo
492       ENDIF
493
494       do islice=1,nslice
495         stot(islice)=stot(islice)+jj(islice)
496       enddo
497
498       enddo
499       enddo
500       write (iout,*) "IPARM",iparm
501       enddo
502
503       if (nslice.eq.1) then
504 #ifdef MPI
505         write (liczba,'(bz,i3.3)') me
506         bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
507      &    prefix(:ilen(prefix))//liczba//".xbin.tmp"
508 #else
509         bprotfile_temp = scratchdir(:ilen(scratchdir))//
510      &     "/"//prefix(:ilen(prefix))//".xbin.tmp"
511 #endif
512         write(iout,*) mm(1)," conformations read",ll(1),
513      &    " conformations written to ", 
514      &    bprotfile_temp(:ilen(bprotfile_temp))
515       else
516         do islice=1,nslice
517           write (liczba1,'(bz,i2.2)') islice
518 #ifdef MPI
519           write (liczba,'(bz,i3.3)') me
520           bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
521      &      prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
522 #else
523           bprotfile_temp = scratchdir(:ilen(scratchdir))//
524      &       "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
525 #endif
526           write(iout,*) mm(islice)," conformations read",ll(islice),
527      &    " conformations written to ", 
528      &    bprotfile_temp(:ilen(bprotfile_temp))
529         enddo
530       endif
531
532 #ifdef MPI
533 c Check if everyone has the same number of conformations
534       call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
535      &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
536       lerr=.false.
537       do i=0,nprocs-1
538         if (i.ne.me) then
539           do islice=1,nslice
540           if (stot(islice).ne.ntot_all(islice,i)) then
541             write (iout,*) "Number of conformations at processor",i,
542      &       " differs from that at processor",me,
543      &       stot(islice),ntot_all(islice,i)," slice",islice
544             lerr = .true.
545           endif
546           enddo
547         endif
548       enddo 
549       if (lerr) then
550         write (iout,*)
551         write (iout,*) "Numbers of conformations read by processors"
552         write (iout,*)
553         do i=0,nprocs-1
554           write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
555         enddo
556         write (iout,*) "Calculation terminated."
557         call flush(iout)
558         return1
559       endif
560       do islice=1,nslice
561         ntot(islice)=stot(islice)
562       enddo
563       return
564 #endif
565  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
566       call flush(iout)
567       return1
568       end
569 c------------------------------------------------------------------------------
570       subroutine card_concat(card,to_upper)
571       implicit none
572       include 'DIMENSIONS.ZSCOPT'
573       include "COMMON.IOUNITS"
574       character*(*) card
575       character*80 karta,ucase
576       logical to_upper
577       integer ilen
578       external ilen
579       read (inp,'(a)') karta
580       if (to_upper) karta=ucase(karta)
581       card=' '
582       do while (karta(80:80).eq.'&')
583         card=card(:ilen(card)+1)//karta(:79)
584         read (inp,'(a)') karta
585         if (to_upper) karta=ucase(karta)
586       enddo
587       card=card(:ilen(card)+1)//karta
588       return
589       end
590 c------------------------------------------------------------------------------
591       subroutine readi(rekord,lancuch,wartosc,default)
592       implicit none
593       character*(*) rekord,lancuch
594       integer wartosc,default
595       integer ilen,iread
596       external ilen
597       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
598       if (iread.eq.0) then
599         wartosc=default
600         return
601       endif
602       iread=iread+ilen(lancuch)+1
603       read (rekord(iread:),*) wartosc
604       return
605       end
606 c----------------------------------------------------------------------------
607       subroutine reada(rekord,lancuch,wartosc,default)
608       implicit none
609       character*(*) rekord,lancuch
610       character*80 aux
611       double precision wartosc,default
612       integer ilen,iread
613       external ilen
614       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
615       if (iread.eq.0) then
616         wartosc=default
617         return
618       endif
619       iread=iread+ilen(lancuch)+1
620       read (rekord(iread:),*) wartosc
621       return
622       end
623 c----------------------------------------------------------------------------
624       subroutine multreadi(rekord,lancuch,tablica,dim,default)
625       implicit none
626       integer dim,i
627       integer tablica(dim),default
628       character*(*) rekord,lancuch
629       character*80 aux
630       integer ilen,iread
631       external ilen
632       do i=1,dim
633         tablica(i)=default
634       enddo
635       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
636       if (iread.eq.0) return
637       iread=iread+ilen(lancuch)+1
638       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
639    10 return
640       end
641 c----------------------------------------------------------------------------
642       subroutine multreada(rekord,lancuch,tablica,dim,default)
643       implicit none
644       integer dim,i
645       double precision tablica(dim),default
646       character*(*) rekord,lancuch
647       character*80 aux
648       integer ilen,iread
649       external ilen
650       do i=1,dim
651         tablica(i)=default
652       enddo
653       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
654       if (iread.eq.0) return
655       iread=iread+ilen(lancuch)+1
656       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
657    10 return
658       end
659 c----------------------------------------------------------------------------
660       subroutine reads(rekord,lancuch,wartosc,default)
661       implicit none
662       character*(*) rekord,lancuch,wartosc,default
663       character*80 aux
664       integer ilen,lenlan,lenrec,iread,ireade
665       external ilen
666       logical iblnk
667       external iblnk
668       lenlan=ilen(lancuch)
669       lenrec=ilen(rekord)
670       iread=index(rekord,lancuch(:lenlan)//"=")
671 c      print *,"rekord",rekord," lancuch",lancuch
672 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
673       if (iread.eq.0) then
674         wartosc=default
675         return
676       endif
677       iread=iread+lenlan+1
678 c      print *,"iread",iread
679 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
680       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
681         iread=iread+1
682 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
683       enddo
684 c      print *,"iread",iread
685       if (iread.gt.lenrec) then
686          wartosc=default
687         return
688       endif
689       ireade=iread+1
690 c      print *,"ireade",ireade
691       do while (ireade.lt.lenrec .and.
692      &   .not.iblnk(rekord(ireade:ireade)))
693         ireade=ireade+1
694       enddo
695       wartosc=rekord(iread:ireade)
696       return
697       end
698 c----------------------------------------------------------------------------
699       subroutine multreads(rekord,lancuch,tablica,dim,default)
700       implicit none
701       integer dim,i
702       character*(*) rekord,lancuch,tablica(dim),default
703       character*80 aux
704       integer ilen,lenlan,lenrec,iread,ireade
705       external ilen
706       logical iblnk
707       external iblnk
708       do i=1,dim
709         tablica(i)=default
710       enddo
711       lenlan=ilen(lancuch)
712       lenrec=ilen(rekord)
713       iread=index(rekord,lancuch(:lenlan)//"=")
714 c      print *,"rekord",rekord," lancuch",lancuch
715 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
716       if (iread.eq.0) return
717       iread=iread+lenlan+1
718       do i=1,dim
719 c      print *,"iread",iread
720 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
721       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
722         iread=iread+1
723 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
724       enddo
725 c      print *,"iread",iread
726       if (iread.gt.lenrec) return
727       ireade=iread+1
728 c      print *,"ireade",ireade
729       do while (ireade.lt.lenrec .and.
730      &   .not.iblnk(rekord(ireade:ireade)))
731         ireade=ireade+1
732       enddo
733       tablica(i)=rekord(iread:ireade)
734       iread=ireade+1
735       enddo
736       end
737 c----------------------------------------------------------------------------
738       subroutine split_string(rekord,tablica,dim,nsub)
739       implicit none
740       integer dim,nsub,i,ii,ll,kk
741       character*(*) tablica(dim)
742       character*(*) rekord
743       integer ilen
744       external ilen
745       do i=1,dim
746         tablica(i)=" "
747       enddo
748       ii=1
749       ll = ilen(rekord)
750       nsub=0
751       do i=1,dim
752 C Find the start of term name
753         kk = 0
754         do while (ii.le.ll .and. rekord(ii:ii).eq." ") 
755           ii = ii+1
756         enddo
757 C Parse the name into TABLICA(i) until blank found
758         do while (ii.le.ll .and. rekord(ii:ii).ne." ") 
759           kk = kk+1
760           tablica(i)(kk:kk)=rekord(ii:ii)
761           ii = ii+1
762         enddo
763         if (kk.gt.0) nsub=nsub+1
764         if (ii.gt.ll) return
765       enddo
766       return
767       end
768 c--------------------------------------------------------------------------------
769       integer function iroof(n,m)
770       ii = n/m
771       if (ii*m .lt. n) ii=ii+1
772       iroof = ii
773       return
774       end