text:
      SUBROUTINE PHO_HARISR(Ihpom,P1,P2,Ipf1,Ipf2,Ipa1,Ipa2,Iv1,Iv2,Q2h,
     &                      Xh1,Xh2,Xhmax1,Xhmax2,Ipb1,Ipb2,Ivo1,Ivo2,
     &                      Xisr1,Xisr2,Irej)
C********************************************************************
C
C     initial state radiation according to DGLAP evolution equations
C     (backward evolution, no spin effects)
C
C     input:    IHPOM     index of hard Pomeron
C                         negative: delete all previous entries
C               P1,P2     4 momenta of hard scattered final partons
C                         (in CMS of hard scattering)
C               IPF1,2    flavours of final partons
C               IPA1,2    flavours of initial partons
C               IV1,2     valence quark labels (0/1)
C               Q2H       momentum transfer (squared, positive)
C               XH1,XH2   x values of initial partons
C               XHMAX1,2  max. x values allowed
C
C     output:   all emitted partons in /POPISR/, final state
C               partons are the first two entries
C               shower evolution traced in /PODGL1/
C               IPB1,2    flavours of new initial partons
C               XISR1,2   x values of new initial partons
C               IVO1,2    valence quark labels (0/1)
C
C     attention: quark numbering according to PDG convention,
C                but 0 for gluons
C
C********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION anorf , b0qcd , beg , cf , cfe , codd , cofd ,
     &                 cog , coh , DEPS , DT_RNDM , ee , ee3 , ee4 ,
     &                 fxp , gam , gb , P1 , P2 , pc
      DOUBLE PRECISION pd1 , pd2 , PHO_PMASS , pje , pjx , pjy , pjz ,
     &                 pm , pn , pp4 , pt2 , pt2new , pt3 , pta1 ,
     &                 ptot1 , px , px3 , py , py3 , pz
      DOUBLE PRECISION pz3 , pz4 , q2 , Q2h , q2new , q2p , qmax , r1 ,
     &                 r3 , RHOMAS , s1 , s3 , scale2 , sfe , sh ,
     &                 shat1 , shz , sidd , sifd , sig
      DOUBLE PRECISION sih , ssh , sshz , theta , TINY , tt , wgdir ,
     &                 wgf , wggap , wgpdf , wgtot , Xh1 , Xh2 , xhma ,
     &                 Xhmax1 , Xhmax2 , xhmi , xi , Xisr1 , Xisr2
      DOUBLE PRECISION xmip , xms4 , xms4m , xnew , xp , zmax , zmin ,
     &                 zz
      INTEGER i , idmo , ifla , iflb , ifsum , ihidx , Ihpom , ii , il ,
     &        indx , ip , ipa , Ipa1 , Ipa2 , ipal , ipb , Ipb1 , Ipb2 ,
     &        ipdfc , Ipf1
      INTEGER Ipf2 , Irej , Iv1 , Iv2 , ival , Ivo1 , Ivo2 , k , kf ,
     &        kk , l , next , ngen , niter , ntry
      SAVE

      PARAMETER (RHOMAS=0.766D0,DEPS=1.D-10,TINY=1.D-10)

      DIMENSION P1(4) , P2(4)

C  input/output channels
      INCLUDE 'inc/poinou'
C  event debugging information
      INCLUDE 'inc/podebg'
C  internal rejection counters
      INCLUDE 'inc/poloop'
C  model switches and parameters
      INCLUDE 'inc/pomdls'
C  data of c.m. system of Pomeron / Reggeon exchange
      INCLUDE 'inc/popcms'
C  some hadron information, will be deleted in future versions
      INCLUDE 'inc/pohdrn'
C  currently activated parton density parametrizations
      INCLUDE 'inc/poppdf'
C  scale parameters for parton model calculations
      INCLUDE 'inc/pohscl'
C  parameters for DGLAP backward evolution in ISR
      INCLUDE 'inc/podgl1'
C  initial state parton radiation (internal part)
      INCLUDE 'inc/point6'
C  some constants
      INCLUDE 'inc/pocons'
C  particles created by initial state evolution
      INCLUDE 'inc/popisr'

      DOUBLE PRECISION PYP , eer , ther , qmaxr
      INTEGER PYK

      DIMENSION xhma(2) , next(2) , pd1(-6:6) , pd2(-6:6) , wggap(-6:6)
     &          , wgpdf(-6:6) , xhmi(2) , gb(4) , pm(4) , pn(4) ,
     &          pc(2,4) , q2(2) , ival(2) , ipal(2) , il(2) , ifsum(2) ,
     &          idmo(2)

      Irej = 0
      ntry = 1000
      niter = 0
C  debug output
      IF ( IDEb(79).GE.10 ) THEN
         IF ( LPRi.GT.4 )
     &         WRITE (LO,'(1X,A,/1X,I10,3I3,5E11.3,2(/5X,4E12.3))')
     &         'PHO_HARISR: KEV,IHPOM,IP1,IP2,Q2H,XH1,XH2,XHM1,XHM2:' ,
     &        KEVent , Ihpom , Ipa1 , Ipa2 , Q2h , Xh1 , Xh2 , Xhmax1 ,
     &        Xhmax2 , P1 , P2
      END IF
      IF ( Ihpom.EQ.0 ) RETURN
C
 100  NACc = 0
      idmo(1) = IDPdg1
      idmo(2) = IDPdg2
C
C  copy final state partons to local fields
      ihidx = ABS(Ihpom)

      IF ( ihidx.GT.MXISR2 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,2I4)')
     &         'PHO_HARISR: no space left in ' ,
     &        '/POPISR/ for hard scattering labels (IHIDX,MXISR2):' ,
     &        ihidx , MXISR2
         Irej = 1
      END IF

      DO k = 1 , 2
         IF ( Ihpom.LT.0 ) IMXisr(k) = 0
         IPOisr(k,1,ihidx) = IMXisr(k) + 1
         ipal(k) = IPOisr(k,1,ihidx)
      END DO
      DO i = 1 , 4
         PHIsr(1,i,ipal(1)) = P1(i)
         PHIsr(2,i,ipal(2)) = P2(i)
      END DO
      IFLisr(1,ipal(1)) = Ipf1
      IFLisr(2,ipal(2)) = Ipf2
C
C  check limitations, initialize /PODGL1/
      IF ( (Q2h.GT.Q2Misr(1)) .AND. (Xh1.LT.Xhmax1) ) THEN
         next(1) = 1
         Q2Sh(1,1) = Q2h
      ELSE
         next(1) = 0
         Q2Sh(1,1) = 0.D0
      END IF
      IF ( (Q2h.GT.Q2Misr(2)) .AND. (Xh2.LT.Xhmax2) ) THEN
         next(2) = 1
         Q2Sh(2,1) = Q2h
      ELSE
         next(2) = 0
         Q2Sh(2,1) = 0.D0
      END IF
C
      ISH(1) = 1
      ISH(2) = 1
      XPSh(1,1) = Xh1
      XPSh(2,1) = Xh2
C
      IFL1(1,1) = Ipa1
      ival(1) = Iv1
      IF ( (Ipa1.EQ.22) .OR. (Ipa1.EQ.990) ) next(1) = 0
      IFL1(2,1) = Ipa2
      ival(2) = Iv2
      IF ( (Ipa2.EQ.22) .OR. (Ipa2.EQ.990) ) next(2) = 0
C
      IF ( LPRi.GT.4 .AND. IDEb(79).GE.17 )
     &      WRITE (LO,'(1X,A,/5X,2I2,3E12.3)')
     &      'PHO_HARISR:INITIAL TESTS (NEXT1,2 Q2H Q21,2)' , next ,
     &     Q2h , Q2Misr
      IF ( next(1)+next(2).EQ.0 ) GOTO 600
C
C  initialize parton shower loop
      b0qcd = (33.D0-2.D0*NFSisr)/6.D0
      AL2isr(1) = PDFlam(1)
      AL2isr(2) = PDFlam(2)
      xhma(1) = Xhmax1
      xhma(2) = Xhmax2
      xhmi(1) = PMIsr(1)/PCMp
      xhmi(2) = PMIsr(2)/PCMp
      ZPSh(1,1) = 1.D0
      ZPSh(2,1) = 1.D0
      shat1 = Xh1*Xh2*ECMp**2
      IF ( IPAmdl(109).EQ.1 ) THEN
         PT2sh(1,1) = Q2h
      ELSE
         PT2sh(1,1) = Q2h*(1.D0-Q2h/shat1)
      END IF
      PT2sh(2,1) = PT2sh(1,1)
      IF ( PT2sh(1,1).LT.Q2Misr(1) ) next(1) = 0
      IF ( PT2sh(2,1).LT.Q2Misr(2) ) next(2) = 0
      THSh(1,1) = 2.D0*SQRT(Q2h/shat1)
      THSh(2,1) = THSh(1,1)
      IFAno(1) = 0
      IFAno(2) = 0
      zz = 1.D0
      IF ( Irej.NE.0 ) GOTO 600
C
C  main generation loop
C -------------------------------------------------
C  choose parton side to become solved
 200  IF ( (next(1)+next(2)).EQ.2 ) THEN
         IF ( Q2Sh(1,ISH(1)).GT.Q2Sh(2,ISH(2)) ) THEN
            ip = 1
         ELSE IF ( Q2Sh(2,ISH(2)).GT.Q2Sh(1,ISH(1)) ) THEN
            ip = 2
         ELSE
            ip = MAX(INT(DT_RNDM(shat1)*2.D0+0.999999D0),1)
         END IF
      ELSE IF ( next(1).EQ.1 ) THEN
         ip = 1
      ELSE IF ( next(2).EQ.1 ) THEN
         ip = 2
      ELSE
         GOTO 600
      END IF
      indx = ISH(ip)
C  INDX now parton position of parton to become solved
C  IP   now side to be treated
      xp = XPSh(ip,indx)
      q2p = Q2Sh(ip,indx)
      pt2 = PT2sh(ip,indx)
      iflb = IFL1(ip,indx)
C  check available x
      xmip = xhmi(ip)
C  cutoff by x limitation: no further development
      IF ( (xhma(ip)-xp).LT.xmip*2.D0 ) THEN
         next(ip) = 0
         Q2Sh(ip,indx) = 0.D0
         IF ( IDEb(79).GE.17 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,/5X,3E12.4,2I3)')
     &            'PHO_HARISR: evolution x-stop (XP,XMIP,XHMA,IP,INDX)'
     &           , xp , xmip , xhma(ip) , ip , indx
         END IF
         GOTO 200
      END IF
C  initial value of evolution variable t
      tt = LOG(AQQali*q2p/AL2isr(ip))
      DO i = -NFSisr , NFSisr
         wggap(i) = 0.D0
         wgpdf(i) = 0.D0
      END DO
C  DGLAP weights
      zmin = xp/xhma(ip)
      zmax = xp/(xp+xmip)
      cf = 4./3.
C  q --> q g, g --> g g
      IF ( iflb.EQ.0 ) THEN
         wggap(0) = cf*((zmax**2-zmin**2)/2.D0-2.D0*(zmax-zmin)
     &              +2.D0*LOG(zmax/zmin))
         DO i = 1 , NFSisr
            wggap(i) = wggap(0)
            wggap(-i) = wggap(0)
         END DO
         wggap(0) = 6.D0*((zmin**3-zmax**3)/3.D0+(zmax**2-zmin**2)
     &              /2.D0-2.D0*(zmax-zmin)
     &              +LOG(zmax/zmin*(1.D0-zmin)/(1.D0-zmax)))
C  q --> g q, g --> q qb
      ELSE IF ( ABS(iflb).LE.6 ) THEN
         wggap(iflb) = cf*((zmin**2-zmax**2)/2.D0-zmax+zmin-2.D0*LOG((
     &                 1.D0-zmax)/(1.D0-zmin)))
         IF ( ival(ip).EQ.0 ) wggap(0)
     &        = 0.5D0*(2./3.*(zmax**3-zmin**3)-zmax**2+zmin**2+zmax-
     &        zmin)
      ELSE
         IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A,I7)')
     &         'PHO_HARISR:ERROR: unsupported particle ID' , iflb
         CALL PHO_ABORT
      END IF
