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