!---------------------------------------------------------------------------!
!- START of Subroutine: write_results
!---------------------------------------------------------------------------!
!! @AUTHOR L. Magnani
!! @VERSION 0.0
!- history:                                                                 !
!---date of creation: 17-09-04
!---modified by     : B.M. Dinelli
!---date            : 11/05/2005
!---Modification    : added lsouthpole and itimeini
!---modified by     : B.M. Dinelli
!---date            : 09/06/2005
!---Modification    : added the possibility of not writing retrieved values
!                     affected by clouds
!------------
!! @SUBJECT writes results files
!------------
!! @DESCRIPTION
!!>
!! The subroutine writes the output files with the initial guess and the
!! retrieved values of all the retrieval targets. 
!! If the user requires it, the retrieved values affected by clouds in the
!! measurements won't be written in the output files
!!<
!---------------------------------------------------------------------------!


      subroutine write_results(nf, orbit, date, rec_time, iseq1, iseq2,
     &     ilimb, ngrid, nret, rlatpar, rlongpar, difpar, qpar,
     &     rxpar, rvcmns, jvetpar, nflag, retr_time,
     &     itarget,indpar,itimeini,lsouthpole,lwritepar,Ppar,Tpar,
     &     l_optest,rapvcm)

      implicit none
      include 'parameters.inc'

!! @INPUT_VARIABLES
      integer*4 nf                   !! file number
      character*11 orbit             !! main orbit numbers
      character*11 date(2)           !! maesurements date
      character*8 rec_time(imxlmb,imxseq) !! sweep recording time
      character*20 retr_time(imxlmb,imxseq) !! time of each retrieved point
      logical lsouthpole   !! flag used to have contiguous North Pole data
      logical lwritepar(imxlmb,imxseq)
      logical l_optest
      integer*4 ilimb(imxseq)        !! number of sweeps in scan iseq
                                        !! index of the nomial OM altitudes used
      integer*4 indpar(imxtrgts,imxlmb) !! in the retrieval for each target
      integer*4 itarget              !! number of retrieval target
      integer*4 iseq1                !! index of the first scan to be analyzed
      integer*4 iseq2                !! index of the second scan to be analyzed
      integer*4 itimeini          !! index of  scan where time originates
      integer*4 ngrid                !! number of grid profiles
      integer*4 nret(imxtrgts)  !! number of retrieval altitudes for each target
      real*8 rlatpar(imxlmb,imxseq)  !! latitude of the grid profiles
      real*8 rlongpar(imxlmb,imxseq) !! longitude of the grid profiles
      real*8 difpar(imxlmb,imxseq)   !! Orbital Coordinate of the retrieval-grid points defined by altitude i of grid profile j
      real*8 Ppar(imxlmb,imxseq)     !! Pressure at retrieved points
      real*8 qpar(imxlmb)            !! altitudes of the retrieval-grid profiles
      real*4 rapvcm(imxtop)          !! vector containing the diagonal elements of a-priori VCM
      real*8 rxpar(imxtop)         !! vector containing all the retrieval parameters
      real*4 rvcmns(imxtop)          !! vector containing the diagonal elements of VCM
      real*8 Tpar(imxlmb,imxseq)     !! Temperature at retrieved points
      integer*4 jvetpar              !! counter of the parameter vector (Input/Output)
      integer*4 nflag                !! reported quantity flag (P=1; T=2; VMR>=3) 

