      subroutine call_mcd(zkey,xz,xlon,xlat,hireskey,
     &      datekey,xdate,localtime,dset,scena,
     &      perturkey,seedin,gwlength,extvarkeys,
     &      pres,dens,temp,zonwind,merwind,
     &      meanvar,extvar,seedout,ier)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Purpose:
c =======
c     call_mcd is a fortran subroutine which extracts and computes
c     meteorogical variables from the Mars Climate Database
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Version:
c =======
c     5.3 07/17 New EUV scenarios and improved gravity waves perturbations,
c               Improved extrapolation above ~250km
c               And added He outputs and a Mars Year 32 scenario
c     5.2 02/15 Added scenarios 24 to 31 for Mars Years 24 to 31
c     5.2 02/15 MCDv5.2 added output of solar zenit angle, Sun-Mars distance,
c               local mean solar time; improved surface layer computations
c     5.2 01/15 MCDv5.2 release with corrected vertical velocity and TEC
c     5.1 05/14 MCDv5.1 release ; improved dust storm high res. with "VL3"
c     5.0 03/14 added near-surface noise perturbation model 
c     5.0 12/12 new scenarios and more outputs and improved EOF perturbations
c     4.3 04/08 new improved large scale (EOF) perturbations model
c     4.2 02/08 improved GW perturbations behaviour near surface (grwpb.F)
c     4.2 02/08 fixed 2 minor bugs (random seed initialization & getsi.F)
c     4.2 02/08 implemented EOF perturbations interpolations
c     4.2 02/08 added zkey=5 option 
c     4.2 06/07 various minor bug fixes
c     4.2 03/07 added interpolation for temp. and winds below first layer
c     4.2 03/07 added altitude-wise RMS
c     4.2 12/06 changed "atmemcd" routine to "call_mcd"; added high resolution
c     4.1 03/06 improved time routines (orbit,mars_ltime,sol2ls,ls2sol)
c     4.1 02/05 add many new output variables.  
c     4.0 12/04 add many new output variables.  
c     4.0 09/04 height calculation in atmemcd (not in test_emcd anymore)      
c     4.0 09/04 choice of vertical coordinates, input output 
c     4.0 09/04 R changes with altitude 
c     4.0 06/04 add variables, add RMS, up and low atmosphere
c     3.2 03/04 saved variables, seasonal interpolation, variable gravity
c     3.1 05/01 bug fixes and improved interpolations
c     3.0 02/01 - Major MCD release based on earlier code
c     2.3 05/09/00
c     2.2 23/03/00
c
c Authors:
c =======
c     v5.0 Implemented new scenarios and possibility to select which
c          extra variables will be computed -- EM
c     v4.3 Implemented improved large scale perturbations model -- EM
c     v4.2 Changed argument list; implemented high resolution and altitude-wise
c          RMS, and many other small changes  -- EM
c     v4.1 Improved Earth to Mars date conversion + local time -- EM
c     v4.1 Final version with further new variables, bug fixed, EOF -- FF
c     v4.0 version with new variables, bug fixed, EOF -- FF
c     v4.0 combining low and up atmosphere, variables added --
c          min, mean and max solar scenario for thermosphere --
c          choice of vertical coordinates, flag zkey, zareoid, zsurface, 
c          zradius,zpressure in input output --       
c          altitude in meters everywhere -- KD, FF
c     v3.2 Variables saved in subroutines, seasonal interpolation added -- SJB
c          Gravity allowed to vary with height in hydrostatic eqn  -- SRL
c
c     v3.1 Bugs fixed, new ls2sol routine and various improvements
c     to interpolation above model top, especially for EOFs -- SRL
c
c     Major updates for v3.0 MCD -- SRL + FF
c       See Changes file for more information.
c
c     Modifications for v2.3 MCD -- SRL
c
c     Derived from code for version 2 by C. HOURTOLLE
c
c
c Arguments (inputs):
c ==================
c   zkey  : <integer>   type of vertical coordinate xz
c                       1 = radius from centre of planet (m)
c                       2 = height above areoid (m) (MOLA zero datum)
c                       3 = height above surface (m)
c                       4 = pressure level (Pa)
c                       5 = altitude above mean Mars Radius(=3396000m) (m)
c   xz    : <real> vertical coordinate (m or Pa, depends on zkey)
c   xlon  : <real> longitude (degrees east)
c   xlat  : <real> latitude (degrees north)
c   hireskey: <integer> flag to switch to high resolution topography
c                    0 = use GCM resolution
c                    1 = switch to high resolution topography
c   datekey: <integer>    type of input date
c                    0 = "Earth time": xdate is given in Julian days
c                                      (localtime must be set to zero)
c                    1 = "Mars date": xdate is the value of Ls
c   xdate : <double precision> date 
c                      IF datekey = 0 : Earth julian date
c                      IF datekey = 1 : Value of Ls (in degrees)
c   localtime : <real> local time (in martian hours) at longitude xlon
c              ONLY USED IF datekey=1 (must be set to 0 if datekey=0)
c   dset  : <character*(*)>    data set
c                              One or more blanks to get the default,
c                              or full directory path of MCD data required
c                              including trailing slash (e.g. '/dir/path/'),
c                              or a link to a directory (e.g. 'MCD_DATA/').
c                              Default is link MCD_DATA in working directory.
c   scena : <integer>          scenario
c                             1 = Climatology ave solar
c                             2 = Climatology min solar
c                             3 = Climatology max solar
c                             4 = dust storm tau=5 (dark dust) min solar           
c                             5 = dust storm tau=5 (dark dust) ave solar           
c                             6 = dust storm tau=5 (dark dust) max solar           
c                             7 = warm scenario - dusty, max solar        
c                             8 = cold scenario - low dust, min solar
c                            24 = Mars Year 24, with associated solar EUV
c                            25 = Mars Year 25, with associated solar EUV
c                            26 = Mars Year 26, with associated solar EUV
c                            27 = Mars Year 27, with associated solar EUV
c                            28 = Mars Year 28, with associated solar EUV
c                            29 = Mars Year 29, with associated solar EUV
c                            30 = Mars Year 30, with associated solar EUV
c                            31 = Mars Year 31, with associated solar EUV
c                            32 = Mars Year 32, with associated solar EUV
c   perturkey: <integer>  perturbation type
c             1: none
c             2: large scale : 'seedin' is seed or signals reseting the
c                EOF perturbations if changed between calls to CALL_MCD
c             3: small scale : 'seedin' is seed or signals reseting the
c                GW perturbations if changed between calls to CALL_MCD
c             4: large and small scale : does both 2 and 3 above
c             5: add 'seedin' times the standard deviation
c                (seedin must not be greater than 4 or less than -4)
c   seedin: <real> 
c           IF perturkey=2,3,4: seed for random number generation
c                               changes in seedin between subsequent calls
c                               trigger reseeding and regeneration of
c                               perturbations
c           IF perturkey=5: coefficient by which standard deviations should
c                           be multiplied before being added to mean values
c   gwlength : <real>  for small scale (ie: gravity wave) perturbations;
c                   gwlength= vertical wavelength (lamda) of gravity wave perturbation (m)
c                             set to zero to get default (16000m). 
c                             horizontal wavelength (lamda_D) is proportional to the
c                             vertical one (lamda_D = 10 x lamda)
c   extvarkeys  : <integer> array output type
c                   extvar(i) = 0 : don't compute and output extvar(i)
c                   extvar(i) = 1 : compute and output extvar(i)
c
c Arguments (outputs):
c ===================
c   pres    : <real> atmospheric pressure (Pa)
c   dens    : <real> atmospheric density (kg/m^3)
c   temp    : <real> atmospheric temperature (K)
c   zonwind : <real> zonal wind component (East-West)
c   merwind : <real> meridional wind component (North-South)
c   meanvar : <real> mean unperturbed values (array of dimension nmeanvar=5)
c               meanvar(1)= mean pressure
c               meanvar(2)= mean density
c               meanvar(3)= mean temperature
c               meanvar(4)= mean zonal wind component
c               meanvar(5)= mean meridional wind component
c   extvar : <real>  extra variables array (of dimension nextvar=100)
c               extvar(1) = Radial distance from planet center (m)
c               extvar(2) = Altitude above areoid (Mars geoid) (m)
c               extvar(3) = Altitude above local surface (m)
c               extvar(4) = orographic height (m) (surface altitude above areoid)
c               extvar(5) = Ls, solar longitude of Mars (deg)
c               extvar(6) = LST local true solar time (hrs)
c               extvar(7) = Universal solar time (LST at lon=0) (hrs)
c               extvar(8) = Air heat capacity Cp (J kg-1 K-1)
c               extvar(9) = gamma=Cp/Cv Ratio of specific heats
c               extvar(10)= density RMS day to day variations (kg/m^3)
c               extvar(11)= LMT: Local mean time (hrs) at sought longitude
c                           (LMT is only computed for Earth date input)
c               extvar(12)= Sun-Mars distance (in Astronomical Unit AU)
c               extvar(13)= scale height H(p) (m)
c               extvar(14)= GCM orography (m)
c               extvar(15)= surface temperature (K)
c               extvar(16)= daily maximum mean surface temperature (K)
c               extvar(17)= daily minimum mean surface temperature (K)
c               extvar(18)= surf. temperature RMS day to day variations (K) 
c               extvar(19)= surface pressure (high resolution if hireskey=1)
c               extvar(20)= GCM surface pressure (Pa)
c               extvar(21)= atmospheric pressure RMS day to day variations (Pa)
c               extvar(22)= surface pressure RMS day to day variations (Pa)
c               extvar(23)= temperature RMS day to day variations (K)
c               extvar(24)= zonal wind RMS day to day variations (m/s)
c               extvar(25)= meridional wind RMS day to day variations (m/s)
c               extvar(26)= vertical wind component (m/s) >0 when downwards!
c               extvar(27)= vertical wind RMS day to day variations (m/s)
c               extvar(28)= small scale perturbation (gravity wave) (kg/m^3)
c               extvar(29)= Surface roughness length z0 (m)
c               extvar(30)= solar flux reflected to space (W/m2)
c               extvar(31)= thermal IR flux to surface (W/m2)
c               extvar(32)= solar flux to surface (W/m2)
c               extvar(33)= thermal IR flux to space (W/m2)
c               extvar(34)= surface H2O (seasonal frost) layer (kg/m2)
c               extvar(35)= surface CO2 ice layer (kg/m2)
c               extvar(36)= DOD: Dust column visible optical depth above surface
c               extvar(37)= DOD RMS day to day variations 
c               extvar(38)= Dust mass mixing ratio (kg/kg)
c               extvar(39)= Dust effective radius (m)
c               extvar(40)= Dust deposition on flat surface (kg m-2 s-1)
c               extvar(41)= Water vapor column (kg/m2)
c               extvar(42)= Water vapor vol. mixing ratio (mol/mol)
c               extvar(43)= Water ice column (kg/m2)
c               extvar(44)= Water ice mixing ratio (mol/mol)
c               extvar(45)= Water ice effective radius (m)
c               extvar(46)= Convective Planetary Boundary Layer (PBL) height (m)
c               extvar(47)= Max. upward convective wind within the PBL (m/s)
c               extvar(48)= Max. downward convective wind within the PBL (m/s)
c               extvar(49)= Convective vertical wind variance at level z (m2/s2)
c               extvar(50)= Convective eddy vertical heat flux at level z (m/s/K)
c               extvar(51)= Surface wind stress (Kg/m/s2)
c               extvar(52)= Surface sensible heat flux (W/m2) (<0 when flux from surf to atm.)
c               extvar(53)= R: Molecular gas constant (J K-1 kg-1)
c               extvar(54)= Air viscosity estimation (N s m-2)
c               extvar(55)= not used (set to zero)
c               extvar(56)= Solar zenith angle (deg)
c               extvar(57)= [CO2] volume mixing ratio (mol/mol)
c               extvar(58)= [N2] volume mixing ratio (mol/mol)
c               extvar(59)= [Ar] volume mixing ratio (mol/mol)
c               extvar(60)= [CO] volume mixing ratio (mol/mol)
c               extvar(61)= [O] volume mixing ratio (mol/mol)
c               extvar(62)= [O2] volume mixing ratio (mol/mol)
c               extvar(63)= [O3] volume mixing ratio (mol/mol)
c               extvar(64)= [H] volume mixing ratio (mol/mol)
c               extvar(65)= [H2] volume mixing ratio (mol/mol)
c               extvar(66)= electron number density (particules/cm3)
c               extvar(67)= CO2 column (kg/m2)
c               extvar(68)= N2 column (kg/m2)
c               extvar(69)= Ar column (kg/m2)
c               extvar(70)= CO column (kg/m2)
c               extvar(71)= O column (kg/m2)
c               extvar(72)= O2 column (kg/m2)
c               extvar(73)= O3 column (kg/m2)
c               extvar(74)= H column (kg/m2)
c               extvar(75)= H2 column (kg/m2)
c               extvar(76)= total electonic content (TEC) (particules/m2)
c               extvar(77)= He column (kg/m2)
c               extvar(78)= [He] volume mixing ratio (mol/mol)
c               extvar(79:100) = unused (set to zero)
c   seedout : <real> current value of the seed of the random number generator
c    ier  : <integer> error flag
c     0 = OK
c     1 = wrong vertical coordinate flag (zkey)
c     2 = wrong dust scenario (scena)
c     3 = wrong value for perturbation flag (perturkey)
c     4 = wrong value for high resolution flag (hireskey)
c     5 = wrong value for date flag (datekey)
c     6 = wrong value for extra variables flag (extvarkey)
c     7 = wrong value for latitude (xlat)
c     8 = inadequate value for gravity wave wavelength (gwlength)
c     9 = wrong value of solar longitude (xdate, if datekey=1)
c    10 = Julian date lies outside [jdate_min:jdate_max] range
c    11 = wrong value of local time (should be in [0:24])
c    12 = Incompatible localtime(.ne.0) and datekey(=0)
c    13 = Unresonable value of 'seedin' (in perturkey=5 case)
c    14 = No dust storm scenario available at such date
c    15 = Could not open a database file (wrong path to database?)
c    16 = Failed loading data from a database file
c    17 = Given (or computed) altitude is underground
c    18 = adding (perturkey=5) perturbation yields unphysical density
c    19 = adding (perturkey=5) perturbation yields unphysical temperature
c    20 = adding (perturkey=5) perturbation yields unphysical pressure
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      include "constants_mcd.inc"

c -------------------------------------------
c     Additional Option (not in the input) / Obsolete as of version 4.2
c     itimint, flag for seasonal interpolation
c     itimint 0: no seasonal interpolation
c     itimint 1: seasonal interpolation
      integer,parameter :: itimint=1
c -------------------------------------------
      integer,parameter :: nextvar=100  ! size of extra variable array extvar()
      integer,parameter :: nmeanvar=5 ! size of meanvar() array
      
c     inputs
c     ******      
      integer,intent(in) :: zkey ! flag for vertical coordinate type
      real,intent(in) :: xz   ! vertical coordinate (m or Pa, depending on zkey)
      real,intent(in) :: xlon ! east longitude (degrees) 
      real,intent(in) :: xlat ! latitude (degrees)
      integer,intent(in) :: hireskey  ! flag: high resolution if =1
      integer,intent(in) :: datekey   ! flag: 0=earth date 1= Mars date
      real*8,intent(in) :: xdate     ! earth julian date (or Ls)
      real,intent(in) :: localtime ! true solar time at longitude lon
! Note: localtime can only be imposed if datekey=1 and should be 0 otherwise
      integer,intent(in) :: perturkey ! flag: perturbation type
      real,intent(in) :: seedin    ! seed for random number generation
! seedin is also used to trigger a reinitialization of perturbation
! (triggered if value of seedin changes between subsequent calls to CALL_EMCD)
      real,intent(in) :: gwlength ! gravity wave perturbation (vertical wavelength of)
      integer,intent(in) :: scena     ! dust scenario
      integer,intent(in) :: extvarkeys(nextvar) ! flag: compute extra variable i
                                         !       if i==1
      character*(*),intent(in) :: dset      ! path to datafiles

c     outputs
c     *******      
      real,intent(out) :: meanvar(nmeanvar) ! array for mean values
      real,intent(out) :: extvar(nextvar)   ! array for extra outputs
      real,intent(out) :: seedout ! current value of 'seed' (used by ran1)
      real,intent(out) :: pres    ! atmospheric pressure
      real,intent(out) :: dens    ! atmospheric density
      real,intent(out) :: temp    ! atmospheric temperature
      real,intent(out) :: zonwind ! zonal (eastward) wind
      real,intent(out) :: merwind ! meridional (northward) wind
      integer,intent(out) :: ier     ! status (error) code

c     local variables
c     *************** 
c     NETCDF file IDs :
c     orographic data
      integer,save :: unet
c     mean fields      
      integer,save :: unetm
c     std dev and rms fields
      integer,save :: unetsd
c     mean up fields (thermosphere)
      integer,save :: unetm_up
c     std dev up and rms fields (thermosphere)
      integer,save :: unetsd_up
c     eof fields
      integer,save :: ueof
c     season2 for seasonal interpolation- mean fields   
      integer,save :: unetm2
c     season2 for seasonal interpolation- std dev fields
      integer,save :: unetsd2
c     season2 for seasonal interpolation- mean fields - up 
      integer,save :: unetm2_up
c     season2 for seasonal interpolation- std dev fields - up
      integer,save :: unetsd2_up
c     dust scenario
      integer           dust
      integer          xdust ! for a second scenario (year-to-year interp.) 
c     seed (and also current index) for random number generator
      integer,save :: seed
!      real              R        ! not used any more, see Rnew
!      parameter        (R=191.2)
      real              Rnew
      real              R_gcm(dimlevs)    ! R at GCM levels
      real              rho_gcm(dimlevs)  ! rho at GCM levels
      real              temp_gcm(dimlevs) ! temperature at GCM levels
c     number of standard deviations to add
      real              nbsig
c   Julian dates that should bracket "xdate" (if given as a Juilan date)
      real*8,parameter :: jdate_min=2378496.5 ! 01-01-1800 00:00
      real*8,parameter :: jdate_max=2524593.5 ! 01-01-2200 00:00
      integer       i,l    ! for loops
      integer       ierr   ! to store return codes from called routines
      integer       levlow ! database level, lower bound (wrt sought alt.)
      integer       levhi  ! database level, upper bound (wrt sought alt.)
      real          areoid     ! distance to center of Mars of reference areoid
      real          areoid_gcm ! GCM areoid (m)
      real          areoid_hr  ! High resolution areoid (m)
      real          height    ! height above areoid (m) 
      real          oroheight ! orographic height (m)
      real          oro_hr    ! MOLA orographic height (m)
      real          oro_gcm   ! GCM orographic height (m)
      real          absheight ! height above surface (m)
      real          ls        ! Solar longitude (deg.)
      real          marsau    ! Sun-Mars distance (in AU)
      real          modelday  ! GCM day [0:668.6]
      real          solzenang ! solar zenith angle (deg)
      real          loctime   ! true solar time at longitude lon
      real          utime     ! universal time (0. to 24. hrs)
                              ! =local true solar time at lon=0
      real          lmeantime ! local meant solar time (01 to 24. hrs) at lon
      real          lon      ! longitude east [-180  180]
      real          lat      ! latitude (degrees)
      real          ps       ! surface pressure
      real          ps_gcm   ! surface pressure (as read from MCD)
      real          ps_hr    ! surface pressure (high res, using MOLA topo.)
      real          ps_psgcm ! high res to GCM surface pressure ratio
      real          p_pgcm(dimlevs) ! high res over GCM pressure ratios
c     variables  (trailing l indicates lower level variable from database)
      real          t,tl     ! atmospheric temperature
      real          p,pl     ! atmospheric pressure
      real          rho,rhol ! density
      real          u,ul,v,vl,w_l ! zonal, meridional and vertical winds
      real          col_h2ovapor,col_h2oice
      real          col_co2,col_co
      real          col_o,col_o2,col_o3
      real          col_h,col_h2
!      real          col_h2o2
      real          col_n2,col_ar
      real          col_elec, col_he
      real          vmr_h2oice,vmr_h2o
      real          vmr_co2,vmr_co
      real          vmr_o,vmr_o2,vmr_o3
      real          vmr_h,vmr_h2
!      real          vmr_h2o2
      real          vmr_n2,vmr_ar,vmr_elec,vmr_he
      real          pertm,pertr,pertrho,pertrhogw
!      real          pertrhoeof
!      data          pertrhoeof /0.0/
      real          pertps
      real          pertu, pertv, pertt
      real          levweight(2) ! weights for vertical interpolation
      real          pratio ! 'correction' parameter (usually = 1.0, unless
                    ! outside of sigma range; it then accounts for
                    ! being above/below lowest/highest sigma level)
      real          lamda   !gravity wave vertical wavelength (m) 
!      real          rdnos(dimnevecs) ! deviates for EOF perturbations
!      save          rdnos
      real, save :: rdeof ! uniform random number in [0,1] for EOF perturb.
      real, save :: dev  ! (random) phase angle for gravity waves
      real, save :: rangw1,rangw2 !random number for gravity wave spread
      real          pscaleheight    ! scale height
      real          tsurf,tsurfmax,tsurfmin
      real          rmstsurf ! RMS of surface temperature
      real          rmsps    ! RMS of surface pressure
      real          co2ice ! CO2 ice at the surface
      real          dod, rmsdod,dust_mmr !,tsddod
      real          dust_reff ! dust effective radius
      real          dust_dep ! dust deposition on flat surface kg/m2/s
      real          profile_reff(dimlevs) ! dust effective radius profile
      real          profile_mmr(dimlevs) ! dust mass mlixing ratio profile
      real          surf_h2o_ice ! surface H2O (seasonal deposit) (kg/m2)
      real          h2oice_reff ! water ice effective radius (m)
!      real          tauh2oice ! water ice cloud opacity (abs. at 825 cm-1)
      real          Cp, gamma, viscosity,Rgas
      real          fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw
      real          rmst, rmsu, rmsv, rmsw
      real          rmsrho    ! rms of density at a given pressure level
      real          altrmsp   ! altitude-wise rms of atmospheric pressure
      real          tmeanl(nmeanvar) ! to temporarily store meanvar()
      real          zradius   ! distance to center of planet
      real          zareoid   ! height above areoid
      real          zsurface  ! height above surface
      real          zpressure ! atmospheric pressure
      real          zmradius  ! altitude above mean Mars radius (3.390E6 m)
      real          sheight(dimlevs) ! altitudes (above surface) of GCM
                                     ! sigma levels
      real          surfstress ! surface stress
      real          sensib_flux ! sensible heat flux
      real          vvv ! vertical velocity variance in PBL
      real          vhf ! vertical eddy heat flux in PBL
      real          wstar ! free convection velocity
      real          zmax ! PBL height
      real,parameter :: g0=3.7257964  ! reference gravitational acceleration
      real,parameter :: a0=3396.E3  ! reference distance at which g=g0

      character*16  name
      character*4   typevar
      character*255 dataset
      integer       lendataset
      character*255 dataset2
c  season numbers
      integer       nums,nums2 ! encompassing 'lower' and 'higher' month #
      integer,save :: numsprec,nums2prec ! previous values of nums & nums2
      integer,save :: numsprec_sd,numsprec_gw
c  seasonal interpolation weights 
      real          wl,wh
c  previous dust scenario
      integer,save :: dustprec 
c  previous wavelength of gravity wave perturbation
      real,save :: prevgwlength
c previous value of seedin
      real,save :: prevseedin
c Large scale EOF perturbed fields
      real          ps_gcm_pert   ! perturbed GCM surface pressure
      real          ps_hr_pert    ! perturbed high res surface pressure
      real          temp_gcm_pert(dimlevs) ! perturbed temperature profile
c Near-surface noise (on Ps, Ts and atmospheric temperatures)
      real,save :: ps_noise ! 'noise' to add to surface pressure
      real,save :: temp_noise ! 'noise' to add to temperature
      real,save :: temp_gcm_noise(dimlevs) ! 'noise' to add to temperatures
                                           ! at GCM levels
      real,save :: ps_noise_dev ! deviate to generate ps_noise
      real,save ::  temp_gcm_noise_dev ! deviate to generate temp_noise
c     function declarations
      real          gasdev,ran1,ls2sol
c     vertical coordinates
      real,save :: aps(dimlevs) ! hybrid coordinate 
      real,save :: bps(dimlevs)       ! hybrid coordinate
      real          sigma(dimlevs)     ! sigma levels to work with
      real          sigma_gcm(dimlevs) ! sigma levels in GCM
      real          sigma_hr(dimlevs)  ! high resolution sigma levels
      real,save :: pseudoalt(dimlevs) ! pseudo altitude
c flag for first call to call_mcd
      logical,save :: firstcall=.true.
c flag to (internally) set EOF coordinates-related interpolation (or not)
      logical,save :: inicoordeof ! see eofpb
c flag to interally add near-surface noise with perturbations
      logical,parameter :: nearsurfnoise=.true.
      logical       air_prop_flag
      logical       rmsoutput
c intermediate variables for large scale (EOF) density perturbation
!      real  lowpress,hipress ! pressures at levhi and levlow
!      real  lowdens,hidens   ! densities at levhi and levlow
!      integer oldlevlow,oldlevhi ! 'unperturbed' levlow and levhi values
!      real  oldlevweight2
!      real  newrhol
      real col_gcm(dimlevs) ! number of mol. between levels in GCM
      real col_hr(dimlevs)  ! number of mol. between levels in hires
      real col_hr2gcm       ! ratio of columns in hires to GCM
      integer flag_ini_profi ! flag of initialization in profi()


c 0. Initializations
      if (firstcall) then ! only for very first call to CALL_MCD
        ! set dummy previous value of seedin
        if (seedin.eq.0) then
          prevseedin=-1.0
        else
          prevseedin=-seedin
        endif
        ! set dummy previous value of gwlength
        if (prevgwlength.eq.0) then
          prevgwlength=1.0
        else
          prevgwlength=-gwlength
        endif
        ! set dummy previous values of months
        numsprec=0
        nums2prec=0
        numsprec_sd=0
        numsprec_gw=0
        ! set dummy previous value of dust scenario
        dustprec=0
        ! set firstcall to false
        firstcall=.false.
      endif

      ! initializations for every call to CALL_MCD
      ier=0
      inicoordeof=.true.
      flag_ini_profi = 0 ! No initialization in profi()

c***********************************************************************
c  1. Check (and eventually convert) input arguments
c***********************************************************************

c  1.1 Path to data set
c     find last non-blank character and length of 'dset' srting
      lendataset=len_trim(dset)
      if (lendataset.eq.0) then  !default data set case
         dataset = 'MCD_DATA/'
         dataset2 = 'MCD_DATA/'
         lendataset=len_trim(dataset)
      else  !symbolic link or full path to datasets is given explicitly
         dataset = dset
         dataset2 = dset
      endif

c  1.2 Check value of dust scenario
      dust=scena
      if ((dust.lt.1).or.((dust.gt.8).and.(dust.lt.24)).or.
     &     (dust.gt.32)) then
        if (output_messages) then
         write(out,*)'CALL_MCD Error: ',scena,' unknown scenario'
        endif
        ier=2
        goto 9999
      endif

c  1.3 Check vertical coordinate type (and set high_res_topo flag)
      if ((zkey.gt.5).or.(zkey.le.0)) then
        if (output_messages) then
          write(out,*)'CALL_MCD Error: zkey=',zkey,
     &                ' unknown vert. coord'
        endif
        ier=1
        goto 9999
      endif


c 1.4 Check various flags
c 'hireskey' flag:
      if ((hireskey.lt.0).or.(hireskey.gt.1)) then
        if (output_messages) then
          write(out,*)'CALL_MCD Error: wrong value for ',
     &                'parameter hireskey'
          write(out,*)'               hireskey=',hireskey
        endif
        ier=4
        goto 9999
      endif

c 'datekey' flag:
      if ((datekey.lt.0).or.(datekey.gt.1)) then
        if (output_messages) then
          write(out,*)'CALL_MCD Error: wrong value for ',
     &                'parameter datekey'
          write(out,*)'               datekey=',datekey
        endif
        ier=5
        goto 9999
      endif

c 'perturkey' flag:
      if ((perturkey.lt.0).or.(perturkey.gt.5)) then
        if (output_messages) then
         write(out,*)'CALL_MCD Error: wrong value for ',
     &               'parameter perturkey'
         write(out,*)'               perturkey=',perturkey
        endif
        ier=3
        goto 9999
      endif

c 'extvarkeys' flag:
      do i=1,nextvar
      if ((extvarkeys(i).lt.0).or.(extvarkeys(i).gt.1)) then
       if (output_messages) then
        write(out,*)'CALL_MCD Error: wrong value for ',
     &              'parameter extvarkeys(',i,')'
        write(out,*)'          extvarkeys(',i,')=',extvarkeys(i)
       endif
       ier=6
       goto 9999
      endif
      enddo

c check that the value of gwlength is resonable
      if (((gwlength.lt.2.E3).and.(gwlength.ne.0.)).or.
     &     (gwlength.gt.30.E3)) then
        if (output_messages) then
         write(out,*)'CALL_MCD Error: wrong value for ',
     &               'parameter gwlength'
         write(out,*)'               gwlength=',gwlength
         write(out,*)'(should be in [2000:30000] range)'
        endif
        ier=8
        goto 9999
      endif

c check that the value of 'seedin' is reasonable (in perturkey=5 case only)
      if ((perturkey.eq.5).and.(abs(seedin).gt.4.0)) then
       if (output_messages) then
        write(out,*)'CALL_MCD Error: wrong value for parameter seedin'
        write(out,*)'               seedin=',seedin
        write(out,*)'(should not be more/less ',
     &              'than +/-4 when perturkey=5)'
       endif
       ier=13
       goto 9999
      endif
      
c 1.4 Check and set latitude and longitude
      if (abs(xlat).gt.90.0) then
        if (output_messages) then
         write(out,*)'CALL_MCD Error: wrong value for latitude=',xlat
        endif
        ier=7
        goto 9999
      else
        lat=xlat
      endif
       
      ! we want longitude to be in [-180:180]
      lon=mod(xlon,360.0)
      if (lon.lt.-180.) then ! in case lon in [-360:-180]
        lon=lon+360.
      endif
      if (lon.gt.180.0) then ! in case lon in [180:360]
        lon=lon-360.
      endif

c 1.5 Perturbations

c 1.5.1 (re-)initialize EOF and GW perturbations, if instructed to do so
! Note: since on first call seedin.ne.prevseedin, the random number generator
! will always be seeded, which is good, since even when no perturbations
! are added, some extra variables need random numbers
      if (seedin.ne.prevseedin) then
        ! generate the seed for random number generation
        ! (seed must then be a strictly negative integer)
        seed=-abs(int(seedin))
        if (seed.eq.0) then
          seed=-1
        endif
        ! compute rdnos() for EOFs
!        do i=1,dimnevecs
!          rdnos(i)=gasdev(seed)
!        enddo
        ! compute rdeof for EOFs
        rdeof=ran1(seed)
        ! compute dev for GW
        dev=ran1(seed)
        rangw1=ran1(seed)
        rangw2=ran1(seed)
        ! compute gaussian deviates for near-surface perturbations
        temp_gcm_noise_dev=gasdev(seed)
        ps_noise_dev=gasdev(seed)
        ! store current seed in seedout
        seedout=seed
        ! set prevseedin to seedin and prevgwlength to gwlength
        prevseedin=seedin
        prevgwlength=gwlength
      endif

c 1.5.2 re-initialize only GW perturbations, if instructed to do so
      if ((gwlength.ne.prevgwlength).and.
     &    ((perturkey.eq.3).or.(perturkey.eq.4))) then
        ! the random number generator necessarily has already been seeded
        ! get a new random phase for the gravity waves
        dev=ran1(seed)
        ! store current value of seed in seedout
        seedout=seed
        ! set prevgwlength to gwlength
        prevgwlength=gwlength
      endif

