_RL q(sNx,sNy,Nrphys,ntracer)

       q(i,j,L,1) = sphy(i,j,L,bi,bj)

     & p(1,1,bi,bj),u,v,t,q,pl(1,1,1,bi,bj),ple(1,1,1,bi,bj),

      call FIZHI_STEP_DIAG(myid,p,u,v,t,q,qq,pkl,dpres,

     &  ptop,p,u,v,t,q,ntracer,ptracer,pl,ple,dpres,pkht,pkl,

      call GWDRAG (myid,p,pl,ple,dpres,pkl,u,v,t,q,phis_var,

      call RAYLEIGH(myid,pl,pkl,pkht,surfz,u,v,t,q,im,jm,lm,

     &    p,pl,ple,dpres,pkht,pkl,u,v,t,q,bi,bj,ntracer,ptracer,qq,

      _RL q(im,jm,lm,ntracer)

     & p,u,v,t,q,pl,ple,dpres,pkht,pkl,surfz,fracland,landtype,radswt,

      subroutine MAXVAL1 (q,im,qval,qmax,imax)

      _RL   q(im), qmax, qval

      if( q(i).gt.qmax ) then

      qmax = q(i)

      subroutine SRCLOUDS (th,q,plk,pl,plke,cloud,cldwat,irun,irise,

      _RL    q(irun,irise)

      rh    = q(i,L)/qs(i,L)

      call CTEI ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )

      subroutine CTEI ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )

      _RL  th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)

          dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)

      if( q(i).lt.minval ) minval = q(i)

      subroutine VQSAT (tt,p,q,dqdt,ldqdt,n)

      _RL tt(n), p(n), q(n), dqdt(n)

      call QSAT ( tt(i),p(i),q(i),dqdt(i),ldqdt )

      function minval (q,im)

      _RL q(im), minval

      if(2.eq.1 )then

      if(2.eq.1 )then

      _RL del_P ,del_th, th, q

      q      = del_th

      th     = th + (1 - 1/sqrt(two))*(del_th - q)

      q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q

      th     = th + (1 + 1/sqrt(two))*(del_th - q)

      q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q

      rv = th + (del_th - two*q)/(two*three)

      _RL del_P ,del_th, th, q

      q      = del_th

      th     = th + (1 - 1/sqrt(two))*(del_th - q)

      q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q

      th     = th + (1 + 1/sqrt(two))*(del_th - q)

      q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q

      rv     = th + (del_th - two*q)/(two*three)

      do q = 1,NUM_INTERP_POINTS

         prof_interp_i(num_file,k,q,bi,bj) = -999

         prof_interp_j(num_file,k,q,bi,bj) = -999

         prof_interp_weights(num_file,k,q,bi,bj) = -999

      do q=1,NUM_INTERP_POINTS

        vec_start2(1)=q

     &  vec_count2, tmp_weights(1,q))

     &  vec_count2, tmp_i(1,q))

     &  vec_count2, tmp_j(1,q))

        do q = 1,NUM_INTERP_POINTS

             prof_interp_weights(num_file,length_for_tile,q,bi,bj)

     &       =tmp_weights(k,q)

             prof_interp_i(num_file,length_for_tile,q,bi,bj)

     &       =tmp_i(k,q)

             prof_interp_j(num_file,length_for_tile,q,bi,bj)

     &       =tmp_j(k,q)

             tmp_sum_weights=tmp_sum_weights+tmp_weights(k,q)

       if ( (tmp_i(k,q).LT.0).OR.(tmp_j(k,q).LT.0)

     & .OR.(tmp_i(k,q).GT.sNx+1).OR.(tmp_j(k,q).GT.sNy+1) ) then

       if ( tmp_weights(k,q) .NE. 0. ) then

       if ( ((tmp_i(k,q).EQ.0).AND.(tmp_j(k,q).EQ.0))

     & .OR.((tmp_i(k,q).EQ.sNx+1).AND.(tmp_j(k,q).EQ.sNy+1))

     & .OR.((tmp_i(k,q).EQ.0).AND.(tmp_j(k,q).EQ.sNy+1))

     & .OR.((tmp_i(k,q).EQ.sNx+1).AND.(tmp_j(k,q).EQ.0)) ) then

       if ( (tmp_weights(k,q).LT.0).OR.(tmp_weights(k,q).GT.1) ) then

      integer varid_intp3, varid_intp4, varid_intp5, q

      do q = 1,NUM_INTERP_POINTS

         prof_interp_i(num_file,k,q,bi,bj) = -999

         prof_interp_j(num_file,k,q,bi,bj) = -999

         prof_interp_weights(num_file,k,q,bi,bj) = -999

        ponderations(q)=tab_coeffs3(q)*weights_cur(q)

        pondsSUM=pondsSUM+ponderations(q)

         do q=1,NUM_INTERP_POINTS

       traj_cur(k)=traj_cur(k)+tab_coeffs1(q)*ponderations(q)/pondsSUM

      integer q,k,kk,kcur,bi,bj

       do q=1,NUM_INTERP_POINTS

               tab_coeffs1(q)=theta(i_cur(q),j_cur(q),k,bi,bj)

               tab_coeffs3(q)=maskC(i_cur(q),j_cur(q),k,bi,bj) 

               tab_coeffs1(q)=salt(i_cur(q),j_cur(q),k,bi,bj) 

               tab_coeffs3(q)=maskC(i_cur(q),j_cur(q),k,bi,bj) 

               tab_coeffs1(q)=uVel(i_cur(q),j_cur(q),k,bi,bj) 

               tab_coeffs3(q)=maskW(i_cur(q),j_cur(q),k,bi,bj)  

               tab_coeffs1(q)=vVel(i_cur(q),j_cur(q),k,bi,bj) 

               tab_coeffs3(q)=maskS(i_cur(q),j_cur(q),k,bi,bj)  

               tab_coeffs1(q)=pTracer(i_cur(q),j_cur(q),k,bi,bj,1) 

               tab_coeffs1(q)=0 

               tab_coeffs3(q)=maskC(i_cur(q),j_cur(q),k,bi,bj) 

               tab_coeffs1(q)=etan(i_cur(q),j_cur(q),bi,bj) 

               tab_coeffs3(q)=maskC(i_cur(q),j_cur(q),1,bi,bj) 

               tab_coeffs1(q)=0.

               tab_coeffs3(q)=0.

      integer q,k,kk,kcur,bi,bj

       do q=1,NUM_INTERP_POINTS

               tab_coeffs1(q)=prof_etan_mean(i_cur(q),j_cur(q),bi,bj) 

               tab_coeffs3(q)=maskC(i_cur(q),j_cur(q),1,bi,bj) 

               tab_coeffs1(q)=0.

               tab_coeffs3(q)=0.

        ponderations(q)=tab_coeffs3(q)*weights_cur(q)

        pondsSUM=pondsSUM+ponderations(q)

         do q=1,NUM_INTERP_POINTS

       traj_cur(k)=traj_cur(k)+tab_coeffs1(q)*ponderations(q)/pondsSUM

      _RL q(sNx,sNy,Nrphys,ntracer)

      _RL del_P ,del_th, th, q

      q      = del_th

      th     = th + (1 - 1/sqrt(two))*(del_th - q)

      q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q

      th     = th + (1 + 1/sqrt(two))*(del_th - q)

      q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q

      SW_temp= th + (del_th - two*q)/(two*three)

      _RL del_P ,del_th, th, q

      q      = del_th

      th     = th + (1 - 1/sqrt(two))*(del_th - q)

      q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q

      th     = th + (1 + 1/sqrt(two))*(del_th - q)

      q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q

      SW_PTMP     = th + (del_th - two*q)/(two*three)

            q=q100(i,j)

     &           q.GT.ssqt(i,j)*0.7 _d 0 ) THEN

     &               + ( (q-ssqt(i,j)*0.7 _d 0)*dm/cheap_pr2 )

            IF (q.GT.0.0214 _d 0 .AND. q.GT.ssqt(i,j)*0.9 _d 0) THEN

              precip = precip + ((q-ssqt(i,j)*0.9 _d 0)*dm/cheap_pr1)

            entrain = cdq(i,j)*q*0.25

        _RL q,precip,ttt,entrain

        q=qair(i,j,bi,bj)

        q=QaR*ssqt

      Dq=qs-q

     &          ( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )

     &         /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) )

     &    ' tsr,tta,q,qsr  =', tsr, tta, q, qsr

       zL=xkar*gravity*zu*( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )

     &                   /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) )

      _RL u,q,Tas,tta,zi,es,qs,tsw

            q=qair(i,j,bi,bj)

            q=QaR * ssq

            delq   = q - ssq

            t0     = ttas*(1. _d 0 + humid_fac*q)

     &                    qstar/(1. _d 0/humid_fac+q))

     &                    qstar/(1. _d 0/humid_fac+q))

     &                    + qstar/(1. _d 0/humid_fac+q))

      _RL ssq, q