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 xx, yy

        xx = xG(i,j,bi,bj) + ddx*( xG(i+1,j,bi,bj) - xG(i,j,bi,bj) )

        xx = xG(i,j,bi,bj) + ddx*( xG(i+1,j,bi,bj) - xG(i,j,bi,bj) )

        xx = xx_ij + ddx*( xx_ip1j - xx_ij )

     I                               xx, yy, bi, bj, myThid )

      _RL xx, yy

           IF ( xG(i,j,bi,bj).LE.xx .AND. xx.LT.xG(i+1,j,bi,bj) ) THEN

             fm = ( xx - xG(i,j,bi,bj) ) / dist

     O                               xx, yy,

            CALL FLT_MAP_IJLOCAL2XY( xx, yy,

            tmp(3) = xx

      _RL ix, jy, i0x, j0y, xx, yy, zz

             CALL FLT_MAP_IJLOCAL2XY( xx, yy,

             tmp(3) = xx

      _RL ix, jy, i0x, j0y, xx, yy, zz

      LOGICAL FUNCTION QUADROOT(aa,bb,cc,xx)

          _RL xx(1:2)

              xx(1) =  - bb + sq

              xx(2) =  - bb - sq

              xx(1) = xx(1) * aa

              xx(2) = xx(2) * aa

              xx(1) =  - cc / bb

              xx(2) =  - cc / bb

      _RL FUNCTION IDEMIX_gofx2(xx)

      _RL x,c,xx

      x=MAX(3.d0,xx)

     I                          xx,n,m,sNx,sNy,x,

      _RL     xx(1:n+1)

          IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR.

     &         ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN

        IF ( x(i,j).LT.xx(2) ) THEN

        ELSE IF ( x(i,j).GE.xx(n) ) THEN

        IF (x(i,j) .GE. xx(n)) THEN

        ELSE IF (x(i,j) .LT. xx(2)) THEN

        ELSE IF ( (x(i,j) .GE. xx(k(i,j)))

     &    .AND.   (x(i,j) .LT. xx(k(i,j)+1)) ) THEN

        ELSE IF (x(i,j) .GE. xx(k(i,j))) THEN

         DO WHILE (x(i,j) .GE. xx(k(i,j)+1))

        ELSE IF (x(i,j) .LT. xx(k(i,j)+1)) THEN

         DO WHILE (x(i,j) .LT. xx(k(i,j)))

      _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(nv)

     I                          xx

      call OPTIM_WRITEDATA( nv, ctrlname, lheaderonly, ff, 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 )