

      subroutine pres0(dset,lat,lon,solar,loctim,pres,alt,ierr)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Purpose:
c =======
c  Pres0 si a subroutine designed to determine accurate 
c  surface pressure (as well as MOLA altitude) on Mars at
c  any given point in latitude, longitude, local time and solar longitude.
c  It uses:
c     1) Reference pressure measurements at Viking Lander 1 site (file VL1.ls)
c     2) High resolution 32 pix/deg MOLA topography (file mola32.nc)
c     3) data issued from the Mars Climate Database (file ps_clim.nc) but is
c  independent from it, so you don't need the full MCD to be installed to
c  use pres0.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Changes and updates:
c ===================
c  -Initial version: F. Forget, Y. Wanherdrick, LMD 2002
c  -Interpolation improvement and code acceleration F.Forget 10/2005 02/2006
c  -Removed spline interpolation; routine improvements E. Millour 11/2006
c  -Improved interpolation (switched to bilinear interpolation of log(pressure))
c  and added 'dset' (path to datafiles) argument E. Millour 12/2006
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Pres0 Return codes: (argument ierr)
c ==================
c ierr=0: No problem
c ierr=1: Bad input argument
c ierr=2: (non fatal) error in routine mola.F
c ierr=3: Failed to read VL1.ls file
c ierr=4: Problem reading data from file ps_clim.nc
c ierr=5: Problem with bilinear interpolation in ps_mcd.F
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

c arguments:
c    inputs:
      character*(*) dset ! Path to MCD datafiles (must end with a '/')
      real lat         ! Latitude coordinate of the point (North degrees)
      real lon         ! Longitude coordinate of the point (East degrees)
      real solar       ! Solar longitude (degrees)
      real loctim      ! Local time at longitude lon (hours)

c   outputs:
      real pres        ! high resolution surface pressure (Pa)
      real alt         ! surface altitude from MOLA (m)
      integer ierr     ! Control variable

c local variables:
      real factcor  ! Pressure correction factor
      real zlon     ! East longitude (deg.) [0:360]
      real H        ! Scale height (m)
      real oro      ! MCD orography at location lon,lat
      real temp     ! Temperature at altitude ~1km
      
! Check input arguments
      if ((lat.gt.90).or.(lat.lt.-90)) then
         write(6,*)"pres0 error: Bad latitude value. lat= ",lat
         ierr=1
         return
      endif

      if ((lon.gt.360).or.(lon.lt.-360)) then
         write(6,*)"pres0 error: Bad longitude value. lon= ",lon
         ierr=1
         return
      endif
      zlon=lon
      if(zlon.lt.0.) zlon =zlon + 360.

      if ((solar.gt.360).or.(solar.lt.0)) then
         write(6,*)"pres0 error: Bad solar longitude value. ls=",
     &        solar
         ierr=1
         return
      endif

      if ((loctim.gt.24).or.(loctim.lt.0)) then
        write(6,*)"pres0 error: Bad local time value. loctim= ",
     &              loctim
         ierr=1
         return
      endif

      ierr=0

! 1. Read MOLA orography
      call mola(dset,lat,zlon,alt,ierr)
      if (ierr.ne.0) return

! 2. Compute a correction factor using VL1 pressure
      call calc_factcor(dset,solar,factcor,ierr)
      if (ierr.ne.0) return


! 3. Get MCD temperature (at ~1km above surface) and surface pressure
      call ps_mcd(dset,lat,zlon,solar,loctim,pres,temp,oro,ierr)
      if (ierr.ne.0) return

! 4. Build Scale height H=R.T/g
! Gas Constant: R (m2.s-2.K-1) = 191
! surface gravity: g (m.s-2) = 3.73
! scale height: H=R.T/g
      H=191*temp/3.73

! 5. Compute topography-corrected surface pressure
      pres=pres*factcor*exp(-(alt-oro)/H)

      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine read_vl(dset,tab,ierr)
