Changeset 2844 for trunk/LMDZ.MARS


Ignore:
Timestamp:
Dec 18, 2022, 11:44:43 PM (2 years ago)
Author:
tpierron
Message:

Mars GCM :
Update xvik program. The program now interpolates surface pressure at 3 locations : Viking Lander 1, Viking Lander 2 and Insight.
Also add as outputs :

  • Diurnal average
  • Harmonics fit (Fourier series) of the diurnal average (6 harmonics)

Columns of each output are described at the top of each file. (Fourier coefficients are also given for harmonics fit)
Total CO2 inventory is also caluclated now by xvik.
TP

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2840 r2844  
    38353835sublimation / condensation in W.m-2
    38363836Also corrected units of shortwave & longwave heat rates in writediagfi
     3837
     3838== 18/12/2022 == TP
     3839Update xvik program. The program now interpolates surface pressure at 3 locations : Viking Lander 1, Viking Lander 2 and Insight.
     3840Also add as outputs :
     3841- Diurnal average
     3842- Harmonics fit (Fourier series) of the diurnal average (6 harmonics)
     3843Columns of each output are described at the top of each file. (Fourier coefficients are also given for harmonics fit)
     3844Total CO2 inventory is also caluclated now by xvik.

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F

    r2824 r2844  
    1010c=======================================================================
    1111c
    12 c  Pression au site Viking
     12c  Pressure at Insight and Viking sites
    1313c
    1414c=======================================================================
    … …  
    4343c   declarations pour les points viking:
    4444c   ------------------------------------
    45       INTEGER ivik(2),jvik(2),ifile(2),iv
     45      INTEGER isite(3),jsite(3),ifile(3),iv
    4646     
    4747      REAL, PARAMETER ::  lonvik1 = -47.95
    … …  
    4949      REAL, PARAMETER ::  lonvik2 =  134.29
    5050      REAL, PARAMETER ::  latvik2 =  47.67
    51      
     51      REAL, PARAMETER ::  loninst =  135.62
     52      REAL, PARAMETER ::  latinst =  4.502
     53         
    5254      REAL, PARAMETER :: phivik1 = -3637
    5355      REAL, PARAMETER :: phivik2 = -4505
    54      
    55      
    56       REAL lonvik(2),latvik(2),phivik(2),phisim(2)
     56      REAL, PARAMETER :: phiinst = -2614
     57     
     58     
     59      REAL lonsite(3),latsite(3),phisite(3),phisim(3)
    5760      REAL unanj
    5861
    … …  
    6467      real t7(iip1,jjp1) ! temperature in 7th atmospheric layer
    6568
    66       REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta
     69      REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,3),zalpha,zbeta
    6770
    6871      LOGICAL firstcal
    … …  
    8588     
    8689      INTEGER Time_unit
     90     
     91      REAL ls2sol
    8792     
    8893
    … …  
    95100      EXTERNAL SSUM
    96101      REAL SSUM
    97      
    98 
     102
     103c   diurnam:
     104c   -------- 
     105      integer di,di_prev
     106      integer k
     107      real ps_gcm_diurnal, ps_diurnal
     108      integer compt_diurn
     109
     110c   harmonics:
     111c   -------- 
     112   
     113      integer, parameter :: nbmax = 999999
     114      integer ndata
     115      real ls_harm
     116      real ls_tab(nbmax)
     117      real ps_diurnal_tab(nbmax),ps_gcm_diurnal_tab(nbmax)
     118      real a_tab(nbmax),b_tab(nbmax)
     119      real a_tab_gcm(nbmax),b_tab_gcm(nbmax)
     120      real a,b
     121      real ps_harm, ps_gcm_harm
     122      integer, parameter :: n_harmo = 6
     123     
     124     
    99125cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    100126
    … …  
    112138     
    113139c-----------------------------------------------------------------------
    114 c   Viking Lander coordinates:
     140c   Viking Lander and Insight coordinates:
    115141c   --------------------------------------------------------------------
    116142
    117       lonvik(1) = lonvik1
    118       latvik(1) = latvik1
    119       lonvik(2) = lonvik2
    120       latvik(2) = latvik2
    121      
    122       phivik(1) = phivik1
    123       phivik(2) = phivik2
    124      
    125       WRITE(*,*) 'Viking coordinates:'
    126       WRITE(*,*) 'latvik:',latvik,' lonvik:',lonvik
    127       WRITE(*,*) 'Phivik:', phivik
    128      
     143      lonsite(1) = lonvik1
     144      latsite(1) = latvik1
     145      lonsite(2) = lonvik2
     146      latsite(2) = latvik2
     147      lonsite(3) = loninst
     148      latsite(3) = latinst     
     149     
     150      phisite(1) = phivik1
     151      phisite(2) = phivik2
     152      phisite(3) = phiinst
     153     
     154      WRITE(*,*) 'Viking1 coordinates:'
     155      WRITE(*,*) 'latvik1:',latvik1,' lonvik1:',lonvik1
     156      WRITE(*,*) 'Phivik1:', phivik1
     157     
     158      WRITE(*,*) 'Viking2 coordinates:'
     159      WRITE(*,*) 'latvik2:',latvik2,' lonvik2:',lonvik2
     160      WRITE(*,*) 'Phivik2:', phivik2
     161     
     162      WRITE(*,*) 'Insight coordinates:'
     163      WRITE(*,*) 'latinst:',latinst,' loninst:',loninst
     164      WRITE(*,*) 'Phiinst:', phiinst
     165                       
    129166      ! convert coordinates to radians
    130       lonvik(1) = lonvik1 * pi/180.
    131       latvik(1) = latvik1 * pi/180.
    132       lonvik(2) = lonvik2 * pi/180.
    133       latvik(2) = latvik2 * pi/180.
    134      
    135 
     167      lonsite(1) = lonvik1 * pi/180.
     168      latsite(1) = latvik1 * pi/180.
     169      lonsite(2) = lonvik2 * pi/180.
     170      latsite(2) = latvik2 * pi/180.
     171      lonsite(3) = loninst * pi/180.
     172      latsite(3) = latinst * pi/180.   
     173     
     174     
     175     
    136176      WRITE(*,*) 'Path to the diagfi files directory'
    137177      READ (*,'(a)')  pathchmp
    … …  
    144184     
    145185     
    146       write (*,*)'>>>>>>>>>>>>>>>>', phivik,g
    147       DO iv=1,2
    148          phivik(iv)=phivik(iv)*3.73
     186      write (*,*)'>>>>>>>>>>>>>>>>', phisite,g
     187      DO iv=1,3
     188         phisite(iv)=phisite(iv)*3.73
    149189      END DO
    150190
    … …  
    154194      ifile(1)=12
    155195      ifile(2)=13
     196      ifile(3)=14
     197     
    156198      kyear=-1
    157199      unitlec=11
    … …  
    171213
    172214c----------------------------------------------------------------------
    173 c   Ouverture des fichiers histoire:
     215c   Open diagfi files :
    174216c----------------------------------------------------------------------
    175217
    … …  
    194236
    195237c----------------------------------------------------------------------
    196 c   Lecture temps :
     238c   Read time :
    197239c----------------------------------------------------------------------
    198240
    … …  
    221263       
    222264                   
    223 c----------------------------------------------------------------------   
    224 c   ponderations pour les 4 points autour de Viking
    225 c----------------------------------------------------------------------
    226 
    227 
    228       DO iv=1,2
     265c------------------------------------------------------   
     266c   Weights for 4 near points at Viking and Insight
     267c------------------------------------------------------
     268
     269      DO iv=1,3
    229270        ! locate index of GCM grid points near VL
    230         do i=1,iim
     271        do i=1,iim
    231272           ! we know longitudes are ordered -180...180
    232            if ((lonvik(iv).ge.rlonu(i)).and.
    233      &         (lonvik(iv).le.rlonu(i+1))) then
    234              ivik(iv)=i
     273           write(*,*) i, lonsite(iv),rlonu(i),rlonu(i+1)
     274           if ((lonsite(iv).ge.rlonu(i)).and.
     275     &         (lonsite(iv).le.rlonu(i+1))) then
     276             isite(iv)=i
    235277             exit
    236278           endif
    237279         enddo
     280
    238281         do j=1,jjm-1
    239282           !we know tha latitudes are ordered 90...-90
    240            if ((latvik(iv).le.rlatv(j)).and.
    241      &         (latvik(iv).ge.rlatv(j+1))) then
    242              jvik(iv)=j
     283           if ((latsite(iv).le.rlatv(j)).and.
     284     &         (latsite(iv).ge.rlatv(j+1))) then
     285             jsite(iv)=j
    243286             exit
    244287           endif
    245288         enddo
    246          zalpha=(lonvik(iv)-rlonu(ivik(iv)))/
    247      s          (rlonu(ivik(iv)+1)-rlonu(ivik(iv)))
    248          zbeta=(latvik(iv)-rlatv(jvik(iv)))/
    249      s          (rlatv(jvik(iv)+1)-rlatv(jvik(iv)))
     289         zalpha=(lonsite(iv)-rlonu(isite(iv)))/
     290     s          (rlonu(isite(iv)+1)-rlonu(isite(iv)))
     291         zbeta=(latsite(iv)-rlatv(jsite(iv)))/
     292     s          (rlatv(jsite(iv)+1)-rlatv(jsite(iv)))
    250293         zw(0,0,iv)=(1.-zalpha)*(1.-zbeta)
    251294         zw(1,0,iv)=zalpha*(1.-zbeta)
    … …  
    255298
    256299c----------------------------------------------------------------------
    257 c   true and model altitude at Viking locations
    258 c----------------------------------------------------------------------
    259 
    260 
    261       DO iv=1,2
     300c   true and model altitude at Viking and Insight locations
     301c----------------------------------------------------------------------
     302
     303
     304      DO iv=1,3
    262305         phisim(iv)=0.
    263306         DO jj=0,1
    264             j=jvik(iv)+jj
     307            j=jsite(iv)+jj
    265308            DO ii=0,1
    266                i=ivik(iv)+ii
     309               i=isite(iv)+ii
    267310               phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j)
    268311            ENDDO
    269312         ENDDO
    270313      ENDDO
    271       PRINT*,'phivik at Viking locations for outputs:',phivik
     314      PRINT*,'phisite at Viking locations for outputs:',phisite
    272315           
    273316
    … …  
    279322
    280323c======================================================================
    281 c   debut de la boucle sur les etats dans un fichier histoire:
     324c   Begin the loop on states in inputs files :
    282325c======================================================================
    283326
    … …  
    391434             IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1
    392435             kyear=iyear
    393              DO ii=1,2
     436             DO ii=1,3
    394437                CLOSE(10+ifile(ii))
     438                CLOSE(20+ifile(ii))             
    395439                CLOSE(2+ifile(ii))
    396440                CLOSE(4+ifile(ii))
    … …  
    404448             ENDDO
    405449             CLOSE(5+ifile(1))
    406              OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted')
    407              OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted')                                 
    408              OPEN(97,file='xprestot'//chr2,form='formatted')
     450             OPEN(ifile(1)+10,file='ps_VL1_year'//chr2,form='formatted')
     451             OPEN(ifile(2)+10,file='ps_VL2_year'//chr2,form='formatted')
     452             OPEN(ifile(3)+10,file='ps_INS_year'//chr2,form='formatted')                                               
     453             OPEN(97,file='prestot_year'//chr2,form='formatted')
     454             
     455
     456c  Sol or ls or both
     457c  Planetary mean surface pressure (Pa)
     458c  Equivalent pressure of CO2 ice at North Polar cap (Pa)
     459c  Equivalent pressure of CO2 ice at South Polar cap (Pa)
     460c  Total amount of CO2 on the planet (Pa)           
     461             
     462             IF (Time_unit == 1) THEN
     463             
     464              WRITE(ifile(1)+10,'(a)') '# Sol , Surface Pressure at VL1 at
     465     &        true (interpolated) altitude (Pa) , 
     466     &        Surface Pressure at VL1 at GCM altitude (Pa)'
     467
     468              WRITE(ifile(2)+10,'(a)') '# Sol , Surface Pressure at VL2 at
     469     &        true (interpolated) altitude (Pa) , 
     470     &        Surface Pressure at VL2 at GCM altitude (Pa)'
     471
     472              WRITE(ifile(3)+10,'(a)') '# Sol , Surface Pressure at Insight at
     473     &        true (interpolated) altitude (Pa) , 
     474     &        Surface Pressure at Insight at GCM altitude (Pa)'
     475
     476              WRITE(97,'(a)') '# Sol , Planetary Mean Surface Pressure (Pa) ,
     477     &        Equivalent pressure of CO2 ice at North Polar cap (Pa) ,
     478     &        Equivalent pressure of CO2 ice at South Polar cap (Pa) ,
     479     &        Total amount of CO2 on the planet (Pa)'     
     480                     
     481             ELSEIF (Time_unit == 2) THEN
     482
     483              WRITE(ifile(1)+10,'(a)') '# Ls (deg) , Surface Pressure at VL1 at
     484     &        true (interpolated) altitude (Pa) , 
     485     &        Surface Pressure at VL1 at GCM altitude (Pa)'
     486
     487              WRITE(ifile(2)+10,'(a)') '# Ls (deg) , Surface Pressure at VL2 at
     488     &        true (interpolated) altitude (Pa) , 
     489     &        Surface Pressure at VL2 at GCM altitude (Pa)'
     490
     491              WRITE(ifile(3)+10,'(a)') '# Ls (deg) , Surface Pressure at Insight at
     492     &        true (interpolated) altitude (Pa) , 
     493     &        Surface Pressure at Insight at GCM altitude (Pa)'
     494     
     495              WRITE(97,'(a)') '# Ls (deg) , Planetary Mean Surface Pressure (Pa) ,
     496     &        Equivalent pressure of CO2 ice at North Polar cap (Pa) ,
     497     &        Equivalent pressure of CO2 ice at South Polar cap (Pa) ,
     498     &        Total amount of CO2 on the planet (Pa)'     
     499             
     500             ELSE           
     501
     502              WRITE(ifile(1)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at VL1 at
     503     &        true (interpolated) altitude (Pa) , 
     504     &        Surface Pressure at VL1 at GCM altitude (Pa)'
     505
     506              WRITE(ifile(2)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at VL2 at
     507     &        true (interpolated) altitude (Pa) , 
     508     &        Surface Pressure at VL2 at GCM altitude (Pa)'
     509
     510              WRITE(ifile(3)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at Insight at
     511     &        true (interpolated) altitude (Pa) , 
     512     &        Surface Pressure at Insight at GCM altitude (Pa)'
     513
     514              WRITE(97,'(a)') '# Sol , Ls (deg) , Planetary Mean Surface Pressure (Pa) ,
     515     &        Equivalent pressure of CO2 ice at North Polar cap (Pa) ,
     516     &        Equivalent pressure of CO2 ice at South Polar cap (Pa) ,
     517     &        Total amount of CO2 on the planet (Pa)'
     518             
     519             ENDIF
     520             
    409521
    410522          ENDIF
    … …  
    442554
    443555
    444 c --------------Write output file xprestot-----------------------
    445 c  Sol ou ls ou les deux
    446 c  Ps_moy_planetaire (Pa)
    447 c  Pequivalente de glace de CO2 au Nord (si entierement sublimee) (Pa)
    448 c  Pequivalente de glace de CO2 au Sud (si entierement sublimee) (Pa)
     556c --------------Write output file prestot-----------------------
     557c  Sol or ls or both
     558c  Planetary mean surface pressure (Pa)
     559c  Equivalent pressure of CO2 ice at North Polar cap (Pa)
     560c  Equivalent pressure of CO2 ice at South Polar cap (Pa)
     561c  Total amount of CO2 on the planet (Pa)
    449562
    450563
    451564         IF(Time_unit == 1) THEN
    452               WRITE(97,'(4e16.6)') dayw,pstot*airtot1
    453      &       , captotN*g*airtot1, captotS*g*airtot1       
     565              WRITE(97,'(5e16.6)') dayw,pstot*airtot1
     566     &       , captotN*g*airtot1, captotS*g*airtot1,       
     567     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
    454568 
    455569         ELSEIF (Time_unit == 2) THEN   
    456               WRITE(97,'(4e16.6)') dayw_ls,pstot*airtot1
    457      &       , captotN*g*airtot1, captotS*g*airtot1
     570              WRITE(97,'(5e16.6)') dayw_ls,pstot*airtot1
     571     &       , captotN*g*airtot1, captotS*g*airtot1,
     572     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
    458573     
    459574         ELSE
    460              WRITE(97,'(5e16.6)') dayw,dayw_ls,pstot*airtot1
    461      &       , captotN*g*airtot1,captotS*g*airtot1
     575             WRITE(97,'(6e16.6)') dayw,dayw_ls,pstot*airtot1
     576     &       , captotN*g*airtot1,captotS*g*airtot1,
     577     &         pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1
    462578     
    463579                   
    … …  
    465581
    466582c----------------------------------------------------------------------
    467 c Loop on Viking sites:
     583c Loop on sites:
    468584c----------------------------------------------------------------------
    469585
    … …  
    474590         IF(.NOT.firstcal) THEN
    475591         
    476           DO iv=1,2
     592          DO iv=1,3
    477593
    478594             zp1=0.
    … …  
    483599             DO jj=0,1
    484600             
    485                 j=jvik(iv)+jj
     601                j=jsite(iv)+jj
    486602               
    487603                DO ii=0,1
    488604               
    489                    i=ivik(iv)+ii
     605                   i=isite(iv)+ii
    490606                   zt=zt+zw(ii,jj,iv)*t7(i,j)
    491607                   zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P)
    … …  
    506622             if (gh.eq.0) then ! if we don't have temperature values
    507623               ! assume a scale height of 10km
    508                zp2=zp1*exp(-(phivik(iv)-phisim(iv))/(3.73*1.e4))
     624               zp2=zp1*exp(-(phisite(iv)-phisim(iv))/(3.73*1.e4))
    509625             else
    510                zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh)
     626               zp2=zp1*exp(-(phisite(iv)-phisim(iv))/gh)
    511627             endif
    512628           
    513           WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh'
    514           WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh
     629          WRITE (*,*) 'iv,pstot,zp2, zp1, phisite(iv),phisim(iv),gh'
     630          WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phisite(iv),phisim(iv),gh
    515631          WRITE(*,*) "------"
    516632             
    517633
    518 c ------Write 2 files (1 for Vl1, 1 for VL2) xpsol ------
     634c ------Write 3 files (for Vl1, VL2, Insight) --------
    519635c  Sol or ls or both
    520 c  Ps site VLi (i=1,2) at GCM altitude (Pa)
    521 c  Ps site VLi (i=1,2) at true (interpolated) altitude (Pa)
     636c  Ps site VLi (i=1,2) or Inisght at true (interpolated) altitude (Pa) (zp2)
     637c  Ps site VLi (i=1,2) or Insight at GCM altitude (Pa) (zp1)
    522638             
    523639             IF(Time_unit == 1) THEN
    … …  
    527643             ELSE   
    528644                WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1
    529              ENDIF     
    530              
     645             ENDIF
     646             
     647         
    531648          ENDDO
    532649
    533650         ENDIF
     651                                 
    534652         
    535653         firstcal=.false.
    … …  
    556674
    557675
     676
    558677c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    559678c   End of loop on the diagfi files
    … …  
    562681      ENDDO
    563682
    564       PRINT*,'altitude of VL1',.001*phis(ivik(1),jvik(1))/g
    565       PRINT*,'altitude of VL2',.001*phis(ivik(2),jvik(2))/g
    566       DO iv=1,2
    567          PRINT*,'Viking',iv,'   i=',ivik(iv),'j  =',jvik(iv)
     683      PRINT*,'altitude of VL1',.001*phis(isite(1),jsite(1))/g
     684      PRINT*,'altitude of VL2',.001*phis(isite(2),jsite(2))/g
     685      PRINT*,'altitude of Ins',.001*phis(isite(3),jsite(3))/g
     686     
     687      DO iv=1,3
     688         PRINT*,'Site',iv,'   i=',isite(iv),'j  =',jsite(iv)
    568689         WRITE(6,7777)
    569      s   (rlonv(i)*180./pi,i=ivik(iv)-1,ivik(iv)+2)
     690     s   (rlonv(i)*180./pi,i=isite(iv)-1,isite(iv)+2)
    570691         print*
    571          DO j=jvik(iv)-1,jvik(iv)+2
     692         DO j=jsite(iv)-1,jsite(iv)+2
    572693            WRITE(6,'(f8.1,10x,5f7.1)')
    573      s   rlatu(j)*180./pi,(phis(i,j)/(g*1000.),i=ivik(iv)-1,ivik(iv)+2)
     694     s   rlatu(j)*180./pi,(phis(i,j)/(g*1000.),i=isite(iv)-1,
     695     &   isite(iv)+2)
    574696         ENDDO
    575697         print*
    … …  
    5837057777  FORMAT ('latitude/longitude',4f7.1)
    584706
    585 
    586 
     707c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     708c   Diurnal Average
     709c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     710
     711      write(*,*) 'ici'
     712      DO i=1,kyear+1
     713       WRITE(chr2(1:1),'(i1)') i
     714       IF(i.GE.9) WRITE(chr2,'(i2)') i
     715       DO iv=1,3
     716        if (iv==1) then
     717         
     718         open(ifile(iv), file = 'ps_VL1_year'//trim(trim(chr2)))
     719         open(ifile(iv)+20, file = 'ps_VL1_year'//trim(chr2)//'_diurnal')
     720         IF(Time_unit == 1) THEN         
     721          write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL1 at true
     722     & (interpolated) altitude (diurnal mean) , PS at VL1 at
     723     &  GCM altitude (diurnal mean)'
     724         ELSEIF (Time_unit == 2) THEN
     725          write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL1 at
     726     & true (interpolated) altitude (diurnal mean) , PS at VL1
     727     & at GCM altitude (diurnal mean)'   
     728         ELSE
     729          write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL1
     730     & at true (interpolated) altitude (diurnal mean) , PS at VL1
     731     & at GCM altitude (diurnal mean)'   
     732         ENDIF   
     733                 
     734        elseif (iv==2) then
     735         
     736         open(ifile(iv), file = 'ps_VL2_year'//trim(chr2))
     737         open(ifile(iv)+20, file = 'ps_VL2_year'//trim(chr2)//'_diurnal')
     738         IF(Time_unit == 1) THEN         
     739          write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL2 at true
     740     & (interpolated) altitude (diurnal mean) , PS at VL1 at
     741     & GCM altitude (diurnal mean)'
     742         ELSEIF (Time_unit == 2) THEN
     743          write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL2
     744     & at true (interpolated) altitude (diurnal mean) , PS at VL1
     745     & at GCM altitude (diurnal mean)'   
     746         ELSE
     747          write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL2
     748     & at true (interpolated) altitude (diurnal mean) , PS at VL1
     749     & at GCM altitude (diurnal mean)'   
     750         ENDIF 
     751                                 
     752        else
     753         open(ifile(iv), file = 'ps_INS_year'//trim(chr2))
     754         open(ifile(iv)+20, file = 'ps_INS_year'//trim(chr2)//'_diurnal')
     755         IF(Time_unit == 1) THEN         
     756          write(ifile(iv)+20,'(a)') '# Sol ,  PS at Insight at true
     757     & (interpolated) altitude (diurnal mean) , PS at VL1
     758     & at GCM altitude (diurnal mean)'
     759         ELSEIF (Time_unit == 2) THEN
     760          write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at Insight at true
     761     & (interpolated) altitude (diurnal mean) , PS at VL1
     762     & at GCM altitude (diurnal mean)'   
     763         ELSE
     764          write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at Insight
     765     & at true (interpolated) altitude (diurnal mean) , PS at VL1
     766     & at GCM altitude (diurnal mean)'   
     767         ENDIF                   
     768        endif   
     769       
     770        READ(ifile(iv),*)
     771        IF(Time_unit == 1) THEN
     772         READ(ifile(iv),*) dayw,zp2,zp1
     773        ELSEIF (Time_unit == 2) THEN   
     774         READ(ifile(iv),*) dayw_ls,zp2,zp1
     775         dayw = ls2sol(dayw_ls)
     776        ELSE   
     777         READ(ifile(iv),*) dayw,dayw_ls,zp2,zp1
     778        ENDIF
     779                             
     780        di=floor(dayw)
     781        di_prev = floor(dayw)
     782               
     783        DO k=1,nbmax
     784         ps_gcm_diurnal = 0.
     785         ps_diurnal = 0.
     786         compt_diurn = 0
     787                   
     788         DO WHILE (di==di_prev)
     789                   
     790          ps_gcm_diurnal = ps_gcm_diurnal + zp1
     791          ps_diurnal = ps_diurnal + zp2
     792          compt_diurn = compt_diurn + 1
     793          IF(Time_unit == 1) THEN
     794           READ(ifile(iv),*,end=777) dayw,zp2,zp1
     795          ELSEIF (Time_unit == 2) THEN   
     796           READ(ifile(iv),*,end=777) dayw_ls,zp2,zp1
     797           dayw=ls2sol(dayw_ls)
     798          ELSE   
     799           READ(ifile(iv),*,end=777) dayw,dayw_ls,zp2,zp1
     800          ENDIF
     801          di=floor(dayw)
     802                                               
     803         ENDDO
     804                       
     805         ps_gcm_diurnal = ps_gcm_diurnal/compt_diurn
     806         ps_diurnal = ps_diurnal/compt_diurn
     807         
     808         IF(Time_unit == 1) THEN         
     809                 
     810          write(ifile(iv)+20,'(i4,2e15.5)') di_prev, ps_diurnal,
     811     &         ps_gcm_diurnal
     812         
     813         ELSEIF (Time_unit == 2) THEN
     814         
     815          write(ifile(iv)+20,'(3e15.5)') dayw_ls, ps_diurnal,
     816     &         ps_gcm_diurnal   
     817         
     818         ELSE
     819         
     820          write(ifile(iv)+20,'(i4,3e15.5)') di_prev, dayw_ls,
     821     &         ps_diurnal, ps_gcm_diurnal       
     822         
     823         ENDIF
     824         
     825                 
     826         di_prev = di
     827                   
     828        ENDDO
     829        close(ifile(iv)+20)
     830        close(ifile(iv))
     831777    ENDDO
     832      ENDDO
     833
     834
     835
     836c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     837c   Harmonics
     838c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     839
     840       
     841      DO i=1,kyear+1
     842       WRITE(chr2(1:1),'(i1)') i
     843       IF(i.GE.9) WRITE(chr2,'(i2)') i     
     844       DO iv=1,3
     845       
     846        if (iv==1) then
     847         open(ifile(iv), file = 'ps_VL1_year'//trim(chr2)//'_diurnal')
     848         open(ifile(iv)+20, file = 'ps_VL1_year'//trim(chr2)//'_harmonics')
     849       
     850         IF(Time_unit == 1) THEN         
     851          write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL1 at true
     852     & (interpolated) altitude (harmonics fit) , PS at VL1 at
     853     &  GCM altitude (harmonics fit)'
     854         ELSEIF (Time_unit == 2) THEN
     855          write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL1 at
     856     & true (interpolated) altitude (harmonics fit) , PS at VL1
     857     & at GCM altitude (harmonics fit)'         
     858         ELSE
     859          write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL1
     860     & at true (interpolated) altitude (harmonics fit) , PS at VL1
     861     & at GCM altitude (harmonics fit)'         
     862         ENDIF                   
     863       
     864        elseif (iv==2) then
     865         open(ifile(iv), file = 'ps_VL2_year'//trim(chr2)//'_diurnal')
     866         open(ifile(iv)+20, file = 'ps_VL2_year'//trim(chr2)//'_harmonics')
     867         
     868         IF(Time_unit == 1) THEN         
     869          write(ifile(iv)+20,'(a)') '# Sol ,  PS at VL2 at true
     870     & (interpolated) altitude (harmonics fit) , PS at VL2 at
     871     &  GCM altitude (harmonics fit)'
     872         ELSEIF (Time_unit == 2) THEN
     873          write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at VL2 at
     874     & true (interpolated) altitude (harmonics fit) , PS at VL2
     875     & at GCM altitude (harmonics fit)'         
     876         ELSE
     877          write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at VL2
     878     & at true (interpolated) altitude (harmonics fit) , PS at VL2
     879     & at GCM altitude (harmonics fit)'         
     880         ENDIF                   
     881       
     882        else
     883         open(ifile(iv), file = 'ps_INS_year'//trim(chr2)//'_diurnal')
     884         open(ifile(iv)+20, file = 'ps_INS_year'//trim(chr2)//'_harmonics')
     885         
     886         IF(Time_unit == 1) THEN         
     887          write(ifile(iv)+20,'(a)') '# Sol ,  PS at Insight at true
     888     & (interpolated) altitude (harmonics fit) , PS at Insight at
     889     &  GCM altitude (harmonics fit)'
     890         ELSEIF (Time_unit == 2) THEN
     891          write(ifile(iv)+20,'(a)') '# Ls (deg) ,  PS at Insight at
     892     & true (interpolated) altitude (harmonics fit) , PS at Insight
     893     & at GCM altitude (harmonics fit)'         
     894         ELSE
     895          write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) ,  PS at Insight
     896     & at true (interpolated) altitude (harmonics fit) , PS at Insight
     897     & at GCM altitude (harmonics fit)'         
     898         ENDIF   
     899                         
     900        endif
     901
     902        READ(ifile(iv),*)
     903       
     904        DO k = 1,nbmax
     905         IF (Time_unit == 1) THEN 
     906          READ(ifile(iv),*,end=99) di_prev, ps_diurnal_tab(k),
     907     &         ps_gcm_diurnal_tab(k)
     908          call sol2ls(real(di_prev),ls_tab(k))
     909         ELSEIF (Time_unit == 2) THEN
     910          READ(ifile(iv),*,end=99) ls_tab(k), ps_diurnal_tab(k),
     911     &         ps_gcm_diurnal_tab(k)     
     912         ELSE 
     913          READ(ifile(iv),*,end=99) di_prev, ls_tab(k),
     914     &    ps_diurnal_tab(k), ps_gcm_diurnal_tab(k)
     915         ENDIF
     916                 
     917        ENDDO
     918       
     919       
     92099      ndata=k-1
     921       
     922        if(ls_tab(ndata).gt.350.) then     
     923 
     924        do k = 0,n_harmo
     925       
     926           call DiscreetFourierHn(ndata,ls_tab,ps_diurnal_tab,k,a,b)
     927           if(modulo(k,2)==0) then
     928             a_tab(k+1)= a
     929             b_tab(k+1)= b
     930           else
     931             a_tab(k+1)= -a
     932             b_tab(k+1)= -b
     933           endif                       
     934         
     935           call DiscreetFourierHn(ndata,ls_tab,ps_gcm_diurnal_tab,k,a,b)
     936           if(modulo(k,2)==0) then
     937             a_tab_gcm(k+1)= a
     938             b_tab_gcm(k+1)= b
     939           else
     940             a_tab_gcm(k+1)= -a
     941             b_tab_gcm(k+1)= -b
     942           endif
     943                   
     944        enddo
     945         
     946         write(ifile(iv)+20,'(a)') 'Fourrier coefficients ak/bk
     947     &(interpolated altitude)'
     948         write(ifile(iv)+20,'(a,7f)') '#',(a_tab(k),k=1,n_harmo+1)
     949         write(ifile(iv)+20,'(a,7f)') '#',(b_tab(k),k=1,n_harmo+1)
     950         write(ifile(iv)+20,'(a)') 'Fourrier coefficients ak/bk
     951     &(GCM altitude)'
     952         write(ifile(iv)+20,'(a,7f)') '#',(a_tab_gcm(k),k=1,n_harmo+1)
     953         write(ifile(iv)+20,'(a,7f)') '#',(b_tab_gcm(k),k=1,n_harmo+1)   
     954         
     955       
     956        do l = 1,669
     957       
     958         call sol2ls(real(l),ls_harm)
     959         ps_harm = a_tab(1)
     960         ps_gcm_harm = a_tab_gcm(1)
     961       
     962         do k = 1,n_harmo
     963       
     964            ps_harm = ps_harm + a_tab(k+1)*
     965     &cos((pi/180.)*k*(ls_harm))
     966     &+ b_tab(k+1)*sin((pi/180.)*k*(ls_harm))
     967 
     968            ps_gcm_harm = ps_gcm_harm + a_tab_gcm(k+1)*
     969     &cos((pi/180.)*k*(ls_harm))
     970     &+ b_tab_gcm(k+1)*sin((pi/180.)*k*(ls_harm))   
     971       
     972         enddo
     973
     974
     975         IF(Time_unit == 1) THEN         
     976                 
     977         write(ifile(iv)+20,'(i4,2e15.5)') l, ps_harm,
     978     &         ps_gcm_harm
     979         
     980         ELSEIF (Time_unit == 2) THEN
     981         
     982         write(ifile(iv)+20,'(3e15.5)') ls_harm, ps_harm,
     983     &         ps_gcm_harm       
     984         
     985         ELSE
     986         
     987         write(ifile(iv)+20,'(i4,3e15.5)') l,ls_harm, ps_harm,
     988     &         ps_gcm_harm       
     989         
     990         ENDIF
     991                       
     992        enddo
     993        close(ifile(iv))
     994        close(ifile(iv)+20)
     995       
     996        else 
     997         write(ifile(iv)+20,'(a)') 'not computed because
     998     & year is not complete'   
     999        endif ! (ls.gt.350)
     1000                   
     1001       ENDDO
     1002      ENDDO       
     1003     
     1004     
     1005     
    5871006      END
    5881007
    … …  
    6791098
    6801099      end subroutine sol2ls
     1100     
     1101     
     1102!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1103      real function ls2sol(ls)
     1104
     1105!  Returns solar longitude, Ls (in deg.), from day number (in sol),
     1106!  where sol=0=Ls=0 at the northern hemisphere spring equinox
     1107
     1108      implicit none
     1109
     1110!  Arguments:
     1111      real, intent(in) :: ls
     1112
     1113!  Local:
     1114      double precision xref,zx0,zteta,zz
     1115!        xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
     1116      double precision year_day
     1117      double precision peri_day,timeperi,e_elips
     1118      double precision pi,degrad
     1119      parameter (year_day=668.6d0) ! number of sols in a amartian year
     1120!      data peri_day /485.0/
     1121      parameter (peri_day=485.35d0) ! date (in sols) of perihelion
     1122!  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
     1123      parameter (timeperi=1.90258341759902d0)
     1124      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
     1125      parameter (pi=3.14159265358979d0)
     1126      parameter (degrad=57.2957795130823d0)
     1127
     1128      if (abs(ls).lt.1.0e-5) then
     1129         if (ls.ge.0.0) then
     1130            ls2sol = 0.0
     1131         else
     1132            ls2sol = real(year_day)
     1133         end if
     1134         return
     1135      end if
     1136
     1137      zteta = ls/degrad + timeperi
     1138      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
     1139      xref = zx0-e_elips*dsin(zx0)
     1140      zz = xref/(2.*pi)
     1141      ls2sol = real(zz*year_day + peri_day)
     1142      if (ls2sol.lt.0.0) ls2sol = ls2sol + real(year_day)
     1143      if (ls2sol.ge.year_day) ls2sol = ls2sol - real(year_day)
     1144
     1145      return
     1146      end
     1147
     1148!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1149
     1150
     1151!****************************************************************
     1152!*   Calculate the Fourier harmonic #n of a periodic discreet   *
     1153!*   function F(x) defined by ndata points.                     *
     1154!* ------------------------------------------------------------ *
     1155!* Inputs:                                                      *
     1156!*            ndata: number of points of discreet function.     *
     1157!*            X    : pointer to table storing xi abscissas.     *
     1158!*            Y    : pointer to table storing yi ordinates.     *
     1159!*                                                              *
     1160!* Outputs:                                                     *
     1161!*            a    : coefficient an of the Fourier series.      *
     1162!*            b:   : coefficient bn of the Fourier series.      *
     1163!* ------------------------------------------------------------ *
     1164!* Reference: "Mathematiques en Turbo-Pascal By Marc Ducamp and *
     1165!*             Alain Reverchon, 1. Analyse, Editions Eyrolles,  *
     1166!*             Paris, 1991" [BIBLI 03].                         *
     1167!*                                                              *
     1168!*                           F90 Version By J-P Moreau, Paris.  *
     1169!*                                  (www.jpmoreau.fr)           *
     1170!****************************************************************
     1171! Note: The Fourier series  of a periodic discreet function F(x)
     1172!       can be written under the form:
     1173!                   n=inf.
     1174!       F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x) 
     1175!                   n=1
     1176! ----------------------------------------------------------------
     1177         
     1178        Subroutine DiscreetFourierHn(ndata, X, Y, n, a, b)
     1179         real X(1:ndata), Y(1:ndata), a, b
     1180         integer ndata,n,i
     1181         real xi,T,om,wa,wb,wc,wd,wg,wh,wi,wl,wm,wn,wp
     1182         real PI
     1183         PI=4.d0*datan(1.d0)
     1184         T=X(ndata)-X(1); xi=(X(ndata)+X(1))/2.d0
     1185         om = 2*PI*n/T; a=0.d0; b=0.d0
     1186         do i=1, ndata-1
     1187             wa=X(i); wb=X(i+1)
     1188             wc=Y(i); wd=Y(i+1)
     1189          if (wa.ne.wb) then
     1190             wg = (wd-wc)/(wb-wa)
     1191             wh = om*(wa-xi); wi=om*(wb-xi)
     1192           if (n==0) then
     1193             a = a + (wb-wa)*(wc+wg/2.d0*(wb-wa))
     1194           else
     1195             wl = cos(wh); wm = sin(wh)
     1196             wn = cos(wi); wp = sin(wi)
     1197             a = a + wg/om*(wn-wl) + wd*wp - wc*wm
     1198             b = b + wg/om*(wp-wm) - wd*wn + wc*wl
     1199           end if
     1200         end if
     1201        end do
     1202        a = a/T; b = b/T
     1203        if (n.ne.0) then
     1204           a = a*2.d0/om; b = b*2.d0/om
     1205        end if
     1206        return
     1207       End     
Note: See TracChangeset for help on using the changeset viewer.