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