!Routine to read the VL1 (or VL2) mean viking pressure files.
      implicit none

c Arguments:
c Inputs:
      character*(*) dset ! Path to MCD datafiles 
      real tab(669,2) ! tab(:,1) contains Ls and tab(:,2) contains pressure
c Output:
      integer ierr ! returned status code (==0 if OK)

c local variables
      integer vik_num,count
      character*10 vlfile ! (smoothed) surface pressure measured by VL
      data vlfile/'VL1.ls'/
      character*140 datafile ! full path+file name

C 1. Open file VL1.ls
      vik_num=18
      datafile=dset//vlfile
      open(vik_num,file=datafile,status="old",iostat=ierr,
     &     form="formatted")

      if (ierr.ne.0) then
         write(6,*)"read_vl error: Could not open file ",vlfile
         ierr=3
         return
      endif

C 2. Read the 669 lines of datafile

      do count=1,669
        read(vik_num,*) tab(count,1),tab(count,2)
      enddo

C 3. Close file and return
      close(vik_num)

      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine calc_factcor(dset,solar,factcor,ierr)
c Compute a correction factor between Viking1 surface pressure data
c and MCD output at the same point
c Note: Position of VL1 (and VL2 and Pathfinder) are taken from
c       Konopliv et al. 2006 and are expressed in the IAU2000 reference system.
c       Corresponding altitudes were computed using bilinear interpolation of
c       MOLA 32 topography data (and using areoid computed from the mgm1025
c       spherical harmonics coefficients) 
      implicit none

! Arguments:
!     inputs:
      character*(*) dset ! Path to MCD datafiles 
      real solar	! Solar longitude (degrees)
!     outputs:
      real factcor	! correction factor
      integer ierr	! status ierr=0 if everything OK

! Local variables
      real vik1_lat, vik1_lon, vik1_alt
      real vik2_lat, vik2_lon, vik2_alt
      real path_lat, path_lon, path_alt
!  Viking1 coordinates
!      parameter(vik1_lat=22.2692, vik1_lon=311.7783, vik1_alt=-3627) !old ref.
      parameter(vik1_lat=22.269628,vik1_lon=312.050381,
     &          vik1_alt=-3637.1396)
!  Viking2 coordinates
!      parameter(vik2_lat=47.6680, vik2_lon=134.0100, vik2_alt=-4505) ! old ref.
      parameter(vik2_lat=47.668093,vik2_lon=134.282101,
     &          vik2_alt=-4497.0474)
!  Pathfinder coordinates
!      parameter(path_lat=19.0949, path_lon=326.4762, path_alt=-3682) ! old ref.
      parameter(path_lat=19.099641, path_lon=326.747364,
     &          path_alt=-3681.9993)

      real presmean,lsinf,lssup,psinf,pssup,presvik !,alt
      real prescalc(12) ! surface pressure at VL1 site, at times 2,4,...24 hours
      real tempcalc(12) ! temperature above (~1km) VL1, at times 2,4,...24 hours
      real orocalc ! MCD orography, at VL1 landing site
      real H ! scale height
      integer iloop
      logical firstcall
      data firstcall/.true./
      save firstcall
      real vik_tab(669,2) ! tab(:,1) contains Ls and tab(:,2) contains pressure
      save vik_tab        ! records by VL1

C 1. If first call, load (smoothed) VL1 measured pressure

      if (firstcall) then
         firstcall=.false.
         call read_vl(dset,vik_tab,ierr)
         if (ierr.ne.0) return
      endif

C 2. Call ps_mcd to get the pressure at viking1 landing site location
c    at given solar longitude and local time from the MCD.

      do iloop=1,12
         call ps_mcd(dset,vik1_lat,vik1_lon,solar,2.0*iloop,
     &        prescalc(iloop),tempcalc(iloop),orocalc,ierr)
         if (ierr.ne.0) return
      enddo

C 3. Compute the correction factor due to difference between orocalc
C    and viking1 altitude

      do iloop=1,12
