c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE writing(iopt)

c     ******************************************************************
c
c     This routine writes outputs
c
c     ******************************************************************


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

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

c---- Argument variable

      integer	iopt


c---- Local variables

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

      character	form_dx*6,
     >		form_dy*6,
     >		form_r*8 ,
     >          form_U*9 ,
     >          form_V*9 ,
     >          form_W*9 ,
     >		form_p*8 ,
     >		form_Z*8

      character	string_r*12,
     >		string_c*12,
     >		string0 *12

      logical	flag

      integer	ijk(itot+1,jtot+1,ktot+1)

      real*8	r_w(itot+1,jtot+1,ktot+1)




c     --------
c     Formats
c     --------

c++
      if (itot.lt.10) then
c++

        write(form_dx,1) itot
1       format('(0',i1,'i3)')

        write(form_dy,2) itot+1
2       format('(0',i1,'i3)')

        write(form_r,3) itot+1
3       format('(0',i1,'f9.3)')

        write(form_U,4) itot+1
4       format('(0',i1,'e11.3)')

        write(form_V,5) itot
5       format('(0',i1,'e11.3)')

        write(form_W,6) itot
6       format('(0',i1,'e11.3)')

        write(form_p,7) itot+1
7       format('(0',i1,'f9.3)')

        write(form_Z,8) itot
8       format('(0',i1,'f9.3)')


c++
      else
c++

        write(form_dx,10) itot
10      format('(',i2,'i3)')

        write(form_dy,20) itot+1
20      format('(',i2,'i3)')

        write(form_r,30) itot+1
30      format('(',i2,'f9.3)')

        write(form_U,40) itot+1
40      format('(',i2,'f8.3)')

        write(form_V,50) itot
50      format('(',i2,'f8.3)')

        write(form_W,60) itot
60      format('(',i2,'f8.3)')

        write(form_p,70) itot+1
70      format('(',i2,'f9.3)')

        write(form_Z,80) itot
80      format('(',i2,'f9.3)')


c++
      endif
c++




c     ---------------------------------------------------
c     Write model set-up, prior xo, and error in prior xo
c     ---------------------------------------------------

c---- Test

c     ####################
      if (iopt.eq.0) then
c     ####################


c---- dx,dy,dz 

      write(9,*)
      write(9,*) 'LONGITUDE AND LATITUDE OF ORIGIN'
      write(9,*)
      write(9,'(2f8.1)') lon_w,lat_s

      write(9,*)
      write(9,*)
      write(9,*) 'FIELD OF DX [deg]'
      write(9,*)

      do j=jtot+1,1,-1
         write(9,form_dx) (nint(dx(i)*deg),i=1,itot)
      enddo

c--

      write(9,*)
      write(9,*)
      write(9,*) 'FIELD OF DY [deg]'
      write(9,*)

      do j=jtot,1,-1
         write(9,form_dy) (nint(dy(j)*deg),i=1,itot+1)
      enddo

c--

      write(9,*)
      write(9,*)
      write(9,*) 'VALUES OF DZ [m]'
      write(9,*)

      do k=ktot,1,-1
         write(9,'(i4)') nint(dz(k))
      enddo


c---- Grid cell depths

      write(9,*)
      write(9,*)
      write(9,*) 'GRID CELL DEPTHS'
      write(9,*)
      do j=jtot,1,-1
         write(9,form_Z) (z(i,j),i=1,itot)
      enddo