C  anomalous/resolved evolution
      ipdfc = 0
      IF ( IPAmdl(110).GE.1 ) THEN
         IF ( (idmo(ip).EQ.22) .AND. (iflb.NE.0) .AND. (iflb.NE.21) )
     &        THEN
            wgdir = 0.D0
            IF ( NQQali.EQ.1 ) THEN
               scale2 = pt2*AQQpd
            ELSE
               scale2 = q2p*AQQpd
            END IF
            CALL PHO_PDF(ip,xp,scale2,0.D0,pd1)
            ipdfc = 1
            CALL PHO_QPMPDF(iflb,xp,scale2,pta1,PVIrtp(ip),wgdir)
            xi = DT_RNDM(xp)*pd1(iflb)
            IF ( wgdir.GT.xi ) THEN
C  debug output
               IF ( LPRi.GT.4 .AND. IDEb(79).GE.17 )
     &               WRITE (LO,'(1X,2A,/5X,4E12.5,I2,I3)')
     &               'PHO_HARISR: ' ,
     &              'direct splitting (WGDIR,WGPDF,X,SCALE2,IP,IFLB)' ,
     &              wgdir , pd1(iflb) , xp , scale2 , ip , iflb
               Q2Sh(ip,indx) = 0.D0
               next(ip) = 0
               IFAno(ip) = indx
               GOTO 200
            END IF
         END IF
      END IF
C
C  rejection loop for z,t sampling
C ------------------------------------
 300  niter = niter + 1
      IF ( niter.GE.ntry ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,2I6)')
     &         'PHO_HARISR: too many rejections' , niter , ntry
         CALL PHO_PREVNT(-1)
C  clean up event
         Irej = 1
         GOTO 100
      END IF
C  PDF weights
      IF ( ipdfc.EQ.0 ) THEN
         IF ( NQQali.EQ.1 ) THEN
            scale2 = pt2*AQQpd
         ELSE
            scale2 = q2p*AQQpd
         END IF
         CALL PHO_PDF(ip,xp,scale2,0.D0,pd1)
      END IF
      ipdfc = 0
C
      wgtot = 0.D0
      DO i = -NFSisr , NFSisr
         wgpdf(i) = pd1(i)/(pd1(iflb)+1.D-12)*5.D0
         wgtot = wgtot + wgpdf(i)*wggap(i)
      END DO
C
C  sample new t value
 400  tt = tt*EXP(MAX(-10.D0,LOG(DT_RNDM(shat1))*b0qcd/wgtot))
      q2new = AL2isr(ip)*EXP(tt)/AQQali
C  debug output
      IF ( LPRi.GT.4 .AND. IDEb(79).GE.20 ) WRITE (LO,'(1X,A,E12.5)')
     &      'PHO_HARISR: pre-selected Q2:' , q2new
C  compare to limits
      IF ( q2new.LT.Q2Misr(ip) ) THEN
         Q2Sh(ip,indx) = 0.D0
         next(ip) = 0
         IF ( LPRi.GT.4 .AND. IDEb(79).GE.17 )
     &         WRITE (LO,'(1X,A,2E10.3,2I3)')
     &         'PHO_HARISR: evolution Q2-stop (Q2,Q2MIN,IP,INDX):' ,
     &        q2new , Q2Misr(ip) , ip , indx
         GOTO 200
      END IF
      Q2Sh(ip,indx) = q2new
      tt = LOG(AQQali*q2new/AL2isr(ip))
C  selection of flavours
      xi = wgtot*DT_RNDM(tt)
      ifla = -NFSisr - 1
 500  ifla = ifla + 1
      xi = xi - wgpdf(ifla)*wggap(ifla)
      IF ( (xi.GT.0.D0) .AND. (ifla.LT.NFSisr) ) GOTO 500
C  debug output
C  selection of z
      IF ( LPRi.GT.4 .AND. IDEb(79).GE.20 ) WRITE (LO,'(1X,A,2I3)')
     &      'PHO_HARISR: pre-selected IFLA (IFLA,IFLB):' , ifla , iflb
      CALL PHO_HARZSP(ifla,iflb,NFSisr,zmin,zmax,zz)
C  debug output
      IF ( LPRi.GT.4 .AND. IDEb(79).GE.20 ) WRITE (LO,'(1X,A,E12.3)')
     &      'PHO_HARISR: pre-selected ZZ' , zz
C  angular ordering
      theta = 4.D0*zz**2*q2new/((ECMp*xp)**2*(1.D0-zz))
      IF ( theta.GT.THSh(ip,indx) ) THEN
         IF ( LPRi.GT.4 .AND. IDEb(79).GE.20 )
     &        WRITE (LO,'(1X,A,2E12.3)')
     &         'PHO_HARISR: reject by angle (NEW/OLD)' , theta ,
     &        THSh(ip,indx)
         GOTO 400
      END IF
C  rejection weight given by new PDFs
      xnew = xp/zz
      pt2new = q2new*(1.D0-zz)
      IF ( NQQali.EQ.1 ) THEN
         scale2 = pt2new*AQQpd
      ELSE
         scale2 = q2new*AQQpd
      END IF
      IF ( scale2.LT.Q2Misr(ip) ) THEN
         Q2Sh(ip,indx) = 0.D0
         next(ip) = 0
         IF ( LPRi.GT.4 .AND. IDEb(79).GE.17 )
     &         WRITE (LO,'(1X,A,2E10.3,2I3)')
     &         'PHO_HARISR: evol.Q2-stop (SCALE2,Q2MIN,IP,INDX):' ,
     &        q2new , Q2Misr(ip) , ip , indx
         GOTO 200
      END IF
      CALL PHO_PDF(ip,xnew,scale2,0.D0,pd2)
      IF ( pd2(ifla).LT.1.D-10 ) GOTO 300
      CALL PHO_PDF(ip,xp,scale2,0.D0,pd1)
      pd1(iflb) = MAX(pd1(iflb),1.D-10)
      wgf = pd2(ifla)/pd1(iflb)/(wgpdf(ifla)+1.D-12)
      IF ( NQQali.EQ.1 ) wgf = wgf*LOG(q2new*AQQali/AL2isr(ip))
     &     /LOG(pt2new*AQQali/AL2isr(ip))
      IF ( (wgf.GT.1.D0) .AND. (IDEb(79).GE.2) ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,E12.3)')
     &         'PHO_HARISR: final weight:' , wgf
         IF ( LPRi.GT.4 ) WRITE (LO,'(6X,A,I7,2I3,3E11.3)')
     &         'EV,IFLA,IFLB,Q2,PT2,Z:' , KEVent , ifla , iflb , q2new ,
     &        pt2new , zz
      END IF
      IF ( wgf.LT.DT_RNDM(xnew) ) GOTO 300

      IF ( IDEb(79).GE.15 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,/3X,3I3,3E11.3)')
     &         'PHO_HARISR: accepted IP,IFLA,IFLB,PT2,Q2,Z:' , ip ,
     &        ifla , iflb , pt2new , q2new , zz
      END IF

      IF ( indx.GE.MXISR3 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,2I4)')
     &         'PHO_HARISR: no space left in ' ,
     &        '/POINT6/ for parton radiation (INDX,MXISR3):' , indx ,
     &        MXISR3
         Irej = 1
         RETURN
      END IF

C  branching accepted, registration
      Q2Sh(ip,indx) = q2new
      PT2sh(ip,indx) = pt2new
      ZPSh(ip,indx) = zz
      IFL2(ip,indx) = ifla - iflb
      Q2Sh(ip,indx+1) = q2new
      PT2sh(ip,indx+1) = PT2sh(ip,indx)
      XPSh(ip,indx+1) = xnew
      THSh(ip,indx+1) = theta
      IFL1(ip,indx+1) = ifla
      ISH(ip) = ISH(ip) + 1

      NACc = NACc + 1

      IF ( NACc.GT.MXISR4 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,2I4)')
     &         'PHO_HARISR: no space left in ' ,
     &        '/POINT6/ for parton radiation (NACC,MXISR4):' , NACc ,
     &        MXISR4
         Irej = 1
         RETURN
      END IF

      SHAt(NACc) = shat1
      IBRa(1,NACc) = ip
      IBRa(2,NACc) = indx
      shat1 = shat1/zz

C  generation of next branching
      IF ( next(1)+next(2).NE.0 ) GOTO 200

C  new initial flavours, x values
 600  Ipb1 = IFL1(1,ISH(1))
      Ipb2 = IFL1(2,ISH(2))
      Xisr1 = XPSh(1,ISH(1))
      Xisr2 = XPSh(2,ISH(2))
      Ivo1 = ival(1)
      Ivo2 = ival(2)
C  valence flavours
      IF ( Ipb1.NE.0 ) THEN
         IF ( ISH(1).GT.1 ) THEN
            CALL PHO_PDF(1,Xisr1,Q2Misr(1),0.D0,pd1)
            IF ( IDPdg1.EQ.22 ) THEN
               CALL PHO_QPMPDF(Ipb1,Xisr1,Q2Misr(1),0.D0,PVIrtp(1),fxp)
               IF ( DT_RNDM(Xisr1)*pd1(Ipb1).GT.pd1(Ipb1)-fxp ) ival(1)
     &              = 1
            ELSE
               CALL PHO_PDF(1,Xisr1,Q2Misr(1),0.D0,pd1)
               IF ( DT_RNDM(Xisr1)*pd1(Ipb1).GT.pd1(-Ipb1) ) ival(1) = 1
            END IF
         END IF
      END IF
      IF ( Ipb2.NE.0 ) THEN
         IF ( ISH(2).GT.1 ) THEN
            CALL PHO_PDF(2,Xisr2,Q2Misr(2),0.D0,pd1)
            IF ( IDPdg2.EQ.22 ) THEN
               CALL PHO_QPMPDF(Ipb2,Xisr2,Q2Misr(2),0.D0,PVIrtp(2),fxp)
               IF ( DT_RNDM(Xisr2)*pd1(Ipb2).GT.pd1(Ipb2)-fxp ) ival(2)
     &              = 1
            ELSE
               IF ( DT_RNDM(Xisr2)*pd1(Ipb2).GT.pd1(-Ipb2) ) ival(2) = 1
            END IF
         END IF
      END IF

C  parton kinematics
      IF ( NACc.GT.0 ) THEN
C  final partons in CMS
         pm(3) = (Xh1-Xh2)*ECMp/2.D0
         pm(4) = (Xh1+Xh2)*ECMp/2.D0
         sh = Xh1*Xh2*ECMp**2
         ssh = SQRT(sh)
         gb(3) = pm(3)/ssh
         gb(4) = pm(4)/ssh
         CALL PHO_ALTRA(gb(4),0.D0,0.D0,-gb(3),P1(1),P1(2),P1(3),P1(4),
     &                  ptot1,PHIsr(1,1,ipal(1)),PHIsr(1,2,ipal(1)),
     &                  PHIsr(1,3,ipal(1)),PHIsr(1,4,ipal(1)))
         CALL PHO_ALTRA(gb(4),0.D0,0.D0,-gb(3),P2(1),P2(2),P2(3),P2(4),
     &                  ptot1,PHIsr(2,1,ipal(2)),PHIsr(2,2,ipal(2)),
     &                  PHIsr(2,3,ipal(2)),PHIsr(2,4,ipal(2)))
         il(1) = 1
         il(2) = 1
         DO i = 1 , NACc
            ipa = IBRa(1,i)
            ipb = 3 - ipa
            il(ipa) = IBRa(2,i)
C  new initial partons in CMS
            sh = SHAt(i)
            ssh = SQRT(sh)
            shz = sh/ZPSh(ipa,il(ipa))
            sshz = SQRT(shz)
            q2(1) = Q2Sh(1,il(1))
            q2(2) = Q2Sh(2,il(2))
            pc(1,1) = 0.D0
            pc(1,2) = 0.D0
            pc(1,3) = SQRT((sh+q2(1)+q2(2))**2-4.D0*q2(1)*q2(2))
     &                /(2.D0*ssh)
            pc(1,4) = (sh-q2(1)+q2(2))/(2.D0*ssh)
            pc(2,1) = 0.D0
            pc(2,2) = 0.D0
            pc(2,3) = -pc(1,3)
            pc(2,4) = ssh - pc(1,4)
            xms4 = PHO_PMASS(IFL2(ipa,il(ipa)),1)**2
            ee3 = (shz-q2(ipa)+q2(ipb)-xms4)/(2.D0*ssh)
            s1 = sh + q2(ipa) + q2(ipb)
            s3 = shz + q2(ipb) + Q2Sh(ipa,il(ipa)+1)
            r1 = SQRT(s1**2-4.D0*q2(ipa)*q2(ipb))
            r3 = SQRT(s3**2-4.D0*q2(ipb)*Q2Sh(ipa,il(ipa)+1))
            IF ( q2(ipb).LT.0.1D0 ) THEN
               xms4m = (q2(ipa)/ZPSh(ipa,il(ipa))-Q2Sh(ipa,il(ipa)+1))
     &                 *(sh/(sh+q2(ipa))-sh/(shz+Q2Sh(ipa,il(ipa)+1)))
            ELSE
               xms4m = (s1*s3-r1*r3)/(2.D0*q2(ipb)) - q2(ipa)
     &                 - Q2Sh(ipa,il(ipa)+1)
            END IF
            ngen = 1
