text:
      program tapet
*
* print the wallpaper for a sos or pos marcs model
*
      character*128 file,modnam
10    print*, 'what model do you want to print?'
      read(5,'(a)') file
      open (1,file=file,status='old',form='unformatted',
     & err=10,convert='big_endian',recl=1848001)
      print*,' output file for listing?'
      read(5,'(a)') file
      open(20,file=file,status='unknown')
      call listmo(1,1)
      end

      SUBROUTINE LISTMO(MO,IARCH)
C        THIS ROUTINE PRINTS A NUMBER (MO) OF MODELS, WHICH ARE STORED ON
C        FORTRAN UNIT IARCH.
C
      include 'parameter.inc'
      include 'lambda.par'
C
      logical varimacro
      character*4 idrab1,idrab2,idrab3,idrab4,idrab5,idrab6
      character*128 modnam
      DIMENSION ABUND(16),TKORRM(NDP),FCORR(NDP),TAU(NDP),TAUS(NDP),
     *T(NDP),PE(NDP),PG(NDP),PRAD(NDP),PTURB(NDP),XKAPR(NDP),RO(NDP),
     *CP(NDP),CV(NDP),AGRAD(NDP),Q(NDP),U(NDP),V(NDP),ANCONV(NDP),
     *PRESMO(30,NDP),FCONV(NDP),RR(NDP),Z(NDP),EMU(NDP),HNIC(NDP),
     *NJ(16),XLR(20+1),IEL(16),ANJON(16,5),PART(16,5),PROV(mkomp,20+1),
     *ABSKA(mkomp),SPRIDA(mkomp),XLB(lpoint),FLUXME(lpoint),
     *  FLUMAG(lpoint),fluxmec(lpoint),
     &  PEP(16),ABNAME(mkomp),SOURCE(mkomp),DUMMY(11*NDP+2)
      DIMENSION W(lpoint),UW(12),BW(21),VW(25),RW(43),XIW(51),XJW(59),
     &  XKW(39)
      doubleprecision drr(ndp)
      CHARACTER*10 DAG,NAME,NAMEP,KLOCK
      CHARACTER*8 ABNAME,SOURCE
      DIMENSION WAVFLX(10),ptio(ndp)
      COMMON /UTPUT/IREAD,IWRIT
      COMMON /MASSE/RELM,RELLUM,RELATM,RELRAD,BOLMAG
      real    abSc,abTi,abV,abMn,abCo
      real inputabund(92)
      character*2 inputaname(92)
      real macrovel(ndp), macrobeta
      integer modtype
*
      common/specnew9/nlp,xlr,taus
* added for osplot *************
      real app(11*ndp),ptot(ndp)
      common /binstruc/bPPR(NDP),bPPT(NDP),bPP(NDP),bGG(NDP),
     & bZZ(NDP),bDD(NDP),
     & bVV(NDP),bFFC(NDP),bPPE(NDP),bTT(NDP),
     & bTAULN(NDP),NbTAU,IbTER,erad(ndp)
      common /struct/ teff,tau,t,z,ptot,prad,pg,pturb,pe,ro,rr,
     &                xkapr
      common /spectr/ nlb,xlb,w,fluxme
      common /spectrc/fluxmec
      common /pressure/ presmo,ptio
* *************
      DATA UW/0.145,0.436,0.910,1.385,1.843,2.126,2.305,2.241,1.270,
     *0.360,0.128,0.028/
      data BW/0.003,0.026,0.179,0.612,1.903,2.615,2.912,
     *3.005,2.990,2.876,2.681,2.388,2.058,1.725,1.416,1.135,0.840,0.568,
     *0.318,0.126,0.019/
      data VW/0.006,0.077,0.434,1.455,2.207,2.703,2.872,
     *2.738,2.505,2.219,1.890,1.567,1.233,0.918,0.680,0.474,0.312,0.200,
     *0.132,0.096,0.069,0.053,0.037,0.022,0.012/
      DATA  RW/0.03,0.06,0.17,0.28,0.39,0.50,0.60,0.69,0.74,0.79,
     &         0.84,0.88,0.91,0.94,0.96,0.98,0.99,1.00,0.97,0.94,
     &         0.90,0.85,0.79,0.73,0.65,0.57,0.50,0.42,0.36,0.31,
     &         0.24,0.17,0.14,0.11,0.08,0.06,0.05,0.04,0.03,0.02,
     &         0.02,0.01,0.01/
      DATA XIW/0.01,0.01,0.09,0.17,0.26,0.36,0.46,0.56,0.66,0.76,
     &         0.86,0.96,0.97,0.98,0.98,0.99,0.99,1.00,0.99,0.98,
     &         0.96,0.93,0.88,0.84,0.78,0.71,0.64,0.58,0.52,0.47,
     &         0.42,0.36,0.32,0.28,0.24,0.20,0.18,0.15,0.12,0.10,
     &         0.09,0.08,0.06,0.05,0.04,0.03,0.03,0.02,0.02,0.01,
     &         0.01/
      DATA XJW/0.01,0.02,0.02,0.03,0.04,0.06,0.11,0.16,0.26,0.35,
     &         0.49,0.62,0.78,0.93,0.89,0.85,0.82,0.78,0.78,0.78,
     &         0.79,0.80,0.82,0.85,0.89,0.93,0.84,0.75,0.70,0.64,
     &         0.64,0.63,0.63,0.63,0.64,0.66,0.67,0.68,0.69,0.70,
     &         0.70,0.70,0.68,0.66,0.63,0.60,0.53,0.46,0.37,0.27,
     &         0.20,0.14,0.12,0.09,0.08,0.06,0.04,0.02,0.01/
      DATA XKW/0.01,0.035,0.17,0.30,0.46,0.65,0.90,0.97,0.995,
     &         1.00,0.995,0.985,0.975,0.965,0.96,0.955,0.95,0.955,
     &         0.965,0.97,0.97,0.965,0.955,0.945,0.94,0.945,0.95,
     &         0.95,0.955,0.95,0.93,0.88,0.785,0.62,0.46,0.27,0.13,
     &         0.05,0.01/
* Ref HL Johnson ApJ 141,923  1965.
*
      DATA NAME/'LOCAL'/,NAMEP/'PARSONS'/
      DATA A,B/.34785485,.65214515/

*      save a,b,name,namep

      IREAD=5
      IWRIT=20
C
      REWIND IARCH
***** added for osplot ***************

        do i=ndp,1,-1
cc          print*,'oldsta; attempting to read binary model with i=',i
          rewind(iarch)
          read(iarch,err=499)
     &        bppr(1:i),bppt(1:i),bpp(1:i),bgg(1:i),bzz(1:i),bdd(1:i),
     &        bvv(1:i),bffc(1:i),bppe(1:i),btt(1:i),btauln(1:i),nbtau,
     &        ibter,erad(1:i)

          if (nbtau.gt.0.and.nbtau.le.i) then
* we should have found the right dimension
* But we have lost the compatibility with old marcs files that
* had no erad variable! BPz 05/07-2000
            print*,'readmo; binary model had ndp= ',i
            print*,'readmo;       we use now ndp= ',ndp
            print*,"readmo: prev model had ",nbtau," depthpoints"
            goto 497
          endif
499       continue
          if (i.eq.1) then
            print*,'ERROR! We have tried to read the input binary model'
            print*,'ERROR! with all possible values of ndp between 1 '
            print*,'ERROR! and ',ndp,'. And nothing works!'
            print*,'ERROR! Length of common block in stored model seems'
            print*,'ERROR! too long compared to NDP.'
            stop
          endif
        enddo
497     continue
*
c      do k=1,nbtau
c        print*,k,bppr(k),bppt(k),bpp(k),bgg(k),bzz(k)
c        print*,k,bdd(k),bvv(k),bffc(k),bppe(k),btt(k)
c        print*,k,btauln(k),erad(k)
c      enddo
      rewind(iarch)
***********************

      DO 1 IMO=1,MO
* first extract the last few records to get the full chemical composition
*  and the long name. BPz 11/01-2008
      do i=1,2
        read(iarch)
      enddo
      read(iarch)
     &     teff,flux,g,palfa,pny,py,pbeta,iline,istral,mihal,
     &            idrab1,idrab2,idrab3,idrab4,idrab5,idrab6,
     &            itmax,nel
      do i=1,2
        read(iarch) jtau
      enddo
      ntpo=0
      do k=1,jtau
        read(iarch) kr,tau(k)
        tauk=alog10(tau(k))+10.01
        ktau=tauk
        if(abs(tauk-ktau).gt.0.02) go to 131
        if(ktau.eq.10) k0=k
        ntpo=ntpo+1
  131   continue
      enddo

      read(iarch)(nj(i),i=1,nel),nlp
      do k=1,ntpo
        do j=1,nlp
          if (j.eq.1) then
            do i=1,nel
              read(iarch)
            enddo
          endif
          read(iarch)
        enddo
      enddo
*****************
c  there are
c 1) old models with nothing after flux record (very old)
c 2) models with longname after flux record (older BPz)
c 3) models with weight record after flux record, followed by long name
c    (Bpz on bohor)
c 4) the v7.3 models with flux record, fluxmec record, abundances record,
c    longname record (BPz v7.3 new models)
c 5) KE models with flux record, fluxmec record, abundance record
c
c flux record
      modtype=0
      read(iarch,err=192) nlb,(xlb(j),fluxme(j),j=1,nlb),(w(j),j=1,nlb)
      read(iarch,end=190,err=191) (fluxmec(j),j=1,nlb)
      read(iarch,err=193,end=193) inputabund,inputaname
      read(iarch,end=199) modnam
      modtype=4
      print*,'Model type 4'
      goto 199
190   modtype=1
      print*,'Model type 1'
      goto 199
191   backspace(iarch)
      read(iarch) modnam
      modtype=2
      print*,'Model type 2'
      goto 199
192   backspace(iarch)
      read(iarch) modnam
      modtype=3
      print*,'Model type 3'
      goto 199
193   modtype=5
      print*,'Model type 5'
      backspace(iarch)
      read(iarch) inputabund
      do j=1,92
        inputaname(j)=' '
      enddo
      goto 199
* we end up here with the right value for modtype
199   continue
      if (modtype.eq.0) then
        print*,'************'
        print*,' could not find model type!'
        print*,'************'
        stop
      endif
      rewind iarch
* Now do the full reading:
      READ(IARCH)
      print*,'******** just before inord record ***********'
      READ(IARCH) INORD,DAG,KLOCK
      print*,'******** just before abund record ***********'
      READ(IARCH,err=999)
     &     TEFF,FLUX,G,PALFA,PNY,PY,PBETA,ILINE,ISTRAL,MIHAL,
     &            IDRAB1,IDRAB2,IDRAB3,IDRAB4,IDRAB5,IDRAB6,
     &            ITMAX,NEL,(ABUND(I),I=1,NEL),abSc,abTi,abV,abMn,abCo,
     &            macrobeta,macrovel
      goto 998
999   continue
      backspace(iarch)
      print*,'******** just before abund record again (no macro)****'
      READ(IARCH)
     &     TEFF,FLUX,G,PALFA,PNY,PY,PBETA,ILINE,ISTRAL,MIHAL,
     &            IDRAB1,IDRAB2,IDRAB3,IDRAB4,IDRAB5,IDRAB6,
     &            ITMAX,NEL,(ABUND(I),I=1,NEL),abSc,abTi,abV,abMn,abCo
      macrobeta=0.
      do i=1,ndp
        macrovel(i)=0.0
      enddo
998   continue
      if (modtype.eq.1.or.modtype.eq.5) then
        write(modnam,195) idrab1,idrab2,idrab3,idrab4,idrab5,idrab6
195     format(6a4)
      endif
* maybe useful for some of KE's models
      print*,' trying double precision for rr(1:ndp)'
      READ(IARCH,err=1950,end=1950)
     &       JTAU,NCORE,DIFLOG,TAUM,RADIUS,(dRR(K),K=1,JTAU)
      do k=1,jtau
        rr(k)=sngl(drr(k))
      enddo
      goto 1951
1950  continue
      backspace(iarch)
      READ(IARCH)
     &       JTAU,NCORE,DIFLOG,TAUM,RADIUS,(RR(K),K=1,JTAU)
1951  continue
      print*,JTAU,NCORE,DIFLOG,TAUM,RADIUS,(RR(K),K=1,JTAU)
      backspace(iarch)
      write(20,219)
      GLOG=ALOG10(G)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,283)
      FNORD=0.1*INORD