c 1.5.3 Check and set GW wavelength, if necessary
      if ((perturkey.eq.3).or.(perturkey.eq.4) ! GW perturbations
!     &     .or.(extvarkey.eq.1)) then ! some extra variables need 'gwlength'
     &   ) then
        ! set wavelength 'lamda'
        if (gwlength.eq.0.0) then ! use a default wavelength
          lamda=16.E3
        else
          lamda=gwlength
          ! Note: it was checked that 2000<gwlength<30000 in 1.4
        endif
      endif

c 1.5.4 Check value of 'seedin' for the 'n times standard dev' perturbation
      nbsig=0 ! dummy initialization (to get rid of compiler warnings)
      if (perturkey.eq.5) then
        ! check that -4<seedin<+4 was done in 1.4
        nbsig=seedin
      endif


c***********************************************************************
c 1.6 Read/evaluate Ls and the corresponding season number
c***********************************************************************

c 1.6.1 Solar longitude and local time
      if (datekey.eq.0) then  ! Earth date
c       (ie: "xdate" contains the Earth Julian date) 

        ! check that "xdate" is a date that is not less than jdate_min
        ! and no more than jdate_max
        if ((xdate.lt.jdate_min).or.(xdate.gt.jdate_max)) then
         if (output_messages) then
          write(out,*)'CALL_MCD error: The given Julian date ',
     &                'xdate=',xdate
          write(out,*)' lies outside of the [jdate_min:jdate_max] ',
     &                'range!!'
         endif
         ier=10
         goto 9999
        endif
        
        ! check that user did not try to impose a localtime
        if (localtime.ne.0.0) then
         if (output_messages) then
          write(out,*)'CALL_MCD error: If using Julian dates, ',
     &              'localtime=',localtime,' must be set to zero!'
         endif
         ier=12
         goto 9999
        endif
        
        ! compute values of ls,marsau,modelday corresponding to "xdate"
        call orbit(xdate,ls,marsau,modelday)
c       call mars_ltime with current longitude for correct local time:
c       localtime is a "local true solar time" or "local apparent time"   
c       which is such that the sun is always highest at noon
c      (differs from mean local time...) 
        call mars_ltime(lon,xdate,ls,loctime,lmeantime)

      else    ! Mars date (ie: "xdate" contains the value of Ls)
        ls=real(xdate)
        ! check that the value of ls makes sense
        if ((ls.lt.0.0).or.(ls.gt.360.0)) then
         if (output_messages) then
           write(out,*)'CALL_MCD error: The given value of ls, xdate=',
     &               xdate
           write(out,*)' should not lie outside the [0:360] range!!'
         endif
         ier=9
         goto 9999
        endif
        
        ! find the modelday which corresponds to ls
        modelday=ls2sol(ls)
        
        ! get localtime from input arguments
        loctime=localtime
        
        ! check that localtime makes sense
        if ((loctime.lt.0.0).or.(loctime.gt.24.0)) then
         if (output_messages) then
           write(out,*)'CALL_MCD error: local time, localtime=',
     &                 localtime
           write(out,*)' should not lie outside the [0:24] range!!'
         endif
         ier=11
         goto 9999
        endif
        
      endif ! of if (datekey.eq.0)

c 1.6.2 Compute the season numbers date corresponds to
      if (((dust.eq.4).or.(dust.eq.5).or.(dust.eq.6)).and.
     & (ls.lt.180.)) then
        if (output_messages) then
         write(out,*)'CALL_MCD Error: no dust storm scenario for Ls=',Ls
        endif
        ier=14
        goto 9999
      endif

      if (itimint.eq.0) then ! no seasonal interpolation
         call season(ls,nums)
      else if (itimint.ge.1) then
c        if seasonal interpolation
         call season2(ls,nums,nums2,wl,wh,dust)
         ! nums, nums2, wl and wh are now set
      endif   

c     Compute utime (true solar time at lon=0)
      call mars_ptime(lon,loctime,utime)

c***********************************************************************
c 2. load appropriate datasets from the database
c***********************************************************************

c     if the scenario changes, reset season numbers to reload all needed arrays
      if (dust.ne.dustprec) then
         numsprec=0
         nums2prec=0
         numsprec_sd=0
         numsprec_gw=0
         dustprec=dust
      endif

c     if the season number changes, mean value has to be reloaded
      if ((nums.ne.numsprec).or.(nums2.ne.nums2prec)) then
      ! Note: both nums and nums2 must be checked because for dust storm
      !       scenarios for limit cases (month 7 and 12) nums=nums2
      !       but moving on to month 8 (or 11), only nums2 (or nums)
      !       has changed
       if (output_messages) then
        write(out,*) 'Loading variables for new season'
       endif

c       open appropriate  NETCDF file(s)
c       *****************************              
        xdust = dust
        ! account for interpolation wrt previous MY, if needed
        if((nums.gt.nums2).AND.(24.le.dust).AND.(dust.le.32)) then 
         if((ls.lt.180).AND.(dust.gt.24)) xdust = dust-1   
        endif    

        call opend(unet,unetm,unetsd,unetm_up,unetsd_up,
     &              nums,xdust,dataset(1:lendataset),ierr)
        if (ierr.ne.0) then
           ier=ierr ! error code is set in opend
           goto 9999
        endif
        
        if (itimint.ge.1) then
         xdust = dust   
         ! account for interpolation wrt next MY, if needed
         if((nums.gt.nums2).AND.(24.le.dust).AND.(dust.le.32)) then
          if((ls.gt.180).AND.(dust.lt.32)) xdust = dust+1   
         endif
         call opend(unet,unetm2,unetsd2,unetm2_up,unetsd2_up,
     &               nums2,xdust,dataset2(1:lendataset),ierr)
        elseif (itimint.eq.0) then   
           call opend(unet,unetm2,unetsd2,unetm2_up,unetsd2_up,
     &     nums,dust,dataset2(1:lendataset),ierr)
        endif   
        if (ierr.ne.0) then
           ier=ierr ! error code is set in opend
           goto 9999
        endif

c       at the first call, load hybrid coordinates and orographic data
c       ************************************************************** 
        if (numsprec.eq.0) then ! true for the very first call to "CALL_MCD"
           ! load hybrid coordinates
           typevar='hybr'
           call loadvar_mcd(unetm,unetm_up,unetm2,unetm2_up,typevar,
     &                  aps,bps,sigma,pseudoalt,ierr)
                ! aps() and bps() (and pseudoalt()) are now set
           ! note that these are 'saved' and need only be read once
           if (ierr.ne.0) then
             ier=ierr ! error value is set in loadvar_mcd
             go to 9999
           endif
           
           ! load orography and areoid (on GCM grid)
           typevar='orog'
           call loadvar_mcd(unet,unet,unet,unet,typevar,
     &                  aps,bps,sigma,pseudoalt,ierr)
           ! taborog(), tabz0() and tabareo() (which are commons in constants_mcd.inc)
           ! are now set
           if (ierr.ne.0) then
              ier=ierr ! error value is set in loadvar_mcd
              go to 9999
           endif
        endif ! of if (numsprec.eq.0)

c       load mean variables      
c       *******************
        if (output_messages) then
          write(out,*) 'Loading new mean variable'
        endif
        typevar='mean'
        call loadvar_mcd(unetm,unetm_up,unetm2,unetm2_up,typevar,
     &                  aps,bps,sigma,pseudoalt,ierr)
        ! var_2d(), var_2d2(), var_3d() and var_3d2() are now set
        if (ierr.ne.0) then
           ier=ierr ! error value is set in loadvar_mcd
           go to 9999
        endif

        ! set numsprec and nums2prec to the current value of nums and nums2
        numsprec=nums
        nums2prec=nums2
        
        flag_ini_profi = 1 ! Force initializations in profi()

      end if !end if ((nums.ne.numsprec).or.(nums2.ne.nums2prec))

c     if large scale perturbations or extra variables are requested, 
c     load data if not yet done    
c     **************************************************************      
!      if ((perturkey.eq.2).or.(perturkey.eq.4).or.(extvarkey.eq.1)) then
      if ((perturkey.eq.2).or.(perturkey.eq.4)) then
         call loadeof(ueof,rdeof,dust,dataset,ierr)
         if (ierr.ne.0) then
            ier=ierr ! error value is set in loadeof
            go to 9999
         endif
      end if

