PROGRAM xvik

      USE filtreg_mod, ONLY: inifilr
      USE comconst_mod, ONLY: dtvr,g,r,pi
      USE comvert_mod, ONLY: pa,preff
      
      IMPLICIT NONE
      
      
c=======================================================================
c
c  Pressure at Insight and Viking sites
c
c=======================================================================


c-----------------------------------------------------------------------
c   declarations:
c-----------------------------------------------------------------------


      include "dimensions.h"
      include "paramet.h"
      include "comdissip.h"
      include "comgeom2.h"
      include "netcdf.inc"      


      INTEGER itau,nbpas,nbpasmx 
      PARAMETER(nbpasmx=1000000)
      REAL temps(nbpasmx)
      INTEGER unitlec
      INTEGER i,j,l,jj
      REAL constR

c   Declarations NCDF:
c   -----------------
      CHARACTER*100  varname
      INTEGER ierr,nid,nvarid,dimid
      INTEGER start_ps(3),start_temp(4),start_co2ice(3)
      INTEGER count_ps(3),count_temp(4),count_co2ice(3)

c   declarations pour les points viking:
c   ------------------------------------
      INTEGER isite(3),jsite(3),ifile(3),iv
      
      REAL, PARAMETER ::  lonvik1 = -47.95
      REAL, PARAMETER ::  latvik1 =  22.27
      REAL, PARAMETER ::  lonvik2 =  134.29
      REAL, PARAMETER ::  latvik2 =  47.67
      REAL, PARAMETER ::  loninst =  135.62
      REAL, PARAMETER ::  latinst =  4.502
          
      REAL, PARAMETER :: phivik1 = -3637
      REAL, PARAMETER :: phivik2 = -4505
      REAL, PARAMETER :: phiinst = -2614
      
      
      REAL lonsite(3),latsite(3),phisite(3),phisim(3)
      REAL unanj

c   variables meteo:
c   ----------------
      REAL vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
      REAL t(iip1,jjp1,llm),ps(iip1,jjp1),pstot, phis(iip1,jjp1)
      REAL co2ice(iip1,jjp1), captotN,captotS
      real t7(iip1,jjp1) ! temperature in 7th atmospheric layer

      REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,3),zalpha,zbeta

      LOGICAL firstcal
      INTEGER*4 day0

      REAL ziceco2(iip1,jjp1)
      REAL day,zt,sollong,sol,dayw,dayw_ls
      REAL airtot1,gh

      INTEGER ii,iyear,kyear

      CHARACTER*2 chr2

       
c   declarations de l'interface avec mywrite:
c   -----------------------------------------

      CHARACTER file*80
      CHARACTER pathchmp*80,pathsor*80,nomfich*80
      
      INTEGER Time_unit
      
      REAL ls2sol
      

c   externe:
c   --------

      EXTERNAL iniconst,inigeom,covcont,mywrite
      EXTERNAL exner,pbar
      EXTERNAL coordij,moy2
      EXTERNAL SSUM
      REAL SSUM

c   diurnam:
c   --------  
      integer di,di_prev
      integer k
      real ps_gcm_diurnal, ps_diurnal
      integer compt_diurn

c   harmonics:
c   --------  
    
      integer, parameter :: nbmax = 999999
      integer ndata
      real ls_harm
      real ls_tab(nbmax)
      real ps_diurnal_tab(nbmax),ps_gcm_diurnal_tab(nbmax)
      real a_tab(nbmax),b_tab(nbmax)
      real a_tab_gcm(nbmax),b_tab_gcm(nbmax)
      real a,b
      real ps_harm, ps_gcm_harm
      integer, parameter :: n_harmo = 6
      
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c-----------------------------------------------------------------------
c   initialisations:
c-----------------------------------------------------------------------
      pi=4.*atan(1.)
      pa=20
      preff=610.

      chr2="0"
      iyear=0
      unanj=669.
      print*,'WARNING!!! Assuming',unanj,'sols/year'
      
c-----------------------------------------------------------------------
c   Viking Lander and Insight coordinates:
c   --------------------------------------------------------------------

      lonsite(1) = lonvik1
      latsite(1) = latvik1
      lonsite(2) = lonvik2
      latsite(2) = latvik2
      lonsite(3) = loninst
      latsite(3) = latinst     
      
      phisite(1) = phivik1
      phisite(2) = phivik2
      phisite(3) = phiinst
     
      WRITE(*,*) 'Viking1 coordinates:'
      WRITE(*,*) 'latvik1:',latvik1,' lonvik1:',lonvik1
      WRITE(*,*) 'Phivik1:', phivik1
      
      WRITE(*,*) 'Viking2 coordinates:'
      WRITE(*,*) 'latvik2:',latvik2,' lonvik2:',lonvik2
      WRITE(*,*) 'Phivik2:', phivik2
      
      WRITE(*,*) 'Insight coordinates:'
      WRITE(*,*) 'latinst:',latinst,' loninst:',loninst
      WRITE(*,*) 'Phiinst:', phiinst
                        
      ! convert coordinates to radians
      lonsite(1) = lonvik1 * pi/180.
      latsite(1) = latvik1 * pi/180.
      lonsite(2) = lonvik2 * pi/180.
      latsite(2) = latvik2 * pi/180.
      lonsite(3) = loninst * pi/180.
      latsite(3) = latinst * pi/180.	
      
      
      
      WRITE(*,*) 'Path to the diagfi files directory'
      READ (*,'(a)')  pathchmp
      WRITE(*,*) 'Path to the dir for outputs'
      READ (*,'(a)')  pathsor
      
      WRITE(*,*) 'Output file time axis in sol (1) '//
     &'in ls (2) ,or both (3)'
      READ (*,*)  Time_unit
      
      
      write (*,*)'>>>>>>>>>>>>>>>>', phisite,g
      DO iv=1,3
         phisite(iv)=phisite(iv)*3.73
      END DO

c-----------------------------------------------------------------------
c   output files:
c-----------------------------------------------------------------------
      ifile(1)=12
      ifile(2)=13
      ifile(3)=14
      
      kyear=-1
      unitlec=11
      
      
      print*,'diagfi file name (without trailing .nc)'
      READ(5,'(a)',err=9999) nomfich
      

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c   loop on the diagfi files:
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      firstcal=.true.
      DO WHILE(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50)
      PRINT *,'>>>  nomfich : ',trim(nomfich)

c----------------------------------------------------------------------
c   Open diagfi files :
c----------------------------------------------------------------------

      file=pathchmp(1:len_trim(pathchmp))//'/'//
     s     nomfich(1:len_trim(nomfich))
      PRINT*,'file.nc: ', file(1:len_trim(file))//'.nc'
      PRINT*,'timestep ',dtvr

      ierr= NF_OPEN(file(1:len_trim(file))//'.nc',NF_NOWRITE,nid)        

c----------------------------------------------------------------------
c   initialise physics:
c----------------------------------------------------------------------

      CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR)

      WRITE (*,*) 'day0 = ' , day0

      CALL conf_gcm( 99, .TRUE. )
      CALL iniconst
      CALL inigeom

c----------------------------------------------------------------------
c   Read time :
c----------------------------------------------------------------------


      ierr= NF_INQ_DIMID (nid,"Time",dimid)
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, 'xvik: Le champ <Time> est absent'
          CALL abort
        ENDIF

      ierr= NF_INQ_DIMLEN (nid,dimid,nbpas)

      ierr = NF_INQ_VARID (nid, "Time", nvarid)
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, temps)
#endif
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, 'xvik: Lecture echouee pour <Time>'
          CALL abort
        ENDIF

        PRINT*,'temps(1:10)',(temps(itau),itau=1,10)
        
        
                    
c------------------------------------------------------   
c   Weights for 4 near points at Viking and Insight 
c------------------------------------------------------

      DO iv=1,3
        ! locate index of GCM grid points near VL
	 do i=1,iim
           ! we know longitudes are ordered -180...180
	   write(*,*) i, lonsite(iv),rlonu(i),rlonu(i+1)
           if ((lonsite(iv).ge.rlonu(i)).and.
     &         (lonsite(iv).le.rlonu(i+1))) then
             isite(iv)=i
             exit
           endif
         enddo

         do j=1,jjm-1
           !we know tha latitudes are ordered 90...-90
           if ((latsite(iv).le.rlatv(j)).and.
     &         (latsite(iv).ge.rlatv(j+1))) then
             jsite(iv)=j
             exit
           endif
         enddo
         zalpha=(lonsite(iv)-rlonu(isite(iv)))/
     s          (rlonu(isite(iv)+1)-rlonu(isite(iv)))
         zbeta=(latsite(iv)-rlatv(jsite(iv)))/
     s          (rlatv(jsite(iv)+1)-rlatv(jsite(iv)))
         zw(0,0,iv)=(1.-zalpha)*(1.-zbeta)
         zw(1,0,iv)=zalpha*(1.-zbeta)
         zw(0,1,iv)=(1.-zalpha)*zbeta
         zw(1,1,iv)=zalpha*zbeta
      ENDDO

c----------------------------------------------------------------------
c   true and model altitude at Viking and Insight locations
c----------------------------------------------------------------------


      DO iv=1,3
         phisim(iv)=0.
         DO jj=0,1
            j=jsite(iv)+jj
            DO ii=0,1
               i=isite(iv)+ii
               phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j)
            ENDDO
         ENDDO
      ENDDO
      PRINT*,'phisite at Viking locations for outputs:',phisite
           

c----------------------------------------------------------------------
c   read variables:
c   -------------------------------------------------------------------

       airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1))

c======================================================================
c   Begin the loop on states in inputs files :
c======================================================================


       count_ps=(/iip1,jjp1,1/)
       count_co2ice=(/iip1,jjp1,1/)
       count_temp=(/iip1,jjp1,llm,1/)
       
       DO itau=1,nbpas

       start_ps=(/1,1,itau/)
       start_co2ice=(/1,1,itau/)
       start_temp=(/1,1,1,itau/)
       
c----------------------------------------------------------------------       
c   read fields:
c----------------------------------------------------------------------


ccccccccc  Load Ps ccccccccccccccccccccccccccc


          ierr = NF_INQ_VARID (nid, "ps", nvarid)
#ifdef NC_DOUBLE
          ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start_ps,count_ps, ps)
#else
          ierr = NF_GET_VARA_REAL(nid, nvarid,start_ps,count_ps, ps)
#endif
          IF (ierr.NE.NF_NOERR) THEN
            PRINT*, 'xvik: Lecture echouee pour <ps>'
            CALL abort
          ENDIF
          
          PRINT*,'ps',ps(iip1/2,jjp1/2)

ccccccccc  Load Temperature ccccccccccccccccccccccccccc


          ierr = NF_INQ_VARID (nid, "temp", nvarid)
#ifdef NC_DOUBLE
          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_temp,count_temp, t)
#else
          ierr = NF_GET_VARA_REAL(nid,nvarid,start_temp,count_temp, t)
#endif
          IF (ierr.NE.NF_NOERR) THEN
            PRINT*, 'xvik: Lecture echouee pour <temp>'
            ! Ehouarn: proceed anyways
            ! CALL abort
            write(*,*)'--> Setting temperature to zero !!!'
            t(1:iip1,1:jjp1,1:llm)=0.0
            write(*,*)'--> looking for temp7 (temp in 7th layer)'
            ierr=NF_INQ_VARID(nid,"temp7", nvarid)
            if (ierr.eq.NF_NOERR) then
            write(*,*) "    OK, found temp7 variable"
#ifdef NC_DOUBLE
            ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start_ps,count_ps,t7)
#else
            ierr=NF_GET_VARA_REAL(nid,nvarid,start_ps,count_ps,t7)
#endif
              if (ierr.ne.NF_NOERR) then
                write(*,*)'xvik: failed loading temp7 !'
                stop
              endif
            else ! no 'temp7' variable
              write(*,*)'  No temp7 variable either !'
              write(*,*)'  Will have to to without ...'
              t7(1:iip1,1:jjp1)=0.0
            endif
          ELSE ! t() was successfully loaded, copy 7th layer to t7()
            t7(1:iip1,1:jjp1)=t(1:iip1,1:jjp1,7)
          ENDIF



ccccccccc  Load co2ice ccccccccccccccccccccccccccc


          ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
#ifdef NC_DOUBLE
          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_co2ice,
     &    count_co2ice,  co2ice)
#else
          ierr = NF_GET_VARA_REAL(nid, nvarid,start_co2ice,
     &    count_co2ice, co2ice)
#endif
          IF (ierr.NE.NF_NOERR) THEN
            PRINT*, 'xvik: Lecture echouee pour <co2ice>'
            CALL abort
          ENDIF

c----------------------------------------------------------------------
c Handle calendar
c ---------------------------------------------------------------------

          day=temps(itau)
          PRINT*,'day ',day
          sol=day+day0
          do while (sol.gt.unanj)
            sol=sol-unanj
          enddo
          WRITE (*,*) 'sol: ',sol,' iyear:',iyear

c----------------------------------------------------------------------
c Open /close files
c ---------------------------------------------------------------------

          IF (iyear.NE.kyear) THEN
             WRITE(chr2(1:1),'(i1)') iyear+1
             WRITE (*,*) 'iyear bis:',iyear
             WRITE (*,*) 'chr2:',trim(chr2)
             IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1
             kyear=iyear
             DO ii=1,3
                CLOSE(10+ifile(ii))
                CLOSE(20+ifile(ii))		
                CLOSE(2+ifile(ii))
                CLOSE(4+ifile(ii))
                CLOSE(6+ifile(ii))
                CLOSE(8+ifile(ii))
                CLOSE(16+ifile(ii))
                CLOSE(12+ifile(ii))
                CLOSE(14+ifile(ii))
                CLOSE(97)
                CLOSE(98)
             ENDDO
             CLOSE(5+ifile(1))
             OPEN(ifile(1)+10,file='ps_VL1_year'//chr2,form='formatted')
             OPEN(ifile(2)+10,file='ps_VL2_year'//chr2,form='formatted')
             OPEN(ifile(3)+10,file='ps_INS_year'//chr2,form='formatted')	                                       
             OPEN(97,file='prestot_year'//chr2,form='formatted')
	     

c  Sol or ls or both 
c  Planetary mean surface pressure (Pa)
c  Equivalent pressure of CO2 ice at North Polar cap (Pa)
c  Equivalent pressure of CO2 ice at South Polar cap (Pa) 
c  Total amount of CO2 on the planet (Pa)	    
	     
             IF (Time_unit == 1) THEN
              
	      WRITE(ifile(1)+10,'(a)') '# Sol , Surface Pressure at VL1 at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at VL1 at GCM altitude (Pa)'

              WRITE(ifile(2)+10,'(a)') '# Sol , Surface Pressure at VL2 at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at VL2 at GCM altitude (Pa)'

              WRITE(ifile(3)+10,'(a)') '# Sol , Surface Pressure at Insight at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at Insight at GCM altitude (Pa)'

              WRITE(97,'(a)') '# Sol , Planetary Mean Surface Pressure (Pa) , 
     &        Equivalent pressure of CO2 ice at North Polar cap (Pa) , 
     &        Equivalent pressure of CO2 ice at South Polar cap (Pa) , 
     &        Total amount of CO2 on the planet (Pa)'     
          	      
	     ELSEIF (Time_unit == 2) THEN 

              WRITE(ifile(1)+10,'(a)') '# Ls (deg) , Surface Pressure at VL1 at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at VL1 at GCM altitude (Pa)'

              WRITE(ifile(2)+10,'(a)') '# Ls (deg) , Surface Pressure at VL2 at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at VL2 at GCM altitude (Pa)'

              WRITE(ifile(3)+10,'(a)') '# Ls (deg) , Surface Pressure at Insight at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at Insight at GCM altitude (Pa)'
     
              WRITE(97,'(a)') '# Ls (deg) , Planetary Mean Surface Pressure (Pa) , 
     &        Equivalent pressure of CO2 ice at North Polar cap (Pa) , 
     &        Equivalent pressure of CO2 ice at South Polar cap (Pa) , 
     &        Total amount of CO2 on the planet (Pa)'      
	      
             ELSE 	     

              WRITE(ifile(1)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at VL1 at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at VL1 at GCM altitude (Pa)'

              WRITE(ifile(2)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at VL2 at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at VL2 at GCM altitude (Pa)'

              WRITE(ifile(3)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at Insight at 
     &        true (interpolated) altitude (Pa) ,  
     &        Surface Pressure at Insight at GCM altitude (Pa)'

              WRITE(97,'(a)') '# Sol , Ls (deg) , Planetary Mean Surface Pressure (Pa) , 
     &        Equivalent pressure of CO2 ice at North Polar cap (Pa) , 
     &        Equivalent pressure of CO2 ice at South Polar cap (Pa) , 
     &        Total amount of CO2 on the planet (Pa)' 
	      
	     ENDIF
	     

          ENDIF
 
          dayw = sol
          call sol2ls(sol,sollong)
          dayw_ls = sollong
          
          
          
c----------------------------------------------------------------------
c Compute average planetary pressure
c ---------------------------------------------------------------------


          pstot=0.
          captotS=0.
          captotN=0.
          DO j=1,jjp1
             DO i=1,iim
                pstot=pstot+aire(i,j)*ps(i,j)
             ENDDO
          ENDDO
 
              DO j=1,jjp1/2
                 DO i=1,iim
                    captotN = captotN  +aire(i,j)*co2ice(i,j)
                 ENDDO
              ENDDO
              DO j=jjp1/2+1, jjp1
                 DO i=1,iim
                    captotS = captotS  +aire(i,j)*co2ice(i,j)
                 ENDDO
              ENDDO


c --------------Write output file prestot----------------------- 
c  Sol or ls or both 
c  Planetary mean surface pressure (Pa)
c  Equivalent pressure of CO2 ice at North Polar cap (Pa)
c  Equivalent pressure of CO2 ice at South Polar cap (Pa) 
c  Total amount of CO2 on the planet (Pa)


         IF(Time_unit == 1) THEN
              WRITE(97,'(5e16.6)') dayw,pstot*airtot1
     &       , captotN*g*airtot1, captotS*g*airtot1,       
     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
 
         ELSEIF (Time_unit == 2) THEN    
              WRITE(97,'(5e16.6)') dayw_ls,pstot*airtot1
     &       , captotN*g*airtot1, captotS*g*airtot1, 
     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
     
         ELSE 
             WRITE(97,'(6e16.6)') dayw,dayw_ls,pstot*airtot1
     &       , captotN*g*airtot1,captotS*g*airtot1,
     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
     
                    
         ENDIF           

c----------------------------------------------------------------------
c Loop on sites:
c----------------------------------------------------------------------

c----------------------------------------------------------------------
c interapolate using temperature in the 7th layer, of surface pressure
c----------------------------------------------------------------------

         IF(.NOT.firstcal) THEN
          
          DO iv=1,3

             zp1=0.
             zp2=0.
             zp2_sm=0.
             zt=0.

             DO jj=0,1
             
                j=jsite(iv)+jj
                
                DO ii=0,1
                
                   i=isite(iv)+ii
                   zt=zt+zw(ii,jj,iv)*t7(i,j)
                   zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P)
                   WRITE (*,*) 'ps around iv',ps(i,j),iv

                ENDDO
             ENDDO
             
             zp1=exp(zp1) ! because of the bilinear interpolation in log(P)
             WRITE (*,*) 'constR ',constR 
             WRITE (*,*) 'zt ',zt
             gh=constR*zt            
             
c---------------------------------------------------------------------- 
c  surface pressure extrapolated using temp. from 7th atmospheric layer
c----------------------------------------------------------------------
           
             if (gh.eq.0) then ! if we don't have temperature values
               ! assume a scale height of 10km
               zp2=zp1*exp(-(phisite(iv)-phisim(iv))/(3.73*1.e4))
             else
               zp2=zp1*exp(-(phisite(iv)-phisim(iv))/gh)
             endif
            
          WRITE (*,*) 'iv,pstot,zp2, zp1, phisite(iv),phisim(iv),gh'
          WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phisite(iv),phisim(iv),gh
          WRITE(*,*) "------"
             

c ------Write 3 files (for Vl1, VL2, Insight) --------
c  Sol or ls or both
c  Ps site VLi (i=1,2) or Inisght at true (interpolated) altitude (Pa) (zp2)
c  Ps site VLi (i=1,2) or Insight at GCM altitude (Pa) (zp1)
              
             IF(Time_unit == 1) THEN
             	WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1
             ELSEIF (Time_unit == 2) THEN    
                WRITE(ifile(iv)+10,'(3e15.5)') dayw_ls,zp2,zp1
             ELSE   
                WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1
             ENDIF 
	      
         
          ENDDO

         ENDIF
	          	 	 
         
         firstcal=.false.


c======================================================================
c   End of loop of variables in the diagfi file
c======================================================================

       if (sol.ge.unanj-1.e-5) then 
         ! end of year reached (with some roundoff margin)
         ! increment iyear
         iyear=iyear+1
       endif

      ENDDO

      ierr= NF_CLOSE(nid)

      PRINT*,'End of file',nomfich
      print*,'Entrer new file name (without trailing .nc)', 
     &" or return to end"
      READ(5,'(a)',err=9999) nomfich



c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c   End of loop on the diagfi files
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      ENDDO

      PRINT*,'altitude of VL1',.001*phis(isite(1),jsite(1))/g
      PRINT*,'altitude of VL2',.001*phis(isite(2),jsite(2))/g
      PRINT*,'altitude of Ins',.001*phis(isite(3),jsite(3))/g
      
      DO iv=1,3
         PRINT*,'Site',iv,'   i=',isite(iv),'j  =',jsite(iv)
         WRITE(6,7777)
     s   (rlonv(i)*180./pi,i=isite(iv)-1,isite(iv)+2)
         print*
         DO j=jsite(iv)-1,jsite(iv)+2
            WRITE(6,'(f8.1,10x,5f7.1)')
     s   rlatu(j)*180./pi,(phis(i,j)/(g*1000.),i=isite(iv)-1,
     &   isite(iv)+2)
         ENDDO
         print*
         print*,'zw'
         write(6,'(2(2f10.4/))') ((zw(ii,jj,iv),ii=0,1),jj=0,1)
         print*,'interpolated altitude (km) ',phisim(iv)/1000./g
      ENDDO
      PRINT*,'R=',r
 9999  PRINT*,'End '

7777  FORMAT ('latitude/longitude',4f7.1)

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c   Diurnal Average
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      write(*,*) 'ici'
      DO i=1,kyear+1
       WRITE(chr2(1:1),'(i1)') i
       IF(i.GE.9) WRITE(chr2,'(i2)') i
       DO iv=1,3
        if (iv==1) then
	 
	 open(ifile(iv), file = 'ps_VL1_year'//trim(trim(chr2)))
	 open(ifile(iv)+20, file = 'ps_VL1_year'//trim(chr2)//'_diurnal')
	 IF(Time_unit == 1) THEN	  
	  write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL1 at true
     & (interpolated) altitude (diurnal mean) , PS at VL1 at
     &  GCM altitude (diurnal mean)' 
         ELSEIF (Time_unit == 2) THEN 
	  write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL1 at
     & true (interpolated) altitude (diurnal mean) , PS at VL1
     & at GCM altitude (diurnal mean)' 	 
	 ELSE
	  write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL1
     & at true (interpolated) altitude (diurnal mean) , PS at VL1
     & at GCM altitude (diurnal mean)' 	 
	 ENDIF   
	 	 
	elseif (iv==2) then
	 
	 open(ifile(iv), file = 'ps_VL2_year'//trim(chr2))
	 open(ifile(iv)+20, file = 'ps_VL2_year'//trim(chr2)//'_diurnal')
	 IF(Time_unit == 1) THEN	  
	  write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL2 at true
     & (interpolated) altitude (diurnal mean) , PS at VL1 at
     & GCM altitude (diurnal mean)' 
         ELSEIF (Time_unit == 2) THEN 
	  write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL2
     & at true (interpolated) altitude (diurnal mean) , PS at VL1
     & at GCM altitude (diurnal mean)' 	 
	 ELSE
	  write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL2
     & at true (interpolated) altitude (diurnal mean) , PS at VL1 
     & at GCM altitude (diurnal mean)' 	 
	 ENDIF  
	 	 	 	 
	else 
	 open(ifile(iv), file = 'ps_INS_year'//trim(chr2))
	 open(ifile(iv)+20, file = 'ps_INS_year'//trim(chr2)//'_diurnal')
	 IF(Time_unit == 1) THEN	  
	  write(ifile(iv)+20,'(a)') '# Sol ,  PS at Insight at true
     & (interpolated) altitude (diurnal mean) , PS at VL1 
     & at GCM altitude (diurnal mean)' 
         ELSEIF (Time_unit == 2) THEN 
	  write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at Insight at true 
     & (interpolated) altitude (diurnal mean) , PS at VL1 
     & at GCM altitude (diurnal mean)' 	 
	 ELSE
	  write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at Insight
     & at true (interpolated) altitude (diurnal mean) , PS at VL1 
     & at GCM altitude (diurnal mean)' 	 
	 ENDIF	 	 	 
	endif	
        
	READ(ifile(iv),*)
        IF(Time_unit == 1) THEN
         READ(ifile(iv),*) dayw,zp2,zp1
        ELSEIF (Time_unit == 2) THEN    
         READ(ifile(iv),*) dayw_ls,zp2,zp1
	 dayw = ls2sol(dayw_ls)
        ELSE   
         READ(ifile(iv),*) dayw,dayw_ls,zp2,zp1
        ENDIF 
	                     
	di=floor(dayw)
        di_prev = floor(dayw)
                
        DO k=1,nbmax
         ps_gcm_diurnal = 0.
	 ps_diurnal = 0.
         compt_diurn = 0
                   
         DO WHILE (di==di_prev)
                   
          ps_gcm_diurnal = ps_gcm_diurnal + zp1
	  ps_diurnal = ps_diurnal + zp2
          compt_diurn = compt_diurn + 1
          IF(Time_unit == 1) THEN
           READ(ifile(iv),*,end=777) dayw,zp2,zp1
          ELSEIF (Time_unit == 2) THEN    
           READ(ifile(iv),*,end=777) dayw_ls,zp2,zp1
	   dayw=ls2sol(dayw_ls)
          ELSE   
           READ(ifile(iv),*,end=777) dayw,dayw_ls,zp2,zp1
          ENDIF 
	  di=floor(dayw)
                                                
         ENDDO
                   	
         ps_gcm_diurnal = ps_gcm_diurnal/compt_diurn
         ps_diurnal = ps_diurnal/compt_diurn
 	 
	 IF(Time_unit == 1) THEN	  
                  
          write(ifile(iv)+20,'(i4,2e15.5)') di_prev, ps_diurnal, 
     &         ps_gcm_diurnal
          
         ELSEIF (Time_unit == 2) THEN 
	 
          write(ifile(iv)+20,'(3e15.5)') dayw_ls, ps_diurnal, 
     &         ps_gcm_diurnal	 
	  
	 ELSE
	 
          write(ifile(iv)+20,'(i4,3e15.5)') di_prev, dayw_ls, 
     &         ps_diurnal, ps_gcm_diurnal	 
	 
	 ENDIF
	  
                 
         di_prev = di
                   
        ENDDO
	close(ifile(iv)+20)
	close(ifile(iv))
777    ENDDO 
      ENDDO



c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c   Harmonics
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        
      DO i=1,kyear+1
       WRITE(chr2(1:1),'(i1)') i
       IF(i.GE.9) WRITE(chr2,'(i2)') i      
       DO iv=1,3
       
        if (iv==1) then
	 open(ifile(iv), file = 'ps_VL1_year'//trim(chr2)//'_diurnal')
	 open(ifile(iv)+20, file = 'ps_VL1_year'//trim(chr2)//'_harmonics')
	
	 IF(Time_unit == 1) THEN	  
	  write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL1 at true
     & (interpolated) altitude (harmonics fit) , PS at VL1 at
     &  GCM altitude (harmonics fit)' 
         ELSEIF (Time_unit == 2) THEN 
	  write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL1 at
     & true (interpolated) altitude (harmonics fit) , PS at VL1
     & at GCM altitude (harmonics fit)' 	 
	 ELSE
	  write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL1
     & at true (interpolated) altitude (harmonics fit) , PS at VL1
     & at GCM altitude (harmonics fit)' 	 
	 ENDIF  	 	 
	
	elseif (iv==2) then
	 open(ifile(iv), file = 'ps_VL2_year'//trim(chr2)//'_diurnal')
	 open(ifile(iv)+20, file = 'ps_VL2_year'//trim(chr2)//'_harmonics')
	 
	 IF(Time_unit == 1) THEN	  
	  write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL2 at true
     & (interpolated) altitude (harmonics fit) , PS at VL2 at
     &  GCM altitude (harmonics fit)' 
         ELSEIF (Time_unit == 2) THEN 
	  write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL2 at
     & true (interpolated) altitude (harmonics fit) , PS at VL2
     & at GCM altitude (harmonics fit)' 	 
	 ELSE
	  write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL2
     & at true (interpolated) altitude (harmonics fit) , PS at VL2
     & at GCM altitude (harmonics fit)' 	 
	 ENDIF	 	 	 
	
	else 
	 open(ifile(iv), file = 'ps_INS_year'//trim(chr2)//'_diurnal')
	 open(ifile(iv)+20, file = 'ps_INS_year'//trim(chr2)//'_harmonics')
	 
	 IF(Time_unit == 1) THEN	  
	  write(ifile(iv)+20,'(a)') '# Sol ,  PS at Insight at true
     & (interpolated) altitude (harmonics fit) , PS at Insight at
     &  GCM altitude (harmonics fit)' 
         ELSEIF (Time_unit == 2) THEN 
	  write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at Insight at
     & true (interpolated) altitude (harmonics fit) , PS at Insight
     & at GCM altitude (harmonics fit)' 	 
	 ELSE
	  write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at Insight
     & at true (interpolated) altitude (harmonics fit) , PS at Insight
     & at GCM altitude (harmonics fit)' 	 
	 ENDIF	 
	 	 	 
	endif

 	READ(ifile(iv),*)
        
	DO k = 1,nbmax
	 IF (Time_unit == 1) THEN  
          READ(ifile(iv),*,end=99) di_prev, ps_diurnal_tab(k), 
     &         ps_gcm_diurnal_tab(k)
          call sol2ls(real(di_prev),ls_tab(k)) 
         ELSEIF (Time_unit == 2) THEN 
          READ(ifile(iv),*,end=99) ls_tab(k), ps_diurnal_tab(k), 
     &         ps_gcm_diurnal_tab(k)	 
	 ELSE  
          READ(ifile(iv),*,end=99) di_prev, ls_tab(k), 
     &    ps_diurnal_tab(k), ps_gcm_diurnal_tab(k)
         ENDIF
                  
        ENDDO
	
        
99      ndata=k-1
        
	if(ls_tab(ndata).gt.350.) then     
  
        do k = 0,n_harmo
        
           call DiscreetFourierHn(ndata,ls_tab,ps_diurnal_tab,k,a,b)
           if(modulo(k,2)==0) then
             a_tab(k+1)= a
             b_tab(k+1)= b
           else 
             a_tab(k+1)= -a
             b_tab(k+1)= -b
           endif                       
          
	   call DiscreetFourierHn(ndata,ls_tab,ps_gcm_diurnal_tab,k,a,b)
           if(modulo(k,2)==0) then
             a_tab_gcm(k+1)= a
             b_tab_gcm(k+1)= b
           else 
             a_tab_gcm(k+1)= -a
             b_tab_gcm(k+1)= -b
           endif
	            
	enddo 
         
         write(ifile(iv)+20,'(a)') 'Fourrier coefficients ak/bk 
     &(interpolated altitude)'
         write(ifile(iv)+20,'(a,7f)') '#',(a_tab(k),k=1,n_harmo+1)
         write(ifile(iv)+20,'(a,7f)') '#',(b_tab(k),k=1,n_harmo+1)
         write(ifile(iv)+20,'(a)') 'Fourrier coefficients ak/bk 
     &(GCM altitude)'
         write(ifile(iv)+20,'(a,7f)') '#',(a_tab_gcm(k),k=1,n_harmo+1)
         write(ifile(iv)+20,'(a,7f)') '#',(b_tab_gcm(k),k=1,n_harmo+1)	 
	 
	
        do l = 1,669
        
         call sol2ls(real(l),ls_harm)
         ps_harm = a_tab(1)
         ps_gcm_harm = a_tab_gcm(1)
	
         do k = 1,n_harmo
        
            ps_harm = ps_harm + a_tab(k+1)*
     &cos((pi/180.)*k*(ls_harm)) 
     &+ b_tab(k+1)*sin((pi/180.)*k*(ls_harm))
 
            ps_gcm_harm = ps_gcm_harm + a_tab_gcm(k+1)*
     &cos((pi/180.)*k*(ls_harm)) 
     &+ b_tab_gcm(k+1)*sin((pi/180.)*k*(ls_harm))    
        
         enddo


	 IF(Time_unit == 1) THEN	  
                  
         write(ifile(iv)+20,'(i4,2e15.5)') l, ps_harm, 
     &         ps_gcm_harm 
          
         ELSEIF (Time_unit == 2) THEN 
	 
         write(ifile(iv)+20,'(3e15.5)') ls_harm, ps_harm, 
     &         ps_gcm_harm	 
	  
	 ELSE
	 
         write(ifile(iv)+20,'(i4,3e15.5)') l,ls_harm, ps_harm, 
     &         ps_gcm_harm	 
	 
	 ENDIF
                       
        enddo
        close(ifile(iv))
	close(ifile(iv)+20)
	
	else  
         write(ifile(iv)+20,'(a)') 'not computed because
     & year is not complete'	
	endif ! (ls.gt.350)
	            
       ENDDO
      ENDDO       
      
      
      
      END

      subroutine sol2ls(sol,Ls)
!==============================================================================
! Purpose: 
! Convert a date/time, given in sol (martian day),
! into solar longitude date/time, in Ls (in degrees),
! where sol=0 is (by definition) the northern hemisphere
!  spring equinox (where Ls=0).
!==============================================================================
! Notes:
! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
! "Ls" will be increased by N*360
! Won't work as expected if sol is negative (then again,
! why would that ever happen?)
!==============================================================================

      implicit none

!==============================================================================
! Arguments:
!==============================================================================
      real,intent(in) :: sol
      real,intent(out) :: Ls

!==============================================================================
! Local variables:
!==============================================================================
      real year_day,peri_day,timeperi,e_elips,twopi,degrad
      data year_day /669./            ! # of sols in a martian year
      data peri_day /485.0/           
      data timeperi /1.9082314/ 
      data e_elips  /0.093358/
      data twopi       /6.2831853/    ! 2.*pi
      data degrad   /57.2957795/      ! pi/180

      real zanom,xref,zx0,zdx,zteta,zz

      integer count_years
      integer iter

!==============================================================================
! 1. Compute Ls
!==============================================================================

      zz=(sol-peri_day)/year_day
      zanom=twopi*(zz-nint(zz))
      xref=abs(zanom)

!  The equation zx0 - e * sin (zx0) = xref, solved by Newton
      zx0=xref+e_elips*sin(xref)
      do iter=1,20 ! typically, 2 or 3 iterations are enough
         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
         zx0=zx0+zdx
         if(abs(zdx).le.(1.e-7)) then
!            write(*,*)'iter:',iter,'     |zdx|:',abs(zdx)
             exit
         endif 
      enddo

      if(zanom.lt.0.) zx0=-zx0

      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
      Ls=zteta-timeperi

      if(Ls.lt.0.) then
         Ls=Ls+twopi
      else
         if(Ls.gt.twopi) then
            Ls=Ls-twopi
         endif
      endif

      Ls=degrad*Ls
! Ls is now in degrees

!==============================================================================
! 1. Account for (eventual) years included in input date/time sol
!==============================================================================

      count_years=0 ! initialize
      zz=sol  ! use "zz" to store (and work on) the value of sol
      do while (zz.ge.year_day)
          count_years=count_years+1
          zz=zz-year_day
      enddo

! Add 360 degrees to Ls for every year
      if (count_years.ne.0) then
         Ls=Ls+360.*count_years
      endif


      end subroutine sol2ls
      
      
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real function ls2sol(ls)

!  Returns solar longitude, Ls (in deg.), from day number (in sol),
!  where sol=0=Ls=0 at the northern hemisphere spring equinox

      implicit none

!  Arguments:
      real, intent(in) :: ls

!  Local:
      double precision xref,zx0,zteta,zz
!        xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
      double precision year_day 
      double precision peri_day,timeperi,e_elips
      double precision pi,degrad 
      parameter (year_day=668.6d0) ! number of sols in a amartian year
!      data peri_day /485.0/
      parameter (peri_day=485.35d0) ! date (in sols) of perihelion
!  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
      parameter (timeperi=1.90258341759902d0)
      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
      parameter (pi=3.14159265358979d0)
      parameter (degrad=57.2957795130823d0)

      if (abs(ls).lt.1.0e-5) then
         if (ls.ge.0.0) then
            ls2sol = 0.0
         else
            ls2sol = real(year_day)
         end if
         return
      end if

      zteta = ls/degrad + timeperi
      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
      xref = zx0-e_elips*dsin(zx0)
      zz = xref/(2.*pi)
      ls2sol = real(zz*year_day + peri_day)
      if (ls2sol.lt.0.0) ls2sol = ls2sol + real(year_day)
      if (ls2sol.ge.year_day) ls2sol = ls2sol - real(year_day)

      return
      end

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


!****************************************************************
!*   Calculate the Fourier harmonic #n of a periodic discreet   *
!*   function F(x) defined by ndata points.                     *
!* ------------------------------------------------------------ *
!* Inputs:                                                      *
!*            ndata: number of points of discreet function.     *
!*            X    : pointer to table storing xi abscissas.     *
!*            Y    : pointer to table storing yi ordinates.     *
!*                                                              *
!* Outputs:                                                     *
!*            a    : coefficient an of the Fourier series.      *
!*            b:   : coefficient bn of the Fourier series.      *
!* ------------------------------------------------------------ *
!* Reference: "Mathematiques en Turbo-Pascal By Marc Ducamp and *
!*             Alain Reverchon, 1. Analyse, Editions Eyrolles,  *
!*             Paris, 1991" [BIBLI 03].                         *
!*                                                              *
!*                           F90 Version By J-P Moreau, Paris.  *
!*                                  (www.jpmoreau.fr)           *
!****************************************************************
! Note: The Fourier series  of a periodic discreet function F(x)
!       can be written under the form:
!                   n=inf. 
!       F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x)  
!                   n=1 
! ----------------------------------------------------------------
         
        Subroutine DiscreetFourierHn(ndata, X, Y, n, a, b)
         real X(1:ndata), Y(1:ndata), a, b
         integer ndata,n,i
         real xi,T,om,wa,wb,wc,wd,wg,wh,wi,wl,wm,wn,wp
         real PI 
         PI=4.d0*datan(1.d0)
         T=X(ndata)-X(1); xi=(X(ndata)+X(1))/2.d0
         om = 2*PI*n/T; a=0.d0; b=0.d0
         do i=1, ndata-1
             wa=X(i); wb=X(i+1)
             wc=Y(i); wd=Y(i+1)
          if (wa.ne.wb) then
             wg = (wd-wc)/(wb-wa)
             wh = om*(wa-xi); wi=om*(wb-xi)
           if (n==0) then
             a = a + (wb-wa)*(wc+wg/2.d0*(wb-wa))
           else 
             wl = cos(wh); wm = sin(wh)
             wn = cos(wi); wp = sin(wi)
             a = a + wg/om*(wn-wl) + wd*wp - wc*wm
             b = b + wg/om*(wp-wm) - wd*wn + wc*wl
           end if
         end if
        end do
        a = a/T; b = b/T
        if (n.ne.0) then
           a = a*2.d0/om; b = b*2.d0/om
        end if
        return
       End