c  fnord=3.0 = apollo version; =4.0 = Sun version
cccc      fnord=3.0
      relrad=log10(radius/6.9598E10)
      if (radius.le.2.) then
        write(20,281) FNORD,DAG,KLOCK
        relm=0.
        relatm=0.
        rellum=0.
        bolMAG=99.
      else
        write(20,282) FNORD,DAG,KLOCK
**************
        print*,'radius, rr(2) ',radius,rr(2)
        relm=10**(glog-4.44+2.*relrad)
CCCCC      RELATM=-3.+ALOG10(TSUN/TEFF/RELM)+0.5*RELLUM
        relatm=log10((rr(2)-radius)/radius)
        RELLUM=4.*ALOG10(TEFF/5780.)+2.*RELRAD
        BOLMAG=4.74-2.5*RELLUM
**************
      endif
      write(20,284)
      write(20,200)
C        CONVERT TO 'PHYSICAL FLUX'
      FLUX=3.14159*FLUX
      write(20,201) TEFF,FLUX,G
      if(radius.le.2.) then
        write(20,2001)
      else
        write(20,2002) 10.**relrad,radius,relm,10.**relatm,10.**rellum,
     ,                 bolmag
      endif
      write(20,2003)PALFA,PNY,PY
      IF(PBETA.LE.0. .and. macrobeta.le.0.) write(20,256)
      IF(PBETA.GT.0.0) write(20,257) PBETA
      IF(macrobeta.GT.0.0) then
        varimacro=.false.
        do i=1,jtau
          if (macrovel(i).ne.macrovel(1)) varimacro=.true.
        enddo
        if (varimacro) then
          write(20,259) macrobeta,macrovel(1)/1.e5
        else
          write(20,258) macrobeta,macrovel(1)/1.e5
        endif
      endif
      IF(ISTRAL.LT.1) write(20,250)
      IF(ISTRAL.GE.1) write(20,251)
      IF(ILINE.LT.1) write(20,252)
      IF(ILINE.GE.1) write(20,253)
      write(20,254) MIHAL,modnam(1:index(modnam,' '))
      DO I=1,NEL
        ABUND(I)=ALOG10(ABUND(I))+12.
      enddo
      if (modtype.eq.1.or.modtype.eq.2.or.modtype.eq.3) then
        write(20,202)
        write(20,203) (ABUND(I),I=1,NEL)
        write(20,2021)
        abSc=alog10(abSc)+12.
        abTi=alog10(abTi)+12.
        abV=alog10(abV)+12.
        abMn=alog10(abMn)+12.
        abCo=alog10(abCo)+12.
        write(20,203) absc,abti,abv,abmn,abco
      else
        j=1
        write(20,1202) (inputaname(i),i=(j-1)*20+1,min(92,(j-1)*20+20))
        write(20,1203) (inputabund(i),i=(j-1)*20+1,min(92,(j-1)*20+20))
        do j=2,5
          write(20,1204) (inputaname(i),i=(j-1)*20+1,
     &                            min(92,(j-1)*20+20))
          write(20,1203) (inputabund(i),i=(j-1)*20+1,
     &                            min(92,(j-1)*20+20))
        enddo
      endif
      write(20,255) ITMAX
      READ(IARCH)JTAU,NCORE,DIFLOG,TAUM,RADIUS,(RR(K),K=1,JTAU)
      READ(IARCH)JTAU,(TKORRM(I),I=1,JTAU),(FCORR(K),K=1,JTAU)
      NTPO=0
      DO 3 K=1,JTAU
      READ(IARCH) KR,TAU(K),TAUS(K),Z(K),T(K),PE(K),PG(K),PRAD(K),
     &            PTURB(K),XKAPR(K),RO(K),EMU(K),CP(K),CV(K),
     &            AGRAD(K),Q(K),U(K),V(K),ANCONV(K),HNIC(K),NMOL,
     &            (PRESMO(J,K),J=1,NMOL),PTIO(K)
      TAUK=ALOG10(TAU(K))+10.01
      KTAU=TAUK
      IF(ABS(TAUK-KTAU).GT.0.02) GO TO 31
      IF(KTAU.EQ.10) K0=K
      NTPO=NTPO+1
   31 CONTINUE
    3 CONTINUE
      write(20,204)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,205)
      DO 4 I=1,JTAU
      FCONV(I)=ANCONV(I)*FLUX
      write(20,206) I,TAU(I),T(I),TKORRM(I),FCONV(I),FCORR(I),I
    4 CONTINUE
***
      READ(IARCH)(NJ(I),I=1,NEL),NLP,(XLR(I),I=1,NLP)
     & ,NPROV,NPROVA,NPROVS,(ABNAME(KP),SOURCE(KP),KP=1,NPROV)
***
      write(20,207)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,208) XLR(NLP)
      Z0=Z(K0)
      DO 5 I=1,JTAU
      Z(I)=Z(I)-Z0
* added for osplot*************
        i1=min(i+1,jtau)
        PTOT(I)=PG(I)+PRAD(I)+.5*(pturb(i)+pturb(i1))
      write(20,209) I,TAU(I),TAUS(I),Z(I),T(I),PE(I),PG(I),PRAD(I),
     &             PTURB(I),XKAPR(I),I
    5 CONTINUE
      write(20,210)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,211)
      DO 6 I=1,JTAU
      write(20,212) I,TAU(I),RO(I),EMU(I),CP(I),CV(I),AGRAD(I),Q(I),
     &     U(I),V(I),ANCONV(I),I
    6 CONTINUE
      write(20,213)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,214)
      DO 7 I=1,JTAU
      HNIC(I)=ALOG10(HNIC(I))
      DO 77 J=1,NMOL
      PRESMO(J,I)=ALOG10(PRESMO(J,I))
   77 CONTINUE
       PTIO(I)=AMAX1(PTIO(I),1.E-30)
       PTIOLG=ALOG10(PTIO(I))
      write(20,215) I,HNIC(I),(PRESMO(JJ,I),JJ=1,13),PTIOLG,I
    7 CONTINUE
      write(20,213)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,216)
      DO I=1,JTAU
        write(20,217) I,(PRESMO(J,I),J=14,16),(PRESMO(J,I),J=18,30),I
      enddo
   10 CONTINUE
cc      write(20,260)
ccc      READ(IARCH)(NJ(I),I=1,NEL),NLP,(XLR(I),I=1,NLP)
ccc     & ,NPROV,NPROVA,NPROVS,(ABNAME(KP),SOURCE(KP),KP=1,NPROV)
      HCK=143922240.
      DO 22 KTAU=1,NTPO
      DO 20 IE=1,NEL
      NJP=NJ(IE)
      READ(IARCH) KR,TAUI,TI,PEI,IEL(IE),ABUND(IE),
     &            (ANJON(IE,JJ),JJ=1,NJP),(PART(IE,JJ),JJ=1,NJP)
   20 CONTINUE
*
************************************************************************
* try to improve computation of electron contribution:
* BPz 31/03-2004
*
* first, check that Pg=P(HI)+P(HII)+P(H2)+P(HeI)+Pe (approx, this is neglecting
* any metals, and supposing He is neutral)
* P(He)= fictpressure(H)*abund(He)
*
         Pcheck=10**hnic(kr)*(1.+anjon(1,2)/anjon(1,1))*(1.+abund(2))+
     &          Pei+
     &          10**presmo(2,kr)*(1+2.*abund(2))
         fictPH=10**hnic(kr)*(1.+anjon(1,2)/anjon(1,1))+
     &          10**presmo(2,kr)*2.
         print*,'(Pgas-pcheck)/pgas = ',(pg(kr)-pcheck)/pg(kr),kr
************************************************************************
*
      DO 21 KLAM=1,NLP
      READ(IARCH) KR,TAUIL,(PROV(J,KLAM),J=1,NPROV),
     &            ABSKA(KLAM),SPRIDA(KLAM)
   21 CONTINUE
      IF(KTAU.LE.1) write(20,218)
      IF(KTAU.GT.1) write(20,219)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,220) TAUI
      write(20,221) TI,PEI
      write(20,222)
      PESUM=0.
* treat hydrogen separately, to account for H2 formation
      pep(1)=10**hnic(kr)*anjon(1,2)/anjon(1,1)
      pesum=pesum+pep(1)
      DO IE=2,NEL
        PEP(IE)=0.
        NJP=NJ(IE)
        if (NJP.ge.2) then
          do JJ=2,NJP
            PEP(IE)=PEP(IE)+fictPH*ABUND(IE)*ANJON(IE,JJ)*(JJ-1)
          enddo
          PESUM=PESUM+PEP(IE)
        endif
      enddo
*
*BPz
      print*,' check  (Pe-Pesum)/pe = ', (pe(kr)-pesum)/pei
*
      DO 19 IE=1,NEL
      NJP=NJ(IE)
CCCC      PEP(IE)=PEP(IE)/PESUM
************************************************************************
* BPz:  renormalize to actual Pe. Missing electron donors will then
*       be accounted for.
*
      PEP(ie)=pep(ie)/pei
*
      if(njp.eq.2) then
        write(20,2232)IEL(IE),ABUND(IE),PEP(IE),(ANJON(IE,JJ),JJ=1,NJP),
     .                (part(ie,jj),jj=1,njp)
      elseif(njp.eq.3) then
        write(20,2233)IEL(IE),ABUND(IE),PEP(IE),(ANJON(IE,JJ),JJ=1,NJP),
     .                (part(ie,jj),jj=1,njp)
      elseif(njp.eq.4) then
        write(20,2234)IEL(IE),ABUND(IE),PEP(IE),(ANJON(IE,JJ),JJ=1,NJP),
     .                (part(ie,jj),jj=1,njp)
      endif
ccccc      write(20,224) (PART(IE,JJ),JJ=1,NJP)
   19 CONTINUE
      write(20,2250)
*
      write(20,225) (ABNAME(KP),KP=1,20)
      HCKT=HCK/TI
      DO 18 KLAM=1,NLP
      ABKLA=ABSKA(KLAM)
      STIM=1.-EXP(-HCKT/XLR(KLAM))
      DO J=1,NPROVA
        PROV(J,KLAM)=PROV(J,KLAM)/ABKLA*STIM
      enddo
      DO J=1,NPROVS
        PROV(NPROVA+J,KLAM)=PROV(NPROVA+J,KLAM)/SPRIDA(KLAM)
      enddo
      write(20,226) XLR(KLAM),ABSKA(KLAM),SPRIDA(KLAM),
     &             (PROV(J,KLAM),J=1,20)
   18 CONTINUE
c
      write(20,2251) (ABNAME(KP),KP=21,NPROV)
      DO 188 KLAM=1,NLP
        write(20,2261) XLR(KLAM),(PROV(J,KLAM),J=21,NPROV)
  188 CONTINUE

   22 CONTINUE
*****************
c  there are
c 1) old models with nothing after flux record (very old)
c 2) models with longname after flux record (older BPz)
c 3) models with weight record after flux record, followed by long name
c    (Bpz on bohor)
c 4) the v7.3 models with flux record, fluxmec record, abundances record,
c    longname record (BPz v7.3 new models)
c 5) KE models with flux record, fluxmec record, abundance record
c
      if (modtype.eq.3) then
        READ(IARCH) NLB,(XLB(J),FLUXME(J),J=1,NLB)
        read(iarch) (W(J),J=1,NLB)
      else
        READ(IARCH) NLB,(XLB(J),FLUXME(J),J=1,NLB),(W(J),J=1,NLB)
      endif
      if (modtype.eq.4.or.modtype.eq.5) then
        read(iarch,end=79,err=79) (fluxmec(j),j=1,nlb)
        print*,(fluxmec(j),j=10000,10010)
      endif
      goto 78
79    print*,'********************************************************'
      print*,'Error while reading fluxmec in model; skipping that part'
      print*,'********************************************************'
78    continue
C        CONVERT TO 'PHYSICAL' FLUXES
      DO J=1,NLB
        FLUXME(J)=3.14159*FLUXME(J)
        fluxmec(j)=3.14159*fluxmec(j)
      enddo
      DO J=1,NLB
        JNORM=J
        IF(XLB(J).GT.5000.) GOTO 26
      enddo
   26 STMAGN=-2.5*ALOG10(FLUXME(JNORM))
      DO I=1,NLB
        FLUMAG(I)=-2.5*ALOG10(FLUXME(I))-STMAGN
      enddo
