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,

      _RL dtts,dt

        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))

             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)

     &        + (Fbot(i,j)+(esurp(i,j)+etop(i,j)+ebot(i,j))/dt)

     &        esurp(i,j)/dt,etop(i,j)/dt,ebot(i,j)/dt

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

     &         (mwater0(i,j)-(rhos*hSnow1(i,j)+rhoi*hIce(i,j)))/dt,

     &         flx2oc(I,J),Fbot(i,j),(etope(i,j)+ebote(i,j))/dt

          freshe=(rhos*hSnow1(i,j)+rhoi*hIce(i,j))/dt

          salte=(rhoi*hIce(i,j)*saltIce)/dt

           flx2oc(i,j)=flx2oc(i,j)+(extend-rqh)/dt

          flx2oc(i,j)=flx2oc(i,j)+extend/dt

      _RL dt

      dt  = thSIce_deltaT

          etop(i,j) = (1. _d 0-frace)*sHeat(i,j) * dt

          etope(i,j) = frace*sHeat(i,j) * dt

          esurp(i,j) = sHeat(i,j) * dt

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

         frzSeaWat(i,j) = rhoi*dhi/dt

     &        - (rhos*hSnow1(i,j)+rhoi*hIce(i,j)))/dt

          hSnow1(i,j) = hSnow1(i,j) + dt*snowP(i,j)/rhos

           IF (evapLoc(i,j)/rhos *dt.GT.hSnow1(i,j)) THEN

            evapLoc(i,j)=evapLoc(i,j)-hSnow1(i,j)*rhos/dt

            hSnow1(i,j) = hSnow1(i,j) - evapLoc(i,j)/rhos *dt

             dhi = evapLoc(i,j)/rhoi*dt

              evapLoc(i,j)=evapLoc(i,j)-hnew(i,j,k)*rhoi/dt

     &          + (rhos*hSnow1(i,j) + rhoi*hIce(i,j))/dt

           frw2oc(i,j) = frw2oc(i,j) + chi*rhoi/dt

      _RL  dt

      dt  = thSIce_dtTemp

          a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +

     &          k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )

     &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )

     &           /(2. _d 0*dt)

     &        - k32*( 4. _d 0*dt*k32*tFrz(i,j)

     &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )

          c10(i,j) = rhoi*Lfresh*hIce(i,j)*Tmlt1 / (2. _d 0*dt)

          tIc2(i,j) = ( 2. _d 0*dt*k32*(tIc1(i,j)+2. _d 0*tFrz(i,j))

     &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))