! Gas Constant: R (m2.s-2.K-1) = 191
! surface gravity: g (m.s-2) = 3.73
         H=191*tempcalc(iloop)/3.73
         prescalc(iloop)=prescalc(iloop)*exp(-(vik1_alt-orocalc)/H)
      enddo

C 4. Compute average (over the day) pressure
      presmean=0.
      do iloop=1,12
         presmean = presmean + prescalc(iloop)/12.
      enddo

C 5. Get encompassing VL1 pressure data indexes and values
      if (solar.le.vik_tab(1,1)) then
         lsinf=vik_tab(669,1)-360
         lssup=vik_tab(1,1)
         psinf=vik_tab(669,2)
         pssup=vik_tab(1,2)
      else if (solar.ge.vik_tab(669,1)) then
         lsinf=vik_tab(669,1)
         lssup=vik_tab(1,1)+360
         psinf=vik_tab(669,2)
         pssup=vik_tab(1,2)
      else
         iloop=1
         do while ((vik_tab(iloop,1).lt.solar).and.(iloop.lt.669))
            iloop=iloop+1
         enddo
         lsinf=vik_tab(iloop-1,1)
         lssup=vik_tab(iloop,1)+360
         psinf=vik_tab(iloop-1,2)
         pssup=vik_tab(iloop,2)
      endif

C 6. Interpolate VL1 pressure data for solar longitude 'solar'

c iloop-1 is the "lower limit" and iloop the "upper limit" in tab_vik
      presvik=psinf+(solar-lsinf)*(pssup-psinf)/(lssup-lsinf)


C 7. Compute the pressure correction factor
      factcor=presvik/presmean

      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine ps_mcd(dset,latin,lonin,solar,time,
     &     ps,temp,altout,ierr)
! Give surface pressure, temperature at ~1km and altitude extracted
! from Mars Climate Database for a point with latitude, longitude, 
! solar longitude and local time coordinates.
! Note: horizontal bilinear interpolation is used for orography and
!       temperature; but for surface pressure, bilinear interpolation
!       of log(pressure) is used.
!       pressure and temperature are interpolated linearily with respect
!       to encompassing monthes and universal time

      implicit none

      include "netcdf.inc"

! Arguments:
!     inputs:
      character*(*) dset ! Path to MCD datafiles 
      real latin        ! north latitude (degrees)
      real lonin	! east longitude (degrees)
      real solar	! solar longitude (degrees)
      real time 	! local time (martian hours)
!     outputs:
      real ps		! interpolated surface pressure
      real altout	! interpolated altitude
      real temp 	! intepolated Temperature at 50m above ground
      integer ierr	! subroutine status (0 if everything OK)

! Local variables

      character*12 psfile ! datafile which contains the MCD data subset
      data psfile/'ps_clim.nc'/
      character*140 datafile   ! full path+datafile
      integer nid,dimid,nvarid ! NetCDF IDs
      integer iloop 
      integer dimlength     ! to verify that dimensions lenghts are fine
      integer iim, jjm, num ! sizes of MCD data arrays