C
C        write(20,FLUXES)
C
**** here we jump over fluxes print out
      goto 1111
*
      write(20,521)
      N=0
      DO 52 I=1,NLB,4
      WAVFLX(N+1)=(XLB(I)+XLB(I+1)+XLB(I+2)+XLB(I+3))/4.
      WAVFLX(N+2)=100.*(A*(FLUXME(I)+FLUXME(I+3))+B*(FLUXME(I+1)+
     &      FLUXME(I+2)))
      N=N+2
      IF(N.LT.10) GO TO 52
      write(20,520) WAVFLX
      N=0
   52 CONTINUE
      IF(N.NE.0) write(20,520)(WAVFLX(I),I=1,N)
  520 FORMAT(1X,5(F12.0,E13.5))
  521 FORMAT('1  BELL''S FLUX              (CENTER WAVE & FLUX)'//)
C
C
      write(20,219)
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,227)
      LLB=NLB
   60 IF(LLB.LT.NLB) write(20,219)
      IF(LLB.LT.NLB) write(20,300) TEFF,GLOG,
     &       modnam(1:(index(modnam,' ')))
      IF(LLB.LT.NLB) write(20,227)
      write(20,228)
      IRAD = 0
      IF(LLB.LT.NLB) GOTO 65
      I=0
      I2=50
      I3=100
      I4=150
      GOTO 650
   65 CONTINUE
      I = I + 149
      I2 = I2 + 149
      I3 = I3 + 149
      I4 = I4 + 149
  650 CONTINUE
      ILI=2
      IF(LLB.GE.51) ILI=3
      IF(LLB.GE.101)ILI=4
      IF(LLB.GE.151)ILI=5
   66 I4=I4+1
   67 I3=I3+1
   68 I2=I2+1
   69 I=I+1
      IRAD = IRAD + 1
      IF(IRAD.GT.50) GOTO 70
      GOTO(70,61,62,63,64),ILI
   64 write(20,229) XLB(I),FLUXME(I),FLUMAG(I),XLB(I2),FLUXME(I2),
     &             FLUMAG(I2),XLB(I3),FLUXME(I3),FLUMAG(I3),
     &             XLB(I4),FLUXME(I4),FLUMAG(I4)
      IF(I4.GE.NLB) ILI=4
      GOTO 66
   63 write(20,229) XLB(I),FLUXME(I),FLUMAG(I),XLB(I2),FLUXME(I2),
     &             FLUMAG(I2),XLB(I3),FLUXME(I3),FLUMAG(I3)
      IF(I3.GE.NLB) ILI=3
      GOTO 67
   62 write(20,229) XLB(I),FLUXME(I),FLUMAG(I),XLB(I2),FLUXME(I2),
     &             FLUMAG(I2)
      IF(I2.GE.NLB) ILI=2
      GOTO 68
   61 write(20,229) XLB(I),FLUXME(I),FLUMAG(I)
      IF(I.GE.NLB) GOTO 70
      GOTO 69
   70 LLB = LLB - 200
      IF(LLB.GT.0) GOTO 60
*
C
C        COMPUTE U, B, V, R, I AND COLOURS
      SUMCOL=0.
      SUMNOR=0.
      write(20,270)
      DO 40 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.3000..OR.XXLB.GE.4200.) GOTO 40
      IINDEX=INT(XXLB/100.)-29
      VIKT=W(I)*UW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   40 CONTINUE
      UFLUX=SUMCOL/SUMNOR
      SUMCOL=0.
      SUMNOR=0.
      DO 41 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.3500..OR.XXLB.GE.5600.) GOTO 41
      IINDEX=INT(XXLB/100.)-34
      VIKT=W(I)*BW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   41 CONTINUE
      BFLUX=SUMCOL/SUMNOR
      SUMCOL=0.
      SUMNOR=0.
      DO 42 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.4700..OR.XXLB.GE.7200.) GOTO 42
      IINDEX=INT(XXLB/100.)-46
      VIKT=W(I)*VW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   42 CONTINUE
      VFLUX=SUMCOL/SUMNOR
* new colours BPz 290890
      SUMCOL=0.
      SUMNOR=0.
      DO 43 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.5300..OR.XXLB.GE.9500.) GOTO 43
      IINDEX=INT(XXLB/100.)-52
      VIKT=W(I)*RW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   43 CONTINUE
      RFLUX=SUMCOL/SUMNOR
      SUMCOL=0.
      SUMNOR=0.
      DO 44 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.6900..OR.XXLB.GE.12000.) GOTO 44
      IINDEX=INT(XXLB/100.)-68
      VIKT=W(I)*XIW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   44 CONTINUE
      XIFLUX=SUMCOL/SUMNOR
      SUMCOL=0.
      SUMNOR=0.
      DO 45 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.9700..OR.XXLB.GE.15600.) GOTO 45
      IINDEX=INT(XXLB/100.)-96
      VIKT=W(I)*XJW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   45 CONTINUE
      XJFLUX=SUMCOL/SUMNOR
      SUMCOL=0.
      SUMNOR=0.
      DO 46 I=1,NLB
      XXLB=XLB(I)
      IF(XXLB.LT.18100..OR.XXLB.GT.25900.) GOTO 46
      IINDEX=INT((XXLB/200.)-89.5)
      VIKT=W(I)*XKW(IINDEX)
      SUMNOR=SUMNOR+VIKT
      SUMCOL=SUMCOL+VIKT*FLUXME(I)
   46 CONTINUE
      XKFLUX=SUMCOL/SUMNOR
      UMAG=-2.5*ALOG10(UFLUX)-STMAGN
      BMAG=-2.5*ALOG10(BFLUX)-STMAGN
      VMAG=-2.5*ALOG10(VFLUX)-STMAGN
      RMAG=-2.5*ALOG10(RFLUX)-STMAGN
      XIMAG=-2.5*ALOG10(XIFLUX)-STMAGN
      XJMAG=-2.5*ALOG10(XJFLUX)-STMAGN
      XKMAG=-2.5*ALOG10(XKFLUX)-STMAGN
      UB=UMAG-BMAG
      BV=BMAG-VMAG
      UV=UMAG-VMAG
      RI=RMAG-XIMAG
      VR=VMAG-RMAG
      VI=VMAG-XIMAG
      VJ=VMAG-XJMAG
      VK=VMAG-XKMAG
      write(20,300) TEFF,GLOG,modnam(1:(index(modnam,' ')))
      write(20,2698)
      write(20,271) UB,BV,UV
      write(20,272) XLB(JNORM),RI,VR,VI,VJ,VK
      write(20,273) UMAG,BMAG,VMAG,RMAG,XIMAG,xjmag
      write(20,274) XKMAG