C  max. virtuality for time-like showers
            qmax = MIN(xms4m,PARmdl(95)*q2(ipa))
            IF ( (IPAmdl(111).GE.1) .AND. (qmax.GT.PARmdl(94)) ) THEN
C  generate time-like parton shower
               kf = IFL2(ipa,il(ipa))
               IF ( kf.EQ.0 ) kf = 21
               eer = MIN(ee3-pc(ipa,4),ECMp)
               ther = 0.

               CALL PY1ENT(1,kf,eer,ther,ther)
               qmaxr = SQRT(qmax)
               CALL PYSHOW(1,0,qmaxr)
C debug output
               IF ( IDEb(79).GE.25 ) THEN
                  IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,/,5X,1P,4E12.4)')
     &                  'PHO_HARISR: ' ,
     &                 'PYSHOW called for EE,QMAX,XMS4M,Q2(IPA)' , eer ,
     &                 qmax , xms4m , q2(ipa)
                  CALL PYLIST(1)
               END IF
               ngen = PYK(0,1)

               IF ( ngen.GT.1 ) THEN
                  pjx = 0.D0
                  pjy = 0.D0
                  pjz = 0.D0
                  pje = 0.D0
                  kk = ipal(ipa)
                  DO k = 3 , ngen

                     IF ( PYK(k,1).LE.4 ) THEN
                        kk = kk + 1

                        IF ( kk.GT.MXISR1 ) THEN
                           IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,2I5)')
     &                         'PHO_HARISR: no space ' ,
     &                        'left in /POPISR/ (KK,MXISR1):' , kk ,
     &                        MXISR1
                           Irej = 1
                           RETURN
                        END IF

                        PHIsr(ipa,1,kk) = PYP(k,1)
                        pjx = pjx + PHIsr(ipa,1,kk)
                        PHIsr(ipa,2,kk) = PYP(k,2)
                        pjy = pjy + PHIsr(ipa,2,kk)
                        PHIsr(ipa,3,kk) = PYP(k,3)
                        pjz = pjz + PHIsr(ipa,3,kk)
                        PHIsr(ipa,4,kk) = PYP(k,4)
                        pje = pje + PHIsr(ipa,4,kk)
                        IFLisr(ipa,kk) = PYK(k,2)

                        IF ( IFLisr(ipa,kk).EQ.21 ) IFLisr(ipa,kk) = 0
                        IF ( IFLisr(ipa,kk).EQ.5 ) IFLisr(ipa,kk) = 3
                        IF ( IFLisr(ipa,kk).EQ.-5 ) IFLisr(ipa,kk) = -3
                     END IF
                  END DO
                  ngen = kk - ipal(ipa)
                  xms4 = (pje+pjz)*(pje-pjz) - pjx**2 - pjy**2
                  pp4 = SQRT(pje**2-xms4)
                  ee3 = (shz-q2(ipa)+q2(ipb)-xms4)/(2.D0*ssh)
C debug output
                  IF ( LPRi.GT.4 .AND. IDEb(79).GE.20 )
     &                  WRITE (LO,'(1X,2A,/,5X,1P,6E12.4)')
     &                  'PHO_HARISR: ' ,
     &                 'time-like shower: PJE,PJX,PJY,PJZ,PP4,XMS4' ,
     &                 pje , pjx , pjy , pjz , pp4 , xms4
               END IF
            END IF
            pz3 = (2.D0*pc(ipa,4)*ee3+q2(ipa)+Q2Sh(ipa,il(ipa)+1)+xms4)
     &            /(2.D0*pc(ipa,3))
            pt3 = (ee3+pz3)*(ee3-pz3) + Q2Sh(ipa,il(ipa)+1)
            IF ( pt3.LT.0.D0 ) THEN
               IF ( LPRi.GT.4 .AND. IDEb(79).GE.5 )
     &               WRITE (LO,'(1X,A,E12.3)')
     &               'PHO_HARISR: rejection due to PT3' , pt3
               GOTO 100
            END IF
            pt3 = SQRT(pt3)
            CALL PHO_SFECFE(sfe,cfe)
            px3 = cfe*pt3
            py3 = sfe*pt3
C
            IF ( ngen.GT.1 ) THEN
C  time-like shower generated
               ee4 = ee3 - pc(ipa,4)
               pz4 = pz3 - pc(ipa,3)
               pp4 = SQRT(pt3**2+pz4**2)
C  Lorentz boost
               gam = (ee4*pje-pp4*pjz)/xms4
               beg = (pje*pp4-ee4*pjz)/xms4
C  rotation angles
               codd = pz4/pp4
               sidd = SQRT(px3**2+py3**2)/pp4
               cofd = 1.D0
               sifd = 0.D0
               IF ( pp4*sidd.GT.1.D-5 ) THEN
                  cofd = px3/(sidd*pp4)
                  sifd = py3/(sidd*pp4)
                  anorf = SQRT(cofd*cofd+sifd*sifd)
                  cofd = cofd/anorf
                  sifd = sifd/anorf
               END IF
C  copy partons back
               kk = ipal(ipa)
               DO k = 1 , ngen
                  kk = kk + 1
                  px = PHIsr(ipa,1,kk)
                  py = PHIsr(ipa,2,kk)
                  pz = PHIsr(ipa,3,kk)
                  coh = PHIsr(ipa,4,kk)
                  ee = gam*coh + beg*pz
                  pz = gam*pz + beg*coh
                  PHIsr(ipa,4,kk) = ee
                  CALL PHO_TRANS(px,py,pz,codd,sidd,cofd,sifd,
     &               PHIsr(ipa,1,kk),PHIsr(ipa,2,kk),PHIsr(ipa,3,kk))
               END DO
               ipal(ipa) = kk
            ELSE
C  no time-like shower generated
               ipal(ipa) = ipal(ipa) + 1
               PHIsr(ipa,1,ipal(ipa)) = px3
               PHIsr(ipa,2,ipal(ipa)) = py3
               PHIsr(ipa,3,ipal(ipa)) = pz3 - pc(ipa,3)
               PHIsr(ipa,4,ipal(ipa)) = ee3 - pc(ipa,4)
               IFLisr(ipa,ipal(ipa)) = IFL2(ipa,il(ipa))
            END IF
            pc(ipa,1) = px3
            pc(ipa,2) = py3
            pc(ipa,3) = pz3
            pc(ipa,4) = ee3
C  boost / rotate into new CMS
            DO k = 1 , 4
               gb(k) = (pc(1,k)+pc(2,k))/sshz
            END DO
            CALL PHO_ALTRA(gb(4),-gb(1),-gb(2),-gb(3),pc(1,1),pc(1,2),
     &                     pc(1,3),pc(1,4),ptot1,pm(1),pm(2),pm(3),
     &                     pm(4))
            cog = pm(3)/ptot1
            sig = SQRT(pm(1)**2+pm(2)**2)/ptot1
            coh = 1.D0
            sih = 0.D0
            IF ( ptot1*sig.GT.1.D-5 ) THEN
               coh = pm(1)/(sig*ptot1)
               sih = pm(2)/(sig*ptot1)
               anorf = SQRT(coh*coh+sih*sih)
               coh = coh/anorf
               sih = sih/anorf
            END IF
            DO k = 1 , 2
               DO l = IPOisr(k,1,ihidx) , ipal(k)
                  CALL PHO_ALTRA(gb(4),-gb(1),-gb(2),-gb(3),PHIsr(k,1,l)
     &               ,PHIsr(k,2,l),PHIsr(k,3,l),PHIsr(k,4,l),ptot1,
     &               pm(1),pm(2),pm(3),pm(4))
                  CALL PHO_TRANI(pm(1),pm(2),pm(3),cog,sig,coh,sih,
     &               pn(1),pn(2),pn(3))
                  CALL PHO_TRANS(pn(1),pn(2),pn(3),1.D0,0.D0,coh,sih,
     &               PHIsr(k,1,l),PHIsr(k,2,l),PHIsr(k,3,l))
                  PHIsr(k,4,l) = pm(4)
               END DO
            END DO
         END DO
C  boost back to global CMS
         pm(3) = (Xisr1-Xisr2)/2.D0
         pm(4) = (Xisr1+Xisr2)/2.D0
         ssh = SQRT(Xisr1*Xisr2)
         gb(3) = pm(3)/ssh
         gb(4) = pm(4)/ssh
         DO k = 1 , 2
            DO l = IPOisr(k,1,ihidx) , ipal(k)
               CALL PHO_ALTRA(gb(4),0.D0,0.D0,gb(3),PHIsr(k,1,l),
     &                        PHIsr(k,2,l),PHIsr(k,3,l),PHIsr(k,4,l),
     &                        ptot1,pm(1),pm(2),pm(3),pm(4))
               PHIsr(k,1,l) = pm(1)
               PHIsr(k,2,l) = pm(2)
               PHIsr(k,3,l) = pm(3)
               PHIsr(k,4,l) = pm(4)
            END DO
         END DO
      END IF
      IPOisr(1,2,ihidx) = ipal(1)
      IPOisr(2,2,ihidx) = ipal(2)
      IMXisr(1) = ipal(1)
      IMXisr(2) = ipal(2)
C
C  debug output
      IF ( IDEb(79).GE.10 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,2I10/,6X,A,2E12.3,2I5)')
     &         'NUMBER OF EMISSIONS' , ISH(1) - 1 , ISH(2) - 1 ,
     &        'NEW X1,X2,IFL1,ILF2' , Xisr1 , Xisr2 , Ipb1 , Ipb2
         IF ( NACc.GT.0 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,2I5,/6X,A)')
     &            'PHO_HARISR: ISR configuration (NITER,NACC)' , niter ,
     &           NACc ,
     &     ' SIDE   NO.   IFLB IFLC     Q2SH    PT2SH     XH         ZZ'
            DO ii = 1 , NACc
               k = IBRa(1,ii)
               i = IBRa(2,ii)
               IF ( LPRi.GT.4 ) WRITE (LO,'(5X,4I5,4E11.3)') k , i ,
     &              IFL1(k,i) , IFL2(k,i) , Q2Sh(k,i) , PT2sh(k,i) ,
     &              XPSh(k,i) , ZPSh(k,i)
            END DO
         END IF
C  check of final configuration
         px3 = 0.D0
         py3 = 0.D0
         pz3 = 0.D0
         ee3 = 0.D0
         ifsum(1) = 0
         ifsum(2) = 0
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A)')
     &         'PHO_HARISR: outgoing partons'
         DO k = 1 , 2
            DO l = IPOisr(k,1,ihidx) , IPOisr(k,2,ihidx)
               IF ( LPRi.GT.4 ) WRITE (LO,'(6X,2I4,I6,4E11.3)') k , l ,
     &              IFLisr(k,l) , PHIsr(k,1,l) , PHIsr(k,2,l) ,
     &              PHIsr(k,3,l) , PHIsr(k,4,l)
               ifsum(k) = ifsum(k) + IFLisr(k,l)
               px3 = px3 + PHIsr(k,1,l)
               py3 = py3 + PHIsr(k,2,l)
               pz3 = pz3 + PHIsr(k,3,l)
               ee3 = ee3 + PHIsr(k,4,l)
            END DO
         END DO
         ifsum(1) = ifsum(1) - Ipb1
         ifsum(2) = ifsum(2) - Ipb2
         pz3 = pz3 - (Xisr1-Xisr2)*ECMp/2.D0
         ee3 = ee3 - (Xisr1+Xisr2)*ECMp/2.D0
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,2I4,4E11.3)')
     &         'CHECK:IFL1,2 PCM(1-4)' , ifsum , px3 , py3 , pz3 , ee3
      END IF
      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      SUBROUTINE W3FI72(ITYPE,FLD,IFLD,IBITL,
     &                  IPFLAG,ID,PDS,
     &                  IGFLAG,IGRID,IGDS,ICOMP,
     &                  IBFLAG,IBMAP,IBLEN,IBDSFL,
     &                  NPTS,KBUF,ITOT,JERR)