!      parameter(iim=73,jjm=38,num=13) ! MCD v3 format
      parameter(iim=65,jjm=49,num=13)  ! MCD v4 format
      real lat(jjm)  ! latitudes (range from 90 to -90)
      real deltalat  ! step in latitudes
      parameter (deltalat=180.0/(jjm-1))
      real lon(iim)  ! longitudes (degrees east) (range from 0 to 360)
      real deltalon  ! step in longitudes
      parameter (deltalon=360.0/(iim-1))
      real ls(num)  ! areocentric longitude (range from 15 to 375)
      real deltals  ! step in solar longitude
      parameter (deltals=30.0)
      real utime   ! universal time corresponding to input (local) 'time'
      real ut(num) ! universal time of day (in hours) (range from 0 to 24)
      real deltalt ! step in universal time
      parameter (deltalt=24.0/(num-1))
      real psmcd(iim,jjm,num,num)   ! surface pressure from MCD
      real tempmcd(iim,jjm,num,num) ! atmospheric temperature from MCD
      real altmcd(iim,jjm)          ! orography from MCD
      save lat,lon,ls,ut,psmcd,tempmcd,altmcd

      integer monthinf,monthsup ! encompassing MCD month indexes
      integer timeinf,timesup ! encompassing MCD time indexes
      integer latp(4) ! latitudes of neighbouring points
      integer latinf,latsup ! indexes (in latitude) of neighbouring points
      integer lonp(4) ! longitudes of neighbouring points

      integer loninf,lonsup ! indexes (in longitude) of neighbouring points
      real colat ! colatitude (=90-latin)
      real dlat,dlon
      real solartmp
      logical firstcall
      data firstcall/.true./
      save firstcall

      real alt(4)     ! to store MCD orography
      real ps1(4,2)   ! to store MCD surface pressure at month(monthinf)
      real temp1(4,2) ! to store MCD temperature at month(monthinf)
      real ps1_tinf   ! pressure at time(timeinf) and month(monthinf)
      real ps1_tsup   ! pressure at time(timesup) and month(monthinf)
      real ps1_m1     ! interpolated pressure for month(monthinf)
      real temp1_tinf ! temperature at time(timeinf) and month(monthinf)
      real temp1_tsup ! temperature at time(timesup) and month(monthinf)
      real temp1_m1   ! interpolated temperature for month(monthinf)
      real ps2(4,2)   ! to store MCD surface pressure at month(monthsup)
      real temp2(4,2) ! to store MCD temperature at month(monthsup)
      real ps2_tinf   ! pressure at time(timeinf) and month(monthsup)
      real ps2_tsup   ! pressure at time(timesup) and month(monthsup)
      real ps2_m2     ! interpolated pressure for month(monthsup)
      real temp2_tinf ! temperature at time(timeinf) and month(monthsup)
      real temp2_tsup ! temperature at time(timesup) and month(monthsup)
      real temp2_m2   ! interpolated temperature for month(monthsup)
      
