call COLUMN(m,n,np,pa,dt,dco2,sco3,scopre,scotem)

        call COLUMN(m,n,np,pa,dt,do3,sco3,scopre,scotem)

        call H2OEXPS(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)

        call CO2EXPS(m,n,np,dco2,pa,dt,co2exp)

        call N2OEXPS(ib,m,n,np,dn2o,pa,dt,n2oexp)

        call CH4EXPS(ib,m,n,np,dch4,pa,dt,ch4exp)

        call COMEXPS(ib,m,n,np,dco2,dt,comexp)

        call CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df11,dt,f11exp)

        call CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df12,dt,f12exp)

        call CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df22,dt,f22exp)

        call B10EXPS(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt

      subroutine COLUMN (m,n,np,pa,dt,sabs0,sabs,spre,stem)

      _RL pa(m,n,np),dt(m,n,np),sabs0(m,n,np)

           stem(i,j,k+1)=stem(i,j,k)+dt(i,j,k)*sabs0(i,j,k)

      subroutine H2OEXPS(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)

      _RL dh2o(m,n,np),pa(m,n,np),dt(m,n,np)

     1        * ( 1.+(aw(ib)+bw(ib)* dt(i,j,k))*dt(i,j,k) )

      subroutine CO2EXPS(m,n,np,dco2,pa,dt,co2exp)

      _RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)

     1             *(1.+(0.0182+1.07e-4*dt(i,j,k))*dt(i,j,k))

     1             *(1.+(0.0042+2.00e-5*dt(i,j,k))*dt(i,j,k))

      subroutine N2OEXPS(ib,m,n,np,dn2o,pa,dt,n2oexp)

      _RL dn2o(m,n,np),pa(m,n,np),dt(m,n,np)

           xc=dn2o(i,j,k)*(1.+(1.9297e-3+4.3750e-6*dt(i,j,k))*dt(i,j,k))

     *        *(1.+(1.3804e-3+7.4838e-6*dt(i,j,k))*dt(i,j,k))

      subroutine CH4EXPS(ib,m,n,np,dch4,pa,dt,ch4exp)

      _RL dch4(m,n,np),pa(m,n,np),dt(m,n,np)

           xc=dch4(i,j,k)*(1.+(1.7007e-2+1.5826e-4*dt(i,j,k))*dt(i,j,k))

     *       *(1.+(5.9590e-4-2.2931e-6*dt(i,j,k))*dt(i,j,k))

      subroutine COMEXPS(ib,m,n,np,dcom,dt,comexp)

      _RL dcom(m,n,np),dt(m,n,np)

           xc=dcom(i,j,k)*(1.+(3.5775e-2+4.0447e-4*dt(i,j,k))*dt(i,j,k))

           xc=dcom(i,j,k)*(1.+(3.4268e-2+3.7401e-4*dt(i,j,k))*dt(i,j,k))

      subroutine CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,dcfc,dt,cfcexp)

      _RL dcfc(m,n,np),dt(m,n,np)

           xf=dcfc(i,j,k)*(1.+(a1+b1*dt(i,j,k))*dt(i,j,k))

           xf=dcfc(i,j,k)*(1.+(a2+b2*dt(i,j,k))*dt(i,j,k))

      subroutine B10EXPS(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt

      _RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)

     1           *(1.+(0.0149+6.20e-5*dt(i,j,k))*dt(i,j,k))

     1           *(1.+(0.0179+1.02e-4*dt(i,j,k))*dt(i,j,k))

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

      _RL pa(m,n,np),dt(m,n,np)

         dt(i,j,k)=ta(i,j,k)-250.0

        call COLUMN(m,n,np,pa,dt,dh2o,sh2o,swpre,swtem)

      _RL dt

      FRAC = 0.5/dt

       tx1(i) = cmb2pa *  (prs(i,ib+1) - prs(i,ib))/(grav*dt)

       cmass(i,L) = rhfrac(i)*cmass(i,L) * dt

       rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i)

      _RL   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)

      parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)

           ft=(log10(taux)-t1)/dt

       _RL dt

       subroutine STEP_FIZHI_CORR (myTime, myIter, myThid, dt)

       dtinv = 1. _d 0 / dt

       _RL dt

        call STEP_PHYSICS(uphy,vphy,thphy,sphy,dt,idim1,idim2,

      print *,' In step fizhi fg, new fizhi fields ',bi,' dt= ',dt

       subroutine STEP_FIZHI_FG (myTime, myIter, myThid, dt)

       _RL dt

         uphy(i,j,L,bi,bj) = uphy(i,j,L,bi,bj) + duphy(i,j,L,bi,bj) * dt

         vphy(i,j,L,bi,bj) = vphy(i,j,L,bi,bj) + dvphy(i,j,L,bi,bj) * dt

         thphy(i,j,L,bi,bj)=thphy(i,j,L,bi,bj) + dthphy(i,j,L,bi,bj)* dt

         sphy(i,j,L,bi,bj) = sphy(i,j,L,bi,bj) + dsphy(i,j,L,bi,bj) * dt

       subroutine STEP_PHYSICS(uphy,vphy,thphy,sphy,dt,

             dt = min(dt,dtts)

               write(0,*)"Time step is ", dt

     *                  (mda(k)*(tda(k)-taa(k)))*dt*recip_drF(k)

     *                  (mda(k)*(sda(k)-saa(k)))*dt*recip_drF(k)

     *              *dt*recip_drF(k2)

     *              *dt*recip_drF(k2)

     *                  dt*recip_drF(kmx+1)

     *                  dt*recip_drF(kmx+1)

        dt = dtts

         IF ( Wd(K2) .NE. 0. _d 0 ) dt = min(dt,drF(k2)/Wd(k2))

         ntime = nint(0.5*int(dtts/dt))

        dt = min(dt,dtts)

     &        'S/R OPPS_CALC: time step = ', dt

     &        (mda(k)*(Pda(k,ktr)-Paa(k,ktr)))*dt*recip_drF(k)

     &           *dt*recip_drF(k2)

     &        mda(kmx)*(Pda(kmx,ktr)-Paa(kmx,ktr))*dt*recip_drF(kmx+1)

      real dtts,dt,sumo,sumn

             dt = dtts

                dt = min(dt,dz(k2)/wd(k2))

                ntime = nint(0.5*int(dtts/dt))

      _RL dtts,dt

      _RL  dt

      dt  = thSIce_deltaT

          etop = (1. _d 0-frace)*sHeating * dt

          etope = frace*sHeating * dt

          esurp = sHeating * dt

      ebot = (flxCnB(i,j)-Fbot) * dt

      fresh = (mwater0 - (rhos*hs + rhoi*hi))/dt

        hs = hs + dt*snowPr/rhos

         IF (evap/rhos *dt.GT.hs) THEN

           evap=evap-hs*rhos/dt

           hs = hs - evap/rhos *dt

               dhi = evap/rhoi*dt

                 evap=evap-hnew(k)*rhoi/dt

         fresh = fresh + (rhos*hs + rhoi*hi)/dt

         fresh = fresh + chi*rhoi/dt

       qleft = qleft + (Fbot+(esurp+etop+ebot)/dt)

     &  'ThSI_CALC_TH: [esurp,etop+ebot]/dt =',esurp/dt,etop/dt,ebot/dt

      fsalt(i,j) = (msalt0 - rhoi*hi*saltIce)/dt

     &   (mwater0-(rhos*hs+rhoi*hi))/dt,evap,fresh,fsalt(i,j)

     &            Qleft,Fbot,(etope+ebote)/dt

         freshe=(rhos*hs+rhoi*hi)/dt

         salte=(rhoi*hi*saltIce)/dt

           qleft=qleft+(extend-rqh)/dt

         qleft=qleft+extend/dt

      _RL  dt

      dt  = thSIce_dtTemp

      a10 = rhoi*cpIce *hi/(2. _d 0*dt) +

     &      k32 * (4. _d 0*dt*k32 + rhoi*cpIce *hi)

     &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi)

     &       /(2. _d 0*dt)

     &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpIce *hi*Tice(2))

     &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi) - fswint

      c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)

      Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)

     &         /(6. _d 0*dt*k32 + rhoi*cpIce *hi)