C $$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:  W3FI72        MAKE A COMPLETE GRIB MESSAGE
C   PRGMMR: FARLEY           ORG: NMC421      DATE:94-11-22
C
C ABSTRACT: MAKES A COMPLETE GRIB MESSAGE FROM A USER SUPPLIED
C   ARRAY OF FLOATING POINT OR INTEGER DATA.  THE USER HAS THE
C   OPTION OF SUPPLYING THE PDS OR AN INTEGER ARRAY THAT WILL BE
C   USED TO CREATE A PDS (WITH W3FI68).  THE USER MUST ALSO
C   SUPPLY OTHER NECESSARY INFO; SEE USAGE SECTION BELOW.
C
C PROGRAM HISTORY LOG:
C   91-05-08  R.E.JONES
C   92-07-01  M. FARLEY    ADDED GDS AND BMS LOGIC.  PLACED EXISTING
C                          LOGIC FOR BDS IN A ROUTINE.
C   92-10-02  R.E.JONES    ADD ERROR EXIT FOR W3FI73
C   93-04-30  R.E.JONES    REPLACE DO LOOPS TO MOVE CHARACTER DATA
C                          WITH XMOVEX, USE XSTORE TO ZERO CHARACTER
C                          ARRAY. MAKE CHANGE SO FLAT FIELD WILL PACK.
C   93-08-06  CAVANAUGH    MODIFIED CALL TO W3FI75
C   93-10-26  CAVANAUGH    ADDED CODE TO RESTORE INPUT FIELD TO ORIGINAL
C                          VALUES IF D-SCALE NOT 0
C   94-01-27  CAVANAUGH    ADDED IGDS ARRAY IN CALL TO W3FI75 TO PROVIDE
C                          INFORMATION FOR BOUSTROPHEDONIC PROCESSING
C   94-03-03  CAVANAUGH    INCREASED SIZE OF GDS ARRAY FOR THIN GRIDS
C   94-05-16  FARLEY       CLEANED UP DOCUMENTATION
C   94-11-10  FARLEY       INCREASED SIZE OF PFLD/IFLD ARRARYS FROM
C                          100K TO 260K FOR .5 DEGREE SST ANAL FIELDS
C   94-12-04  R.E.JONES    CHANGE DOCUMENT FOR IPFLAG.
C   95-10-31  IREDELL      REMOVED SAVES AND PRINTS
C   98-05-19  Gilbert      Increased array dimensions to handle grids
C                          of up to 500,000 grid points.
C   95-10-31  IREDELL      GENERALIZED WORD SIZE
C   98-12-21  Gilbert      Replaced Function ICHAR with mov_a2i.
C   99-02-01  Gilbert      Changed the method of zeroing out array KBUF.
C                          the old method, using W3FI01 and XSTORE was
C                          incorrect with 4-byte integers and 8-byte reals.
C 2001-06-07  Gilbert      Removed calls to xmovex.
C                          changed IPFLD from integer to character.
C   10-02-19  GAYNO        FIX ALLOCATION OF ARRAY BMS
C
C USAGE:  CALL W3FI72(ITYPE,FLD,IFLD,IBITL,
C        &            IPFLAG,ID,PDS,
C        &            IGFLAG,IGRID,IGDS,ICOMP,
C        &            IBFLAG,IBMAP,IBLEN,IBDSFL,
C        &            IBDSFL,
C        &            NPTS,KBUF,ITOT,JERR)
C
C   INPUT ARGUMENT LIST:
C     ITYPE    - 0 = FLOATING POINT DATA SUPPLIED IN ARRAY 'FLD'
C                1 = INTEGER DATA SUPPLIED IN ARRAY 'IFLD'
C     FLD      - REAL ARRAY OF DATA (AT PROPER GRIDPOINTS) TO BE
C                CONVERTED TO GRIB FORMAT IF ITYPE=0.
C                SEE REMARKS #1 & 2.
C     IFLD     - INTEGER ARRAY OF DATA (AT PROPER GRIDPOINTS) TO BE
C                CONVERTED TO GRIB FORMAT IF ITYPE=1.
C                SEE REMARKS #1 & 2.
C     IBITL    - 0 = COMPUTER COMPUTES LENGTH FOR PACKING DATA FROM
C                    POWER OF 2 (NUMBER OF BITS) BEST FIT OF DATA
C                    USING 'VARIABLE' BIT PACKER W3FI58.
C                8, 12, ETC. COMPUTER RESCALES DATA TO FIT INTO THAT
C                    'FIXED' NUMBER OF BITS USING W3FI59.
C                SEE REMARKS #3.
C
C     IPFLAG   - 0 = MAKE PDS FROM USER SUPPLIED ARRAY (ID)
C                1 = USER SUPPLYING PDS
C                NOTE: IF PDS IS GREATER THAN 30, USE IPLFAG=1.
C                THE USER COULD CALL W3FI68 BEFORE HE CALLS
C                W3FI72. THIS WOULD MAKE THE FIRST 30 BYTES OF
C                THE PDS, USER THEN WOULD MAKE BYTES AFTER 30.
C     ID       - INTEGER ARRAY OF  VALUES THAT W3FI68 WILL USE
C                TO MAKE AN EDITION 1 PDS IF IPFLAG=0.  (SEE THE
C                DOCBLOCK FOR W3FI68 FOR LAYOUT OF ARRAY)
C     PDS      - CHARACTER ARRAY OF VALUES (VALID PDS SUPPLIED
C                BY USER) IF IPFLAG=1. LENGTH MAY EXCEED 28 BYTES
C                (CONTENTS OF BYTES BEYOND 28 ARE PASSED
C                THROUGH UNCHANGED).
C
C     IGFLAG   - 0 = MAKE GDS BASED ON 'IGRID' VALUE.
C                1 = MAKE GDS FROM USER SUPPLIED INFO IN 'IGDS'
C                    AND 'IGRID' VALUE.
C                SEE REMARKS #4.
C     IGRID    - #   = GRID IDENTIFICATION (TABLE B)
C                255 = IF USER DEFINED GRID; IGDS MUST BE SUPPLIED
C                      AND IGFLAG MUST =1.
C     IGDS     - INTEGER ARRAY CONTAINING USER GDS INFO (SAME
C                FORMAT AS SUPPLIED BY W3FI71 - SEE DOCKBLOCK FOR
C                LAYOUT) IF IGFLAG=1.
C     ICOMP    - RESOLUTION AND COMPONENT FLAG FOR BIT 5 OF GDS(17)
C                0 = EARTH ORIENTED WINDS
C                1 = GRID ORIENTED WINDS
C
C     IBFLAG   - 0 = MAKE BIT MAP FROM USER SUPPLIED DATA
C                # = BIT MAP PREDEFINED BY CENTER
C                SEE REMARKS #5.
C     IBMAP    - INTEGER ARRAY CONTAINING BIT MAP
C     IBLEN    - LENGTH OF BIT MAP WILL BE USED TO VERIFY LENGTH
C                OF FIELD (ERROR IF IT DOESN'T MATCH).
C
C     IBDSFL   - INTEGER ARRAY CONTAINING TABLE 11 FLAG INFO
C                BDS OCTET 4:
C                (1) 0 = GRID POINT DATA
C                    1 = SPHERICAL HARMONIC COEFFICIENTS
C                (2) 0 = SIMPLE PACKING
C                    1 = SECOND ORDER PACKING
C                (3) ... SAME VALUE AS 'ITYPE'
C                    0 = ORIGINAL DATA WERE FLOATING POINT VALUES
C                    1 = ORIGINAL DATA WERE INTEGER VALUES
C                (4) 0 = NO ADDITIONAL FLAGS AT OCTET 14
C                    1 = OCTET 14 CONTAINS FLAG BITS 5-12
C                (5) 0 = RESERVED - ALWAYS SET TO 0
C         BYTE 6 OPTION 1 NOT AVAILABLE (AS OF 5-16-93)
C                (6) 0 = SINGLE DATUM AT EACH GRID POINT
C                    1 = MATRIX OF VALUES AT EACH GRID POINT
C         BYTE 7 OPTION 0 WITH SECOND ORDER PACKING N/A (AS OF 5-16-93)
C                (7) 0 = NO SECONDARY BIT MAPS
C                    1 = SECONDARY BIT MAPS PRESENT
C                (8) 0 = SECOND ORDER VALUES HAVE CONSTANT WIDTH
C                    1 = SECOND ORDER VALUES HAVE DIFFERENT WIDTHS
C
C   OUTPUT ARGUMENT LIST:
C     NPTS     - NUMBER OF GRIDPOINTS IN ARRAY FLD OR IFLD
C     KBUF     - ENTIRE GRIB MESSAGE ('GRIB' TO '7777')
C                EQUIVALENCE TO INTEGER ARRAY TO MAKE SURE IT
C                IS ON WORD BOUNARY.
C     ITOT     - TOTAL LENGTH OF GRIB MESSAGE IN BYTES
C     JERR     - = 0, COMPLETED MAKING GRIB FIELD WITHOUT ERROR
C                  1, IPFLAG NOT 0 OR 1
C                  2, IGFLAG NOT 0 OR 1
C                  3, ERROR CONVERTING IEEE F.P. NUMBER TO IBM370 F.P.
C                  4, W3FI71 ERROR/IGRID NOT DEFINED
C                  5, W3FK74 ERROR/GRID REPRESENTATION TYPE NOT VALID
C                  6, GRID TOO LARGE FOR PACKER DIMENSION ARRAYS
C                     SEE AUTOMATION DIVISION FOR REVISION!
C                  7, LENGTH OF BIT MAP NOT EQUAL TO SIZE OF FLD/IFLD
C                  8, W3FI73 ERROR, ALL VALUES IN IBMAP ARE ZERO
C
C   OUTPUT FILES:
C     FT06F001 - STANDARD FORTRAN OUTPUT PRINT FILE
C
C   SUBPROGRAMS CALLED:
C     LIBRARY:
C       W3LIB    - W3FI58, W3FI59, W3FI68, W3FI71, W3FI73, W3FI74
C                  W3FI75, W3FI76
C       FORTRAN 90 INTRINSIC - BIT_SIZE
C
C REMARKS:
C   1)  IF BIT MAP TO BE INCLUDED IN MESSAGE, NULL DATA SHOULD
C       BE INCLUDED IN FLD OR IFLD.  THIS ROUTINE WILL TAKE CARE
C       OF 'DISCARDING' ANY NULL DATA BASED ON THE BIT MAP.
C   2)  UNITS MUST BE THOSE IN GRIB DOCUMENTATION:  NMC O.N. 388
C       OR WMO PUBLICATION 306.
C   3)  IN EITHER CASE, INPUT NUMBERS WILL BE MULTIPLIED BY
C       '10 TO THE NTH' POWER FOUND IN ID(25) OR PDS(27-28),
C       THE D-SCALING FACTOR, PRIOR TO BINARY PACKING.
C   4)  ALL NMC PRODUCED GRIB FIELDS WILL HAVE A GRID DEFINITION
C       SECTION INCLUDED IN THE GRIB MESSAGE.  ID(6) WILL BE
C       SET TO '1'.
C       - GDS WILL BE BUILT BASED ON GRID NUMBER (IGRID), UNLESS
C         IGFLAG=1 (USER SUPPLYING IGDS).  USER MUST STILL SUPPLY
C         IGRID EVEN IF IGDS PROVIDED.
C   5)  IF BIT MAP USED THEN ID(7) OR PDS(8) MUST INDICATE THE
C       PRESENCE OF A BIT MAP.
C   6)  ARRAY KBUF SHOULD BE EQUIVALENCED TO AN INTEGER VALUE OR
C       ARRAY TO MAKE SURE IT IS ON A WORD BOUNDARY.
C   7)  SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C $$
C
      REAL            FLD(*)
C
      INTEGER         IBDSFL(*)
      INTEGER         IBMAP(*)
      INTEGER         ID(*)
      INTEGER         IFLD(*)
      INTEGER         IGDS(*)
      INTEGER         IB(4)
      INTEGER         NLEFT, NUMBMS
