text:
Subroutine hijsft(jp, jt, jout, ierror)
  Parameter (maxstr=150001)
  Common /hjcrdn/yp(3, 300), yt(3, 300)
  Common /hparnt/hipr1(100), ihpr2(50), hint1(100), ihnt2(50)
  Common /hijdat/hidat0(10, 10), hidat(10)
  Common /rndf77/nseed
  Common /hjjet1/npj(300), kfpj(300, 500), pjpx(300, 500), pjpy(300, 500), pjpz(300, 500), pjpe(300, 500), pjpm(300, 500), ntj(300), kftj(300, 500), pjtx(300, 500), pjty(300, 500), pjtz(300, 500), pjte(300, 500), pjtm(300, 500)
  Common /hstrng/nfp(300, 15), pp(300, 15), nft(300, 15), pt(300, 15)
  Common /dpmcm1/jjp, jjt, amp, amt, apx0, atx0, ampn, amtn, amp0, amt0, nfdp, nfdt, wp, wm, sw, xremp, xremt, dpkc1, dpkc2, pp11, pp12, pt11, pt12, ptp2, ptt2
  Common /dpmcm2/ndpm, kdpm(20, 2), pdpm1(20, 5), pdpm2(20, 5)
  Save
  ierror = 0
  jjp = jp
  jjt = jt
  ndpm = 0
  If (jp>ihnt2(1) .Or. jt>ihnt2(3)) Return
  epp = pp(jp, 4) + pp(jp, 3)
  epm = pp(jp, 4) - pp(jp, 3)
  etp = pt(jt, 4) + pt(jt, 3)
  etm = pt(jt, 4) - pt(jt, 3)
  wp = epp + etp
  wm = epm + etm
  sw = wp*wm
  If (wp<0.0 .Or. wm<0.0) Goto 1000
  If (jout==0) Then
    If (epp<0.0) Goto 1000
    If (epm<0.0) Goto 1000
    If (etp<0.0) Goto 1000
    If (etm<0.0) Goto 1000
    If (epp/(epm+0.01)<=etp/(etm+0.01)) Return
  End If
  ihnt2(11) = jp
  ihnt2(12) = jt
  miss = 0
  pkc1 = 0.0
  pkc2 = 0.0
  pkc11 = 0.0
  pkc12 = 0.0
  pkc21 = 0.0
  pkc22 = 0.0
  dpkc11 = 0.0
  dpkc12 = 0.0
  dpkc21 = 0.0
  dpkc22 = 0.0
  If (nfp(jp,10)==1 .Or. nft(jt,10)==1) Then
    If (nfp(jp,10)==1) Then
      phi1 = ulangl(pp(jp,10), pp(jp,11))
      ppjet = sqrt(pp(jp,10)**2+pp(jp,11)**2)
      pkc1 = ppjet
      pkc11 = pp(jp, 10)
      pkc12 = pp(jp, 11)
    End If
    If (nft(jt,10)==1) Then
      phi2 = ulangl(pt(jt,10), pt(jt,11))
      ptjet = sqrt(pt(jt,10)**2+pt(jt,11)**2)
      pkc2 = ptjet
      pkc21 = pt(jt, 10)
      pkc22 = pt(jt, 11)
    End If
    If (ihpr2(4)>0 .And. ihnt2(1)>1 .And. ihnt2(3)>1) Then
      If (nfp(jp,10)==0) Then
        phi = -phi2
      Else If (nft(jt,10)==0) Then
        phi = phi1
      Else
        phi = (phi1+phi2-hipr1(40))/2.0
      End If
      bx = hint1(19)*cos(hint1(20))
      by = hint1(19)*sin(hint1(20))
      xp0 = yp(1, jp)
      yp0 = yp(2, jp)
      xt0 = yt(1, jt) + bx
      yt0 = yt(2, jt) + by
      r1 = max(1.2*ihnt2(1)**0.3333333, sqrt(xp0**2+yp0**2))
      r2 = max(1.2*ihnt2(3)**0.3333333, sqrt((xt0-bx)**2+(yt0-by)**2))
      If (abs(cos(phi))<1.0E-5) Then
        dd1 = r1
        dd2 = r1
        dd3 = abs(by+sqrt(r2**2-(xp0-bx)**2)-yp0)
        dd4 = abs(by-sqrt(r2**2-(xp0-bx)**2)-yp0)
        Goto 5
      End If
      bb = 2.0*sin(phi)*(cos(phi)*yp0-sin(phi)*xp0)
      cc = (yp0**2-r1**2)*cos(phi)**2 + xp0*sin(phi)*(xp0*sin(phi)-2.0*yp0*cos(phi))
      dd = bb**2 - 4.0*cc
      If (dd<0.0) Goto 10
      xx1 = (-bb+sqrt(dd))/2.0
      xx2 = (-bb-sqrt(dd))/2.0
      dd1 = abs((xx1-xp0)/cos(phi))
      dd2 = abs((xx2-xp0)/cos(phi))
      bb = 2.0*sin(phi)*(cos(phi)*(yt0-by)-sin(phi)*xt0) - 2.0*bx
      cc = (bx**2+(yt0-by)**2-r2**2)*cos(phi)**2 + xt0*sin(phi)*(xt0*sin(phi)-2.0*cos(phi)*(yt0-by)) - 2.0*bx*sin(phi)*(cos(phi)*(yt0-by)-sin(phi)*xt0)
      dd = bb**2 - 4.0*cc
      If (dd<0.0) Goto 10
      xx1 = (-bb+sqrt(dd))/2.0
      xx2 = (-bb-sqrt(dd))/2.0
      dd3 = abs((xx1-xt0)/cos(phi))
      dd4 = abs((xx2-xt0)/cos(phi))
      5 dd1 = min(dd1, dd3)
      dd2 = min(dd2, dd4)
      If (dd1<hipr1(13)) dd1 = 0.0
      If (dd2<hipr1(13)) dd2 = 0.0
      If (nfp(jp,10)==1 .And. ppjet>hipr1(11)) Then
        dp1 = dd1*hipr1(14)/2.0
        dp1 = min(dp1, ppjet-hipr1(11))
        pkc1 = ppjet - dp1
        dpx1 = cos(phi1)*dp1
        dpy1 = sin(phi1)*dp1
        pkc11 = pp(jp, 10) - dpx1
        pkc12 = pp(jp, 11) - dpy1
        If (dp1>0.0) Then
          cthep = pp(jp, 12)/sqrt(pp(jp,12)**2+ppjet**2)
          dpz1 = dp1*cthep/sqrt(1.0-cthep**2)
          dpe1 = sqrt(dpx1**2+dpy1**2+dpz1**2)
          eppprm = pp(jp, 4) + pp(jp, 3) - dpe1 - dpz1
          epmprm = pp(jp, 4) - pp(jp, 3) - dpe1 + dpz1
          If (eppprm<=0.0 .Or. epmprm<=0.0) Goto 15
          epp = eppprm
          epm = epmprm
          pp(jp, 10) = pkc11
          pp(jp, 11) = pkc12
          npj(jp) = npj(jp) + 1
          kfpj(jp, npj(jp)) = 21
          pjpx(jp, npj(jp)) = dpx1
          pjpy(jp, npj(jp)) = dpy1
          pjpz(jp, npj(jp)) = dpz1
          pjpe(jp, npj(jp)) = dpe1
          pjpm(jp, npj(jp)) = 0.0
          pp(jp, 3) = pp(jp, 3) - dpz1
          pp(jp, 4) = pp(jp, 4) - dpe1
        End If
      End If
      15 If (nft(jt,10)==1 .And. ptjet>hipr1(11)) Then
        dp2 = dd2*hipr1(14)/2.0
        dp2 = min(dp2, ptjet-hipr1(11))
        pkc2 = ptjet - dp2
        dpx2 = cos(phi2)*dp2
        dpy2 = sin(phi2)*dp2
        pkc21 = pt(jt, 10) - dpx2
        pkc22 = pt(jt, 11) - dpy2
        If (dp2>0.0) Then
          cthet = pt(jt, 12)/sqrt(pt(jt,12)**2+ptjet**2)
          dpz2 = dp2*cthet/sqrt(1.0-cthet**2)
          dpe2 = sqrt(dpx2**2+dpy2**2+dpz2**2)
          etpprm = pt(jt, 4) + pt(jt, 3) - dpe2 - dpz2
          etmprm = pt(jt, 4) - pt(jt, 3) - dpe2 + dpz2
          If (etpprm<=0.0 .Or. etmprm<=0.0) Goto 16
          etp = etpprm
          etm = etmprm
          pt(jt, 10) = pkc21
          pt(jt, 11) = pkc22
          ntj(jt) = ntj(jt) + 1
          kftj(jt, ntj(jt)) = 21
          pjtx(jt, ntj(jt)) = dpx2
          pjty(jt, ntj(jt)) = dpy2
          pjtz(jt, ntj(jt)) = dpz2
          pjte(jt, ntj(jt)) = dpe2
          pjtm(jt, ntj(jt)) = 0.0
          pt(jt, 3) = pt(jt, 3) - dpz2
          pt(jt, 4) = pt(jt, 4) - dpe2
        End If
      End If
      16 dpkc11 = -(pp(jp,10)-pkc11)/2.0
      dpkc12 = -(pp(jp,11)-pkc12)/2.0
      dpkc21 = -(pt(jt,10)-pkc21)/2.0
      dpkc22 = -(pt(jt,11)-pkc22)/2.0
      wp = epp + etp
      wm = epm + etm
      sw = wp*wm
    End If
  End If
  10 ptp02 = pp(jp, 1)**2 + pp(jp, 2)**2
  ptt02 = pt(jt, 1)**2 + pt(jt, 2)**2
  amq = max(pp(jp,14)+pp(jp,15), pt(jt,14)+pt(jt,15))
  amx = hipr1(1) + amq
  amp0 = amx
  dpm0 = amx
  nfdp = 0
  If (nfp(jp,5)<=2 .And. nfp(jp,3)/=0) Then
    amp0 = ulmass(nfp(jp,3))
    nfdp = nfp(jp, 3) + 2*nfp(jp, 3)/abs(nfp(jp,3))
    dpm0 = ulmass(nfdp)
    If (dpm0<=0.0) Then
      nfdp = nfdp - 2*nfdp/abs(nfdp)
      dpm0 = ulmass(nfdp)
    End If
  End If
  amt0 = amx
  dtm0 = amx
  nfdt = 0
  If (nft(jt,5)<=2 .And. nft(jt,3)/=0) Then
    amt0 = ulmass(nft(jt,3))
    nfdt = nft(jt, 3) + 2*nft(jt, 3)/abs(nft(jt,3))
    dtm0 = ulmass(nfdt)
    If (dtm0<=0.0) Then
      nfdt = nfdt - 2*nfdt/abs(nfdt)
      dtm0 = ulmass(nfdt)
    End If
  End If
  ampn = sqrt(amp0**2+ptp02)
  amtn = sqrt(amt0**2+ptt02)
  snn = (ampn+amtn)**2 + 0.001
  If (sw<snn+0.001) Goto 4000
  swptn = 4.0*(max(amp0,amt0)**2+max(ptp02,ptt02))
  swptd = 4.0*(max(dpm0,dtm0)**2+max(ptp02,ptt02))
  swptx = 4.0*(amx**2+max(ptp02,ptt02))
  If (sw<=swptn) Then
    pkcmx = 0.0
  Else If (sw>swptn .And. sw<=swptd .And. npj(jp)==0 .And. ntj(jt)==0) Then
    pkcmx = sqrt(sw/4.0-max(amp0,amt0)**2) - sqrt(max(ptp02,ptt02))
  Else If (sw>swptd .And. sw<=swptx .And. npj(jp)==0 .And. ntj(jt)==0) Then
    pkcmx = sqrt(sw/4.0-max(dpm0,dtm0)**2) - sqrt(max(ptp02,ptt02))
  Else If (sw>swptx) Then
    pkcmx = sqrt(sw/4.0-amx**2) - sqrt(max(ptp02,ptt02))
  End If
  If (nfp(jp,10)==1 .Or. nft(jt,10)==1) Then
    If (pkc1>pkcmx) Then
      pkc1 = pkcmx
      pkc11 = pkc1*cos(phi1)
      pkc12 = pkc1*sin(phi1)
      dpkc11 = -(pp(jp,10)-pkc11)/2.0
      dpkc12 = -(pp(jp,11)-pkc12)/2.0
    End If
    If (pkc2>pkcmx) Then
      pkc2 = pkcmx
      pkc21 = pkc2*cos(phi2)
      pkc22 = pkc2*sin(phi2)
      dpkc21 = -(pt(jt,10)-pkc21)/2.0
      dpkc22 = -(pt(jt,11)-pkc22)/2.0
    End If
    dpkc1 = dpkc11 + dpkc21
    dpkc2 = dpkc12 + dpkc22
    nfp(jp, 10) = -nfp(jp, 10)
    nft(jt, 10) = -nft(jt, 10)
    Goto 40
  End If
  isng = 0
  If (ihpr2(13)/=0 .And. ranart(nseed)<=hidat(4)) isng = 1
  If ((nfp(jp,5)==3 .Or. nft(jt,5)==3) .Or. (npj(jp)/=0 .Or. nfp(jp,10)/=0) .Or. (ntj(jt)/=0 .Or. nft(jt,10)/=0)) isng = 0
  If (ihpr2(5)==0) Then
    pkc = hipr1(2)*sqrt(-alog(1.0-ranart(nseed)*(1.0-exp(-pkcmx**2/hipr1(2)**2))))
    Goto 30
  End If
  xminhi = 0.0
  xmaxhi = pkcmx**2
  pkc = hirnd2(3, xminhi, xmaxhi)
  pkc = sqrt(pkc)
  If (pkc>hipr1(20)) pkc = hipr1(2)*sqrt(-alog(exp(-hipr1(20)**2/hipr1(2)**2)-ranart(nseed)*(exp(-hipr1(20)**2/hipr1(2)**2)-exp(-pkcmx**2/hipr1(2)**2))))
  If (isng==1) pkc = 0.65*sqrt(-alog(1.0-ranart(nseed)*(1.0-exp(-pkcmx**2/0.65**2))))
  30 phi0 = 2.0*hipr1(40)*ranart(nseed)
  pkc11 = pkc*sin(phi0)
  pkc12 = pkc*cos(phi0)
  pkc21 = -pkc11
  pkc22 = -pkc12
  dpkc1 = 0.0
  dpkc2 = 0.0
  40 pp11 = pp(jp, 1) + pkc11 - dpkc1
  pp12 = pp(jp, 2) + pkc12 - dpkc2
  pt11 = pt(jt, 1) + pkc21 - dpkc1
  pt12 = pt(jt, 2) + pkc22 - dpkc2
  ptp2 = pp11**2 + pp12**2
  ptt2 = pt11**2 + pt12**2
  ampn = sqrt(amp0**2+ptp2)
  amtn = sqrt(amt0**2+ptt2)
  snn = (ampn+amtn)**2 + 0.001
  wp = epp + etp
  wm = epm + etm
  sw = wp*wm
  If (sw<snn) Then
    miss = miss + 1
    If (miss<=100) Then
      pkc = 0.0
      Goto 30
    End If
    If (ihpr2(10)/=0) Write (6, *) 'Error occured in Pt kick section of HIJSFT'
    Goto 4000
  End If
  ampd = sqrt(dpm0**2+ptp2)
  amtd = sqrt(dtm0**2+ptt2)
  ampx = sqrt(amx**2+ptp2)
  amtx = sqrt(amx**2+ptt2)
  dpn = ampn**2/sw
  dtn = amtn**2/sw
  dpd = ampd**2/sw
  dtd = amtd**2/sw
  dpx = ampx**2/sw
  dtx = amtx**2/sw
  spntd = (ampn+amtd)**2
  spntx = (ampn+amtx)**2
  spdtn = (ampd+amtn)**2
  spxtn = (ampx+amtn)**2
  spdtx = (ampd+amtx)**2
  spxtd = (ampx+amtd)**2
  sdd = (ampd+amtd)**2
  sxx = (ampx+amtx)**2
  If (sw>sxx+0.001) Then
    If (isng==0) Then
      d1 = dpx
      d2 = dtx
      nfp3 = 0
      nft3 = 0
      Goto 400
    Else
      If ((nfp(jp,5)==3 .And. nft(jt,5)==3) .Or. (npj(jp)/=0 .Or. nfp(jp,10)/=0) .Or. (ntj(jt)/=0 .Or. nft(jt,10)/=0)) Then
        d1 = dpx
        d2 = dtx
        nfp3 = 0
        nft3 = 0
        Goto 400
      End If
      If (ranart(nseed)>0.5 .Or. (nft(jt,5)>2 .Or. ntj(jt)/=0 .Or. nft(jt,10)/=0)) Then
        d1 = dpn
        d2 = dtx
        nfp3 = nfp(jp, 3)
        nft3 = 0
        Goto 220
      Else
        d1 = dpx
        d2 = dtn
        nfp3 = 0
        nft3 = nft(jt, 3)
        Goto 240
      End If
    End If
  Else If (sw>max(spdtx,spxtd)+0.001 .And. sw<=sxx+0.001) Then
    If (((npj(jp)==0 .And. ntj(jt)==0 .And. ranart(nseed)>0.5) .Or. (npj(jp)==0 .And. ntj(jt)/=0)) .And. nfp(jp,5)<=2) Then
      d1 = dpd
      d2 = dtx
      nfp3 = nfdp
      nft3 = 0
      Goto 220
    Else If (ntj(jt)==0 .And. nft(jt,5)<=2) Then
      d1 = dpx
      d2 = dtd
      nfp3 = 0
      nft3 = nfdt
      Goto 240
    End If
    Goto 4000
  Else If (sw>min(spdtx,spxtd)+0.001 .And. sw<=max(spdtx,spxtd)+0.001) Then
    If (spdtx<=spxtd .And. npj(jp)==0 .And. nfp(jp,5)<=2) Then
      d1 = dpd
      d2 = dtx
      nfp3 = nfdp
      nft3 = 0
      Goto 220
    Else If (spdtx>spxtd .And. ntj(jt)==0 .And. nft(jt,5)<=2) Then
      d1 = dpx
      d2 = dtd
      nfp3 = 0
      nft3 = nfdt
      Goto 240
    End If
    If (((npj(jp)==0 .And. ntj(jt)==0 .And. ranart(nseed)>0.5) .Or. (npj(jp)==0 .And. ntj(jt)/=0)) .And. nfp(jp,5)<=2) Then
      d1 = dpn
      d2 = dtx
      nfp3 = nfp(jp, 3)
      nft3 = 0
      Goto 220
    Else If (ntj(jt)==0 .And. nft(jt,5)<=2) Then
      d1 = dpx
      d2 = dtn
      nfp3 = 0
      nft3 = nft(jt, 3)
      Goto 240
    End If
    Goto 4000
  Else If (sw>max(spntx,spxtn)+0.001 .And. sw<=min(spdtx,spxtd)+0.001) Then
    If (((npj(jp)==0 .And. ntj(jt)==0 .And. ranart(nseed)>0.5) .Or. (npj(jp)==0 .And. ntj(jt)/=0)) .And. nfp(jp,5)<=2) Then
      d1 = dpn
      d2 = dtx
      nfp3 = nfp(jp, 3)
      nft3 = 0
      Goto 220
    Else If (ntj(jt)==0 .And. nft(jt,5)<=2) Then
      d1 = dpx
      d2 = dtn
      nfp3 = 0
      nft3 = nft(jt, 3)
      Goto 240
    End If
    Goto 4000
  Else If (sw>min(spntx,spxtn)+0.001 .And. sw<=max(spntx,spxtn)+0.001) Then
    If (spntx<=spxtn .And. npj(jp)==0 .And. nfp(jp,5)<=2) Then
      d1 = dpn
      d2 = dtx
      nfp3 = nfp(jp, 3)
      nft3 = 0
      Goto 220
    Else If (spntx>spxtn .And. ntj(jt)==0 .And. nft(jt,5)<=2) Then
      d1 = dpx
      d2 = dtn
      nfp3 = 0
      nft3 = nft(jt, 3)
      Goto 240
    End If
    Goto 4000
  Else If (sw<=min(spntx,spxtn)+0.001 .And. (npj(jp)/=0 .Or. ntj(jt)/=0)) Then
    Goto 4000
  Else If (sw<=min(spntx,spxtn)+0.001 .And. nfp(jp,5)>2 .And. nft(jt,5)>2) Then
    Goto 4000
  Else If (sw>sdd+0.001 .And. sw<=min(spntx,spxtn)+0.001) Then
    d1 = dpd
    d2 = dtd
    nfp3 = nfdp
    nft3 = nfdt
    Goto 100
  Else If (sw>max(spntd,spdtn)+0.001 .And. sw<=sdd+0.001) Then
    If (ranart(nseed)>0.5) Then
      d1 = dpd
      d2 = dtn
      nfp3 = nfdp
      nft3 = nft(jt, 3)
      Goto 100
    Else
      d1 = dpn
      d2 = dtd
      nfp3 = nfp(jp, 3)
      nft3 = nfdt
      Goto 100
    End If
  Else If (sw>min(spntd,spdtn)+0.001 .And. sw<=max(spntd,spdtn)+0.001) Then
    If (spntd>spdtn) Then
      d1 = dpd
      d2 = dtn
      nfp3 = nfdp
      nft3 = nft(jt, 3)
      Goto 100
    Else
      d1 = dpn
      d2 = dtd
      nfp3 = nfp(jp, 3)
      nft3 = nfdt
      Goto 100
    End If
  Else If (sw<=min(spntd,spdtn)+0.001) Then
    d1 = dpn
    d2 = dtn
    nfp3 = nfp(jp, 3)
    nft3 = nft(jt, 3)
    Goto 100
  End If
  Write (6, *) ' Error in HIJSFT: There is no path to here'
  Return
  100 nfp5 = max(2, nfp(jp,5))
  nft5 = max(2, nft(jt,5))
  bb1 = 1.0 + d1 - d2
  bb2 = 1.0 + d2 - d1
  If (bb1**2<4.0*d1 .Or. bb2**2<4.0*d2) Then
    miss = miss + 1
    If (miss>100 .Or. pkc==0.0) Goto 3000
    pkc = pkc*0.5
    Goto 30
  End If
  If (ranart(nseed)<0.5) Then
    x1 = (bb1-sqrt(bb1**2-4.0*d1))/2.0
    x2 = (bb2-sqrt(bb2**2-4.0*d2))/2.0
  Else
    x1 = (bb1+sqrt(bb1**2-4.0*d1))/2.0
    x2 = (bb2+sqrt(bb2**2-4.0*d2))/2.0
  End If
  ihnt2(13) = 2
  Goto 600
  220 nfp5 = max(2, nfp(jp,5))
  nft5 = 3
  If (nfp3==0) nfp5 = 3
  bb2 = 1.0 + d2 - d1
  If (bb2**2<4.0*d2) Then
    miss = miss + 1
    If (miss>100 .Or. pkc==0.0) Goto 3000
    pkc = pkc*0.5
    Goto 30
  End If
  xmin = (bb2-sqrt(bb2**2-4.0*d2))/2.0
  xmax = (bb2+sqrt(bb2**2-4.0*d2))/2.0
  miss4 = 0
  222 x2 = hirnd2(6, xmin, xmax)
  x1 = d1/(1.0-x2)
  If (x2*(1.0-x1)<(d2+1.E-4/sw)) Then
    miss4 = miss4 + 1
    If (miss4<=1000) Goto 222
    Goto 5000
  End If
  ihnt2(13) = 2
  Goto 600
  240 nfp5 = 3
  nft5 = max(2, nft(jt,5))
  If (nft3==0) nft5 = 3
  bb1 = 1.0 + d1 - d2
  If (bb1**2<4.0*d1) Then
    miss = miss + 1
    If (miss>100 .Or. pkc==0.0) Goto 3000
    pkc = pkc*0.5
    Goto 30
  End If
  xmin = (bb1-sqrt(bb1**2-4.0*d1))/2.0
  xmax = (bb1+sqrt(bb1**2-4.0*d1))/2.0
  miss4 = 0
  242 x1 = hirnd2(6, xmin, xmax)
  x2 = d2/(1.0-x1)
  If (x1*(1.0-x2)<(d1+1.E-4/sw)) Then
    miss4 = miss4 + 1
    If (miss4<=1000) Goto 242
    Goto 5000
  End If
  ihnt2(13) = 2
  Goto 600
  400 nfp5 = 3
  nft5 = 3
  bb1 = 1.0 + d1 - d2
  bb2 = 1.0 + d2 - d1
  If (bb1**2<4.0*d1 .Or. bb2**2<4.0*d2) Then
    miss = miss + 1
    If (miss>100 .Or. pkc==0.0) Goto 3000
    pkc = pkc*0.5
    Goto 30
  End If
  xmin1 = (bb1-sqrt(bb1**2-4.0*d1))/2.0
  xmax1 = (bb1+sqrt(bb1**2-4.0*d1))/2.0
  xmin2 = (bb2-sqrt(bb2**2-4.0*d2))/2.0
  xmax2 = (bb2+sqrt(bb2**2-4.0*d2))/2.0
  miss4 = 0
  410 x1 = hirnd2(4, xmin1, xmax1)
  x2 = hirnd2(4, xmin2, xmax2)
  If (nfp(jp,5)==3 .Or. nft(jt,5)==3) Then
    x1 = hirnd2(6, xmin1, xmax1)
    x2 = hirnd2(6, xmin2, xmax2)
  End If
  If (abs(nfp(jp,1)*nfp(jp,2))>1000000 .Or. abs(nfp(jp,1)*nfp(jp,2))<100) Then
    x1 = hirnd2(5, xmin1, xmax1)
  End If
  If (abs(nft(jt,1)*nft(jt,2))>1000000 .Or. abs(nft(jt,1)*nft(jt,2))<100) Then
    x2 = hirnd2(5, xmin2, xmax2)
  End If
  If (abs(nfp(jp,1)*nfp(jp,2))>1000000) x1 = 1.0 - x1
  xxp = x1*(1.0-x2)
  xxt = x2*(1.0-x1)
  If (xxp<(d1+1.E-4/sw) .Or. xxt<(d2+1.E-4/sw)) Then
    miss4 = miss4 + 1
    If (miss4<=1000) Goto 410
    Goto 5000
  End If
  ihnt2(13) = 3
  600 Continue
  If (x1*(1.0-x2)<(ampn**2-1.E-4)/sw .Or. x2*(1.0-x1)<(amtn**2-1.E-4)/sw) Then
    miss = miss + 1
    If (miss>100 .Or. pkc==0.0) Goto 2000
    pkc = 0.0
    Goto 30
  End If
  epp = (1.0-x2)*wp
  epm = x1*wm
  etp = x2*wp
  etm = (1.0-x1)*wm
  pp(jp, 3) = (epp-epm)/2.0
  pp(jp, 4) = (epp+epm)/2.0
  If (epp*epm-ptp2<0.0) Goto 6000
  pp(jp, 5) = sqrt(epp*epm-ptp2)
  nfp(jp, 3) = nfp3
  nfp(jp, 5) = nfp5
  pt(jt, 3) = (etp-etm)/2.0
  pt(jt, 4) = (etp+etm)/2.0
  If (etp*etm-ptt2<0.0) Goto 6000
  pt(jt, 5) = sqrt(etp*etm-ptt2)
  nft(jt, 3) = nft3
  nft(jt, 5) = nft5
  pp(jp, 1) = pp11 - pkc11
  pp(jp, 2) = pp12 - pkc12
  kcdip = 1
  kcdit = 1
  If (abs(nfp(jp,1)*nfp(jp,2))>1000000 .Or. abs(nfp(jp,1)*nfp(jp,2))<100) Then
    kcdip = 0
  End If
  If (abs(nft(jt,1)*nft(jt,2))>1000000 .Or. abs(nft(jt,1)*nft(jt,2))<100) Then
    kcdit = 0
  End If
  If ((kcdip==0 .And. ranart(nseed)<0.5) .Or. (kcdip/=0 .And. ranart(nseed)<0.5/(1.0+(pkc11**2+pkc12**2)/hipr1(22)**2))) Then
    pp(jp, 6) = (pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0 + pp(jp, 6)
    pp(jp, 7) = (pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0 + pp(jp, 7)
    pp(jp, 8) = (pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0 + pp(jp, 8) + pkc11
    pp(jp, 9) = (pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0 + pp(jp, 9) + pkc12
  Else
    pp(jp, 8) = (pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0 + pp(jp, 8)
    pp(jp, 9) = (pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0 + pp(jp, 9)
    pp(jp, 6) = (pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0 + pp(jp, 6) + pkc11
    pp(jp, 7) = (pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0 + pp(jp, 7) + pkc12
  End If
  pp(jp, 1) = pp(jp, 6) + pp(jp, 8)
  pp(jp, 2) = pp(jp, 7) + pp(jp, 9)
  pt(jt, 1) = pt11 - pkc21
  pt(jt, 2) = pt12 - pkc22
  If ((kcdit==0 .And. ranart(nseed)<0.5) .Or. (kcdit/=0 .And. ranart(nseed)<0.5/(1.0+(pkc21**2+pkc22**2)/hipr1(22)**2))) Then
    pt(jt, 6) = (pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0 + pt(jt, 6)
    pt(jt, 7) = (pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0 + pt(jt, 7)
    pt(jt, 8) = (pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0 + pt(jt, 8) + pkc21
    pt(jt, 9) = (pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0 + pt(jt, 9) + pkc22
  Else
    pt(jt, 8) = (pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0 + pt(jt, 8)
    pt(jt, 9) = (pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0 + pt(jt, 9)
    pt(jt, 6) = (pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0 + pt(jt, 6) + pkc21
    pt(jt, 7) = (pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0 + pt(jt, 7) + pkc22
  End If
  pt(jt, 1) = pt(jt, 6) + pt(jt, 8)
  pt(jt, 2) = pt(jt, 7) + pt(jt, 9)
  If (npj(jp)/=0) nfp(jp, 5) = 3
  If (ntj(jt)/=0) nft(jt, 5) = 3
  If (epp/(epm+0.0001)<etp/(etm+0.0001) .And. abs(nfp(jp,1)*nfp(jp,2))<1000000) Then
    Do jsb = 1, 15
      psb = pp(jp, jsb)
      pp(jp, jsb) = pt(jt, jsb)
      pt(jt, jsb) = psb
      nsb = nfp(jp, jsb)
      nfp(jp, jsb) = nft(jt, jsb)
      nft(jt, jsb) = nsb
    End Do
  End If
  Return
  1000 ierror = 1
  If (ihpr2(10)==0) Return
  Write (6, *) '     Fatal HIJSFT start error,abandon this event'
  Write (6, *) '     PROJ E+,E-,W+', epp, epm, wp
  Write (6, *) '     TARG E+,E-,W-', etp, etm, wm
  Write (6, *) '     W+*W-, (APN+ATN)^2', sw, snn
  Return
  2000 ierror = 0
  If (ihpr2(10)==0) Return
  Write (6, *) '     (2)energy partition fail,'
  Write (6, *) '     HIJSFT not performed, but continue'
  Write (6, *) '     MP1,MPN', x1*(1.0-x2)*sw, ampn**2
  Write (6, *) '     MT2,MTN', x2*(1.0-x1)*sw, amtn**2
  Return
  3000 ierror = 0
  If (ihpr2(10)==0) Return
  Write (6, *) '     (3)something is wrong with the pt kick, '
  Write (6, *) '     HIJSFT not performed, but continue'
  Write (6, *) '     D1=', d1, ' D2=', d2, ' SW=', sw
  Write (6, *) '     HISTORY NFP5=', nfp(jp, 5), ' NFT5=', nft(jt, 5)
  Write (6, *) '     THIS COLLISON NFP5=', nfp5, ' NFT5=', nft5
  Write (6, *) '     # OF JET IN PROJ', npj(jp), ' IN TARG', ntj(jt)
  Return
  4000 ierror = 0
  If (ihpr2(10)==0) Return
  Write (6, *) '     (4)unable to choose process, but not harmful'
  Write (6, *) '     HIJSFT not performed, but continue'
  Write (6, *) '     PTP=', sqrt(ptp2), ' PTT=', sqrt(ptt2), ' SW=', sw
  Write (6, *) '     AMCUT=', amx, ' JP=', jp, ' JT=', jt
  Write (6, *) '     HISTORY NFP5=', nfp(jp, 5), ' NFT5=', nft(jt, 5)
  Return
  5000 ierror = 0
  If (ihpr2(10)==0) Return
  Write (6, *) '     energy partition failed(5),for limited try'
  Write (6, *) '     HIJSFT not performed, but continue'
  Write (6, *) '     NFP5=', nfp5, ' NFT5=', nft5
  Write (6, *) '     D1', d1, ' X1(1-X2)', x1*(1.0-x2)
  Write (6, *) '     D2', d2, ' X2(1-X1)', x2*(1.0-x1)
  Return
  6000 pkc = 0.0
  miss = miss + 1
  If (miss<100) Goto 30
  ierror = 1
  If (ihpr2(10)==0) Return
  Write (6, *) ' ERROR OCCURED, HIJSFT NOT PERFORMED'
  Write (6, *) ' Abort this event'
  Write (6, *) 'MTP,PTP2', epp*epm, ptp2, '  MTT,PTT2', etp*etm, ptt2
  Return
End Subroutine hijsft

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
!
! This file is released under terms of BSD license
! See LICENSE file for more information
!
! Test case for issue #598
module sca49
  contains
  subroutine sub1()
    integer :: i, j
    integer, parameter :: ny = 10, nx = 10

    !$claw define dimension k(1:n_term) &
    !$claw sca
    do j=1, ny
      do i=1, nx
      end do
    end do
  end subroutine sub1
end module

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text:
subroutine dff_test2
use enforce
use cluster
use paracluster

implicit none

integer :: i,n,k
real, allocatable :: dft(:)
real :: arg,freq

allocate(dft(npas))

do i = 1,npas
	dft(i) = 0.d0
enddo

do k=1,npas
	freq = 1.d0/(k*tstep)
	arg = -(2.d0*pi*dble(1.d0/freq)/npas)
	do n =1,npas
		dft(k) = dft(k) + cos(arg*dble(n))*vel_act_est(n)
	enddo
	write(60,*) freq,dft(k)*dft(k)
enddo

end subroutine dff_test2

!======================================================

subroutine dff_test
use ENFORCE!, only: vel_act_est
use CLUSTER
use PARACLUSTER

implicit none

!variables needed:      - N = npas + 1
!			- k = index for N
!			- n = index for chosen frequency
!PLAN:
!	INITIALISE -- size of dft

integer :: N,N_i,freq,a,b,i
real, allocatable :: dft(:,:)
real :: arg

!N = npas + 1!SIZE(vel_acf_est)
!do i = 0,npas
!	write(*,*) vel_act_est(i)
!enddo

allocate(dft(0:npas,2))

do a = 0,npas
	dft(a,1) = 0.d0
	dft(a,2) = 0.d0
end do

!open(11,file = 'fort.59', status = 'old' , action = 'read')

!do b = 0,N
!	read(11,*) dft(b)
!end do
!open(60, file = 'fort.60', status = 'new', action = 'write')

do freq = 1,npas
	do N_i = 1,npas
	arg = -(2*pi*dble(N_i)/npas)
	!write(13,*) N_i,arg
	dft(freq,1) = dft(freq,1) + cos(arg*dble(freq)*tstep)*vel_act_est(N_i)!/npas
	!write(14,*) N_i ,dft(N_i,1) , cos(arg) ,vel_act_est(freq)
	!dft(N_i,2) = dft(N_i,2) + sin(arg*dble(freq))*vel_act_est(freq)!/npas
	!write(15,*) N_i,dft(N_i,2) , sin(arg) ,vel_act_est(freq)
	end do
write(60,*) freq , dft(freq,1)*dft(freq,1) !+ dft(N_i,2)*dft(N_i,2)
end do

!do N_i = 0,npas
	!dft(N_i,1)=dft(N_i,1)/npas
	!dft(N_i,2)=dft(N_i,2)/npas
	!write(60,*) N_i , dft(N_i,1)*dft(N_i,1) + dft(N_i,2)*dft(N_i,2)
!end do

end subroutine dff_test

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS


text:
!
! LBLRTM_r2p1p1_Module
!
! Module containing procedures for LBLRTM Record r2p1p1.
!
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, 24-Dec-2012
!                     paul.vandelst@noaa.gov
!

MODULE LBLRTM_r2p1p1_Module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE File_Utility   , ONLY: File_Open
  USE Message_Handler, ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  ! Line-by-line model parameters
  USE LBL_Parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Visibility
  ! ----------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: LBLRTM_r2p1p1_type
  ! Procedures
  PUBLIC :: LBLRTM_r2p1p1_Write

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Message string length
  INTEGER, PARAMETER :: ML = 256
  ! The record I/O format
  CHARACTER(*), PARAMETER :: LBLRTM_R2P1P1_FMT = '(f10.4,f10.4,f10.4,a3,i2,1x,f7.2,f8.3,f7.2,f7.2,f8.3,f7.2)'

  ! -------------
  ! Derived types
  ! -------------
  TYPE :: LBLRTM_r2p1p1_type
    REAL(fp)     :: pave   = 0.0_fp  !  Average pressure of layer
    REAL(fp)     :: tave   = 0.0_fp  !  Average temperature of layer
    REAL(fp)     :: secntk = 0.0_fp  !  Column amount scale factor for layer
    CHARACTER(3) :: ityl   = ''      !  Overrides control of DV ratio between layers
    INTEGER      :: ipath  = 0       !  Merge direction
    REAL(fp)     :: altzb  = 0.0_fp  !  Altitude of BOTTOM of current layer
    REAL(fp)     :: pzb    = 0.0_fp  !  Pressure at BOTTOM of current layer
    REAL(fp)     :: tzb    = 0.0_fp  !  Temperature at BOTTOM of current layer
    REAL(fp)     :: altz   = 0.0_fp  !  Altitude of TOP of current layer
    REAL(fp)     :: pz     = 0.0_fp  !  Pressure at TOP of current layer
    REAL(fp)     :: tz     = 0.0_fp  !  Temperature at TOP of current layer
  END TYPE LBLRTM_r2p1p1_type

CONTAINS

  FUNCTION LBLRTM_r2p1p1_Write(r2p1p1,fid) RESULT(err_stat)

    ! Arguments
    TYPE(LBLRTM_r2p1p1_type), INTENT(IN) :: r2p1p1
    INTEGER                 , INTENT(IN) :: fid
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_r2p1p1_Write'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Setup
    err_stat = SUCCESS
    ! ...Check unit is open
    IF ( .NOT. File_Open(fid) ) THEN
      msg = 'File unit is not connected'
      CALL Cleanup(); RETURN
    END IF

    ! Write the record
    WRITE( fid,FMT=LBLRTM_R2P1P1_FMT,IOSTAT=io_stat,IOMSG=io_msg) r2p1p1
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing record - '//TRIM(io_msg)
      CALL Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE CleanUp()
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE CleanUp

  END FUNCTION LBLRTM_r2p1p1_Write

END MODULE LBLRTM_r2p1p1_Module

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
      subroutine chabund
c*****************************************************************************
c  this routine changes the input abundances according to the users request
c  and the re-runs the synthesis with the new abundances OVERWRITING the
c  original output files, but using the same model, linelist, etc
c*****************************************************************************

      implicit real*8 (a-h,o-z)
      include 'Atmos.com'
      include 'Factor.com'
      include 'Mol.com'
      include 'Pstuff.com'
      include 'Multistar.com'
      real*8 xnum
      real*4 shortnum
      character choice2

c*****set the local temporary parameters to what they were for the last synth
      newnumiso = numiso
      newnumisosyn = numisosyn
      newnumpecatom = numpecatom
      newnumatomsyn = numatomsyn
      if (newnumatomsyn .eq. 0) newnumatomsyn = 1
      do i=3,95
         newpec(i) = pec(i)
         do j=1,newnumatomsyn
            newpecabund(i,j) = pecabund(i,j)
         enddo
      enddo
      do i=1,newnumiso
         if(newnumiso .ge. 1) then
            newisotope(i)=isotope(i)
            do j=1,newnumisosyn
               newisoabund(i,j)=isoabund(i,j)
            enddo
         endif
      enddo
      if (newnumatomsyn .le. 0) newnumatomsyn = 1

c*****SPECIAL CASE
c     for the "binary" driver, present a limited set of options; only
c     the abundances of the elements already chosen for variation in
c     the parameter file can be altered
11    if (control .eq. 'binary ') then
         istat = ivcleof(4,1)
         if (syncount .eq. 1) then
            array = 'Abundance Alterations for the PRIMARY STAR:   '
         else
            array = 'Abundance Alterations for the SECONDARY STAR: '
         endif
         istat = ivwrite(5,5,array,46)
         array = 'Atom   Abundance (relative to initial model)'
         istat = ivwrite(6,1,array,49)
         line = 6
         do i=3,95
            if(newpec(i) .eq. 1) then
               line = line + 1
               write (array,1001) i,(newpecabund(i,j),j=1,newnumatomsyn)
               istat = ivwrite(line,1,array,57)
            endif
         enddo
         array = 'Options: c = change abundances '
         istat = ivwrite(line+1,1,array,41)
         array = '         x = exit no changes    q = rerun synth'
         istat = ivwrite(line+2,1,array,47)
         array = 'What is your choice? '
         nchars = 21
         call getasci (nchars,line+4)
         choice = chinfo(1:1)
         if (choice.ne.'c' .and. choice.ne.'x' .and.
     .       choice.ne.'q') go to 11
      endif

c*****for the "synth" driver, present the user with the current abundance
c     alterations, and the options; find out what is desired
1     if (control .eq. 'synth'  .or.
     .    control .eq. 'isotop' .or.
     .    control .eq. 'synpop') then
         line = 4
         istat = ivcleof(line,1)
         write (array,1006)
         line = line + 1
         istat = ivwrite(6,1,array,72)
         line = line + 1
         do i=3,95
            if(newpec(i) .eq. 1) then
               line = line + 1
               write (array,1001) i,(newpecabund(i,j),j=1,newnumatomsyn)
               istat = ivwrite(line,1,array,57)
            endif
         enddo
         if (newnumiso .gt. 0) then
            do i=1,newnumiso
               line = line + 1
               write (array,1002) i, newisotope(i),
     .                            (newisoabund(i,j),j=1,newnumisosyn)
               istat = ivwrite(line,1,array,68)
            enddo
         endif
         write (array,1007)
         istat = ivwrite(line+1,1,array,61)
         write (array,1008)
         istat = ivwrite(line+2,1,array,72)
         array = 'What is your choice? '
         nchars = 21
         call getasci (nchars,line+3)
         choice = chinfo(1:1)
         if (choice.ne.'c' .and. choice.ne.'i' .and.
     .       choice.ne.'n' .and. choice.ne.'x' .and.
     .       choice.ne.'q') go to 1
      endif

c*****for "synth", change elemental abundances
      if     (choice .eq. 'c') then
20       istat = ivcleof(line+1,1)
         array = 'Which element to change? '
         nchars = 25
         call getnum (nchars,line+1,xnum,shortnum)
         if (xnum.lt.2.0 .or. xnum.gt.95.) go to 20
         j = nint(xnum)
         if (newpec(j) .eq. 1) then
25          array = 'n = new abundances, or z = zero offsets? '
            nchars = 41
            call getasci (nchars,line+2)
            choice2 = chinfo(1:1)
            if     (choice2 .eq. 'z') then
               newpec(j) = 0
               do i=1,5
                  newpecabund(j,i) = 0.0
               enddo
               newnumpecatom = newnumpecatom - 1
            elseif (choice2 .eq. 'n') then
               write (array,1004)
               istat = ivwrite(line+3,1,array,40)
               read (*,*) (newpecabund(j,i),i=1,newnumatomsyn)
            else
               go to 25
            endif
         else
            newpec(j) = 1
            write (array,1004)
            istat = ivwrite(line+2,1,array,40)
            read (*,*) (newpecabund(j,i),i=1,newnumatomsyn)
            newnumpecatom = newnumpecatom + 1
         endif

c*****for "synth", change isotopic factors
      elseif (choice .eq. 'i') then
         istat = ivcleof(line+1,1)
         array = 'Options: c = change an isotopic factor'
         istat = ivwrite(line+1,1,array,49)
         array = '         n = enter a new isotope'
         istat = ivwrite(line+2,1,array,48)
30       array = 'What is your choice? '
         nchars = 22
         call getasci (nchars,line+3)
         choice2 = chinfo(1:1)
         if (choice2 .eq. 'c') then
35          istat = ivcleof(line+1,1)
            array =  'Which isotope number from the list? '
            nchars = 36
            call getnum (nchars,line+1,xnum,shortnum)
            j = nint(xnum)
            if (j.lt.1 .or. j.gt.newnumiso) go to 35
            istat = ivcleof(line+2,1)
            array = 'What are the new division factors? '
            istat = ivwrite(line+2,1,array,35)
            read (*,*) (newisoabund(j,i),i=1,newnumisosyn)
         elseif (choice2 .eq. 'n') then
            newnumiso = newnumiso + 1
            istat = ivcleof(line+1,1)
            array = 'What is the new isotope designation? '
            nchars = 37
            call getnum (nchars,line+4,xnum,shortnum)
            newisotope(newnumiso) = xnum
            array = 'What are its division factors? '
            istat = ivwrite(line+2,1,array,31)
            read (*,*) (newisoabund(newnumiso,i),i=1,newnumisosyn)
         else
            go to 30
         endif

c*****for "synth", change the number of syntheses
      elseif (choice .eq. 'n') then
55       array = 'How many synths? '
         nchars = 17
         call getnum (nchars,line+5,xnum,shortnum)
         if (xnum .gt. 5.) go to 55
         newnumatomsyn = nint(xnum)
         go to 1

c*****for "synth", exit the routine without changing anything
      elseif (choice .eq. 'x') then
         return

c*****for "synth", make the proposed alterations permanent; then
c     return to the calling routine, which allegedly will
c     redo the syntheses.
      elseif (choice .eq. 'q') then
         numiso = newnumiso
         numisosyn = newnumisosyn
         numpecatom = newnumpecatom
         numatomsyn = newnumatomsyn
         do i=3,95
            pec(i) = newpec(i)
            do j=1,numatomsyn
               pecabund(i,j) = newpecabund(i,j)
            enddo
         enddo
         do i=1,numiso
            if(numiso .ge. 1) then
               isotope(i) = newisotope(i)
               do j=1,numisosyn
                  isoabund(i,j)=newisoabund(i,j)
               enddo
            endif
         enddo
         return
      endif

c*****loop back and print out the main menu again
      if (control .eq. 'synth  ') then
         go to 1
      else
         go to 11
      endif

c*****format statements
1001  format (i3, 4x, 5(f8.3, 2x))
1002  format (i2, 2x, f10.5, 5(3x, f8.3))
1004  format ('Enter the new offsets on the line below:')
1006     format ('element, abundance offsets   OR   ',
     .            'isotope number, isotope name, factors')
1007     format ('Options: c = change abundance       ',
     .           'i = change isotopic ratio')
1008     format('         n = change # syntheses     q = rerun ',
     .          'syntheses        x = exit')

      end

ext: f


!*==zpbcon.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b ZPBCON
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZPBCON + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
!                          RWORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          UPLO
!       INTEGER            INFO, KD, LDAB, N
!       DOUBLE PRECISION   ANORM, RCOND
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   RWORK( * )
!       COMPLEX*16         AB( LDAB, * ), WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZPBCON estimates the reciprocal of the condition number (in the
!> 1-norm) of a complex Hermitian positive definite band matrix using
!> the Cholesky factorization A = U**H*U or A = L*L**H computed by
!> ZPBTRF.
!>
!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
!> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] UPLO
!> \verbatim
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangular factor stored in AB;
!>          = 'L':  Lower triangular factor stored in AB.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!> \endverbatim
!>
!> \param[in] KD
!> \verbatim
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
!> \endverbatim
!>
!> \param[in] AB
!> \verbatim
!>          AB is COMPLEX*16 array, dimension (LDAB,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**H*U or A = L*L**H of the band matrix A, stored in the
!>          first KD+1 rows of the array.  The j-th column of U or L is
!>          stored in the j-th column of the array AB as follows:
!>          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).
!> \endverbatim
!>
!> \param[in] LDAB
!> \verbatim
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD+1.
!> \endverbatim
!>
!> \param[in] ANORM
!> \verbatim
!>          ANORM is DOUBLE PRECISION
!>          The 1-norm (or infinity-norm) of the Hermitian band matrix A.
!> \endverbatim
!>
!> \param[out] RCOND
!> \verbatim
!>          RCOND is DOUBLE PRECISION
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
!>          estimate of the 1-norm of inv(A) computed in this routine.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is COMPLEX*16 array, dimension (2*N)
!> \endverbatim
!>
!> \param[out] RWORK
!> \verbatim
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!> \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 Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup complex16OTHERcomputational
!
!  =====================================================================
      SUBROUTINE ZPBCON(Uplo,N,Kd,Ab,Ldab,Anorm,Rcond,Work,Rwork,Info)
      IMPLICIT NONE
!*--ZPBCON136
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Uplo
      INTEGER Info , Kd , Ldab , N
      DOUBLE PRECISION Anorm , Rcond
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION Rwork(*)
      COMPLEX*16 Ab(Ldab,*) , Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ONE , ZERO
      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
!     ..
!     .. Local Scalars ..
      LOGICAL upper
      CHARACTER normin
      INTEGER ix , kase
      DOUBLE PRECISION ainvnm , scale , scalel , scaleu , smlnum
      COMPLEX*16 zdum
!     ..
!     .. Local Arrays ..
      INTEGER isave(3)
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      INTEGER IZAMAX
      DOUBLE PRECISION DLAMCH
      EXTERNAL LSAME , IZAMAX , DLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL XERBLA , ZDRSCL , ZLACN2 , ZLATBS
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , DBLE , DIMAG
!     ..
!     .. Statement Functions ..
      DOUBLE PRECISION CABS1
!     ..
!     .. Statement Function definitions ..
      CABS1(zdum) = ABS(DBLE(zdum)) + ABS(DIMAG(zdum))
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      Info = 0
      upper = LSAME(Uplo,'U')
      IF ( .NOT.upper .AND. .NOT.LSAME(Uplo,'L') ) THEN
         Info = -1
      ELSEIF ( N<0 ) THEN
         Info = -2
      ELSEIF ( Kd<0 ) THEN
         Info = -3
      ELSEIF ( Ldab<Kd+1 ) THEN
         Info = -5
      ELSEIF ( Anorm<ZERO ) THEN
         Info = -6
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('ZPBCON',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      Rcond = ZERO
      IF ( N==0 ) THEN
         Rcond = ONE
         RETURN
      ELSEIF ( Anorm==ZERO ) THEN
         RETURN
      ENDIF
!
      smlnum = DLAMCH('Safe minimum')
!
!     Estimate the 1-norm of the inverse.
!
      kase = 0
      normin = 'N'
      DO
         CALL ZLACN2(N,Work(N+1),Work,ainvnm,kase,isave)
         IF ( kase/=0 ) THEN
            IF ( upper ) THEN
!
!           Multiply by inv(U**H).
!
               CALL ZLATBS('Upper','Conjugate transpose','Non-unit',    &
     &                     normin,N,Kd,Ab,Ldab,Work,scalel,Rwork,Info)
               normin = 'Y'
!
!           Multiply by inv(U).
!
               CALL ZLATBS('Upper','No transpose','Non-unit',normin,N,  &
     &                     Kd,Ab,Ldab,Work,scaleu,Rwork,Info)
            ELSE
!
!           Multiply by inv(L).
!
               CALL ZLATBS('Lower','No transpose','Non-unit',normin,N,  &
     &                     Kd,Ab,Ldab,Work,scalel,Rwork,Info)
               normin = 'Y'
!
!           Multiply by inv(L**H).
!
               CALL ZLATBS('Lower','Conjugate transpose','Non-unit',    &
     &                     normin,N,Kd,Ab,Ldab,Work,scaleu,Rwork,Info)
            ENDIF
!
!        Multiply by 1/SCALE if doing so will not cause overflow.
!
            scale = scalel*scaleu
            IF ( scale/=ONE ) THEN
               ix = IZAMAX(N,Work,1)
               IF ( scale<CABS1(Work(ix))*smlnum .OR. scale==ZERO ) EXIT
               CALL ZDRSCL(N,scale,Work,1)
            ENDIF
            CYCLE
         ENDIF
!
!     Compute the estimate of the reciprocal condition number.
!
         IF ( ainvnm/=ZERO ) Rcond = (ONE/ainvnm)/Anorm
         EXIT
      ENDDO
!
!
!
!     End of ZPBCON
!
      END SUBROUTINE ZPBCON


text:
c      @(#)rw_main.f	1.2 11/23/94
c
c      file: rw_main.f
       subroutine rw_main
       return
       end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
MODULE params_obs

USE common, ONLY: r_size

IMPLICIT NONE

PUBLIC

INTEGER,SAVE :: nobs
INTEGER,PARAMETER :: nid_obs=8
INTEGER,PARAMETER :: id_u_obs=2819
INTEGER,PARAMETER :: id_v_obs=2820
INTEGER,PARAMETER :: id_t_obs=3073
INTEGER,PARAMETER :: id_s_obs=5521      !(OCEAN)
INTEGER,PARAMETER :: id_ssh_obs=5526    !(OCEAN)
INTEGER,PARAMETER :: id_sst_obs=5525    !(OCEAN)
INTEGER,PARAMETER :: id_sss_obs=5522    !(OCEAN)
INTEGER,PARAMETER :: id_eta_obs=5351    !(OCEAN)
INTEGER,PARAMETER :: id_sic_obs=5283    !(SEAICE) sea ice concentration (fraction)
INTEGER,PARAMETER :: nid_dobs=3
INTEGER,PARAMETER :: id_x_obs=1111   !(OCEAN) (DRIFTERS) !STEVE: may want to change this depending on type of drifters
INTEGER,PARAMETER :: id_y_obs=2222   !(OCEAN) (DRIFTERS) !STEVE: may want to change this depending on type of drifters
INTEGER,PARAMETER :: id_z_obs=3333   !(OCEAN) (DRIFTERS) !STEVE: may want to change this depending on type of drifters
INTEGER,PARAMETER :: id_q_obs=3330
INTEGER,PARAMETER :: id_rh_obs=3331
INTEGER,PARAMETER :: id_ps_obs=14593
INTEGER,PARAMETER :: id_rain_obs=9999
INTEGER,PARAMETER :: nid_sfcflxobs=14       !(DO_SFCFLUXES)
INTEGER,PARAMETER :: id_atm_u_obs=5287      !(OCEAN) (ATMOS)
INTEGER,PARAMETER :: id_atm_v_obs=5284      !(OCEAN) (ATMOS)
INTEGER,PARAMETER :: id_atm_t_obs=5285      !(OCEAN) (ATMOS)
INTEGER,PARAMETER :: id_atm_q_obs=5281      !(OCEAN) (ATMOS)
INTEGER,PARAMETER :: id_atm_ps_obs=5280     !(OCEAN) (ATMOS)

!ISSUE: make these namelist controllable:
INTEGER :: obs1nrec = 6                     ! The number of records in the obs1-formatted file (previous 6, 7 adds a time record).
INTEGER :: obs2nrec = 9                     ! The number of records in the obs2-formatted file (previous 8, 9 adds a time record).

! Remove all observations above 65ºN due to tripolar grid
LOGICAL :: DO_REMOVE_65N = .false.
!STEVE: I made this a command line adjustable argument in obsop.f90 :: -rm65N .true.

LOGICAL            :: DO_POTTEMP_to_INSITU = .true.  ! Conversion to observation space. This is needed if the
                                                     ! observations aren't converted to potential temperature
                                                     ! (as is done by most - NCEP, SODA, NASA/GMAO, etc.). But
                                                     ! unlike that approach, this does not require synthetic salinity
                                                     ! observations to be constructed from climatologies.
                                                     ! This approach is theoretically better, but investigation must
                                                     ! be done to ensure model biases do not cause significant errors.
                                                     ! (a warning from J. Carton of potential difficulty)
                                                     !
                                                     ! Only one can be true, this one takes prioirty
                                                     !
LOGICAL            :: DO_INSITU_to_POTTEMP = .false. ! Technically, this would require matching an observed salinity
                                                     ! measurement with each observed in situ temperature measurement
                                                     ! and using it to compute the potential temperature. The opposite
                                                     ! process is quite a bit easier.

END MODULE params_obs

max_stars_repo_name: GEOS-ESM/Ocean-LETKF

max_issues_repo_name: GEOS-ESM/Ocean-LETKF

max_forks_repo_name: GEOS-ESM/Ocean-LETKF