xx=global2local_i(xpart(ip,bi,bj),bi,bj,mythid)

               call FLT_BILINEAR  (xx,yy,uu,kp,uVel,  2,bi,bj)

               call FLT_BILINEAR  (xx,yy,vv,kp,vVel,  3,bi,bj)

               call FLT_BILINEAR2D(xx,yy,pp,   cg2d_x,1,bi,bj)

               call FLT_BILINEAR  (xx,yy,tt,kp,theta, 1,bi,bj)

               call FLT_BILINEAR  (xx,yy,ss,kp,salt,  1,bi,bj)

      _RL xx, yy, xlo, xhi, ylo, yhi

      INTEGER bi, bj, xx, yy

                     xx=global2local_i(xpart(ip,bi,bj),bi,bj,mythid)

                     scalex=recip_dxF(INT(xx),INT(yy),bi,bj)

                     scaley=recip_dyF(INT(xx),INT(yy),bi,bj)

                        call FLT_TRILINEAR(xx,yy,zz,uu,uVel,2,bi,bj)

                        call FLT_TRILINEAR(xx,yy,zz,vv,vVel,3,bi,bj)

                        call FLT_BILINEAR(xx,yy,uu,kp,uVel,2,bi,bj)

                        call FLT_BILINEAR(xx,yy,vv,kp,vVel,3,bi,bj)

                     xt=xx+0.5*deltaTmom*uu*scalex

      _RL xx, yy, xt, yt

            xx=global2local_i(xpart(ip,bi,bj),bi,bj,mythid)

               call FLT_BILINEAR2D(xx,yy,pp,cg2d_x,1,bi,bj)

               call FLT_BILINEAR  (xx,yy,uu,k,uVel,  2,bi,bj)

               call FLT_BILINEAR  (xx,yy,vv,k,vVel,  3,bi,bj)

               call FLT_BILINEAR  (xx,yy,tt,k,theta, 1,bi,bj)

               call FLT_BILINEAR  (xx,yy,ss,k,salt,  1,bi,bj)

      _RL xx, yy, xlo, xhi, ylo, yhi

      _RL xx, yy

      _RL xx, yy, zz, phi, scalex, scaley, scalez

      _RL xx, yy, phi, scalex, scaley

      _RL xx, yy, phi, scalex, scaley

         xx=trant(i,j)*dlayer(i,j,k2)

         flcu(i,j,k1) =flcu(i,j,k1)+xx

         flxu(i,j,k1) =flxu(i,j,k1)+xx*fclr(i,j)

         xx=trant(i,j)*dlayer(i,j,k1)

         flcd(i,j,k2) =flcd(i,j,k2)+xx

         flxd(i,j,k2) =flxd(i,j,k2)+xx*fclr(i,j)

      _RL xx,xx1,xx2,xx3

           xx=dh2o(i,j,k)*(pa(i,j,k)/500.0)

           h2oexp(i,j,k,1)=exp(-xx*0.10624)

           xx=h2oexp(i,j,k,1)*h2oexp(i,j,k,1)

           xx=xx*xx

           h2oexp(i,j,k,2)=xx*xx

           xx=h2oexp(i,j,k,2)*h2oexp(i,j,k,2)

           xx=xx*xx

           h2oexp(i,j,k,3)=xx*xx

           xx=h2oexp(i,j,k,3)*h2oexp(i,j,k,3)

           xx=xx*xx

           h2oexp(i,j,k,4)=xx*xx

           xx=h2oexp(i,j,k,4)*h2oexp(i,j,k,4)

           xx=xx*xx

           h2oexp(i,j,k,5)=xx*xx

           xx=h2oexp(i,j,k,5)*h2oexp(i,j,k,5)

           xx=xx*xx

           h2oexp(i,j,k,6)=xx*xx

           xx=dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5

           co2exp(i,j,k,1,1)=exp(-xx*2.656e-5)

           xx=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)

           xx=xx*xx

           co2exp(i,j,k,2,1)=xx*xx

           xx=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)

           xx=xx*xx

           co2exp(i,j,k,3,1)=xx*xx

           xx=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)

           xx=xx*xx

           co2exp(i,j,k,4,1)=xx*xx

           xx=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)

           xx=xx*xx

           co2exp(i,j,k,5,1)=xx*xx

           xx=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)

           xx=xx*xx

           co2exp(i,j,k,6,1)=xx*xx

           xx=dn2o(i,j,k)*(1.+(1.4476e-3+3.6656e-6*dt(i,j,k))*dt(i,j,k))

           n2oexp(i,j,k,1)=exp(-xx*0.25238)

           xx=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)

           xx1=xx*xx

           n2oexp(i,j,k,2)=xx*xx1*xx2*xx3

      _RL xx

           xx=   0.3153*th2o(i,j,1)

           xx=xx+0.4604*th2o(i,j,2)

           xx=xx+0.1326*th2o(i,j,3)

           xx=xx+0.0798*th2o(i,j,4)

           xx=xx+0.0119*th2o(i,j,5)

           tran(i,j)=tran(i,j)*xx

           xx=    0.2673*tco2(i,j,1,1)

           xx=xx+ 0.2201*tco2(i,j,2,1)

           xx=xx+ 0.2106*tco2(i,j,3,1)

           xx=xx+ 0.2409*tco2(i,j,4,1)

           xx=xx+ 0.0196*tco2(i,j,5,1)

           xx=xx+ 0.0415*tco2(i,j,6,1)

           tran(i,j)=tran(i,j)*xx

           xx=   0.970831*tn2o(i,j,1)

           xx=xx+0.029169*tn2o(i,j,2)

           tran(i,j)=tran(i,j)*(xx-1.0)

      _RL dp,xx,p1,dwe,dpe,a1,b1,fk1,a2,b2,fk2

         xx=pa(i,j,k)*0.001618*wa(i,j,k)*wa(i,j,k)*dp

         dcont(i,j,k) = xx*exp(1800./ta(i,j,k)-6.081)+1.e-10

      _RL  zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2,

                xx  = one-ff*ssc

                taup= tau*xx

                sscp= ssc*(one-ff)/xx

                xx  =  three*gp 

                gm1 =  (seven - sscp*(four+xx))*fourth

                gm2 = -(one   - sscp*(four-xx))*fourth

                gm3 =  (two   -        zth*xx )*fourth

                xx  = gm1-gm2

                akk = sqrt((gm1+gm2)*xx)

                alf1 = gm1 - gm3 * xx

                alf2 = gm2 + gm3 * xx

                xx  = akk * two

                all1 = (gm3 - alf2 * zth    )*xx*td

                bll = (one - gm3 + alf1*zth)*xx

                xx  = akk * zth

                st7 = one - xx

                st8 = one + xx

                xx  = akk * gm3

                cll = (alf2 + xx) * st7

                dll = (alf2 - xx) * st8

                xx  = akk * (one-gm3)

                fll = (alf1 + xx) * st8

                ell = (alf1 - xx) * st7

      _RL  xx,uuu,ttt,emt,up1,um1,st1

             xx  = one-ssc*g0

             uuu = sqrt( xx/(one-ssc))

             ttt = sqrt( xx*(one-ssc)*three )*tau

             xx  = um1*emt

             st1 = one / ((up1+xx) * (up1-xx))

      _RL  denm,xx

         xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is)

         fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm

         fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm

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

          xx=1./.3

          ic=int( (clog1+3.15)*xx+1.)

          iw=int( (wlog+4.15)*xx+1.)

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

          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

      _RL  zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2,

                xx  = one-ff*ssc

                taup= tau*xx

                sscp= ssc*(one-ff)/xx

                xx  =  three*gp 

                gm1 =  (seven - sscp*(four+xx))*fourth

                gm2 = -(one   - sscp*(four-xx))*fourth

                gm3 =  (two   -        zth*xx )*fourth

                xx  = gm1-gm2

                akk = sqrt((gm1+gm2)*xx)

                alf1 = gm1 - gm3 * xx

                alf2 = gm2 + gm3 * xx

                xx  = akk * two

                all1 = (gm3 - alf2 * zth    )*xx*td

                bll = (one - gm3 + alf1*zth)*xx

                xx  = akk * zth

                st7 = one - xx

                st8 = one + xx

                xx  = akk * gm3

                cll = (alf2 + xx) * st7

                dll = (alf2 - xx) * st8

                xx  = akk * (one-gm3)

                fll = (alf1 + xx) * st8

                ell = (alf1 - xx) * st7

      _RL  xx,uuu,ttt,emt,up1,um1,st1

             xx  = one-ssc*g0

             uuu = sqrt( xx/(one-ssc))

             ttt = sqrt( xx*(one-ssc)*three )*tau

             xx  = um1*emt

             st1 = one / ((up1+xx) * (up1-xx))

      _RL  denm,xx

         xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is)

         fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm

         fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm

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

          xx=1./.3

          ic=int( (clog1+3.15)*xx+1.)

          iw=int( (wlog+4.15)*xx+1.)

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

          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

      _RL     xx(nv)

     I                          xx

      call OPTIM_WRITEDATA( nv, ctrlname, lheaderonly, ff, xx )

      _RL   xx(nn)

     O                    xx

      call LSOPT_TOP( nn, xx, objf, adxx

      _RL   xx(nn)

      _RL   xx

      pointer (pxx,xx(1))

      _RL   xx(nmax)

      call OPTIM_INITMOD( nn, xx )

      _RL   xx(nn)

     I                  xx,

      call OPTIM_READDATA( nn, ctrlname, lheaderonly, objf,   xx )

      print *, ' leaving simul with xx, adxx = ', xx(1), adxx(1)

            xx(i) = xdiff(i)

            z = xx(i) + tact*dd(i)

            if ((z.ne.xx(i)) .and. (z.ne.xdiff(i))) go to 950

               xx(i) = xx(i) + tg*dd(i)

            xdiff(i) = xx(i) + tact*dd(i)

         xx(i) = xx(i) + tg*dd(i)

      double precision xx(nn), dd(nn), gg(nn), xdiff(nn)

     &     , dd, gg, xx, xdiff

         print *, 'pathei-lsopt xx(1), gg(1) ', xx(1), gg(1)

         call SIMUL( indic, nn, xx, ff, gg )

         print *, 'pathei: nach simul:   xx(1), gg(1) = ', xx(1), gg(1)

         call DOSTORE( nn, xx, .true., 1 )

      subroutine LSOPT_TOP( nn, xx, ff, gg, simul, optline

         call DOSTORE( nn, xx, .false., 1 )

     $           , DNRM2( nn, xx, 1 ), 0.

     &     , gg, dd, xx, xdiff

     &     , dd, gg, xx, xdiff

      call DOSTORE( nn, xx, .true., 1 )

     $     , DNRM2( nn, xx, 1 ), tact*DNRM2( nn, dd, 1 )

     &        , gg, dd, xx, xdiff

      double precision  xx(nn), ff, gg(nn), epsx, fmin, epsg

     $        , '  norm of x...................', DNRM2( nn, xx, 1 )

         xdiff(i) = xx(i) + tact*dd(i)

      double precision    gg(nn), dd(nn), xx(nn), xdiff(nn)

     &     , gg, dd, xx, xdiff

      real*4 xx(NN), gg(NN)

      subroutine LSWRI( isize, phniter, NN, xx, gg, lphprint )