C
      CHARACTER * 1   BDS11(11)
      CHARACTER * 1   KBUF(*)
      CHARACTER * 1   PDS(*)
      CHARACTER * 1   GDS(200)
      CHARACTER(1),ALLOCATABLE:: BMS(:)
      CHARACTER(1),ALLOCATABLE:: PFLD(:)
      CHARACTER(1),ALLOCATABLE:: IPFLD(:)
      CHARACTER * 1   SEVEN
      CHARACTER * 1   ZERO
C
C
C   ASCII REP OF  /'G', 'R', 'I', 'B'/
C
      DATA  IB    / 71,  82,  73,  66/
C
      IER    = 0
      IBERR  = 0
      JERR   = 0
      IGRIBL = 8
      IPDSL  = 0
      LENGDS = 0
      LENBMS = 0
      LENBDS = 0
      ITOSS  = 0
C
C            1.0   PRODUCT DEFINITION SECTION(PDS).
C
C   SET ID(6) TO 1 ...OR... MODIFY PDS(8) ...
C      REGARDLESS OF USER SPECIFICATION...
C   NMC GRIB FIELDS WILL ALWAYS HAVE A GDS
C
      IF (IPFLAG .EQ.0) THEN
        ID(6) = 1
        CALL W3FI68(ID,PDS)
      ELSE IF (IPFLAG .EQ. 1) THEN
        IF (IAND(mov_a2i(PDS(8)),64) .EQ. 64) THEN
C         BOTH GDS AND BMS
          PDS(8) = CHAR(192)
        ELSE IF (mov_a2i(PDS(8)) .EQ. 0) THEN
C         GDS ONLY
          PDS(8) = CHAR(128)
        END IF
        CONTINUE
      ELSE
C       PRINT *,' W3FI72 ERROR, IPFLAG IS NOT 0 OR 1 IPFLAG = ',IPFLAG
        JERR = 1
        GO TO 900
      END IF
C
C     GET LENGTH OF PDS
C
      IPDSL = mov_a2i(PDS(1)) * 65536 + mov_a2i(PDS(2)) * 256 +
     &        mov_a2i(PDS(3))
C
C            2.0   GRID DEFINITION SECTION (GDS).
C
C     IF IGFLAG=1 THEN USER IS SUPPLYING THE IGDS INFORMATION
C
      IF (IGFLAG .EQ. 0) THEN
        CALL W3FI71(IGRID,IGDS,IGERR)
        IF (IGERR .EQ. 1) THEN
C         PRINT *,' W3FI71 ERROR, GRID TYPE NOT DEFINED...',IGRID
          JERR = 4
          GO TO 900
        END IF
      END IF
      IF (IGFLAG .EQ. 0  .OR.  IGFLAG .EQ.1) THEN
        CALL W3FI74(IGDS,ICOMP,GDS,LENGDS,NPTS,IGERR)
        IF (IGERR .EQ. 1) THEN
C         PRINT *,' W3FI74 ERROR, GRID REP TYPE NOT VALID...',IGDS(3)
          JERR = 5
          GO TO 900
        ELSE
        END IF
      ELSE
C       PRINT *,' W3FI72 ERROR, IGFLAG IS NOT 0 OR 1 IGFLAG = ',IGFLAG
        JERR = 2
        GO TO 900
      END IF
C
C            3.0   BIT MAP SECTION (BMS).
C
C     SET ITOSS=1 IF BITMAP BEING USED.  W3FI75 WILL TOSS DATA
C     PRIOR TO PACKING.  LATER CODING WILL BE NEEDED WHEN THE
C     'PREDEFINED' GRIDS ARE FINALLY 'DEFINED'.
C
      IF (mov_a2i(PDS(8)) .EQ. 64 .OR.
     &    mov_a2i(PDS(8)) .EQ. 192)   THEN
        ITOSS = 1
        IF (IBFLAG .EQ. 0) THEN
          IF (IBLEN .NE. NPTS) THEN
C           PRINT *,' W3FI72 ERROR, IBLEN .NE. NPTS = ',IBLEN,NPTS
            JERR = 7
            GO TO 900
          END IF
          IF (MOD(IBLEN,16).NE.0) THEN
             NLEFT  = 16 - MOD(IBLEN,16)
          ELSE
             NLEFT  = 0
          END IF
          NUMBMS = 6 + (IBLEN+NLEFT) / 8
          ALLOCATE(BMS(NUMBMS))
          ZERO = CHAR(00)
          BMS = ZERO
          CALL W3FI73(IBFLAG,IBMAP,IBLEN,BMS,LENBMS,IER)
          IF (IER .NE. 0) THEN
C           PRINT *,' W3FI73 ERROR, IBMAP VALUES ARE ALL ZERO'
            JERR = 8
            GO TO 900
          END IF
        ELSE
C         PRINT *,'   BIT MAP PREDEFINED BY CENTER, IBFLAG = ',IBFLAG
        END IF
      END IF
C
C            4.0   BINARY DATA SECTION (BDS).
C
C            4.1   SCALE THE DATA WITH D-SCALE FROM PDS(27-28)
C
      JSCALE = mov_a2i(PDS(27)) * 256 + mov_a2i(PDS(28))
      IF (IAND(JSCALE,32768).NE.0) THEN
        JSCALE = - IAND(JSCALE,32767)
      END IF
      SCALE  = 10.0 ** JSCALE
      IF (ITYPE .EQ. 0) THEN
        DO 410 I = 1,NPTS
          FLD(I) = FLD(I) * SCALE
  410   CONTINUE
      ELSE
        DO 411 I = 1,NPTS
          IFLD(I) = NINT(FLOAT(IFLD(I)) * SCALE)
  411   CONTINUE
      END IF
C
C            4.2   CALL W3FI75 TO PACK DATA AND MAKE BDS.
C
      ALLOCATE(PFLD(NPTS*4))
C
      IF(IBDSFL(2).NE.0) THEN
        ALLOCATE(IPFLD(NPTS*4))
        IPFLD=char(0)
      ELSE
        ALLOCATE(IPFLD(1))
      ENDIF
C
      CALL W3FI75(IBITL,ITYPE,ITOSS,FLD,IFLD,IBMAP,IBDSFL,
     &         NPTS,BDS11,IPFLD,PFLD,LEN,LENBDS,IBERR,PDS,IGDS)
C
      IF(IBDSFL(2).NE.0) THEN
C        CALL XMOVEX(PFLD,IPFLD,NPTS*4)
         do ii = 1, NPTS*4
            PFLD(ii) = IPFLD(ii)
         enddo
      ENDIF
        DEALLOCATE(IPFLD)
C
        IF (IBERR .EQ. 1) THEN
          JERR = 3
          GO TO 900
        END IF
C            4.3   IF D-SCALE NOT 0, RESCALE INPUT FIELD TO
C                   ORIGINAL VALUE
C
      IF (JSCALE.NE.0) THEN
          DSCALE = 1.0 / SCALE
          IF (ITYPE.EQ.0) THEN
              DO 412 I = 1, NPTS
                  FLD(I)  = FLD(I) * DSCALE
  412         CONTINUE
          ELSE
              DO 413 I = 1, NPTS
                  FLD(I)  = NINT(FLOAT(IFLD(I)) * DSCALE)
  413         CONTINUE
          END IF
      END IF
C
C            5.0   OUTPUT SECTION.
C
C            5.1   ZERO OUT THE OUTPUT ARRAY KBUF.
C
      ZERO    = CHAR(00)
      ITOT    = IGRIBL + IPDSL + LENGDS + LENBMS + LENBDS + 4
C     PRINT *,'IGRIBL  =',IGRIBL
C     PRINT *,'IPDSL   =',IPDSL
C     PRINT *,'LENGDS  =',LENGDS
C     PRINT *,'LENBMS  =',LENBMS
C     PRINT *,'LENBDS  =',LENBDS
C     PRINT *,'ITOT    =',ITOT
      KBUF(1:ITOT)=ZERO
C
C            5.2   MOVE SECTION 0 - 'IS' INTO KBUF (8 BYTES).
C
      ISTART  = 0
      DO 520 I = 1,4
        KBUF(I) = CHAR(IB(I))
  520 CONTINUE
C
      KBUF(5) = CHAR(MOD(ITOT / 65536,256))
      KBUF(6) = CHAR(MOD(ITOT /   256,256))
      KBUF(7) = CHAR(MOD(ITOT        ,256))
      KBUF(8) = CHAR(1)
C
C            5.3   MOVE SECTION 1 - 'PDS' INTO KBUF (28 BYTES).
C
      ISTART  = ISTART + IGRIBL
      IF (IPDSL.GT.0) THEN
C        CALL XMOVEX(KBUF(ISTART+1),PDS,IPDSL)
         do ii = 1, IPDSL
            KBUF(ISTART+ii) = PDS(ii)
         enddo
      ELSE
C       PRINT *,'LENGTH OF PDS LESS OR EQUAL 0, IPDSL = ',IPDSL
      END IF
C
C            5.4   MOVE SECTION 2 - 'GDS' INTO KBUF.
C
      ISTART  = ISTART + IPDSL
      IF (LENGDS .GT. 0) THEN
C        CALL XMOVEX(KBUF(ISTART+1),GDS,LENGDS)
         do ii = 1, LENGDS
            KBUF(ISTART+ii) = GDS(ii)
         enddo
      END IF
C
C            5.5   MOVE SECTION 3 - 'BMS' INTO KBUF.
C
      ISTART  = ISTART + LENGDS
      IF (LENBMS .GT. 0) THEN
C        CALL XMOVEX(KBUF(ISTART+1),BMS,LENBMS)
         do ii = 1, LENBMS
            KBUF(ISTART+ii) = BMS(ii)
         enddo
      END IF
C
C            5.6   MOVE SECTION 4 - 'BDS' INTO KBUF.
C
C                  MOVE THE FIRST 11 OCTETS OF THE BDS INTO KBUF.
C
      ISTART  = ISTART + LENBMS
C      CALL XMOVEX(KBUF(ISTART+1),BDS11,11)
      do ii = 1, 11
         KBUF(ISTART+ii) = BDS11(ii)
      enddo
C
C                  MOVE THE PACKED DATA INTO THE KBUF
C
      ISTART  = ISTART + 11
      IF (LEN.GT.0) THEN
C        CALL XMOVEX(KBUF(ISTART+1),PFLD,LEN)
         do ii = 1, LEN
            KBUF(ISTART+ii) = PFLD(ii)
         enddo
      END IF
C
C                  ADD '7777' TO END OFF KBUF
C   NOTE THAT THESE 4 OCTETS NOT INCLUDED IN ACTUAL SIZE OF BDS.
C
      SEVEN  = CHAR(55)
      ISTART = ITOT - 4
      DO 562 I = 1,4
        KBUF(ISTART+I) = SEVEN
 562  CONTINUE
C
 900  CONTINUE
      IF(ALLOCATED(BMS)) DEALLOCATE(BMS)
      IF(ALLOCATED(PFLD)) DEALLOCATE(PFLD)
      RETURN
      END

ext: f

max_stars_repo_name: martinremy/wps

max_issues_repo_name: martinremy/wps

max_forks_repo_name: WRF-CMake/WPS


text:
 Note: Infernal Affairs Commission is simply a reference to the Internal Affairs Commission IAC, when said Commission members vote in a manner that others find galling.  The name is obviously a humorous play on Internal Affairs Commission.

People applying to Internal Affairs Commission IAC: be advised that members of the student body may take issue with your voting decisions.