***** here we come if we don't print the fluxes nor the colours
1111  continue
    1 CONTINUE
  200 FORMAT('M o d e l   P a r a m e t e r s   e t c.'/)
  281 FORMAT(' * The following model was computed by  MARCS35 PP  '
     & ,F5.1,2X,A10,2X,A8,'  *')
  282 FORMAT(' * The following model was computed by  MARCS35 SPH '
     & ,F5.1,2X,A10,2X,A8,'  *')
  283 FORMAT(1X,81('*'))
  284 FORMAT(1X,81('*')////)
  201 FORMAT('Effective temperature',6x,F12.0,'  K'/' Total ',
     +'flux (= sigma*Teff**4)',1PE11.3,'  erg/s/cm**2'/' Acceleration',
     +' of gravity',0PF16.3,'  cm/s**2')
 2001 format(' Plane-parallel model')
 2002 format(' Spherical model with  Radius',1p,e11.3,'  Solar radii',
     ,  8x,'or ',3x,e11.3,' cm',/23x,
     , 'Mass  ',0p,f11.2,'  Solar masses'/23x,'Atm/Radius',2p,f7.2,'  %',
     , /23x,'Lumin.',1p,e11.2,'  Solar luminosities',' or M(bol)',
     , 0p,f8.3)
 2003 format(/' Convection Parameters'/' PALF',
     *'A (L/Hp)=',F5.2,',  PNY (ny)=',F5.2,',  PY (y)=',F5.3/)
  256 FORMAT(' Turbulence pressure is neglected'/)
  257 FORMAT(' Turbulence pressure from convection is included and ',
     ,'PBETA=',F4.2/)
  258 FORMAT(' Macroturbulence pressure ',
     & 'is included and BETA=',F4.2,' Vmacro= ',f5.2,
     &   ' km/s   constant with depth'/)
  259 FORMAT(' Macroturbulence pressure ',
     & 'is included and BETA=',F4.2,' Vmacro(1)= ',f5.2,
     &   ' km/s   variable with depth'/)
  250 FORMAT(' Convection has been included in this model')
  251 FORMAT(' Convection has  NOT  been included in this model')
  253 FORMAT(' Line blanketing has been included in this model')
  252 FORMAT(' Line blanketing has  NOT  been included in this model')
  254 FORMAT(' The Stroemgren equation has been used for the uppermost',
     *I4,' points'/' Identification    ',A)
  255 FORMAT(//' The following model was obtained after',I3,' itera',
     *'tion(s)')
c...v....1....v....2....v....3....v....4....v....5....v....6....v....7....v....8
  202 FORMAT(//'  Log. Abundances used in model calculations'/2X,'H'
     *,5X,'He',4X,'C',5X,'N',5X,'O',5X,'Ne',4X,'Na',4X,'Mg',4X,'Al',4X,
     *'Si',4X,'S',5X,'K',5X,'Ca',4X,'Cr',4X,'Fe',4X,'Ni')
 2021 format(/2x,'Sc',4x,'Ti',4x,'V',5x,'Mn',4x,'Co')
  203 FORMAT (16F6.2)
 1202 format(//'  Log. Abundances used in model calculations'/
     &         20(2x,a2,2x))
 1204 format(20(2x,a2,2x))
 1203 FORMAT (20F6.2)
c...v....1....v....2....v....3....v....4....v....5....v....6....v....7....v....8
  204 FORMAT(////,'C o r r e c t i o n s   i n   t h e   l a s t   i t e
     * r a t i o n')
  205 FORMAT(/3X,'K',4X,'TauRoss',6X,'T',7X,'Delta(T)',5X,'Fconv',
     *4X,'Delta(Fconv)',2X,'K')
  206 FORMAT(I4,F13.7,F9.2,F11.2,1PE12.3,E12.2,I5)
  207 FORMAT(////,'M o d e l   A t m o s p h e r e     (cgs units)')
  208 FORMAT(/4X,'K',5X,'TauRoss',5X,'Tau(',f6.0,')',1X,
     & 'Geom. Depth',
     & 5X,
     &'T',8X,'Pe',10X,'Pg',9X,'Prad',8X,'Pturb',4X,'KappaRoss',5X,'K')
  209 FORMAT(I5,F12.7,F13.7,1PE13.3,0PF9.0,5(1PE12.3),I6)
c...v....1....v....2....v....3....v....4....v....5....v....6....v....7....v....8
  210 FORMAT(////'T h e r m o d y n a m i c a l   q u a n t i t i e s
     *and   C o n v e c t i o n   (cgs units)')
  211 FORMAT(/3X,'K',4X,'TauRoss',4X,'Density',5X,'Mu',9X,'Cp',10X,
     *'Cv',7X,'AdGrad',9X,'Q',7X,'Sound vel.',2X,'Conv. vel.',1X,'Fconv/
     *F',3X,'K')
  212 FORMAT(I4,F12.5,1PE11.3,0PF8.3,6(1PE12.3),0PF9.5,I4)
  213 FORMAT(////'L o g a r i t h m i c   M o l e c u l a r   p a r t i
     *a l   p r e s s u r e s    (cgs units)')
  214 FORMAT(/'  K',3X,'P(H)',2X,'P(H-)',2X,'P(H2)',2X,'P(H2+)',1X,
     *'P(H2O)',1X,'P(OH)',2X,'P(CH)',2X,'P(CO)',2X,'P(CN)',2X,'P(C2)',2X
     *,'P(N2)',2X,'P(O2)',2X,'P(NO)',2X,'P(NH)',1X,'P(TiO)',4X,'K')
  215 FORMAT(I3,15F7.2,2X,I3)
  216 FORMAT(//'  K',3X,'C2H2',3X,'HCN',4X,'C2H',4X,'HS',5X,'SiH',4X
     * ,'C3H',4X,'C3',5X,'CS',5X,'SiC',4X,'SiC2',3X,'NS',5X,'SiN',4X
     * ,'SiO',4X,'SO',5X,'S2',5X,'SiS',5X,'K')
  217 FORMAT(I3,16F7.2,2X,I3)
  218 FORMAT(////'I o n i z a t i o n   c o n d i t i o n s   and   A b
     *s o r p t i o n  m e c h a n i s m s')
  219 FORMAT(////)
  220 FORMAT(/' *******************'/' * Tau=',F11.7,' *'/' ************
     ********')
  221 FORMAT(/' T=',F7.0,'  Pe=',1PE9.2)
  222 FORMAT(/'  Element  Abundance'
     *                      ,' ElCntb Ionization fractions',17X,'Partiti
     *on functions'/30X,'I',7X,'II',6X,'III',5X,'IV',12X,'I',9X,'II',8X,
     *'III',7X,'IV')
 2232 FORMAT(6X,A2,1PE12.3,0PF6.3,2(F8.4),16x,4x,1p,2e10.2)
 2233 FORMAT(6X,A2,1PE12.3,0PF6.3,3(F8.4),8x,4x,1p,3e10.2)
 2234 FORMAT(6X,A2,1PE12.3,0PF6.3,4(F8.4),4x,1p,4e10.2)
  223 FORMAT(6X,A2,1PE12.3,0PF6.3,4(F8.4))
  224 FORMAT(62X,4(1PE10.2))
c...v....1....v....2....v....3....v....4....v....5....v....6....v....7....v....8
 2250 FORMAT(//16X,'P e r c e n t a g e s   o f   C o n t i n u o u s'
     * ,'   A b s o r p t i o n   a n d   S c a t t e r i n g')
  225 FORMAT(/'  Wavel     Abs     Scat   ',20(A5))
 2251 format(/'  Wavel     ',24(a5))
  226 FORMAT(F7.0,2(1pE9.2),2p,20F5.1)
 2261 format(f7.0,3x,2p,24f5.1)
  227 FORMAT(' F l u x e s    (Physical fluxes in ergs/s/cm**2/Angstrom)
     *'/)
  228 FORMAT(4(6X,'LAMBDA',5X,'FLUX',7X,'MAGN'))
  229 FORMAT(4(F12.0,1PE11.2,0PF9.3))
  260 FORMAT(///' Molecules other than H2 and H2+ have not been consider
     *ed in the  Temperature - Electron pressure - Pressure  balance')
 2698 FORMAT(//' (UBV Transmission functions after Matthews + Sandage.
     *  Air mass = 1)')
  270 FORMAT(////'G I A N T  L I N E  C O L O U R S')
  271 FORMAT(/8X,'U - B =',F8.3//9X,'B - V =',F8.3//9X,'U - V =',F8.3
     *////)
  272 FORMAT(' TENTATIVE CONTINUUM COLOURS  (R AT',F6.0,'A AND I AT 9000
     *.A)'//9X,'R - I =',F8.3,//9X,'V - R =',F8.3//9X,'V - I =',F8.3
     &//9X,'V - J =',F8.3//9X,'V - K =',F8.3)
  273 FORMAT(////' U =',F10.3,'  B =',F10.3,'  V =',F10.3,'  R =',F10.3,
     *'  I =',F10.3,'  J =',F10.3)
  274 FORMAT(//' K =',F10.3,'  L =',F10.3,'  M =',F10.3,'  N =',F10.3)
  300 FORMAT('Teff=',F6.0,' log g=',F5.2,1X,a)
        print*,nlb
        print*,(xlb(j),j=10000,10010)
        print*,(w(j),j=10000,10010)
        print*,(fluxme(j),j=10000,10010)
        print*,(fluxmec(j),j=10000,10010)
        print*,fluxme(107855), lpoint,fluxme(lpoint)
        print*,'*** returning from readmo'
      RETURN
         END

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
      SUBROUTINE NEWSTR (START, NEWID)

C  NEWSTR searches for a non-existent structure identifier, starting
C  with START, and returns it, if found, in NEWID.

C structure status indicator
      INTEGER    PSNOEX,     PSEMPT,     PSNEMP
      PARAMETER (PSNOEX = 0, PSEMPT = 1, PSNEMP = 2)

      INTEGER     START, NEWID, IX, STRSTI, ERRIND

      DO 100 IX = START, START+100
         CALL PQSTST (IX, ERRIND, STRSTI)
         CALL CHKINQ ('pqstst', ERRIND)
         IF (STRSTI .EQ. PSNOEX) THEN
            NEWID = IX
            RETURN
         ENDIF
100   CONTINUE
      CALL UNMSG ('NEWSTR cannot find a new structure identifier.')

      END

ext: f


FUNCTION gser_s(a,x,gln)
  USE mo_kind, only: sp, i4
  USE mo_nrutil, ONLY : nrerror
  use mo_nr, ONLY : gammln
  IMPLICIT NONE
  REAL(SP), INTENT(IN) :: a,x
  REAL(SP), OPTIONAL, INTENT(OUT) :: gln
  REAL(SP) :: gser_s
  INTEGER(I4), PARAMETER :: ITMAX=100
  REAL(SP), PARAMETER :: EPS=epsilon(x)
  INTEGER(I4) :: n
  REAL(SP) :: ap,del,summ
  if (x == 0.0) then
     gser_s=0.0
     RETURN
  end if
  ap=a
  summ=1.0_sp/a
  del=summ
  do n=1,ITMAX
     ap=ap+1.0_sp
     del=del*x/ap
     summ=summ+del
     if (abs(del) < abs(summ)*EPS) exit
  end do
  if (n > ITMAX) call nrerror('a too large, ITMAX too small in gser_s')
  if (present(gln)) then
     gln=gammln(a)
     gser_s=summ*exp(-x+a*log(x)-gln)
  else
     gser_s=summ*exp(-x+a*log(x)-gammln(a))
  end if
END FUNCTION gser_s

FUNCTION gser_v(a,x,gln)
  USE mo_kind, only: sp, i4, lgt
  USE mo_nrutil, ONLY : assert_eq,nrerror
  use mo_nr, ONLY : gammln
  IMPLICIT NONE
  REAL(SP), DIMENSION(:), INTENT(IN) :: a,x
  REAL(SP), DIMENSION(:), OPTIONAL, INTENT(OUT) :: gln
  REAL(SP), DIMENSION(size(a)) :: gser_v
  INTEGER(I4), PARAMETER :: ITMAX=100
  REAL(SP), PARAMETER :: EPS=epsilon(x)
  INTEGER(I4) :: n
  REAL(SP), DIMENSION(size(a)) :: ap,del,summ
  LOGICAL(LGT), DIMENSION(size(a)) :: converged,zero
  n=assert_eq(size(a),size(x),'gser_v')
  zero=(x == 0.0)
  where (zero) gser_v=0.0
  ap=a
  summ=1.0_sp/a
  del=summ
  converged=zero
  do n=1,ITMAX
     where (.not. converged)
        ap=ap+1.0_sp
        del=del*x/ap
        summ=summ+del
        converged = (abs(del) < abs(summ)*EPS)
     end where
     if (all(converged)) exit
  end do
  if (n > ITMAX) call nrerror('a too large, ITMAX too small in gser_v')
  if (present(gln)) then
     if (size(gln) < size(a)) call &
          nrerror('gser_v: Not enough space for gln')
     gln=gammln(a)
     where (.not. zero) gser_v=summ*exp(-x+a*log(x)-gln)
  else
     where (.not. zero) gser_v=summ*exp(-x+a*log(x)-gammln(a))
  end if
END FUNCTION gser_v

FUNCTION dgser_s(a,x,gln)
  USE mo_kind, only: dp, i4
  USE mo_nrutil, ONLY : nrerror
  use mo_nr, ONLY : gammln
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: a,x
  REAL(DP), OPTIONAL, INTENT(OUT) :: gln
  REAL(DP) :: dgser_s
  INTEGER(I4), PARAMETER :: ITMAX=100
  REAL(DP), PARAMETER :: EPS=epsilon(x)
  INTEGER(I4) :: n
  REAL(DP) :: ap,del,summ
  if (x == 0.0_dp) then
     dgser_s=0.0_dp
     RETURN
  end if
  ap=a
  summ=1.0_dp/a
  del=summ
  do n=1,ITMAX
     ap=ap+1.0_dp
     del=del*x/ap
     summ=summ+del
     if (abs(del) < abs(summ)*EPS) exit
  end do
  if (n > ITMAX) call nrerror('a too large, ITMAX too small in dgser_s')
  if (present(gln)) then
     gln=gammln(a)
     dgser_s=summ*exp(-x+a*log(x)-gln)
  else
     dgser_s=summ*exp(-x+a*log(x)-gammln(a))
  end if
END FUNCTION dgser_s

FUNCTION dgser_v(a,x,gln)
  USE mo_kind, only: dp, i4, lgt
  USE mo_nrutil, ONLY : assert_eq,nrerror
  use mo_nr, ONLY : gammln
  IMPLICIT NONE
  REAL(DP), DIMENSION(:), INTENT(IN) :: a,x
  REAL(DP), DIMENSION(:), OPTIONAL, INTENT(OUT) :: gln
  REAL(DP), DIMENSION(size(a)) :: dgser_v
  INTEGER(I4), PARAMETER :: ITMAX=100
  REAL(DP), PARAMETER :: EPS=epsilon(x)
  INTEGER(I4) :: n
  REAL(DP), DIMENSION(size(a)) :: ap,del,summ
  LOGICAL(LGT), DIMENSION(size(a)) :: converged,zero
  n=assert_eq(size(a),size(x),'dgser_v')
  zero=(x == 0.0_dp)
  where (zero) dgser_v=0.0_dp
  ap=a
  summ=1.0_dp/a
  del=summ
  converged=zero
  do n=1,ITMAX
     where (.not. converged)
        ap=ap+1.0_dp
        del=del*x/ap
        summ=summ+del
        converged = (abs(del) < abs(summ)*EPS)
     end where
     if (all(converged)) exit
  end do
  if (n > ITMAX) call nrerror('a too large, ITMAX too small in dgser_v')
  if (present(gln)) then
     if (size(gln) < size(a)) call &
          nrerror('dgser_v: Not enough space for gln')
     gln=gammln(a)
     where (.not. zero) dgser_v=summ*exp(-x+a*log(x)-gln)
  else
     where (.not. zero) dgser_v=summ*exp(-x+a*log(x)-gammln(a))
  end if
END FUNCTION dgser_v


text:
Jump to Timeline #Navigation Navigation

  Flood hits
  Tule House destroyed
  State Senator  O. B. Powers
  Sheriff and tax collector  C. H. Gray

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
-makelib ies_lib/xil_defaultlib -sv \
  "C:/Xilinx/Vivado/2017.4/data/ip/xpm/xpm_fifo/hdl/xpm_fifo.sv" \
  "C:/Xilinx/Vivado/2017.4/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
  "C:/Xilinx/Vivado/2017.4/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
-endlib
-makelib ies_lib/xpm \
  "C:/Xilinx/Vivado/2017.4/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib ies_lib/axi_infrastructure_v1_1_0 \
  "../../../../AXI2TCM_sim_VHDL.srcs/sources_1/bd/design_1/ipshared/ec67/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib ies_lib/smartconnect_v1_0 -sv \
  "../../../../AXI2TCM_sim_VHDL.srcs/sources_1/bd/design_1/ipshared/02c8/hdl/sc_util_v1_0_vl_rfs.sv" \
-endlib
-makelib ies_lib/axi_protocol_checker_v2_0_1 -sv \
  "../../../../AXI2TCM_sim_VHDL.srcs/sources_1/bd/design_1/ipshared/3b24/hdl/axi_protocol_checker_v2_0_vl_rfs.sv" \
-endlib
-makelib ies_lib/xil_defaultlib -sv \
  "../../../bd/design_1/ip/design_1_axi_vip_0_0/sim/design_1_axi_vip_0_0_pkg.sv" \
-endlib
-makelib ies_lib/axi_vip_v1_1_1 -sv \
  "../../../../AXI2TCM_sim_VHDL.srcs/sources_1/bd/design_1/ipshared/a16a/hdl/axi_vip_v1_1_vl_rfs.sv" \
-endlib
-makelib ies_lib/xil_defaultlib -sv \
  "../../../bd/design_1/ip/design_1_axi_vip_0_0/sim/design_1_axi_vip_0_0.sv" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/design_1/sim/design_1.vhd" \
  "../../../bd/design_1/ip/design_1_AXI2TCM_v1_0_0_0/sim/design_1_AXI2TCM_v1_0_0_0.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  glbl.v
-endlib


text:
***
      SUBROUTINE kick(kw,m1,m1n,m2,ecc,sep,jorb,v,ffb,ECS,mfin)
      implicit none
*
      integer kw,k
      INTEGER idum
      COMMON /VALUE3/ idum
      INTEGER idum2,iy,ir(32)
      COMMON /RAND3/ idum2,iy,ir
      real*8 m1,m2,m1n,ecc,sep,jorb,ecc2,mrem,mean_ej
      real*8 pi,twopi,gmrkm,yearsc,rsunkm
      parameter(yearsc=3.1557d+07,rsunkm=6.96d+05)
      real*8 mm,em,dif,der,del,r
      real*8 u1,u2,vk,v(4),s,theta,phi,ss
      real*8 sphi,cphi,stheta,ctheta,salpha,calpha
      real*8 vr,vr2,vk2,vn2,hn2
      real*8 mu,cmu,vs(3),v1,v2,mx1,mx2
      real*8 maxwellian
      real*8 sigma1,sigma2,mxns
      integer bhflag
      COMMON /VALUE4/ sigma1,sigma2,mxns,bhflag
      real ran3,xx
      external ran3
*
      INTEGER ceflag,tflag,ifflag,nsflag,wdflag,piflag
      COMMON /FLAGS/ ceflag,tflag,ifflag,nsflag,wdflag,piflag
*
      INTEGER ECS
      real*8 ffb,mfin
*
      do k = 1,3
         vs(k) = 0.d0
      enddo
*     if(kw.eq.14.and.bhflag.eq.0) goto 95
*
* To take into account the ECS
      if(ECS .eq. 1)then
* Check if the supernova is driven by electron-cupture...
          maxwellian = sigma2
      else
* ... or iron core collapse.
          maxwellian = sigma1
      endif
*
      pi = ACOS(-1.d0)
      twopi = 2.d0*pi
* Conversion factor to ensure velocities are in km/s using mass and
* radius in solar units.
      gmrkm = 1.906125d+05
*
* Find the initial separation by randomly choosing a mean anomaly.
      if(sep.gt.0.d0.and.ecc.ge.0.d0)then
         xx = RAN3(idum)
         mm = xx*twopi
         em = mm
 2       dif = em - ecc*SIN(em) - mm
         if(ABS(dif/mm).le.1.0d-04) goto 3
         der = 1.d0 - ecc*COS(em)
         del = dif/der
         em = em - del
         goto 2
 3       continue
         r = sep*(1.d0 - ecc*COS(em))
*
* Find the initial relative velocity vector.
         salpha = SQRT((sep*sep*(1.d0-ecc*ecc))/(r*(2.d0*sep-r)))
         calpha = (-1.d0*ecc*SIN(em))/SQRT(1.d0-ecc*ecc*COS(em)*COS(em))
         vr2 = gmrkm*(m1+m2)*(2.d0/r - 1.d0/sep)
         vr = SQRT(vr2)
      else
         vr = 0.d0
         vr2 = 0.d0
         salpha = 0.d0
         calpha = 0.d0
      endif
*
* Generate Kick Velocity using Maxwellian Distribution (Phinney 1992).
* Use Henon''s method for pairwise components (Douglas Heggie 22/5/97).
      do 20 k = 1,2
         u1 = RAN3(idum)
         u2 = RAN3(idum)
* Generate two velocities from polar coordinates S & THETA.
* way to avoid negative values in LOG
         ss = 1.d0 - u1
         if(ss.le.0.d0)then
            ss = 1.d-6
         endif
         s = maxwellian*SQRT(-2.d0*LOG(ss))
         theta = twopi*u2
         v(2*k-1) = s*COS(theta)
         v(2*k) = s*SIN(theta)
 20   continue
      vk2 = v(1)**2 + v(2)**2 + v(3)**2
      vk = SQRT(vk2)
      v(4) = vk
*
***
* Generate Kick Velocity for Black Holes...
* Standard case
      if(kw.eq.14.and.bhflag.eq.1)then
         mrem = m1n
         vk = vk*(1.d0 - ffb)
         vk2 = vk*vk
*     ... if they are generated via direct collapse vk = 0
***
* Generate Kick Velocity for Black Holes and Neutron Stars
* Fallback case
      elseif((kw.eq.14.or.kw.eq.13).and.(bhflag.eq.2))then
         mrem = m1n
         vk = vk*(1.d0 - ffb)
         vk2 = vk*vk
***
* Generate Kick Velocity for Black Holes and Neutron Stars
* Mass Ejected case:
*   vk = vk*(mej/<mej>)*(1.2/mrem)
      elseif((kw.eq.14.or.kw.eq.13).and.(bhflag.eq.3))then
         mrem = m1n
*
         if(kw.eq.14 .and. ffb.eq.1.d0)then
            mfin = mrem
         endif
* mean mass ejected <mej> in Msun (in a population of 10^6 single stars)
* model 1
         if(nsflag.eq.2) mean_ej = 9.d0
         if(nsflag.eq.3) mean_ej = 7.5d0
         vk = vk*((mfin-mrem)/mean_ej)*(1.2d0/mrem)
         vk2 = vk*vk
***
* Generate Kick Velocity for Black Holes and Neutron Stars
* FULL KICK FOR BOTH NS AND BH
      elseif((kw.eq.14.or.kw.eq.13).and.(bhflag.eq.4))then
         mrem = m1n
         vk = vk
         vk2 = vk*vk
***
      elseif((kw.eq.14.and.bhflag.eq.0).or.kw.lt.0)then
         vk2 = 0.d0
         vk = 0.d0
         if(kw.lt.0) kw = 13
*
      else
         write(66,*)'We have a problem in kicks.f'
      endif
***
*     rescale base on new vk
      do k = 1,3
         v(k) = v(k)*vk/v(4)
      end do
      v(4) = vk
*
*  Correction (Noticed by Peter & Mirek)
      u1 = RAN3(idum)
      u2 = RAN3(idum)
      theta = twopi*u2
***
      sphi = -1.d0 + 2.d0*u1
      phi = ASIN(sphi)
      cphi = COS(phi)
      stheta = SIN(theta)
      ctheta = COS(theta)
      if(sep.le.0.d0.or.ecc.lt.0.d0) goto 90
*
* Determine the magnitude of the new relative velocity.
* (NG 12/17: there was a wrong sign in the brackets)
      vn2 = vk2+vr2-2.d0*vk*vr*(ctheta*cphi*salpha+stheta*cphi*calpha)
* Calculate the new semi-major axis.
      sep = 2.d0/r - vn2/(gmrkm*(m1n+m2))
      sep = 1.d0/sep
*     if(sep.le.0.d0)then
*        ecc = 1.1d0
*        goto 90
*     endif
* Determine the magnitude of the cross product of the separation vector
* and the new relative velocity.
      v1 = vk2*sphi*sphi
      v2 = (vk*ctheta*cphi-vr*salpha)**2
      hn2 = r*r*(v1 + v2)
* Calculate the new eccentricity.
      ecc2 = 1.d0 - hn2/(gmrkm*sep*(m1n+m2))
      ecc2 = MAX(ecc2,0.d0)
      ecc = SQRT(ecc2)
* Calculate the new orbital angular momentum taking care to convert
* hn to units of Rsun^2/yr.
      jorb = (m1n*m2/(m1n+m2))*SQRT(hn2)*(yearsc/rsunkm)
* Determine the angle between the new and old orbital angular
* momentum vectors.
      cmu = (vr*salpha-vk*ctheta*cphi)/SQRT(v1 + v2)
      mu = ACOS(cmu)
* Calculate the components of the velocity of the new centre-of-mass.
 90   continue
      if(ecc.le.1.0)then
* Calculate the components of the velocity of the new centre-of-mass.
         mx1 = vk*m1n/(m1n+m2)
         mx2 = vr*(m1-m1n)*m2/((m1n+m2)*(m1+m2))
         vs(1) = mx1*ctheta*cphi + mx2*salpha
         vs(2) = mx1*stheta*cphi + mx2*calpha
         vs(3) = mx1*sphi
      else
* Calculate the relative hyperbolic velocity at infinity (simple method).
         sep = r/(ecc-1.d0)
*        cmu = SQRT(ecc-1.d0)
*        mu = ATAN(cmu)
         mu = ACOS(1.d0/ecc)
         vr2 = gmrkm*(m1n+m2)/sep
         vr = SQRT(vr2)
         vs(1) = vr*SIN(mu)
         vs(2) = vr*COS(mu)
         vs(3) = 0.d0
         ecc = MIN(ecc,99.99d0)
      endif
*
 95   continue
* Save the information about the colculation of the natal kick velocity
* Save even the ecc to know what systems are still bound after the SN
*      if(kw.eq.13)then
*         WRITE(33,*)vk,phi,theta,ecc,sep,mfin,mrem,m1n
*      elseif(kw.eq.14)then
*         WRITE(44,*)vk,phi,theta,ecc,sep,mfin,mrem,m1n
*      endif
*
      RETURN
      END
***

max_stars_repo_name: michelamapelli/PeTar

max_issues_repo_name: michelamapelli/PeTar

max_forks_repo_name: michelamapelli/PeTar


text:
module type_function_status
#include <messenger.h>
    use mod_constants,                  only: BOUNDARY_ADVECTIVE_FLUX
    use type_mesh,                  only: mesh_t
    use type_equationset_function_data, only: equationset_function_data_t
    use type_function_status_data,      only: function_status_data_t
    use type_face_info,                 only: face_info_t
    use type_function_info,             only: function_info_t
    implicit none

    !> Class is used to manage the contributions and linearizations of functions.
    !!
    !! Domain component contains arrays of logicals that indicate the status of a given function.
    !!
    !! So, in this way, we can compute a flux on a face and add the contribution to both elements
    !! by registering that the function contribution was added. This container holds the registration
    !! information as arrays of logicals.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    type, public :: function_status_t

       type(function_status_data_t), allocatable :: dom(:)

    contains

        procedure   :: init

        procedure   :: compute_function             !< Check if a function has been executed
        procedure   :: linearize_function           !< Check if a function has been linearized

        procedure   :: compute_function_equation    !< Check whether a function has contributed to a particular equation
        procedure   :: linearize_function_equation  !< Check whether a function contribution to a particualar equation has been linearized

        procedure   :: register_function_computed   !< Register a function as computed
        procedure   :: register_function_linearized !< Register a function as linearized

        procedure   :: clear                        !< Reset function status data

    end type function_status_t
    !*****************************************************************************************

contains

    !>  Initialize class storage.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  mesh            Array of mesh_t instances.
    !!  @param[in]  function_data   Contains the number of each function type for every domain
    !!
    !-----------------------------------------------------------------------------------------
    subroutine init(self, mesh, function_data)
        class(function_status_t),           intent(inout)   :: self
        type(mesh_t),                   intent(in)      :: mesh
        type(equationset_function_data_t),  intent(in)      :: function_data(:)

        integer :: idom, ndom, ierr

        !
        ! Deallocate storage if necessary in case this is being called as a
        ! reinitialization routine.
        !
        if (allocated(self%dom)) deallocate(self%dom)

        ndom = mesh%ndomains()
        allocate(self%dom(ndom), stat=ierr)
        if (ierr /= 0) call AllocationError

        do idom = 1,ndom
            call self%dom(idom)%init( mesh%domain(idom), function_data(idom) )
        end do

    end subroutine init
    !*****************************************************************************************

    !>  Call to determine if a function needs computed for a given face, element.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  face_info   Container for face information. Location in mesh, type, etc.
    !!  @param[in]  flux_info   Container for function information. idiff, idonor, iflux, type
    !!  @return     res         True if function needs computed. False if already computed.
    !!
    !-----------------------------------------------------------------------------------------
    function compute_function(self, face_info, function_info) result(res)
        class(function_status_t),   intent(in)              :: self
        type(face_info_t),          intent(in)              :: face_info
        type(function_info_t),      intent(in)              :: function_info

        logical :: function_already_computed
        logical :: res

        associate(idomain_l => face_info%idomain_l,  ielement_l => face_info%ielement_l, iface => face_info%iface, &
                  type => function_info%type, ifcn => function_info%ifcn )

        function_already_computed = self%dom(idomain_l)%function_computed(type,ielement_l,iface,ifcn)

        if ( function_already_computed ) then
            res = .false.
        else
            res = .true.
        end if

        end associate

    end function compute_function
    !*****************************************************************************************

    !>  Call to determine if a particular equation has been contributed to from a given function.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  face_info   Container for face information. Location in mesh, type, etc.
    !!  @param[in]  flux_info   Container for function information. idiff, idonor, iflux, type
    !!  @param[in]  ieqn        Equation index within the function definition
    !!  @return     res         True if function needs computed. False if already computed.
    !!
    !-----------------------------------------------------------------------------------------
    function compute_function_equation(self, face_info, function_info, ieqn) result(res)
        class(function_status_t),   intent(in)              :: self
        type(face_info_t),          intent(in)              :: face_info
        type(function_info_t),      intent(in)              :: function_info
        integer(ik),                intent(in)              :: ieqn

        logical :: function_already_computed
        logical :: res

        associate(idomain_l => face_info%idomain_l,  ielement_l => face_info%ielement_l, iface => face_info%iface, &
                  type => function_info%type, ifcn => function_info%ifcn )

        function_already_computed = self%dom(idomain_l)%function_equation_computed(type,ielement_l,iface,ifcn,ieqn)

        if ( function_already_computed ) then
            res = .false.
        else
            res = .true.
        end if

        end associate

    end function compute_function_equation
    !*****************************************************************************************

    !>  Call to determine if a function needs linearized for a given face, element, seed.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  face_info   Container for face information. Location in mesh, type, etc.
    !!  @param[in]  flux_info   Container for function information. idiff, idonor, iflux, type
    !!  @return     res         True if function needs linearized. False if already linearized.
    !!
    !-----------------------------------------------------------------------------------------
    function linearize_function(self, face_info, function_info) result(res)
        class(function_status_t),   intent(in)  :: self
        type(face_info_t),          intent(in)  :: face_info
        type(function_info_t),      intent(in)  :: function_info

        logical :: function_already_linearized
        logical :: differentiating, res

        associate(idomain_l => face_info%idomain_l,  ielement_l => face_info%ielement_l, iface => face_info%iface, &
                  type => function_info%type, ifcn => function_info%ifcn,  idiff => function_info%idiff )

        differentiating = (idiff /= 0)

        if (differentiating) then

            function_already_linearized = self%dom(idomain_l)%function_linearized(type,ielement_l,iface,ifcn,idiff)
            if ( function_already_linearized ) then
                res = .false.
            else
                res = .true.
            end if

        else

            res = .false.

        end if

        end associate

    end function linearize_function
    !*****************************************************************************************

    !>  Call to determine if a particular equation has been linearized from a given function.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  face_info   Container for face information. Location in mesh, type, etc.
    !!  @param[in]  flux_info   Container for function information. idiff, idonor, iflux, type
    !!  @return     res         True if function needs linearized. False if already linearized.
    !!
    !-----------------------------------------------------------------------------------------
    function linearize_function_equation(self, face_info, function_info, ieqn) result(res)
        class(function_status_t),   intent(in)  :: self
        type(face_info_t),          intent(in)  :: face_info
        type(function_info_t),      intent(in)  :: function_info
        integer(ik),                intent(in)  :: ieqn

        logical :: function_already_linearized
        logical :: differentiating, res

        associate(idomain_l => face_info%idomain_l,  ielement_l => face_info%ielement_l, iface => face_info%iface, &
                  type => function_info%type, ifcn => function_info%ifcn,  idiff => function_info%idiff )

        differentiating = (idiff /= 0)

        if (differentiating) then

            function_already_linearized = self%dom(idomain_l)%function_equation_linearized(type,ielement_l,iface,ifcn,idiff,ieqn)
            if ( function_already_linearized ) then
                res = .false.
            else
                res = .true.
            end if

        else

            res = .false.

        end if

        end associate

    end function linearize_function_equation
    !*****************************************************************************************

    !>  Register a function as having been computed for a given element/face.
    !!
    !!  This sets both the 'function_computed' and 'function_equation_computed' logicals.
    !!  So, in this way, because this gets called for every equation, the 'function_computed' value
    !!  will get set multiple times. Once for every equation. That is okay.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  face_info       Contains face indices
    !!  @param[in]  function_info   Contains function indices, type, linearization index
    !!  @param[in]  ieqn            Index of equation being registered
    !!
    !!
    !------------------------------------------------------------------------------------------
    subroutine register_function_computed(self,face_info,function_info,ieqn)
        class(function_status_t),   intent(inout)   :: self
        type(face_info_t),          intent(in)      :: face_info
        type(function_info_t),      intent(in)      :: function_info
        integer(ik),                intent(in)      :: ieqn

        associate(idomain_l => face_info%idomain_l,  ielement_l => face_info%ielement_l, iface => face_info%iface, &
                  type => function_info%type, ifcn => function_info%ifcn)

            self%dom(idomain_l)%function_computed(type,ielement_l,iface,ifcn) = .true.
            self%dom(idomain_l)%function_equation_computed(type,ielement_l,iface,ifcn,ieqn) = .true.

        end associate
    end subroutine register_function_computed
    !******************************************************************************************

    !>  Register a function as having been linearized for a given element/face/blk.
    !!
    !!  This sets both the 'function_linearized' and 'function_equation_linearized' logicals.
    !!  So, in this way, because this gets called for every equation, the 'function_linearized' value
    !!  will get set multiple times. Once for every equation. That is okay.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  face_info       Contains face indices
    !!  @param[in]  function_info   Contains function indices, type, linearization index
    !!  @param[in]  ieqn            Index of equation being registered
    !!
    !------------------------------------------------------------------------------------------
    subroutine register_function_linearized(self,face_info,function_info,ieqn)
        class(function_status_t),   intent(inout)   :: self
        type(face_info_t),          intent(in)      :: face_info
        type(function_info_t),      intent(in)      :: function_info
        integer(ik),                intent(in)      :: ieqn

        associate(idomain_l => face_info%idomain_l,  ielement_l => face_info%ielement_l, iface => face_info%iface, &
                  type => function_info%type, ifcn => function_info%ifcn, idiff => function_info%idiff )

            self%dom(idomain_l)%function_linearized(type,ielement_l,iface,ifcn,idiff) = .true.
            self%dom(idomain_l)%function_equation_linearized(type,ielement_l,iface,ifcn,idiff,ieqn) = .true.

        end associate
    end subroutine register_function_linearized
    !*****************************************************************************************

    !>  Call clear on each domain local function status data. Should be called before every
    !!  evaluation of the right-hand side.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!
    !----------------------------------------------------------------------------------------
    subroutine clear(self)
        class(function_status_t),   intent(inout)   :: self

        integer :: idom

        do idom = 1,size(self%dom)
            call self%dom(idom)%clear()
        end do

    end subroutine clear
    !****************************************************************************************

end module type_function_status

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
PROGRAM test

  USE hdf5_wrapper
  IMPLICIT NONE
  CHARACTER(LEN=256)          :: fname
  INTEGER                     :: ifile
  INTEGER*4                   :: dozen,intattr
  integer*4, dimension(5,3)   :: i4attr
  real, dimension(5,3)        :: r4attr
  real*8, dimension(3,100000) :: r8arr
  integer*4, dimension(10,20) :: i4arr
  INTEGER :: i,j,k
  CHARACTER(LEN=20) :: str(5,5)
  integer, dimension(7) :: start,count
  integer, dimension(10) :: i1d
  character(len=50) :: dtype
  integer :: size
  integer :: nsets
  character(len=200), dimension(100) :: set_name

  fname = "./test.hdf5"

  HDF_VERBOSITY=2

  CALL hdf5_open_file(ifile, fname, readonly=.true.)

  call hdf5_read_data(ifile,"Dozen",Dozen)
  write(*,*)'One Dozen is : ',dozen

  call hdf5_read_data(ifile,"/Arrays/Real/R8_Array",r8arr(1:1,:))
  write(*,*)'R8_Array(1:10) = ',r8arr(1:1,1:10)

  call hdf5_read_attribute(ifile,"Arrays/I4_Attribute",i4attr)
  write(*,*)'The 2d integer*4 attribute is : ',i4attr

  call hdf5_read_attribute(ifile,"Arrays/I4_Scalar",intattr)
  write(*,*)'The scalar integer*4 attribute is : ',intattr

  call hdf5_read_attribute(ifile,"Arrays/R4_Attribute",r4attr)
  write(*,*)'The 2d real*4 attribute is : ',r4attr

  call hdf5_read_data(ifile,"Arrays/Test-String",str)
  DO i=1,2
     DO j=1,2
        WRITE(*,*)'String: ',i,j,TRIM(str(i,j))
     END DO
  END DO

  call hdf5_read_data(ifile,"/Arrays/I4_Array",i4arr)
  write(*,*)'I4_Array:'
  write(*,*)i4arr(1,1),i4arr(2,1),i4arr(3,1),i4arr(4,1),i4arr(5,1)
  write(*,*)i4arr(1,2),i4arr(2,2),i4arr(3,2),i4arr(4,2),i4arr(5,2)
  write(*,*)i4arr(1,3),i4arr(2,3),i4arr(3,3),i4arr(4,3),i4arr(5,3)

  i4arr=i4arr*0
  start(1)=3
  start(2)=2
  count(1)=3
  count(2)=2
  call hdf5_read_data(ifile,"/Arrays/I4_Array",i4arr,start=start,count=count)
  write(*,*)'Subsection from I4_Array:'
  write(*,*)i4arr(1,1),i4arr(2,1),i4arr(3,1),i4arr(4,1),i4arr(5,1)
  write(*,*)i4arr(1,2),i4arr(2,2),i4arr(3,2),i4arr(4,2),i4arr(5,2)
  write(*,*)i4arr(1,3),i4arr(2,3),i4arr(3,3),i4arr(4,3),i4arr(5,3)

  start(1)=4
  start(2)=2
  count(1)=1
  count(2)=1
  call hdf5_read_data(ifile,"/Arrays/I4_Array",i,start=start,count=count)
  write(*,*)'One element from I4_Array:'
  write(*,*)i

  start(1)=2
  start(2)=2
  count(1)=3
  count(2)=1
  i4arr=i4arr*0
  call hdf5_read_data(ifile,"/Arrays/I4_Array",i4arr(3,1:3),start=start,count=count)
  write(*,*)'Part of a row from I4_Array read into the third column of a 2D array:'
  write(*,*)i4arr(1,1),i4arr(2,1),i4arr(3,1),i4arr(4,1),i4arr(5,1)
  write(*,*)i4arr(1,2),i4arr(2,2),i4arr(3,2),i4arr(4,2),i4arr(5,2)
  write(*,*)i4arr(1,3),i4arr(2,3),i4arr(3,3),i4arr(4,3),i4arr(5,3)

  start(1)=2
  start(2)=1
  count(1)=1
  count(2)=3
  i4arr=i4arr*0
  call hdf5_read_data(ifile,"/Arrays/I4_Array",i4arr(1:3,3),start=start,count=count)
  write(*,*)'A column from I4_Array read into the third row of a 2D array:'
  write(*,*)i4arr(1,1),i4arr(2,1),i4arr(3,1),i4arr(4,1),i4arr(5,1)
  write(*,*)i4arr(1,2),i4arr(2,2),i4arr(3,2),i4arr(4,2),i4arr(5,2)
  write(*,*)i4arr(1,3),i4arr(2,3),i4arr(3,3),i4arr(4,3),i4arr(5,3)

  CALL hdf5_get_type(ifile,"Dozen",dtype,size)
  write(*,*)'Type and size of "Dozen" are ',trim(dtype),size

  CALL hdf5_get_type(ifile,"/Arrays/Real/R8_Array",dtype,size)
  write(*,*)'Type and size of "R8_Array" are ',trim(dtype),size

  write(*,*)'Datasets in root group:'
  call hdf5_list_datasets(ifile,"/",nsets,set_name)
  do i = 1, nsets, 1
     write(*,*)"  ",trim(set_name(i))
  end do

  write(*,*)'Datasets in group "Arrays":'
  call hdf5_list_datasets(ifile,"/Arrays",nsets,set_name)
  do i = 1, nsets, 1
     write(*,*)"  ",trim(set_name(i))
  end do

  write(*,*)'Groups in root group:'
  call hdf5_list_groups(ifile,"/",nsets,set_name)
  do i = 1, nsets, 1
     write(*,*)"  ",trim(set_name(i))
  end do

  write(*,*)'Attributes in group "Arrays":'
  call hdf5_list_attributes(ifile,"/Arrays",nsets,set_name)
  do i = 1, nsets, 1
     write(*,*)"  ",trim(set_name(i))
  end do

  CALL hdf5_close_file(ifile)

END PROGRAM test

max_stars_repo_name: galtay/urchin

max_issues_repo_name: galtay/urchin

max_forks_repo_name: galtay/urchin


text:
*> \brief \b DLAQZ4
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAQZ4 + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE DLAQZ4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
*     $    NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ,
*     $    QC, LDQC, ZC, LDZC, WORK, LWORK, INFO )
*      IMPLICIT NONE
*
*      Function arguments
*      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
*      INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
*     $    NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
*
*      DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
*     $    Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ),
*     $    WORK( * ), SR( * ), SI( * ), SS( * )
*
*      INTEGER, INTENT( OUT ) :: INFO
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAQZ4 Executes a single multishift QZ sweep
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ILSCHUR
*> \verbatim
*>          ILSCHUR is LOGICAL
*>              Determines whether or not to update the full Schur form
*> \endverbatim
*> \param[in] ILQ
*> \verbatim
*>          ILQ is LOGICAL
*>              Determines whether or not to update the matrix Q
*> \endverbatim
*>
*> \param[in] ILZ
*> \verbatim
*>          ILZ is LOGICAL
*>              Determines whether or not to update the matrix Z
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices A, B, Q, and Z.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*> \endverbatim
*>
*> \param[in] NSHIFTS
*> \verbatim
*>          NSHIFTS is INTEGER
*>          The desired number of shifts to use
*> \endverbatim
*>
*> \param[in] NBLOCK_DESIRED
*> \verbatim
*>          NBLOCK_DESIRED is INTEGER
*>          The desired size of the computational windows
*> \endverbatim
*>
*> \param[in] SR
*> \verbatim
*>          SR is DOUBLE PRECISION array. SR contains
*>          the real parts of the shifts to use.
*> \endverbatim
*>
*> \param[in] SI
*> \verbatim
*>          SI is DOUBLE PRECISION array. SI contains
*>          the imaginary parts of the shifts to use.
*> \endverbatim
*>
*> \param[in] SS
*> \verbatim
*>          SS is DOUBLE PRECISION array. SS contains
*>          the scale of the shifts to use.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, N)
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, N)
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[in,out] QC
*> \verbatim
*>          QC is DOUBLE PRECISION array, dimension (LDQC, NBLOCK_DESIRED)
*> \endverbatim
*>
*> \param[in] LDQC
*> \verbatim
*>          LDQC is INTEGER
*> \endverbatim
*>
*> \param[in,out] ZC
*> \verbatim
*>          ZC is DOUBLE PRECISION array, dimension (LDZC, NBLOCK_DESIRED)
*> \endverbatim
*>
*> \param[in] LDZC
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK.  LWORK >= max(1,N).
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Thijs Steel, KU Leuven
*
*> \date May 2020
*
*> \ingroup doubleGEcomputational
*>
*  =====================================================================
      SUBROUTINE DLAQZ4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
     $                   NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q,
     $                   LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK,
     $                   INFO )
      IMPLICIT NONE

*     Function arguments
      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
      INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
     $         NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC

      DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
     $                  Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ),
     $                  ZC( LDZC, * ), WORK( * ), SR( * ), SI( * ),
     $                  SS( * )

      INTEGER, INTENT( OUT ) :: INFO

*     Parameters
      DOUBLE PRECISION :: ZERO, ONE, HALF
      PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )

*     Local scalars
      INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
     $           ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
      DOUBLE PRECISION :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
*
*     External functions
      EXTERNAL :: XERBLA, DGEMM, DLAQZ1, DLAQZ2, DLASET, DLARTG, DROT,
     $            DLACPY

      INFO = 0
      IF ( NBLOCK_DESIRED .LT. NSHIFTS+1 ) THEN
         INFO = -8
      END IF
      IF ( LWORK .EQ.-1 ) THEN
*        workspace query, quick return
         WORK( 1 ) = N*NBLOCK_DESIRED
         RETURN
      ELSE IF ( LWORK .LT. N*NBLOCK_DESIRED ) THEN
         INFO = -25
      END IF

      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DLAQZ4', -INFO )
         RETURN
      END IF

*     Executable statements

      IF ( NSHIFTS .LT. 2 ) THEN
         RETURN
      END IF

      IF ( ILO .GE. IHI ) THEN
         RETURN
      END IF

      IF ( ILSCHUR ) THEN
         ISTARTM = 1
         ISTOPM = N
      ELSE
         ISTARTM = ILO
         ISTOPM = IHI
      END IF

*     Shuffle shifts into pairs of real shifts and pairs
*     of complex conjugate shifts assuming complex
*     conjugate shifts are already adjacent to one
*     another

      DO I = 1, NSHIFTS-2, 2
         IF( SI( I ).NE.-SI( I+1 ) ) THEN
*
            SWAP = SR( I )
            SR( I ) = SR( I+1 )
            SR( I+1 ) = SR( I+2 )
            SR( I+2 ) = SWAP

            SWAP = SI( I )
            SI( I ) = SI( I+1 )
            SI( I+1 ) = SI( I+2 )
            SI( I+2 ) = SWAP

            SWAP = SS( I )
            SS( I ) = SS( I+1 )
            SS( I+1 ) = SS( I+2 )
            SS( I+2 ) = SWAP
         END IF
      END DO

*     NSHFTS is supposed to be even, but if it is odd,
*     then simply reduce it by one.  The shuffle above
*     ensures that the dropped shift is real and that
*     the remaining shifts are paired.

      NS = NSHIFTS-MOD( NSHIFTS, 2 )
      NPOS = MAX( NBLOCK_DESIRED-NS, 1 )

*     The following block introduces the shifts and chases
*     them down one by one just enough to make space for
*     the other shifts. The near-the-diagonal block is
*     of size (ns+1) x ns.

      CALL DLASET( 'FULL', NS+1, NS+1, ZERO, ONE, QC, LDQC )
      CALL DLASET( 'FULL', NS, NS, ZERO, ONE, ZC, LDZC )

      DO I = 1, NS, 2
*        Introduce the shift
         CALL DLAQZ1( A( ILO, ILO ), LDA, B( ILO, ILO ), LDB, SR( I ),
     $                SR( I+1 ), SI( I ), SS( I ), SS( I+1 ), V )

         TEMP = V( 2 )
         CALL DLARTG( TEMP, V( 3 ), C1, S1, V( 2 ) )
         CALL DLARTG( V( 1 ), V( 2 ), C2, S2, TEMP )

         CALL DROT( NS, A( ILO+1, ILO ), LDA, A( ILO+2, ILO ), LDA, C1,
     $              S1 )
         CALL DROT( NS, A( ILO, ILO ), LDA, A( ILO+1, ILO ), LDA, C2,
     $              S2 )
         CALL DROT( NS, B( ILO+1, ILO ), LDB, B( ILO+2, ILO ), LDB, C1,
     $              S1 )
         CALL DROT( NS, B( ILO, ILO ), LDB, B( ILO+1, ILO ), LDB, C2,
     $              S2 )
         CALL DROT( NS+1, QC( 1, 2 ), 1, QC( 1, 3 ), 1, C1, S1 )
         CALL DROT( NS+1, QC( 1, 1 ), 1, QC( 1, 2 ), 1, C2, S2 )

*        Chase the shift down
         DO J = 1, NS-1-I

            CALL DLAQZ2( .TRUE., .TRUE., J, 1, NS, IHI-ILO+1, A( ILO,
     $                   ILO ), LDA, B( ILO, ILO ), LDB, NS+1, 1, QC,
     $                   LDQC, NS, 1, ZC, LDZC )

         END DO

      END DO

*     Update the rest of the pencil

*     Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
*     from the left with Qc(1:ns+1,1:ns+1)'
      SHEIGHT = NS+1
      SWIDTH = ISTOPM-( ILO+NS )+1
      IF ( SWIDTH > 0 ) THEN
         CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
     $               A( ILO, ILO+NS ), LDA, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ILO,
     $                ILO+NS ), LDA )
         CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
     $               B( ILO, ILO+NS ), LDB, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ILO,
     $                ILO+NS ), LDB )
      END IF
      IF ( ILQ ) THEN
         CALL DGEMM( 'N', 'N', N, SHEIGHT, SHEIGHT, ONE, Q( 1, ILO ),
     $               LDQ, QC, LDQC, ZERO, WORK, N )
         CALL DLACPY( 'ALL', N, SHEIGHT, WORK, N, Q( 1, ILO ), LDQ )
      END IF

*     Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
*     from the right with Zc(1:ns,1:ns)
      SHEIGHT = ILO-1-ISTARTM+1
      SWIDTH = NS
      IF ( SHEIGHT > 0 ) THEN
         CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, A( ISTARTM,
     $               ILO ), LDA, ZC, LDZC, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ISTARTM,
     $                ILO ), LDA )
         CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, B( ISTARTM,
     $               ILO ), LDB, ZC, LDZC, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ISTARTM,
     $                ILO ), LDB )
      END IF
      IF ( ILZ ) THEN
         CALL DGEMM( 'N', 'N', N, SWIDTH, SWIDTH, ONE, Z( 1, ILO ), LDZ,
     $               ZC, LDZC, ZERO, WORK, N )
         CALL DLACPY( 'ALL', N, SWIDTH, WORK, N, Z( 1, ILO ), LDZ )
      END IF

