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

        call STEP_PHYSICS(uphy,vphy,thphy,sphy,dt,im1,im2,jm1,jm2,

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

       _RL dt

       dt = 1.

        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

         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,