{{{
Dear Internal Affairs Commissioners,

I was shocked to learn that you voted down the Pledge of Allegiance bill.  I
read the bill, it is wellwritten and perfectly fair.  I cannot understand
how all of you, in your right minds, could vote against a bill to say the
Pledge of Allegiance to the U.S. flag.

Your no vote just doesnt make any sense.  Especially when considering
that UC Davis, as a public research university, relies on funding from the
STATE and FEDERAL levels of government.
Your decision rings with ingratitude for the deluge of money that the
state and federal governments rain down on this school.

Every time the ASUCD Senate has their meeting, many different opinions
bounce around the room.  Isnt it wonderful that the U.S. Constitution
allows for that freedom of speech, without any fear of retribution from the
government?  Voting no on this bill is a disservice to the principles that
even allow ASUCD to exist.

Any and every reputable political body in the entire nation says the Pledge
of Allegiance before their meetings begin.  You are all incredibly immature
for voting down this bill, especially since reciting the Pledge, or even
standing up to recognize the flag, would be optional.

And for those of you who are worried about selfpreservation...
Congratulations on having this in your permanent political record.  I hope
you realize that your chance for ever getting elected to any public office
is ruined, for the rest of your lives.

No discerning citizen would vote for a man/woman to represent them, if that
man/woman is on RECORD saying that the Pledge of Allegiance is not good
enough to be said, aloud, in a public university.  Especially since the
recitation of the Pledge, in this and in every situation, is optional.

Im not trying to be meanspirited about this particular point, but lets
face it.  If you ever run for any public office, your chances of victory are
very slim if the media somehow discovers your no vote.

You have also slimmed your chances of getting hired by any
nongovernmental organization or business that would take issue with your
vote.  All companies do background checks, as Im sure you are well aware.

So I wish you the best of luck in your future ventures, because with this
embarassing affair on record, youll need it.

Shame on you, those who voted no.  Youre not even elected officials.  Let
the Senators, elected by popular vote, do the legislating in ASUCD,
especially when it comes to bills of such incredible importance as this one.

When it comes down to principles, we as people dont bleed blue and gold, or
red and yellow.  Like it or not, we are granted many freedoms by the U.S.
Constitution, and saying the Pledge before a public meeting is simply a way
of recognizing those freedoms.  Thats all.

Sincerely,
Cameron Menezes  UC Davis class of 2005
}}}

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Please see ASUCD Senate Bill #84, ASUCD Information Systems webpage (http://asucd.ucdavis.edu/government/govdocuments/senatebills/20052006)
Please see ASUCD Senate Bill #12, ASUCD Information Systems webpage (http://asucd.ucdavis.edu/government/govdocuments/senatebills/20062007)

