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