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