Internal Affairs Commissioners voting to have a voluntary recitation of the Pledge of Allegiance at Senate meetings:
Joseph Bleckman  UC Davis class of 2009  (Senate Bill#12, fall 2006)
Tom Woelfel  UC Davis class of 2006        (Senate Bill#84, spring 2006)

Internal Affairs Commissioners voting NOT to have a voluntary recitation of the Pledge of Allegiance at Senate meetings:
Derick Lennox  UC Davis class of 2009  (Senate Bill#12, fall 2006 & Senate Bill#84, spring 2006)
Mike Rivera  UC Davis class of 2008     (Senate Bill#12, fall 2006 & Senate Bill#84, spring 2006)
Matt Schock  UC Davis class of 2006   (Senate Bill#12, fall 2006 & Senate Bill#84, spring 2006)
Users/MatthewShannon Matthew Shannon  UC Davis class of 2008  (Senate Bill#84, spring 2006)
Users/JillWeinstein Jill Weinstein  UC Davis class of 2007   (Senate Bill#12, fall 2006 & Senate Bill#84, spring 2006)
Users/MaxMikalonis Max Mikalonis  UC Davis class of 2009  (Senate Bill#12, fall 2006)
Alix Giroux  UC Davis class of 2009    (Senate Bill#12, fall 2006)
Users/PaulHarms Paul Harms  UC Davis class of 2009    (Senate Bill#84, spring 2006)

Internal Affairs Commissioners ABSTAINING to vote on whether to have a voluntary recitation of the Pledge of Allegiance at Senate meetings:
Jack Shannon  UC Davis class of 2008  (Senate Bill#84, spring 2006)
Eric Wu      UC Davis class of 2006  (Senate Bill#84, spring 2006)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module map_non_coincident_grids_mod

use precision_mod
use coords_mod
use field_section_mod
use pointer_mod
use subfield_mod

implicit none

type, public :: bounds
  real(kind=double_precision) :: west_extreme_lon,east_extreme_lon
  real(kind=double_precision) :: north_extreme_lat,south_extreme_lat
end type

type, public,abstract :: vertex_coords
end type

type, public, extends(vertex_coords) :: unstructured_grid_vertex_coords
  real(kind=double_precision), dimension(:), pointer :: vertex_lats
  real(kind=double_precision), dimension(:), pointer :: vertex_lons
end type

interface unstructured_grid_vertex_coords
  procedure :: unstructured_grid_vertex_coords_constructor
end interface unstructured_grid_vertex_coords

type, public, abstract :: non_coincident_grid_mapper
  class(field_section), pointer  :: mask
  class(field_section), pointer  :: cell_numbers
  class(field_section), pointer  :: pixel_center_lats
  class(field_section), pointer  :: pixel_center_lons
  class(field_section), pointer  :: cell_vertex_coords
  class(subfield), pointer  :: area_to_consider_mask
  class(subfield), pointer  :: secondary_area_to_consider_mask
  class(coords), pointer :: coarse_cell_coords
  class(bounds), pointer :: cell_bounds
  class(section_coords),pointer :: fine_grid_shape
  integer, dimension(:), pointer :: section_min_lats
  integer, dimension(:), pointer :: section_min_lons
  integer, dimension(:), pointer :: section_max_lats
  integer, dimension(:), pointer :: section_max_lons
  logical :: display_progress  = .true.
  contains
    procedure :: set_cell_numbers
    procedure :: generate_pixels_in_cell_mask
    procedure :: check_if_pixel_is_in_cell
    procedure :: generate_cell_numbers
    procedure :: process_cell
    procedure :: generate_limits
    procedure :: offset_limits
    procedure(process_pixel_for_limits), deferred :: process_pixel_for_limits
    procedure(generate_cell_bounds), deferred :: generate_cell_bounds
    procedure(generate_areas_to_consider), deferred :: generate_areas_to_consider
    procedure(check_if_pixel_center_is_in_bounds), deferred :: &
      check_if_pixel_center_is_in_bounds
    procedure(create_new_mask), deferred :: create_new_mask
    procedure(assign_cell_numbers), deferred :: assign_cell_numbers
    procedure(print_progress), deferred :: print_progress
end type

abstract interface
  subroutine generate_cell_bounds(this)
    import non_coincident_grid_mapper
    implicit none
    class(non_coincident_grid_mapper), intent(inout) :: this
  end subroutine generate_cell_bounds

  subroutine generate_areas_to_consider(this)
    import non_coincident_grid_mapper
    implicit none
    class(non_coincident_grid_mapper), intent(inout) :: this
  end subroutine generate_areas_to_consider

  function check_if_pixel_center_is_in_bounds(this,pixel_center_lat,pixel_center_lon) &
      result(is_in_bounds)
    import non_coincident_grid_mapper
    import double_precision
    implicit none
    class(non_coincident_grid_mapper), intent(inout) :: this
    real(kind=double_precision), intent(inout) :: pixel_center_lat, pixel_center_lon
    logical :: is_in_bounds
  end function check_if_pixel_center_is_in_bounds

  subroutine create_new_mask(this)
    import non_coincident_grid_mapper
    implicit none
    class(non_coincident_grid_mapper), intent(inout) :: this
  end subroutine create_new_mask

  subroutine assign_cell_numbers(this)
    import non_coincident_grid_mapper
    implicit none
    class(non_coincident_grid_mapper), intent(inout) :: this
  end subroutine

  subroutine process_pixel_for_limits(this,coords_in)
    import non_coincident_grid_mapper
    import coords
    implicit none
    class(non_coincident_grid_mapper), intent(inout) :: this
    class(coords), pointer, intent(inout) :: coords_in
  end subroutine

  subroutine print_progress(this)
    import non_coincident_grid_mapper
    implicit none
    class(non_coincident_grid_mapper), intent(in) :: this
  end subroutine print_progress
end interface

type, extends(non_coincident_grid_mapper), public :: &
  icon_icosohedral_cell_latlon_pixel_ncg_mapper
  logical :: longitudal_range_centered_on_zero
contains
  procedure :: icon_icosohedral_cell_get_vertex_coords
  procedure :: latlon_pixel_generate_area_to_consider
  procedure, nopass :: icon_icosohedral_cell_calculate_line
  procedure :: generate_cell_bounds => icon_icosohedral_cell_generate_cell_bounds
  procedure :: generate_areas_to_consider => latlon_pixel_generate_areas_to_consider
  procedure :: check_if_pixel_center_is_in_bounds => &
    icon_icosohedral_cell_check_if_pixel_center_is_in_bounds
  procedure :: create_new_mask => latlon_pixel_create_new_mask
  procedure :: assign_cell_numbers => &
    icon_icosohedral_cell_latlon_pixel_assign_cell_numbers
  procedure :: process_pixel_for_limits => latlon_process_pixel_for_limits
  procedure :: init_icon_icosohedral_cell_latlon_pixel_ncg_mapper
  procedure :: icon_icosohedral_cell_latlon_pixel_ncg_mapper_destructor
  procedure :: print_progress => icon_icosohedral_cell_print_progress
end type

interface icon_icosohedral_cell_latlon_pixel_ncg_mapper
    procedure icon_icosohedral_cell_latlon_pixel_ncg_mapper_constructor
end interface icon_icosohedral_cell_latlon_pixel_ncg_mapper

contains

  function unstructured_grid_vertex_coords_constructor(vertex_lats_in,&
                                                         vertex_lons_in) &
                                                         result(constructor)

    real(kind=double_precision), dimension(:), target, intent(in) :: &
        vertex_lats_in
    real(kind=double_precision), dimension(:), target, intent(in) :: &
        vertex_lons_in
    type(unstructured_grid_vertex_coords) :: constructor
      constructor%vertex_lats => vertex_lats_in
      constructor%vertex_lons => vertex_lons_in
  end function unstructured_grid_vertex_coords_constructor

  subroutine generate_pixels_in_cell_mask(this,cell_coords)
    class(non_coincident_grid_mapper), intent(inout) :: this
    class(coords), pointer, intent(in) :: cell_coords
    !Despite the name this field section actually covers the entire field
      call this%create_new_mask()
      this%coarse_cell_coords => cell_coords
      call this%generate_cell_bounds()
      call this%generate_areas_to_consider()
      call this%area_to_consider_mask%for_all(check_if_pixel_is_in_cell_wrapper,this)
      if (associated(this%secondary_area_to_consider_mask)) then
        call this%secondary_area_to_consider_mask%for_all(check_if_pixel_is_in_cell_wrapper,&
                                                          this)
      end if
      if (this%display_progress) call this%print_progress()
  end subroutine generate_pixels_in_cell_mask

  subroutine check_if_pixel_is_in_cell_wrapper(this,coords_in)
    class(*), intent(inout) :: this
    class(coords), pointer,intent(inout) :: coords_in
      select type (this)
      class is (non_coincident_grid_mapper)
        call this%check_if_pixel_is_in_cell(coords_in)
      end select
  end subroutine check_if_pixel_is_in_cell_wrapper

  subroutine check_if_pixel_is_in_cell(this,coords_in)
    class(non_coincident_grid_mapper), intent(inout) :: this
    class(coords), pointer,intent(inout) :: coords_in
    class(*), pointer :: pixel_center_lat_ptr, pixel_center_lon_ptr
    real(kind=double_precision) :: pixel_center_lat, pixel_center_lon
      pixel_center_lat_ptr => this%pixel_center_lats%get_value(coords_in)
      select type (pixel_center_lat_ptr)
        type is (real(kind=double_precision))
          pixel_center_lat = pixel_center_lat_ptr
      end select
      deallocate(pixel_center_lat_ptr)
      pixel_center_lon_ptr => this%pixel_center_lons%get_value(coords_in)
      select type (pixel_center_lon_ptr)
        type is (real(kind=double_precision))
          pixel_center_lon = pixel_center_lon_ptr
      end select
      deallocate(pixel_center_lon_ptr)
      call this%mask%set_value(coords_in,&
                               this%check_if_pixel_center_is_in_bounds(pixel_center_lat,&
                                                                       pixel_center_lon))
      deallocate(coords_in)
  end subroutine check_if_pixel_is_in_cell

  function generate_cell_numbers(this) &
      result(cell_numbers)
    class(non_coincident_grid_mapper), intent(inout) :: this
    class(field_section), pointer  :: cell_numbers
      call this%cell_vertex_coords%for_all_section(process_cell_wrapper,this)
      cell_numbers =>this%cell_numbers
  end function generate_cell_numbers

  subroutine set_cell_numbers(this,cell_numbers_in)
    class(non_coincident_grid_mapper), intent(inout) :: this
    class(field_section), pointer  :: cell_numbers_in
      this%cell_numbers => cell_numbers_in
  end subroutine set_cell_numbers

  subroutine process_cell_wrapper(this,coords_in)
    class(*), intent(inout) :: this
    class(coords), pointer,intent(inout) :: coords_in
      select type (this)
        class is (non_coincident_grid_mapper)
          call this%process_cell(coords_in)
      end select
  end subroutine process_cell_wrapper

  subroutine process_cell(this,coords_in)
    class(non_coincident_grid_mapper), intent(inout) :: this
    class(coords), pointer,intent(inout) :: coords_in
      call this%generate_pixels_in_cell_mask(coords_in)
      call this%assign_cell_numbers()
      call this%mask%deallocate_data()
      deallocate(this%mask)
      deallocate(this%cell_bounds)
      call this%area_to_consider_mask%destructor()
      deallocate(this%area_to_consider_mask)
      if (associated(this%secondary_area_to_consider_mask)) then
      call this%secondary_area_to_consider_mask%destructor()
       deallocate(this%secondary_area_to_consider_mask)
      end if
      deallocate(coords_in)
  end subroutine process_cell

  subroutine generate_limits(this)
    class(non_coincident_grid_mapper), intent(inout) :: this
      call this%cell_numbers%for_all_section(process_pixel_for_limits_wrapper,this)
  end subroutine generate_limits

  subroutine offset_limits(this,lat_offset)
    class(non_coincident_grid_mapper), intent(inout) :: this
    integer, intent(in) :: lat_offset
      this%section_min_lats(:) = this%section_min_lats(:) + lat_offset
      this%section_max_lats(:) = this%section_max_lats(:) + lat_offset
  end subroutine offset_limits

  subroutine process_pixel_for_limits_wrapper(this,coords_in)
    class(*), intent(inout) :: this
    class(coords), pointer, intent(inout) :: coords_in
      select type(this)
        class is (non_coincident_grid_mapper)
          call this%process_pixel_for_limits(coords_in)
      end select
  end subroutine process_pixel_for_limits_wrapper

  subroutine latlon_process_pixel_for_limits(this,coords_in)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
    class(coords), pointer,intent(inout) :: coords_in
    integer :: i,j
    integer :: working_cell_number
    class(*), pointer :: working_cell_number_ptr
      select type (coords_in)
      type is (latlon_coords)
        i = coords_in%lat
        j = coords_in%lon
      end select
      working_cell_number_ptr => this%cell_numbers%get_value(coords_in)
      select type (working_cell_number_ptr)
      type is (integer)
        working_cell_number = working_cell_number_ptr
      end select
      if (this%section_min_lats(working_cell_number) > i) then
        this%section_min_lats(working_cell_number) = i
      end if
      if (this%section_max_lats(working_cell_number) < i) then
        this%section_max_lats(working_cell_number) = i
      end if
      if (this%section_min_lons(working_cell_number) > j) then
        this%section_min_lons(working_cell_number) = j
      end if
      if (this%section_max_lons(working_cell_number) < j) then
        this%section_max_lons(working_cell_number) = j
      end if
      deallocate(working_cell_number_ptr)
      deallocate(coords_in)
  end subroutine latlon_process_pixel_for_limits

  subroutine latlon_pixel_generate_areas_to_consider(this)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
    real(kind=double_precision) :: rotated_west_extreme_lon, rotated_east_extreme_lon
    real(kind=double_precision) :: temp_for_swapping_lons
    real(kind=double_precision) :: fine_grid_zero_line
      select type (latlon_fine_grid_shape => this%fine_grid_shape)
      type is (latlon_section_coords)
        fine_grid_zero_line = latlon_fine_grid_shape%zero_line
      end select
      rotated_west_extreme_lon = this%cell_bounds%west_extreme_lon
      rotated_east_extreme_lon = this%cell_bounds%east_extreme_lon
      if (rotated_west_extreme_lon < 0  + fine_grid_zero_line) then
        rotated_west_extreme_lon = rotated_west_extreme_lon + 360.0
      end if
      if (rotated_east_extreme_lon < 0 + fine_grid_zero_line) then
        rotated_east_extreme_lon = rotated_east_extreme_lon + 360.0
      end if
      if (rotated_west_extreme_lon > 360.0  + fine_grid_zero_line) then
        rotated_west_extreme_lon = rotated_west_extreme_lon - 360.0
      end if
      if (rotated_east_extreme_lon > 360.0 + fine_grid_zero_line) then
        rotated_east_extreme_lon = rotated_east_extreme_lon - 360.0
      end if
      if (rotated_east_extreme_lon - rotated_west_extreme_lon > 180.0) then
        temp_for_swapping_lons = rotated_west_extreme_lon
        rotated_west_extreme_lon = rotated_east_extreme_lon
        rotated_east_extreme_lon = temp_for_swapping_lons
      end if
      if (rotated_west_extreme_lon < rotated_east_extreme_lon .or. &
          rotated_east_extreme_lon == 0.0 + fine_grid_zero_line) then
        if(rotated_east_extreme_lon == 0.0 + fine_grid_zero_line) then
          rotated_east_extreme_lon = 360.0 + fine_grid_zero_line
        end if
        this%area_to_consider_mask => &
          this%latlon_pixel_generate_area_to_consider(rotated_west_extreme_lon,&
                                                      rotated_east_extreme_lon)
        this%secondary_area_to_consider_mask => null()
      else if (rotated_west_extreme_lon > rotated_east_extreme_lon) then
        this%area_to_consider_mask => &
          this%latlon_pixel_generate_area_to_consider(0.0_double_precision  + fine_grid_zero_line,&
                                                      rotated_east_extreme_lon)
        this%secondary_area_to_consider_mask => &
          this%latlon_pixel_generate_area_to_consider(rotated_west_extreme_lon, &
                                                      360.0_double_precision + fine_grid_zero_line)
      else
        stop "Error - cell appears to have no width"
      end if
  end subroutine  latlon_pixel_generate_areas_to_consider

  function latlon_pixel_generate_area_to_consider(this,area_min_lon, area_max_lon) &
      result(area_to_consider)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
    class(latlon_section_coords), pointer :: area_to_consider_section_coords
    class(*),pointer,dimension(:,:) :: mask
    class(subfield), pointer  :: area_to_consider
    real(kind=double_precision) :: area_min_lon, area_max_lon
    type(latlon_coords),pointer :: pixel_in_row_or_column
    class(*), pointer :: pixel_center_lat_ptr,pixel_center_lon_ptr
    real(kind=double_precision)  :: pixel_center_lat,pixel_center_lon
    real(kind=double_precision) :: working_min_lon, working_max_lon
    integer :: area_min_lat_index, area_max_lat_index
    integer :: area_min_lon_index, area_max_lon_index
    integer :: nlat,nlon
    integer :: i
    logical :: minimum_found
      area_min_lat_index = 0
      area_max_lat_index = 0
      area_min_lon_index = 0
      area_max_lon_index = 0
      working_min_lon = 99999.0
      working_max_lon = -99999.0
      select type (latlon_fine_grid_shape => this%fine_grid_shape)
        type is (latlon_section_coords)
          nlat = latlon_fine_grid_shape%section_width_lat
          nlon = latlon_fine_grid_shape%section_width_lon
      end select
      minimum_found = .FALSE.
      allocate(pixel_in_row_or_column)
      do i=1,nlat
        pixel_in_row_or_column = latlon_coords(i,1)
        pixel_center_lat_ptr => this%pixel_center_lats%get_value(pixel_in_row_or_column)
        select type (pixel_center_lat_ptr)
          type is (real(kind=double_precision))
            pixel_center_lat = pixel_center_lat_ptr
        end select
        deallocate(pixel_center_lat_ptr)
        if ( pixel_center_lat <= this%cell_bounds%north_extreme_lat &
            .and. .not. minimum_found) then
          area_min_lat_index = i
          minimum_found = .TRUE.
        end if
        area_max_lat_index = i
        if (pixel_center_lat < &
            this%cell_bounds%south_extreme_lat) exit
      end do
      minimum_found = .FALSE.
      do i=1,nlon
        pixel_in_row_or_column = latlon_coords(1,i)
        pixel_center_lon_ptr => this%pixel_center_lons%get_value(pixel_in_row_or_column)
        select type (pixel_center_lon_ptr)
          type is (real(kind=double_precision))
            pixel_center_lon = pixel_center_lon_ptr
        end select
        deallocate(pixel_center_lon_ptr)
        if (pixel_center_lon >= area_min_lon .and. &
            pixel_center_lon < working_min_lon) then
          area_min_lon_index = i
          working_min_lon = pixel_center_lon
        end if
        if (pixel_center_lon < area_max_lon .and. &
            pixel_center_lon > working_max_lon) then
          area_max_lon_index = i
          working_max_lon = pixel_center_lon
        end if
      end do
      deallocate(pixel_in_row_or_column)
      allocate(area_to_consider_section_coords,&
               source=latlon_section_coords(area_min_lat_index, &
                                            area_min_lon_index, &
                                            area_max_lat_index - area_min_lat_index + 1, &
                                            area_max_lon_index - area_min_lon_index + 1))
      allocate(logical::mask(1:area_max_lat_index - area_min_lat_index + 1, &
                             1:area_max_lon_index - area_min_lon_index + 1))
      select type (mask)
      type is (logical)
        mask = .FALSE.
      end select
      area_to_consider=> latlon_subfield(mask,area_to_consider_section_coords)
      deallocate(area_to_consider_section_coords)
  end function latlon_pixel_generate_area_to_consider

  subroutine latlon_pixel_create_new_mask(this)
      class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
      class(*),pointer,dimension(:,:) :: new_mask
        select type(fine_grid_shape => this%fine_grid_shape)
        type is (latlon_section_coords)
          allocate(logical::new_mask(1:fine_grid_shape%section_width_lat, &
                                     1:fine_grid_shape%section_width_lon))
          select type(new_mask)
          type is (logical)
            new_mask = .False.
          end select
          this%mask => latlon_field_section(new_mask,fine_grid_shape)
        end select
  end subroutine latlon_pixel_create_new_mask

  subroutine icon_icosohedral_cell_print_progress(this)
      class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(in) :: this
      character*20 :: output
        select type (cell_vertex_coords => this%cell_vertex_coords)
        type is (icon_single_index_field_section)
          select type (coarse_cell_coords => this%coarse_cell_coords)
          type is (generic_1d_coords)
            if(mod(coarse_cell_coords%index,cell_vertex_coords%get_num_points()/10) == 0) then
              write(output,*) 100*coarse_cell_coords%index/cell_vertex_coords%get_num_points()
              write(*,*) trim(output) // " % complete"
            end if
          end select
        end select
  end subroutine icon_icosohedral_cell_print_progress

  subroutine icon_icosohedral_cell_generate_cell_bounds(this)
      class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
      real(kind=double_precision) :: vertex_one_lat, vertex_two_lat, vertex_three_lat
      real(kind=double_precision) :: vertex_one_lon, vertex_two_lon, vertex_three_lon
        vertex_one_lat = this%icon_icosohedral_cell_get_vertex_coords(this%coarse_cell_coords,1,.TRUE.)
        vertex_two_lat = this%icon_icosohedral_cell_get_vertex_coords(this%coarse_cell_coords,2,.TRUE.)
        vertex_three_lat = this%icon_icosohedral_cell_get_vertex_coords(this%coarse_cell_coords,3,.TRUE.)
        vertex_one_lon = this%icon_icosohedral_cell_get_vertex_coords(this%coarse_cell_coords,1,.FALSE.)
        vertex_two_lon = this%icon_icosohedral_cell_get_vertex_coords(this%coarse_cell_coords,2,.FALSE.)
        vertex_three_lon = this%icon_icosohedral_cell_get_vertex_coords(this%coarse_cell_coords,3,.FALSE.)
        if (vertex_one_lat == 90.0 .or. vertex_one_lat == -90.0) then
          vertex_one_lon = vertex_two_lon
        else if (vertex_two_lat == 90.0 .or. vertex_two_lat == -90.0) then
          vertex_two_lon = vertex_three_lon
        else if (vertex_three_lat == 90.0 .or. vertex_three_lat == -90.0) then
          vertex_three_lon = vertex_one_lon
        end if
        if (vertex_two_lon - vertex_one_lon > 180.0) then
          vertex_one_lon = vertex_one_lon + 360.0
        else if (vertex_two_lon - vertex_one_lon < -180.0) then
          vertex_two_lon = vertex_two_lon + 360.0
        end if
        if (vertex_three_lon - vertex_one_lon > 180.0) then
          vertex_one_lon = vertex_one_lon + 360.0
          vertex_two_lon = vertex_two_lon + 360.0
        else if (vertex_three_lon - vertex_one_lon < -180.0) then
          vertex_three_lon = vertex_three_lon + 360.0
        end if
        allocate(this%cell_bounds)
        this%cell_bounds%west_extreme_lon = &
                         min(min(vertex_one_lon,vertex_two_lon),vertex_three_lon)
        this%cell_bounds%east_extreme_lon = &
                         max(max(vertex_one_lon,vertex_two_lon),vertex_three_lon)
        this%cell_bounds%south_extreme_lat = &
                         min(min(vertex_one_lat,vertex_two_lat),vertex_three_lat)
        this%cell_bounds%north_extreme_lat = &
                         max(max(vertex_one_lat,vertex_two_lat),vertex_three_lat)
  end subroutine icon_icosohedral_cell_generate_cell_bounds

  function icon_icosohedral_cell_check_if_pixel_center_is_in_bounds(this,pixel_center_lat,pixel_center_lon) &
        result(is_in_bounds)
      class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
      real(kind=double_precision), intent(inout) :: pixel_center_lat, pixel_center_lon
      real(kind=double_precision) :: vertex_one_lat, vertex_one_lon
      real(kind=double_precision) :: vertex_two_lat, vertex_two_lon
      real(kind=double_precision) :: vertex_three_lat, vertex_three_lon
      real(kind=double_precision) :: vertex_temp_lat, vertex_temp_lon
      real(kind=double_precision) :: lat_difference_at_pixel_center_lon
      real(kind=double_precision) :: lat_difference_at_vertex_three_lon
      logical :: is_in_bounds
      integer :: i
        select type (coarse_cell_coords => this%coarse_cell_coords)
        type is (generic_1d_coords)
          vertex_one_lat = &
            this%icon_icosohedral_cell_get_vertex_coords(coarse_cell_coords,1,.TRUE.)
          vertex_one_lon = &
            this%icon_icosohedral_cell_get_vertex_coords(coarse_cell_coords,1,.FALSE.)
          vertex_two_lat = &
            this%icon_icosohedral_cell_get_vertex_coords(coarse_cell_coords,2,.TRUE.)
          vertex_two_lon = &
            this%icon_icosohedral_cell_get_vertex_coords(coarse_cell_coords,2,.FALSE.)
          vertex_three_lat = &
            this%icon_icosohedral_cell_get_vertex_coords(coarse_cell_coords,3,.TRUE.)
          vertex_three_lon = &
            this%icon_icosohedral_cell_get_vertex_coords(coarse_cell_coords,3,.FALSE.)
        end select
        if (vertex_two_lon - vertex_one_lon > 180.0) then
          vertex_two_lon = vertex_two_lon - 360.0
        else if(vertex_two_lon - vertex_one_lon < -180.0) then
          vertex_two_lon = vertex_two_lon + 360.0
        end if
        if (vertex_three_lon - vertex_one_lon > 180.0) then
          vertex_three_lon = vertex_three_lon - 360.0
        else if(vertex_three_lon - vertex_one_lon < -180.0) then
          vertex_three_lon = vertex_three_lon + 360.0
        end if
        if (pixel_center_lon - vertex_one_lon > 180.0) then
          pixel_center_lon = pixel_center_lon - 360.0
        else if(pixel_center_lon - vertex_one_lon < -180.0) then
          pixel_center_lon = pixel_center_lon + 360.0
        end if
        is_in_bounds = .TRUE.
        do i = 1,3
          vertex_temp_lon = vertex_three_lon
          vertex_three_lon = vertex_two_lon
          vertex_two_lon = vertex_one_lon
          vertex_one_lon = vertex_temp_lon
          vertex_temp_lat = vertex_three_lat
          vertex_three_lat = vertex_two_lat
          vertex_two_lat = vertex_one_lat
          vertex_one_lat = vertex_temp_lat
          if (vertex_one_lon-vertex_two_lon == 0 .or. &
              vertex_two_lat == 90.0 .or. vertex_two_lat == -90.0) then
            is_in_bounds = is_in_bounds .and. &
                           SIGN(1.0_double_precision,pixel_center_lon - vertex_one_lon) == &
                           SIGN(1.0_double_precision,vertex_three_lon - vertex_one_lon)
          else if (vertex_two_lat - vertex_one_lat == 0  ) then
            is_in_bounds = is_in_bounds .and. &
                           SIGN(1.0_double_precision,pixel_center_lat - vertex_one_lat) == &
                           SIGN(1.0_double_precision,vertex_three_lat - vertex_one_lat)
          else if(vertex_one_lat == 90.0 .or. vertex_one_lat == -90.0) then
            is_in_bounds = is_in_bounds .and. &
                           SIGN(1.0_double_precision,pixel_center_lon - vertex_two_lon) == &
                           SIGN(1.0_double_precision,vertex_three_lon - vertex_two_lon)
          else
            lat_difference_at_pixel_center_lon = pixel_center_lat - &
              this%icon_icosohedral_cell_calculate_line(pixel_center_lon, vertex_one_lon, &
                                                        vertex_two_lon, vertex_one_lat, &
                                                        vertex_two_lat)
            lat_difference_at_vertex_three_lon = vertex_three_lat - &
              this%icon_icosohedral_cell_calculate_line(vertex_three_lon, vertex_one_lon, &
                                                        vertex_two_lon, vertex_one_lat, &
                                                        vertex_two_lat)
            if (lat_difference_at_pixel_center_lon == 0) then
              is_in_bounds = is_in_bounds .and. &
                             (lat_difference_at_vertex_three_lon >= 0)
            else
              is_in_bounds = is_in_bounds .and. &
                              SIGN(1.0_double_precision, &
                                   lat_difference_at_pixel_center_lon) == &
                              SIGN(1.0_double_precision, &
                                   lat_difference_at_vertex_three_lon)
            end if
          end if
        end do
  end function icon_icosohedral_cell_check_if_pixel_center_is_in_bounds

  function icon_icosohedral_cell_calculate_line(x,x1,x2,y1,y2) result(y)
    real(kind=double_precision), intent(in) :: x,x1,x2,y1,y2
    real(kind=double_precision) :: y
      y = ((y2 - y1)/(x2 - x1))*(x-x1) + y1
  end function icon_icosohedral_cell_calculate_line

  function icon_icosohedral_cell_get_vertex_coords(this,coords_in,vertex_num_in,return_lat) result (coordinate)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(in) :: this
    class(coords), intent(in) ::  coords_in
    integer, intent(in) :: vertex_num_in
    logical, intent(in) :: return_lat
    real(kind=double_precision) :: coordinate
    real(kind=double_precision) :: lat_coordinate
    class(*), pointer :: vertex_positions_ptr
    integer :: vertex_num
      vertex_num = vertex_num_in
      vertex_positions_ptr => this%cell_vertex_coords%get_value(coords_in)
      select type (vertex_positions_ptr)
       type is (unstructured_grid_vertex_coords)
       lat_coordinate = vertex_positions_ptr%vertex_lats(vertex_num)
       if (return_lat) then
          coordinate = lat_coordinate
       else
          !If at pole return one of the other coordinates instead
          if (lat_coordinate == 90.0 .or. lat_coordinate == -90.0) then
            if (vertex_num == 1) then
              vertex_num = 2
            else
              vertex_num = 1
            end if
          end if
          coordinate = vertex_positions_ptr%vertex_lons(vertex_num)
          if (this%longitudal_range_centered_on_zero) then
            if (coordinate < 0.0) then
              coordinate = coordinate + 360.0
            end if
          end if
       end if
     end select
     deallocate(vertex_positions_ptr)
  end function icon_icosohedral_cell_get_vertex_coords

  subroutine icon_icosohedral_cell_latlon_pixel_assign_cell_numbers(this)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
      select type (mask=>this%mask)
      type is (latlon_field_section)
        select type (data=>mask%data)
        type is (logical)
          select type(cell_numbers => this%cell_numbers)
          type is (latlon_field_section)
            select type (cell_numbers_data => cell_numbers%data)
            type is (integer)
              select type (coarse_cell_coords=>this%coarse_cell_coords)
              type is (generic_1d_coords)
                  where(data) cell_numbers_data = coarse_cell_coords%index
              end select
            end select
          end select
        end select
      end select
  end subroutine icon_icosohedral_cell_latlon_pixel_assign_cell_numbers

  subroutine init_icon_icosohedral_cell_latlon_pixel_ncg_mapper(this,&
                                                                pixel_center_lats,&
                                                                pixel_center_lons,&
                                                                cell_vertex_coords,&
                                                                fine_grid_shape, &
                                                                longitudal_range_centered_on_zero_in)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
    class(latlon_field_section), pointer, intent(in)  :: pixel_center_lats
    class(latlon_field_section), pointer, intent(in)  :: pixel_center_lons
    class(icon_single_index_field_section), pointer, intent(inout)  :: cell_vertex_coords
    class(latlon_section_coords),pointer,intent(in) :: fine_grid_shape
    class(*),pointer,dimension(:,:) :: new_cell_numbers
    logical, optional :: longitudal_range_centered_on_zero_in
      if (present(longitudal_range_centered_on_zero_in)) then
        this%longitudal_range_centered_on_zero = longitudal_range_centered_on_zero_in
      else
        this%longitudal_range_centered_on_zero = .false.
      end if
      this%pixel_center_lats => pixel_center_lats
      this%pixel_center_lons => pixel_center_lons
      this%cell_vertex_coords => cell_vertex_coords
      this%fine_grid_shape => fine_grid_shape
      allocate(integer::new_cell_numbers(1:fine_grid_shape%section_width_lat, &
                                         1:fine_grid_shape%section_width_lon))
      allocate(this%section_min_lats(cell_vertex_coords%get_num_points()))
      allocate(this%section_min_lons(cell_vertex_coords%get_num_points()))
      allocate(this%section_max_lats(cell_vertex_coords%get_num_points()))
      allocate(this%section_max_lons(cell_vertex_coords%get_num_points()))
      this%section_min_lats(:) = fine_grid_shape%section_width_lat + 1
      this%section_min_lons(:) = fine_grid_shape%section_width_lon + 1
      this%section_max_lats(:) = 0
      this%section_max_lons(:) = 0
      select type(new_cell_numbers)
      type is (integer)
        new_cell_numbers = 0
      end select
      this%cell_numbers => latlon_field_section(new_cell_numbers,fine_grid_shape)
  end subroutine init_icon_icosohedral_cell_latlon_pixel_ncg_mapper

  subroutine icon_icosohedral_cell_latlon_pixel_ncg_mapper_destructor(this)
    class(icon_icosohedral_cell_latlon_pixel_ncg_mapper), intent(inout) :: this
      deallocate(this%section_min_lats)
      deallocate(this%section_min_lons)
      deallocate(this%section_max_lats)
      deallocate(this%section_max_lons)
      call this%cell_numbers%deallocate_data()
      deallocate(this%cell_numbers)
  end subroutine icon_icosohedral_cell_latlon_pixel_ncg_mapper_destructor

  function icon_icosohedral_cell_latlon_pixel_ncg_mapper_constructor(pixel_center_lats,&
                                                                     pixel_center_lons,&
                                                                     cell_vertex_coords,&
                                                                     fine_grid_shape, &
                                                                     longitudal_range_centered_on_zero_in) &
                                                                     result(constructor)
        type(icon_icosohedral_cell_latlon_pixel_ncg_mapper), allocatable :: constructor
        class(latlon_field_section), pointer, intent(in)  :: pixel_center_lats
        class(latlon_field_section), pointer, intent(in)  :: pixel_center_lons
        class(icon_single_index_field_section), pointer, intent(inout)  :: cell_vertex_coords
        class(latlon_section_coords),pointer, intent(in) :: fine_grid_shape
        logical, optional :: longitudal_range_centered_on_zero_in
        logical :: longitudal_range_centered_on_zero
            if (present(longitudal_range_centered_on_zero_in)) then
              longitudal_range_centered_on_zero = longitudal_range_centered_on_zero_in
            else
              longitudal_range_centered_on_zero = .false.
            end if
            allocate(constructor)
            call constructor%init_icon_icosohedral_cell_latlon_pixel_ncg_mapper(pixel_center_lats,&
                                                                                pixel_center_lons,&
                                                                                cell_vertex_coords,&
                                                                                fine_grid_shape,&
                                                                                longitudal_range_centered_on_zero)
  end function icon_icosohedral_cell_latlon_pixel_ncg_mapper_constructor

end module map_non_coincident_grids_mod

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD