_RL wsqr,oldflux,newflux,entrainrate

         entrainrate = (newflux - oldflux)/newflux

     &         - entrainrate*(pmix - Pd(k2,ktr))

     &    'S/R OPPS_CALC: Dd,De,entr,k ',Dd(k2),De(k2),entrainrate,k2

         wsqr = wsqr - wsqr*abs(entrainrate)+ gravity*

     &          'S/R OPPS_CALC: ', wsqr, entrainrate,

      real pdensity,wsqr,oldflux,newflux,entrainrate,adtemp

                 entrainrate = (newflux - oldflux)/newflux

                  sd(k2+1)=sd(k2) - entrainrate*(smix - sd(k2))

                  td(k2+1)=td(k2) - entrainrate*(thmix - td(k2))

                  write(0,*)"Dd,De,entr,k ",Dd(k2),De(k2),entrainrate,k2

                 wsqr = wsqr - wsqr*abs(entrainrate)+ 9.81*