*     The following block chases the shifts down to the bottom
*     right block. If possible, a shift is moved down npos
*     positions at a time

      K = ILO
      DO WHILE ( K < IHI-NS )
         NP = MIN( IHI-NS-K, NPOS )
*        Size of the near-the-diagonal block
         NBLOCK = NS+NP
*        istartb points to the first row we will be updating
         ISTARTB = K+1
*        istopb points to the last column we will be updating
         ISTOPB = K+NBLOCK-1

         CALL DLASET( 'FULL', NS+NP, NS+NP, ZERO, ONE, QC, LDQC )
         CALL DLASET( 'FULL', NS+NP, NS+NP, ZERO, ONE, ZC, LDZC )

*        Near the diagonal shift chase
         DO I = NS-1, 0, -2
            DO J = 0, NP-1
*              Move down the block with index k+i+j-1, updating
*              the (ns+np x ns+np) block:
*              (k:k+ns+np,k:k+ns+np-1)
               CALL DLAQZ2( .TRUE., .TRUE., K+I+J-1, ISTARTB, ISTOPB,
     $                      IHI, A, LDA, B, LDB, NBLOCK, K+1, QC, LDQC,
     $                      NBLOCK, K, ZC, LDZC )
            END DO
         END DO

*        Update rest of the pencil