c---- Vorticity flags 

      do k=ktot,1,-1

         write(99,*)
         write(99,*)
         write(99,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(99,*)
         do j=jtot,1,-1
            write(99,*) (flag_vor(i,j,k),i=1,itot)
         enddo

      enddo


c---- Prior density 

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

         ijk(i,j,k) = 0
         if (flag_r(i,j,k)) ijk(i,j,k) = 1

          if (flag_r(i,j,k)) then
            r_w(i,j,k) = r(i,j,k) + rho_0
          else
            r_w(i,j,k) = r(i,j,k)
          endif

      enddo
      enddo
      enddo

c--

      do k=ktot+1,1,-1

         write(40,*)
         write(40,100) z_l(k)/1000. 
         write(40,*) 
         do j=jtot+1,1,-1
            write(40,form_r) (r_w(i,j,k),i=1,itot+1)
         enddo

      enddo


c---- Prior transport component U

      do k=ktot,1,-1 

         write(41,*)
         write(41,200) z_l(k+1)/1000.,z_l(k)/1000.
         write(41,*)
         do j=jtot,1,-1
            write(41,form_U) (U(i,j,k),i=1,itot+1)
         enddo 

      enddo


c---- Prior transport component V 

      do k=ktot,1,-1

         write(42,*)
         write(42,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(42,*)
         do j=jtot+1,1,-1
            write(42,form_V) (V(i,j,k),i=1,itot)
         enddo

      enddo


c---- Prior transport component W 

      do k=ktot+1,1,-1

         write(43,*)
         write(43,100) z_l(k)/1000. 
         write(43,*)
         do j=jtot,1,-1
            write(43,form_W) (W(i,j,k),i=1,itot)
         enddo

      enddo


c---- Prior tracer values 

      do k=ktot+1,1,-1

         write(44,*)
         write(44,100) z_l(k)/1000. 
         write(44,*) 
         do j=jtot+1,1,-1
            write(44,form_p) (p(i,j,k),i=1,itot+1)
         enddo

      enddo


c---- Error in prior density 

      do k=ktot+1,1,-1

         write(81,*)
         write(81,100) z_l(k)/1000.
         write(81,*)
         do j=jtot+1,1,-1
            write(81,form_r) (e_r(i,j,k),i=1,itot+1)
         enddo

      enddo


c---- Error in prior transport component U

      do k=ktot,1,-1 

         write(82,*)
         write(82,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(82,*)
         do j=jtot,1,-1
            write(82,form_U) (e_U(i,j,k),i=1,itot+1)
         enddo 

      enddo


c---- Error in prior transport component V

      do k=ktot,1,-1

         write(83,*)
         write(83,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(83,*)
         do j=jtot+1,1,-1
            write(83,form_V) (e_V(i,j,k),i=1,itot)
         enddo

      enddo


c---- Error in prior transport component W

      do k=ktot+1,1,-1

         write(84,*)
         write(84,100) z_l(k)/1000. 
         write(84,*)
         do j=jtot,1,-1
            write(84,form_W) (e_W(i,j,k),i=1,itot)
         enddo

      enddo


c---- Error in prior tracer values

      do k=ktot+1,1,-1

         write(85,*)
         write(85,100) z_l(k)/1000.
         write(85,*)
         do j=jtot+1,1,-1
            write(85,form_p) (e_p(i,j,k),i=1,itot+1)
         enddo

      enddo


c---- Exit the routine

      return


c---- End of test

c     ######
      endif
c     ######




c     ----------------
c     Rewind all files
c     ----------------

c---- r,U,V,W,p

      do ifile=10,14
         rewind(ifile)
      enddo


c---- e_r,e_U,e_V,e_W,e_p

      do ifile=60,64
         rewind(ifile)
      enddo


c---- res_mas,res_twu,res_twv,res_vor,res_den,res_pos,res_p

      do ifile=27,33
         rewind(ifile)
      enddo


c---- mis_r,mis_U,mis_V,mis_W,mis_p

      do ifile=22,26
         rewind(ifile)
      enddo


c---- P

      rewind (97)




c     -----------------
c     Write posterior x 
c     -----------------

c---- Posterior density 

      do k=ktot+1,1,-1

         write(10,*)
         write(10,100) z_l(k)/1000.
         write(10,*) 
         do j=jtot+1,1,-1
            write(10,form_r) (r(i,j,k),i=1,itot+1)
         enddo

      enddo

100   format('FIELD AT ',f3.1,' KM')


c---- Posterior transport component U

      do k=ktot,1,-1 

         write(11,*)
         write(11,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(11,*)
         do j=jtot,1,-1
            write(11,form_U) (U(i,j,k),i=1,itot+1)
         enddo 

      enddo

200   format('TRANSPORT BETWEEN ',f3.1,'-',f3.1,' KM')


c---- Posterior transport component V 

      do k=ktot,1,-1

         write(12,*)
         write(12,200) z_l(k+1)/1000.,z_l(k)/1000.
         write(12,*)
         do j=jtot+1,1,-1
            write(12,form_V) (V(i,j,k),i=1,itot)
         enddo

      enddo


c---- Posterior transport component W 

      do k=ktot+1,1,-1

         write(13,*)
         write(13,100) z_l(k)/1000.
         write(13,*)
         do j=jtot,1,-1
            write(13,form_W) (W(i,j,k),i=1,itot)
         enddo

      enddo


c---- Posterior tracer values 

      do k=ktot+1,1,-1

         write(14,*)
         write(14,100) z_l(k)/1000.
         write(14,*) 
         do j=jtot+1,1,-1
            write(14,form_p) (p(i,j,k),i=1,itot+1)
         enddo

      enddo




c     ---------------------------
c     Write errors in posterior x
c     ---------------------------

c---- Error in posterior density 

      do k=ktot+1,1,-1

         write(60,*)
         write(60,100) z_l(k)/1000. 
         write(60,*)
         do j=jtot+1,1,-1
            write(60,form_r) (e_r(i,j,k),i=1,itot+1)
         enddo

      enddo


c---- Error in posterior transport component U

      do k=ktot,1,-1 

         write(61,*)
         write(61,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(61,*)
         do j=jtot,1,-1
            write(61,form_U) (e_U(i,j,k),i=1,itot+1)
         enddo 

      enddo


c---- Error in posterior transport component V

      do k=ktot,1,-1

         write(62,*)
         write(62,200) z_l(k+1)/1000.,z_l(k)/1000. 
         write(62,*)
         do j=jtot+1,1,-1
            write(62,form_V) (e_V(i,j,k),i=1,itot)
         enddo

      enddo


c---- Error in posterior transport component W

      do k=ktot+1,1,-1

         write(63,*)
         write(63,100) z_l(k)/1000.
         write(63,*)
         do j=jtot,1,-1
            write(63,form_W) (e_W(i,j,k),i=1,itot)
         enddo

      enddo


c---- Error in posterior tracer values 

      do k=ktot+1,1,-1

         write(64,*)
         write(64,100) z_l(k)/1000. 
         write(64,*)
         do j=jtot+1,1,-1
            write(64,form_p) (e_p(i,j,k),i=1,itot+1)
         enddo

      enddo




c     ---------------------
c     Write model residuals
c     ---------------------

c---- Residuals for mass conservation
 
      kcomp_w = 0
      kcomp   = 0

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

         kcomp_w = kcomp_w + 1
         if (flag_mas(i,j,k)) then
           kcomp = kcomp + 1
           write(27,*) kcomp_w,res_mas(kcomp),i,j,k 
         endif

      enddo
      enddo
      enddo


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

      kcomp_w = 0
      kcomp   = 0

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

         kcomp_w = kcomp_w + 1
         if (flag_twu(i,j,k)) then 
             kcomp = kcomp + 1
             write(28,*) kcomp_w,res_twu(kcomp),i,j,k 
         endif

      enddo
      enddo
      enddo


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

      kcomp_w = 0
      kcomp   = 0

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

         kcomp_w = kcomp_w + 1
         if (flag_twv(i,j,k)) then
             kcomp = kcomp + 1
             write(29,*) kcomp_w,res_twv(kcomp),i,j,k 
         endif

      enddo
      enddo
      enddo


c---- Residuals for vorticity balance

      kcomp_w = 0
      kcomp   = 0

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

         kcomp_w = kcomp_w + 1
         if (flag_vor(i,j,k)) then
           kcomp = kcomp + 1
           write(30,*) kcomp_w,res_vor(kcomp),i,j,k 
         endif

      enddo
      enddo
      enddo


c---- Residuals for density equation

      kcomp_w = 0
      kcomp   = 0
     
      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         kcomp_w = kcomp_w + 1
         if (flag_den(i,j,k)) then
           kcomp = kcomp + 1
           write(31,*) kcomp_w,res_den(kcomp),i,j,k 
         endif

      enddo
      enddo
      enddo


c---- Residuals for tracer equation

      kcomp_w = 0
      kcomp   = 0

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

         kcomp_w = kcomp_w + 1
         if (flag_pos(i,j,k)) then
           kcomp = kcomp + 1
           write(32,*) kcomp_w,res_pos(kcomp),i,j,k
         endif

      enddo
      enddo
      enddo




c     ------------------
c     Write data misfits
c     ------------------

c---- Misfits for density

      kcomp = 0

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

         if (flag_r(i,j,k)) then
           kcomp = kcomp + 1
           write(22,*) kcomp,(x(kcomp)-xo(kcomp))*sqrt(Co_1(kcomp))
         endif 
    
      enddo
      enddo
      enddo 


c---- Misfits for transport component U

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

         if (flag_U(i,j,k)) then
           kcomp     = kcomp + 1
           write(23,*) kcomp,(x(kcomp)-xo(kcomp))*sqrt(Co_1(kcomp))
         endif

      enddo
      enddo
      enddo


c---- Misfits for transport component V

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

         if (flag_V(i,j,k)) then
           kcomp     = kcomp + 1
           write(24,*) kcomp,(x(kcomp)-xo(kcomp))*sqrt(Co_1(kcomp)) 
         endif

      enddo
      enddo
      enddo


c---- Misfits for transport component W

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

         if (flag_W(i,j,k)) then
           kcomp     = kcomp + 1
           write(25,*) kcomp,(x(kcomp)-xo(kcomp))*sqrt(Co_1(kcomp)) 
         endif

      enddo
      enddo
      enddo


c---- Misfits for tracer 

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

         if (flag_p(i,j,k)) then
           kcomp = kcomp + 1
           write(26,*) kcomp,(x(kcomp)-xo(kcomp))*sqrt(Co_1(kcomp))
         endif 
    
      enddo
      enddo
      enddo 




c     ---------------------------------------------
c     Write error covariance matrix for posterior x
c     ---------------------------------------------

      do n1=1,n_o
      do n2=1,n_o

         string_r = string(n1)
         string_c = string(n2)

         flag = .false.

         if (string_r(1:2).eq.'VV' .and. 
     >	     string_c(1:2).eq.'VV') then 

           if (string_r(7:8).eq.'08' .and. string_c(7:8).eq.'08' .or.
     >	       string_r(7:8).eq.'36' .and. string_c(7:8).eq.'36' .or.
     >	       string_r(7:8).eq.'42' .and. string_c(7:8).eq.'42') then

               flag = .true.

           endif

         endif

         if (iflag.eq.2 .or.
     >	     iflag.eq.3) flag = .false.

         if (flag) then
            write(97,111) string(n1),' ',string(n2),' ',Coa(n1,n2) 
         endif

      enddo
      enddo
 
111   format(a12,a1,a12,a1,f9.3)




c     ------------------
c     End of subroutine 
c     ------------------

      return
 
      end

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

