_RL q(sNx,sNy,Nrphys,ntracer)

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

     . xlons,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

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

      integer k,l,m,q,bi,bj,iG,jG, myThid,num_file,length_for_tile

        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

     & k,q,tmp_i(k,q),tmp_j(k,q),tmp_weights(k,q)

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

     & k,q,tmp_weights(k,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,i,j,k,kk,kcur,iG,jG,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,i,j,k,kk,kcur,iG,jG,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 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)

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