*        Update A(k+1:k+ns+np, k+ns+np:istopm) and
*        B(k+1:k+ns+np, k+ns+np:istopm)
*        from the left with Qc(1:ns+np,1:ns+np)'
         SHEIGHT = NS+NP
         SWIDTH = ISTOPM-( K+NS+NP )+1
         IF ( SWIDTH > 0 ) THEN
            CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC,
     $                  LDQC, A( K+1, K+NS+NP ), LDA, ZERO, WORK,
     $                  SHEIGHT )
            CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( K+1,
     $                   K+NS+NP ), LDA )
            CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC,
     $                  LDQC, B( K+1, K+NS+NP ), LDB, ZERO, WORK,
     $                  SHEIGHT )
            CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( K+1,
     $                   K+NS+NP ), LDB )
         END IF
         IF ( ILQ ) THEN
            CALL DGEMM( 'N', 'N', N, NBLOCK, NBLOCK, ONE, Q( 1, K+1 ),
     $                  LDQ, QC, LDQC, ZERO, WORK, N )
            CALL DLACPY( 'ALL', N, NBLOCK, WORK, N, Q( 1, K+1 ), LDQ )
         END IF

*        Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
*        from the right with Zc(1:ns+np,1:ns+np)
         SHEIGHT = K-ISTARTM+1
         SWIDTH = NBLOCK
         IF ( SHEIGHT > 0 ) THEN
            CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE,
     $                  A( ISTARTM, K ), LDA, ZC, LDZC, ZERO, WORK,
     $                  SHEIGHT )
            CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
     $                   A( ISTARTM, K ), LDA )
            CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE,
     $                  B( ISTARTM, K ), LDB, ZC, LDZC, ZERO, WORK,
     $                  SHEIGHT )
            CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
     $                   B( ISTARTM, K ), LDB )
         END IF
         IF ( ILZ ) THEN
            CALL DGEMM( 'N', 'N', N, NBLOCK, NBLOCK, ONE, Z( 1, K ),
     $                  LDZ, ZC, LDZC, ZERO, WORK, N )
            CALL DLACPY( 'ALL', N, NBLOCK, WORK, N, Z( 1, K ), LDZ )
         END IF

         K = K+NP

      END DO