c     if small scale perturbations or extra variables are requested, 
c     load data if not yet done
c     **************************************************************      
!      if (((perturkey.eq.3).or.(perturkey.eq.4).or.(extvarkey.eq.1))
      if (((perturkey.eq.3).or.(perturkey.eq.4))
     &     .and.(nums.ne.numsprec_gw)) then
        if (output_messages) then
          write(out,*) 'Loading new grwp variable'
        endif
        typevar='grwp'
        call loadvar_mcd(unet,unet,unet,unet,typevar,
     &                  aps,bps,sigma,pseudoalt,ierr)
        if (ierr.ne.0) then
           ier=ierr ! error value is set in loadvar_mcd
           go to 9999
        endif
        numsprec_gw=nums
        
        flag_ini_profi = 1 !  Force initializations in profi()

      end if
      
c     if n std dev perturbations or extra variables are requested, 
c     load data if not yet done
c     ************************************************************ 
!      if (((perturkey.eq.5).or.(extvarkey.eq.1))
      if ((extvarkeys(10).eq.1).or.(extvarkeys(18).eq.1).or.
     &    (extvarkeys(21).eq.1).or.(extvarkeys(22).eq.1).or.
     &    (extvarkeys(23).eq.1).or.(extvarkeys(24).eq.1).or.
     &    (extvarkeys(25).eq.1).or.(extvarkeys(27).eq.1).or.
     &    (extvarkeys(38).eq.1)) then
        rmsoutput=.true.
      else
        rmsoutput=.false.
      endif
      if (((perturkey.eq.5).or.(rmsoutput))
     &     .and.(nums.ne.numsprec_sd)) then
       ! Note: we load both kinds of RMS since (a rare but possible
       !       eventuality) 'zkey' could change from one call to call_mcd
       !       to the next
        ! Pressure-wise RMS
        typevar='rms'
        call loadvar_mcd(unetsd,unetsd_up,unetsd2,unetsd2_up,typevar,
     &                  aps,bps,sigma,pseudoalt,ierr)
        if (ierr.ne.0) then
           ier=ierr ! error value is set in loadvar_mcd
           go to 9999
        endif
        ! Altitude-wise RMS
        typevar='arms'
        call loadvar_mcd(unetsd,unetsd_up,unetsd2,unetsd2_up,typevar,
     &                  aps,bps,sigma,pseudoalt,ierr)
        if (ierr.ne.0) then
           ier=ierr ! error value is set in loadvar_mcd
           go to 9999
        endif
        numsprec_sd=nums
      
        flag_ini_profi = 1 !  Force initializations in profi()
      
      end if

c First "dummy" call of profi() from call_mcd().
c Does not produce the profile, does only initialization (if necessary).
      if(flag_ini_profi.eq.1) then
c      name='rho'
       call profi(rho_gcm,lon,lat,utime,name,ierr,itimint,wl,wh,1,0,0)
      endif

c***********************************************************************
c  3.  Retrieve main meteorological variables
c***********************************************************************

c  3.1 Get GCM orography for location lon,lat
      name='orography'
      call var2d(oro_gcm,lon,lat,1.0,name,ierr,itimint,wl,wh,1.0)
      ! oro_gcm is now set to the value of GCM orography at (lon,lat)

      name='areoid'
      call var2d(areoid_gcm,lon,lat,1.0,name,ierr,itimint,wl,wh,1.0)
      ! areoid_gcm is now set to the value of GCM areoid at (lon,lat)
      
      name='z0'
      call var2d(z_0,lon,lat,1.0,name,ierr,itimint,wl,wh,1.0)
      ! z_0 is now set to the value of GCM roughness length at (lon,lat)
      
      if (hireskey.eq.1) then
      ! get high resolution areoid at (lon,lat)
        call molareoid(dataset(1:lendataset),lon,lat,areoid_hr)
      endif
      
c  3.2 Get surface pressure at desired lon,lat and utime
      ! low resolution (ie: GCM) surface pressure
      name='ps'
      call var2d(ps_gcm,lon,lat,utime,name,ierr,itimint,wl,wh,1.0)

      ! high resolution (ie: using MOLA topography) surface pressure
      if (hireskey.eq.1) then
        call pres0(dataset(1:lendataset),dust,lat,lon,ls,utime,
     &             ps_gcm,oro_gcm,wl,wh,ps_hr,oro_hr,ierr)
        if (ierr.ne.0) then
          ier=ierr ! error value set in pres0
          goto 9999
        endif
        ! ps_hr and oro_hr are now set
        oroheight=oro_hr
        areoid=areoid_hr
        ps=ps_hr
        ps_psgcm=ps_hr/ps_gcm
      else
        oroheight=oro_gcm
        areoid=areoid_gcm
        ps=ps_gcm
        ps_psgcm=1.
      endif ! of if (hireskey.eq.1)

c  3.3 Build sigma levels
      ! GCM sigma levels
      do l=1,dimlevs
        sigma_gcm(l)=aps(l)/ps_gcm+bps(l)
      enddo

      ! high res sigma levels
      if (hireskey.eq.1) then
        call build_sigma_hr(sigma_gcm,ps_gcm,ps_hr,sigma_hr,p_pgcm)
      endif
      
      ! sigma levels that will be used further on
      if (hireskey.eq.1) then
        do l=1,dimlevs
          sigma(l)=sigma_hr(l)
        enddo
      else
        do l=1,dimlevs
          sigma(l)=sigma_gcm(l)
          p_pgcm(l)=1. ! No high_res, so set p_pgcm(:) to 1
        enddo
      endif
      
c  3.4 Get GCM temperature profile for location lon,lat and time utime
      name='temp'
      call profi(temp_gcm,lon,lat,utime,name,ierr,itimint,wl,wh,0,0,0)

c  3.5 Get GCM density profile for location lon,lat and time utime
      name='rho'
      call profi(rho_gcm,lon,lat,utime,name,ierr,itimint,wl,wh,0,0,0)

c  3.6 Compute R profile at location lon,lat and utime
      do i=1,dimlevs
        R_gcm(i)=(aps(i)+bps(i)*ps_gcm)/(rho_gcm(i)*temp_gcm(i))
      enddo

      if (hireskey.eq.1) then
        ! assign col_gcm() and col_hr() arrays
        call colint(aps,bps,oroheight,ps_gcm,ps_hr,R_gcm,temp_gcm,
     &              col_gcm,col_hr)
      endif

c  3.7 Find database levels encompassing sought altitude and determine
c      corresponding weights for vertical interpolation
      call getsi(xz,zkey,lon,lat,oroheight,areoid,ps,sigma,utime,dset,
     &           levhi,levlow,levweight,pratio,ierr,itimint,wl,wh,
     &           R_gcm,temp_gcm,zareoid,zradius,zsurface,zpressure,
     &           zmradius,sheight)
      ! levhi,levlow,levweight(),pratio are now set
      ! zareoid,zradius,zsurface,zpressure and sheight() are also set
      
      if (ierr.ne.0) then
! Note: Error message and error value are given in getsi
         ier=ierr
        goto 9999
      end if

      height=zareoid              ! height above areoid
      absheight=height-oroheight  ! height above surface
      
c      if underground : stop
      if (absheight.lt.-0.01) then
        if (output_messages) then
          write(out,*)'CALL_MCD Error: underground object '
          write(out,*)'                absheight=',absheight
        endif
        ier=17
        goto 9999
      endif

c  3.8 Compute mean field values (winds,density,temperature,pressure)
c      at chosen location

      name='u'
      call var3d(ul,lon,lat,zsurface,levhi,levlow,levweight,pratio,
     &      sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,itimint,wl,wh)
      name='v' 
      call var3d(vl,lon,lat,zsurface,levhi,levlow,levweight,pratio,
     &      sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,itimint,wl,wh)
      name='temp'
      call var3d(tl,lon,lat,zsurface,levhi,levlow,levweight,pratio,
     &      sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,itimint,wl,wh)
      name='rho'
      call var3d(rhol,lon,lat,zsurface,levhi,levlow,levweight,pratio,
     &      sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,itimint,wl,wh)

!     Mean pressure:
      tmeanl(1)=pratio*ps*(sigma(levlow)
     &     +(sigma(levhi)-sigma(levlow))*levweight(2))
      tmeanl(2)=rhol ! mean density
      tmeanl(3)=tl   ! mean temperature
      tmeanl(4)=ul   ! mean zonal wind
      tmeanl(5)=vl   ! mean meridional wind
      
      ! store these mean values for output
      do i=1,nmeanvar
         meanvar(i)=tmeanl(i)
      enddo

c***********************************************************************
c 4. Add perturbations, if required
c***********************************************************************

c 4.1 Add large scale variability, if required
c*********************************************      
      if ((perturkey.eq.2).or.(perturkey.eq.4)) then
      
      ! 4.1.0 Compute and add near-surface perturbations
        if (nearsurfnoise.and.(perturkey.gt.1)) then
          call nearsurfacenoise(ps,ps_noise_dev,ps_noise,
     &     temp_gcm_noise_dev,temp_gcm_noise,sheight)

          do l=1,dimlevs
            temp_gcm(l)=temp_gcm(l)+temp_gcm_noise(l)
          enddo
        else
          ps_noise=0
          temp_noise=0
        endif ! of if (nearsurfnoise.and.(perturkey.gt.1))

      ! 4.1.1 Build perturbed surface pressure
         name='ps'
         call eofpb(inicoordeof,scena,pertm,pertr,rdeof,lon,lat,
     &        levhi,levlow,levweight,pratio,modelday,name,ierr)
         ps_gcm_pert=ps_gcm+pertr+ps_noise
         ! recompute 'high resolution' surface pressure, if required
         if (hireskey.eq.1) then
           call pres0(dataset(1:lendataset),dust,lat,lon,ls,utime,
     &                ps_gcm_pert,oro_gcm,wl,wh,ps_hr_pert,oro_hr,ierr)
           if (ierr.ne.0) then
             ier=ierr ! error value set in pres0
             goto 9999
           endif
           ps=ps_hr_pert
           ps_psgcm=ps_hr_pert/ps_gcm_pert
         else
           ps=ps_gcm_pert
           ps_psgcm=1.
         endif ! if (hireskey.eq.1)

      ! 4.1.2 Build corresponding perturbed sigma levels
         ! perturbed GCM sigma levels
         do l=1,dimlevs
           sigma_gcm(l)=aps(l)/ps_gcm_pert+bps(l)
         enddo
         ! perturbed high res. sigma levels
         if (hireskey.eq.1) then
           call build_sigma_hr(sigma_gcm,ps_gcm_pert,
     &                         ps_hr_pert,sigma_hr,p_pgcm)
         endif         
         ! sigma levels which will be used further on
         if (hireskey.eq.1) then
           do l=1,dimlevs
             sigma(l)=sigma_hr(l)
           enddo
         else
           do l=1,dimlevs
             sigma(l)=sigma_gcm(l)
             p_pgcm(l)=1  ! no high res, p_pgcm(:)=1
           enddo
         endif
         
      ! 4.1.3 Build perturbed temperature profile
         name='temp'
         ! NB: imposed perturbed temperatures on the profile up to
         !     levhi+2 (since what happens above won't matter when
         !     integrating the hydrostatic equation in getsi)
         call profi_eof(inicoordeof,scena,temp_gcm,rdeof,lon,lat,
     &                  modelday,levhi+2,name,temp_gcm_pert,ierr)

         ! store 'unperturbed' values of levlow and levhi
!         oldlevlow=levlow
!         oldlevhi=levhi
!         oldlevweight2=levweight(2)

      ! 4.1.4 Compute new 'levhi,levlow,levweight,pratio'
         call getsi(xz,zkey,lon,lat,oroheight,areoid,ps,sigma,utime,
     &           dset,levhi,levlow,levweight,pratio,ierr,itimint,wl,wh,
     &          R_gcm,temp_gcm_pert,zareoid,zradius,zsurface,zpressure,
     &          zmradius,sheight)
      ! levhi,levlow,levweight(),pratio are now set
      ! zareoid,zradius,zsurface,zpressure and sheight are also set
         if (ierr.ne.0) then
! Note: Error message and error value are given in getsi
           ier=ierr
           goto 9999
         endif

         height=zareoid              ! height above areoid
         absheight=height-oroheight  ! height above surface

c      if underground : stop
         if (absheight.lt.-0.01) then
           if (output_messages) then
             write(out,*)'CALL_MCD Error: underground object '
             write(out,*)'                absheight=',absheight
           endif
           ier=17
           goto 9999
         endif

      ! 4.1.5 Get new mean values and add EOF perturbations
         name='u'
!         call var3d(ul,lon,lat,zsurface,levhi,levlow,levweight,pratio,
!     &              sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
!     &              itimint,wl,wh)
         if (zsurface.ge.z_0) then
         ! if above aerodynamic roughness length, compute and add perturbation
           call eofpb(inicoordeof,scena,pertm,pertr,rdeof,lon,lat,
     &              levhi,levlow,levweight,pratio,modelday,name,ierr)
           if (ierr.ne.0) then
             ier=ierr
             goto 9999
           endif
           ul=ul+pertr
c         else ! below aerodynamic roughness length, no perturbation
         endif

         name='v'
!         call var3d(vl,lon,lat,zsurface,levhi,levlow,levweight,pratio,
!     &              sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
!     &              itimint,wl,wh)
         if (zsurface.ge.z_0) then
         ! if above aerodynamic roughness length, compute and add perturbation
           call eofpb(inicoordeof,scena,pertm,pertr,rdeof,lon,lat,
     &              levhi,levlow,levweight,pratio,modelday,name,ierr)
           if (ierr.ne.0) then
             ier=ierr
             goto 9999
           endif
           vl=vl+pertr
c         else ! below aerodynamic roughness length, no perturbation
         endif
         
         name='temp'
!         call var3d(tl,lon,lat,zsurface,levhi,levlow,levweight,pratio,
!     &              sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
!     &              itimint,wl,wh)
         call eofpb(inicoordeof,scena,pertm,pertr,rdeof,lon,lat,
     &              levhi,levlow,levweight,pratio,modelday,name,ierr)
         if (ierr.ne.0) then
           ier=ierr
           goto 9999
         endif
         if (nearsurfnoise) then
           temp_noise=temp_gcm_noise(levlow)+
     &       (temp_gcm_noise(levhi)-temp_gcm_noise(levlow))*levweight(1)
         else
           temp_noise=0.0
         endif
         tl=tl+pertr+temp_noise

! straightforward (model 0) version:
!        get density as 'mean field' from perturbed state
!         name='rho'
!         call var3d(rhol,lon,lat,zsurface,levhi,levlow,levweight,pratio,
!     &              sheight,p_pgcm,ps_psgcm,utime,name,ierr,
!     &              itimint,wl,wh)
          ! deviation to mean density is thus:
!          pertrhoeof=rhol-meanvar(2)

! new version: (model 1) compute density from other variables
          ! first compute value of R using linear interpolation
!!          rhol=R_gcm(levlow)+(R_gcm(levhi)-R_gcm(levlow))*levweight(1)
          ! use unperturbed value of R=P/(rho*T)
          rhol=tmeanl(1)/(tmeanl(2)*tmeanl(3))
          ! density is then:
          rhol=zpressure/(rhol*tl)
          ! test: use lower value of R
!          rhol=zpressure/(R_gcm(levlow)*tl)
          ! test: use old value of pressure
!          rhol=tmeanl(1)/(rhol*tl)
!           rhol=tmeanl(1)/(rhol*tmeanl(3))
          ! deviation to mean density is thus:
!          pertrhoeof=rhol-meanvar(2)
! Note: computing R using linear interpolation yields 'arcs' in reconstructed
!       density deviation (wrt mean value) vertical profiles
!       Using the value of R corresponding to unperturbed P,T,rho yields
!       occasionnal 'kicks' (due to vertical interval 'mismatches' between
!       perturbed & unperturbed profiles) in density deviations.

! other new version: (model 2)
!   compute rho at levlow and levhi and then interpolate
!   but using 'old' mean level indexes & weights (otherwise there are
!   'kicks' in density deviation (wrt mean value) vertical
!   profiles when oldlevlow /= levlow and oldlevhi /= levhi).
          ! first compute pressures of encompassing levels:
!          lowpress=ps*sigma(oldlevlow)*pratio
!          hipress=ps*sigma(oldlevhi)*pratio
          ! now compute corresponding densities:
!          lowdens=lowpress/(R_gcm(oldlevlow)*temp_gcm_pert(oldlevlow))
!          hidens=hipress/(R_gcm(oldlevhi)*temp_gcm_pert(oldlevhi))
          ! interpolate density using lowdens and hidens
!          rhol=lowdens*p_pgcm(oldlevlow)+oldlevweight2*
!     &         (hidens*p_pgcm(oldlevhi)-lowdens*p_pgcm(oldlevlow))
!          rhol=rhol*pratio ! correction for densities out of sigma range

          ! deviation to mean density is thus:
!          pertrhoeof=rhol-meanvar(2)

!      else if(extvarkey.eq.1) then
      ! just compute for density as an extra variable
!        name='rho'
!        call eofpb(scena,pertm,pertrhoeof,rdeof,lon,lat,
!     &             levhi,levlow,levweight,pratio,modelday,name,ierr)
!        pertrhoeof=0.0
      endif ! if ((perturkey.eq.2).or.(perturkey.eq.4)) 

c 4.2 Add small scale (gravity wave) variability, if required
c************************************************************      
      if ((perturkey.eq.3).or.(perturkey.eq.4)) then
         name='u'
         absheight=height-oroheight
         call grwpb(dset,pertu,rangw1,rangw2,dev,lamda,lon,lat,
     &              absheight,tmeanl(2),tmeanl(4),tmeanl(5),ps,sigma,
     &              utime,name,ierr,itimint,wl,wh,
     &              R_gcm,temp_gcm,sheight,p_pgcm,ps_psgcm)
         name='v'
         call grwpb(dset,pertv,rangw1,rangw2,dev,lamda,lon,lat,
     &              absheight,tmeanl(2),tmeanl(4),tmeanl(5),ps,sigma,
     &              utime,name,ierr,itimint,wl,wh,
     &              R_gcm,temp_gcm,sheight,p_pgcm,ps_psgcm)
         name='temp'
         call grwpb(dset,pertt,rangw1,rangw2,dev,lamda,lon,lat,
     &              absheight,tmeanl(2),tmeanl(4),tmeanl(5),ps,sigma,
     &              utime,name,ierr,itimint,wl,wh,
     &              R_gcm,temp_gcm,sheight,p_pgcm,ps_psgcm)
         name='rho'
         call grwpb(dset,pertrho,rangw1,rangw2,dev,lamda,lon,lat,
     &              absheight,tmeanl(2),tmeanl(4),tmeanl(5),ps,sigma,
     &              utime,name,ierr,itimint,wl,wh,
     &              R_gcm,temp_gcm,sheight,p_pgcm,ps_psgcm)
         pertrhogw = pertrho
         ul = ul + pertu
         vl = vl + pertv
         tl = tl + pertt
         ! Check that the resulting density is not unphysical
         ! (was reported to happen in rare cases), and thus limit density
         ! to never be less than 5% of what it was before adding pert.
         ! (5% is quite extreme, tests show large decreases most often
         ! do not lead to less than 50% of value)
         rhol=max(rhol+pertrho,0.05*rhol)
      else if(extvarkeys(28).eq.1) then
      ! just compute for density to output as an extra variable
        name='rho'
        if (gwlength.eq.0.0) then ! use a default wavelength
          lamda=16.E3
        else
          lamda=gwlength
          ! Note: it was checked that 2000<gwlength<30000 in 1.4
        endif
        call grwpb(dset,pertrhogw,rangw1,rangw2,dev,lamda,lon,lat,
     &             absheight,tmeanl(2),tmeanl(4),tmeanl(5),ps,sigma,
     &             utime,name,ierr,itimint,wl,wh,
     &             R_gcm,temp_gcm,sheight,p_pgcm,ps_psgcm)
      endif

c 4.3 Add n sigmas, if required
c****************************** 
      if (perturkey.eq.5) then
        if (zkey.eq.4) then ! Get pressure-wise RMSs
          name='rmsu'
          call var3d(pertu,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
          name='rmsv'
          call var3d(pertv,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
          name='rmsrho'
          call var3d(pertrho,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
          name='rmstemp'
          call var3d(pertt,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
        else ! get altitude-wise RMSs
          name='armsu'
          call var3d(pertu,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
          name='armsv'
          call var3d(pertv,lon,lat,zsurface,levhi,levlow,levweight,
     &         pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
          name='armsrho'
          call var3d(pertrho,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
          name='armstemp'
          call var3d(pertt,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &               itimint,wl,wh)
        endif ! of if (zkey.eq.4)
         name='rmsps'
         call var2d(pertps,lon,lat,1.0,name,ierr,itimint,wl,wh,ps_psgcm)
         ul = ul + (nbsig*pertu)
         vl = vl + (nbsig*pertv)
         rhol = rhol + (nbsig*pertrho)
         tl = tl + (nbsig*pertt)
         ps = ps + (nbsig*pertps)
         ! check that this has not led to unphysical values
         if (rhol.lt.0.) then
           if (output_messages) then
             write(out,*)'CALL_MCD Error: unphysical density:',rhol
             write(out,*)'  due to addition of ',nbsig,
     &                 ' times the day to day variability'
           endif
           ier=18
           goto 9999
         endif
         if (tl.lt.0.) then
           if (output_messages) then
             write(out,*)'CALL_MCD Error: unphysical temperature:',tl
             write(out,*)'  due to addition of ',nbsig,
     &                 ' times the day to day variability'
           endif
           ier=19
           goto 9999
         endif
         if (ps.lt.0.) then
           if (output_messages) then
             write(out,*)'CALL_MCD Error: unphysical pressure:',ps
             write(out,*)'  due to addition of ',nbsig,
     &                 ' times the day to day variability'
           endif
           ier=20
           goto 9999
         endif
      endif ! of if (perturkey.eq.5)

c 4.4 Store atmospheric fields (pressure, temperature, density, winds)
c     (re-)compute atmospheric pressure
      pl=ps*(sigma(levlow)+(sigma(levhi)-sigma(levlow))*levweight(2))
      pl=pratio*pl

      t=tl      ! temperature
      p=pl      ! pressure
      rho=rhol  ! density
      u=ul      ! zonal wind
      v=vl      ! meridional wind

c 4.5 (re)compute extra variables, if required
c*********************************************   
! Preliminary stuff:
         do i=1,nextvar 
            extvar(i)=0.0
         enddo

! call air_properties under many conditions:
       air_prop_flag=.false.
       if ((extvarkeys(49).eq.1).or.(extvarkeys(50).eq.1).or.
     &     (extvarkeys(8).eq.1).or.(extvarkeys(9).eq.1)) then
         ! prerequisites to air_properties:
         name='vmr_co2'
         call var3d(vmr_co2,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
         name='vmr_n2'
         call var3d(vmr_n2,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
         name='vmr_ar'
         call var3d(vmr_ar,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
         name='vmr_o'
         call var3d(vmr_o,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
         name='vmr_co'
         call var3d(vmr_co,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
         call air_properties(t,vmr_co2,vmr_n2,vmr_ar,vmr_o,vmr_co,
     &                       Cp,gamma,viscosity,Rgas)
         air_prop_flag=.true.
       endif

         ! the first few extvar() are always computed
         extvar(1)=zradius
         extvar(2)=zareoid
         extvar(3)=zsurface
         extvar(4)=oroheight
         extvar(5)=ls
         extvar(6)=loctime
         extvar(7)=utime

         if (extvarkeys(8).eq.1) extvar(8)=Cp
         if (extvarkeys(9).eq.1) extvar(9)=gamma
         
         ! Density RMS day to day variations (kg/m3)
         if (extvarkeys(10).eq.1) then
           if (zkey.eq.4) then   ! vertical coordinate is pressure
             extvar(10)=0
           else
             name='armsrho'
             call var3d(rmsrho,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
             extvar(10)=rmsrho
           endif
         endif
         
         ! Local mean solar time (hrs) at sought longitude
         if (extvarkeys(11).eq.1) then
           ! only computed if using earth date as input
           if (datekey.eq.0) extvar(11)=lmeantime
         endif
         
         ! Sun-Mars distance:
         if (extvarkeys(12).eq.1) then
           if (datekey.eq.0) then ! it has already been computed
             extvar(12)=marsau
           else
             call sunmarsdistance(ls,marsau)
             extvar(12)=marsau
           endif
         endif
         
         ! Atmospheric scale height
         if (extvarkeys(13).eq.1) then
           Rnew=p/(rho*t)
           pscaleheight=Rnew*t/(g0*a0**2/(a0+real(zareoid))**2)
           extvar(13)=pscaleheight
         endif
         
         ! GMC orography (always computed)
         if (extvarkeys(14).eq.1) extvar(14)=oro_gcm 
         
         ! Surface temperature (K)
         if (extvarkeys(15).eq.1) then
           name='tsurf'
           call var2d(tsurf,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(15)=tsurf
         endif
         
         if (extvarkeys(16).eq.1) then
           call max2d(tsurfmax,lon,lat,name,ierr,itimint,wl,wh)
           extvar(16)=tsurfmax
         endif
         
         if (extvarkeys(17).eq.1) then
           call min2d(tsurfmin,lon,lat,name,ierr,itimint,wl,wh)
           extvar(17)=tsurfmin
         endif
         
         if (extvarkeys(18).eq.1) then
           name='rmstsurf'
           call var2d(rmstsurf,lon,lat,1.0,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(18)=rmstsurf
         endif
         
         if (extvarkeys(19).eq.1) then
           if (hireskey.eq.1) then
             extvar(19)=ps_hr     ! high res Ps
           else
             extvar(19)=ps_gcm
           endif
         endif
         
         extvar(20)=ps_gcm ! ps_gcm is always computed
         
         ! Atmospheric pressure RMS day to day variations
         if (extvarkeys(21).eq.1) then
           if (zkey.eq.4) then   ! vertical coordinate is pressure
             extvar(21)=0.0
           else
             name='armspressure'
             call var3d(altrmsp,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
             extvar(21)=altrmsp
           endif
         endif
         
         ! Surface pressure RMS day to day variation
         if (extvarkeys(22).eq.1) then
           name='rmsps'
           call var2d(rmsps,lon,lat,1.0,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(22)=rmsps
         endif
         
         ! Atmospheric temperature RMS day to day variations (K)
         if (extvarkeys(23).eq.1) then
           if (zkey.eq.4) then ! pressure-wise RMS
             name='rmstemp'
             call var3d(rmst,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           else ! altitude-wise RMS
             name='armstemp'
             call var3d(rmst,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           endif
           extvar(23)=rmst
         endif
         
         ! Zonal wind RMS day to day variations (m/s)
         if (extvarkeys(24).eq.1) then
           if (zkey.eq.4) then ! pressure-wise RMS
             name='rmsu'
             call var3d(rmsu,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           else ! altitude-wise RMS
             name='armsu'
             call var3d(rmsu,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           endif
           extvar(24)=rmsu
         endif
         
         ! Meridional wind RMS day to day variations (m/s)
         if (extvarkeys(25).eq.1) then
           if (zkey.eq.4) then ! pressure-wise RMS
             name='rmsv'
             call var3d(rmsv,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           else ! altitude-wise RMS
             name='armsv'
             call var3d(rmsv,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           endif
           extvar(25)=rmsv
         endif
         
         ! Vertical wind component (m/s) >0 when downwards!
         if (extvarkeys(26).eq.1) then
           name='w' 
           call var3d(w_l,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(26)=w_l
         endif
         
         ! Vertical wind RMS day to day variations (m/s)
         if (extvarkeys(27).eq.1) then
           if (zkey.eq.4) then ! pressure-wise RMS
             name='rmsw'
             call var3d(rmsw,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           else ! altitude-wise RMS
             name='armsw'
             call var3d(rmsw,lon,lat,zsurface,levhi,levlow,levweight,
     &               pratio,sheight,p_pgcm,sigma,ps_psgcm,1.0,name,ierr,
     &                itimint,wl,wh)
           endif
           extvar(27)=rmsw
         endif
         
         ! Gravity wave density perturbation
         if (extvarkeys(28).eq.1) extvar(28)=pertrhogw

         ! Surface roughness length (m) (always computed)
         if (extvarkeys(29).eq.1) extvar(29)=z_0
         
         ! Solar flux reflected to space (W/m2)
         if (extvarkeys(30).eq.1) then
           name='fluxtop_sw'
           call var2d(fluxtop_sw,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(30)=fluxtop_sw
         endif

         if (extvarkeys(31).eq.1) then
           name='fluxsurf_lw'
           call var2d(fluxsurf_lw,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(31)=fluxsurf_lw
         endif
         
         if (extvarkeys(32).eq.1) then
           name='fluxsurf_sw'
           call var2d(fluxsurf_sw,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(32)=fluxsurf_sw
         endif
         
         if (extvarkeys(33).eq.1) then
           name='fluxtop_lw'
           call var2d(fluxtop_lw,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(33)=fluxtop_lw
         endif
         
         ! Surface H2O (seasonal frost),
         ! set to 0.5 kg/m2 for perrenial H2O ice deposits
         if (extvarkeys(34).eq.1) then
           name='h2oice'
           call var2d(surf_h2o_ice,lon,lat,utime,name,ierr,itimint,wl,
     &                wh,ps_psgcm)
           extvar(34)=surf_h2o_ice
         endif
         
         ! Surface CO2 ice layer (kg/m2)
         if (extvarkeys(35).eq.1) then
           name='co2ice'
           call var2d(co2ice,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(35)=co2ice
         endif
         
         if (extvarkeys(36).eq.1) then
           ! DOD: Dust column visible optical depth
           name='tauref'
           call var2d(dod,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! NB: tauref is dod at 610 Pa, so rescale it to current ps
           dod=dod*(ps/610.)
           extvar(36)=dod
         endif
         
         ! DOD RMS day to day variations 
         if (extvarkeys(37).eq.1) then
           name='rmstauref'
           call var2d(rmsdod,lon,lat,1.0,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(37)=rmsdod
         endif
         
         ! dust mass mixing ratio (kg/kg_air)
         if ((extvarkeys(38).eq.1)) then
!             call dustmix(scena,lat,Ls,dod,ps,p,dust_mmr) ! dustmix for MCD4
           name='dustq'
           call var3d(dust_mmr,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(38)=dust_mmr
         endif
         
         ! dust effective radius (m)
         if ((extvarkeys(39).eq.1)) then
           name="reffdust"
           call var3d(dust_reff,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(39)=dust_reff
         endif
         
         ! Dust deposition on flat surface (kg / m2 / s)
         if (extvarkeys(40).eq.1) then
           ! get the dust effective radius profile (we need profile_reff(1))
           name="reffdust"
           call profi(profile_reff,lon,lat,utime,name,ier,itimint,
     &                  wl,wh,0,1,1)
           ! get the dust mass mixing ratio profile (we need profile_mmr(1))
           name='dustq'
           call profi(profile_mmr,lon,lat,utime,name,ier,itimint,
     &                  wl,wh,0,1,1)
           ! compute dust deposition
           call dust_deposition(ps,temp_gcm(1),rho_gcm(1),
     &                       profile_mmr(1),profile_reff(1),dust_dep)
           extvar(40)=dust_dep
         endif
         
         ! Water vapour column (kg/m2)
         if (extvarkeys(41).eq.1) then
           name='col_h2ovapor'
           call var2d(col_h2ovapor,lon,lat,utime,name,ierr,itimint,wl,
     &                wh,ps_psgcm)
           extvar(41)=col_h2ovapor
           
          if(hireskey.eq.1) then 
           name='vmr_h2ovapor'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(41)=extvar(41)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(41).eq.1)
         
         ! Water vapor vol. mixing ratio (mol/mol)
         if (extvarkeys(42).eq.1) then
           name='vmr_h2ovapor'
           call var3d(vmr_h2o,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(42)=vmr_h2o
         endif
         
         ! Water ice column (kg/m2)
         if (extvarkeys(43).eq.1) then
           name='col_h2oice'
           call var2d(col_h2oice,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(43)=col_h2oice
          
          if(hireskey.eq.1) then 
           name='vmr_h2oice'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,      
     &                  col_gcm,col_hr) 
           extvar(43)=extvar(43)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(43).eq.1)
         
         ! Water ice mixing ratio (mol/mol)
         if (extvarkeys(44).eq.1) then
           name='vmr_h2oice'
           call var3d(vmr_h2oice,lon,lat,zsurface,levhi,levlow,
     &     levweight,pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,
     &              ierr,itimint,wl,wh)
           extvar(44)=vmr_h2oice
         endif
         
         ! Water ice effective radius (m)
         if (extvarkeys(45).eq.1) then
           name='reffice'
           call var3d(h2oice_reff,lon,lat,zsurface,levhi,levlow,
     &     levweight,pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,
     &              ierr,itimint,wl,wh)
           ! Corrections to remove extreme unphysical values (due to
           ! undelying averaging and limit with low ice content/high radii)
           if (extvarkeys(44).eq.0) then ! in case vmr_h2oice is not known
             ! get water ice mixing ratio
             name='vmr_h2oice'
             call var3d(vmr_h2oice,lon,lat,zsurface,levhi,levlow,
     &       levweight,pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,
     &              ierr,itimint,wl,wh)
           endif
           if (vmr_h2oice.le.1.e-6) then ! very low ice content
             h2oice_reff=0. ! set radius to zero
           endif
           if (h2oice_reff.ge.1.e-4) then ! unphysically high reff
             h2oice_reff=1.e-4 ! allow it to be at most 100 microns
           endif
           extvar(45)=h2oice_reff
         endif
         
         ! PBL height (m)
         if ((extvarkeys(46).eq.1).or.(extvarkeys(49).eq.1)
     &                            .or.(extvarkeys(50).eq.1)) then
           name='zmax' 
           call var2d(zmax,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           extvar(46)=zmax
         endif

         ! Max. upward convective wind within the PBL (m/s)
         if (extvarkeys(47).eq.1) then
           name='wstar'
           call var2d(wstar,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! wstar is the convective vertical velocity scale 
           ! multiply it by 2.75 to obtain maximum updraft velocity
           extvar(47)=wstar*2.75
         endif
         
         ! Max. downward convective wind within the PBL (m/s)
         if (extvarkeys(48).eq.1) then
           name='wstar'
           call var2d(wstar,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! wstar is the convective vertical velocity scale 
           ! multiply it by 1.75 to obtain maximum downdraft velocity
           extvar(48)=wstar*1.75
         endif

         ! Vertical velocity variance at z (m2/s2)
         if (extvarkeys(49).eq.1) then
           ! This has only a meaning inside the PBL; set it to zero otherwise
           if (zsurface.le.zmax) then
             name='vvv'
             call var3d(vvv,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           else
             vvv=0.
           endif
           extvar(49)=vvv
         endif 
         
         ! Vertical eddy heat flux at z (m/s/K)
         if (extvarkeys(50).eq.1) then
           ! This has only a meaning inside the PBL; set it to zero otherwise
           if (zsurface.le.zmax) then
             name='vhf'
             call var3d(vhf,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           else
             vhf=0.
           endif
           extvar(50)=vhf
         endif

         ! Surface stress (Kg/m/s2)
         if (extvarkeys(51).eq.1) then
           name='surfstress' ! is a 2D var but pbl_parameter is handled in var3d
           ! note that high zsurface values lead to numerical overflows in
           ! pbl_parameters, so we call var3d with zsurface=10.
           call var3d(surfstress,lon,lat,10.,levhi,levlow,
     &       levweight,pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,
     &         ierr,itimint,wl,wh)
           extvar(51)=surfstress
         endif
         
         ! Sensible heat flux (W/m2)
         !     (negative when flux is from surface to atmosphere)
         if (extvarkeys(52).eq.1) then
           name='sensib_flux' ! is a 2D var but pbl_parameter is handled in var3d
           ! note that high zsurface values lead to numerical overflows in
           ! pbl_parameters, so we call var3d with zsurface=10.
           call var3d(sensib_flux,lon,lat,10.,levhi,levlow,
     &       levweight,pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,
     &         ierr,itimint,wl,wh)
           extvar(52)=sensib_flux
         endif
         
         if (extvarkeys(53).eq.1) extvar(53)=Rgas
         if (extvarkeys(54).eq.1) extvar(54)=viscosity
         
         ! Solar zenith angle (deg)
         if (extvarkeys(56).eq.1) then
           call solarzenithangle(lat,ls,loctime,solzenang)
           extvar(56)=solzenang
         endif
         
         ! CO2 vol. mixing ratio (mol/mol)
         if (extvarkeys(57).eq.1) then
           if (.not.air_prop_flag) then
             name='vmr_co2'
             call var3d(vmr_co2,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           endif
           extvar(57)=vmr_co2
         endif
         
         ! N2 vol. mixing ratio (mol/mol)
         if (extvarkeys(58).eq.1) then
           if (.not.air_prop_flag) then
             name='vmr_n2'
             call var3d(vmr_n2,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           endif
           extvar(58)=vmr_n2
         endif
         
         !Ar vol. mixing ratio (mol/mol)
         if (extvarkeys(59).eq.1) then
           if (.not.air_prop_flag) then
             name='vmr_ar'
             call var3d(vmr_ar,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           endif
           extvar(59)=vmr_ar
         endif
         
         ! CO vol. mixing ratio (mol/mol)
         if (extvarkeys(60).eq.1) then
           if (.not.air_prop_flag) then
             name='vmr_co'
             call var3d(vmr_co,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           endif
           extvar(60)=vmr_co
         endif
         
         ! O vol. mixing ratio (mol/mol)
         if (extvarkeys(61).eq.1) then
           if (.not.air_prop_flag) then
             name='vmr_o'
             call var3d(vmr_o,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           endif
           extvar(61)=vmr_o
         endif
         
         !O2 vol. mixing ratio (mol/mol)
         if (extvarkeys(62).eq.1) then
           name='vmr_o2'
           call var3d(vmr_o2,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(62)=vmr_o2
         endif
                  
         ! O3 ozone vol. mixing ratio (mol/mol)
         if (extvarkeys(63).eq.1) then
           name='vmr_o3'
           call var3d(vmr_o3,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(63)=vmr_o3
         endif
         
         !H vol. mixing ratio (mol/mol)
         if (extvarkeys(64).eq.1) then
           name='vmr_h'
           call var3d(vmr_h,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(64)=vmr_h
         endif
         
         !H2 vol. mixing ratio (mol/mol)
         if (extvarkeys(65).eq.1) then
           name='vmr_h2'
           call var3d(vmr_h2,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(65)=vmr_h2
         endif
         
         ! electron number density (particules/cm3)
         if (extvarkeys(66).eq.1) then
           ! Give an electron Volume Mixing Ratio only below 5.e-6 Pa
           ! (above the "chemistry ionosphere", one would needs to model 
           !  the full ionosphere dynamics, which isn't the case here)
           if (p.ge.5.e-6) then
            name='vmr_elec'
            call var3d(vmr_elec,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
            ! vmr_elec is in mol/mol , convert to electronic number
            ! density in particules/cm3
            !c_e=vmr_e*P_air[Pa]*6.023e23[part/mol]/(8.314*1e6[Pa*cm3/(mol*K)]*T[K])->part/cm3
            vmr_elec=vmr_elec*6.022e23*pl/(8.314e6*tl)
           else
             vmr_elec=0.
           endif
           extvar(66)=vmr_elec
         endif
         
         ! CO2 column (kg/m2)
         if (extvarkeys(67).eq.1) then
           name="c_co2"
           call var2d(col_co2,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_co2 is given in mol cm-2; convert to kg/m2
           col_co2=col_co2*(44./6.022e22)
           extvar(67)=col_co2

          if(hireskey.eq.1) then 
           name='vmr_co2'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(67)=extvar(67)*col_hr2gcm
          endif

         endif ! of if (extvarkeys(67).eq.1)
         
         ! N2 column (kg/m2)
         if (extvarkeys(68).eq.1) then
           name="c_n2"
           call var2d(col_n2,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_n2 is given in mol cm-2; convert to kg/m2
           col_n2=col_n2*(28./6.022e22)
           extvar(68)=col_n2
         
          if(hireskey.eq.1) then 
           name='vmr_n2'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(68)=extvar(68)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(68).eq.1)
         
         ! Ar column (kg/m2)
         if (extvarkeys(69).eq.1) then
           name="c_ar"
           call var2d(col_ar,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_n2 is given in mol cm-2; convert to kg/m2
           col_ar=col_ar*(40./6.022e22)
           extvar(69)=col_ar
           
          if(hireskey.eq.1) then 
           name='vmr_ar'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(69)=extvar(69)*col_hr2gcm
          endif
           
         endif ! of if (extvarkeys(69).eq.1)
         
         ! CO column (kg/m2)
         if (extvarkeys(70).eq.1) then
           name="c_co"
           call var2d(col_co,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_co is given in mol cm-2; convert to kg/m2
           col_co=col_co*(28./6.022e22)
           extvar(70)=col_co
           
          if(hireskey.eq.1) then 
           name='vmr_co'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(70)=extvar(70)*col_hr2gcm
          endif
           
         endif ! of if (extvarkeys(70).eq.1)
         
         ! O column (kg/m2)
         if (extvarkeys(71).eq.1) then
           name="c_o"
           call var2d(col_o,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_o is given in mol cm-2; convert to kg/m2
           col_o=col_o*(16./6.022e22)
           extvar(71)=col_o
         
          if(hireskey.eq.1) then 
           name='vmr_o'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(71)=extvar(71)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(71).eq.1)
         
         ! O2 column (kg/m2)
         if (extvarkeys(72).eq.1) then
           name="c_o2"
           call var2d(col_o2,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_o2 is given in mol cm-2; convert to kg/m2
           col_o2=col_o2*(32./6.022e22)
           extvar(72)=col_o2
         
          if(hireskey.eq.1) then 
           name='vmr_o2'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,      
     &                  col_gcm,col_hr) 
           extvar(72)=extvar(72)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(72).eq.1)
         
         ! O3 column (kg/m2)
         if (extvarkeys(73).eq.1) then
           name="c_o3"
           call var2d(col_o3,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_o3 is given in mol cm-2; convert to kg/m2
           col_o3=col_o3*(48./6.022e22)
           extvar(73)=col_o3
           
          if(hireskey.eq.1) then 
           name='vmr_o3'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(73)=extvar(73)*col_hr2gcm
          endif
           
         endif ! of if (extvarkeys(73).eq.1)
         
         ! H column (kg/m2)
         if (extvarkeys(74).eq.1) then
           name="c_h"
           call var2d(col_h,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_h is given in mol cm-2; convert to kg/m2
           col_h=col_h*(1./6.022e22)
           extvar(74)=col_h
         
          if(hireskey.eq.1) then 
           name='vmr_h'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                  col_gcm,col_hr) 
           extvar(74)=extvar(74)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(74).eq.1)
         
         ! H2 column (kg/m2)
         if (extvarkeys(75).eq.1) then
           name="c_h2"
           call var2d(col_h2,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_h2 is given in mol cm-2; convert to kg/m2
           col_h2=col_h2*(2./6.022e22)
           extvar(75)=col_h2
         
          if(hireskey.eq.1) then 
           name='vmr_h2'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,      
     &                  col_gcm,col_hr) 
           extvar(75)=extvar(75)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(75).eq.1)
         
         ! total electonic content (TEC) (particules/m2)
         if (extvarkeys(76).eq.1) then
           name="c_elec"
           call var2d(col_elec,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_elec is given in mol cm-2; convert to particules/m2
           col_elec=col_elec*((1./1822.89)/(6.022e22*9.109e-31))
           extvar(76)=col_elec
         endif

         ! He column (kg/m2)
         if (extvarkeys(77).eq.1) then
           name="c_he"
           call var2d(col_he,lon,lat,utime,name,ierr,itimint,wl,wh,
     &                ps_psgcm)
           ! c_he is given in mol cm-2; convert to kg/m2
           col_he=col_he*(4./6.022e22)
           extvar(77)=col_he
         
          if(hireskey.eq.1) then 
           name='vmr_he'           
           call colcor(col_hr2gcm,lon,lat,utime,name,ierr,itimint,wl,wh,      
     &                  col_gcm,col_hr) 
           extvar(77)=extvar(77)*col_hr2gcm
          endif
         
         endif ! of if (extvarkeys(77).eq.1)

         !He vol. mixing ratio (mol/mol)
         if (extvarkeys(78).eq.1) then
           name='vmr_he'
           call var3d(vmr_he,lon,lat,zsurface,levhi,levlow,levweight,
     &             pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,name,ierr,
     &              itimint,wl,wh)
           extvar(78)=vmr_he
         endif

c     copy to output
      pres=p
      dens=rho
      temp=t
      zonwind=u
      merwind=v

      return

c     Error handling : all the outputs are set to a missing data value
 9999 pres=-999.0
      dens=-999.0
      temp=-999.0
      zonwind=-999.0
      merwind=-999.0
      do i=1,nmeanvar
         meanvar(i)=-999.0
      end do
      do i=1,nextvar
         extvar(i)=-999.0
      enddo

      end ! end of call_mcd

!=======================================================================
      SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0, 
     & zzlay,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out,n_out,
     & T_out,u_out,ustar,tstar,vhf,vvv)
      IMPLICIT NONE
!=======================================================================
!
!   Anlysis of the PBL from input temperature, wind field and thermals outputs.
!
!   -------  
!
!   Author: Arnaud Colaitis 09/01/12
!   MCD VERSION from gcm svn revision 662
!   -------
!
!   Arguments:
!   ----------
!
!   inputs:
!   ------
!     ngrid            size of the horizontal grid
!     nlay             size of the vertical grid
!     pz0(ngrid)       surface roughness length
!     zzlay(ngrid,nlay)   height of mid-layers
!     pu(ngrid,nlay)   u component of the wind
!     pv(ngrid,nlay)   v component of the wind
!     wstar_in(ngrid)  free convection velocity in PBL
!     hfmax(ngrid)     maximum vertical turbulent heat flux in thermals
!     zmax(ngrid)      height reached by the thermals (pbl height)
!     pts(ngrid)       surface temperature
!     ph(ngrid,nlay)   potential temperature T*(p/ps)^kappa
!     z_out(n_out)     heights of interpolation
!     n_out            number of points for interpolation
!
!   outputs:
!   ------
!
!     Teta_out(ngrid,n_out)  interpolated teta
!     u_out(ngrid,n_out)     interpolated u
!     ustar(ngrid)     friction velocity
!     tstar(ngrid)     friction temperature
!
!
!=======================================================================
!
!-----------------------------------------------------------------------
!   Declarations:
!   -------------

!   Arguments:
!   ----------

      INTEGER, INTENT(IN) :: ngrid,nlay,n_out
      REAL, INTENT(IN) :: pz0(ngrid),ps(ngrid),pplay(ngrid,nlay)
      REAL, INTENT(IN) :: zzlay(ngrid,nlay)
      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
      REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid)
      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
      REAL, INTENT(IN) :: z_out(n_out)

!    Outputs:
!    --------

      REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out)
      REAL Teta_out(ngrid,n_out)
      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid)

!   Local:
!   ------

      INTEGER ig,k,n
      REAL karman,nu,rcp,pg
      DATA karman,nu,rcp,pg/.41,0.001,0.2567930,3.72/
      SAVE karman,nu,rcp,pg

!    Local(2):
!    ---------

      REAL zout
      REAL rib(ngrid)  ! Bulk Richardson number
      REAL fm(ngrid) ! stability function for momentum
      REAL fh(ngrid) ! stability function for heat
      REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
          ! phim = 1+betam*zeta   or   (1-bm*zeta)**am
          ! phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
      REAL betam, betah, alphah, bm, bh, lambda
          ! ah and am are assumed to be -0.25 and -0.5 respectively
      REAL cdn(ngrid),chn(ngrid)  ! neutral momentum and heat drag coefficient
      REAL pz0t        ! initial thermal roughness length. (local)
      REAL ric         ! critical richardson number
      REAL reynolds(ngrid)    ! reynolds number for UBL
      REAL prandtl(ngrid)     ! prandtl number for UBL
      REAL pz0tcomp(ngrid)     ! computed z0t
      REAL ite
      REAL residual!,zcd0,z1
      REAL pcdv(ngrid),pcdh(ngrid)
      REAL zu2(ngrid)                  ! Large-scale wind at first layer
      REAL pbl_teta(ngrid)             ! mixed-layer potential temperature
      INTEGER pbl_height_index(ngrid)  ! index of nearest vertical grid point for zmax
      REAL dteta(ngrid,nlay),x(ngrid)  ! potential temperature gradient and z/zi
      REAL dvhf(ngrid),dvvv(ngrid)     ! dimensionless vertical heat flux and 
                                       ! dimensionless vertical velocity variance
      REAL vhf(ngrid),vvv(ngrid)       ! vertical heat flux and vertical velocity variance
      INTEGER ii(1)
      LOGICAL dummy

!------------------------------------------------------------------------
!------------------------------------------------------------------------
! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
!------------------------------------------------------------------------
!------------------------------------------------------------------------
      zout=0 ! dummy initialization to get rid of compiler warning
      dummy=.false.
      IF (z_out(1) .le. 0.) THEN
      dummy=.true.
      ENDIF

      DO n=1,n_out

c Initialisation :

      ustar(:)=0.
      tstar(:)=0.
      IF (dummy) THEN
      zout=1.
      ELSE
      zout=z_out(n)
      ENDIF
      reynolds(:)=0.
      pz0t = 0.
      pz0tcomp(:) = 0.1*pz0(:)
      rib(:)=0.
      pcdv(:)=0.
      pcdh(:)=0.

      ! this formulation assumes alphah=1., implying betah=betam
      ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :

      bm=16.            !UBL
      bh=16.            !UBL
      alphah=1.
      betam=5.         !SBL
      betah=5.         !SBL
      lambda=(sqrt(bh/bm))/alphah
      ric=betah/(betam**2)
      DO ig=1,ngrid
       ite=0.
       residual=abs(pz0tcomp(ig)-pz0t)

       zu2(ig)=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
     &     + (log(1.+0.7*wstar_in(ig) + 2.3*wstar_in(ig)**2))**2

       DO WHILE((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))

         pz0t=pz0tcomp(ig)
         IF (zu2(ig) .ne. 0.) THEN
            ! Richardson number formulation proposed by D.E. England et al. (1995)
          rib(ig) = (pg/pts(ig))
     &        *sqrt(zzlay(ig,1)*pz0(ig))
     &        *(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t)))
     &        *(ph(ig,1)-pts(ig))/(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))
         ELSE
!            print*,'warning, infinite Richardson at surface'
!            print*,pu(ig,1),pv(ig,1)
            rib(ig) = ric
         ENDIF

         z1z0=zzlay(ig,1)/pz0(ig)
         z1z0t=zzlay(ig,1)/pz0t

         cdn(ig)=karman/log(z1z0)
         cdn(ig)=cdn(ig)*cdn(ig)
         chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) 

         ! STABLE BOUNDARY LAYER :
         IF (rib(ig) .gt. 0.) THEN
            ! From D.E. England et al. (95)
            prandtl(ig)=1.
            if(rib(ig) .lt. ric) then
               ! Assuming alphah=1. and bh=bm for stable conditions :
               fm(ig)=((ric-rib(ig))/ric)**2
               fh(ig)=((ric-rib(ig))/ric)**2
            else
               ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
               fm(ig)=0.
               fh(ig)=0.
            endif
         ! UNSTABLE BOUNDARY LAYER :
         ELSE
            ! From D.E. England et al. (95)
            fm(ig)=sqrt(1.-lambda*bm*rib(ig))
            fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*
     &                     (1.-lambda*bm*rib(ig))**0.25
            prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/
     &             ((1.-lambda*bh*rib(ig))**0.5)
         ENDIF
 
        reynolds(ig)=karman*sqrt(fm(ig))
     &              *sqrt(zu2(ig))
     &              *pz0(ig)/(log(z1z0)*nu)
        pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3*
     &              (reynolds(ig)**0.25)*(prandtl(ig)**0.5))
        residual = abs(pz0t-pz0tcomp(ig))
        ite = ite+1

       ENDDO  ! of while
       pz0t=0.

! Drag computation:

         pcdv(ig)=cdn(ig)*fm(ig)
         pcdh(ig)=chn(ig)*fh(ig)
       
      ENDDO ! of ngrid

!------------------------------------------------------------------------
!------------------------------------------------------------------------
! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
!------------------------------------------------------------------------
!------------------------------------------------------------------------

! u* theta* computation

      DO ig=1,ngrid
         IF (rib(ig) .ge. ric) THEN
           ustar(ig)=0.
           tstar(ig)=0.
         ELSE
           ustar(ig)=sqrt(pcdv(ig))
     &        *sqrt(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))
           tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
     &        /sqrt(pcdv(ig))
         ENDIF
      ENDDO

! Interpolation:

      DO ig=1,ngrid
        IF(zout .lt. pz0tcomp(ig)) THEN
           u_out(ig,n)=0.
           Teta_out(ig,n)=pts(ig)

        ELSE
          IF (rib(ig) .ge. ric) THEN ! ustar=tstar=0  (and fm=fh=0)
           u_out(ig,n)=0
           Teta_out(ig,n)=pts(ig)
          ELSE
           u_out(ig,n)= ustar(ig)*log(zout/pz0(ig))/
     &(karman*sqrt(fm(ig)))

           Teta_out(ig,n)=pts(ig)+(tstar(ig)*sqrt(fm(ig))*log(zout/
     & (pz0tcomp(ig)))/
     &(karman*fh(ig)))
          ENDIF
        ENDIF

        IF (zout .lt. pz0(ig)) THEN
           u_out(ig,n)=0.
        ENDIF 

      ENDDO

! when using convective adjustment without thermals, a vertical potential temperature
! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
! interpolated at any height in the surface layer is theta at the first level.

!      IF ((.not.calltherm).and.(calladj)) THEN
!       Teta_out(:,n)=ph(:,1)
!       u_out(:,n)=(sqrt(cdn(:))*sqrt(pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1))
!     &                                /karman)*log(zout/pz0(:))
!      ENDIF
              T_out(:,n) = Teta_out(:,n)*(exp(
     &   (zout/zzlay(:,1))*(log(pplay(:,1)/ps))
     &                             )
     &                         )**rcp

      ENDDO   !of n=1,n_out


!------------------------------------------------------------------------
!------------------------------------------------------------------------
! PART III : WSTAR COMPUTATION
!------------------------------------------------------------------------
!------------------------------------------------------------------------

! Detection of the mixed-layer potential temperature
! ------------

! Nearest index for the pbl height

!      IF (calltherm) THEN

      pbl_height_index(:)=1

      DO k=1,nlay-1
         DO ig=1, ngrid
            IF (abs(zmax(ig)-zzlay(ig,k)) .lt. 
     &              abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN
               pbl_height_index(ig)=k
            ENDIF
         ENDDO
      ENDDO

! Potential temperature gradient

      dteta(:,nlay)=0.
      DO k=1,nlay-1
         DO ig=1, ngrid
         dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k))
         ENDDO
      ENDDO

! Computation of the pbl mixed layer temperature

      DO ig=1, ngrid
         ii=MINLOC(abs(dteta(ig,1:pbl_height_index(ig))))
         pbl_teta(ig) = ph(ig,ii(1))
      ENDDO


!------------------------------------------------------------------------
!------------------------------------------------------------------------
! PART IV : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
!------------------------------------------------------------------------
!------------------------------------------------------------------------

! We follow Spiga et. al 2010 (QJRMS)
! ------------

      DO ig=1, ngrid
         IF (zmax(ig) .gt. 0.) THEN
            x(ig) = zout/zmax(ig)
         ELSE
            x(ig) = 999.
         ENDIF
      ENDDO

      DO ig=1, ngrid
         ! dimensionless vertical heat flux
         IF (x(ig) .le. 0.3) THEN
            dvhf(ig) = ((-3.85/log(x(ig)))+0.07*log(x(ig)))
     &                                       *exp(-4.61*x(ig))
         ELSEIF (x(ig) .le. 1.) THEN
            dvhf(ig) = -1.52*x(ig) + 1.24
         ELSE
            dvhf(ig) = 0.
         ENDIF
         ! dimensionless vertical velocity variance
         IF (x(ig) .le. 1.) THEN
            dvvv(ig) = 2.05*(x(ig)**(2./3.))*(1.-0.64*x(ig))**2
         ELSE
            dvvv(ig) = 0.
         ENDIF
      ENDDO

      vhf(:) = dvhf(:)*hfmax(:)
      vvv(:) = dvvv(:)*(wstar_in(:))**2

!      ENDIF ! of if calltherm

      IF (dummy) THEN
      vhf(:)=0.
      vvv(:)=0.
      T_out(:,1)=pts(:)
      u_out(:,1)=0.
      !ustar is computed
      !star is computed
      ENDIF

      RETURN
      END

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine eofpb(inicoord,scena,pertm,pertr,rdeof,lon,lat,
     &                 levhi,levlow,levweight,pratio,day,name,ier)
      
c     Calculate an EOF perturbation on the 3d variable=name at 
c     longitude=lon, latitude=lat and day=otherday.
c     where otherday=day+(rdeof-0.5)*eof_twindow , with rdeof in [0:1]
c
c     Improved variability model...EOFs calculated in longitude-sigma plane
c
c         Updated feb-2008: added lon-lat-time interpolation of perturbations
c                           + fixed "dust storm case" (where data is only
c                            available over half of the year). E.M.
      implicit none

      include "constants_mcd.inc"
      
c     inputs
      logical       inicoord ! flag to (re-)set evaluation of lon/lat/time
                             ! index of encompassing neighbours
      integer       scena   ! dust scenario (4,5 or 6 == dust storm)
!      real          rdnos(dimnevecs) !Uniform deviates for a single profile
      real          rdeof            ! uniform deviate
      real          lon              !longitude east of point 
      real          lat              !latitude of point
      integer       levhi            !database level upper bound
      integer       levlow           !database level lower bound
      real          levweight(2)     !level weight for interpolation
c     (1) for linear in height (2) for linear in pressure
      real          pratio           !ratio of pressure to extreme value if out of range
      real          day              !model day [0:668.6]
      character*16  name             !name of variable
      
c     outputs
      integer   ier     ! error flag (0=OK, 1=NOK)
      real      pertm   ! perturbation corresponding to trend (not used!!)
      real      pertr   ! perturbation corresponding to random component
      
c     local variables
      logical firstcall              ! flag to signal initializations
      data firstcall/.true./
      save firstcall
      real lat_eof(dimlateo)         ! latitudes, along EOF grid
      real lon_eof(dimloneo)        ! longitudes, along EOF grid
      save lat_eof,lon_eof
      
      integer indlat(2) ! indexes of encompassing latitudes
      data indlat /0,0/ ! dummy initialization to get rid of compiler warning
      integer indlon(2) ! indexes of encompassing longitudes
      data indlon /0,0/ ! dummy initialization to get rid of compiler warning
      integer indtime(2) ! indexes of encompassing dates
      data indtime /0,0/ ! dummy initialization to get rid of compiler warning
      save indlat,indlon,indtime
      
      real rel_lat ! relative position of lat between encompassing latitudes
      data rel_lat /0.0/ ! dummy initialization to get rid of compiler warning
      real rel_lon ! relative position of lon between encompassing longitudes
      data rel_lon /0.0/ ! dummy initialization to get rid of compiler warning
      real rel_time ! relative position of time between encompassing times
      data rel_time /0.0/ ! dummy initialization to get rid of compiler warning
      save rel_lat,rel_lon,rel_time
      
      real marsyear
      parameter (marsyear=668.5921) ! number of sols in a martian year
      ! Time window over which EOF perturbations are picked:
      real eof_twindow ! window size (sols)
      parameter (eof_twindow=30.0)
      
      real otherday  ! model day at which values are sought
      ! values at 4 encompassing nodes and 2 encompassing times
      real values(4,2) !,smthvalues(4,2)
      ! bilinearly interpolated values, at 2 encompassing times
      real val(2) !,smthval(2)
      ! final values, interpolated in time
      real fval !,fsmthval
      
      integer       i,inode,ilon,ilat,itim
      real     norm    ! normalization factor
!      real     average !  average value
!      real     pc(dimnevecs) ! principal components
!      real     pcsmth(dimnevecs) ! smoothed principal components
!      real     pcpert(dimnevecs) ! 'perturbation' principal components
      real     pcpert(dimnevecs,4,2) ! 'perturbation' principal components
      save     pcpert
      real          evecs(dimnevecs)

!      write(out,*) ""
!      write(out,*) "Entering eofpb, name:",name

      ! 1. Initializations
      
      ! first call only initializations
      if (firstcall) then
       ! build latitudes of EOF grid
       ! (NB: dimlateo,latmaxeo,... are known from constants_mcd.inc)
       do i=1,dimlateo
         lat_eof(i)=latmaxeo-(i-1)*deltalateo
       enddo
       ! build longitudes of EOF grid
       do i=1,dimloneo
         lon_eof(i)=lonmineo+(i-1)*deltaloneo
       enddo
        firstcall=.false.
      endif

      ! initialize return code
      ier = 0
      

      ! 2. find indexes of encompassing points
      ! NB: Only need to do this once every call to call_mcd
      if (inicoord) then
      ! 2.1. find encompassing longitudes
       if (lon.ge.lon_eof(dimloneo)) then ! wraparound
        indlon(1)=dimloneo
        indlon(2)=1
        rel_lon=(lon-lon_eof(dimloneo))/deltaloneo
       else
        indlon(1)=1
        do i=1,dimloneo-1
          if ((lon.ge.lon_eof(i)).and.(lon.lt.lon_eof(i+1))) then
            indlon(1)=i
          endif
        enddo
        indlon(2)=indlon(1)+1
        rel_lon=(lon-lon_eof(indlon(1)))/deltaloneo
       endif ! of if (lon.ge.lon_eof(dimloneo))

      ! 2.2 Find encompassing latitudes
       if (lat.ge.lat_eof(1)) then ! beyond northernest latitude
        indlat(1)=1
        indlat(2)=1
        rel_lat=0
       elseif (lat.le.lat_eof(dimlateo)) then ! beyond southernest latitude
        indlat(1)=dimlateo
        indlat(2)=dimlateo
        rel_lat=0
       else ! find encompassing latitudes
        indlat(1)=1
        do i=1,dimlateo-1
          if ((lat.le.lat_eof(i)).and.(lat.gt.lat_eof(i+1))) then
          indlat(1)=i
          endif
        enddo
        indlat(2)=indlat(1)+1
        rel_lat=(lat_eof(indlat(1))-lat)/deltalateo
       endif ! of if (lat.ge.lat_eo(1)) elseif 

      ! 2.3 compute time at which values will be sought
       otherday=mod(marsyear+day+(rdeof-0.5)*eof_twindow,marsyear)
!      otherday=day ! temporary, for tests
!      write(out,*)'eofpb: day=',day,' otherday=',otherday

      ! 2.4 Find encompassing times
       if ((scena.ge.4).and.(scena.le.6)) then
        ! specific to dust storm scenarios (data in Ls 180-360)
        if (otherday.lt.373) then
         ! first day of series; no interpolation
         indtime(1)=1
         indtime(2)=1
         rel_time=0.0
        elseif (otherday.ge.(dimeoday-1)) then
         ! last day of series; no interpolation
         indtime(1)=dimeoday-372
         indtime(2)=dimeoday-372
         rel_time=0.0
        else
         indtime(1)=1
         do i=373,dimeoday-1
          if ((otherday.ge.i).and.(otherday.lt.(i+1))) then
            indtime(1)=i-372 ! because dust storm data starts at month #7
          endif
         enddo
         indtime(2)=indtime(1)+1
         rel_time=(otherday-372.0)-indtime(1)
        endif ! of if (otherday.lt.373) elseif (otherday.ge.(dimeoday-1))
       else ! data is known all year round (MY24,cold & warm scenarios)
        if (otherday.lt.1.0) then ! wraparound
         indtime(1)=dimeoday
         indtime(2)=1
         rel_time=otherday
        elseif (otherday.ge.(dimeoday-1)) then
        ! specific workaround to handle case of
        ! last day of year (which is not of unity length, in the true calendar)
         indtime(1)=dimeoday-1
         indtime(2)=dimeoday
         rel_time=(otherday-int(otherday))*1./0.6 ! so that rel_time in [0:1]
        else
         indtime(1)=1
         do i=1,dimeoday-1
          if ((otherday.ge.i).and.(otherday.lt.(i+1))) then
            indtime(1)=i
          endif
         enddo
         indtime(2)=indtime(1)+1
         rel_time=otherday-indtime(1)
        endif ! of if (otherday.lt.1.0)
       endif ! of if (scena.ge.4).and.(scena.le.6)
      
!        write(out,*) "eofpb: day=",day,' otherday=',otherday
!        write(out,*) "eofpb: indtime(1)=",indtime(1)
!        write(out,*) "eofpb: rel_time=",rel_time

      ! 2.5 Build pcpert(dimnevecs,4,2)
       do itim=1,2
        do inode=1,4
          if (inode.le.2) then ! node # 1 or # 2
            ilat=indlat(1)
          else ! node # 3 or # 4
            ilat=indlat(2)
          endif
          ! read PCs
!          do i=1,dimnevecs
!            pc(i)=tabpc(ilat,indtime(itim),i)
!          enddo
          ! read smoothed PCs
!          do i=1,dimnevecs
!            pcsmth(i)=tabpcsmth(ilat,indtime(itim),i)
!          enddo
          ! build pcpert, 'perturbation' component of principal component
          do i=1,dimnevecs
            ! Note: we multiply by 1.205 to preserve overall variance
!            pcpert(i,inode,itim)=1.205*(pc(i)-pcsmth(i))
            pcpert(i,inode,itim)=1.205*(tabpc(ilat,indtime(itim),i)
     &                                 -tabpcsmth(ilat,indtime(itim),i))
          enddo
        enddo ! of do inode=1,4
       enddo ! of do itim=1,2
      
      endif ! of if (inicoord)


      ! 3. compute values at encompassing nodes and time
      ! values are organized as follows:  v(1,i) v(2,i)
      !                                   v(4,i) v(3,i)

      do itim=1,2 ! loop on encompassing time
        do inode=1,4 ! loop on encompassing nodes
         if (inode.eq.1) then
          ilon=indlon(1)
          ilat=indlat(1)
         elseif (inode.eq.2) then
          ilon=indlon(2)
          ilat=indlat(1)
         elseif (inode.eq.3) then
          ilon=indlon(2)
          ilat=indlat(2)
         else
          ilon=indlon(1)
          ilat=indlat(2)
         endif
c     read normalisation factor, averages and eofs
c     normalisation factor is no longer the standard deviation
c
c     EOF perturbations are stored up to the last model level
         if (levlow.lt.dimlevs) then  ! we are entirely within EOF range
          if (name.eq.'u') then
            norm=tabeonormu(ilat)
            do i=1,dimnevecs
            evecs(i)= tabeou(ilon,ilat,levlow,i)
     &             +  (tabeou(ilon,ilat,levhi,i)
     &             -  tabeou(ilon,ilat,levlow,i))*levweight(1)
            enddo
!            average=tabeouave(ilon,ilat,levlow)
!     &             + (tabeouave(ilon,ilat,levhi)-
!     &                tabeouave(ilon,ilat,levlow))*levweight(1)
          elseif (name.eq.'v') then
            norm=tabeonormv(ilat)
            do i=1,dimnevecs
            evecs(i)= tabeov(ilon,ilat,levlow,i)
     &             +  (tabeov(ilon,ilat,levhi,i)
     &             -  tabeov(ilon,ilat,levlow,i))*levweight(1)
            enddo
!            average=tabeovave(ilon,ilat,levlow)
!     &             + (tabeovave(ilon,ilat,levhi)-
!     &                tabeovave(ilon,ilat,levlow))*levweight(1)
          elseif (name.eq.'temp') then
            norm=tabeonormt(ilat)
            if (levweight(1).eq.0.0) then ! happens a lot when
              ! eofpb is called from profi_eof
              do i=1,dimnevecs
                evecs(i)= tabeot(ilon,ilat,levlow,i)
              enddo
            else ! standard case do a vertical interpolation
              do i=1,dimnevecs
                evecs(i)= tabeot(ilon,ilat,levlow,i)
     &                 +  (tabeot(ilon,ilat,levhi,i)
     &                 -  tabeot(ilon,ilat,levlow,i))*levweight(1)
             enddo
            endif
!            average=tabeotave(ilon,ilat,levlow)
!     &             + (tabeotave(ilon,ilat,levhi)-
!     &                tabeotave(ilon,ilat,levlow))*levweight(1)
!          elseif (name.eq.'rho') then
!            norm=tabeonormr(ilat)
!!        write(out,*)'norm:',norm
!            do i=1,dimnevecs
!            evecs(i)= tabeorho(ilon,ilat,levlow,i)
!     &             +  (tabeorho(ilon,ilat,levhi,i)
!     &             -  tabeorho(ilon,ilat,levlow,i))*levweight(2)
c         write(out,*)'evecs(',i,'):',evecs(i)
!            enddo
          elseif (name.eq.'ps') then
            norm=tabeonormp(ilat)
            do i=1,dimnevecs
               evecs(i)= tabeops(ilon,ilat,i)
            enddo
!            average=tabeopsave(ilon,ilat)
          else
           if (output_messages) then
            write(out,*)'EOFPB Error: ',name,' unknown variable name'
           endif
           ier=1
           goto 9999
          endif
         else  ! we are above the top EOF level, most perturbations constant
          if (name.eq.'u') then
            norm=tabeonormu(ilat)
            do i=1,dimnevecs
              evecs(i)= tabeou(ilon,ilat,dimlevs,i)
            enddo
!            average=tabeouave(ilon,ilat,levlow)
          elseif (name.eq.'v') then
            norm=tabeonormv(ilat)
            do i=1,dimnevecs
              evecs(i)= tabeov(ilon,ilat,dimlevs,i)
            enddo
!            average=tabeovave(ilon,ilat,levlow)
          elseif (name.eq.'temp') then
            norm=tabeonormt(ilat)
            do i=1,dimnevecs
              evecs(i)= tabeot(ilon,ilat,dimlevs,i)
            enddo
!            average=tabeotave(ilon,ilat,levlow)
!          elseif (name.eq.'rho') then
!            norm=tabeonormr(ilat)
!            do i=1,dimnevecs
!              evecs(i)= tabeorho(ilon,ilat,dimlevs,i)
!            enddo
          elseif (name.eq.'ps') then
            norm=tabeonormp(ilat)
            do i=1,dimnevecs
              evecs(i)= tabeops(ilon,ilat,i)
            enddo
!            average=tabeopsave(ilon,ilat)
          else
            if (output_messages) then
             write(out,*)'EOFPB Error: ',name,' unknown variable name'
            endif
            ier=1
            goto 9999
          endif ! of if (name.eq.'u') elseif ...
         endif ! of if (levlow.lt.dimlevs)
c
c     calculate perturbation
c
!          pertm=0.0
!          pertr=0.0
          values(inode,itim)=0.0 ! initialization
!          smthvalues(inode,itim)=0.0 ! initialization
          do i=1,dimnevecs !72 !dimnevecs
!            smthvalues(inode,itim)=smthvalues(inode,itim)+
!     &                         pcsmth(i)*evecs(i)
            ! here values(:,:) is used to store the perturbation
            values(inode,itim)=values(inode,itim)+
     &                         pcpert(i,inode,itim)*evecs(i)
!            pertm=pertm+pcsmth(i)*evecs(i)
!       write(out,*)'i:',i,' pcvar:',pcvar(i),' rdnos:',rdnos(i),
!     &            ' evecs:',evecs(i)
!            pertr=pertr+pcvar(i)*rdnos(i)*evecs(i)
!            write(out,*)' pertr:',pertr
          enddo
c
c     renormalise and add average
c
!          pertm=pertm*norm
!          pertr=pertr*norm
!          smthvalues(inode,itim)=average+smthvalues(inode,itim)*norm
          values(inode,itim)=values(inode,itim)*norm
          ! values(:,:) now contains the perturbed field
        enddo ! do inode=1,4
        
        ! bilinear interpolation using four neighbours
!        if ((name.eq.'rho').or.(name.eq.'ps')) then
!        if (name.eq.'ps') then
          ! interpolate the log() ! NO, because these are deviations to mean
!          val(itim)=(1.-rel_lon)*(1.-rel_lat)*log(values(1,itim))
!     &              +rel_lon*(1.-rel_lat)*log(values(2,itim))
!     &              +rel_lon*rel_lat*log(values(3,itim))
!     &              +(1.-rel_lon)*rel_lat*log(values(4,itim))
!          val(itim)=exp(val(itim))
!        else ! bilinear interpolation of the variable
          val(itim)=(1.-rel_lon)*(1.-rel_lat)*values(1,itim)
     &              +rel_lon*(1.-rel_lat)*values(2,itim)
     &              +rel_lon*rel_lat*values(3,itim)
     &              +(1.-rel_lon)*rel_lat*values(4,itim)
!        endif
      enddo ! do itim=1,2

      ! 4. linear interpolation in time
      fval=(1.0-rel_time)*val(1)+rel_time*val(2)
!      fsmthval=(1.0-rel_time)*smthval(1)+rel_time*smthval(2)

      pertr=fval

!      write(out,*)'eofpb: pertr=',pertr

      ! set inicoord to false, since rel_lat, rel_lon, rel_time ... are
      ! now known (and won't change aduring current call to call_mcd).
      inicoord=.false.

      return

c     error handling
 9999 pertm = 0.0
      pertr = 0.0
      return

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine profi_eof(inicoord,scena,profile,rdeof,lon,lat,day,
     &                     toplevpert,name,profile_pert,ier)

c     Add EOF perturbations to all vertical levels of a (known) profile

      implicit none

      include "constants_mcd.inc"

c     inputs
      logical inicoord ! flag to (re-)set evaluation of lon/lat/time
                       ! passed on to eofpb
      integer scena ! dust scenario 
      real    profile(dimlevs) ! the GCM profile at GCM levels
!      real    rdnos(dimnevecs) ! uniform deviates
      real    rdeof ! uniform deviate, for EOF perturbations
      real    lon   ! east longitude (degrees)
      real    lat   ! latitude (degrees)
      real    day   ! model day, in [0:668.6]
      integer toplevpert ! model layer up to which perturbation will be added
      character*16 name !name of variable
      
c     outputs
      real    profile_pert(dimlevs) ! perturbed profile at GCM levels
      integer ier     ! error flag (=0 if OK)

c local variables
      integer i
      real pertm ! perturbation corresponding to trend (from eofpb)
      real pertr ! perturbation corresponding to random component (from eofpb)
      integer levlo ! local value for call to eofpb
      integer levhi ! local value for call to eofpb
      real levweight(2) ! local value for call to eofpb
      real pratio ! local (unused) value for call to eofpb
      
      integer toppertlev ! level up to which perturbation will be added
      
      if (toplevpert.gt.dimlevs) then
      ! in case input toplevpert is greater than dimlevs
        toppertlev=dimlevs
      else
        toppertlev=toplevpert
      endif
      
      if (name.eq.'temp') then
        ! add perturbation to levels up to toppertlev
        do i=1,toppertlev
         if (i.lt.dimlevs) then ! general case
          levweight(1)=0.0 ! all the weight on levlo
          levlo=i
          levhi=i+1
          ! get EOF perturbation for level levlo
          call eofpb(inicoord,scena,pertm,pertr,rdeof,lon,lat,
     &               levhi,levlo,levweight,pratio,day,name,ier)
          ! add EOF perturbation
          profile_pert(levlo)=profile(levlo)+pertr
         else ! special case if last level is also perturbed
          levweight(1)=1.0 ! i.e. all the weight on levhi
          levlo=dimlevs-1
          levhi=dimlevs
          ! get EOF perturbation for level levhi
          call eofpb(inicoord,scena,pertm,pertr,rdeof,lon,lat,
     &               levhi,levlo,levweight,pratio,day,name,ier)
          ! add EOF perturbation
          profile_pert(levhi)=profile(levhi)+pertr
         endif ! of if (i.lt.dimlevs)
        enddo
        
        ! levels above toppertlev are not perturbed
        do i=toppertlev+1,dimlevs
          profile_pert(i)=profile(i)
        enddo
      else
        if (output_messages) then
          write(out,*) 'PROFI_EOF Error: unknowm name:',name
        endif
        ier=1
      endif ! of if (name.eq.'temp')
      
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      function gasdev(idum)

c     Return a Gaussian deviate with zero mean and unit standard deviation.
c     Uses function ran1()
c     If idum is negative, then deviates are (re-)initialized.
c Note: When initialized, the function actually computes 2 deviates,
c       'gset' and 'gasdev'. It then returns 'gasdev' and saves 'gset'
c       which will be returned at the next call of the function.
      implicit none

c     inputs
      integer idum    ! Seed for random number generator
                      ! is altered if calls to function ran1() occur

      real    gasdev

c     local variables
      integer iset    ! flag (triggers the generation of 2 deviates if =0)
      data    iset/0/
      save    iset
      real    v1      ! random number (between -1 and 1)
      real    v2      ! random number (between -1 and 1)
      real    r
      real    fac
      real    gset    ! extra deviate (saved for next call to function)
      save    gset

c called function
      real    ran1

      if (idum.lt.0) iset=0  ! Reinitialize deviates generation
      if (iset.eq.0) then ! compute 2 deviates, gasdev and gset
1       v1=2.*ran1(idum)-1.
        v2=2.*ran1(idum)-1.
        r=v1*v1+v2*v2
        if ((r.ge.1.).or.(r.eq.0.)) goto 1
        fac=sqrt(-2.*log(r)/r)
        gset=v1*fac
        gasdev=v2*fac
        iset=1  ! set flag to signal there is an extra deviate at hand
      else   ! return saved deviate 'gset'
        gasdev=gset
        iset=0  ! set flag to signal there is no extra deviate at hand
      endif

      return
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine getsi(xz,zkey,lon,lat,oroheight,areoid,
     &                 ps,sigma,utime,dset,
     &               levhi,levlow,levweight,pratio,ier,itimint,wl,wh,
     &    R_gcm,temp_gcm,zareoid,zradius,zsurface,zpressure,zmradius,
     &    sheight)

c     Find the nearest datafile levels and vertical interpolation weights at
c     longitude=lon, latitude=lat and vertical position xz (which may be
c     any of 'above areoid' or 'above surface'
c      or 'radius from center of planet' 
c      or 'altitude above mean radius'
c      or 'Pressure' coordinate)
c     Last Updated feb-2008: added 'zmradius' input
c                           +improved Tmean interpolation   E.M. 
      implicit none

      include "constants_mcd.inc"
      
c     inputs

      real          xz ! altitude (see zkey)
      integer       zkey             
c  zkey    :                   switch to choose which z variable
c                               is input and used to find other three :
c           zkey=1  zradius (m)  -> zareoid, zsurface, zpressure, zmradius
c           zkey=2  zareoid (m)  -> zradius, zsurface, zpressure, zmradius
c           zkey=3  zsurface (m) -> zradius, zareoid, zpressure, zmradius
c           zkey=4  zpressure (Pa) -> zradius, zareoid, zsurface, zmradius
c           zkey=5  zmradius (m) -> zradius, zareoid, zsurface, zpressure
c           zkey=-5  case of getsi call by grwpb 
      real          lon            !east longitude
      real          lat            !latitude
      real          oroheight      !height of surface above reference areoid (m)
      real          areoid         ! reference areoid (m)
      real          ps             ! surface pressure (Pa)
      real          sigma(dimlevs) ! sigma levels
      real          utime          !Universal time (0 to 24 hrs)=local time at lon=0
      character*(*) dset           ! Path to Dataset
      integer       itimint        !seasonal interpolation flag
      real          wl, wh         !seasonal interpolation weights
      real          R_gcm(dimlevs) ! R at GCM levels
      real          temp_gcm(dimlevs) ! temperature at GCM levels

c     outputs
      integer       levhi          !database level upper bound
      integer       levlow         !database level lower bound
      real          levweight(2)   !level weight for interpolation
c     (1) for linear in height (2) for linear in pressure
      real          pratio         !p/p(top, bottom) for extrapolation,
                                   ! 1 if in range
      integer       ier            !error flag (0=OK, not 0 =NOK)

      real zradius     ! distance to center of planet
      real zareoid     ! height above areoid
      real zsurface    ! height above surface
      real zpressure   ! pressure
      real zmradius    ! altitude above mean radius of planet
      real sheight(dimlevs) ! altitude above surface of sigma levels

c     local variables
      real marsradius ! mean Mars radius (m)
      parameter (marsradius=3.396E6)
      real          height         !height of point above areoid (m)
      real          absheight         !height of point above surface (m)
      integer       l
      real          zsig
      
c     Gravity on mean areoid
c     The areoid is defined as a surface of constant gravitational plus
c     rotational potential. The inertial rotation rate of Mars is assumed
c     to be 0.70882187E-4 rad/s. This potential is the mean value at the
c     equator at a radius of 3396.000 km, namely 12652804.7 m^2/s^2,
c     calculated from Goddard Mars Gravity Model mgm1025
c     [LEMOINEETAL2001] evaluated to degree and order 50.
      real          g, g0, a0
      parameter (a0=3396.E3,g0=3.7257964)
c      real          R 
c      data          R /191.1/
c     R replaced by Rnew  (varies with altitude)
!      real          Rnew(dimlevs)
      real          Tmean          ! "mean" temperature of a layer
      real          Rogct(dimlevs) 
      real          vmr_gcm(dimlevs)  ! vmr at GCM levels
      real          vmr(nbvar3d)
      integer       i,igas,ierr
      real          dz,dz_max,dz_min,xsig,z

      ier=0


c     Calculate altitude above the surface of each model layer: sheight(l)
c     integrate hydrostatic equation
c     Rogct is the R/g variable depending on dimlevs 
      g = g0*(a0/(a0+oroheight))**2
      Rogct(1) = R_gcm(1)/g          !Rnew(1)/g
      Tmean = temp_gcm(1)            !t(1)
      sheight(1)= - Rogct(1)*temp_gcm(1)*log(sigma(1))

      do l=2, dimlevs
        if (temp_gcm(l).ne.temp_gcm(l-1)) then
         Tmean = real((dble(temp_gcm(l)-temp_gcm(l-1))) /
     &                 dlog(dble(temp_gcm(l))/dble(temp_gcm(l-1))))
         ! NB: double precision must be used here in case
         ! temp_gcm(l) and temp_gcm(l-1) are almost equal  
        else
         Tmean = temp_gcm(l)
        end if
        g = g0*(a0/(a0+oroheight+sheight(l-1)))**2
        Rogct(l) = R_gcm(l)/g
        sheight(l) = sheight(l-1) - Rogct(l)*Tmean*
     &               log(sigma(l)/sigma(l-1))
      end do

cccccccccccccccccccccccccccccccccccccccccccccccc     
c     height calculation      
      if (zkey.eq.1) then ! input xz is distance to center of planet
              zradius=xz
      elseif (zkey.eq.2) then ! input xz is height above areoid
              zareoid=xz
      elseif (zkey.eq.3) then ! input xz is altitude above surface
            zsurface=xz
      elseif (zkey.eq.4) then ! input xz is atmospheric pressure
          zpressure=xz
c         compute zsurface
          if(zpressure.gt.ps) then
            if (output_messages) then
              write(out,*)'GETSI Error: underground object '
            endif
            ier=17
            goto 9999
          end if
          
          zsig = real(zpressure)/ps
          if (zsig.gt. sigma(1)) then
             zsurface= -log(zsig)*R_gcm(1)*temp_gcm(1)/g0
          else if (zsig.lt.sigma(dimlevs) ) then
             g = g0*(a0/(a0+oroheight+sheight(dimlevs)))**2
!             zsurface= sheight(dimlevs)
!     &        -log(zsig/sigma(dimlevs))*R_gcm(dimlevs)
!     &                                 *temp_gcm(dimlevs)/g
             ! Compute zsurface
             do igas=1,nbvar3d
              ! extract the vmr of various gases
              if(ma3d(igas).lt.0) cycle   
              call profi(vmr_gcm,lon,lat,utime,varname3d(igas),
     &            ierr,itimint,wl,wh,0,dimlevs,dimlevs)
              vmr(igas) = vmr_gcm(dimlevs)
             enddo 
c
c Search for the altitude corresponding to the given pressure
c by the bisection method
c
             dz_max = 100000.0e+3 ! max possible delta
             dz_min = 0.0   

             do i=0,100 ! max number of iterations
              dz = 0.5*(dz_max+dz_min)
              z  = (a0+oroheight+sheight(dimlevs))*dz/
     &             (a0+oroheight+sheight(dimlevs)+dz) ! z=r0*dz/(r0+dz)
              xsig = 0.0

              do igas=1,nbvar3d
               if(ma3d(igas).lt.0) cycle   
               xsig = xsig + vmr(igas)*sigma(dimlevs)*
     &                exp(-g*z/(8314.4598/ma3d(igas)*temp_gcm(dimlevs)))               
              enddo

              if(abs(xsig-zsig)/zsig.lt.0.05) goto 5555
              if(xsig.gt.zsig) then 
               dz_min = dz
              else 
               dz_max = dz
              endif 
             enddo !i
   
             ! If the loop repeat to the end this means that 
             ! for a given range and precision the conversion of altitude
             ! has failed. Further the closest altitude is used.
             !    
             if (output_messages) then
              write(out,*)'CALL_MCD Error: conversion of altitude',
     &               'failed in getsi().'
             endif   
   
 5555        continue    
             zsurface= sheight(dimlevs)+dz
          else
            do  l=1,dimlevs-1
              if ((zsig.ge.sigma(l+1)).and.(zsig.le.sigma(l))) then
                zsurface= sheight(l) + (sheight(l+1) -sheight(l)) 
     &            * log(zsig/sigma(l)) / log(sigma(l+1)/sigma(l)) 
                exit
              end if
            end do
          end if
      elseif (zkey.eq.5) then ! input xz is altitude above mean Mars radius
        zmradius=xz
      endif      
      
c absheight = zsurface for the first call of getsi       
c absheight = xz for call to getsi from  grwpb      
      if (zkey.eq.-5) then ! call from grwpb
        absheight=xz
      else
        ! compute required zradius,zareoid,zsurface and/or zmradius
        if (zkey.eq.1.) then ! zradius is known
          zareoid=zradius-areoid
          zsurface=zareoid-oroheight
          zmradius=zradius-marsradius
        elseif (zkey.eq.2) then ! zareoid is known
          zradius=zareoid+areoid
          zsurface=zareoid-oroheight
          zmradius=(areoid-marsradius)+zareoid
        elseif ((zkey.eq.3).or.(zkey.eq.4)) then ! zsurface is known
          zareoid=zsurface+oroheight
          zradius=zareoid+areoid
          zmradius=(areoid-marsradius)+zareoid
        elseif (zkey.eq.5) then ! zmradius is known
          zradius=marsradius+zmradius
          zareoid=zmradius-(areoid-marsradius)
          zsurface=zareoid-oroheight
        endif
        
        height=zareoid              ! height above areoid
        absheight=height-oroheight  ! height above surface
      endif

ccccccccccccccccccccccccccccccccccccccccccccccccccc      
      
c     find levhi, levlow and levweight
c     compute g (at height=zareoid)
      g = g0*(a0/(a0+oroheight+absheight))**2
      if (absheight.lt.sheight(1)) then
c       below the lowest layer
        levhi=1
        levlow=1
        levweight(1)=1.
        levweight(2)=0.
        pratio=exp((sheight(1)-absheight)*g/(R_gcm(1)*temp_gcm(1)))
      elseif (absheight.ge.sheight(dimlevs)) then
c       above the top layer
        levhi=dimlevs
        levlow=dimlevs
        levweight(1)=0.
        levweight(2)=1.
!        pratio=exp((sheight(dimlevs)-absheight)*g/(R_gcm(dimlevs)
!     &                                             *temp_gcm(dimlevs)))
c z = (r1/r0)**2*r0*dz/(r0+dz) where r1 is alt where g was calculated 
        z = (a0+oroheight+absheight)**2/(a0+oroheight+sheight(dimlevs))*
     &      (absheight-sheight(dimlevs))/(a0+oroheight+absheight)
        pratio = 0.0
        do igas=1,nbvar3d
         ! Pressure is computed assuming hydrostatic for each species
         if(ma3d(igas).lt.0) cycle   
         call profi(vmr_gcm,lon,lat,utime,varname3d(igas),
     &            ierr,itimint,wl,wh,0,dimlevs,dimlevs)
         pratio = pratio+vmr_gcm(dimlevs)*
     &          exp(-z*g/(8314.4598/ma3d(igas)*temp_gcm(dimlevs)))
        enddo        
      else
c       general case: sheight(1)<=absheight<=sheight(dimlevs)
       do l=1,dimlevs-1
        if ((absheight.ge.sheight(l))
     &     .and.(absheight.lt.sheight(l+1))) then
          levhi=l+1
          levlow=l
          levweight(1)=(absheight-sheight(levlow))
     &                /(sheight(levhi)-sheight(levlow))
          levweight(2)=(1.-(sigma(levhi)/sigma(levlow))**levweight(1))
     &                /(1.-(sigma(levhi)/sigma(levlow)))
          pratio=1.
        endif
       enddo
      endif

c     Compute pressure
      zpressure=ps*(sigma(levlow)+(sigma(levhi)-sigma(levlow))
     &                              *levweight(2))*pratio

      return

c     Errror handling
 9999 levhi=0
      levlow=0
      levweight(1)=0.
      levweight(2)=0.
      pratio=0.

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine grid4(lon,lat,dlon,dlat,t,u)

c     Given longitude=lon and latitude=lat find the nearest 4 horizontal
c     gridpoints in the database and bilinear interpolation weights.
c
c     Grid points are arranged as follows:
c              4 3
c              1 2 
      
      implicit none

      include "constants_mcd.inc"

c     inputs
      real    lon      !east longitude of point
      real    lat      !latitude of point

c     outputs
      integer dlon(4)  !index, along longitudes, of database points
      integer dlat(4)  !index, along latitudes, of database points
      real    t        !weight (normalized longitudinal distance to point 1)
      real    u        !weight (normalized latitudinal distance to point 1)

c     local variables
      real    alon,alat
      integer i

c     wraparound
      if (lon.ge.lonmax) then
         dlon(1)=dimlon
         dlon(2)=1
         dlon(3)=1
         dlon(4)=dimlon
         t=(lon-lonmax)/deltalon
      else
         alon=lonmin
         do i=1,dimlon-1
            if ((lon.ge.alon).and.(lon.lt.alon+deltalon)) then
               dlon(4)=i
               dlon(3)=i+1
               dlon(1)=i
               dlon(2)=i+1
               t=(lon-alon)/deltalon
            endif
            alon=alon+deltalon
         enddo
      endif

      if (lat.le.latmin) then
         dlat(4)=1
         dlat(3)=1
         dlat(1)=1
         dlat(2)=1
         u=0.
      elseif (lat.ge.latmax) then
         dlat(4)=dimlat
         dlat(3)=dimlat
         dlat(1)=dimlat
         dlat(2)=dimlat
         u=0.
      else
         alat=latmin
         do i=1,dimlat-1
            if ((lat.ge.alat).and.(lat.lt.alat+deltalat)) then
               dlat(4)=i+1
               dlat(3)=i+1
               dlat(1)=i
               dlat(2)=i
               u=(lat-alat)/deltalat
            endif
            alat=alat+deltalat
         enddo       
      endif

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine grwpb(dset,pert,rangw1,rangw2,dev,lamda,lon,lat,
     &                 absheight,rho,u,v,ps,sigma,utime,name,ier,
     &                 itimint,wl,wh,R_gcm,temp_gcm,sheight,
     &                 p_pgcm,ps_psgcm)

c     Small scale gravity wave perturbation model.
c     Computes perturbation on u, v, theta or rho, with wave phase factor.
c
c     Last updated Feb. 2008: improved behaviour when near surface (in order
c                             to get rid of 'jumps' in perturbations when run
c                             at very fine resolution)
      implicit none

      include "constants_mcd.inc"

c     inputs
      character*(*)     dset     ! dataset directory
      real         dev     !phase of perturbation (rand number between 0 and 1)
      real         rangw1,rangw2
      real         lamda   !vertical wavelength of g.w. (m)
      real         lon     !east longitude of point
      real         lat     !latitude of point
      real         absheight  !height (above surface) of point
      real         rho     !density at height=height
      real         u       !zonal wind at height=height
      real         v       !meridional wind at height=height
      real         ps      ! surface pressure
      real         sigma(dimlevs) ! sigma levels 
      real         utime   !Universal time (0. to 24. hrs) = local time at lon=0
      integer      itimint !seasonal interpolation flag
      real         wl,wh   !seasonal interpolation weights
      real         R_gcm(dimlevs) ! R at GCM levels
      real         temp_gcm(dimlevs) ! temperature at GCM levels
      real         sheight(dimlevs) ! altitude of GCM sigma levels
      real         p_pgcm(dimlevs) ! high res to GCM pressure ratios
      real         ps_psgcm ! High res to GCM surface pressure ratio 
      character*16 name    !name of variable to perturb

c     outputs
      integer      ier     !error flag (0=OK, 1=NOK)
      real         pert    !perturbation

c     local variables
      ! zareoid,zpressure,zradius,zsurface only local
      ! we do not care of those variables here
      real         xz,zareoid,zpressure,zradius,zsurface,zmradius   
      integer      zkey
      integer      ierr
      real         a, aprime
      real         rho0,rho1
      real         u0,u1
      real         v0,v1
      real         tmpheight
      integer      levlow,levhi
      real         levweight(2)
      real         pratio
      real         dz,dz_max
      real         sig
      real         umag
      real         hmax    !for perturbations above hmax, use amplitude at hmax (m)
      parameter   (hmax=100.E3) 
      real         usat    !for windspeeds below usat, wave amplitude saturates
      Parameter   (usat=0.5)
      real         dalr    !dry adiabatic lapse rate (K/m)
      parameter   (dalr=4.5E-3)
      real         pi 
      parameter   (pi=3.14159265359)
      character*16 tmpname
      real         h,Dist
      real         lamda_D
      real, save :: prevrangw1=-999.
      real, save :: lat0,lon0
      real, parameter :: degtorad=pi/180.
      
!      write(out,*)""
!      write(out,*) "Entering grwpb, name:",name
!      write(out,*) "                dev=",dev
!      write(out,*) "                lamda=",lamda
!      write(out,*) "                absheight=",absheight
!      write(out,*) "                rho=",rho
!      write(out,*) "                u=",u
!      write(out,*) "                v=",v

      ier=0

c     1. get rho0, u0, v0
c        use mean data from level 1, not height=0.0
      levlow=1
      levhi=1
      levweight(1)=1.
      levweight(2)=1.
      pratio=1.0
      tmpname='u'
      call var3d(u0,lon,lat,sheight(1),levhi,levlow,levweight,pratio,
     &           sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &           itimint,wl,wh)
      tmpname='v'
      call var3d(v0,lon,lat,sheight(1),levhi,levlow,levweight,pratio,
     &           sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &           itimint,wl,wh)
      tmpname='rho'
      call var3d(rho0,lon,lat,sheight(1),levhi,levlow,levweight,pratio,
     &           sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &           itimint,wl,wh)

!      write(out,*) "grwpb: u0=",u0," v0=",v0," rho0=",rho0 
      
c     get sub-grid scale variance
      tmpname='substd'
      call var2d(sig,lon,lat,1.0,tmpname,ierr,itimint,wl,wh,1.0)

c     2. Compute delta z, the vertical displacement
c        (note that here we assume N=N0, i.e. no change in Brunt Vaisala 
c          frequency between altitude 'absheight' and 1st layer)
      if (absheight.le.hmax) then
         umag=sqrt(u**2+v**2)
c        the wave amplitude becomes large as the wind speed becomes small
c        and the wave should saturate
         if (umag.lt.usat) umag=usat
c        *********************               
c         theoretical version:  
          dz=sig*sqrt((rho0*sqrt(u0**2+v0**2))/(rho*umag))
c        *********************               
c         Version 3 :  
c        dz=1.E3*sqrt((rho0*sqrt(u0**2+v0**2)*sig*1.E-3)/(rho*umag))
c        *********************               
         tmpheight=absheight
      else
c         perturbation above hmax km:
c         Using the same amplitude and same perturbation
c         as at hmax km for T,u,v. The perturbation for rho is
c         scaled on the T perturbation for rho at height=height     
         tmpheight=hmax          ! set height for perturbation calculation
c         get info on variable at hmax m:
         xz=tmpheight
         zkey=-5
         call getsi(xz,zkey,lon,lat,0.0,0.0,ps,sigma,utime,dset,
     &              levhi,levlow,levweight,pratio,ierr,itimint,wl,wh,
     &              R_gcm,temp_gcm,zareoid,zradius,zsurface,zpressure,
     &              zmradius,sheight)
         tmpname='rho'
         call var3d(rho1,lon,lat,absheight,levhi,levlow,levweight,
     &          pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &              itimint,wl,wh)
         tmpname='u'
         call var3d(u1,lon,lat,absheight,levhi,levlow,levweight,
     &          pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &              itimint,wl,wh)
         tmpname='v'
         call var3d(v1,lon,lat,absheight,levhi,levlow,levweight,
     &          pratio,sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &              itimint,wl,wh)
         umag=sqrt(u1**2+v1**2)
c        the wave amplitude becomes large as the wind speed becomes small
c        and the wave should saturate
         if (umag.lt.usat) umag=usat
c        dz=(rho0*sqrt(u0**2+v0**2)*sig)/(rho1*umag)
c        dz=sqrt(dz)
          dz=sig*sqrt((rho0*sqrt(u0**2+v0**2))/(rho1*umag))
      end if
!     lat0,lon0 coordinates from which the distance Dist is calculated
!     remain the same as long as the seed is not changed
      if (rangw1.ne.prevrangw1) then
         lat0=lat
         lon0=lon
         prevrangw1=rangw1
      endif      
!      Horizontal wavelength proportional to the vertical one
      lamda_D=10.*lamda
      h=(sin((lat-lat0)*degtorad/2.))**2+cos(lat0*degtorad)
     &  *cos(lat*degtorad)*(sin((lon-lon0)*degtorad/2.))**2
      Dist=(2.*3390.E3)*asin(sqrt(h))
c     apply simple test for saturation (require theta_z > 0)
c     and compute wave
!      write(*,*) "grwpb: lamda=",lamda
      dz_max=lamda/(2.*pi)
!     Random saturation to make the maximum amplitude of the wave vary
      dz_max=dz_max*(1.+3.*rangw2)       
      if (dz.gt.dz_max) then
         dz=dz_max
      endif
      dz=dz*sin(2.*pi*dev+(2.*pi*absheight)/lamda
     &                   +(2.*pi*Dist)/lamda_D)  

c     3. Find perturbation from change delta z in mean profile
c     (more accurate than taking first derivative).
c     Using an oroheight of 0km in getsi here implies an insignificant
c     error in the calculated value of g used to find the distance to
c     the perturbation, very much smaller than other assumptions made
c     (interpolation, adiabatic, constant g, etc.).

c     don't allow dz to take us below the surface
!      if ((tmpheight+dz).le.0.) then
      if ((tmpheight+dz).le.0) then
! old version rescale underground to near-surface
!         dz = -0.99*tmpheight
! new version: mirror underground depths to above surface heights
        dz=-tmpheight-(tmpheight+dz)
      endif

! But, avoid being taken too close to the surface
      if ((tmpheight+dz).le.1000.) then
! map dz so that (tmpheight+dz) spans [100:1000] instead of [0:1000]
        dz=(((1000.0-100.0)/1000.0)*(tmpheight+dz)+100.0)-tmpheight
      endif

!      write(*,*) "grwpb: finally dz=",dz," tmpheight+dz=",tmpheight+dz

      tmpname = name
      if (name.eq.'rho') then
        ! density perturbation is computed from temperature perturbation
        ! (see below)
        tmpname= 'temp'
      endif
      
      ! get 'a': value of variable at altitude 'tmpheight'
      xz=tmpheight
      zkey=-5
      call getsi(xz,zkey,lon,lat,0.0,0.0,ps,sigma,utime,dset,
     &           levhi,levlow,levweight,pratio,ierr,itimint,wl,wh,
     &           R_gcm,temp_gcm,zareoid,zradius,zsurface,zpressure,
     &           zmradius,sheight)
      call var3d(a,lon,lat,xz,levhi,levlow,levweight,pratio,
     &           sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &           itimint,wl,wh)
      
      ! get 'aprime': value of variable at altitude 'tmpheight+dz'
      xz=tmpheight+dz
      call getsi(xz,zkey,lon,lat,0.0,0.0,ps,sigma,utime,dset,
     &           levhi,levlow,levweight,pratio,ierr,itimint,wl,wh,
     &           R_gcm,temp_gcm,zareoid,zradius,zsurface,zpressure,
     &           zmradius,sheight)
      call var3d(aprime,lon,lat,xz,levhi,levlow,levweight,pratio,
     &           sheight,p_pgcm,sigma,ps_psgcm,utime,tmpname,ierr,
     &           itimint,wl,wh)
      
      ! perturbation is just the difference between values at altitudes
      ! tmpheight and (tmpheight+dz) (with corrections for temperature
      ! and density)
      pert = aprime - a


c     correction for perturbation to potential temperature
      if (tmpname.eq.'temp') then
        pert = pert + dalr*dz
      endif

      ! reduce 'pert' if below 2000.0m (so that it smoothly goes to zero
      ! at surface)
      if (absheight.le.2000.0) then
        pert=pert*sin((absheight/2000.0)*(pi/2.))**2
      endif

c     correction for perturbation on density
c     pert(rho) = -rho*pert(T)/(T+pert(T))
      if (name.eq.'rho') then
         pert = -rho*pert/(a+pert)
      end if 

      ! if near the surface then quench perturbations

!      write(out,*) "grwpb: pert=",pert
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine loadvar_mcd(nf,nf_up,nf2,nf_up2,typevar,
     &                   aps,bps,sigma,pseudoalt,ier)

c     load arrays corresponding to the variable typevar

      implicit none
      
      include "netcdf.inc"
      include "constants_mcd.inc"

c variable typevar to load :
c 2D LOW mean = tsurf,ps,co2ice,fluxsurf_lw,fluxtop_lw,fluxsurf_sw,
c               fluxtop_sw,dod,col_h2ovapor,col_h2oice! 2D values
c               q2,vmr_h2ovapor,vmr_h2oice,vmr_o3,temp,u,v,rho ! 3D low (nf=unem)
c UP     mean = temp,u,v,w,rho,vmr_o,vmr_co2,_vmr_co !3D up (nf=unem_up)
c        sigm = sigma values of mean field file (nf=unetm)
c        orog = orography values of mountain file (nf=unet)
c        eofs = normu,normu,normt,normp,normr,u,v,t,rho,ps,pcsmth,pcvar
c               of eof field file (nf=ueof) 
c        grwp = substd values of mountain file (nf=unet)
c 2D&LOW stdv = tsdtsurf,tsdps,tsddod, !2D
c               tsdtemp,tsdsu,tsdsv,tsdrho !3D low values of
c               standard deviation field file (nf=unetsd)
c 2D&LOW rms  = rmstsurf,rmsps,rmsco2ice, !2D
c               rmstemp,rmsu,rmsv,rmsrho !3D low values of
c               rms field file (nf=unetsd)
c UP     stdv = tsdtemp,tsdsu,tsdsv,tsdw,tsdrho !3D up values of 
c               standard deviation field file (nf=unetsd_up)
c UP     rms  = rmstemp,rmssu,rmssv,rmsw,rmsrho !3D up values of 
c               rms field file (nf=unetsd_up)

c     inputs
      integer      nf,nf_up,nf2,nf_up2 !NetCDF file IDs
      character(len=*) :: typevar ! variable 'type'
      character*16 name    ! variable name

c     outputs      
      real aps,bps,sigma,pseudoalt
      dimension aps(dimlevs)        ! hybrid coordinate
      dimension bps(dimlevs)        ! hybrid coordinate
      dimension sigma(dimlevs)    ! sigma levels
      dimension pseudoalt(dimlevs)  ! pseudo altitude
      integer      ier              ! error flag (0=OK, 1=NOK)

c     local variables     
      integer      ierr
      integer      varid ! NetCDF variable ID
      integer      l,k,iloop,jloop
      real         aps_low(dimlevs),bps_low(dimlevs) ! hybrid coord.
      real         aps_up(dimlevs) ! hybrid coord. upper atm.

      ! 2D variables with RMS
      character*50 sd2d(nbsd2d)  /  'tsurf','ps','tauref'/
      ! 3D variables with RMS (pressure-wise and altitude-wise)
      character*50 sd3d(nbsd3d+1) / 'temp','u','v','rho','w','pressure'/
      
      real temp2d(dimlon,dimlat) ! temporary array

      if (typevar.eq.'mean') then

c MEAN 2d 
        do k=1,nbvar2d
          call get_2d(nf,k,varname2d(k),1)
          call get_2d(nf2,k,varname2d(k),2)
        enddo
c MEAN  only in low file! up data extrapolated 
        do k=1,nbvarlow-nbcom
          call get_3d(nf,k,varname3d(k),low,1)
          call get_3d(nf2,k,varname3d(k),low,2)
          call extrapol3d(typevar,k,up,pseudoalt)
        enddo
c MEAN up and low COMMON
        do k=nbvarlow-nbcom+1,nbvarlow
          call get_3d(nf,k,varname3d(k),low,1)
          call get_3d(nf2,k,varname3d(k),low,2)
          call get_3d(nf_up,k,varname3d(k),up,1)
          call get_3d(nf_up2,k,varname3d(k),up,2)
        enddo
c MEAN up - low data extrapolated 
        do k=nbvarlow+1,nbvar3d
          call get_3d(nf_up,k,varname3d(k),up,1)
          call get_3d(nf_up2,k,varname3d(k),up,2)
          call extrapol3d(typevar,k,low,pseudoalt)
        enddo

      elseif (typevar.eq.'hybr') then

        name='aps'
        ierr = NF_INQ_VARID(nf_up,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf_up,varid,aps_up)
        if (ierr.ne.nf_noerr) goto 9999

c if low atmosphere, need aps and bps         
        name='aps'
        ierr = NF_INQ_VARID(nf,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf,varid,aps_low)
        if (ierr.ne.nf_noerr) goto 9999
        name='bps'
        ierr = NF_INQ_VARID(nf,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf,varid,bps_low)
        if (ierr.ne.nf_noerr) goto 9999

c       hybrid coordinates: aps() and bps()
        do l=1,low
          aps(l)=aps_low(l)
          bps(l)=bps_low(l)
        enddo
        do l=low+1,dimlevs
          aps(l)=aps_up(l-low)
          bps(l)=0
        enddo

c       pseudoalt:
        do l=1,dimlevs
          pseudoalt(l)=-10*log(aps(l)/610.+bps(l))
        enddo
             
c       sigma Thermosphere : !!! not used any more
c        sigma(dimlevs+1) = sigma(dimlevs)*exp(-20./7.)  
         
      elseif (typevar.eq.'orog') then

        name='orography'
        ierr = NF_INQ_VARID(nf,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf,varid,temp2d)
        if (ierr.ne.nf_noerr) goto 9999
        do jloop=1,dimlat
          do iloop=1,dimlon ! taborog() is a common in constants_mcd.inc
           taborog(iloop,jloop,1)=temp2d(iloop,dimlat+1-jloop) 
          enddo
        enddo
        
        name='areoid'
        ierr = NF_INQ_VARID(nf,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf,varid,temp2d)
        if (ierr.ne.nf_noerr) goto 9999
        do jloop=1,dimlat
          do iloop=1,dimlon ! tabareo() is a common in constants_mcd.inc
           tabareo(iloop,jloop,1)=temp2d(iloop,dimlat+1-jloop) 
          enddo
        enddo

        name='z0'
        ierr = NF_INQ_VARID(nf,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf,varid,temp2d)
        if (ierr.ne.nf_noerr) goto 9999
        do jloop=1,dimlat
          do iloop=1,dimlon ! tabareo() is a common in constants_mcd.inc
           tabz0(iloop,jloop,1)=temp2d(iloop,dimlat+1-jloop) 
          enddo
        enddo

      elseif (typevar.eq.'grwp') then

        name='substd'
        ierr = NF_INQ_VARID(nf,name,varid) 
        if (ierr.ne.nf_noerr) goto 9999
        ierr = NF_GET_VAR_REAL(nf,varid,temp2d)
        if (ierr.ne.nf_noerr) goto 9999
        do jloop=1,dimlat
          do iloop=1,dimlon ! tabsubstd() is a common in constants_mcd.inc
           tabsubstd(iloop,jloop,1)=temp2d(iloop,dimlat+1-jloop) 
          enddo
        enddo

      elseif (typevar.eq.'rms') then

c RMS 2d 
        do k=1,nbsd2d
          call getsd_2d(nf,typevar,k,"rms"//sd2d(k),1)
          call getsd_2d(nf2,typevar,k,"rms"//sd2d(k),2)
        enddo
c RMS low - up data extrapolated
c      call getsd_3d(nf,typevar,1,nbsdlow-nbcomsd,"tsd"//sd3d(k),
c     $ low) 
c RMS up and low 
        do k=nbsdlow-nbcomsd+1,nbsdlow  
          call getsd_3d(nf,typevar,k,"rms"//sd3d(k),low,1)
          call getsd_3d(nf2,typevar,k,"rms"//sd3d(k),low,2)
          call getsd_3d(nf_up,typevar,k,"rms"//sd3d(k),up,1)
          call getsd_3d(nf_up2,typevar,k,"rms"//sd3d(k),up,2)
        enddo
c RMS up - low data extrapolated 
c        WARNING : uncomment these lines if to have a rms only in the
c        low atm (not the case in version 4)
c        do k=nbsdlow+1,nbsd3d
c         call getsd_3d(nf_up,typevar,k,"rms"//sd3d(k),up,1)
c         call getsd_3d(nf_up2,typevar,k,"rms"//sd3d(k),up,2)
c         call extrapol3d(typevar,k,low)
c        enddo

      elseif (typevar.eq.'arms') then
c Altitude-wise RMS
        do k=nbsdlow-nbcomsd+1,nbsdlow+1 ! +1 because we also want pressure
          call getsd_3d(nf,typevar,k,"arms"//sd3d(k),low,1)
          call getsd_3d(nf2,typevar,k,"arms"//sd3d(k),low,2)
          call getsd_3d(nf_up,typevar,k,"arms"//sd3d(k),up,1)
          call getsd_3d(nf_up2,typevar,k,"arms"//sd3d(k),up,2)
!          write(out,*)'k=',k,'vararms3d(2,48,2,k)=',vararms3d(2,48,2,k)
        enddo

      endif ! of if (typevar.eq.'mean') elseif ...

      return

c     Error handling
 9999 ier=16
      if (output_messages) then
        write(out,*)"loadvar_mcd Error : impossible to load ",
     &            name," from ",typevar, "value file"
      endif
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine loadeof(nf,rdeof,dust,dataset,ier)

c     load EOF arrays (which are commons in 'constants_mcd.inc')
      implicit none
     
      include "netcdf.inc"
      include "constants_mcd.inc"

c     inputs
      integer   nf !NetCDF file ID (file may have been previously opened)
      real      rdeof ! uniform random number in [0,1] for EOF perturb.
      integer   dust  ! dust scenario
      character(len=*) :: dataset ! path to datafiles
c     outputs
      integer      ier                 !error flag (0=OK, 1=NOK)

      character*16 name ! name of variable to read
      integer ierr      ! NetCDF routine (returned) status
      integer varid     ! NeCDF identifier
      character(len=45) :: datafile
      integer :: EOFyear ! EOF scenario
      integer, save :: EOFyear_prev=0 ! previous EOF scenario
      real :: rdeof2 ! random (between 0 and 1) generated from rdeof

      ! initialize ier to zero
      ier=0
      
      ! When using the dust=1-3 scenarios, there is now the possibility that we
      ! will swich between EOF files
      EOFyear=dust ! general case, the EOF scenario is the same as the dust scenario
      
      if ((dust.ge.1).and.(dust.le.3)) then
        ! build rdeof2 from rdeof, as the third (and following) digits
        ! of rdeof, and have it be between 0 and 1.
        rdeof2=100.*rdeof-floor(100.*rdeof)
!        write(*,*) "loadeof: rdeof=",rdeof," rdeof2=",rdeof2
        ! depending on rdeof2, we swich between EOF scenarios
        if (rdeof2.lt.1./7.) then
          EOFyear=2 ! MY24 scenario
        elseif (rdeof2.lt.2./7.) then
          EOFyear=4 ! MY26 scenario
        elseif (rdeof2.lt.3./7.) then
          EOFyear=5 ! MY27 scenario
        elseif (rdeof2.lt.4./7.) then
          EOFyear=7 ! MY29 scenario
        elseif (rdeof2.lt.5./7.) then
          EOFyear=8 ! MY30 scenario
        elseif (rdeof2.lt.6./7.) then
          EOFyear=9 ! MY31 scenario
        else
          EOFyear=10 ! MY32 scenario
        endif
      endif

      if (EOFyear.ne.EOFyear_prev) then
        ! we need to load some EOF data:
        ierr=NF_CLOSE(nf) ! close a possibly previously opened file
        ! which EOF datafile ?
        if ((dust.ge.1).and.(dust.le.3)) then ! clim scenarios
          if (EOFyear.eq.2) datafile="clim_aveEUV/MY24_all_var_eo.nc"
          if (EOFyear.eq.3) datafile="clim_aveEUV/MY25_all_var_eo.nc"
          if (EOFyear.eq.4) datafile="clim_aveEUV/MY26_all_var_eo.nc"
          if (EOFyear.eq.5) datafile="clim_aveEUV/MY27_all_var_eo.nc"
          if (EOFyear.eq.6) datafile="clim_aveEUV/MY28_all_var_eo.nc"
          if (EOFyear.eq.7) datafile="clim_aveEUV/MY29_all_var_eo.nc"
          if (EOFyear.eq.8) datafile="clim_aveEUV/MY30_all_var_eo.nc"
          if (EOFyear.eq.9) datafile="clim_aveEUV/MY31_all_var_eo.nc"
          if (EOFyear.eq.10) datafile="clim_aveEUV/MY32_all_var_eo.nc"
        elseif (dust.eq.4) then
          datafile="strm/strm_all_min_eo.nc"
        elseif (dust.eq.5) then
          datafile="strm/strm_all_ave_eo.nc"
        elseif (dust.eq.6) then
          datafile="strm/strm_all_max_eo.nc"
        elseif (dust.eq.7) then
          datafile="warm/warm_all_max_eo.nc"
        elseif (dust.eq.8) then
          datafile="cold/cold_all_min_eo.nc"
        else if ((dust.ge.24).and.(dust.le.32)) then
          datafile=""
          ! NB these files are stored in clim_aveEUV
          write(datafile,'(a14,i2,a14)')"clim_aveEUV/MY",
     &                                       dust,"_all_var_eo.nc"
        else
          ! no corresponding EOF file!
          goto 9999
        endif
        
        ! open EOF datafile:
        if (output_messages) then
          write(out,*) "Opening ",trim(dataset)//trim(datafile)
        endif
        ierr=NF_OPEN(trim(dataset)//trim(datafile),NF_NOWRITE,nf)
        if ((ierr.ne.nf_noerr).and.output_messages) then
          write(*,*) "Failed to open ",trim(dataset)//trim(datafile)
          write(*,*) NF_STRERROR(ierr)
!        else
!          write(*,*) "OK: opened ",trim(dataset)//trim(datafile)
        endif

        if (ierr.ne.nf_noerr) goto 9999
        
        ! load dataset
         ! Norm verctor for u
         name='normu'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeonormu)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! Norm vector for v
         name='normv'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeonormv)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! Norm vector for temperature
         name='normt'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeonormt)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! Norm vector for surface pressure
         name='normp'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeonormp)
         if (ierr.ne.nf_noerr) goto 9999
       
!         name='normr'
!         ierr = NF_INQ_VARID(nf,name,varid)
!         if (ierr.ne.nf_noerr) goto 9999
!         ierr = NF_GET_VAR_REAL(nf,varid,tabeonormr)
!         if (ierr.ne.nf_noerr) goto 9999
       
         ! zonal velocity
         name='u'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeou)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! meridional velocity
         name='v'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeov)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! atmospheric tempeerature
         name='temp' 
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeot)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! density
!         name='rho'
!         ierr = NF_INQ_VARID(nf,name,varid)
!         if (ierr.ne.nf_noerr) goto 9999
!         ierr = NF_GET_VAR_REAL(nf,varid,tabeorho)
!         if (ierr.ne.nf_noerr) goto 9999
       
         ! surface pressure
         name='ps'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabeops)
         if (ierr.ne.nf_noerr) goto 9999
       
         ! principal components
         name='pc'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabpc)
         if (ierr.ne.nf_noerr) goto 9999

         ! smoothed principal components
         name='pcsmth'
         ierr = NF_INQ_VARID(nf,name,varid)
         if (ierr.ne.nf_noerr) goto 9999
         ierr = NF_GET_VAR_REAL(nf,varid,tabpcsmth)
         if (ierr.ne.nf_noerr) goto 9999
       
!         name='pcvar'
!         ierr = NF_INQ_VARID(nf,name,varid)
!         if (ierr.ne.nf_noerr) goto 9999
!         ierr = NF_GET_VAR_REAL(nf,varid,tabpcvar)
!         if (ierr.ne.nf_noerr) goto 9999

        EOFyear_prev=EOFyear
      endif ! of if (EOFyear.ne.EOFyear_prev)

        return 
c     Error handling
 9999 ier=16
      if (output_messages) then
        write(out,*)"LOADEOF Error : failed to load ",
     &       trim(name)," from eof file ",trim(datafile)
      endif
      return

        end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mars_ltime(lon,jdate,ls,localtime,lmeantime)

c     Compute local true solar time at longitude=lon (in deg. east)
c     at julian date jdate, also given in Ls.

      implicit none

      real,intent(in) :: lon        !east longitude (in deg.)
      real*8,intent(in) :: jdate    !julian date
      real,intent(in) :: ls       !aerocentric longitude (in deg.)
      real,intent(out) :: localtime !local true solar time (in martian hours)
      real,intent(out) :: lmeantime ! local mean solar time (martian hours)

c     local variables
      double precision  LMT_0,E_T,LTST_0
      double precision Mars_LMT0
c marsday: number of seconds in a sol (matian day)
c      double precision     marsday
c      parameter (marsday=88775.245d0)
c  earthday: number of seconds in a day: 24*60*60
c      double precision     earthday
c      parameter (earthday=86400.d0)
c julian date of reference date: 01-01-1976 at 00:00:00
c      double precision     jdate_ref
c      parameter (jdate_ref=2442778.5d0) 
c LMT_ref: mean solar time at 0 deg. longitude at date jdate_ref
c      double precision    LMT_ref ! in martian hours
c      parameter (LMT_ref=16.1725d0) ! 01-01-1976 at 00:00:00
c orbital eccentricity
      double precision    eccentricity
      parameter (eccentricity=0.0934d0)
c Ls of perihelion (in degrees)
      double precision    Ls_peri
      parameter (Ls_peri=250.99d0)
c obliquity of equator to orbit (in deg.)
      double precision    obliquity
      parameter (obliquity=25.1919d0)
      double precision    pi,degtorad
      parameter (pi=3.14159265358979d0)
      parameter (degtorad=pi/180.0d0)

c compute LMT, local mean time, (in martian hours) at longitude 0
c      LMT_0=LMT_ref+24.0*(jdate-jdate_ref)*(earthday/marsday)
c      do while(LMT_0.le.0.d0) ! in case jdate<jdate_ref
c        LMT_0=LMT_0+24.0d0
c      enddo
c      LMT_0=dmod(LMT_0,24.0d0)
! a slightly more accurate computation of LMT_0
       LMT_0=Mars_LMT0(jdate)

c compute equation of time (see eq. 10.17, p.419 in
c "satellites orbits and missions", by M. Capderou),
c with a 24/(2*pi) factor to express it in martian hours
      E_T=(2.d0*eccentricity*dsin((ls-Ls_peri)*degtorad)
     &     -dtan(obliquity*degtorad/2.0)
     &      *dtan(obliquity*degtorad/2.0)*dsin(2.0*ls*degtorad))
     &       *24.0d0/(2.0*pi)


c compute true solar time at longitude 0
      LTST_0=LMT_0-E_T

c compute local true solar time at longitude lon
      localtime=real(LTST_0) + (lon/15.0)
      localtime=mod(24.0+localtime,24.0)

c compute local mean solar time at longitude lon
      lmeantime=real(LMT_0) + (lon/15.0)
      lmeantime=mod(24.0+lmeantime,24.0)
      
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function Mars_LMT0(djul)
      ! Compute local mean solar time on Mars at lon=0 in Mars hours (and
      ! decimal fraction thereof) for a given Julian date
      !
      ! Form M. Capderou, using references M. Allison & M. McEwen (Mars24)
      ! and computations from J. Tauber, Observatoire de Paris
      !
      ! Uses the fact that Julian date 2451549.5 (midnight at Greenwich)
      ! corresponds to almost midnight (23:59:39) LMT at Airy0 crater
      ! With 44796 added to avoid negative values (start of MSD)
      !
      ! Factor a1 from Allison & McEwen.
      ! Factors a2 and a3 added by M. Capderou
      ! for a better match with Mars24 outputs in the "modern era".
      ! Accurate to the second or less for years 1950 to 2050 (wrt Mars24).

      implicit none
      double precision, intent(in) :: djul ! julian date
      double precision, parameter :: dj2000=2451545.0d0
      double precision, parameter :: coeff=1.02749125d0 ! length of sol in days
      double precision :: a1,a2,a3,atot
      double precision :: x,y,xt,xanterre
      
      a1=7.2d-4
      xt=djul-(dj2000+4.5)
      x=xt/coeff
      x=x+44796.0d0-a1
      y=x-int(x)
      y=y*24.d0  ! Martian hours
      
      ! Small corrections to better fit to Mars24
      a2=41.0d0
      xanterre=(djul-dj2000)/365.25 ! Earth years wrt 01/01/2000
      if (xanterre.lt.0) then
        a3=28.5/45.*(0.-xanterre)+0.5
      else
        if (xanterre.lt.20.) then
          a3=3./20.*(20.-xanterre)-2.5
        else
          a3=-1
        endif
      endif
      atot=a2-a3
      atot=atot/3600.d0 ! Martian hours
      
      Mars_LMT0=y+atot
      
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mars_ptime(lon,lt,ut)

c     Convert local true solar time, lt (0..24) at east longitude (in deg.) 
c     into "universal time"  ut (0..24), the local true solar time at lon=0

      implicit none
      
c     inputs
      real lon    !longitude east (in degrees)
      real lt     !local time (in martian hours)
      
c     output
      real ut     !universal time  (local true solar time at lon=0)


      ut=mod(24.0+lt-lon/15.0,24.0)

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mcd_time(ut,itime,w)
      
c     From universal time ut (0..24 hrs) , find the 
c     2 nearest timestep in the database (1 - 12) 
c     and give the linear interpolation weight of itime(1)
c 
c     Note that it is midnight universal time (and thus at longitude=0.0)
c     at itime = 12

      implicit none

c     input
      real ut       !universal time (0. to 24. hrs) =local time at lon=0

c     outputs
      integer  itime(2)     ! 2 nearest timestep in database 
      real  w            ! linear interpolation weight of itime(1)
                         ! ie: w=1 if ut=itime(1) and w=0 if ut=itime(2)

      itime(1) = int(ut/2.)
      if (itime(1).eq.0) itime(1) = 12
      itime(2) = itime(1) +1
      if (itime(2).ge.13.) itime(2) = itime(2) -12

      w = 1 - ut/2. + int(ut/2.)
      
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine opend(unet,unetm,unetsd,unetm_up,unetsd_up,
     &            num,dust,dset,ier)

c     Open the appropriate NETCDF file corresponding to Ls, dust scenario  
c     and data set.
c     
c     Sept. 2012. Scenario files are now in individual subdirectories
c     April 2008. Added search for 'datafile.gz' in same directory in case
c                'datafile' is not found. EM

      implicit none

      include "netcdf.inc"
      include "constants_mcd.inc"

c     inputs
c NETCDF channels number      
      integer       unet       !for mountain fields
      integer       unetm      !for mean fields
      integer       unetsd     !for STD and RMS deviation fields
      integer       unetm_up   !for mean fields in thermosphere 
      integer       unetsd_up  !for STD and RMS fields in thermosphere 

      integer       num        !season num
      integer       dust       !Dust senario
      character*(*) dset       !Dataset

c     output
      integer       ier        !Error flag (0=OK, 1=NOK)

c     local variables
      integer       ierr     ! return code from NetCDF functions
      character*258 datfile  ! directory where data is to be found
      character*4   scen     ! dust scenario name ('clim','cold', 'warm' or 'strm')
      character*12  scenlow  ! dust scenario directory
                             ! for lower atm ('cold','warm',etc.)
      character*12  scenup   ! dust scenario directory
                             ! for upper atm ('cold','warm',etc.)
      character*3   solar    ! EUV scenario ('min' or 'ave' or 'max')
      character*2   saison   ! month # (from '01' to '12')
      character*2   typ      ! type of datafile ('me' or 'sd' or 'eo')
      character*258 gzfile   ! 'datfile'//'.gz'
      
      ier=0
c     close files 
      ierr=NF_CLOSE(unetm)
      ierr=NF_CLOSE(unetsd)
      ierr=NF_CLOSE(unetm_up)
      ierr=NF_CLOSE(unetsd_up)
      ierr=NF_CLOSE(unet)

c     dust scenario
      if (dust.eq.1.or.dust.eq.2.or.dust.eq.3) then
         scen='clim'
         scenlow='clim_aveEUV'
         if (dust.eq.1) then
           scenup='clim_aveEUV'
           solar='ave' 
         endif
         if (dust.eq.2) then
           scenup='clim_minEUV'
           solar='min' 
         endif
         if (dust.eq.3) then
           scenup='clim_maxEUV'
           solar='max'
         endif 
      elseif ((dust.eq.4).or.(dust.eq.5).or.(dust.eq.6)) then
         scen='strm'
         scenlow='strm'
         scenup='strm'
         if (dust.eq.4) solar='min' 
         if (dust.eq.5) solar='ave' 
         if (dust.eq.6) solar='max' 
      elseif (dust.eq.7) then
         scen='warm'
         scenlow='warm'
         scenup='warm'
         solar='max'
      elseif (dust.eq.8) then
         scen='cold'
         scenlow='cold'
         scenup='cold'
         solar='min'
      elseif ((dust.ge.24).and.(dust.le.32)) then
         write(scen,'(a2,i2)') "MY",dust
         write(scenlow,'(a2,i2)') "MY",dust
         write(scenup,'(a2,i2)') "MY",dust
         solar="var"
      else
         if (output_messages) then
           write(out,*) 'pb in opend with dust= ', dust
           write(out,*) '   no such scenario!'
         endif
         stop
      endif

c     season number
      write(saison,'(i2.2,1x)') num 

c     mean file
      typ='me'
      datfile=trim(dset)//trim(scenlow)//'/'//trim(scen)//
     &        '_'//saison//'_'//typ//'.nc'
      if (output_messages) then
        write(out,*) "Opening ",trim(datfile)
      endif
      ierr=NF_OPEN(datfile,NF_NOWRITE,unetm)
      if (ierr.ne.NF_NOERR) then ! failed to open file
        ! check for a '.gz' version of the datafile
        gzfile=trim(dset)//trim(scenlow)//'/'//trim(scen)//
     &         '_'//saison//'_'//typ//'.nc.gz'
        open(40,file=gzfile,iostat=ierr,status="old")
        if (ierr.ne.NF_NOERR) then ! failed to open file.gz
          goto 9999
        else ! there is a file.gz around
          goto 8888
        endif
      endif

c     standard deviation file
      typ='sd'
      datfile=trim(dset)//trim(scenlow)//'/'//trim(scen)//
     &        '_'//saison//'_'//typ//'.nc'
      if (output_messages) then
        write(out,*) "Opening ",trim(datfile)
      endif
      ierr=NF_OPEN(datfile,NF_NOWRITE,unetsd)
      if (ierr.ne.NF_NOERR) then ! failed to open file
        ! check for a '.gz' version of the datafile
        gzfile=trim(dset)//trim(scenlow)//'/'//trim(scen)//
     &         '_'//saison//'_'//typ//'.nc.gz'
        open(40,file=gzfile,iostat=ierr,status="old")
        if (ierr.ne.NF_NOERR) then ! failed to open file.gz
          goto 9999
        else ! there is a file.gz around
          goto 8888
        endif
      endif
      
c     thermo Mean file
      typ='me'
      datfile=trim(dset)//trim(scenup)//'/'//trim(scen)//
     &        '_'//saison//'_thermo_'//solar//'_'//typ//'.nc'
      if (output_messages) then
        write(out,*) "Opening ",trim(datfile)
      endif
      ierr=NF_OPEN(datfile,NF_NOWRITE,unetm_up)
      if (ierr.ne.NF_NOERR) then ! failed to open file
        ! check for a '.gz' version of the datafile
        gzfile=trim(dset)//trim(scenup)//'/'//trim(scen)//
     &         '_'//saison//'_thermo_'//solar//'_'//typ//'.nc.gz'
        open(40,file=gzfile,iostat=ierr,status="old")
        if (ierr.ne.NF_NOERR) then ! failed to open file.gz
          goto 9999
        else ! there is a file.gz around
          goto 8888
        endif
      endif

c     thermo standard deviation file
      typ='sd'
      datfile=trim(dset)//trim(scenup)//'/'//trim(scen)//
     &        '_'//saison//'_thermo_'//solar//'_'//typ//'.nc'
      if (output_messages) then
        write(out,*) "Opening ",trim(datfile)
      endif
      ierr=NF_OPEN(datfile,NF_NOWRITE,unetsd_up)
      if (ierr.ne.NF_NOERR) then ! failed to open file
        ! check for a '.gz' version of the datafile
        gzfile=trim(dset)//trim(scenup)//'/'//trim(scen)//
     &         '_'//saison//'_thermo_'//solar//'_'//typ//'.nc.gz'
        open(40,file=gzfile,iostat=ierr,status="old")
        if (ierr.ne.NF_NOERR) then ! failed to open file.gz
          goto 9999
        else ! there is a file.gz around
          goto 8888
        endif
      endif
      
c     mountain file
      datfile=dset//'mountain.nc'
      ierr=NF_OPEN(datfile,NF_NOWRITE,unet)
      if (ierr.ne.NF_NOERR) goto 9999

      return

c     Error handling

 8888 ier=21    ! could not open file but found file.gz instead
      if (output_messages) then
        write(out,*)"Error in opend: cannot open file ",trim(datfile)
        write(out,*)"  but found file ",trim(gzfile)
        write(out,*)"  which should be uncompressed using third-party"//
     &              " software (e.g. gunzip on Unix or Winzip on "//
     &              " Windows) to produce sought file."
      endif
      return
      
 9999 ier=15    ! could not open file (and no file.gz around)
      if (output_messages) then
        write(out,*)"Error in opend: cannot open file ",trim(datfile)
      endif
      return

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine orbit(date,ls,marsau,outmodelday)

c     Given Julian date,compute Ls, sun-Mars distance
c     in AU and model day

      implicit none

c     input
      real*8     date   !Julian date

c     outputs
      real       ls     !Ls: Aerocentric longitude (in deg.)
      real       marsau !Sun-Mars distance (in AU)
      real       outmodelday

c     local variables
      double precision     modelday
c marsday: number of seconds in a sol (matian day)
      double precision     marsday
      parameter (marsday=88775.245d0)
c  earthday: number of seconds in a day: 24*60*60
      double precision     earthday
      parameter (earthday=86400.d0)
      double precision     marsyear
      parameter (marsyear=668.5921d0) ! number of sols in a martian year
c julian date for 19-12-1975 at 4:00:00, at which Ls=0.0
      double precision jdate_ref
      parameter (jdate_ref=2442765.667d0)
      double precision pi,degtorad
      parameter (pi=3.14159265358979d0)
      parameter (degtorad=pi/180.0d0) 

      real       dummy

c     convert Julian day to model day
      modelday=(date-jdate_ref)*earthday/marsday
       do while(modelday.le.0.d0) !in case date<jdate_ref
         modelday=modelday+marsyear
       enddo
      modelday=dmod(modelday,marsyear)
      dummy=real(modelday)
      outmodelday=real(modelday)

      call sol2ls(dummy,ls)

c compute sun-mars distance (in AU)
      call sunmarsdistance(ls,marsau)

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine sunmarsdistance(ls,marsau)
      ! compute Sun-Mars distance (in AU)
      implicit none
      real,intent(in) :: ls ! solar longitude (deg)
      real,intent(out) :: marsau ! Sun-Mars distance (AU)
      
      ! semi-major axis of orbit (in AU)
      double precision,parameter :: sma=1.52368d0
      ! orbital eccentricity
      double precision,parameter :: eccentricity=0.09340d0
      ! Ls of perihelion
      double precision,parameter :: Ls_peri=250.99
      double precision,parameter :: pi=3.14159265358979d0
      double precision,parameter :: degtorad=pi/180.0d0

      marsau=real(sma*(1.0d0-eccentricity*eccentricity)/
     &       (1.0d0+eccentricity*dcos((dble(ls)-Ls_peri)*degtorad)))

      end
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine solarzenithangle(lat,ls,localtime,solzenang)
      ! compute solar zenith angle
      implicit none
      
      real,intent(in) :: lat ! latitude (deg)
      real,intent(in) :: ls ! solar longitude (deg)
      real,intent(in) :: localtime ! local true solar time (hours)
      real,intent(out) :: solzenang ! (deg) solar zenith angle
      
      ! local variables:
      double precision :: declin ! declination of the Sun (rad)
      double precision,parameter :: obliquity=25.1919d0
      double precision,parameter :: pi=3.14159265358979d0
      double precision,parameter :: degtorad=pi/180.0d0
      double precision,parameter :: radtodeg=180.0d0/pi
      
      double precision :: mu0 ! cos of solar zenith angle
      
      ! Compute Sun's declination
      declin=asin(sin(ls*degtorad)*sin(obliquity*degtorad))
      
      ! Compute mu0
      mu0=sin(lat*degtorad)*sin(declin)+
     &    cos(lat*degtorad)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
      
      solzenang=ACOS(mu0)*radtodeg
      
      end
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine season(ls,numsaison)

c     compute the season number for a given areocentric longitude

      implicit none

c     input
      real    ls

c     output 
      integer numsaison

      if ((ls.ge.0.0).and.(ls.le.30)) then
         numsaison=1
      endif
      if ((ls.ge.30.0).and.(ls.le.60)) then
         numsaison=2
      endif
      if ((ls.ge.60.0).and.(ls.le.90)) then
         numsaison=3
      endif
      if ((ls.ge.90.0).and.(ls.le.120)) then
         numsaison=4
      endif
      if ((ls.ge.120.0).and.(ls.le.150)) then
         numsaison=5
      endif
      if ((ls.ge.150.0).and.(ls.le.180)) then
         numsaison=6
      endif
      if ((ls.ge.180.0).and.(ls.le.210)) then
         numsaison=7
      endif
      if ((ls.ge.210.0).and.(ls.le.240)) then
         numsaison=8
      endif
      if ((ls.ge.240.0).and.(ls.le.270)) then
         numsaison=9
      endif
      if ((ls.ge.270.0).and.(ls.le.300)) then
         numsaison=10
      endif
      if ((ls.ge.300.0).and.(ls.le.330)) then
         numsaison=11
      endif
      if ((ls.ge.330.0).and.(ls.le.360)) then
         numsaison=12
      endif
      
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine season2(ls,numsaison,nums2,wl,wh,dust)

c     compute the 2 season numbers(nums,nums2) from areocentric longitude
c     for when itimint = 1, ie for seasonal interpolation
c     wl= lower season weighting, wh=higher season weighting

      implicit none

c     inputs
      real    ls
      integer dust

c     outputs 
      integer numsaison,nums2
      real    wl,wh

c     local variable
      real    lsc      

c     compute lower season
      numsaison=int((ls/30)+0.5)
      if (numsaison.eq.0) numsaison=12
c     compute higher season
      nums2=numsaison+1   
      if (nums2.eq.13) nums2=1
c     compute centre of season to find weightings
      lsc=((numsaison-1)*30)+15
      if (ls.lt.15) lsc=-15
c     compute weightings
      wh=(ls-lsc)/30.
      wl=1-wh

c     Dust storm scenario Special case for Ls near 180 and Ls near 360
      if ((dust.eq.4).or.(dust.eq.5).or.(dust.eq.6))then 
          if(numsaison.eq.6) numsaison=7
          if(nums2.eq.1) nums2=12
      endif 

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine sol2ls(sol,ls)

c  convert a given martian day number (sol)
c  into corresponding solar longitude, Ls (in degr.),
c  where sol=0=Ls=0 is the
c  northern hemisphere spring equinox.

      implicit none

c     input 
      real    sol
c     output 
      real    ls

c     local variables
      double precision    year_day,peri_day,timeperi,e_elips
      double precision    pi,radtodeg
c number of martian days (sols) in a martian year
      parameter (year_day=668.6d0)
c perihelion date (in sols)
      parameter (peri_day=485.35d0)
c orbital eccentricity
      parameter (e_elips=0.09340d0)
      parameter (pi=3.14159265358979d0)
c  radtodeg: 180/pi
      parameter (radtodeg=57.2957795130823d0)
c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
      parameter (timeperi=1.90258341759902d0)

      double precision    zanom,xref,zx0,zdx,zteta,zz
c  xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly        
      integer iter

      zz=(sol-peri_day)/year_day
      zanom=2.*pi*(zz-nint(zz))
      xref=dabs(zanom)

c  The equation zx0 - e * sin (zx0) = xref, solved by Newton
      zx0=xref+e_elips*dsin(xref)
      do 110 iter=1,10
         zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
         if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough
           goto 120
         endif
         zx0=zx0+zdx
  110 continue
  120 continue
      zx0=zx0+zdx
      if(zanom.lt.0.) zx0=-zx0

c compute true anomaly zteta, now that eccentric anomaly zx0 is known
      zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.))

c compute Ls
      ls=real(zteta-timeperi)
      if(ls.lt.0.) ls=ls+2.*real(pi)
      if(ls.gt.2.*pi) ls=ls-2.*real(pi)
c convert Ls in deg.
      ls=real(radtodeg)*ls

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real function ls2sol(ls)

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

      implicit none

c  Arguments:
      real ls

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

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

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

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine var2d(a,lon,lat,utime,name,ier,itimint,wl,wh,ps_psgcm)

c     Retrieve the value of 2-d variable 'name' at longitude=lon, latitude=lat 
c     and universal time=utime.
c     - Bilinear interpolation is used to interpolate variable to user
c       specified longitude and latitude. For surface pressure (name='ps'),
c       bilinear interpolation of log(pressure) is used.
c     - Linear interpolation is used to interpolate variable to user
c       specified time (of day) 'utime'.
c     - If (itimint>=1) then linear seasonal interpolation (with wl=weight 
c       of lower season and wh=weight of higher season) is then used.

      implicit none
      
      include "constants_mcd.inc"
      
c     inputs
      real         lon     !longitude (east) of point
      real         lat     !latitude of point
      real         utime   !Universal time (0. to 24. hrs) =local time at lon=0
      character*16 name    !Name of variable
      integer      itimint !seasonal interpolation flag
                           !  (==1 if seasonal interpolation is required)
      real         wl,wh   !seasonal interpolation weights
      real         ps_psgcm ! High res to GCM surface pressure ratio 

c     outputs
      integer      ier     !error flag (0=OK, 1=NOK)
      real         a       !value of variable

c     local variables
      real         at(2)    !seasonal interpolation variable  
      integer      i, j
      integer      k,sd,rms ! "flags" associated to variable 'name'
      integer      itime(2) ! nearest MCD time indexes
      integer      iut
      integer      dlon(4) ! nearest grid point longitude indexes
      integer      dlat(4) ! nearest grid point latitude indexes
      real         y(2,4)    ! to store nearest grid point values
      real         x(2,2)  ! bilinearly interpolated values ("lower" season)
      real         t,u   ! longitude and latitude-wise interpolation weights
      real         w     ! time-wise interpolation weight
      integer      itype,it
      
      ! rmsname(:) same ordering as  sd2d(:) in call_mcd with "rms" prefix
      character*12 rmsname(nbsd2d) /'rmstsurf','rmsps', 'rmstauref'/
      
      ier = 0

!      write(*,*)"var2d: name=",trim(name)

c     find nearest 4 grid points
      call grid4(lon,lat,dlon,dlat,t,u)
c     dlon() and dlat() now contain longitude and latitude indexes
c     and t and u (interpolation weights) are consequently set
!      write(*,*)"var2d: dlon=",dlon," dlat=",dlat
!      write(*,*)"var2d: t=",t," u=",u

c     find nearest 2 timestep :
      call mcd_time(utime,itime,w)
!      write(*,*)"var2d: itime=",itime," w=",w

c     flags sd and rms initialized to 999
      k=0
      sd=999 
      rms=999

c     Associate numbers with variable name to read loaded variable
c     (note these should be compatible with what is stated in 'loadvar_mcd')
c     ------------------------------------------------------------
      do itype=0,1
       if(itype.eq.0) then
        do k=1,nbvar2d
         if(name.eq.varname2d(k)) then
          sd = 0   
          goto 1111   
         endif     
        enddo    
       else
        do k=1,nbsd2d
         if(name.eq.rmsname(k)) then
          rms = 1   
          goto 1111   
         endif    
        enddo   
       endif   
      enddo      
      
      if((name.ne.'orography').and.(name.ne.'z0').and.
     &   (name.ne.'areoid').and.(name.ne.'substd')) then
c        CASE of an unexpected name (orography, areoid, z0, and substd
c             are treated below)
       if(output_messages) then
         write(out,*) 'problem using subroutine var2d : the name ',name
         write(out,*) 'is not recognized'
       endif
       stop
      endif       
     
1111  continue

      do it=1,2
 
c it=1 Retrieving variable for "lower" season
c it=2 Retrieving variable for "higher" season (if seasonal interpolation)
c     --------------------------------------
c Note: Arrays taborog(), tabsubstd(), var_2d(),....
c       are known from "constants_mcd.inc" and have been previously filled
c     loop on the 2 nearest timestep :
      do j = 1 , 2
         iut = itime(j)

c     retrieve the four values at the nearest grid points from the array
          if (name.eq.'orography') then
            do i=1,4
               y(it,i)=taborog(dlon(i),dlat(i),1)
            enddo
          elseif (name.eq.'areoid') then
            do i=1,4
               y(it,i)=tabareo(dlon(i),dlat(i),1)
            enddo
          elseif (name.eq.'substd') then
            do i=1,4
               y(it,i)=tabsubstd(dlon(i),dlat(i),1)
            enddo
          elseif (name.eq.'z0') then
            do i=1,4
               y(it,i)=tabz0(dlon(i),dlat(i),1)
            enddo
          elseif (sd.eq.0) then
           if(it.eq.1) then  
            do i=1,4
              y(it,i)=var_2d(dlon(i),dlat(i),iut,k)
            enddo
           else
            do i=1,4
              y(it,i)=var_2d2(dlon(i),dlat(i),iut,k)
            enddo   
           endif 
          elseif (rms.eq.1) then
            if(it.eq.1) then  
             do i=1,4
               y(it,i)=varrms2d(dlon(i),dlat(i),k)
             enddo
            else
              do i=1,4
               y(it,i)=varrms2d2(dlon(i),dlat(i),k)
             enddo  
            endif    
          endif
          
          if ((name.eq.'ps').or.
     &        (name.eq.'rmsps')) then
c         bilinear interpolation of log(surface pressure)
            x(it,j)=(1.-t)*(1.-u)*log(y(it,1))+t*(1.-u)*log(y(it,2))
     &          +t*u*log(y(it,3))+(1.-t)*u*log(y(it,4))
            x(it,j)=exp(x(it,j))
          else
c         bilinear interpolation in space
            x(it,j)=(1.-t)*(1.-u)*y(it,1)+t*(1.-u)*y(it,2)
     &          +t*u*y(it,3)+(1.-t)*u*y(it,4)
          endif
      end do ! of do j=1,2 loop

c     linear interpolation in time:
      at(it) = w*x(it,1) + (1-w)*x(it,2)
            
      if (itimint.lt.1) exit
      
      enddo ! of do it=1,2
      
      if (itimint.ge.1) then
c     Season interpolation between the two seasons
c     --------------------------------------------        
       a = wl*at(1) + wh*at(2)
      else
       a = at(1)
      endif 

c     Multiply (ie: rescale) by ps_psgcm for some variables
      ! NB: 'tauref' is dod (dust optical depth) at 610 Pa.
      !    so we don't need to rescale it
!      if ((name.eq.'tauref').or.(name.eq.'rmstauref').or.
!     &    (name.eq.'rmsps').or.
!     &    (name.eq.'col_h2ovapor')) then
      if (name.eq.'rmsps') then
        a=a*ps_psgcm
      endif

c     error handling (unused yet)
! 9999 a=0.
!      return
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine interpol(k,up_low,vtype,pseudoalt)

c     calculate the vertical interpolation between low and up atmosphere    
c     when data are not in both files (low and thermo)    
      
      implicit none
      
      include "constants_mcd.inc"

c     inputs      
      integer    k,up_low
      character*4 vtype
      real       pseudoalt
      dimension  pseudoalt(dimlevs)

c     local variables
      integer     startk,endk
      integer     iloop,kloop,jloop,it

      if (up_low.eq.up) then 
      startk=low+1
      endk=low+3
      else
      startk=low-2        
      endk=low
      endif

c     3D extrapolation for the first 3 layers of the thermosphere          
      if (vtype.eq.'mean') then
       do it=1,dimuti
        if (k.ne.10) then
         do kloop=startk,endk
           do jloop=1,dimlat
            do iloop=1,dimlon
       var_3d(iloop,jloop,kloop,it,k)=var_3d(iloop,jloop,startk-1,it,k)+
     $ (pseudoalt(kloop)-pseudoalt(startk-1))*
     $ (var_3d(iloop,jloop,endk+1,it,k)-
     $ var_3d(iloop,jloop,startk-1,it,k))
     $ /(pseudoalt(endk+1)-pseudoalt(startk-1))

       var_3d2(iloop,jloop,kloop,it,k)=
     $ var_3d2(iloop,jloop,startk-1,it,k)+
     $ (pseudoalt(kloop)-pseudoalt(startk-1))*
     $ (var_3d2(iloop,jloop,endk+1,it,k)-
     $ var_3d2(iloop,jloop,startk-1,it,k))
     $ /(pseudoalt(endk+1)-pseudoalt(startk-1))
            enddo
           enddo
          enddo
      elseif(k.eq.10) then   ! "o" case
         do kloop=2,30      ! interpolation for the all low atmosphere
           do jloop=1,dimlat
            do iloop=1,dimlon
        var_3d(iloop,jloop,kloop,it,k)=
     $ exp(
     $ log(var_3d(iloop,jloop,1,it,k))+
     $ (pseudoalt(kloop)-pseudoalt(1))*
     $ (log(max(var_3d(iloop,jloop,31,it,k),1.E-30))-
     $ log(var_3d(iloop,jloop,1,it,k)))
     $ /(pseudoalt(31)-pseudoalt(1))
     $ )

       var_3d2(iloop,jloop,kloop,it,k)=
     $ exp(
     $ log(var_3d2(iloop,jloop,1,it,k))+
     $ (pseudoalt(kloop)-pseudoalt(1))*
     $ (log(max(var_3d2(iloop,jloop,31,it,k),1.E-30))-
     $ log(var_3d2(iloop,jloop,1,it,k)))
     $ /(pseudoalt(31)-pseudoalt(1))
     $ )
            
            enddo !iloop
           enddo !jloop
          enddo !kloop
         endif! "o" case 
        enddo   
      elseif (vtype.eq.'rms') then
c    2D extrapolation for the first 3 layers of the thermosphere          
         do kloop=startk,endk
           do jloop=1,dimlat
            do iloop=1,dimlon
       varrms3d(iloop,jloop,kloop,k)=
     $ varrms3d(iloop,jloop,startk-1,k)+
     $ (pseudoalt(kloop)-pseudoalt(startk-1))*
     $ (varrms3d(iloop,jloop,endk+1,k)-
     $ varrms3d(iloop,jloop,startk-1,k))
     $ /(pseudoalt(endk+1)-pseudoalt(startk-1))

       varrms3d2(iloop,jloop,kloop,k)=
     $ varrms3d2(iloop,jloop,startk-1,k)+
     $ (pseudoalt(kloop)-pseudoalt(startk-1))*
     $ (varrms3d2(iloop,jloop,endk+1,k)-
     $ varrms3d2(iloop,jloop,startk-1,k))
     $ /(pseudoalt(endk+1)-pseudoalt(startk-1))
            enddo
           enddo
          enddo
       endif !vtype
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine extrapol3d(vtype,k,up_low,pseudoalt)

c  extrapolate 3d data when not in the initial netcdf files 
c  call subroutine interpol 

      implicit none
      
      include "constants_mcd.inc"
      include "netcdf.inc"

c     inputs      
      integer     k,up_low
      character*4 vtype
      real        pseudoalt
      dimension   pseudoalt(dimlevs)

c     local variables      
      integer     iloop,kloop,jloop,it
      real        var_ini(nbvarup-nbcom)

c     initial values in low atmosphere ! not valid for MCDv5 !!!
!      var_ini(1)=0.2E-11!o : only for surface
!      var_ini(2)=95.32E-2 !co2
!      var_ini(3)=8.E-4 !co
!      var_ini(4)=0.027!n2

      if (up_low.eq.up) then !extrapolation in the thermosphere 

       if (vtype.eq.'mean') then
        do it=1,dimuti
         do kloop=low+4,dimlevs
           do jloop=1,dimlat
            do iloop=1,dimlon
                 var_3d(iloop,jloop,kloop,it,k)=0.
                 var_3d2(iloop,jloop,kloop,it,k)=0.
            enddo
           enddo
          enddo
         enddo

       elseif (vtype.eq.'stdv'.or.vtype.eq.'rms')    then 
         do kloop=low+4,dimlevs
          do jloop=1,dimlat
           do iloop=1,dimlon
             if (vtype.eq.'rms') then
                varrms3d(iloop,jloop,kloop,k)=0.
                varrms3d2(iloop,jloop,kloop,k)=0.
!             else if (vtype.eq.'stdv') then
!                varsd3d(iloop,jloop,kloop,k)=0.
!                varsd3d2(iloop,jloop,kloop,k)=0.
             endif
            enddo
           enddo
          enddo

       endif    

c extrapolation for the first 3 layers of the thermosphere          
c *********************************************************        
      call  interpol(k,up,vtype,pseudoalt)

      else  ! extrapolation in the lower atmosphere   

       if (vtype.eq.'mean') then
         do it=1,dimuti
          do kloop=1,low-3
           do jloop=1,dimlat
            do iloop=1,dimlon
               var_3d(iloop,jloop,kloop,it,k)=var_ini(k-nbvarlow)
               var_3d2(iloop,jloop,kloop,it,k)=var_ini(k-nbvarlow)
             enddo
            enddo
           enddo
           enddo

       elseif (vtype.eq.'stdv'.or.vtype.eq.'rms')    then 
         do kloop=1,low-3
           do jloop=1,dimlat
            do iloop=1,dimlon
              if (vtype.eq.'rms') then
               varrms3d(iloop,jloop,kloop,k)=0.
               varrms3d2(iloop,jloop,kloop,k)=0.
!              else if (vtype.eq.'stdv') then
!               varsd3d(iloop,jloop,kloop,k)=0.
!               varsd3d2(iloop,jloop,kloop,k)=0.
              endif
            enddo
           enddo
          enddo
        endif    

c extrapolation for the last 3 layers of the low atmosphere          
c *********************************************************        
       call  interpol(k,low,vtype,pseudoalt)
        
       endif! up_low          

      return
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine get_2d(nf,k,varname,order)

c  get_2d reads the 2d mean variables in low netcdf files.
c  Fills var_2d() and var_2d2() (defined in constants_mcd.inc)  

      implicit none

      include "constants_mcd.inc"
      include "netcdf.inc"

c     inputs      
      integer      nf,k,order
      character*50 varname

c     local variables     
      integer      ierr,iloop,jloop,it
      integer      var2didin(nbvar2d)
      real         temp2d(dimlon,dimlat,dimuti)

c     write(out,*),'varname(k) k',k,varname
c     NETCDF reading 2d mean variable      
         ierr=NF_INQ_VARID(nf,varname,var2didin(k))
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GET_2D Error:"
             write(out,*) NF_STRERROR(ierr), varname
           endif
           stop
         endif
         ierr = NF_GET_VAR_REAL(nf,var2didin(k),temp2d)
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GET_2D Error:"
             write(out,*) NF_STRERROR(ierr)
           endif
           stop
         endif

         do it=1,dimuti
           do jloop=1,dimlat
            do iloop=1,dimlon
                  if (order.eq.1) then
                  var_2d(iloop,jloop,it,k)=
     &            temp2d(iloop,dimlat+1-jloop,it) 
                  elseif (order.eq.2) then 
                  var_2d2(iloop,jloop,it,k)=
     &            temp2d(iloop,dimlat+1-jloop,it) 
                  endif
            enddo
           enddo
          enddo

      return
      end 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine get_3d(nf,k,varname,up_low,order)

c  getsd_3D reads the 3D mean variables in up or low netcdf files.
c  Fills var_3d() and var_3d2() (defined in constants_mcd.inc)  

      implicit none

      include "constants_mcd.inc"
      include "netcdf.inc"

c     inputs      
      integer      ierr,nf,k,up_low,order
      character*50 varname

c     local variables      
      integer      iloop,kloop,jloop,it
      integer      var3didin(nbvar3d)
      real         temp3dlow(dimlon,dimlat,low,dimuti)
      real         temp3dup(dimlon,dimlat,up,dimuti)

c     write(out,*),'varname(k) k',k,varname
      ierr=NF_INQ_VARID(nf,varname,var3didin(k))
c         ierr=NF_INQ_VARID(nf,varname,var3didin2(k))
      if(ierr.ne.NF_NOERR) then
        if (output_messages) then
          write(out,*) "GET_3D Error:"
          write(out,*) NF_STRERROR(ierr),varname
        endif
        stop 
      endif
      if (up_low.eq.up) then
         ierr = NF_GET_VAR_REAL(nf,var3didin(k),temp3dup)
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GET_3D Error:"
             write(out,*) NF_STRERROR(ierr)
           endif
           stop
         endif

         do it=1,dimuti
          do kloop=low+1,dimlevs
           do jloop=1,dimlat
            do iloop=1,dimlon
             if (order.eq.1) then
                  var_3d(iloop,jloop,kloop,it,k)=
     &            temp3dup(iloop,dimlat+1-jloop,kloop-low,it) 
             elseif(order.eq.2) then      
                  var_3d2(iloop,jloop,kloop,it,k)=
     &            temp3dup(iloop,dimlat+1-jloop,kloop-low,it) 
             endif
            enddo
           enddo
          enddo
         enddo 

      else    

        ierr = NF_GET_VAR_REAL(nf,var3didin(k),temp3dlow)
        if(ierr.ne.NF_NOERR) then
          if (output_messages) then
            write(out,*) "GET_3D Error:"
            write(out,*) NF_STRERROR(ierr)
          endif
          stop
        endif

        do it=1,dimuti
         do kloop=1,low
          do jloop=1,dimlat
           do iloop=1,dimlon
            if (order.eq.1) then
                  var_3d(iloop,jloop,kloop,it,k)=
     &            temp3dlow(iloop,dimlat+1-jloop,kloop,it) 
            elseif(order.eq.2) then      
                  var_3d2(iloop,jloop,kloop,it,k)=
     &            temp3dlow(iloop,dimlat+1-jloop,kloop,it) 
            endif     
             enddo
            enddo
          enddo
         enddo

      endif ! of if (up_low.eq.up)

      return
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine getsd_2d(nf,vtype,k,varname,order)

c  getsd_2D reads the 2D RMS variables from netcdf files and
c fills corresponding varrms2d(k) and varrms2d2(k) arrays (which are defined
c as common in constants_mcd.inc)

      implicit none

      include "constants_mcd.inc"
      include "netcdf.inc"

c     inputs      
      integer      nf      ! Input NetCDF file ID
      character*4  vtype
      integer      k
      character*50 varname
      integer      order

c     local variables      
      integer      ierr,iloop,jloop
      integer      sd2didin(nbsd2d)    ! to store ID # of variable
      real         temp2d(dimlon,dimlat,dimuti) ! to temporarely store variable

c      write(out,*),'varname(k) k',k,varname

c  Read variable from NetCDF file      
         ierr=NF_INQ_VARID(nf,varname,sd2didin(k))
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GETSD_2D Error:"
             write(out,*) NF_STRERROR(ierr),varname
           endif
           stop 
         endif
         ierr = NF_GET_VAR_REAL(nf,sd2didin(k),temp2d)
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GETSD_2D Error:"
             write(out,*) NF_STRERROR(ierr) , 'stop in getsd3d'
           endif
           stop 
         endif

c Copy variable
           do jloop=1,dimlat
            do iloop=1,dimlon

             if (vtype.eq.'rms') then
                  if (order.eq.1) then
                  varrms2d(iloop,jloop,k)=
     &            temp2d(iloop,dimlat+1-jloop,1) 
                  elseif (order.eq.2) then
                  varrms2d2(iloop,jloop,k)=
     &            temp2d(iloop,dimlat+1-jloop,1)
                  endif 
             endif
                               
            enddo
           enddo

      return
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine getsd_3d(nf,vtype,k,varname,up_low,order)

c getsd_3D reads the 3D RMS and ARMS variables in up or low netcdf files 
c and fills corresponding varrms3d(k) and vararms3d(k) arrays (which are 
c defined as common in constants_mcd.inc)
      implicit none
      
      include "constants_mcd.inc"
      include "netcdf.inc"
    
c     inputs      
      integer      nf      ! Input NetCDF file ID
      character(len=*) ::  vtype
      integer      k
      character(len=*) :: varname
      integer      up_low
      integer      order

c     local variables
      integer      iloop,kloop,jloop
      integer      ierr
      integer      sd3didin(nbsd3d+1)     ! to store ID # of variable
      real         temp3dlow(dimlon,dimlat,low,dimuti)
! temp3dlow(:,:,:,:) to temporarely store variable from "low" file
      real         temp3dup(dimlon,dimlat,up,dimuti)
! temp3dup(:,:,:,:) to temporarely store variable from "up" file

!      write(out,*)'getsd_3d: k varname',k,trim(varname)

c  Get ID of variable from NetCDF file      
         ierr=NF_INQ_VARID(nf,varname,sd3didin(k))
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GETSD_3D Error:"
             write(out,*) NF_STRERROR(ierr),varname
           endif
           stop 
         endif

      if (up_low.eq.up) then
      ! Read variable from "up" file
         ierr = NF_GET_VAR_REAL(nf,sd3didin(k),temp3dup)
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GETSD_3D Error:"
             write(out,*) NF_STRERROR(ierr)
           endif
           stop
         endif

         ! Copy variable
          do kloop=low+1,dimlevs
           do jloop=1,dimlat
            do iloop=1,dimlon

             if (vtype.eq.'rms') then
                  if (order.eq.1) then
                    varrms3d(iloop,jloop,kloop,k)=
     &                temp3dup(iloop,dimlat+1-jloop,kloop-low,1) 
                  elseif (order.eq.2) then
                    varrms3d2(iloop,jloop,kloop,k)=
     &                temp3dup(iloop,dimlat+1-jloop,kloop-low,1) 
                  endif
             else if (vtype.eq.'arms') then
                  if (order.eq.1) then
                    vararms3d(iloop,jloop,kloop,k)=
     &                temp3dup(iloop,dimlat+1-jloop,kloop-low,1) 
                  elseif (order.eq.2) then
                    vararms3d2(iloop,jloop,kloop,k)=
     &                temp3dup(iloop,dimlat+1-jloop,kloop-low,1) 
                  endif
             endif
             
            enddo
           enddo
          enddo

      else ! Read variable from "low" file
         ierr = NF_GET_VAR_REAL(nf,sd3didin(k),temp3dlow)
         if(ierr.ne.NF_NOERR) then
           if (output_messages) then
             write(out,*) "GETSD_3D Error:"
             write(out,*) NF_STRERROR(ierr)
           endif
           stop
         endif

         ! Copy variable (and permute latitude index)
          do kloop=1,low
           do jloop=1,dimlat
            do iloop=1,dimlon

             if (vtype.eq.'rms') then
                  if (order.eq.1) then
                  varrms3d(iloop,jloop,kloop,k)=
     &              temp3dlow(iloop,dimlat+1-jloop,kloop,1) 
                  elseif (order.eq.2) then
                  varrms3d2(iloop,jloop,kloop,k)=
     &              temp3dlow(iloop,dimlat+1-jloop,kloop,1) 
                  endif
             else if (vtype.eq.'arms') then
                  if (order.eq.1) then
                  vararms3d(iloop,jloop,kloop,k)=
     &              temp3dlow(iloop,dimlat+1-jloop,kloop,1) 
                  elseif (order.eq.2) then
                  vararms3d2(iloop,jloop,kloop,k)=
     &              temp3dlow(iloop,dimlat+1-jloop,kloop,1) 
                  endif
             endif
              
            enddo
           enddo
          enddo
       endif! of if (up_low.eq.up)         

      
      return
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine var3d(a,lon,lat,zsurface,levhi,levlow,levweight,pratio,
     &                 sheight,p_pgcm,sigma,ps_psgcm,utime,name,ier,
     &                 itimint,wl,wh)

c     Retrieve the value of 3-d variable=name at longitude=lon, latitude=lat
c     and Universal time=utime. The height (above local surface) of the
c     variable is given by zsurface, although only levhi, levlow and levweight
c     are needed for most vertical interpolations.
c     Bilinear interpolation (see called routine profi) is used to get user
c     specified longitude and latitude.
c     Linear interpolation of variables with height is done in the vertical
c     direction, exept for density (and similar fields), which is
c     interpolated linearly with pressure.
c     When below first atmospheric level, specific "interpolations" are used
c     for temperature and horizontal winds (Monin-Obukhov)
c     For this reason, var3d now needs sigma in input
c     
      implicit none

      include "constants_mcd.inc"

c     inputs
      real         lon         ! longitude (east) of point (deg.)
      real         lat         ! latitude of point (deg.)
      real         zsurface    ! altitude above surface of point (m)
      integer      levhi       ! database level upper bound
      integer      levlow      ! database level lower bound
      real         levweight(2)! level weight for vertical interpolation
c     (1) for linear in height (2) for linear in pressure
      real         pratio   !ratio of pressure to extreme value if out of range
      real         sheight(dimlevs) ! altitude of GCM sigma levels
      real         p_pgcm(dimlevs) ! high res to GCM pressure ratios
      real         sigma(dimlevs) ! sigma levels in MCD grid
      real         ps_psgcm ! High res to GCM surface pressure ratio 
      real         utime    !Universal time (0. to 24. hrs)=local time at lon=0
      integer      itimint  !seasonal interpolation flags
      real         wl,wh    !seasonal interpolation weights 
      character(len=*):: name     !Name of variable

c     outputs
      integer       ier        !error flag (0=OK, 1=NOK)
      real          a          !value of variable

c     local variables
      real    profile(dimlevs) ! profile in MCD grid 
      character*16 tmpname     ! Name of variable (for call to var2d)
      real    tsurf,ps ! surface temperature and surface pressure
      real    zu(dimlevs),zv(dimlevs) ! winds (used for temperature interpolation)
      real    hfmax,zmax,wstar_in ! pbl parameters from the thermals for interpolation
      real    pplev(dimlevs),ph(dimlevs) !level pressure and potential temp
      real    T_out,u_out,ustar,tstar,vvv,vhf ! outputs of pbl_parameters
      real    rcp ! ratio of r over cp
!      real    pz0 !surface roughness length interpolated at output location
      real    norm1 ! norm of velocity vector at first latitude level
      integer ierr
      real          vmr_gcm(dimlevs)  ! vmr at GCM levels
      real          rho_gcm(dimlevs)  ! density at GCM levels
      real          temp_gcm(dimlevs) ! temperature at GCM levels
      real          coef              ! coef accounting composition
      real          xratio,x,xi,xsum  ! coef for extrapolation
      integer       igas
      
      real          g, g0, a0
      parameter (a0=3396.E3,g0=3.7257964)      

      ! calls to 2D (and 3D) variables from pbl_parameter are handled in var3d to avoid needless complexity in var2d.
      if ((name.eq.'surfstress').or.(name.eq.'sensib_flux')
     &    .or.(name.eq.'vvv').or.(name.eq.'vhf')) then
          goto 1000
      endif

      ier = 0

      ! get horizontally interpolated GCM vertical profile of variable 'name'
      call profi(profile,lon,lat,utime,name,ier,
     &           itimint,wl,wh,0,levlow,levhi)
      if (ier.ne.0) then
        if (output_messages) then
         write(out,*)'VAR3D Error : profile not available for variable'
     &        ,name
        endif
        ier=1
        goto 9999
      endif


      if ((name.eq.'rho').or.
     &    (name.eq.'rmsrho').or.(name.eq.'armsrho').or.
     &    (name.eq.'armspressure')) then
c     interpolate densities linearly in pressure
         a=profile(levlow)*p_pgcm(levlow)+levweight(2)*
     &    (profile(levhi)*p_pgcm(levhi)-profile(levlow)*p_pgcm(levlow))
c     correction for densities which are out of sigma range
c         a=pratio*a
         if(.not.(levlow.eq.levhi.and.levlow.eq.dimlevs)) then
          a=pratio*a   
         else    
          g = g0*(a0/(a0+sheight(dimlevs)))**2
          call profi(temp_gcm,lon,lat,utime,'temp',ierr,itimint,wl,wh,
     &               0,dimlevs,dimlevs)          
          tmpname = 'ps'
          call var2d(ps,lon,lat,utime,tmpname,ierr,
     &               itimint,wl,wh,ps_psgcm)
          xratio = 0.0
          do igas=1,nbvar3d
           if(ma3d(igas).lt.0) cycle   
           call profi(vmr_gcm,lon,lat,utime,varname3d(igas),
     &                ierr,itimint,wl,wh,0,dimlevs,dimlevs)
           coef = ma3d(igas)*sigma(dimlevs)*ps/
     &          (8314.4598*temp_gcm(dimlevs)*a)  ! a is rho or pressure
           if(name.eq.'armspressure') coef = 1.0
c           xratio = xratio+coef*vmr_gcm(dimlevs)*
c     &              exp((sheight(dimlevs)-zsurface)*g/
c     &              (8314.4598/ma3d(igas)*temp_gcm(dimlevs)))
c           
           xratio = xratio+coef*vmr_gcm(dimlevs)*
     &        exp(-(a0+sheight(dimlevs))*(zsurface-sheight(dimlevs))/
     &             (a0+zsurface)*
     &              g/(8314.4598/ma3d(igas)*temp_gcm(dimlevs)))           
          enddo
          a=xratio*a
         endif ! of if(.not.(levlow.eq.levhi.and.levlow.eq.dimlevs))
      elseif ((name.eq.'vmr_h2ovapor').or.
     &        (name.eq.'vmr_co2').or.(name.eq.'vmr_co').or.
     &        (name.eq.'vmr_n2').or.(name.eq.'vmr_ar').or.
     &        (name.eq.'vmr_o').or.(name.eq.'vmr_o2').or.
     &        (name.eq.'vmr_o3').or.(name.eq.'vmr_he').or.
     &        (name.eq.'vmr_h').or.(name.eq.'vmr_h2')) then
        ! Default linear interpolation between levels:
         a=profile(levlow)+(profile(levhi)-profile(levlow))*levweight(1)
        ! If abowe last atmospheric level
         if ((levlow.eq.levhi).and.(levlow.eq.dimlevs)) then
          g = g0*(a0/(a0+sheight(dimlevs)))**2
          call profi(temp_gcm,lon,lat,utime,'temp',ierr,itimint,wl,wh,
     &               0,dimlevs,dimlevs)
          xsum = 0.0
          do igas=1,nbvar3d
           if(ma3d(igas).lt.0) cycle   
           call profi(vmr_gcm,lon,lat,utime,varname3d(igas),
     &                ierr,itimint,wl,wh,0,dimlevs,dimlevs)
c           x = vmr_gcm(dimlevs)*
c     &              exp((sheight(dimlevs)-zsurface)*g/
c     &              (8314.4598/ma3d(igas)*temp_gcm(dimlevs)))
c
           x = vmr_gcm(dimlevs)*
     &        exp(-(a0+sheight(dimlevs))*(zsurface-sheight(dimlevs))/
     &             (a0+zsurface)*
     &              g/(8314.4598/ma3d(igas)*temp_gcm(dimlevs)))
           
           if(name.eq.varname3d(igas)) xi = x
           xsum = xsum + x
          enddo            
          a = xi/xsum             
         endif ! of if ((levlow.eq.levhi).and.(levlow.eq.dimlevs))
      elseif ((name.eq.'temp')) then
        ! Default linear interpolation between levels:
        a=profile(levlow)+(profile(levhi)-profile(levlow))*levweight(1)
        ! If below first atmospheric level, use linear
        ! interpolation with surface temperature
        if ((levlow.eq.levhi).and.(levlow.eq.1)) then
           goto 1000
        endif ! of if ((levlow.eq.levhi).and.(levlow.eq.1))
      elseif ((name.eq.'u').or.(name.eq.'v')) then
        ! Default linear interpolation between levels:
        a=profile(levlow)+(profile(levhi)-profile(levlow))*levweight(1)
        ! If below first atmopheric level, mimic log boundary layer
        ! with roughness length z_0; below z_0, set wind velocity to zero
        if ((levlow.eq.levhi).and.(levlow.eq.1)) then
           goto 1000
        endif
      elseif ((name.eq.'vvv').or.(name.eq.'vhf')) then
           goto 1000
      else ! most variables are interpolated linearly in height
        a=profile(levlow)+(profile(levhi)-profile(levlow))*levweight(1)
      endif
c

c ---------------------------------------------------------------------
c Special case : Monin Obukhov interpolation
c for variables : u,v, theta below z1, vhf, vvv, surfstress, sensibflux
c----------------------------------------------------------------------

      return

1000  continue
          ! get surface temperature at same location
          ! in the PBL, rcp is constant
          rcp = 0.2567930
          tmpname="tsurf"
          call var2d(tsurf,lon,lat,utime,tmpname,
     &               ierr,itimint,wl,wh,ps_psgcm)
          ! get surface pressure at same location
          tmpname="ps"
          call var2d(ps,lon,lat,utime,tmpname,
     &               ierr,itimint,wl,wh,ps_psgcm)
          ! get wind profile at same location
          tmpname="u"
          call profi(zu,lon,lat,utime,tmpname,ierr,itimint,wl,wh,0,0,0)
          tmpname="v"
          call profi(zv,lon,lat,utime,tmpname,ierr,itimint,wl,wh,0,0,0)
          ! get thermals parameters at same location
          tmpname="zmax"
          call var2d(zmax,lon,lat,utime,tmpname,
     &               ierr,itimint,wl,wh,ps_psgcm)
          tmpname="hfmax"
          call var2d(hfmax,lon,lat,utime,tmpname,
     &               ierr,itimint,wl,wh,ps_psgcm)
          tmpname="wstar"
          call var2d(wstar_in,lon,lat,utime,tmpname,
     &               ierr,itimint,wl,wh,ps_psgcm)
          ! get temp profile at same location
          tmpname="temp"
          call profi(profile,lon,lat,utime,tmpname,ier,itimint,wl,wh,
     &               0,0,0)
          ! compute pressure profile at same location
          pplev(:) = sigma(:)*ps
          ! compute corresponding potential temperature profile
          ph(:) = profile(:)*((ps/pplev(:))**rcp)
          ! interpolate
          call pbl_parameters(1,dimlevs,ps,pplev,z_0,
     & sheight,zu,zv,wstar_in,hfmax,zmax,tsurf,ph,zsurface,1,
     & T_out,u_out,ustar,tstar,vhf,vvv)

      if (name.eq.'vvv') then
          a=vvv
      elseif (name.eq. 'vhf') then
          a=vhf
      elseif ((name.eq.'u').or.(name.eq.'v')) then
          ! N.B: u_out is the norm of the velocity vector
          norm1=sqrt(zu(1)*zu(1)+zv(1)*zv(1)) ! norm of velocity at atm level 1
          if (norm1.eq.0) then
            ! should almost never happen, but just in case...
            a=0
          else
            if (u_out.gt.norm1) then
              ! should never be the case, but might, because of roundoff
              ! effects, when norm1 is very near (but not exactly) zero
              u_out=norm1
            endif
            if (name.eq.'u') then
              a=zu(1)*(u_out/norm1)
            else
              a=zv(1)*(u_out/norm1)
            endif
          endif ! of if (norm1.eq.0)
      elseif (name.eq.'temp') then
          a=T_out
      elseif (name.eq.'surfstress') then
          tmpname="rho" 
          call profi(profile,lon,lat,utime,tmpname,ier,itimint,
     &               wl,wh,0,1,1)
          a= profile(1)*p_pgcm(1)*ustar*ustar
      elseif (name.eq.'sensib_flux') then
          ! at surface, cp is constant == 744.499
          tmpname="rho"
          call profi(profile,lon,lat,utime,tmpname,ier,itimint,
     &               wl,wh,0,1,1)
          a= profile(1)*p_pgcm(1)*744.499*ustar*tstar
      endif

      return

c     Error handling
 9999 a=0.
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine max2d(a,lon,lat,name,ier,itimint,wl,wh)

c     Retrieve the daily maximum value of 2-d variable=name at 
c     longitude=lon, latitude=lat

      implicit none

      include "constants_mcd.inc"

c     inputs
      real         lon     !west longitude of point
      real         lat     !latitude of point
      character*16 name    !Name of variable
      integer      itimint !seasonal interpolation flag 
      real         wl,wh   !seasonal interpolation weightings

c     outputs
      integer      ier     !Error flag (0=OK, 1=NOK)
      real         a       !maximum of variable

c     local variables
      integer      i
      real         x
      real         xmax
      real         utime
      
      ier=0
      xmax=-1.e20
      do i=0,22,2
        utime=float(i)
        call var2d(x,lon,lat,utime,name,ier,itimint,wl,wh,1.0)
        if(ier.ne.0) then
           ier=1
           if (output_messages) then
             write(out,*)"MAX2D Error : impossible to find ",name,
     &                 " maximum"
           endif
           goto 9999
        endif
        if (x.gt.xmax) xmax=x
      enddo
      a=xmax

      return

c     Error handling
 9999 a=0.
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine min2d(a,lon,lat,name,ier,itimint,wl,wh)

c     Retrieve the daily minimum value of 2-d variable=name at 
c     longitude=lon, latitude=lat

      implicit none

      include "constants_mcd.inc"

c     inputs
      real         lon     !longitude of point
      real         lat     !latitude of point
      character*16 name    !Name of variable

c     output
      integer      ier     !Error flag (0=OK, 1=NOK)
      real         a       !minimum of variable
      real         wl,wh   !seasonal interpolation weightings
      integer      itimint !seasonal interpolation flag

c     local variables
      integer      i
      real         x
      real         xmin
      real         utime

      ier=0
      xmin=1.e20
      do i=0,22,2
        utime=float(i)
        call var2d(x,lon,lat,utime,name,ier,itimint,wl,wh,1.0)
        if(ier.ne.0) then
           ier=1
           if (output_messages) then
             write(out,*)"MIN2D Error : impossible to find ",name,
     &                 " minimum"
           endif
           goto 9999
        endif
        if (x.lt.xmin) xmin=x
      enddo
      a=xmin
      return

c     Error handling
 9999 a=0.
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      function ran1(idum)
c     "Minimal" random number generator of Park and Miller with Bays-Durham
c     shuffle and added safeguards 
c     Return a uniform random deviate between 0.0 and 1.0 (exclusive of the
c     endpoints values). Call with idum a negative integer to initialize;
c     thereafter, do not alter idum between successive deviates in a sequence.
c     RNMX should approximate the largest floating value that is less than 1.

      implicit none

c input:
      integer idum ! Seed for random number generator (if negative)
                   ! the value of idum is altered by the function

c returned value:
      real ran1

c local variables:
      integer IA,IM,IQ,IR,NTAB,NDIV
      real AM,EPS,RNMX
      parameter (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      integer j,k
      integer iv(NTAB),iy
      save iv,iy
      data iv /NTAB*0/, iy /0/

      if ((idum.le.0).or.(iy.eq.0)) then
         idum=max(-idum,1)
         do j=NTAB+8,1,-1
            k=idum/IQ
            idum=IA*(idum-k*IQ)-IR*k
            if (idum.lt.0) idum=idum+IM
            if (j.le.NTAB) iv(j)=idum
         end do
         iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!      subroutine dustmix(scena,lat,Ls,dod,psurf,p,dust_mmr)
!     EM 2012: not needed for MCDv5
c     Subroutine used to reconstruct the dust mixing ratio 
c     assumed in the GCM run, at a given pressure level. 
c     F. Forget 2005

!      implicit none

!      include "constants_mcd.inc"

c     INPUT
!      integer scena  ! scenario
!      real lat  ! latitude  (degree) 
!      real Ls       ! solar longitude (deg)
!      real dod      ! dust optical depth
!      real psurf      ! surface pressure (Pa)
!      real p        ! local pressure (Pa)
c     OUTPUT
!      real dust_mmr  ! dust mass mixing ratio (kg/kg)
c     Local Variable
!      real radius, Qext, rho_dust
!      real zlsconst, topdust,topdust0,qextrhor,g
!      real zp, expfactor, xlat
!      logical firstcall
!      data firstcall /.true./
!      save firstcall,qextrhor,g
!      real pi
!      parameter   (pi=3.14159265359)


!      topdust=0. ! dummy initialization to get rid of compiler warnings

c     initialisation
c     --------------
!      if (firstcall) then
!         radius=1.8e-6  ! effective radius of dust (m)
!         Qext=3.  ! dust visible single scattering extinction coeff.
!         rho_dust=2500.  ! Mars dust density (kg.m-3)
!         qextrhor= (3./4.)*Qext/(rho_dust*radius)
!         g=3.72
!         firstcall = .false. 
!      end if
!      xlat = lat*pi/180.

c     Altitude of the top of the dust layer
c     -------------------------------------

!      zlsconst=SIN(ls*pi/180. - 2.76)
      ! TEMPORARILY REMOVE TESTS BELOW; MAKE NO SENSE WITH NEW DUST SCENARIOS
!      if (scena.eq.8) then
!             topdust= 30.        ! constant dust layer top = 30 km
!      else if (scena.eq.7) then      ! "Viking" scenario
!            topdust0=60. -22.*SIN(xlat)**2
!            topdust=topdust0+18.*zlsconst
!      else if(scena.le.6) then         !"MGS" scenario (MY24, storm)
!            topdust=60.+18.*zlsconst
!     &               -(32+18*zlsconst)*sin(xlat)**4
!     &               - 8*zlsconst*(sin(xlat))**5
!
!      else
!      topdust=70.
!      if (scena.gt.10) then
!        if (output_messages) then
!          write(out,*) "DUST_MMR error: problem with scena"
!        endif
!        stop 
!      endif

c     Computing dust mass mixing ratio
c     --------------------------------
!      if(p.gt.700./(988.**(topdust/70.))) then
!           zp=(700./p)**(70./topdust)
!           expfactor=max(exp(0.007*(1.-max(zp,1.))),1.E-10)
!      else
!           expfactor=1.E-10
!      endif
c     qextrhor is  (3/4)*Qext/(rho*reff)
!      dust_mmr = expfactor* dod * g / (psurf * qextrhor)
!      return
!      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine dust_deposition(ps,temp,rho,mmr,reff,dep)
      
      implicit none
      
      real,intent(in) :: ps ! surface pressure
      real,intent(in) :: temp ! atm. temperature in first layer
      real,intent(in) :: rho ! atm. density in first layer
      real,intent(in) :: mmr ! dust mass mixing ratio
      real,intent(in) :: reff ! dust effective radius
      real,intent(out) :: dep ! dust deposition on flat surface kg / m2 / s
      
      ! local variables
      logical,save :: firstcall=.true.
      real,save :: a,b,xnu,nueff
      real :: r_sed
      real :: vstokes
      ! physical constants:
      real,parameter :: visc=1.e-5 ! Molecular viscosity of CO2 gas (N.s.m-2)
      real,parameter :: molrad=2.2e-10 ! effective CO2 gas molecular radius (m)
      real,parameter :: rhodust=2500. ! Mars dust density (kg.m-3)
      real,parameter :: g=3.72 ! gravity
      
      IF (firstcall) THEN

c       Preliminary calculations for sedimentation velocity :
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c       Constant to compute stokes speed simple formulae
c       (Vstokes =  b * rho* r**2   with   b= (2/9) * rho * g / visc
         b = 2./9. * g / visc

c       - Constant  to compute gas mean free path
c        l= (T/P)*a, with a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
         a = 0.707*8.31/(4*3.1416* molrad**2  * 6.023e23)

c       - Correction to account for non-spherical shape (Murphy et al.  1990)
c   (correction = 0.85 for irregular particles, 0.5 for disk shaped particles)
         a = a  * 0.5

c       Correction to compute "sedimentation" effective radius r_sed 
c       from effective radius reff with the lognormal distribution
c       (r_sed = xnu * reff)
        nueff= 0.5   ! dust size distribution effective variance in GCM runs 
        xnu =  (1 + nueff)**2

      ENDIF ! of IF (firstcall)
      
      ! compute dust deposition on flat surface
      r_sed=xnu*reff
      vstokes = b * rhodust * r_sed**2 *
     &             (1 + 1.333* ( a*temp/ps )/r_sed)
      dep = mmr * rho * vstokes
      
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine air_properties(T,vmr_co2,vmr_n2,vmr_ar,vmr_o,vmr_co,
     &                       Cp,gamma,viscosity,Rgas)

      implicit none
                                                                                
c     Subroutine used to compute some useful gas properties
c     from the climate database data
c     F. Forget 2005 / E. Millour 2012
                                                                                
c     INPUT
      real,intent(in) ::  T ! temperature   (K)
      ! gas volume mixing ratios (mol/mol)
      real,intent(in) :: vmr_co2,vmr_n2,vmr_ar,vmr_o,vmr_co 
      
c     OUTPUT
      real,intent(out) :: cp   ! Heat capacity at constant pressure (J kg-1? K-1)
      real,intent(out) :: gamma   ! Cp/Cv Ratio of specific heats 
      real,intent(out) :: viscosity    ! air viscosity (N s m-2)
      real,intent(out) :: Rgas   ! gas constant  (J kg-1 K-1)
      
c     Local Variable
c     --------------

!      real vmr_ar   ! Argon vmr
c     Molecular heat capacity (J mol-1 K-1)
      real Cp_co2, Cp_co, Cp_n2, Cp_ar, Cp_O
c     Coefficient to compute conductivity = A*T**0.69  (W m-1 K-1)
!      real A_co2, A_co, A_n2, A_ar, A_O
c     Molecular mass
      real M_co2, M_co, M_n2, M_ar, M_O
!      real conduct,
      real meanmolmass,svmr
      logical firstcall
      data firstcall /.true./
      save firstcall
    
      save Cp_co2, Cp_co, Cp_n2, Cp_ar, Cp_O
!      save A_co2, A_co, A_n2, A_ar, A_O
      save M_co2, M_co, M_n2, M_ar, M_O

      integer i,j
      real,save :: M(5)
      real,save :: Phi(5,5)
      real :: nu(5),mfrac(5)
      real denom

c     initialisation
c     --------------
      if (firstcall) then
c           Molecular mass (kg mol-1)
            M_co2=44.01e-3
            M_co=28.01e-3
            M_n2=28.013e-3
            M_Ar=39.95e-3
            M_o=16.e-3
c           Molecular heat capacity (J mol-1 K-1)
            Cp_co2 = 37.
            Cp_co  = 29
            Cp_n2 = 29
            Cp_ar = 20
            Cp_o = 21
c           Coeff for conductivity :
!            A_co2=3.07e-4
!            A_o=7.6e-4
!            A_n2=5.6e-4
!            A_co=4.87e-4
!            A_ar=3.4e-4
! Build coefficients for viscosity: 1: CO2, 2: CO, 3: N2, 4: Ar, 5: O
             ! molar masses
             M(1)=44 ; M(2)=28 ; M(3)=28 ; M(4)=39.95 ; M(5)=16
             !Phi(i,j) coefficients:
             do i=1,5
               do j=1,5
                 phi(i,j)=sqrt(M(j)/M(i))
               enddo
             enddo

            firstcall=.false.
      end if

c     Estimating Argon mixing ratio ! Not for MCDv5 where vmr_ar is known
c     -----------------------------   
!      vmr_ar = 0.6 *vmr_n2

! 03/2011: improvement, use S. Lebonnois' formulation for Cp(T) for CO2
      Cp_co2=1.e3*M_co2*(T/460.)**0.35

c     Mean Molecular mass
c     ------------------
      Svmr = vmr_co2+vmr_n2+vmr_o+vmr_co+vmr_ar
      meanmolmass= (M_co2*vmr_co2+M_n2*vmr_n2+M_o*vmr_o
     &      + M_co*vmr_co+M_ar*vmr_ar)/Svmr
      Rgas = 8.3144/meanmolmass

c     Computing Cp (Neglecting H2 => not valid if very high)
c     ------------
      Cp = (Cp_co2*vmr_co2+Cp_n2*vmr_n2+Cp_o*vmr_o
     &      + Cp_co*vmr_co+Cp_ar*vmr_ar)/Svmr


c     Conversion from J mol-1 K-1 to J kg-1 K-1
      Cp= Cp / meanmolmass

c     Computing gamma
      gamma = 1./(1-Rgas/Cp)

c     Computing viscosity
c     -------------------
c     Therma molecular Conductivity (W m-1 K-1)
!      conduct = (A_co2*vmr_co2+A_n2*vmr_n2+A_o*vmr_o
!     &        + A_co*vmr_co + A_ar*vmr_ar)/Svmr
!      conduct=T**0.69 * conduct
!      viscosity = conduct / (0.25*(4*Cp + 5*Rgas))  
  
! 03/2011: use (temperature dependent) viscosity of gases
      ! CO2: Sutherland formula
      nu(1)=1.572e-6*(T**1.5)/(T+240.)
      ! CO: Sutherland formula
      nu(2)=1.428e-6*(T**1.5)/(T+118.)
      ! N2: Sutherland formula
      nu(3)=1.407e-6*(T**1.5)/(T+111.)
      ! Ar: Van Itterbeek & Van Paemel, Physica V,no 10, 1938
      nu(4)=7.592e-6*(T/90.)**0.883
      ! O: Dalgarno & Smith Planet. Space Sci., Vol 9, pp. 1-2, 1962
      nu(5)=3.34e-7*(T)**0.71
      ! 
      mfrac(1)=vmr_co2
      mfrac(2)=vmr_co
      mfrac(3)=vmr_n2
      mfrac(4)=vmr_ar
      mfrac(5)=vmr_o
      
      viscosity=0
      do i=1,5
       denom=0
       do j=1,5
         denom=denom+mfrac(j)*phi(i,j)
       enddo
       viscosity=viscosity+(mfrac(i)*nu(i))/denom
      enddo
      
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine build_sigma_hr(sigma_gcm,ps_gcm,ps_hr,sigma_hr,p_pgcm)
      
      implicit none
      include "constants_mcd.inc"

c     inputs
      real sigma_gcm(dimlevs) ! GCM sigma levels
      real ps_gcm             ! GCM surface pressure
      real ps_hr              ! High res surface pressure
c     outputs
      real sigma_hr(dimlevs)  ! High res sigma levels
      real p_pgcm(dimlevs)    ! high res to GCM pressure ratios
      
c     local variables
      integer l
      real x  ! lower layer compression (-0.9<x<0) or dilatation (0.<x<0.9)
      real rp     ! surface pressure ratio ps_hr/ps_gcm
      real deltaz   ! corresponding pseudo-altitude difference (km)
      double precision f  ! coefficient f= p_hr / p_gcm
      double precision z  ! altitude of transition of p_hr toward p_gcm (km)

c 1. Coefficients
      rp=ps_hr/ps_gcm
      deltaz=-10.*log(rp)
      x = min(max(0.12*(abs(deltaz)-1.),0.),0.8)
      if(deltaz.gt.0) x=-x
      z=max(deltaz + 3.,3.)

      do l=1,dimlevs
        f=rp*sigma_gcm(l)**x
c        f=f+(1-f)*0.5*(1+tanh(6.*(-10.*log(sigma_gcm(l))-z)/z))
c        sigma_hr(l)=f*sigma_gcm(l)/rp
        p_pgcm(l)=f+(1.d0-f)*0.5d0*
     &            (1.d0+tanh(6.d0*(-10.d0*log(dble(sigma_gcm(l)))-z)/z))
        sigma_hr(l)=p_pgcm(l)*sigma_gcm(l)/rp
      enddo

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine pres0(dset,dust,lat,lon,solar,utime,ps_MCD,oro_MCD,
     &                 wl,wh,pres,alt,ierr)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Purpose:
c =======
c  Pres0 si a subroutine designed to yield high resolution MOLA topography
c  and recompute surface pressure accordingly. 
c  It uses:
c     1) Reference pressure measurements at Viking Lander 1 site (file VL1.ls)
c     2) High resolution MOLA topography (file mola32.nc)
c     3) Data from the Mars Climate Database
c
c Built from the v4.1 pres0 tool + improvements
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      include "constants_mcd.inc"

c inputs:
      character*(*) dset ! Path to MCD datafiles 
      integer dust  ! dust scenario
      real lon      ! Longitude coordinate of the point (East degrees)
      real lat      ! Latitude coordinate of the point (North degrees)
      real solar    ! Solar longitude (degrees)
      real utime    ! Universal time (hours) (time at lon=0)
      real ps_MCD   ! surface pressure (from MCD)
      real oro_MCD  ! orography from MCD
      real wl,wh    ! season-wise interpolation weights

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 temp     ! Temperature at altitude ~1km

      character*16 name    !Name of variable to retrieve from MCD
      integer ier
      real profile(dimlevs) ! to temporarily store temperature profile
      integer ref_alt ! reference altitude index for temperature retrieval
      parameter (ref_alt=7) ! 7th level is ~ 1km above surface
      
! Ensure that (local) longitude zlon is in [0:360]
! since given "lon" is in [-180:180]

      zlon=lon
      if(zlon.lt.0.) zlon =zlon + 360.


      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,dust,solar,wl,wh,factcor,ierr)
      if (ierr.ne.0) return

! 3. Get MCD temperature (at ~1km above surface)
      name='temp'
      call profi(profile,lon,lat,utime,name,ier,
     &           1,wl,wh,0,ref_alt,ref_alt)
      ! keep temperature of ref_alt layer (~1km)
        temp=profile(ref_alt)
      ! Ehouarn test:

! 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
      H=191*temp/3.73

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

      end
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine read_vl(dset,tab,ierr)
!Routine to read the VL1 and VL3 mean viking pressure files.
      implicit none
      include "constants_mcd.inc"

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

c local variables
      integer vik_num,count
      character*140 vlfile ! (smoothed) surface pressure measured by VL
c      data vlfile/'MCD_DATA/VL1.ls'/

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

      if (ierr.ne.0) then
        if (output_messages) then
          write(out,*) "READ_VL Error: Failed opening file ",
     &                 trim(vlfile)
        endif
        return
      endif

C 2. Read the 669 lines of file
      do count=1,669
        read(vik_num,*) tab(count,1),tab(count,2)
      enddo

C 3. Close file and return
      close (vik_num)

c 4. Open file VL3.ls
      vik_num=18
      vlfile=dset//'VL3.ls'
      open(vik_num,file=vlfile,status="old",iostat=ierr,
     &     form="formatted")

      if (ierr.ne.0) then
        if (output_messages) then
          write(out,*) "READ_VL Error: Failed opening file ",
     &                 trim(vlfile)
        endif
        return
      endif

C 2. Read the 669 lines of file
      do count=1,669
        read(vik_num,*) tab(count,1),tab(count,3)
      enddo

C 3. Close file and return
      close (vik_num)

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine calc_factcor(dset,dust,solar,wl,wh,factcor,ierr)
c Calculate a correction factor between Viking1 surface pressure data
c and MCD output at the same point
c Note: Position of VL1 is taken from Konopliv et al. 2006 
c       and is expressed in the IAU2000 reference system.
c       Corresponding altitude was computed using bilinear interpolation of
c       MOLA 32 topography data (and using areoid computed from the mgm1025
c       spherical harmonics coefficients) 
!     In case of dust storm scenario, use VL3 input not VL1
      implicit none
      include "constants_mcd.inc"

! Arguments:
!     inputs:
      character*(*) dset ! path to MCD datafiles
      integer dust       ! dust scenario
      real solar         ! Solar longitude (degrees)
      real wl,wh         ! weights for interpolation
!     outputs:
      real factcor        ! correction factor
      integer ierr        ! status ierr=0 if everything OK

! Local variables
      ! Viking 1 coordinates
      real vik1_lat, vik1_lon, vik1_alt
      parameter(vik1_lat=22.269628,vik1_lon=312.050381,
     &          vik1_alt=-3637.1396)
      ! Vinking 3 coordinates
      real,parameter :: vik3_lat=-45.
      real,parameter :: vik3_lon=360.
      real,parameter :: vik3_alt=1322.000
      real :: vik_lat, vik_lon, vik_alt
      real vik_MCD_oro ! MCD orography at VL coordinates
      real prescalc(12) ! surface pressure at VL1 site, at times 2,4,...24 hours
      real tempcalc(12) ! temperature above VL1, at times 2,4,...24 hours
      integer ref_alt ! reference altitude index for temperature retrieval
      parameter (ref_alt=7)
      real H        ! scale height
      real presmean ! mean value (over a day) of surface pressure
      real lsinf,lssup,psinf,pssup,presvik !,alt
      integer iloop
      logical firstcall
      data firstcall/.true./
      save firstcall
      real,save :: vik_tab(669,3)
      integer vik_ind ! index for vik_tab() column to use; 2:VL1 or 3:VL3

      character*16 name    ! Name of variable to retrieve from MCD
      integer ier
      real profile(dimlevs) ! to temporarily store temperature profile
      
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. Read surface pressure and temperature (1km above surface) from MCD, 
C    at VL1 site, at universal times (of day) 2, 4, ... 24 hours.

! depending on dust scenario, use VL1 or VL3 as reference
      if ((dust.ge.4).and.(dust.le.6)) then
        vik_lat=vik3_lat
        vik_lon=vik3_lon
        vik_alt=vik3_alt
        vik_ind=3 ! index for vik_tab() column to use
      else ! regular case use VL1
        vik_lat=vik1_lat
        vik_lon=vik1_lon
        vik_alt=vik1_alt
        vik_ind=2 ! index for vik_tab() column to use
      endif

C 2.1 Read pressure and temperature from MCD at utime 2, 4, ... 24
      do iloop=1,12
        name='ps'
        ! N.B. vard2d expects longitudes in [-180:180]
        call var2d(prescalc(iloop),vik_lon-360.,vik_lat,real(2*iloop),
     &              name,ier,1,wl,wh,1.0)
        ! retrieve temperature profile
        name='temp'
        call profi(profile,vik_lon-360.,vik_lat,2.0*iloop,
     &              name,ier,1,wl,wh,0,ref_alt,ref_alt)
        ! keep temperature of ref_alt layer (~1km)
        tempcalc(iloop)=profile(ref_alt)
c        write(out,*)'CALC_FACTCOR: i=',iloop,
c     &            ' prescalc(i)=',prescalc(iloop),
c     &            ' tempcalc(i)=',tempcalc(iloop)
      enddo  

C 2.2 Get MCD orography corresponding to VL1 site
      name='orography'
      call var2d(vik_MCD_oro,vik_lon-360.,vik_lat,0.0,
     &           name,ier,1,wl,wh,1.0)
c      write(out,*) 'CALC_FACTCOR: vik_MCD_oro=',vik_MCD_oro


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(-(vik_alt-vik_MCD_oro)/H)
      enddo

C 4. Compute average (over the day) pressure
      presmean=0.
      do iloop=1,12
         presmean = presmean + prescalc(iloop)/12.
      enddo
c      write(out,*)'CALC_FACTCOR: presmean=',presmean

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,vik_ind)
         pssup=vik_tab(1,vik_ind)
      else if (solar.ge.vik_tab(669,1)) then
         lsinf=vik_tab(669,1)
         lssup=vik_tab(1,1)+360
         psinf=vik_tab(669,vik_ind)
         pssup=vik_tab(1,vik_ind)
      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,vik_ind)
         pssup=vik_tab(iloop,vik_ind)
      endif

C 6. Interpolate VL1 pressure data to 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      write(out,*)'CALC_FACTCOR: presvik=',presvik

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

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine nearsurfacenoise(ps,ps_noise_dev,ps_noise,
     &     temp_gcm_noise_dev,temp_gcm_noise,sheight)
      
! Routine to build perturbations (noise) for surface pressure and
! atmospheric temperatures near the surface. The amplitudes and
! range over which these perturbations are computed are hard coded
! in the routine (see below). E.M. 03/2009.
      implicit none
      include "constants_mcd.inc"
! arguments -inputs-:
      real ps ! surface pressure
      real ps_noise_dev ! gaussian deviate (of unit standard deviation)
!      real temp_gcm(dimlevs) ! atmospheric temperatures
      real temp_gcm_noise_dev ! gaussian deviate 
      real sheight(dimlevs) ! altitude of layers
! arguments -outputs-:
      real ps_noise ! perturbation to add to surface pressure
      real temp_gcm_noise(dimlevs) ! perturbations to add to temperatures

! local variables:
      real old_ps_noise_dev ! to store value from previous calls
      data old_ps_noise_dev /-77./ ! dummy initial value
      save old_ps_noise_dev 

      real ps_noise_level ! standard deviation of (relative) pert. in pressure
      parameter (ps_noise_level=0.01) !  e.g. 0.02=2%
      
      real T_noise_level ! standard deviation of pert in temperature (K)
      parameter (T_noise_level=3) ! e.g. 5K
      
      real T_noise_mid ! altitude (m) at which temp pert. will be half of max
      parameter (T_noise_mid=6000)
      
      real T_noise_delta ! range (m) over which noise drops from max to zero
      parameter (T_noise_delta=4000)

      integer lay
      real coeff ! altitude-dependent coefficient for temperature pert.
      
      if (ps_noise_dev.ne.old_ps_noise_dev) then ! compute perturbations
      
        ! 1. surface pressure perturbation
        ps_noise=ps*ps_noise_dev*ps_noise_level
!        write(*,*) 'ps_noise_dev=',ps_noise_dev,
!     &             ' temp_gcm_noise_dev=',temp_gcm_noise_dev
!        write(*,*) 'ps=',ps,
!     &             ' ps_noise=',ps_noise
        
        ! 2. atmospheric temperature perturbations
        do lay=1,dimlevs
          coeff=0.5*
     &          (1.-tanh(6.*(sheight(lay)-T_noise_mid)/T_noise_delta))
          temp_gcm_noise(lay)=T_noise_level*coeff*temp_gcm_noise_dev
!          write(*,*) 'l=',lay,'sheight(l)=',sheight(lay),
!     &                 'temp_gcm_noise(l)=',temp_gcm_noise(lay)
        enddo
        
        ! reset old_ps_moise_dev
        old_ps_noise_dev=ps_noise_dev
      endif ! of if (ps_noise_dev.ne.old_ps_noise_dev)
      
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine colint(aps,bps,oroheight,ps_gcm,ps_hr,R_gcm,temp_gcm,
     &                  col_gcm,col_hr)
c
c     Computes number of mol. between levels for hi and GCM resolution
c     Used for correction of column in hires. 
c
      implicit none

      include "constants_mcd.inc"
      
c     inputs
      real,intent(in) :: aps(dimlevs)   ! hybrid coordinate 
      real,intent(in) :: bps(dimlevs)   ! hybrid coordinate      
      real,intent(in) :: oroheight      !height of surface above reference areoid (m)
      real,intent(in) :: ps_gcm         ! surface pressure (Pa)
      real,intent(in) :: ps_hr          ! surface pressure (Pa)
      real,intent(in) :: R_gcm(dimlevs) ! R at GCM levels
      real,intent(in) :: temp_gcm(dimlevs) ! temperature at GCM levels

c     outputs
      real,intent(out) :: col_gcm(dimlevs) ! number of mol. between levels !! actual dimension used is: dimlevs-1 !!
      real,intent(out) :: col_hr(dimlevs)  ! number of mol. between levels !! actual dimension used is: dimlevs-1 !!      

c     local variables
      integer       l
c     Gravity on mean areoid
c     The areoid is defined as a surface of constant gravitational plus
c     rotational potential. The inertial rotation rate of Mars is assumed
c     to be 0.70882187E-4 rad/s. This potential is the mean value at the
c     equator at a radius of 3396.000 km, namely 12652804.7 m^2/s^2,
c     calculated from Goddard Mars Gravity Model mgm1025
c     [LEMOINEETAL2001] evaluated to degree and order 50.
      real          g, g0, a0
      parameter (a0=3396.E3,g0=3.7257964)

      real          Tmean            ! "mean" temperature of a layer
      real          Rogct
      real          sig_gcm(dimlevs) ! sigma levels
      real          sig_hr(dimlevs)  ! sigma levels
      real          x(dimlevs)       ! the stub for compatibility
      real          h_low,h_hi       ! height above surface     

      ! GCM sigma levels
      do l=1,dimlevs
        sig_gcm(l)=aps(l)/ps_gcm+bps(l)
      enddo
      ! high res sigma levels
      call build_sigma_hr(sig_gcm,ps_gcm,ps_hr,sig_hr,x)
c      
c     Calculate altitude above the surface of each model layer: h(l)
c     integrate hydrostatic equation
c
      g     = g0*(a0/(a0+oroheight))**2
      Rogct = R_gcm(1)/g
      h_low =-Rogct*temp_gcm(1)*log(sig_gcm(1))
      h_hi  =-Rogct*temp_gcm(1)*log(sig_hr(1))

      do l=2, dimlevs
       if (temp_gcm(l).ne.temp_gcm(l-1)) then
        Tmean = real((dble(temp_gcm(l)-temp_gcm(l-1))) /
     &                dlog(dble(temp_gcm(l))/dble(temp_gcm(l-1))))
       else
        Tmean = temp_gcm(l)
       end if
       ! For GCM 
       g = g0*(a0/(a0+oroheight+h_low))**2
       Rogct = R_gcm(l)/g
       h_low = h_low - Rogct*Tmean*log(sig_gcm(l)/sig_gcm(l-1))
       col_gcm(l-1)=-Rogct/1.380648e-23*(sig_gcm(l)-sig_gcm(l-1))*ps_gcm
       ! For hires
       g = g0*(a0/(a0+oroheight+h_hi))**2
       Rogct = R_gcm(l)/g
       h_hi = h_hi - Rogct*Tmean*log(sig_hr(l)/sig_hr(l-1))
       col_hr(l-1) =-Rogct/1.380648e-23*(sig_hr(l)-sig_hr(l-1))*ps_hr
      end do

      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine colcor(a,lon,lat,utime,name,ier,itimint,wl,wh,      
     &                  col_gcm,col_hr)
c
c     Calculate the ratio of column in hi res to column in GCM res. 
c     Used for correction of the column in hires.
c
      implicit none

      include "constants_mcd.inc"

c     inputs
      real,intent(in) :: lon     ! east longitude of point
      real,intent(in) :: lat     ! latitude of point
      real,intent(in) :: utime   ! Universal time (0. to 24. hrs)=local time at lon=0
      character(len=*),intent(in) :: name    ! Name of variable
      integer,intent(in) :: itimint ! seasonal interpolation flag (1==yes,0==no)
      real,intent(in) :: wl,wh   ! weights for seasonal interpolation      
      real,intent(in) :: col_gcm(dimlevs) !
      real,intent(in) :: col_hr(dimlevs) !

c     outputs
      real,intent(out) :: a      ! col_hr/col_gcm
      integer,intent(out) :: ier    ! error flag (0=OK, not 0 =NOK)

c     local variables
      integer       l
      real          x_gcm,x_hr,vmr,cor
      real          vmr_gcm(dimlevs)      
           
      call profi(vmr_gcm,lon,lat,utime,name,ier,itimint,wl,wh,0,0,0)
      
      x_gcm = 0.0
      x_hr  = 0.0
      do l=2,dimlevs
       vmr = 0.5*(vmr_gcm(l-1)+vmr_gcm(l))
       x_gcm = x_gcm + col_gcm(l-1)*vmr
       x_hr  = x_hr  + col_hr( l-1)*vmr
      enddo
      a = x_hr/x_gcm 
          
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine profi(a,lon,lat,utime,name,ier,itimint,wl,wh,flag_ini,
     &                 ib,ie)

c     Retrieve a vertical profile of 3-d variable=name at east longitude=lon,
c     latitude=lat and universal time=utime.
c     horizontal bilinear interpolation is used for all variables exept
c     density, for which the horizontal interpolation is bilinear
c     interpolation of log(rho)
c     Note: see routine loadvar_mcd.F about the way data arrays are organized.
! IMPORTRANT: order of variables here must match order in "loadvar_mcd"
c
c flag_ini = 0 : NO initialization, use previously computed values (if possible)
c flag_ini = 1 : DO initialization for ALL variables and exit
c flag_ini > 1 : DO initialization for particular variable and extract the profile
c 
c The profile is extracted for indexes from ib to ie.
c If ib=0 and ie=0 the complete profile is extracted (from 1 to dimlevs)
c For any other ie,ib the profile is extracted for 1<=index_begin<=dimlevs, 1<=index_end<=dimlevs.
c If ib>ie then error.
c If lon=var_lon AND lat=var_lon AND utime=var_tuime AND  var_ib<=ib AND ie<=var_ie  the stored profile is taken.
c
      implicit none

      include "constants_mcd.inc"

c     inputs
      real,intent(in) :: lon     ! east longitude of point
      real,intent(in) :: lat     ! latitude of point
      real,intent(in) :: utime   ! Universal time (0. to 24. hrs)=local time at lon=0
      character(len=*),intent(in) :: name    ! Name of variable
      integer,intent(in) :: itimint ! seasonal interpolation flag (1==yes,0==no)
      real,intent(in) :: wl,wh   ! weights for seasonal interpolation  
      integer,intent(in) :: flag_ini ! flag of initialization (in the case of new season or scenario)
                            ! 0= no initialization; 
                            ! 1= do initialization and exit; 
                            ! 1> do initialization and profile
      integer,intent(in) :: ib,ie   ! begin end indexes
c     outputs
      integer,intent(out) ::  ier        ! error flag (0=OK, 1=Not OK)
      real,intent(out) :: a(dimlevs) ! Interpolated profile 

c     local variables
      integer  i,j,l    ! for loops
c      integer  k,sd,rms,arms ! flags/index of variable
      integer  itime(2) ! nearest MCD time indexes (1-12)
      integer  iut      ! temporary variable to store itime()
      integer  dlon(4)  ! longitude indexes of nearby MCD grid points
      integer  dlat(4)  ! latitude indexes of nearby MCD grid points
      integer  idlon ! to store dlon(i) outside i loop to save time (?)
      integer  idlat ! to store dlat(i) outside i loop to save time (?)
      real     t,u ! weights for horizontal (bilinear) interpolation
      real     w ! weight for time-of-day interpolation
      real     y(2,dimlevs,4) ! low/hi, MCD profiles at nearby grid points
      real     x(2,dimlevs,2) ! low/hi, horizontally interpolated y() at times itime()
      real     b(2,dimlevs)
      
      integer  ivar,it,k,var_type
      integer  indxb,indxe ! begin, end index for the profile extraction
      integer,parameter :: nvar = nbvar3d+nbsd3d+(nbsd3d+1)
c var_a() is a profile at var_lon(),var_lat(), var_utime() for ivar
      real, save :: var_lon(nvar)      
      real, save :: var_lat(nvar)
      real, save :: var_utime(nvar)
      real, save :: var_a(nvar,dimlevs)
      integer, save :: var_ib(nvar)
      integer, save :: var_ie(nvar)    
      
      character*12 rmsname(nbsd3d)    /'rmstemp', 'rmsu', 'rmsv',
     &                                  'rmsrho', 'rmsw'/
      character*12 armsname(nbsd3d+1)/'armstemp','armsu','armsv',
     &                                 'armsrho','armsw','armspressure'/

      ier=0    ! initialize error flag to 0 (== OK)
c
c flag_ini=1 
c Initialization for ALL variables (e.g. in case of a new season or scenario)
c by some dummy "unrealistic" values
c 
      if(flag_ini.eq.1) then
       do ivar=1,nvar
        var_lon(ivar)   = 1000.0
        var_lat(ivar)   = 1000.0
        var_utime(ivar) = 1000.0
        var_ib(ivar)    =-1
        var_ie(ivar)    =-1
       enddo
       return 
      endif

c associate variable with name
c var_type: 0=sd; 1=rms; 2=arms
c k       : index in the array
c ivar    : index of the profile 
c
      do var_type=0,2
       if(var_type.eq.0) then
        do k=1,nbvar3d
         if(varname3d(k).eq.name) then          
          ivar = k
          goto 1111
         endif 
        enddo
       elseif(var_type.eq.1) then
        do k=1,nbsd3d
         if(rmsname(k).eq.name) then          
          ivar = nbvar3d+k
          goto 1111
         endif   
        enddo
       elseif(var_type.eq.2) then
        do k=1,(nbsd3d+1)
         if(armsname(k).eq.name) then          
          ivar = nbvar3d+nbsd3d+k
          goto 1111
         endif   
        enddo   
       endif    
      enddo

c     CASE of an unexpected name
      if(output_messages) then
       write(out,*) 'problem using subroutine profi : the name ',
     &                  trim(name)
       write(out,*) 'is not recognized'
      endif
      ier=1
      stop
      
1111  continue

c
c flag_ini>1
c Initialization for a particular variable and proceed (i.e. as in "old" version) 
c
      if(flag_ini.gt.1) then
       var_lon(ivar)   = 1000.0
       var_lat(ivar)   = 1000.0
       var_utime(ivar) = 1000.0
       var_ib(ivar)    =-1
       var_ie(ivar)    =-1
      endif
c
c if ib=0 and ie=0 extract full profile: from 1 to dimlevs
c else for any ie,ib: 1<=lb<=dimlevs, 1<=le<=dimlevs 
c if lb>le: error
c      
      if(ib.eq.0.and.ie.eq.0) then
       indxb = 1
       indxe = dimlevs
      else
       indxb = ib
       indxe = ie
       ! check beginning and ending index coherence
       if(indxb.lt.1.or.indxe.gt.dimlevs.or.indxb.gt.indxe) then
        if (output_messages) then
         write(out,*) 'profi() : problem with indexes ibegin=',indxb,
     &                ' iend=',indxe
        endif
        ier=1
        stop
       endif ! indxb>indxe 
      endif ! else if ib=0 and ie=0    
c          
c
c
      if(utime.eq.var_utime(ivar)) then
       if(lon.eq.var_lon(ivar)) then
        if(lat.eq.var_lat(ivar)) then
         if(indxb.ge.var_ib(ivar).and.indxe.le.var_ie(ivar)) then   
          do l=indxb,indxe
           a(l) = var_a(ivar,l)          
          enddo
          return
         endif
        endif    
       endif    
      endif    
     
c     Retrieving variable for  "lower" season
c     ------------------------------------------------
c     find the 4 neighbouring MCD grid points
      call grid4(lon,lat,dlon,dlat,t,u)
c     find nearest 2 timestep :
      call mcd_time(utime,itime,w) 
c
c it=1 Retrieving variable for "lower season"
c it=2 Retrieving variable for "higher" season (if seasonal interpolation)
c 
      do it=1,2
c     --------------------------------------
c     For seasonal interpolation
c     loop on the 2 nearest timestep :
      do j = 1 , 2
         iut = itime(j) ! MCD time
c        get 4 profiles at those points
c    MEAN variables
       if (var_type.eq.0) then
            do i=1,4
              idlon=dlon(i)
              idlat=dlat(i)
               do l=indxb,indxe
                 if(it.eq.1) then  
                  y(it,l,i)=var_3d(idlon,idlat,l,iut,k)
                 else
                  y(it,l,i)=var_3d2(idlon,idlat,l,iut,k)   
                 endif    
               enddo
            enddo
       endif!sd=0
c    RMS variables
       if (var_type.eq.1) then
           do i=1,4
             idlon=dlon(i)
             idlat=dlat(i)
              do l=indxb,indxe
               if(it.eq.1) then   
                y(it,l,i)=varrms3d(idlon,idlat,l,k)
               else
                y(it,l,i)=varrms3d2(idlon,idlat,l,k)   
               endif    
              enddo
           enddo
       endif
c    ARMS variables
       if (var_type.eq.2) then
           do i=1,4
             idlon=dlon(i)
             idlat=dlat(i)
              do l=indxb,indxe
               if(it.eq.1) then   
                y(it,l,i)=vararms3d(idlon,idlat,l,k)
               else
                y(it,l,i)=vararms3d2(idlon,idlat,l,k)   
               endif    
              enddo
           enddo
       endif

       if ((name.eq.'rho').or.
     &     (name.eq.'tsdrho').or.
     &     (name.eq.'rmsrho').or.
     &     (name.eq.'armsrho').or.
     &     (name.eq.'armspressure')) then
c    bilinear interpolation of log(rho)
         do l=indxb,indxe
          x(it,l,j)=(1.-t)*(1.-u)*log(y(it,l,1))+t*(1.-u)*log(y(it,l,2))
     &           +t*u*log(y(it,l,3))+(1.-t)*u*log(y(it,l,4))
          x(it,l,j)=exp(x(it,l,j))
         enddo
       else
c    bilinear interpolation
         do l=indxb,indxe
           x(it,l,j)=(1.-t)*(1.-u)*y(it,l,1)+t*(1.-u)*y(it,l,2)
     &           +t*u*y(it,l,3)+(1.-t)*u*y(it,l,4)
         enddo
       endif
      enddo! loop on 2 nearest MCD timestep

c    linear interpolation in (hour-of-day) time :
      do l=indxb,indxe
        b(it,l) = w*x(it,l,1) + (1-w)*x(it,l,2)
c      if(k.eq.4) write(out,*)'l,t ',l,a(l)
      enddo

      if(itimint.lt.1) exit ! if NO seasonal interpolation
      
      enddo ! it
c      
c     Season interpolation between the two seasons
c     --------------------------------------------
      if(itimint.lt.1) then      
       do l=indxb,indxe
        a(l) = b(1,l)
       enddo 
      else
       do l=indxb,indxe 
        a(l) = (wl)*b(1,l) + (wh)*b(2,l) 
       enddo
      endif     
c
c store the profile (from indxb to indxe)
c
      do l=indxb,indxe
       var_a(ivar,l) = a(l)   
      enddo    
      var_lat(ivar)   = lat
      var_lon(ivar)   = lon
      var_utime(ivar) = utime
      var_ib(ivar)    = indxb
      var_ie(ivar)    = indxe

      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
