c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE residua(opt)

c     ******************************************************************
c
c     This routine computes the residual for each model equation
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Argument variable

      character*3 opt


c---- Local variables

      integer	i,j,k,kcomp,kcomp0

      real*8	S,imbal

      real*8	conse_m,
     >		therm_U,
     >		therm_V,
     >		conse_r,
     >		conse_p




c     -------------------------------
c     Residuals for mass conservation
c     -------------------------------

      kcomp = 0

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (flag_mas(i,j,k)) then
           kcomp          = kcomp + 1
c          -----------------------------------
           res_mas(kcomp) = conse_m('R',i,j,k)
c          ----------------------------------- 
           if (opt.eq.'nor') res_mas(kcomp) =
     >			     res_mas(kcomp) / e_mas
         endif

      enddo
      enddo
      enddo




c     ------------------------------
c     Residuals for thermal wind (U)
c     ------------------------------

      kcomp = 0

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot-1

         if (flag_twu(i,j,k)) then
             kcomp          = kcomp + 1
c            ------------------------------------------------------
             res_twu(kcomp) = U(i,j,k+1)/dz(k+1) - 
     >			      U(i,j,k  )/dz(k  ) - therm_U(i,j,k+1) 
c            ------------------------------------------------------
             res_twu(kcomp) =
     >	     res_twu(kcomp) * (dz(k)+dz(k+1))/2.
             if (opt.eq.'nor') res_twu(kcomp) =
     >			       res_twu(kcomp) / e_twu
         endif

      enddo
      enddo
      enddo     




c     ------------------------------
c     Residuals for thermal wind (V)
c     ------------------------------

      kcomp = 0

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot-1

         if (flag_twv(i,j,k)) then
             kcomp          = kcomp + 1
c            -----------------------------------------------------
             res_twv(kcomp) = V(i,j,k+1)/dz(k+1) - 
     >			      V(i,j,k  )/dz(k  ) - therm_V(i,j,k+1)
c            ------------------------------------------------------
             res_twv(kcomp) =
     >	     res_twv(kcomp) * (dz(k)+dz(k+1))/2.
             if (opt.eq.'nor') res_twv(kcomp) =
     >			       res_twv(kcomp) / e_twv
         endif

      enddo
      enddo
      enddo    




c     -------------------------------
c     Residuals for vorticity balance
c     -------------------------------

      kcomp = 0

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (flag_vor(i,j,k)) then

           kcomp          = kcomp + 1

           res_vor(kcomp) = (sin(lat(j+1))-sin(lat(j)))/2. 
           res_vor(kcomp) =
     >	   res_vor(kcomp) * (V(i,j  ,k)/sin(lat(j  ))+
     >			     V(i,j+1,k)/sin(lat(j+1)))
           res_vor(kcomp) =
     >	   res_vor(kcomp) - W(i,j,k+1) + W(i,j,k) 

           if (opt.eq.'nor') res_vor(kcomp) =
     >			     res_vor(kcomp) / e_vor

         endif

      enddo
      enddo
      enddo




c     ------------------------------
c     Residuals for density equation
c     ------------------------------

      kcomp = 0

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (flag_den(i,j,k)) then
           kcomp          = kcomp + 1
c          -----------------------------------
           res_den(kcomp) = conse_r('R',i,j,k)
c          -----------------------------------
           if (opt.eq.'nor') res_den(kcomp) =
     >			     res_den(kcomp) / abs(e_den)
         endif 

      enddo
      enddo
      enddo




c     -----------------------------
c     Residuals for tracer equation
c     -----------------------------

      kcomp = 0

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (flag_pos(i,j,k)) then
           kcomp          = kcomp + 1
c          -----------------------------------
           res_pos(kcomp) = conse_p('R',i,j,k)
c          -----------------------------------
           if (opt.eq.'nor') res_pos(kcomp) =
     >                       res_pos(kcomp) / abs(e_pos)
         endif

      enddo
      enddo
      enddo




c     -------------------
c     End of the routine
c     -------------------

      return

      end

c     ******************************************************************