*     The following block removes the shifts from the bottom right corner
*     one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).

      CALL DLASET( 'FULL', NS, NS, ZERO, ONE, QC, LDQC )
      CALL DLASET( 'FULL', NS+1, NS+1, ZERO, ONE, ZC, LDZC )

*     istartb points to the first row we will be updating
      ISTARTB = IHI-NS+1
*     istopb points to the last column we will be updating
      ISTOPB = IHI

      DO I = 1, NS, 2
*        Chase the shift down to the bottom right corner
         DO ISHIFT = IHI-I-1, IHI-2
            CALL DLAQZ2( .TRUE., .TRUE., ISHIFT, ISTARTB, ISTOPB, IHI,
     $                   A, LDA, B, LDB, NS, IHI-NS+1, QC, LDQC, NS+1,
     $                   IHI-NS, ZC, LDZC )
         END DO

      END DO

*     Update rest of the pencil

*     Update A(ihi-ns+1:ihi, ihi+1:istopm)
*     from the left with Qc(1:ns,1:ns)'
      SHEIGHT = NS
      SWIDTH = ISTOPM-( IHI+1 )+1
      IF ( SWIDTH > 0 ) THEN
         CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
     $               A( IHI-NS+1, IHI+1 ), LDA, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
     $                A( IHI-NS+1, IHI+1 ), LDA )
         CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
     $               B( IHI-NS+1, IHI+1 ), LDB, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
     $                B( IHI-NS+1, IHI+1 ), LDB )
      END IF
      IF ( ILQ ) THEN
         CALL DGEMM( 'N', 'N', N, NS, NS, ONE, Q( 1, IHI-NS+1 ), LDQ,
     $               QC, LDQC, ZERO, WORK, N )
         CALL DLACPY( 'ALL', N, NS, WORK, N, Q( 1, IHI-NS+1 ), LDQ )
      END IF