C 1. If first call, load MCD data 
      if (firstcall) then

         firstcall=.false.

         datafile=dset//psfile
         
         ierr=NF_OPEN(datafile,NF_NOWRITE,nid)
         if (ierr.ne.NF_NOERR) then
            write(6,*)'ps_mcd error: could not open file ',
     &                trim(datafile)
            ierr=4
            return
         endif

         ierr=NF_INQ_DIMID(nid,"longitude",dimid)
         dimlength=0
         ierr=NF_INQ_DIMLEN(nid,dimid,dimlength)
         if (dimlength.ne.iim) then
           write(6,*)'ps_mcd error: unexpected longitude lenght !!!'
           write(6,*)'     expected: ',iim,'   found: ',dimlength
           ierr=4
           return
         endif

         ierr=NF_INQ_DIMID(nid,"latitude",dimid)
         dimlength=0
         ierr=NF_INQ_DIMLEN(nid,dimid,dimlength)
         if (dimlength.ne.jjm) then
           write(6,*)'ps_mcd error: unexpected latitude lenght !!!'
           write(6,*)'     expected: ',jjm,'   found: ',dimlength
           ierr=4
           return
         endif

         ierr=NF_INQ_DIMID(nid,"ls",dimid)
         dimlength=0
         ierr=NF_INQ_DIMLEN(nid,dimid,dimlength)
         if (dimlength.ne.num) then
           write(6,*)'ps_mcd error: unexpected ls lenght !!!'
           write(6,*)'     expected: ',num,'   found: ',dimlength
           ierr=4
           return
         endif

         ierr=NF_INQ_DIMID(nid,"ut",dimid)
         dimlength=0
         ierr=NF_INQ_DIMLEN(nid,dimid,dimlength)
         if (dimlength.ne.num) then
           write(6,*)'ps_mcd error: unexpected ut lenght !!!'
           write(6,*)'     expected: ',num,'   found: ',dimlength
           ierr=4
           return
         endif

         ierr = NF_INQ_VARID (nid, "latitude",nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <latitude> not found"
            ierr=4
            return
         endif

         ierr = NF_INQ_VARID (nid, "longitude",nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <longitude> not found"
            ierr=4
            return
         endif

         ierr = NF_INQ_VARID (nid, "ls",nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, ls)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <ls> not found"
            ierr=4
            return
         endif

         ierr = NF_INQ_VARID (nid, "ut",nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, ut)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <ut> not found"
            ierr=4
            return
         endif

         ierr = NF_INQ_VARID (nid, "ps", nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, psmcd)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <ps> not found"
            ierr=4
            return
         endif

         ierr = NF_INQ_VARID (nid, "temp", nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, tempmcd)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <temp> not found"
            ierr=4
            return
         endif

         ierr = NF_INQ_VARID (nid, "alt", nvarid)
         ierr = NF_GET_VAR_REAL(nid, nvarid, altmcd)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"ps_mcd error: <alt> not found"
            ierr=4
            return
         endif

         ierr = NF_CLOSE(nid)

      endif  !  of if(firstcall)

C 2. Interpolate MCD data

C 2.1 Prepare things; find neighbouring points and month and time indexes

C 2.1.1 Find encompassing months indexes

      if (solar.lt.(deltals/2.0)) then
         solartmp=solar+360.0     !for interpolation along Ls
      else
         solartmp=solar
      endif
      monthinf=int(((deltals/2.0)+solartmp)/deltals)
      monthsup=monthinf+1

C 2.1.2 Find encompassing universal time indexes

! compute universal time from local time
      utime=time-lonin/15.0

! universal time should be >=0 and <=24
      if (utime.lt.0.0) then
        utime=utime+24.0
      endif
      utime=mod(utime,24.0)
      
      timeinf=1+int(utime/deltalt)
      timesup=timeinf+1

C 2.1.3 Find the four nearest MCD lon-lat grid points
      ! latitudes
      colat=90-latin

      if (colat.eq.0) then
         latsup=1
         latinf=1
      else if (colat.eq.180) then
         latsup=jjm
         latinf=jjm
      else
         latsup=1+int(colat/deltalat)
         latinf=latsup+1
      endif
      
      ! longitudes
      loninf=1+int(mod(lonin,360.0)/deltalon)
      lonsup=loninf+1
      

! Four nearest points are arranged as follows:  1 2
!                                               4 3
      latp(1)=latsup
      latp(2)=latsup
      latp(3)=latinf
      latp(4)=latinf

      lonp(1)=loninf
      lonp(2)=lonsup
      lonp(3)=lonsup
      lonp(4)=loninf

C 2.2 Interpolations

C 2.2.1 Orography (time independent, so only bilinear interpolation
C                  on lon-lat grid is required) 
      do iloop=1,4
        alt(iloop)=altmcd(lonp(iloop),latp(iloop))
      enddo

C Bilinear interpolation based on the four neighbours
! Note: dlat is the "distance" to lat(latsup)
      if ((colat.eq.0).or.(colat.eq.180)) then
         dlat=0
      else if (colat.lt.deltalat) then
         dlat=colat/deltalat
      else 
         dlat=(lat(latsup)-latin)/deltalat
      endif

! Note: dlon is the "distance" to lon(loninf)
      dlon=(lonin-(loninf-1)*deltalon)/deltalon

      if ((dlat.gt.1).or.(dlon.gt.1).or.(dlat.lt.0).or.(dlon.lt.0)) then
         write(6,*)"ps_mcd error: on dlat or dlon in ps_mcd" 
         write(6,*)"dlat: ",dlat
         write(6,*)"lat: ",latin
         write(6,*)"dlon: ",dlon
         write(6,*)"lon: ",lonin
         ierr=5
         return
      endif

      altout=(1.-dlon)*(1.-dlat)*alt(1)+(1.-dlat)*dlon*alt(2)+
     &        dlat*(1.-dlon)*alt(4)+dlat*dlon*alt(3)

C 2.2.2 Surface pressure and atmospheric temperature

C 2.2.2.1 For month(monthinf)

C Get pressures and temperatures at times timeinf and timesup
      do iloop=1,4
        ps1(iloop,1)=psmcd(lonp(iloop),latp(iloop),timeinf,monthinf)
        ps1(iloop,2)=psmcd(lonp(iloop),latp(iloop),timesup,monthinf)
        temp1(iloop,1)=tempmcd(lonp(iloop),latp(iloop),timeinf,monthinf)
        temp1(iloop,2)=tempmcd(lonp(iloop),latp(iloop),timesup,monthinf)
      enddo

C bilinear interpolations of log(pressure)
      ps1_tinf=(1.-dlon)*(1.-dlat)*log(ps1(1,1))
     &        +(1.-dlat)*dlon*log(ps1(2,1))
     &        +dlat*(1.-dlon)*log(ps1(4,1))
     &        +dlat*dlon*log(ps1(3,1))
      ps1_tinf=exp(ps1_tinf)
      ps1_tsup=(1.-dlon)*(1.-dlat)*log(ps1(1,2))
     &        +(1.-dlat)*dlon*log(ps1(2,2))
     &        +dlat*(1.-dlon)*log(ps1(4,2))
     &        +dlat*dlon*log(ps1(3,2))
      ps1_tsup=exp(ps1_tsup)

C bilinear interpolations for temperature
      temp1_tinf=(1.-dlon)*(1.-dlat)*temp1(1,1)+
     &           (1.-dlat)*dlon*temp1(2,1)+
     &           dlat*(1.-dlon)*temp1(4,1)+
     &           dlat*dlon*temp1(3,1)
      temp1_tsup=(1.-dlon)*(1.-dlat)*temp1(1,2)+
     &           (1.-dlat)*dlon*temp1(2,2)+
     &           dlat*(1.-dlon)*temp1(4,2)+
     &           dlat*dlon*temp1(3,2)

C linear interpolation in universal time
      ps1_m1=ps1_tinf+((utime-ut(timeinf))/deltalt)*(ps1_tsup-ps1_tinf)
      temp1_m1=temp1_tinf+
     &         ((utime-ut(timeinf))/deltalt)*(temp1_tsup-temp1_tinf)

C 2.2.2.2 For month(monthsup)

C Get pressures and temperatures at times timeinf and timesup
      do iloop=1,4
        ps2(iloop,1)=psmcd(lonp(iloop),latp(iloop),timeinf,monthsup)
        ps2(iloop,2)=psmcd(lonp(iloop),latp(iloop),timesup,monthsup)
        temp2(iloop,1)=tempmcd(lonp(iloop),latp(iloop),timeinf,monthsup)
        temp2(iloop,2)=tempmcd(lonp(iloop),latp(iloop),timesup,monthsup)
      enddo

C bilinear interpolation of log(surface pressure)
      ps2_tinf=(1.-dlon)*(1.-dlat)*log(ps2(1,1))+
     &         (1.-dlat)*dlon*log(ps2(2,1))+
     &         dlat*(1.-dlon)*log(ps2(4,1))+
     &         dlat*dlon*log(ps2(3,1))
      ps2_tinf=exp(ps2_tinf)
      ps2_tsup=(1.-dlon)*(1.-dlat)*log(ps2(1,2))+
     &         (1.-dlat)*dlon*log(ps2(2,2))+
     &         dlat*(1.-dlon)*log(ps2(4,2))+
     &         dlat*dlon*log(ps2(3,2))
      ps2_tsup=exp(ps2_tsup)

c bilinear interpolation for temperature
      temp2_tinf=(1.-dlon)*(1.-dlat)*temp2(1,1)+
     &           (1.-dlat)*dlon*temp2(2,1)+
     &           dlat*(1.-dlon)*temp2(4,1)+
     &           dlat*dlon*temp2(3,1)
      temp2_tsup=(1.-dlon)*(1.-dlat)*temp2(1,2)+
     &           (1.-dlat)*dlon*temp2(2,2)+
     &           dlat*(1.-dlon)*temp2(4,2)+
     &           dlat*dlon*temp2(3,2)

C linear interpolation in universal time
      ps2_m2=ps2_tinf+((utime-ut(timeinf))/deltalt)*(ps2_tsup-ps2_tinf)
      temp2_m2=temp2_tinf+
     &         ((utime-ut(timeinf))/deltalt)*(temp2_tsup-temp2_tinf)

C 2.2.2.3 linear interpolation in Ls

      ps=ps1_m1+((solartmp-ls(monthinf))/deltals)*(ps2_m2-ps1_m1)

      temp=temp1_m1+
     &      ((solartmp-ls(monthinf))/deltals)*(temp2_m2-temp1_m1)

      
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine mola(dset,latitude,longitude,alt,ierr)
!Give the MOLA altitude (alt), given lat and lon coordinates
!Using bilinear interpolation from 32 pixels/degree MOLA file

      implicit none

      include "netcdf.inc"

! Arguments
! inputs
      character*(*) dset ! Path to MCD datafiles 
      real latitude ! north latitude (degrees)
      real longitude ! east longitude (degrees)
! outputs
      real alt     ! above areoid altitude of surface (meters)
      integer ierr ! returned status code (==0 if OK)

! Local variables

      logical firstcall
      data firstcall/.true./
      save firstcall
      character*140 datafile ! full path+file name
      character*12 molafile ! MOLA datafile 
c      data molafile/'mola_32.2.nc'/
c      data molafile/'mola16.nc'/
      data molafile/'mola32.nc'/
      real resol
c      parameter(resol=16) ! MOLA pixel/degree resolution
      parameter(resol=32)
      integer jjm, iim   ! # of longitude and latitude MOLA data values
      parameter(jjm=180*resol, iim=2*jjm)
      integer*2 altmola(iim,jjm) ! MOLA topography dataset
      save altmola               ! saved so that it needs to be read only once
      integer nid,nvarid         ! NetCDF file and variable IDs

      integer*2 mintopo_check,maxtopo_check ! known min and max of MOLA dataset
c      parameter(mintopo_check=-8156,maxtopo_check=21191) ! mola_32.2.nc
c      parameter(mintopo_check=-8177,maxtopo_check=21191) ! mola16.nc
      parameter(mintopo_check=-8206,maxtopo_check=21191) ! mola32.nc
      integer*2 mintopo, maxtopo ! min and max of read dataset

      real topo(4) ! neighboring MOLA points (for bilinear interpolation)
      integer latsup,latinf,loninf,lonsup ! indexes of neighboring points
      real dlat    ! normalized "latitudinal distance" to point 3 
      real dlon    ! normalized "longitudinal distance" to point 3
      real colat   ! colatitude
      real lat,lon ! latitude and longitude, local values (degrees)
      integer i,j  ! for loops

C 1. Load MOLA dataset upon first call
C    (dataset is stored in a 'save' array)

      if (firstcall) then
         firstcall=.false.

C 1.1. Open MOLA file
         datafile=dset//molafile
         write(6,*) "Loading MOLA topography from file ",
     &              trim(datafile)

         ierr=NF_OPEN(datafile,NF_NOWRITE,nid)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"Error in mola: Could not open file ",
     &               trim(datafile)
            stop
         endif

C 1.2. Load data
         ierr = NF_INQ_VARID (nid, "alt", nvarid)
	 ! note that MOLA "alt" are given as "short" (16 bits integers)
         ierr = NF_GET_VAR_INT2(nid, nvarid, altmola)
         if (ierr.ne.NF_NOERR) then
            write(6,*)"Error in mola: <alt> not found"
            stop
         endif

c 1.3. Close MOLA file
         ierr = NF_CLOSE(nid)

c 1.4 Check that the MOLA dataset was correctly loaded

         mintopo=mintopo_check
         maxtopo=maxtopo_check
         do i=1,iim
            do j=1,jjm
               mintopo=min(mintopo,altmola(i,j))
               maxtopo=max(maxtopo,altmola(i,j))
            enddo
         enddo
         if ((mintopo.ne.mintopo_check).or.
     &       (maxtopo.ne.maxtopo_check)) then
            write(6,*)"Error in mola:"
            write(6,*)" Mola file ",trim(molafile)," was not well read"
            write(6,*)" Minimum found: ", mintopo
            write(6,*)" Minimum should be:",mintopo_check
            write(6,*)" Maximum found: ", maxtopo
            write(6,*)" Maximum should be:",maxtopo_check
            ierr=2
            return
         endif
         write(6,*)"Done loading MOLA topography"

      endif ! End of if(firstcall)

C 2. Check that input longitude and latitude make sense
      lat=latitude
      if((lat.gt.90).or.(lat.lt.-90)) then
        write(6,*)"Error in mola: Wrong value for latitude"
        stop
      endif

c longitude must range from 0 to 360
      lon=longitude
      do while(lon.gt.360.)
        lon=lon-360.
      enddo
      
      do while(lon.lt.0.)
        lon=lon+360.
      enddo
      
C 3. Identify the four neighboring points from MOLA dataset
!    These points are arranged as follows:  1 2
!                                           3 4

      colat=90-lat

      if (colat.lt.1./(2.*resol)) then
         latsup=1 
         latinf=1
         dlat=0
      else if (colat.gt.180.-1./(2.*resol)) then
         latsup=jjm
         latinf=jjm
         dlat=0
      else
         latsup=1+int((colat-1./(2.*resol))*resol)
         latinf=latsup+1
         dlat=1-(colat-(1./(2.*resol)+(latsup-1)/resol))*resol
      endif
! Note: dlat is the (normalized) "latitudinal distance" to point 3 
!       ie: dlat=0 if lat=latitude of point 3
!           dlat=1 if lat=latitude of point 1

      if ((lon.lt.1./(2.*resol)).or.(lon.ge.(360-1./(2.*resol)))) then
         loninf=iim
         lonsup=1
         if (lon.lt.1./(2.*resol)) then
            dlon=lon*resol+0.5
         else
            dlon=(lon-(1./(2.*resol)+(loninf-1)/resol))*resol
         endif
      else
         if (((lon-1./(2.*resol))*resol).ge.0) then
            loninf=1+int((lon-1./(2.*resol))*resol)
         else
            loninf=int((lon-1./(2.*resol))*resol)
         endif
         lonsup=loninf+1
         dlon=(lon-(1./(2.*resol)+(loninf-1)/resol))*resol
      endif
! Note: dlon is the (normalized) "longitudinal distance" to point 3
!       ie: dlon=0 if lon=longitude of point 3
!           dlon=1 if lon=longitude of point 4

      if((dlat.gt.1).or.(dlon.gt.1).or.(dlat.lt.0).or.(dlon.lt.0)) then
         write(6,*)"Error in mola: on dlat or dlon" 
         write(6,*)"dlat: ",dlat
         write(6,*)"lat: ",lat
         write(6,*)"dlon: ",dlon
         write(6,*)"lon: ",lon
         ierr=2
         return
      endif

C 4. Interpolation

! Four nearest points are arranged as follows:  1 2
!                                               3 4
      topo(1)=real(altmola(loninf,latsup))
      topo(2)=real(altmola(lonsup,latsup))
      topo(3)=real(altmola(loninf,latinf))
      topo(4)=real(altmola(lonsup,latinf))

! Use bilinear interpolation to evaluate alt

      alt=(1.-dlon)*(1.-dlat)*topo(3)+(1.-dlat)*dlon*topo(4)
     &     +dlat*(1.-dlon)*topo(1)+dlat*dlon*topo(2)

      end

