program check_co2cycle ! ******************************************************** ! Program to interpolate GCM surface pressure at Viking ! Landers 1/2 and Insight locations, compute diurnal ! average and make a harmonic reconstruction. ! Inputs: diagfi.nc, concat.nc, Xhistins.nc ! Outputs: prestot_yearX ! ps_SITE_yearX ! ps_SITE_yearX_diurnal ! ps_SITE_yearX_harmonics ! Author: J. Mauxion, adapted from xvik.F ! ******************************************************** implicit none !======================================================================= ! Pressure at Insight and Viking sites !======================================================================= !----------------------------------------------------------------------- ! declarations: !----------------------------------------------------------------------- include "netcdf.inc" real :: PI=4.0*atan(1.0) ! PI value integer :: itau ! Current timestep integer :: i,j,l,jj,ii ! Indexes for loops on lon, lat, ls and neighbour cells of sites real :: constR, g, unanj, day0, rad ! Constant reduced gas constant, gravtiational field, number of days in a year, ! starting day and Mars radius ! Declarations NETCDF: ! -------------------- character (len=100) :: varname ! Temporary varname integer :: ierr, nid, varid, dimid ! Error handler, file index and tmp file index, temporary var and dim id integer, dimension(3) :: start2D, count2D ! Start indices and number elements for ps and co2ice integer, dimension(4) :: start3D, count3D ! Same as above for temp ! Declarations Viking/Insight points: ! ----------------------------------- integer, dimension(3) :: isite(3),jsite(3),ifile(3) ! Sites indices (i/j=lon/lat), sites file units integer :: iv ! Loop index on sites character (len=10), dimension(3) :: namelist=(/'VL1','VL2','INS'/) !An array for name tags ! Nine reals below: lon, lat, phi(alt) declarations for the three data sites 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, dimension(3) :: lonsite, latsite, phisite, phisim ! lon/lat/phi(alt) arrays for the 3 sites location, corresponding phi in simulated grid ! Grid related variables: ! ----------------------- integer :: iend, jend ! # of grid points along lon/lat to be looped over integer :: lonlength ! # of grid points along longitude real, dimension(:), allocatable :: lon ! longitude integer :: latlength ! # of grid points along latitude real, dimension(:), allocatable :: lat ! latitude integer :: altlength ! # of grid point along altitude integer :: timelength ! # of grid point along time real, dimension(:), allocatable :: time ! time real :: ctrl(100) ! controle, whose size is know to be 100 real, dimension(:,:), allocatable :: phis ! To read phis (once) real, dimension(:,:), allocatable :: area ! To read area (once) ! Fields to read and compute: ! --------------------------- real, dimension(:,:,:), allocatable :: t ! To read temperature real, dimension(:,:), allocatable :: ps ! To read surface pressure real, dimension(:,:), allocatable :: co2ice ! To read co2 surface ice real, dimension(:,:), allocatable :: t7 ! To read temperature from 7th atmospheric layer (if relevant) real :: pstot, captotN, captotS ! To store calculated total surface pressure and North/South total co2 caps real :: zp1,zp2,zt ! To store interpolated GCM and True surface pressure and interpolated temperature real :: zw(0:1,0:1,3) ! Weight of each of the four near points real :: zalpha,zbeta ! Site relative lon/lat inside matching grid cell logical :: interpolated ! Is it an interpolated grid (using xios) ? real day,sol,dayw,dayw_ls ! To store current day inside diagfi, actual sol since beginning of the year ! and sol/ls to be written real airtot1,gh ! Inverse of total surface, RT product (from PV=nRT) integer iyear,kyear ! Current year, base year (equal as long as ongoing year is not over) character(len=2) :: chr2 ! Some character for int-to-char conversion ! Input/output related variables ! ------------------------------ character(len=80) :: file ! Complete input file path character(len=80) :: pathchmp, pathsor, nomfich ! Input, output dir, generic diagfi file name integer :: time_unit ! Time unit of output real :: ls2sol ! To store output of function ls2sol ! Diurnal average related variables: ! ---------------------------------- integer di,di_prev ! Current day, base day (equal as long as ongoing day is not over) integer k ! Index for loop and some counting real ps_gcm_diurnal, ps_diurnal ! To store GCM/true pressure for diurnal average integer compt_diurn ! Counter for number of timesteps per day for average computing ! Harmonics reconstruction related variables: ! ------------------------------------------- integer, parameter :: nbmax = 999999 ! Some size limit JM: should be cpt on the fly integer ndata ! Final data size real ls_harm ! To store temporary ls value real ls_tab(nbmax) ! To store tabulated ls value for harmonics reconstruction real ps_diurnal_tab(nbmax),ps_gcm_diurnal_tab(nbmax) ! Associated GCM/true pressure real a_tab(nbmax),b_tab(nbmax) ! Tabulated Fourrier coefficients, true pressure real a_tab_gcm(nbmax),b_tab_gcm(nbmax) ! Tabulated Fourrier coefficiens, GCM pressure real a,b ! Temporary Fourrier coefficients real ps_harm, ps_gcm_harm ! To store GCM/true pressure for harmonics fit computation integer, parameter :: n_harmo = 6 ! Number of harmonics of the fit !----------------------------------------------------------------------- ! initialisations: !----------------------------------------------------------------------- unanj=669. ! Mars year length in sols write(*,*) 'WARNING!!! Assuming ',unanj,' sols/year' !----------------------------------------------------------------------- ! Viking Lander and Insight coordinates: !-------------------------------------------------------------------- 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 write(*,*) 'Path to the diagfi files directory ?' read (*,*) pathchmp write(*,*) 'Path to the dir for outputs ?' read (*,*) pathsor write(*,*) 'Diagfi file name (without trailing .nc) ?' read(5,'(a)') nomfich write(*,*) 'Output file time axis in sol (1) in ls (2) ,or both (3) ?' read (*,*) time_unit !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Reading first diagfi to load the grid, topography and physi!-related variables !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% file=trim(pathchmp)//'/'//trim(nomfich)//'.nc' ! Current working file name write(*,*) "Complete file path is:", file write(*,*) "Loading first diagfi to read grid properties,",& " topography and physical constants..." ierr= NF_OPEN(file,NF_NOWRITE,nid) if (ierr.ne.NF_NOERR) then stop "ERROR: Can't open file." endif !!! ! Identifying whether it is interpolated from dynamico or native lonlat !!! ierr=NF_INQ_DIMID(nid,"time_counter",dimid) ! Get dimid if (ierr.eq.NF_NOERR) then write(*,*) "INFO: Found dimension 'time_counter'. ",& "I'll assume the grid was interpolated from a dynamico grid." interpolated=.true. else write(*,*) "INFO: Dimension 'time_counter' is missing. ",& "I'll assume the grid is a native lonlat grid." interpolated=.false. endif !!! ! Getting axis lengths and arrays !!! ! Longitude ierr=NF_INQ_DIMID(nid,"longitude",dimid) ! Get dimid if (ierr.ne.NF_NOERR) then ierr=NF_INQ_DIMID(nid,"lon",dimid) ! Try both spellings if (ierr.ne.NF_NOERR) then ! xios interpolation dumps lon_ ierr=NF_INQ_DIMID(nid,"lon_dom_out",dimid) if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get lon dimid." endif endif endif ierr=NF_INQ_DIMLEN(nid,dimid,lonlength) ! Get dim len if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to read lon length." endif ! If not interpolated, there is a redundant longitude in the diagfi if (.not.interpolated) then iend = lonlength-1 else iend = lonlength endif allocate(lon(lonlength)) ierr = NF_INQ_VARID(nid,"longitude",varid) ! Get varid if (ierr.ne.NF_NOERR) then ierr = NF_INQ_VARID(nid,"lon",varid) ! Try both spellings if (ierr.ne.NF_NOERR) then ! xios interpolation dumps lon_ ierr=NF_INQ_VARID(nid,"lon_dom_out",varid) if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get lon varid." endif endif endif ierr = NF_GET_VAR_REAL(nid,varid,lon) ! Read var if (ierr.ne.NF_NOERR) then stop " ERROR: Failed to read lon array." endif ! Latitude ierr=NF_INQ_DIMID(nid,"latitude",dimid) ! Get dimid if (ierr.ne.NF_NOERR) then ierr=NF_INQ_DIMID(nid,"lat",dimid) ! Try both spellings if (ierr.ne.NF_NOERR) then ! xios interpolation dumps lat_ ierr=NF_INQ_DIMID(nid,"lat_dom_out",dimid) if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get lat dimid." endif endif endif ierr=NF_INQ_DIMLEN(nid,dimid,latlength) ! Get dim len if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to read lat length." endif jend = latlength ! here there is a redundancy no matter what allocate(lat(latlength)) ierr = NF_INQ_VARID(nid,"latitude",varid) ! Get varid if (ierr.ne.NF_NOERR) then ierr = NF_INQ_VARID(nid,"lat",varid) ! Try both spellings if (ierr.ne.NF_NOERR) then ! xios interpolation dumps lat_ ierr = NF_INQ_VARID(nid,"lat_dom_out",varid) if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get lat varid." endif endif endif ierr = NF_GET_VAR_REAL(nid,varid,lat) ! Read var if (ierr.ne.NF_NOERR) then stop " ERROR: Failed to read lat array." endif ! Altitude ierr=NF_INQ_DIMID(nid,"altitude",dimid) ! Get dimid if (ierr.ne.NF_NOERR) then ierr=NF_INQ_DIMID(nid,"alt",dimid) ! Try both spellings if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get alt dimid." endif endif ierr=NF_INQ_DIMLEN(nid,dimid,altlength) ! Get dim len if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to read alt length." endif ! Phis allocate(phis(lonlength,latlength)) ierr = NF_INQ_VARID(nid,"phisinit",varid) ! Get varid if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get phisinit varid." endif ierr = NF_GET_VAR_REAL(nid,varid,phis) ! Read var if (ierr.ne.NF_NOERR) then stop " ERROR: Failed to read phisinit array." endif write(*,*) "Fields lon, lat, alt and phis passed." write(*,*) "Reading controle and loading relevant variables." ! Controle (we assumed length of 100) ierr = NF_INQ_VARID(nid,"controle",varid) if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get controle varid." endif ! Get the array ierr = NF_GET_VAR_REAL(nid,varid,ctrl) if (ierr.ne.NF_NOERR) then stop " ERROR: Failed to read controle array." endif write(*,*) "Field controle passed." write(*,*) "INFO:" !!! ! Controle is loaded, we can store relevant planet-related parameter ! and convert site altitude to surface geopeotential !!! day0 = ctrl(3) ! controle(3) = day_ini = starting day of simulation write (*,*) '>>> day0 = ' , day0 write (*,*) 'NOTE: Resetting day0 to zero in order to have each file ranging from sol=0 to 669.' write (*,*) ' This assumes you provided at least one year of integration with potential offset.' day0 = 0.0 ! Resetting to zero, otherwise one year+offset yields fractioned year over multiple files... rad = ctrl(5) ! controle(5) = rad = Mars radius write (*,*) '>>> rad = ' , rad g = ctrl(7) ! controle(7) = g = Mars gravitational field write (*,*) '>>> g = ' , g constR =1000.*8.314511/ctrl(8) ! controle(8)=mugaz = molar mass (g.mol-1) of atm. ! 1000 is conversion factor, 8.314... universal gas constant write (*,*) '>>> constR = ' , constR write(*,*) "Reading area." ! Area allocate(area(lonlength,latlength)) ierr = NF_INQ_VARID(nid,"area",varid) ! Get varid if (ierr.ne.NF_NOERR) then ierr = NF_INQ_VARID(nid,"aire",varid) ! Try both spellings if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get area varid." endif endif ierr = NF_GET_VAR_REAL(nid,varid,area) ! Read var if (ierr.ne.NF_NOERR) then stop " ERROR: Failed to read area array." endif ierr = NF_CLOSE(nid) ! No more initial reading write(*,*) "Field area passed." write(*,*) "Loading first diagfi to read grid properties,",& " topography and physical constants... Passed." !------------------------------------------------------------------- ! Convert altitudes to true surface geopotential now that we have g !------------------------------------------------------------------- write (*,*)'>>> Site altitudes', phisite do iv=1,3 phisite(iv)=phisite(iv)*g enddo !------------------------------------------------------------------------- ! Compute the inverse totale surface area now that we have the local area !------------------------------------------------------------------------- airtot1 = 0. do j=1,jend do i=1,iend ! Handle any lon except redundant one airtot1 = airtot1 + area(i,j) enddo enddo write (*,*) '>>> Total surface area', airtot1 airtot1 = 1./airtot1 !------------------------------------------------------ ! Weights for 4 near points at Viking and Insight !------------------------------------------------------ ! NOTE: We don't care about periodic completion ! (i.e. testing iend <= lonsite <= iend+1) because ! 1) Redundant point does not always exist in the file ! 2) the obs coordinates are hard-coded to values that ! are neither lon=pm 180 nor lat=pm 90 ! JM do iv=1,3 ! Initialize anyway isite(iv) = -1 jsite(iv) = -1 ! Locate index of GCM grid points near VL et al. do i=1,iend-1 ! We know longitudes are ordered -180...180 write(*,*) i, lonsite(iv),lon(i),lon(i+1) if ((lonsite(iv).ge.lon(i)).and.(lonsite(iv).le.lon(i+1))) then isite(iv)=i exit endif enddo do j=1,jend-1 ! We know that latitudes are ordered 90...-90 if ((latsite(iv).le.lat(j)).and.(latsite(iv).ge.lat(j+1))) then jsite(iv)=j exit endif enddo ! Check something was found if (isite(iv).eq.-1) stop "ERROR: Site longitude not found in grid." if (jsite(iv).eq.-1) stop "ERROR: Site latitude not found in grid." ! Relative site position inside identified grid cell zalpha=(lonsite(iv)-lon(isite(iv)))/(lon(isite(iv)+1)-lon(isite(iv))) zbeta=(latsite(iv)-lat(jsite(iv)))/(lat(jsite(iv)+1)-lat(jsite(iv))) ! Weight of each near point 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 !--------------------------------------------------------- ! True and model altitude at Viking and Insight locations !--------------------------------------------------------- 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 write(*,*) 'Real phisite at Viking/Insight locations: ', phisite write(*,*) 'Interpolated phisim at Viking/Insight locations: ', phisim write(*,*) 'Starting to loop on diagfi files...' !!! ! Init reading shapes and allocate reading arrays !!! count2D=(/lonlength,latlength,1/) count3D=(/lonlength,latlength,altlength,1/) allocate(ps(lonlength,latlength)) allocate(co2ice(lonlength,latlength)) allocate(t7(lonlength,latlength)) allocate(t(lonlength,latlength,altlength)) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Initializations for starting first year, opening files, write headers etc... !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Some random, different units for each site files ifile(1)=12 ifile(2)=13 ifile(3)=14 chr2="1" ! Its the char equivalent of year, starting at 1 iyear=1 ! See above kyear=iyear ! Base year is iyear, as we're starting first year OPEN(ifile(1)+10,file=trim(pathsor)//'/ps_VL1_year'//chr2,form='formatted') OPEN(ifile(2)+10,file=trim(pathsor)//'/ps_VL2_year'//chr2,form='formatted') OPEN(ifile(3)+10,file=trim(pathsor)//'/ps_INS_year'//chr2,form='formatted') OPEN(97,file=trim(pathsor)//'/prestot_year'//chr2,form='formatted') if (time_unit.eq.1) then write(ifile(1)+10,*) '# Sol , Surface Pressure at VL1 at true (interpolated)',& ' (interpolated) altitude (Pa) , Surface Pressure at VL1 at GCM altitude (Pa)' write(ifile(2)+10,*) '# Sol , Surface Pressure at VL2 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL2 at GCM altitude (Pa)' write(ifile(3)+10,*) '# Sol , Surface Pressure at Insight at true',& ' (interpolated) altitude (Pa) , Surface Pressure at Insight at GCM altitude (Pa)' write(97,*) '# 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)' else if (time_unit.eq.2) then write(ifile(1)+10,*) '# Ls (deg) , Surface Pressure at VL1 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL1 at GCM altitude (Pa)' write(ifile(2)+10,*) '# Ls (deg) , Surface Pressure at VL2 at true',& ' (interpolated) altitude (Pa) Surface Pressure at VL2 at GCM altitude (Pa)' write(ifile(3)+10,*) '# Ls (deg) , Surface Pressure at Insight at true',& ' (interpolated) altitude (Pa) , Surface Pressure at Insight at GCM altitude (Pa)' write(97,*) '# 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,*) '# Sol , Ls (deg) , Surface Pressure at VL1 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL1 at GCM altitude (Pa)' write(ifile(2)+10,*) '# Sol , Ls (deg) , Surface Pressure at VL2 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL2 at GCM altitude (Pa)' write(ifile(3)+10,*) '# Sol , Ls (deg) , Surface Pressure at Insight at true',& ' (interpolated) altitude (Pa) , Surface Pressure at Insight at GCM altitude (Pa)' write(97,*) '# 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 ! of if time_unit !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Loop on the diagfi files: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% do while(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50) write(*,*) '>>> nomfich : ',trim(nomfich) !---------------------------------------------------------------------- ! Open diagfi files : !---------------------------------------------------------------------- file=trim(pathchmp)//'/'//trim(nomfich)//'.nc' write(*,*) 'Opening file.nc: ', file ierr= NF_OPEN(file,NF_NOWRITE,nid) if (ierr.ne.NF_NOERR) then stop "ERROR: Can't open file." endif !---------------------------------------------------------------------- ! Read time : !---------------------------------------------------------------------- ierr = NF_INQ_DIMID (nid,"Time",dimid) ! Get dimid if (ierr.NE.NF_NOERR) then ierr = NF_INQ_DIMID (nid,"time_counter",dimid) ! Try both spellings if (ierr.NE.NF_NOERR) then stop "ERROR: Failed to get time dimid." endif endif ierr=NF_INQ_DIMLEN(nid,dimid,timelength) ! Get dim len if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to read time length." endif allocate(time(timelength)) ierr = NF_INQ_VARID(nid,"Time",varid) ! Get varid if (ierr.ne.NF_NOERR) then ierr = NF_INQ_VARID(nid,"time_instant",varid) ! Try both spellings if (ierr.ne.NF_NOERR) then stop "ERROR: Failed to get time varid." endif endif ierr = NF_GET_VAR_REAL(nid,varid,time) ! Read var if (ierr.ne.NF_NOERR) then stop " ERROR: Failed to read time array." endif write(*,*) 'time(1:10) ', time(1:10) !====================================================================== ! Begin the loop on states in inputs files : !====================================================================== do itau=1,timelength !---------------------------------------------------------------------- ! Open/close files if a new year was set at the end of last loop !--------------------------------------------------------------------- if (iyear.ne.kyear) then ! kyear keeps base year in memory, while iyear=kyear ! the ongoing year is not finished yet, otherwise iyear ! was updated and we need new files write(chr2(1:1),'(i1)') iyear ! Convert new year int to character if(iyear.gt.9) then ! If next year needs two digits, adjust chr2 write(chr2,'(i2)') iyear endif write (*,*) '>>> Finished year ',kyear ! This is finishing year write (*,*) '>>> Starting year ',iyear ! This is starting year write (*,*) '>>> chr2:',trim(chr2) ! Check int-to-char conversion kyear=iyear ! New base year ! Close previous year files, open next year files do ii=1,3 CLOSE(10+ifile(ii)) enddo CLOSE(97) OPEN(ifile(1)+10,file=trim(pathsor)//'/ps_VL1_year'//chr2,form='formatted') OPEN(ifile(2)+10,file=trim(pathsor)//'/ps_VL2_year'//chr2,form='formatted') OPEN(ifile(3)+10,file=trim(pathsor)//'/ps_INS_year'//chr2,form='formatted') OPEN(97,file=trim(pathsor)//'/prestot_year'//chr2,form='formatted') ! Write headers: ! Sol or ls or both ! 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) if (time_unit.eq.1) then write(ifile(1)+10,*) '# Sol , Surface Pressure at VL1 at true (interpolated)',& ' (interpolated) altitude (Pa) , Surface Pressure at VL1 at GCM altitude (Pa)' write(ifile(2)+10,*) '# Sol , Surface Pressure at VL2 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL2 at GCM altitude (Pa)' write(ifile(3)+10,*) '# Sol , Surface Pressure at Insight at true',& ' (interpolated) altitude (Pa) , Surface Pressure at Insight at GCM altitude (Pa)' write(97,*) '# 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)' else if (time_unit.eq.2) then write(ifile(1)+10,*) '# Ls (deg) , Surface Pressure at VL1 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL1 at GCM altitude (Pa)' write(ifile(2)+10,*) '# Ls (deg) , Surface Pressure at VL2 at true',& ' (interpolated) altitude (Pa) Surface Pressure at VL2 at GCM altitude (Pa)' write(ifile(3)+10,*) '# Ls (deg) , Surface Pressure at Insight at true',& ' (interpolated) altitude (Pa) , Surface Pressure at Insight at GCM altitude (Pa)' write(97,*) '# 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,*) '# Sol , Ls (deg) , Surface Pressure at VL1 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL1 at GCM altitude (Pa)' write(ifile(2)+10,*) '# Sol , Ls (deg) , Surface Pressure at VL2 at true',& ' (interpolated) altitude (Pa) , Surface Pressure at VL2 at GCM altitude (Pa)' write(ifile(3)+10,*) '# Sol , Ls (deg) , Surface Pressure at Insight at true',& ' (interpolated) altitude (Pa) , Surface Pressure at Insight at GCM altitude (Pa)' write(97,*) '# 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 ! of if time_unit endif ! of if iyear!=kyear ! Adapt reading starting point start2D=(/1,1,itau/) start3D=(/1,1,1,itau/) !---------------------------------------------------------------------- ! Read fields: !---------------------------------------------------------------------- ! Load ps ierr = NF_INQ_VARID (nid, "ps", varid) ierr = NF_GET_VARA_REAL(nid,varid,start2D,count2D,ps) if (ierr.NE.NF_NOERR) then stop "ERROR: Failed to read ." endif write(*,*) 'ps ', ps(lonlength/2,latlength/2) ! Load Temperature ierr = NF_INQ_VARID (nid, "temp", varid) ierr = NF_GET_VARA_REAL(nid,varid,start3D,count3D,t) if (ierr.NE.NF_NOERR) then write(*,*) 'INFO: Failed to read .' ! Ehouarn: proceed anyways write(*,*)'--> Setting temperature to zero.' t(:,:,:)=0.0 write(*,*)'--> Looking for temp7 (temp in 7th layer).' ierr=NF_INQ_VARID(nid,"temp7", varid) if (ierr.eq.NF_NOERR) then write(*,*) " OK, found temp7 variable." ierr=NF_GET_VARA_REAL(nid,varid,start2D,count2D,t7) if (ierr.ne.NF_NOERR) then stop 'ERROR: Failed to read temp7.' endif else ! no 'temp7' variable write(*,*)' No temp7 variable either.' write(*,*)' I will have to do without, setting t7 to zero.' t7(:,:)=0.0 endif else ! t() was successfully loaded, copy 7th layer to t7() t7=t(:,:,7) endif write(*,*) 't7 ', t7(lonlength/2,latlength/2) ! Load co2ice ierr = NF_INQ_VARID (nid, "co2ice", varid) ierr = NF_GET_VARA_REAL(nid, varid,start2D,count2D, co2ice) if (ierr.NE.NF_NOERR) then stop 'ERROR: Failed to read .' endif write(*,*) 'co2ice ', co2ice(lonlength/2,latlength/2) !---------------------------------------------------------------------- ! Handle calendar !--------------------------------------------------------------------- day=time(itau) ! Total #sols since beginning of current diagfi print*,'>>> Current diagfi file: day = ',day sol=day+day0 do while (sol.gt.unanj) sol=sol-unanj ! With offset day0 it can be more than 669, remove full years enddo write(*,*) '>>> Current year: sol = ',sol,', iyear = ',iyear dayw = sol ! Sol to write in file call sol2ls(sol,dayw_ls) ! dayw_ls=Ls to write in file !---------------------------------------------------------------------- ! Compute average planetary pressure !--------------------------------------------------------------------- ! Reset monitoring variables pstot=0. captotS=0. captotN=0. ! Total surface pressure do j=1,jend do i=1,iend pstot=pstot+area(i,j)*ps(i,j) enddo enddo ! Total northern CO2 ice cap do j=1,jend/2+1 do i=1,iend captotN = captotN+area(i,j)*co2ice(i,j) enddo enddo ! Total southern CO2 ice cap do j=jend/2+1,jend do i=1,iend captotS = captotS+area(i,j)*co2ice(i,j) enddo enddo !--------------Write output file prestot----------------------- ! Sol or ls or both ! 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) if (time_unit.eq.1) then write(97,'(5e16.6)') dayw,pstot*airtot1,captotN*g*airtot1,& captotS*g*airtot1,pstot*airtot1+captotN*g*airtot1+captotS*g*airtot1 else if (time_unit.eq.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 !---------------------------------------------------------------------- ! Loop on sites: !---------------------------------------------------------------------- do iv=1,3 ! Resetting for current calculation zp1=0. zp2=0. zt=0. !---------------------------------------------------------------------- ! Interpolate using temperature in the 7th layer, of surface pressure !---------------------------------------------------------------------- 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 ! Basically gravity times scale height !---------------------------------------------------------------------- ! Surface pressure extrapolated using temp. from 7th atmospheri!layer !---------------------------------------------------------------------- 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))/(g*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(*,*) "------" !------Write 3 files (for VL1, VL2, Insight) -------- ! Sol or ls or both ! Ps site VLi (i=1,2) or Inisght at true (interpolated) altitude (Pa) (zp2) ! Ps site VLi (i=1,2) or Insight at GCM altitude (Pa) (zp1) if(time_unit.eq.1) then write(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1 else if (time_unit.eq.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 ! of do iv !====================================================================== ! End of reading/interpolating variables in the diagfi file !====================================================================== if (sol.ge.(unanj-1.e-5)) then ! End of year reached (with some roundoff margin) ! Increment iyear iyear=iyear+1 endif enddo ! of do itau ierr= NF_CLOSE(nid) day0 = day0 + day ! Update day0 at the end of each diagfi file write(*,*) "End of file ", nomfich write(*,*) "New day0 is ", day0 write(*,*) "Enter new file name (without trailing .nc) or return to end ?" read(5,'(a)') nomfich !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! End of loop on the diagfi files !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% deallocate(time) ! deallocate time for next diagfi (in case different len) enddo ! It's over, let's close the files do ii=1,3 CLOSE(10+ifile(ii)) enddo CLOSE(97) ! Some diagnostics output write(*,*) 'Altitude of VL1 ',.001*phis(isite(1),jsite(1))/g write(*,*) 'Altitude of VL2 ',.001*phis(isite(2),jsite(2))/g write(*,*) 'Altitude of Ins ',.001*phis(isite(3),jsite(3))/g do iv=1,3 write(*,*) '===> Site ',iv,' i=',isite(iv),'j =',jsite(iv) write(*,*) 'Four nearest points: ' do jj=0,1 j=jsite(iv)+jj do ii=0,1 i=isite(iv)+ii write(*,*) '>>> lon/lat/alt: ', lon(i),lat(j),phis(i,j)/(g*1000.) write(*,*) '>>> interpolation weight: ', zw(ii,jj,iv) enddo enddo write(*,*) '===> Site interpolated altitude (km): ', phisim(iv)/g/1000. enddo write(*,*) "Total CO2 monitoring and surface pressure interpolation: DONE." !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Diurnal Average !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% write(*,*) 'Computing diurnal average...' if (iyear.ne.kyear) then ! If they are different then we actually covered iyear-1 iyear = iyear-1 endif do i=1,iyear ! At this stage, iyear is the total number of year previously processed (including incomplete). write(chr2(1:1),'(i1)') i if(i.gt.9) then ! i2 if 2 digits. write(chr2,'(i2)') i endif do iv=1,3 ! Browse the three sites if (iv==1) then ! Open corresponding site/year for reading open(ifile(iv), file = trim(pathsor)//'/ps_VL1_year'//trim(chr2)) ! Open a diurnal file once average is performed open(ifile(iv)+20, file = trim(pathsor)//'/ps_VL1_year'//trim(chr2)//'_diurnal') if(time_unit.eq.1) then ! Again, different treatment depending on time unit write(ifile(iv)+20,*) '# Sol , PS at VL1 at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' else if (time_unit.eq.2) then write(ifile(iv)+20,*) '# Ls (deg) , PS at VL1 at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' else write(ifile(iv)+20,*) '# Sol , Ls (deg) , PS at VL1 at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' endif else if (iv==2) then open(ifile(iv), file = trim(pathsor)//'/ps_VL2_year'//trim(chr2)) open(ifile(iv)+20, file = trim(pathsor)//'/ps_VL2_year'//trim(chr2)//'_diurnal') if(time_unit.eq.1) then write(ifile(iv)+20,*) '# Sol , PS at VL2 at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' else if (time_unit.eq.2) then write(ifile(iv)+20,*) '# Ls (deg) , PS at VL2 at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' else write(ifile(iv)+20,*) '# 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 = trim(pathsor)//'/ps_INS_year'//trim(chr2)) open(ifile(iv)+20, file = trim(pathsor)//'/ps_INS_year'//trim(chr2)//'_diurnal') if(time_unit.eq.1) then write(ifile(iv)+20,*) '# Sol , PS at Insight at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' else if (time_unit.eq.2) then write(ifile(iv)+20,*) '# Ls (deg) , PS at Insight at true (interpolated)',& ' altitude (diurnal mean) , PS at VL1 at GCM altitude (diurnal mean)' else write(ifile(iv)+20,*) '# 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),*) ! Skip first line, then read first data line if(time_unit.eq.1) then read(ifile(iv),*,iostat=ierr) dayw,zp2,zp1 else if (time_unit.eq.2) then read(ifile(iv),*,iostat=ierr) dayw_ls,zp2,zp1 dayw = ls2sol(dayw_ls) else read(ifile(iv),*,iostat=ierr) dayw,dayw_ls,zp2,zp1 endif ! Make sure first line reading went well (and also initiate for next do while ierr loop) if (ierr.ne.0) then if (ierr.le.0) then write(*,*) "ERROR: Hitting EOF before entering time loop.",& "Check your ps_"//trim(namelist(iv))//"_year_"//trim(chr2)//" file." stop else write(*,*) "ERROR: Something went wrong when reading first line of data.",& "Check your ps_"//trim(namelist(iv))//"_year_"//trim(chr2)//" file." stop endif endif di=floor(dayw) ! First day (0<=dayw<1) di_prev = floor(dayw) ! Initializing previous day as first day for starting do while (ierr.eq.0) ! Loop as long as we do not reach EOF or an error at reading ! Reset for current day ps_gcm_diurnal = 0. ps_diurnal = 0. compt_diurn = 0 do while (di.eq.di_prev) ! Sum over each time step within a day ps_gcm_diurnal = ps_gcm_diurnal + zp1 ps_diurnal = ps_diurnal + zp2 compt_diurn = compt_diurn + 1 if(time_unit.eq.1) then read(ifile(iv),*,iostat=ierr) dayw,zp2,zp1 else if (time_unit.eq.2) then read(ifile(iv),*,iostat=ierr) dayw_ls,zp2,zp1 dayw=ls2sol(dayw_ls) else read(ifile(iv),*,iostat=ierr) dayw,dayw_ls,zp2,zp1 endif ! Check that reading went well if (ierr.gt.0) then write(*,*) "ERROR: Something went wrong when reading data.",& "Check your ps_"//trim(namelist(iv))//"_year_"//trim(chr2)//" file." stop else if (ierr.lt.0) then ! We reached EOF, if the day isn't over we need to break the loop (otherwise infinite ?) exit endif di=floor(dayw) enddo ! of while di==di_prev ! Compute average and write down for current day ps_gcm_diurnal = ps_gcm_diurnal/compt_diurn ps_diurnal = ps_diurnal/compt_diurn if(time_unit.eq.1) then write(ifile(iv)+20,'(i4,2e15.5)') di_prev, ps_diurnal,ps_gcm_diurnal else if (time_unit.eq.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 ! of do while ierr.eq.0 ! Year is done, close files close(ifile(iv)+20) close(ifile(iv)) enddo ! pf do iv=1,3 enddo ! of do i=1,kyear+1 write(*,*) 'Computing diurnal average: DONE.' !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Harmonics !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% write(*,*) 'Computing harmonics...' do i=1,iyear ! Same, iyear should have good value write(chr2(1:1),'(i1)') i if(i.gt.9) then write(chr2,'(i2)') i endif k=0 ! Reset year-lasting counter do iv=1,3 if (iv==1) then open(ifile(iv),file=trim(pathsor)//'/ps_VL1_year'//trim(chr2)//'_diurnal') open(ifile(iv)+20,file=trim(pathsor)//'/ps_VL1_year'//trim(chr2)//'_harmonics') if(time_unit.eq.1) then write(ifile(iv)+20,*) '# Sol , PS at VL1 at true (interpolated) altitude (harmonics fit) ,',& ' PS at VL1 at GCM altitude (harmonics fit)' else if (time_unit.eq.2) then write(ifile(iv)+20,*) '# Ls (deg) , PS at VL1 at true (interpolated) altitude (harmonics fit) ,',& ' PS at VL1 at GCM altitude (harmonics fit)' else write(ifile(iv)+20,*) '# Sol , Ls (deg) , PS at VL1 at true (interpolated) altitude (harmonics fit) ,',& ' PS at VL1 at GCM altitude (harmonics fit)' endif else if (iv==2) then open(ifile(iv),file=trim(pathsor)//'/ps_VL2_year'//trim(chr2)//'_diurnal') open(ifile(iv)+20,file=trim(pathsor)//'/ps_VL2_year'//trim(chr2)//'_harmonics') if(time_unit.eq.1) then write(ifile(iv)+20,*) '# Sol , PS at VL2 at true (interpolated) altitude (harmonics fit) ,',& ' PS at VL2 at GCM altitude (harmonics fit)' else if (time_unit.eq.2) then write(ifile(iv)+20,*) '# Ls (deg) , PS at VL2 at true (interpolated) altitude (harmonics fit) ,',& ' PS at VL2 at GCM altitude (harmonics fit)' else write(ifile(iv)+20,*) '# 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=trim(pathsor)//'/ps_INS_year'//trim(chr2)//'_diurnal') open(ifile(iv)+20,file=trim(pathsor)//'/ps_INS_year'//trim(chr2)//'_harmonics') if(time_unit.eq.1) then write(ifile(iv)+20,*) '# Sol , PS at Insight at true (interpolated) altitude (harmonics fit) ,',& ' PS at Insight at GCM altitude (harmonics fit)' else if (time_unit.eq.2) then write(ifile(iv)+20,*) '# Ls (deg) , PS at Insight at true (interpolated) altitude (harmonics fit) ,',& ' PS at Insight at GCM altitude (harmonics fit)' else write(ifile(iv)+20,*) '# 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),*,iostat=ierr) ! Skip header and initiate ierr ! Make sure first line reading went well if (ierr.ne.0) then if (ierr.le.0) then write(*,*) "ERROR: Hitting EOF before entering time loop.",& "Check your ps_"//trim(namelist(iv))//"_year_"//trim(chr2)//"_diurnal file." stop else write(*,*) "ERROR: Something went wrong when reading first line of data.",& "Check your ps_"//trim(namelist(iv))//"_year_"//trim(chr2)//"_diurnal file." stop endif endif do while (ierr.eq.0) k=k+1 ! Update counter (skips header line) if (time_unit.eq.1) then read(ifile(iv),*,iostat=ierr) di_prev,ps_diurnal_tab(k),ps_gcm_diurnal_tab(k) call sol2ls(real(di_prev),ls_tab(k)) else if (time_unit.eq.2) then read(ifile(iv),*,iostat=ierr) ls_tab(k),ps_diurnal_tab(k),ps_gcm_diurnal_tab(k) else read(ifile(iv),*,iostat=ierr) di_prev,ls_tab(k),ps_diurnal_tab(k), ps_gcm_diurnal_tab(k) endif ! Check reading went well if (ierr.gt.0) then write(*,*) "ERROR: Something went wrong when reading data.",& "Check your ps_"//trim(namelist(iv))//"_year_"//trim(chr2)//"_diurnal file." stop endif enddo ndata=k-1 ! Cause updated counter one too far if(ls_tab(ndata).ge.(360.-1.e-3)) then ! We need a full year, roundoff margin do k = 0,n_harmo ! Loop over the chosen nb of harmonics ! Compute Fourier coefficient for sites interpolated 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 ! Compute Fourier coefficient for GCM native 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 Fourrier coefficients write(ifile(iv)+20,*) '# Fourrier coefficients ak/bk (interpolated altitude)' write(ifile(iv)+20,'(a,7E12.4)') '# ',(a_tab(k),k=1,n_harmo+1) write(ifile(iv)+20,'(a,7E12.4)') '# ',(b_tab(k),k=1,n_harmo+1) write(ifile(iv)+20,*) '# Fourrier coefficients ak/bk (GCM altitude)' write(ifile(iv)+20,'(a,7E12.4)')'# ',(a_tab_gcm(k),k=1,n_harmo+1) write(ifile(iv)+20,'(a,7E12.4)')'# ',(b_tab_gcm(k),k=1,n_harmo+1) ! Reconstruct from Fourier coefficients, discrete temporal axis of 1 sol ! (but converted in Ls) 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 ! Write down the reconstructed field at day l if(time_unit.eq.1) then write(ifile(iv)+20,'(i4,2e15.5)') l,ps_harm,ps_gcm_harm else if (time_unit.eq.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 ! of do l=1,669 else write(*,*) 'INFO: not computed because year ',i,' is not complete.' endif ! of ls.ge.360 close(ifile(iv)) close(ifile(iv)+20) enddo ! of do iv enddo ! of do kyear write(*,*) 'Computing harmonics: DONE.' 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 !============================================================================== 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: eccentri!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) endif return endif 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 function ls2sol !**************************************************************** !* Calculate the Fourier harmoni!#n of a periodi!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 Mar!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 periodi!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 subroutine DiscreetFourierHn