*     Update A(istartm:ihi-ns,ihi-ns:ihi)
*     from the right with Zc(1:ns+1,1:ns+1)
      SHEIGHT = IHI-NS-ISTARTM+1
      SWIDTH = NS+1
      IF ( SHEIGHT > 0 ) THEN
         CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, A( ISTARTM,
     $               IHI-NS ), LDA, ZC, LDZC, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ISTARTM,
     $                IHI-NS ), LDA )
         CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, B( ISTARTM,
     $               IHI-NS ), LDB, ZC, LDZC, ZERO, WORK, SHEIGHT )
         CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ISTARTM,
     $                IHI-NS ), LDB )
      END IF
      IF ( ILZ ) THEN
         CALL DGEMM( 'N', 'N', N, NS+1, NS+1, ONE, Z( 1, IHI-NS ), LDZ,
     $               ZC, LDZC, ZERO, WORK, N )
         CALL DLACPY( 'ALL', N, NS+1, WORK, N, Z( 1, IHI-NS ), LDZ )
      END IF

      END SUBROUTINE

max_stars_repo_name: msk/lapack

max_issues_repo_name: msk/lapack

max_forks_repo_name: msk/lapack


text:
        subroutine marqcoeferr(y,ndata,a,m,mfit,beta)

	implicit real*4 (a-h,o-z)

	include 'PARAMETER' !por kn,kl y kld
        parameter (mmax=16*kn+2,kld4=4*kld,kldn=mmax*kld4)
        dimension y(ndata),beta(*)
        dimension dyda(mmax),a(*),ymod(kld4),ymodobs(kld4),dvcal(kldn)
        real*4 sigreal(kld4)
	integer m(*),npos(kld4)

	common/uvesalida/ymodobs
        common/derivchi/dvcal  ! deriv. chi^2 para los errores
        common/primera2/ntotal,ntotal4,ists
        common/posiciones/npos
	common/sigrealchi/sigreal,chireal,sumsq

      do 12 j=1,mfit
         beta(j)=0.
12    continue

      call fperfil2err(m,a,ymod,dvcal)
      do i=1,ndata
         ymodobs(i)=ymod(npos(i))
      end do

      sig0=0.	!Diferencia cuadratica entre perfiles obs. y cal.

      do i=1,ndata
         dy=(y(i)-ymodobs(i))/sigreal(i)
         sig0=sig0+dy*dy
      enddo

      do 15 i=1,ndata

	 do j=1,mfit
	    jj=(j-1)*ntotal4+npos(i)
	    dyda(j)=dvcal(jj)
	 end do

         do 14 j=1,mfit
            ww=dyda(j)/sigreal(i)
            beta(j)=beta(j)+ww*ww
14       continue

15    continue

      do 16 j=1,mfit
         beta(j)=beta(j)/sig0
16    continue

      return
      end

max_stars_repo_name: cdiazbas/SIRcode

max_issues_repo_name: cdiazbas/SIRcode

max_forks_repo_name: cdiazbas/SIRcode


text:
module matrix_mod
implicit none

contains

    subroutine get_matrix_dimensions(file_name, matrix_dset_name, &
                                     dims)
        use hdf5
        use, intrinsic :: iso_fortran_env, only : dp => REAL64, error_unit
        implicit none
        integer, parameter :: matrix_rank = 2
        character(len=*), intent(in) :: file_name, matrix_dset_name
        integer, dimension(matrix_rank) :: dims
        integer(hid_t) :: file_id, matrix_dset_id, dset_prop_id, &
                          matrix_mspace_id, matrix_fspace_id
        integer(hsize_t), dimension(matrix_rank) :: matrix_dims, &
                                                max_matrix_dims
        integer :: error
        call h5open_f(error)

        ! open the HDF5 file, overwriting it if it alread exists
        call h5fopen_f(file_name, H5F_ACC_RDONLY_F, file_id, error)
        if (error /= 0) then
            write (unit=error_unit, fmt="(2A)") "# error opening file ", &
                                                file_name
            stop
        end if

        ! open the matrix dataset, and get the matrix file dataspace
        call h5dopen_f(file_id, matrix_dset_name, matrix_dset_id, error)
        call h5dget_space_f(matrix_dset_id, matrix_fspace_id, error)

        ! retrieve the dimensions of the matrix, and allocate the matrix
        call h5sget_simple_extent_dims_f(matrix_fspace_id, matrix_dims, &
                                         max_matrix_dims, error)
        ! close the matrix dataset
        call h5dclose_f(matrix_dset_id, error)

        ! close HDF5 file
        call h5fclose_f(file_id, error)

        ! finalize Fortran HDF5 interface
        call h5close_f(error)
        dims = int(matrix_dims)
    end subroutine get_matrix_dimensions

    subroutine read_matrix(file_name, matrix_dset_name, matrix)
        use hdf5
        use, intrinsic :: iso_fortran_env, only : dp => REAL64, error_unit
        implicit none
        character(len=*), intent(in) :: file_name, matrix_dset_name
        real(kind=dp), allocatable, dimension(:, :), intent(out) :: matrix
        integer(hid_t) :: file_id, matrix_dset_id, dset_prop_id, &
                          matrix_mspace_id, matrix_fspace_id
        integer, parameter :: matrix_rank = 2
        integer(hsize_t), dimension(matrix_rank) :: matrix_dims, &
                                                    max_matrix_dims, &
                                                    chunk_dims
        integer :: error
        ! initialize Fortran interface
        call h5open_f(error)

        ! open the HDF5 file, overwriting it if it alread exists
        call h5fopen_f(file_name, H5F_ACC_RDONLY_F, file_id, error)
        if (error /= 0) then
            write (unit=error_unit, fmt="(2A)") "# error opening file ", &
                                                file_name
            stop
        end if

        ! open the matrix dataset, and get the matrix file dataspace
        call h5dopen_f(file_id, matrix_dset_name, matrix_dset_id, error)
        call h5dget_space_f(matrix_dset_id, matrix_fspace_id, error)

        ! retrieve the dimensions of the matrix, and allocate the matrix
        call h5sget_simple_extent_dims_f(matrix_fspace_id, matrix_dims, &
                                         max_matrix_dims, error)
        allocate(matrix(matrix_dims(1), matrix_dims(2)), stat=error)
        if (error /= 0) then
            write (unit=error_unit, fmt="(A)") &
                "matrix chunk can not be allocated"
            stop
        end if

        ! create the matrix memory dataspace
        call h5screate_simple_f(matrix_rank, matrix_dims, &
                                matrix_mspace_id, error)

        ! read the data
        call h5dread_f(matrix_dset_id, H5T_NATIVE_DOUBLE, matrix, &
                       matrix_dims, error, mem_space_id=matrix_mspace_id, &
                       file_space_id=matrix_fspace_id)

        ! close dataspaces
        call h5sclose_f(matrix_mspace_id, error)
        call h5sclose_f(matrix_fspace_id, error)

        ! close the matrix dataset
        call h5dclose_f(matrix_dset_id, error)

        ! close HDF5 file
        call h5fclose_f(file_id, error)

        ! finalize Fortran HDF5 interface
        call h5close_f(error)
    end subroutine read_matrix

    subroutine print_matrix(matrix)
        use, intrinsic :: iso_fortran_env, only : dp => REAL64
        implicit none
        real(kind=dp), dimension(:, :), intent(in) :: matrix
        integer i
        character(len=50) :: format_str
        write (format_str, "('(', I0, 'F8.3)')") size(matrix, 2)
        do i = 1, size(matrix, 1)
            print format_str, matrix(i, :)
        end do
    end subroutine print_matrix

end module matrix_mod

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


! A proxy which divides many calls over few instances.
!
! Put this component between a driver and a set of models, or between a
! macro model and a set of micro models. It will let the driver or macro-
! model submit as many calls as it wants, and divide them over the available
! (micro)model instances in a round-robin fashion.
!
! Assumes a fixed number of micro-model instances.

program load_balancer
    use ymmsl
    use libmuscle
    implicit none

    type(LIBMUSCLE_PortsDescription) :: ports
    type(LIBMUSCLE_Instance) :: instance

    type(LIBMUSCLE_Message) :: rmsg
    type(LIBMUSCLE_DataConstRef) :: rdata

    integer :: started, done, num_calls, num_workers

    ports = LIBMUSCLE_PortsDescription_create()
    call LIBMUSCLE_PortsDescription_add(ports, YMMSL_Operator_F_INIT, 'front_in[]')
    call LIBMUSCLE_PortsDescription_add(ports, YMMSL_Operator_O_I, 'back_out[]')
    call LIBMUSCLE_PortsDescription_add(ports, YMMSL_Operator_S, 'back_in[]')
    call LIBMUSCLE_PortsDescription_add(ports, YMMSL_Operator_O_F, 'front_out[]')
    instance = LIBMUSCLE_Instance_create(ports)
    call LIBMUSCLE_PortsDescription_free(ports)

    do while (LIBMUSCLE_Instance_reuse_instance(instance, .false.))
        ! F_INIT
        started = 0
        done = 0

        num_calls = LIBMUSCLE_Instance_get_port_length(instance, 'front_in')
        num_workers = LIBMUSCLE_Instance_get_port_length(instance, 'back_out')

        call LIBMUSCLE_Instance_set_port_length(instance, 'front_out', num_calls)

        do while (done < num_calls)
            do while ((started - done < num_workers) .and. (started < num_calls))
                rmsg = LIBMUSCLE_Instance_receive_with_settings_on_slot(instance, 'front_in', started)
                call LIBMUSCLE_Instance_send(instance, 'back_out', rmsg, mod(started, num_workers))
                call LIBMUSCLE_Message_free(rmsg)
                started = started + 1
            end do

            rmsg = LIBMUSCLE_Instance_receive_with_settings_on_slot(instance, 'back_in', mod(done, num_workers))
            call LIBMUSCLE_Instance_send(instance, 'front_out', rmsg, done)
            call LIBMUSCLE_Message_free(rmsg)
            done = done + 1
        end do
    end do

    call LIBMUSCLE_Instance_free(instance)

end program load_balancer