!! @INTERNAL_VARIABLES
      integer*4 j
      integer*4 k
      integer*4 l
      integer*4 jlimb
      integer*4 nf1
      character*110 sheader
      character*110 sheader1
      real*8 rq_marco
      real*8 rq_new
      real*8 fact
      real*8 rtemp
      real*4 resd

      resd=0.0e0
      fact=-1.0D0/(dlog(2.0D0))
      nf1=nf+1


      write(nf,'(a)')'Results file'
      write(nf,'(a)') orbit
      if(lsouthpole)write(nf,'(a)')
     &  'Orbit built to have data contiguously acquired at North Pole'
      write(nf,'(a,a11,3x,a11)')'Orbit dates ',date(1),date(2)
      write(nf,'(a,a8)')'First scan recorded at ',rec_time(1,itimeini)
      if(itimeini.eq.iseq1)then
      write(nf,'(a,a8)')'Last scan recorded at ',
     &            rec_time(ilimb(iseq2),iseq2)
      else
      write(nf,'(a,a8)')'Last scan recorded at ',
     &            rec_time(ilimb(itimeini-1),itimeini-1)
      end if
      write(nf,'(//)')
      write(nf,'(a)')'Number of retrieved profiles'
      call skop(nf)
      write(nf,'(i4)') ngrid

      if (nflag.eq.1) then
              sheader = '     LAT     LONG      O.C.     ALT.'//
     &'      P        ESD              Time           Inf.'
             sheader1 = '    (deg)   (deg)     (deg)     (km)'//
     &'    (hPa)     (hPa)             (UTC)         (bit)'
      elseif (nflag.eq.2) then
              sheader = '     LAT     LONG      O.C.     ALT.'//
     &'      T        ESD              Time           Inf.'
             sheader1 = '    (deg)    (deg)    (deg)     (km)'//
     &'     (K)       (K)             (UTC)          (bit)'
      else
              sheader = '     LAT     LONG      O.C.     ALT.'//
     &'     VMR       ESD             Time           Temp     Press'//
     &'     Inf'
             sheader1 = '    (deg)    (deg)    (deg)     (km)'//
     &'    (ppmV)    (ppmV)           (UTC)           (K)     (hPa)'//
     &'    (bit)'
      endif
      if(l_optest)then
       write(nf1,'(a)') orbit
       write(nf1,'(//)')
       write(nf1,'(a)')'Number of retrieved profiles'
       call skop(nf1)
       write(nf1,'(i4)') ngrid
      end if

      do k=1,ngrid
         jlimb=0
         do j=1,nret(itarget)
          l=indpar(itarget,j)
          if(lwritepar(l,k))jlimb=jlimb+1
         end do
         write(nf,'(/,a,i4)')'Profile: ',k
         write(nf,'(a)')'Number of altitudes'
         call skop(nf)
         write(nf,'(i4,/)') jlimb
         write(nf,'(a)') sheader
         write(nf,'(a)') sheader1
         call skop(nf)
         if(l_optest)then
          write(nf1,'(/,a,i4)')'Profile: ',k
          write(nf1,'(a)')'Number of altitudes'
          call skop(nf1)
          write(nf1,'(i4,/)') jlimb
          write(nf1,'(a)') sheader
          write(nf1,'(a)') sheader1
          call skop(nf1)
         end if
         do j=1,nret(itarget)
           l=indpar(itarget,j)
           if(rvcmns(jvetpar).eq.0.0e0.and.l_optest)then
             resd=sqrt(rapvcm(jvetpar))
           else
             resd=sqrt(rvcmns(jvetpar))
           end if
           if(nflag.le.2)then
           if(lwritepar(l,k))then
            if(l_optest)then
             if(rvcmns(jvetpar).ne.0.0e0)then
             rtemp=log(sqrt(rvcmns(jvetpar)/rapvcm(jvetpar)))
             rq_marco=fact*rtemp
             else
              rq_marco=0.0e0
             end if
             write(nf1,'(4f9.2,1pe11.4,2(1pe10.2),0pf10.1)')
     &           rlatpar(l,k),rlongpar(l,k),difpar(l,k),
     &           qpar(l),rxpar(jvetpar),sqrt(rvcmns(jvetpar)),
     &           sqrt(rapvcm(jvetpar)),rq_marco
             write(nf,'(4f9.2,1pe11.4,1pe10.2,a23,0pf6.1)')
     &           rlatpar(l,k),rlongpar(l,k),difpar(l,k),
     &           qpar(l),rxpar(jvetpar),resd,
     &           retr_time(l,k),rq_marco
            else
             write(nf,'(4f9.2,1pe11.4,1pe10.2,a23)')
     &           rlatpar(l,k),rlongpar(l,k),difpar(l,k),
     &           qpar(l),rxpar(jvetpar),resd,
     &           retr_time(l,k)
            end if
           end if
           else
            if(lwritepar(l,k))then
             if(l_optest)then
             if(rvcmns(jvetpar).ne.0.0e0)then
              rtemp=log(sqrt(rvcmns(jvetpar)/rapvcm(jvetpar)))
              rq_marco=fact*rtemp
             else
              rq_marco=0.0e0
             end if
              write(nf1,
     &         '(4f9.2,1pe11.4,2(1pe10.2),0pf10.1)')
     &           rlatpar(l,k),rlongpar(l,k),difpar(l,k),
     &           qpar(l),rxpar(jvetpar),sqrt(rvcmns(jvetpar)),
     &           sqrt(rapvcm(jvetpar)),rq_marco
              write(nf,
     &        '(4f9.2,1pe11.4,1pe10.2,a23,0pf7.2,1pe11.4,0pf6.1)')
     &           rlatpar(l,k),rlongpar(l,k),difpar(l,k),
     &           qpar(l),rxpar(jvetpar),resd,
     &           retr_time(l,k),Tpar(l,k),Ppar(l,k),rq_marco
             else
              write(nf,
     &        '(4f9.2,1pe11.4,1pe10.2,a23,0pf7.2,1pe11.4)')
     &           rlatpar(l,k),rlongpar(l,k),difpar(l,k),
     &           qpar(l),rxpar(jvetpar),resd,
     &           retr_time(l,k),Tpar(l,k),Ppar(l,k)
             end if
            end if
           end if
           jvetpar=jvetpar+1
         end do
      end do
      
      return
      end
