c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE A_modxx

c     ******************************************************************
c
c     This routine constructs the matrix A of the linear system Ax=b
c     corresponding to the linear equations of the geostrophic model:
c
c     - mass conservation
c     - thermal wind (u)
c     - thermal wind (v)
c     - linear vorticity balance
c
c     ******************************************************************

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

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

c---- Local variables

      integer	i,j,k,n,n1,n2,kcomp

      character*12 stringx

      character*12 stringr_s,stringr_n,stringr_w,stringr_e,
     >		   stringU_l,stringU_u,stringU_w,stringU_e,
     >		   stringV_l,stringV_u,stringV_s,stringV_n,
     >		   stringW_l,stringW_u

      real*8	dz_m




c     --------------------------------
c     Initializing counting for A rows
c     --------------------------------

c---- Number of linear constraint 

      kcomp = 0


c---- Number of each linear constraint

      n_mas = 0
      n_twu = 0
      n_twv = 0
      n_vor = 0      


c---- Matrix of linear constraint

      do n1=1,idim_n_li
      do n2=1,idim_n_o
         A_mod(n1,n2) = 0.d0
      enddo 
      enddo




c     ----------------------------
c     Mass conservation statements
c     ----------------------------

c---- General loop on statements or grid boxes

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


c---- Flag test

c++
      if (flag_mas(i,j,k)) then
c++


c---- Statement counter

      kcomp = kcomp + 1
      n_mas = n_mas + 1


c---- Variables required for grid box (i,j,k)

c     -------------------------------------
      stringU_w = stringx('UU',i  ,j  ,k  )
      stringU_e = stringx('UU',i+1,j  ,k  )
      stringV_s = stringx('VV',i  ,j  ,k  )
      stringV_n = stringx('VV',i  ,j+1,k  )
      stringW_l = stringx('WW',i  ,j  ,k  )
      stringW_u = stringx('WW',i  ,j  ,k+1)
c     -------------------------------------


c---- Elements in the kcomp-th row of A

c++
      do n=1,n_o
c++

         if (string(n).eq.stringU_w .or.
     >	     string(n).eq.stringV_s .or.
     >	     string(n).eq.stringW_l) then

             A_mod(kcomp,n) = -1.d0

         else if (string(n).eq.stringU_e .or.
     >		  string(n).eq.stringV_n .or.
     >		  string(n).eq.stringW_u) then

             A_mod(kcomp,n) = +1.d0

         else

             A_mod(kcomp,n) =  0.d0

         endif

c++
      enddo
c++


c---- End of flag test

c++
      endif
c++ 


c---- End of loop 

c     *****
      enddo
      enddo
      enddo
c     *****




c     -------------------------
c     Thermal wind U statements 
c     -------------------------

c---- General loop on statements or meridional strips i,j,k -> i,j+1,k 

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


c---- Flag test 

c++
      if (flag_twu(i,j,k)) then
c++


c---- Statement counter

      kcomp = kcomp + 1
      n_twu = n_twu + 1


c---- Variables required for meridional strip i,j,k+1 -> i,j+1,k+1 

c     -----------------------------------
      stringU_u = stringx('UU',i,j  ,k+1)
      stringU_l = stringx('UU',i,j  ,k  )
      stringr_s = stringx('rr',i,j  ,k+1)
      stringr_n = stringx('rr',i,j+1,k+1)
c     -----------------------------------


c---- Elements of kcomp-th row of A

c++
      do n=1,n_o
c++

         dz_m = (dz(k)+dz(k+1)) / 2.

         if (string(n).eq.stringU_u) then

           A_mod(kcomp,n) = +dz_m / dz(k+1)

         else if (string(n).eq.stringU_l) then

           A_mod(kcomp,n) = -dz_m / dz(k)

         else if (string(n).eq.stringr_s) then

           A_mod(kcomp,n) = -sv*g*dy(j) / (2*omega*rho_0)
           A_mod(kcomp,n) =
     >	   A_mod(kcomp,n) / (cos(lat(j+1))-cos(lat(j)))
           A_mod(kcomp,n) =
     >	   A_mod(kcomp,n) * dz_m**2

         else if (string(n).eq.stringr_n) then

           A_mod(kcomp,n) = +sv*g*dy(j) / (2*omega*rho_0)
           A_mod(kcomp,n) =
     >     A_mod(kcomp,n) / (cos(lat(j+1))-cos(lat(j)))
           A_mod(kcomp,n) =
     >     A_mod(kcomp,n) * dz_m**2 

         else

           A_mod(kcomp,n) = 0.d0

         endif

