x2=xG(i-offset+1,j-offset,bi,bj)

                    if ((x2-x1).gt.180) x2=x2-360

                    if ((x1-x2).gt.180) x2=x2+360

                    dx=0.5*(x3+x4-x1-x2)

                    dx=0.5*(x2+x4-x1-x3)

                    x2=xG(i-offset+1,j-offset,bi,bj)

                    if ((x2-x1).gt.180) x2=x2-360

                    if ((x1-x2).gt.180) x2=x2+360

                    dx=0.5*(x3+x4-x1-x2)

                    dx=0.5*(x2+x4-x1-x3)

      _RL     fac, x1, x2, x3, x4, y1, y2, y3, y4, dx, dy

            x2=x*x

            x3=x2*x

            a = x3 + k1plocal*x2 + k12p*x + k123p

            da = 3.0*x2 + 2.0*k1plocal*x + k12p

            b = x2 + k1local*x + k12

     &        kwlocal/x2 +

     &        stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +

     &        ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -

     &        pt*x2*(3.0*a-x*da)/a2

              x2=x*x

              x3=x2*x

              a = x3 + k1plocal*x2 + k12p*x + k123p

              da = 3.0*x2 + 2.0*k1plocal*x + k12p

              b = x2 + k1local*x + k12

        _RL  x2

        _RL  x2

        _RL  x2

      _RL x1,x2,x3,we,pe,fw,fp,pa,pb,pc,ax,ba,bb,t1,ca,cb,t2

        x2=(spre(i,j,k2)-spre(i,j,k1))/x1

        pe=(log10(x2)-p1)/dpe

      _RL aa1,aa2,aa3,aa4,aa5,pp,x2,x3,x4,x5

      _RL  xx,clog1,wlog,dc,dw,x1,x2,y2

          x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw

          y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc

      real*8 x1, x2, xs, t           

         x2=2.0*port_rand(seed)-1.0

         xs=x1**2+x2**2

      _RL  xx,clog1,wlog,dc,dw,x1,x2,y2

          x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw

          y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc