Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F	(revision 2832)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F	(revision 2844)
@@ -10,5 +10,5 @@
 c=======================================================================
 c
-c  Pression au site Viking
+c  Pressure at Insight and Viking sites
 c
 c=======================================================================
@@ -43,5 +43,5 @@
 c   declarations pour les points viking:
 c   ------------------------------------
-      INTEGER ivik(2),jvik(2),ifile(2),iv
+      INTEGER isite(3),jsite(3),ifile(3),iv
       
       REAL, PARAMETER ::  lonvik1 = -47.95
@@ -49,10 +49,13 @@
       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 lonvik(2),latvik(2),phivik(2),phisim(2)
+      REAL, PARAMETER :: phiinst = -2614
+      
+      
+      REAL lonsite(3),latsite(3),phisite(3),phisim(3)
       REAL unanj
 
@@ -64,5 +67,5 @@
       real t7(iip1,jjp1) ! temperature in 7th atmospheric layer
 
-      REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta
+      REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,3),zalpha,zbeta
 
       LOGICAL firstcal
@@ -85,4 +88,6 @@
       
       INTEGER Time_unit
+      
+      REAL ls2sol
       
 
@@ -95,6 +100,27 @@
       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
 
@@ -112,26 +138,40 @@
       
 c-----------------------------------------------------------------------
-c   Viking Lander coordinates:
+c   Viking Lander and Insight coordinates:
 c   --------------------------------------------------------------------
 
-      lonvik(1) = lonvik1
-      latvik(1) = latvik1
-      lonvik(2) = lonvik2
-      latvik(2) = latvik2
-      
-      phivik(1) = phivik1
-      phivik(2) = phivik2
-      
-      WRITE(*,*) 'Viking coordinates:'
-      WRITE(*,*) 'latvik:',latvik,' lonvik:',lonvik
-      WRITE(*,*) 'Phivik:', phivik
-      
+      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
-      lonvik(1) = lonvik1 * pi/180.
-      latvik(1) = latvik1 * pi/180.
-      lonvik(2) = lonvik2 * pi/180.
-      latvik(2) = latvik2 * pi/180.
-      
-
+      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
@@ -144,7 +184,7 @@
       
       
-      write (*,*)'>>>>>>>>>>>>>>>>', phivik,g
-      DO iv=1,2
-         phivik(iv)=phivik(iv)*3.73
+      write (*,*)'>>>>>>>>>>>>>>>>', phisite,g
+      DO iv=1,3
+         phisite(iv)=phisite(iv)*3.73
       END DO
 
@@ -154,4 +194,6 @@
       ifile(1)=12
       ifile(2)=13
+      ifile(3)=14
+      
       kyear=-1
       unitlec=11
@@ -171,5 +213,5 @@
 
 c----------------------------------------------------------------------
-c   Ouverture des fichiers histoire:
+c   Open diagfi files :
 c----------------------------------------------------------------------
 
@@ -194,5 +236,5 @@
 
 c----------------------------------------------------------------------
-c   Lecture temps :
+c   Read time :
 c----------------------------------------------------------------------
 
@@ -221,31 +263,32 @@
         
                     
-c----------------------------------------------------------------------   
-c   ponderations pour les 4 points autour de Viking
-c----------------------------------------------------------------------
-
-
-      DO iv=1,2
+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
+	 do i=1,iim
            ! we know longitudes are ordered -180...180
-           if ((lonvik(iv).ge.rlonu(i)).and.
-     &         (lonvik(iv).le.rlonu(i+1))) then
-             ivik(iv)=i
+	   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 ((latvik(iv).le.rlatv(j)).and.
-     &         (latvik(iv).ge.rlatv(j+1))) then
-             jvik(iv)=j
+           if ((latsite(iv).le.rlatv(j)).and.
+     &         (latsite(iv).ge.rlatv(j+1))) then
+             jsite(iv)=j
              exit
            endif
          enddo
-         zalpha=(lonvik(iv)-rlonu(ivik(iv)))/
-     s          (rlonu(ivik(iv)+1)-rlonu(ivik(iv)))
-         zbeta=(latvik(iv)-rlatv(jvik(iv)))/
-     s          (rlatv(jvik(iv)+1)-rlatv(jvik(iv)))
+         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)
@@ -255,19 +298,19 @@
 
 c----------------------------------------------------------------------
-c   true and model altitude at Viking locations
-c----------------------------------------------------------------------
-
-
-      DO iv=1,2
+c   true and model altitude at Viking and Insight locations
+c----------------------------------------------------------------------
+
+
+      DO iv=1,3
          phisim(iv)=0.
          DO jj=0,1
-            j=jvik(iv)+jj
+            j=jsite(iv)+jj
             DO ii=0,1
-               i=ivik(iv)+ii
+               i=isite(iv)+ii
                phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j)
             ENDDO
          ENDDO
       ENDDO
-      PRINT*,'phivik at Viking locations for outputs:',phivik
+      PRINT*,'phisite at Viking locations for outputs:',phisite
            
 
@@ -279,5 +322,5 @@
 
 c======================================================================
-c   debut de la boucle sur les etats dans un fichier histoire:
+c   Begin the loop on states in inputs files :
 c======================================================================
 
@@ -391,6 +434,7 @@
              IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1
              kyear=iyear
-             DO ii=1,2
+             DO ii=1,3
                 CLOSE(10+ifile(ii))
+                CLOSE(20+ifile(ii))		
                 CLOSE(2+ifile(ii))
                 CLOSE(4+ifile(ii))
@@ -404,7 +448,75 @@
              ENDDO
              CLOSE(5+ifile(1))
-             OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted')
-             OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted')                                  
-             OPEN(97,file='xprestot'//chr2,form='formatted')
+             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
@@ -442,22 +554,26 @@
 
 
-c --------------Write output file xprestot----------------------- 
-c  Sol ou ls ou les deux 
-c  Ps_moy_planetaire (Pa)
-c  Pequivalente de glace de CO2 au Nord (si entierement sublimee) (Pa)
-c  Pequivalente de glace de CO2 au Sud (si entierement sublimee) (Pa) 
+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,'(4e16.6)') dayw,pstot*airtot1
-     &       , captotN*g*airtot1, captotS*g*airtot1       
+              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,'(4e16.6)') dayw_ls,pstot*airtot1
-     &       , captotN*g*airtot1, captotS*g*airtot1
+              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,'(5e16.6)') dayw,dayw_ls,pstot*airtot1
-     &       , captotN*g*airtot1,captotS*g*airtot1
+             WRITE(97,'(6e16.6)') dayw,dayw_ls,pstot*airtot1
+     &       , captotN*g*airtot1,captotS*g*airtot1,
+     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
      
                     
@@ -465,5 +581,5 @@
 
 c----------------------------------------------------------------------
-c Loop on Viking sites:
+c Loop on sites:
 c----------------------------------------------------------------------
 
@@ -474,5 +590,5 @@
          IF(.NOT.firstcal) THEN
           
-          DO iv=1,2
+          DO iv=1,3
 
              zp1=0.
@@ -483,9 +599,9 @@
              DO jj=0,1
              
-                j=jvik(iv)+jj
+                j=jsite(iv)+jj
                 
                 DO ii=0,1
                 
-                   i=ivik(iv)+ii
+                   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)
@@ -506,18 +622,18 @@
              if (gh.eq.0) then ! if we don't have temperature values
                ! assume a scale height of 10km
-               zp2=zp1*exp(-(phivik(iv)-phisim(iv))/(3.73*1.e4))
+               zp2=zp1*exp(-(phisite(iv)-phisim(iv))/(3.73*1.e4))
              else
-               zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh)
+               zp2=zp1*exp(-(phisite(iv)-phisim(iv))/gh)
              endif
             
-          WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh'
-          WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh
+          WRITE (*,*) 'iv,pstot,zp2, zp1, phisite(iv),phisim(iv),gh'
+          WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phisite(iv),phisim(iv),gh
           WRITE(*,*) "------"
              
 
-c ------Write 2 files (1 for Vl1, 1 for VL2) xpsol ------
+c ------Write 3 files (for Vl1, VL2, Insight) --------
 c  Sol or ls or both
-c  Ps site VLi (i=1,2) at GCM altitude (Pa)
-c  Ps site VLi (i=1,2) at true (interpolated) altitude (Pa)
+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
@@ -527,9 +643,11 @@
              ELSE   
                 WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1
-             ENDIF     
-              
+             ENDIF 
+	      
+         
           ENDDO
 
          ENDIF
+	          	 	 
          
          firstcal=.false.
@@ -556,4 +674,5 @@
 
 
+
 c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 c   End of loop on the diagfi files
@@ -562,14 +681,17 @@
       ENDDO
 
-      PRINT*,'altitude of VL1',.001*phis(ivik(1),jvik(1))/g
-      PRINT*,'altitude of VL2',.001*phis(ivik(2),jvik(2))/g
-      DO iv=1,2
-         PRINT*,'Viking',iv,'   i=',ivik(iv),'j  =',jvik(iv)
+      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=ivik(iv)-1,ivik(iv)+2)
+     s   (rlonv(i)*180./pi,i=isite(iv)-1,isite(iv)+2)
          print*
-         DO j=jvik(iv)-1,jvik(iv)+2
+         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=ivik(iv)-1,ivik(iv)+2)
+     s   rlatu(j)*180./pi,(phis(i,j)/(g*1000.),i=isite(iv)-1,
+     &   isite(iv)+2)
          ENDDO
          print*
@@ -583,6 +705,303 @@
 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
 
@@ -679,2 +1098,110 @@
 
       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     