c++
      enddo
c++


c---- End of flag test 

c++
      endif
c++


c---- End of loop 

c     ******
      enddo
      enddo
      enddo
c     ******




c     -------------------------
c     Thermal wind V statements 
c     -------------------------

c---- General loop on statements or zonal strips i,j,k+1 -> i+1,j,k+1 

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


c---- Flag test 

c++
      if (flag_twv(i,j,k)) then 
c++


c---- Statement counter

      kcomp = kcomp + 1
      n_twv = n_twv + 1


c---- Variables required for zonal strip i,j,k -> i,j+1,k 

c     -----------------------------------
      stringV_u = stringx('VV',i  ,j,k+1)
      stringV_l = stringx('VV',i  ,j,k  )
      stringr_w = stringx('rr',i  ,j,k+1)
      stringr_e = stringx('rr',i+1,j,k+1)
c     -----------------------------------


c---- Elements of kcomp-th row of A

c++
      do n=1,n_o
c++

         dz_m = (dz(k)+dz(k+1)) / 2.

         if (string(n).eq.stringV_u) then

           A_mod(kcomp,n) = +dz_m / dz(k+1)

         else if (string(n).eq.stringV_l) then

           A_mod(kcomp,n) = -dz_m / dz(k)

         else if (string(n).eq.stringr_w) then

           A_mod(kcomp,n) = -sv*g / (2*omega*rho_0*sin(lat(j)))
           A_mod(kcomp,n) =
     >	   A_mod(kcomp,n) * dz_m**2 

         else if (string(n).eq.stringr_e) then

           A_mod(kcomp,n) = +sv*g / (2*omega*rho_0*sin(lat(j)))
           A_mod(kcomp,n) =
     >     A_mod(kcomp,n) * dz_m**2 

         else

           A_mod(kcomp,n) =  0.d0

         endif

c++
      enddo
c++


c---- End of flag test

c++
      endif
c++


c---- End of loop 

c     *****
      enddo
      enddo
      enddo
c     *****




c     ----------------------------
c     Vorticity balance statements 
c     ----------------------------

c---- General loop on statements or grid boxes 

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


c---- Flag test

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


c---- Statement counter

      kcomp = kcomp + 1
      n_vor = n_vor + 1


c---- Variables required for grid box (i,j,k)

c     -----------------------------------
      stringV_s = stringx('VV',i,j  ,k  )
      stringV_n = stringx('VV',i,j+1,k  )
      stringW_l = stringx('WW',i,j  ,k  )
      stringW_u = stringx('WW',i,j  ,k+1)
c     -----------------------------------


c---- Elements of the kcomp-th row of A

c++
      do n=1,n_o
c++

         if (string(n).eq.stringV_s) then

           A_mod(kcomp,n) =  (  sin(lat(j+1))-sin(lat(j))) /
     >		  	     (2*              sin(lat(j))) 

         else if (string(n).eq.stringV_n) then

           A_mod(kcomp,n) =  (  sin(lat(j+1))-sin(lat(j))) /
     >                       (2*sin(lat(j+1            )))

         else if (string(n).eq.stringW_l) then

           A_mod(kcomp,n) = +1.d0 

         else if (string(n).eq.stringW_u) then

           A_mod(kcomp,n) = -1.d0 

         else

           A_mod(kcomp,n) =  0.d0

         endif

c++
      enddo
c++


c---- End of flag test

c++
      endif
c++


c---- End of loop 

c     *****
      enddo
      enddo
      enddo
c     *****




c     --------------------------------------------------------
c     Number of rows of A_mod (= # of linear model statements)
c     --------------------------------------------------------

      n_li = kcomp




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

      return

      end

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

