Date: Mon, 24 Nov 2003 17:16:49 +0100
From: Jean-Jacques Valette <valette@cls.fr>
To: urs Hugentobler <urs.hugentobler@aiub.unibe.ch>
Cc: petr.stepanek@cls.fr, martine Feissel-Vernier <feissel@ensg.ign.fr>,
    Karine Le Bail <lebail@ensg.ign.fr>,
    gilles Tavernier <Gilles.Tavernier@cnes.fr>,
    jean-pierre.granier@cnes.fr, laurent Soudarin <Laurent.Soudarin@cls.fr>,
    Biancale Richard <Richard.Biancale@cnes.fr>
Subject: DORIS and DTM1994

Dear Urs and Petr,

Following our DORIS meeting in Berne, I have asked to Richard about the DTM 
air Drag model applied in the GINS software.

The reference of the original work from F. Barlier et al. is:
F. Barlier, C. Berger, J.L. Falin, G. Kokckarts and G. Thuillier. A 
Thermospheric model based on the satellite drag data. Ann. Geophys., t. 34, 
1 , 1978, p 9-24.

You will find the fortran 90 code in the attached file. It is the one 
Laurent routinely uses for the LEGOS/CLS analysis center.

Richard mentioned some recent refinements of this model, so-called DTM2000 
based on the CHAMP accelorometer data. It is also available and may be 
asked to Richard or Sean Bruinsma (Sean.Bruinsma@cnes.fr).

Many thanks to Richard and GRGS team and hope it will help in your analysis.

Best regards,

JJacques
------------------------------------------------------------------------------
    Jean-Jacques VALETTE (Jean-Jacques.Valette@cls.fr)
    COLLECTE LOCALISATION SATELLITES
    8-10 Rue Hermès
    Parc Technologique du Canal
    31526 Ramonville-Saint-Agne Cedex
    FRANCE
    Tel (+33) 05.61.39.47.62
    Fax (+33) 05.61.39.48.06
    <http://www.cls.fr>http://www.cls.fr
-------------------------------------------------------------------------------
  [ Part 2: "Attached Text" ]


subroutine dtm94(day,f,fbar,ap,alti,hl,alat,xlon,tz,tinf,ro,wmm,ddh)
  !
  !***********************************************************************
  !
  ! HISTORIQUE
  ! ----------
  !
  ! ROLE
  ! ----
  !*auteur
  !*version decembre 94
  !*but calcul de la temperature, concentrations, densite totale
  !*par entrees
  !     day=jour de l annee
  !     f=flux instantane a t - 1j
  !     fbar=flux moyen a t
  !     akp= kp tri-horaire (avec un delai de 6-abs(alat)*0.033 en heure)
  !     alti=altitude en km superieure a 120 km
  !     hl=heure locale + 12h (en radian)
  !     aphi= latitude du point en radian
  !     alon= longitude du point en radian
  !*par sorties
  !     tz=temperature a l altitude alti
  !     tinf=temperature exospherique
  !     dbc(1)=concentration hydrogene atomique
  !     dbc(2)=concentration helium
  !     dbc(3)=concentration oxygene atomique
  !     dbc(4)=concentration azote moleculaire
  !     dbc(5)=concentration oxygene moleculaire
  !     ro=densite totale en g/cm3
  !     wmm=masse moleculaire moyenne
  !     ddh=derivee de la densite par rapport a l'altitude en km
  !
  ! ARGUMENTS
  ! ---------
  !*ARG  day   E  jour de l annee
  !*ARG  f     E  flux instantane a t - 1j et derive
  !*ARG  fbar  E  flux moyen a t et derive
  !*ARG  ap    E? ap moyen (DTM), a t, -1j et derivees
  !*ARG  alti  E  altitude en km superieue a 120 km
  !*ARG  hl    E  heure locale + 12h
  !*ARG  alat  E  latitude du point en radian
  !*ARG  xlon  E  longitude du point en radian
  !*ARG  tz    S  temperature a l altitude alti
  !*ARG  tinf  S  temperature exospherique
  !*ARG  ro    S  densite totale en g/cm3
  !*ARG  wmm   S  masse moleculaire moyenne
  !*ARG  ddh   S  derivee de la densite par rapport a l'altitude en km
  !      
  ! MODULES
  ! -------
  !
  ! APPELES
  ! -------
  !
  ! APPELANTS
  ! ---------
  !
  !***********************************************************************
  !
  use f90_kind
  use common_plgdtm
  use common_hlocal
  use common_cons
  use gldtm94_interface
  use bint_interface
  ! pour inline bint
  use egalite_interface
  ! 
  !.. Implicit Declarations .. 
  implicit none
  ! 
  !.. Parameters .. 
  integer, parameter :: nlatm94 = 39
  ! 
  !.. Formal Arguments .. 
  real(SINGLE), intent(in) :: day
  real(SINGLE), intent(in) :: f
  real(SINGLE), intent(in) :: fbar
  real(SINGLE), dimension(7), intent(in) :: ap
  real(SINGLE), intent(in) :: alti
  real(SINGLE), intent(in) :: hl
  real(SINGLE), intent(in) :: alat
  real(SINGLE), intent(in) :: xlon
  real(SINGLE), intent(out) :: tz
  real(SINGLE), intent(out) :: tinf
  real(SINGLE), intent(out) :: ro
  real(SINGLE), intent(out) :: wmm
  real(SINGLE), intent(out) :: ddh
  ! 
  !.. Local Scalars .. 
  integer :: i,ikp
  real(SINGLE) :: cpmg = .19081_S
  real(SINGLE) :: gsurf = 980.665_S
  real(SINGLE) :: re = 6356.77_S
  real(SINGLE) :: rgas = 831.4_S
  real(SINGLE) :: spmg = .98163_S
  real(SINGLE) :: xlmg = -1.2392_S
  real(SINGLE) :: zlb = 120._S
  real(SINGLE) :: akp,apt,c,c2,c4,clmlmg,cmg,cmg2,cmg4,dzeta,expsz,fbm150, &
                  fmfb,gamma,gdelaz2,gdelh,gdelhe,gdelo,gdelo2,gdelt,glb,sbc,s2, &
                  sigma,sigzeta,sp,t120,t120tz,tinftz,tp120,upapg,zeta
  ! 
  !.. Local Arrays .. 
  integer, dimension(6) :: ma
  real(SINGLE), dimension(6) :: alefa,cc,dbc,dbase,fz,vma
  real(SINGLE), dimension(nlatm94) :: &
    az2,daz2,dh,dhe,do2,dtt,h,he,o,o2,t0,tp,tt,xdo
  ! 
  !.. Intrinsic Functions .. 
  intrinsic abs, asin, cos, exp, sin
  ! 
  !.. Data Declarations .. 
  data alefa/-0.40_S,-0.38_S,0._S,0._S,0._S,0._S/
  data ma/1,4,16,28,32,14/
  data vma/ &
       1.6606e-24_S,6.6423e-24_S,26.569e-24_S,46.4958e-24_S,53.1381e-24_S,23.2479e-24_S/
  data (tt(i), i = 1,39)/ &
       1000.5_S,0.94610e-02_S,0.42671e-01_S,0.17948e-02_S,-0.79900e-05_S,0.33674e-02_S, &
       0.22629e-01_S,0.37865e-01_S,-0.19230e-01_S,-0.92415e-02_S,-0.21066e03_S, &
       0.10326e-01_S,0.28856e-01_S,-0.76311e02_S,-0.18498e00_S,-0.20306e-01_S, &
       0.14473e-01_S,-0.36306e01_S,-0.28940e-01_S,-0.17367e03_S,-0.10819e00_S, &
       -0.19993e-02_S,0.33966e-02_S,-0.16132e-01_S,-0.96095e-02_S,-0.10455e00_S, &
       0.45752e-02_S,0.45844e-02_S,0.17764e-01_S,-0.42270e-02_S,-0.35667e-02_S, &
       -0.36061e-03_S,0.10493e-01_S,0.45712e-02_S,-0.21751e-04_S,0.15110e-02_S, &
       0.21666e-02_S,-0.11423e-05_S,0.65823e-04_S/
  data (h(i), i = 1,39)/ &
       1.761e+05_S,-1.33700e-01_S,0._S,-1.24600e-02_S,0._S,-1.93000e-02_S,-6.00000e-02_S, &
       -0.20000e-02_S,5.87800e-02_S,0._S,9.22700e01_S,0._S,0._S,0._S,3.30100e-01_S, &
       1.04500e-01_S,0._S,-1.47700e01_S,-9.06500e-02_S,-7.20000e01_S,2.09400e-01_S, &
       2.83000e-02_S,0._S,8.57100e-02_S,-2.47500e-02_S,3.83000e-01_S,2.94100e-02_S,0._S, &
       -3.97400e-03_S,4.35600e-02_S,0._S,0._S,0._S,0._S,0._S,0._S,0._S,0._S,0._S/
  data (he(i), i = 1,39)/ &
       2.791e+07_S,0.10965e00_S,-0.19084e00_S,-0.20772e-03_S,0.48346e-05_S,0.21123e-02_S, &
       0.22119e-03_S,-0.16174e00_S,-0.92223e-01_S,-0.83007e-02_S,0.21353e03_S, &
       0.23504e00_S,-0.79051e-01_S,0.11040e03_S,-0.12677e01_S,-0.18512e-01_S, &
       0.66625e-01_S,-0.18703e03_S,-0.42149e-01_S,-0.21671e03_S,-0.12779e00_S, &
       -0.61824e-02_S,-0.17450e-01_S,-0.43600e-01_S,-0.52850e-01_S,0.31204e00_S, &
       -0.21372e-01_S,-0.24529e-01_S,-0.36727e-02_S,-0.88330e-01_S,0.33991e-01_S, &
       0.51314e-02_S,-0.14740e-01_S,0.85563e-02_S,0.19571e-02_S,-0.45034e-02_S, &
       -0.10241e00_S,-0.16867e-04_S,-0.14265e-02_S/
  data (o(i), i = 1,39)/ &
       0.8472e+11_S,-0.66447e-01_S,-0.97415e-01_S,0.12284e-02_S,0.44976e-05_S, &
       0.53580e-02_S,0.25573e-02_S,-0.98221e-01_S,0.10091e00_S,0.62609e-02_S,0.11565e02_S, &
       0.17600e00_S,-0.71284e-01_S,0.10639e03_S,0.33295e00_S,-0.11448e00_S,-0.42425e-02_S, &
       -0.67511e00_S,-0.41712e-01_S,0.13443e03_S,-0.65932e-01_S,-0.19344e-01_S, &
       -0.85846e-02_S,0.86347e-01_S,0.85679e-01_S,0.45486e-01_S,-0.38719e-01_S, &
       -0.12872e-01_S,-0.91778e-01_S,0.42617e-01_S,0.43386e-01_S,0.56894e-02_S, &
       0.66512e-02_S,-0.10174e-01_S,0.27058e-02_S,-0.47104e-02_S,-0.14639e-01_S, &
       -0.64057e-05_S,0.15178e-02_S/
  data (az2(i), i = 1,39)/ &
       3.2045e+11_S,-0.14019e00_S,0.57220e-01_S,0.11260e-02_S,-0.20779e-05_S, &
       0.34239e-02_S,-0.11588e-01_S,+0.55158e-01_S,-0.10049e-01_S,-0.43905e-01_S, &
       0.19594e03_S,0.32483e-01_S,0.56500e-01_S,0.88202e02_S,0.28812e00_S,-0.31427e-01_S, &
       0.00000e00_S,-0.20015e03_S,0.60700e-01_S,0.51153e02_S,-0.46104e-01_S, &
       -0.97000e-02_S,0.34910e-02_S,0._S,0._S,-0.73329e-01_S,0.21159e-01_S,0.70334e-02_S,0._S, &
       0._S,-0.69199e-02_S,0._S,-0.72899e-02_S,0._S,-0.29200e-02_S,0.27047e-02_S, &
       -0.59939e-02_S,0.00000e00_S,0.82127e-02_S/
  data o2(1)/4.775e+10_S/
  data t0(1)/380.0_S/
  data tp(1)/14.348_S/
  ! 
  ! ... Executable Statements ...
  ! 
  !
  ddh = 0._S
  ro = 0._S
  !
  !   calcul des polynomes de legendre
  c = sin(alat)
  c2 = c * c
  c4 = c2 * c2
  sbc = cos(alat)
  s2 = sbc * sbc
  p10 = c
  p20 = 1.5_S*c2 - 0.5_S
  p30 = c * (2.5_S*c2-1.5_S)
  p40 = 4.375_S*c4 - 3.75_S*c2 + 0.375_S
  p50 = c * (7.875_S*c4-8.75_S*c2+1.875_S)
  p11 = sbc
  p21 = 3._S * c * sbc
  p31 = sbc * (7.5_S*c2-1.5_S)
  p41 = c * sbc * (17.5_S*c2-7.5_S)
  p51 = sbc * (39.375_S*c4-26.25_S*c2+1.875_S)
  p22 = 3._S * s2
  p32 = 15._S * c * s2
  p42 = s2 * (52.5_S*c2-7.5_S)
  p52 = 3._S*c*p42 - 2._S*p32
  p33 = 15._S * sbc * s2
  !
  !   calcul des polynomes de legendre / pole magnetique (79n,71w)
  clmlmg = cos(xlon-xlmg)
  sp = sbc*cpmg*clmlmg + c*spmg
  if (sp > 1._S) then
    sp = 1.0_S
  end if
  if (sp < -1._S) then
    sp = -1._S
  end if
  !
  cmg = sp ! pole magnetique
  cmg2 = cmg * cmg
  cmg4 = cmg2 * cmg2
  p10mg = cmg
  p20mg = 1.5_S*cmg2 - 0.5_S
  p40mg = 4.375_S*cmg4 - 3.75_S*cmg2 + 0.375_S
  !
  !   heure locale
  hl0 = hl + pi
  ch = cos(hl0)
  sh = sin(hl0)
  c2h = ch*ch - sh*sh
  s2h = 2._S * ch * sh
  c3h = c2h*ch - s2h*sh
  s3h = s2h*ch + c2h*sh
  !
  !   flux
  fmfb = f - fbar
  fbm150 = fbar - 150._S
  !
  !   kp
  apt = ap(4) + 2._S*abs(asin(cmg))/pi*(ap(3)-ap(4))
  akp = bint(apt)
  !
  !   calcul de la fonction g(l) / tinf, t120, tp120
  ikp = 1
  call gldtm94(fmfb,fbm150,akp,day,tt,gdelt,1._S,ikp)
  dtt(1) = 1._S + gdelt
  tinf = tt(1) * dtt(1)
  t120 = t0(1)
  tp120 = tp(1)
  !  Attention BUG: pour certaines configurations tinf est inferieur
  !  a t120 ce qui est impossible.
  !  on impose donc tinf =t120+1 si tinf <ou= t120
  if (tinf <= t120) then
    tinf = t120 + 1._S
  end if
  !
  !   calcul des concentrations n(z): h, he, o, n2, o2
  sigma = tp120 / (tinf-t120)
  dzeta = (re+zlb) / (re+alti)
  zeta = (alti-zlb) * dzeta
  sigzeta = sigma * zeta
  expsz = exp(-sigzeta)
  tz = tinf - (tinf-t120)*expsz
  !
  !   calcul de la fonction g(l) / h, he, o, n2, o2
  ikp = 2
  call gldtm94(fmfb,fbm150,akp,day,h,gdelh,0._S,ikp)
  dh(1) = exp(gdelh)
  dbase(1) = h(1) * dh(1)
  call gldtm94(fmfb,fbm150,akp,day,he,gdelhe,0._S,ikp)
  dhe(1) = exp(gdelhe)
  dbase(2) = he(1) * dhe(1)
  call gldtm94(fmfb,fbm150,akp,day,o,gdelo,1._S,ikp)
  xdo(1) = exp(gdelo)
  dbase(3) = o(1) * xdo(1)
  call gldtm94(fmfb,fbm150,akp,day,az2,gdelaz2,1._S,ikp)
  daz2(1) = exp(gdelaz2)
  dbase(4) = az2(1) * daz2(1)
  call gldtm94(fmfb,fbm150,akp,day,o2,gdelo2,1._S,ikp)
  do2(1) = exp(gdelo2)
  dbase(5) = o2(1) * do2(1)
  !
  glb = gsurf / (1._S+zlb/re)**2
  glb = glb / (sigma*rgas*tinf)
  t120tz = t120 / tz
  tinftz = tinf / tz
  do i = 1,5
    gamma = REAL(ma(i),SINGLE) * glb
    upapg = 1._S + alefa(i) + gamma
    fz(i) = t120tz**upapg * exp(-sigzeta*gamma)
    !   concentrations en h, he, o, n2, o2
    cc(i) = dbase(i) * fz(i)
    !   densites en h, he, o, n2, o2
    dbc(i) = cc(i) * vma(i)
    !
    !   densite totale
    ro = ro + dbc(i)
    !   derivee partielle / altitude (km)
    ddh = ddh + dbc(i)*sigma*dzeta*(1._S+alefa(i)-upapg*tinftz)
  end do
  !
  !   masse moleculaire moyenne
  !
  wmm = ro / (vma(1)*(cc(1)+cc(2)+cc(3)+cc(4)+cc(5)))
end subroutine dtm94 

 module common_pardtm 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! modif: 02/2001 s. bruinsma 
   ! 
   ! ROLE 
   ! ---- 
   ! 
   ! VARIABLES PUBLIQUES 
   ! ------------------- 
   !*VAR  ar       
   !*VAR  az       coefficients du calcul de g(l) de dtm pour l'azote 
   !*VAR  az2      coefficients du calcul de g(l) de dtm pour l'azote moleculaire 
   !*VAR  co       2      
   !*VAR  co2      
   !*VAR  h        2       
   !*VAR  h2       
   !*VAR  he       coefficients du calcul de g(l) de dtm pour l'hydrogene moleculaire 
   !*VAR  o        coefficients du calcul de g(l) de dtm pour l'oxygene 
   !*VAR  o2       coefficients du calcul de g(l) de dtm pour l'oxygene moleculaire 
   !*VAR  t0       coefficients du calcul de g(l) de dtm pour la temperature(z0) 
   !*VAR  t0_dtm   
   !*VAR  tp       coefficients du calcul de g(l) de dtm pour le gradient de temperature(z0)  
   !*VAR  tt       coefficients du calcul de g(l) de dtm pour la temperature exospherique   
   !*VAR  tau      opacite (IR) max. de l'atmosphere martienne pendant tempete 
   !*VAR  xlsdeb   Ls du debut de la tempete   (deg) 
   !*VAR  xlsfin   Ls de la fin de la tempete  (deg) 
   !*VAR  taubg    opacite de l'atmosphere martienne moyenne (background dust opacity) 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use parameters_divers 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   real(SINGLE), dimension(nlatm) :: az,az2,h,he,o,o2,t0_dtm,tp,tt,t0,co2,co, & 
     ar,h2 
   real(SINGLE) :: tau,taubg,xlsdeb,xlsfin 
 end module common_pardtm 

 module common_plgdtm 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   ! 
   ! VARIABLES PUBLIQUES 
   ! ------------------- 
   !*VAR  p10     A COMPLETER 
   !*VAR  p10mg   A COMPLETER 
   !*VAR  p11     A COMPLETER 
   !*VAR  p20     A COMPLETER 
   !*VAR  p20mg   A COMPLETER 
   !*VAR  p21     A COMPLETER 
   !*VAR  p22     A COMPLETER 
   !*VAR  p30     A COMPLETER 
   !*VAR  p31     A COMPLETER 
   !*VAR  p32     A COMPLETER 
   !*VAR  p33     A COMPLETER 
   !*VAR  p40     A COMPLETER 
   !*VAR  p40mg   A COMPLETER 
   !*VAR  p41     A COMPLETER 
   !*VAR  p42     A COMPLETER 
   !*VAR  p50     A COMPLETER 
   !*VAR  p51     A COMPLETER 
   !*VAR  p52     A COMPLETER 
   !*VAR  p60     A COMPLETER 
   !*VAR  p62     A COMPLETER 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   real(SINGLE) :: p10,p10mg,p11,p20,p20mg,p21,p22,p30,p31,p32,p33,p40,p40mg, & 
                   p41,p42,p50,p51,p52,p60,p62 
 end module common_plgdtm 

  
 subroutine dtm(day,f,fbar,ap,alti,hl,alat,xlon,tz,tinf,ro,dpra,wmm,ddh) 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*aut R.Biancale 
   !*ver dec 1995 
   !*rol calcul de la temperature,densite totale 
   !*par entrees 
   !     day=jour de l annee 
   !     f=flux instantane a t - 1j et derive 
   !     fbar=flux moyen a t et derive 
   !     ap(4)= ap moyen (DTM), a t, -1j et derivees 
   !     alti=altitude en km superieure a 120 km 
   !     hl=heure locale + 12h 
   !     alat= latitude du point en radian 
   !     xlon= longitude du point en radian 
   !*par sorties 
   !     tz=temperature a l altitude alti 
   !     tinf=temperature exospherique 
   !     dbc(1)=concentration hydrogene atomique 
   !     dbc(2)=concentration helium 
   !     dbc(3)=concentration oxygene atomique 
   !     dbc(4)=concentration azote moleculaire 
   !     dbc(5)=concentration oxygene moleculaire 
   !     dbc(6)=concentration azote atomique 
   !     ro=densite totale en g/cm3 
   !     wmm=masse moleculaire moyenne 
   !     ddh=derivee de la densite par rapport a l'altitude en km 
   !     dpra(..,9)=derivees partielles de la densite / parametres de la 
   !         temperature de la thermopause et des concentrations en h, 
   !         HE, O, N2, O2, T120, TP120 
   ! 
   ! ARGUMENTS 
   ! --------- 
   !*ARG  day   E  jour de l annee 
   !*ARG  f     E? flux instantane a t - 1j et derive 
   !*ARG  fbar  E? flux moyen a t et derive 
   !*ARG  ap    E? ap moyen (DTM), a t, -1j et derivees 
   !*ARG  alti  E  altitude en km superieure a 120 km 
   !*ARG  hl    E  heure locale + 12h 
   !*ARG  alat  E  latitude du point en radian 
   !*ARG  xlon  E  longitude du point en radian 
   !*ARG  tz    S  temperature a l altitude alti 
   !*ARG  tinf  S  temperature exospherique 
   !*ARG  ro    S  densite totale en g/cm3 
   !*ARG  dpra  S? derivee de la densite par rapport a l'altitude en km 
   !*ARG  wmm   S  masse moleculaire moyenne 
   !*ARG  ddh   S  derivee de la densite par rapport a l'altitude en km 
   ! 
   ! MODULES 
   ! ------- 
   !*MOD  datmo    E   npara 
   !               E?  libatm 
   !*MOD  hlocal   S   c2h, c3h, ch, hl0, s2h, s3h, sh 
   !*MOD  pardtm   E?  az, az2, h, he, o, o2, t0_dtm, tp, tt 
   !*MOD  plgdtm   S   p10, p10mg, p11, p20, p20mg, p21, p22, p30, p31, p32, 
   !               S   p33, p40, p40mg, p41, p42, p50, p51, p52 
   !*MOD  pourzpf  E   natmo 
   !               NU  coeff_corr, nasp, nb_lib_stoc, nbaro, ncdga, 
   !               NU  ncdgc, ncdgd, ncdgl, nfroto, ngrav, ngravd, nmar, 
   !               NU  nnbias, npolh, npolk, npolt, npolx, npoly, 
   !               NU  nporto, npreso, nrbias, nredif, nst, nst_tra, 
   !               NU  ntbias, ntherm, ntopo, nxbias, nybias, nzbias, 
   !               NU  tfrot, tnbias, tpres, trbias, ttbias, txbias, 
   !               NU  tybias, tzbias 
   ! 
   ! APPELES 
   ! ------- 
   ! bint 
   ! egalite 
   ! gldtm 
   ! 
   ! APPELANTS 
   ! --------- 
   ! frotem 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use parameters_divers 
   use common_pourzpf 
   use common_plgdtm 
   use common_pardtm 
   use common_hlocal 
   use common_datmo 
   use common_cons 
   use common_fic 
   use gldtm_interface 
   use bint_interface 
   use egalite_interface 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
   ! 
   !.. Formal Arguments .. 
   real(SINGLE), intent(in) :: day 
   real(SINGLE), dimension(2), intent(inout) :: f 
   real(SINGLE), dimension(2), intent(in) :: fbar 
   real(SINGLE), dimension(7), intent(in) :: ap 
   real(SINGLE), intent(in) :: alti 
   real(SINGLE), intent(in) :: hl 
   real(SINGLE), intent(in) :: alat 
   real(SINGLE), intent(in) :: xlon 
   real(SINGLE), intent(out) :: tz 
   real(SINGLE), intent(out) :: tinf 
   real(SINGLE), intent(out) :: ro 
   real(SINGLE), dimension(nlatm,9), intent(out) :: dpra 
   real(SINGLE), intent(out) :: wmm 
   real(SINGLE), intent(out) :: ddh 
   ! 
   !.. Local Scalars .. 
   integer :: i,kleq 
   real(SINGLE) :: cpmg = .19081_S 
   real(SINGLE) :: gsurf = 980.665_S 
   real(SINGLE) :: re = 6356.77_S 
   real(SINGLE) :: rgas = 831.4_S 
   real(SINGLE) :: spmg = .98163_S 
   real(SINGLE) :: xlmg = -1.2392_S 
   real(SINGLE) :: zlb = 120._S 
   real(SINGLE) :: akp2,akp4,c,c2,c4,clmlmg,cmg,cmg2,cmg4,dt120,dtinf,dtp120, & 
                   dzeta,dzeta2,expsz,gamma,gdelaz,gdelaz2,gdelh,gdelhe,gdelo, & 
                   gdelo2,gdelt,gdelt0,gdeltp,glb,sbc,s2,sigma,sigzeta,sp,t120, & 
                   t120tt,t120tz,tinftz,tp120,upapg,xlog,zeta 
   ! 
   !.. Local Arrays .. 
   integer, dimension(6) :: ma 
   real(SINGLE), dimension(4) :: akp 
   real(SINGLE), dimension(6) :: alefa,cc,dbc,dbase,fz,vma 
   real(SINGLE), dimension(nlatm) :: daz,daz2,dh,dhe,do2,dt0,dtp,dtt,xdo 
   ! 
   !.. Intrinsic Functions .. 
   intrinsic cos, exp, sin 
   ! 
   !.. Data Declarations .. 
   data alefa/-0.40_S,-0.38_S,0._S,0._S,0._S,0._S/ 
   data ma/1,4,16,28,32,14/ 
   data vma/ & 
        1.6606e-24_S,6.6423e-24_S,26.569e-24_S,46.4958e-24_S,53.1381e-24_S,23.2479e-24_S/ 
   ! 
   ! ... Executable Statements ... 
   ! 
   ! 
   do i = 1,nlatm 
     dtt(i) = 0._S 
     dh(i) = 0._S 
     dhe(i) = 0._S 
     xdo(i) = 0._S 
     daz2(i) = 0._S 
     do2(i) = 0._S 
     daz(i) = 0._S 
     dt0(i) = 0._S 
     dtp(i) = 0._S 
   end do 
   ro = 0._S 
   ddh = 0._S 
   dtinf = 0._S 
   dt120 = 0._S 
   dtp120 = 0._S 
   ! 
   !   calcul des polynomes de legendre 
   c = sin(alat) 
   c2 = c * c 
   c4 = c2 * c2 
   sbc = cos(alat) 
   s2 = sbc * sbc 
   p10 = c 
   p20 = 1.5_S*c2 - 0.5_S 
   p30 = c * (2.5_S*c2-1.5_S) 
   p40 = 4.375_S*c4 - 3.75_S*c2 + 0.375_S 
   p50 = c * (7.875_S*c4-8.75_S*c2+1.875_S) 
   p60 = (5.5_S*c*p50 - 2.5_S*p40)/3.0_S 
   p11 = sbc 
   p21 = 3._S * c * sbc 
   p31 = sbc * (7.5_S*c2-1.5_S) 
   p41 = c * sbc * (17.5_S*c2-7.5_S) 
   p51 = sbc * (39.375_S*c4-26.25_S*c2+1.875_S) 
   p22 = 3._S * s2 
   p32 = 15._S * c * s2 
   p42 = s2 * (52.5_S*c2-7.5_S) 
   p52 = 3._S*c*p42 - 2._S*p32 
   p62 = 2.75_S*c*p52 - 1.75_S*p42 
   p33 = 15._S * sbc * s2 
   ! 
   !   calcul des polynomes de legendre / pole magnetique (79n,71w) 
   clmlmg = cos(xlon-xlmg) 
   sp = sbc*cpmg*clmlmg + c*spmg 
   ! 
   cmg = sp ! pole magnetique 
   cmg2 = cmg * cmg 
   cmg4 = cmg2 * cmg2 
   !     smg=sqrt(cpcl*cpcl+cpsl*cpsl) 
   !     smg2=smg*smg 
   p10mg = cmg 
   p20mg = 1.5_S*cmg2 - 0.5_S 
   p40mg = 4.375_S*cmg4 - 3.75_S*cmg2 + 0.375_S 
   ! 
   !   heure locale (hl0=hl+pi) 
   hl0 = hl + pi 
   ch = cos(hl0) 
   sh = sin(hl0) 
   c2h = ch*ch - sh*sh 
   s2h = 2._S * ch * sh 
   c3h = c2h*ch - s2h*sh 
   s3h = s2h*ch + c2h*sh 
   ! 
   ! modif slb: il n'y a pas toujours le flux instantane dans acsol2 
   if(f(1) < 1.0_S) f(1)=fbar(1) 
   ! 
   !   kp 
   akp(1) = bint(ap(3)) ! a -3h 
   akp2 = bint(ap(4)) ! a -6h 
   !     calcul de la variation horaire 
   akp(2) = (akp(1)-akp2) * .333333333_S 
   ! 
   !   kp moyen a -1j 
   akp(3) = bint(ap(6)) ! moyen a -1j 
   akp4 = bint(ap(7)) ! a -2j 
   !     calcul de la variation moyenne journaliere 
   akp(4) = akp(3) - akp4 
   ! 
   !   calcul de la fonction g(l) / tinf, t120, tp120 
   kleq=1 
   call gldtm(f,fbar,akp,day,tt,dtt,gdelt,1._S,xlon,kleq) 
   dtt(1) = 1._S + gdelt 
   tinf = tt(1) * dtt(1) 
  
   kleq=0                                        !conditions d'equinox 
   if(day < 59.0_S .or. day   284.0_S) kleq=-1   !             hiver 
   if(day   99.0_S .and. day < 244.0_S) kleq= 1  !             ete 
   call gldtm(f,fbar,akp,day,t0_dtm,dt0,gdelt0,0._S,xlon,kleq) 
   dt0(1) = (t0_dtm(1) + gdelt0)/t0_dtm(1) 
   t120 = t0_dtm(1) + gdelt0 
   call gldtm(f,fbar,akp,day,tp,dtp,gdeltp,0.0_S,xlon,kleq) 
   dtp(1) = (tp(1) + gdeltp)/tp(1) 
   tp120 = tp(1) + gdeltp 
   !  Attention BUG: pour certaines configurations tinf est inferieur 
   !  a t120 ce qui est impossible. 
   !  on impose donc tinf =t120+1 si tinf <ou= t120 
   if (tinf <= t120) then 
     tinf = t120 + 1._S 
   end if 
   ! 
   !   calcul des concentrations n(z): H, HE, O, N2, O2, N 
   sigma = tp120 / (tinf-t120) 
   dzeta = (re+zlb) / (re+alti) 
   zeta = (alti-zlb) * dzeta 
   dzeta2 = dzeta * dzeta 
   sigzeta = sigma * zeta 
   expsz = exp(-sigzeta) 
   tz = tinf - (tinf-t120)*expsz 
   kleq=1 
   ! 
   !   calcul de la fonction g(l) / H, HE, O, N2, O2, N 
   call gldtm(f,fbar,akp,day,h,dh,gdelh,0._S,xlon,kleq) 
   dh(1) = exp(gdelh) 
   dbase(1) = h(1) * dh(1) 
   call gldtm(f,fbar,akp,day,he,dhe,gdelhe,0._S,xlon,kleq) 
   dhe(1) = exp(gdelhe) 
   dbase(2) = he(1) * dhe(1) 
   call gldtm(f,fbar,akp,day,o,xdo,gdelo,1._S,xlon,kleq) 
   xdo(1) = exp(gdelo) 
   dbase(3) = o(1) * xdo(1) 
   call gldtm(f,fbar,akp,day,az2,daz2,gdelaz2,1._S,xlon,kleq) 
   daz2(1) = exp(gdelaz2) 
   dbase(4) = az2(1) * daz2(1) 
   call gldtm(f,fbar,akp,day,o2,do2,gdelo2,1._S,xlon,kleq) 
   do2(1) = exp(gdelo2) 
   dbase(5) = o2(1) * do2(1) 
   if (egalite(az(1),0._S)) then 
     daz(1) = 1._S 
     dbase(6) = 0._S 
   else 
     call gldtm(f,fbar,akp,day,az,daz,gdelaz,1._S,xlon,kleq) 
     daz(1) = exp(gdelaz) 
     dbase(6) = az(1) * daz(1) 
   end if 
   ! 
   glb = gsurf / (1._S+zlb/re)**2 
   glb = glb / (sigma*rgas*tinf) 
   t120tz = t120 / tz 
 !CS 
   if(tp120   21.0_S .or. tp120 < 6.0_S) then 
     print*, gdelt0,dt0(1),t0_dtm(1),gdeltp 
     print*, tinf,tz,t120,tp120 
   endif 
   xlog = LOG(t120tz) 
   tinftz = tinf / tz 
   t120tt = t120 / (tinf-t120) 
   do i = 1,6 
     gamma = REAL(ma(i),SINGLE) * glb 
     upapg = 1._S + alefa(i) + gamma 
     fz(i) = t120tz**upapg * exp(-sigzeta*gamma) 
     !   concentrations en H, HE, O, N2, O2, N 
     cc(i) = dbase(i) * fz(i) 
     !   densites en H, HE, O, N2, O2, N 
     dbc(i) = cc(i) * vma(i) 
     ! 
     !   densite totale 
     ro = ro + dbc(i) 
     !   derivee partielle / altitude (km) 
     ddh = ddh + dbc(i)*sigma*dzeta2*(1._S+alefa(i)-upapg*tinftz) 
     !   derivee partielle / tinf, t120, tp120 
     if (npara/=0 .and. natmo/=0) then 
       ! pour iter. sup. 
       dtinf = dtinf + & 
               dbc(i)* & 
               (gamma*(sigzeta+t120tt*xlog)+ & 
                upapg*tinftz*(expsz*(1._S+sigzeta)-1._S))/tinf 
       dt120 = dt120 + & 
               dbc(i)*(upapg*(1._S-t120tz*expsz*(1._S+sigzeta))-gamma*t120tt-xlog)/ & 
               t120 
       dtp120 = dtp120 - dbc(i)*(zeta*upapg*expsz/tz+gamma*xlog/tp120) 
     end if 
   end do 
   ! 
   !   masse moleculaire moyenne 
   wmm = ro / (vma(1)*(cc(1)+cc(2)+cc(3)+cc(4)+cc(5)+cc(6))) 
   ! 
   if (npara==0 .or. natmo==0) return 
   !   derivees partielles / parametres de temperature exospherique et 
   !                         des concentrations en H, HE, O, N2, O2, N 
   !                         de T120, TP120 
   if (libatm(1,1)  = 1) then 
     dpra(1,1) = dtinf * dtt(1) 
   end if 
   if (libatm(1,2)  = 1) then 
     dpra(1,2) = vma(1) * fz(1) * dh(1) 
   end if 
   if (libatm(1,3)  = 1) then 
     dpra(1,3) = vma(2) * fz(2) * dhe(1) 
   end if 
   if (libatm(1,4)  = 1) then 
     dpra(1,4) = vma(3) * fz(3) * xdo(1) 
   end if 
   if (libatm(1,5)  = 1) then 
     dpra(1,5) = vma(4) * fz(4) * daz2(1) 
   end if 
   if (libatm(1,6)  = 1) then 
     dpra(1,6) = vma(5) * fz(5) * do2(1) 
   end if 
   if (libatm(1,7)  = 1) then 
     dpra(1,7) = vma(6) * fz(6) * daz(1) 
   end if 
   if (libatm(1,8)  = 1) then 
     dpra(1,8) = dt120 * dt0(1) 
   end if 
   if (libatm(1,9)  = 1) then 
     dpra(1,9) = dtp120 * dtp(1) 
   end if 
   do i = 2,nlatm 
     if (libatm(i,1)  = 1) then 
       dpra(i,1) = tt(1) * dtinf * dtt(i) 
     end if 
     if (libatm(i,2)  = 1) then 
       dpra(i,2) = dbc(1) * dh(i) 
     end if 
     if (libatm(i,3)  = 1) then 
       dpra(i,3) = dbc(2) * dhe(i) 
     end if 
     if (libatm(i,4)  = 1) then 
       dpra(i,4) = dbc(3) * xdo(i) 
     end if 
     if (libatm(i,5)  = 1) then 
       dpra(i,5) = dbc(4) * daz2(i) 
     end if 
     if (libatm(i,6)  = 1) then 
       dpra(i,6) = dbc(5) * do2(i) 
     end if 
     if (libatm(i,7)  = 1) then 
       dpra(i,7) = dbc(6) * daz(i) 
     end if 
     if (libatm(i,8)  = 1) then 
       dpra(i,8) = t0_dtm(1) * dt120 * dt0(i) 
     end if 
     if (libatm(i,9)  = 1) then 
       dpra(i,9) = tp(1) * dtp120 * dtp(i) 
     end if 
   end do 
 end subroutine dtm 

  
 subroutine dtm94(day,f,fbar,ap,alti,hl,alat,xlon,tz,tinf,ro,wmm,ddh) 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*auteur 
   !*version decembre 94 
   !*but calcul de la temperature, concentrations, densite totale 
   !*par entrees 
   !     day=jour de l annee 
   !     f=flux instantane a t - 1j 
   !     fbar=flux moyen a t 
   !     akp= kp tri-horaire (avec un delai de 6-abs(alat)*0.033 en heure) 
   !     alti=altitude en km superieure a 120 km 
   !     hl=heure locale + 12h (en radian) 
   !     aphi= latitude du point en radian 
   !     alon= longitude du point en radian 
   !*par sorties 
   !     tz=temperature a l altitude alti 
   !     tinf=temperature exospherique 
   !     dbc(1)=concentration hydrogene atomique 
   !     dbc(2)=concentration helium 
   !     dbc(3)=concentration oxygene atomique 
   !     dbc(4)=concentration azote moleculaire 
   !     dbc(5)=concentration oxygene moleculaire 
   !     ro=densite totale en g/cm3 
   !     wmm=masse moleculaire moyenne 
   !     ddh=derivee de la densite par rapport a l'altitude en km 
   ! 
   ! ARGUMENTS 
   ! --------- 
   !*ARG  day   E  jour de l annee 
   !*ARG  f     E  flux instantane a t - 1j et derive 
   !*ARG  fbar  E  flux moyen a t et derive 
   !*ARG  ap    E? ap moyen (DTM), a t, -1j et derivees 
   !*ARG  alti  E  altitude en km superieue a 120 km 
   !*ARG  hl    E  heure locale + 12h 
   !*ARG  alat  E  latitude du point en radian 
   !*ARG  xlon  E  longitude du point en radian 
   !*ARG  tz    S  temperature a l altitude alti 
   !*ARG  tinf  S  temperature exospherique 
   !*ARG  ro    S  densite totale en g/cm3 
   !*ARG  wmm   S  masse moleculaire moyenne 
   !*ARG  ddh   S  derivee de la densite par rapport a l'altitude en km 
   !      
   ! MODULES 
   ! ------- 
   ! 
   ! APPELES 
   ! ------- 
   ! 
   ! APPELANTS 
   ! --------- 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use common_plgdtm 
   use common_hlocal 
   use common_cons 
   use gldtm94_interface 
   use bint_interface 
   ! pour inline bint 
   use egalite_interface 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
   integer, parameter :: nlatm94 = 39 
   ! 
   !.. Formal Arguments .. 
   real(SINGLE), intent(in) :: day 
   real(SINGLE), intent(in) :: f 
   real(SINGLE), intent(in) :: fbar 
   real(SINGLE), dimension(7), intent(in) :: ap 
   real(SINGLE), intent(in) :: alti 
   real(SINGLE), intent(in) :: hl 
   real(SINGLE), intent(in) :: alat 
   real(SINGLE), intent(in) :: xlon 
   real(SINGLE), intent(out) :: tz 
   real(SINGLE), intent(out) :: tinf 
   real(SINGLE), intent(out) :: ro 
   real(SINGLE), intent(out) :: wmm 
   real(SINGLE), intent(out) :: ddh 
   ! 
   !.. Local Scalars .. 
   integer :: i,ikp 
   real(SINGLE) :: cpmg = .19081_S 
   real(SINGLE) :: gsurf = 980.665_S 
   real(SINGLE) :: re = 6356.77_S 
   real(SINGLE) :: rgas = 831.4_S 
   real(SINGLE) :: spmg = .98163_S 
   real(SINGLE) :: xlmg = -1.2392_S 
   real(SINGLE) :: zlb = 120._S 
   real(SINGLE) :: akp,apt,c,c2,c4,clmlmg,cmg,cmg2,cmg4,dzeta,expsz,fbm150, & 
                   fmfb,gamma,gdelaz2,gdelh,gdelhe,gdelo,gdelo2,gdelt,glb,sbc,s2, & 
                   sigma,sigzeta,sp,t120,t120tz,tinftz,tp120,upapg,zeta 
   ! 
   !.. Local Arrays .. 
   integer, dimension(6) :: ma 
   real(SINGLE), dimension(6) :: alefa,cc,dbc,dbase,fz,vma 
   real(SINGLE), dimension(nlatm94) :: & 
     az2,daz2,dh,dhe,do2,dtt,h,he,o,o2,t0,tp,tt,xdo 
   ! 
   !.. Intrinsic Functions .. 
   intrinsic abs, asin, cos, exp, sin 
   ! 
   !.. Data Declarations .. 
   data alefa/-0.40_S,-0.38_S,0._S,0._S,0._S,0._S/ 
   data ma/1,4,16,28,32,14/ 
   data vma/ & 
        1.6606e-24_S,6.6423e-24_S,26.569e-24_S,46.4958e-24_S,53.1381e-24_S,23.2479e-24_S/ 
   data (tt(i), i = 1,39)/ & 
        1000.5_S,0.94610e-02_S,0.42671e-01_S,0.17948e-02_S,-0.79900e-05_S,0.33674e-02_S, & 
        0.22629e-01_S,0.37865e-01_S,-0.19230e-01_S,-0.92415e-02_S,-0.21066e03_S, & 
        0.10326e-01_S,0.28856e-01_S,-0.76311e02_S,-0.18498e00_S,-0.20306e-01_S, & 
        0.14473e-01_S,-0.36306e01_S,-0.28940e-01_S,-0.17367e03_S,-0.10819e00_S, & 
        -0.19993e-02_S,0.33966e-02_S,-0.16132e-01_S,-0.96095e-02_S,-0.10455e00_S, & 
        0.45752e-02_S,0.45844e-02_S,0.17764e-01_S,-0.42270e-02_S,-0.35667e-02_S, & 
        -0.36061e-03_S,0.10493e-01_S,0.45712e-02_S,-0.21751e-04_S,0.15110e-02_S, & 
        0.21666e-02_S,-0.11423e-05_S,0.65823e-04_S/ 
   data (h(i), i = 1,39)/ & 
        1.761e+05_S,-1.33700e-01_S,0._S,-1.24600e-02_S,0._S,-1.93000e-02_S,-6.00000e-02_S, & 
        -0.20000e-02_S,5.87800e-02_S,0._S,9.22700e01_S,0._S,0._S,0._S,3.30100e-01_S, & 
        1.04500e-01_S,0._S,-1.47700e01_S,-9.06500e-02_S,-7.20000e01_S,2.09400e-01_S, & 
        2.83000e-02_S,0._S,8.57100e-02_S,-2.47500e-02_S,3.83000e-01_S,2.94100e-02_S,0._S, & 
        -3.97400e-03_S,4.35600e-02_S,0._S,0._S,0._S,0._S,0._S,0._S,0._S,0._S,0._S/ 
   data (he(i), i = 1,39)/ & 
        2.791e+07_S,0.10965e00_S,-0.19084e00_S,-0.20772e-03_S,0.48346e-05_S,0.21123e-02_S, & 
        0.22119e-03_S,-0.16174e00_S,-0.92223e-01_S,-0.83007e-02_S,0.21353e03_S, & 
        0.23504e00_S,-0.79051e-01_S,0.11040e03_S,-0.12677e01_S,-0.18512e-01_S, & 
        0.66625e-01_S,-0.18703e03_S,-0.42149e-01_S,-0.21671e03_S,-0.12779e00_S, & 
        -0.61824e-02_S,-0.17450e-01_S,-0.43600e-01_S,-0.52850e-01_S,0.31204e00_S, & 
        -0.21372e-01_S,-0.24529e-01_S,-0.36727e-02_S,-0.88330e-01_S,0.33991e-01_S, & 
        0.51314e-02_S,-0.14740e-01_S,0.85563e-02_S,0.19571e-02_S,-0.45034e-02_S, & 
        -0.10241e00_S,-0.16867e-04_S,-0.14265e-02_S/ 
   data (o(i), i = 1,39)/ & 
        0.8472e+11_S,-0.66447e-01_S,-0.97415e-01_S,0.12284e-02_S,0.44976e-05_S, & 
        0.53580e-02_S,0.25573e-02_S,-0.98221e-01_S,0.10091e00_S,0.62609e-02_S,0.11565e02_S, & 
        0.17600e00_S,-0.71284e-01_S,0.10639e03_S,0.33295e00_S,-0.11448e00_S,-0.42425e-02_S, & 
        -0.67511e00_S,-0.41712e-01_S,0.13443e03_S,-0.65932e-01_S,-0.19344e-01_S, & 
        -0.85846e-02_S,0.86347e-01_S,0.85679e-01_S,0.45486e-01_S,-0.38719e-01_S, & 
        -0.12872e-01_S,-0.91778e-01_S,0.42617e-01_S,0.43386e-01_S,0.56894e-02_S, & 
        0.66512e-02_S,-0.10174e-01_S,0.27058e-02_S,-0.47104e-02_S,-0.14639e-01_S, & 
        -0.64057e-05_S,0.15178e-02_S/ 
   data (az2(i), i = 1,39)/ & 
        3.2045e+11_S,-0.14019e00_S,0.57220e-01_S,0.11260e-02_S,-0.20779e-05_S, & 
        0.34239e-02_S,-0.11588e-01_S,+0.55158e-01_S,-0.10049e-01_S,-0.43905e-01_S, & 
        0.19594e03_S,0.32483e-01_S,0.56500e-01_S,0.88202e02_S,0.28812e00_S,-0.31427e-01_S, & 
        0.00000e00_S,-0.20015e03_S,0.60700e-01_S,0.51153e02_S,-0.46104e-01_S, & 
        -0.97000e-02_S,0.34910e-02_S,0._S,0._S,-0.73329e-01_S,0.21159e-01_S,0.70334e-02_S,0._S, & 
        0._S,-0.69199e-02_S,0._S,-0.72899e-02_S,0._S,-0.29200e-02_S,0.27047e-02_S, & 
        -0.59939e-02_S,0.00000e00_S,0.82127e-02_S/ 
   data o2(1)/4.775e+10_S/ 
   data t0(1)/380.0_S/ 
   data tp(1)/14.348_S/ 
   ! 
   ! ... Executable Statements ... 
   ! 
   ! 
   ddh = 0._S 
   ro = 0._S 
   ! 
   !   calcul des polynomes de legendre 
   c = sin(alat) 
   c2 = c * c 
   c4 = c2 * c2 
   sbc = cos(alat) 
   s2 = sbc * sbc 
   p10 = c 
   p20 = 1.5_S*c2 - 0.5_S 
   p30 = c * (2.5_S*c2-1.5_S) 
   p40 = 4.375_S*c4 - 3.75_S*c2 + 0.375_S 
   p50 = c * (7.875_S*c4-8.75_S*c2+1.875_S) 
   p11 = sbc 
   p21 = 3._S * c * sbc 
   p31 = sbc * (7.5_S*c2-1.5_S) 
   p41 = c * sbc * (17.5_S*c2-7.5_S) 
   p51 = sbc * (39.375_S*c4-26.25_S*c2+1.875_S) 
   p22 = 3._S * s2 
   p32 = 15._S * c * s2 
   p42 = s2 * (52.5_S*c2-7.5_S) 
   p52 = 3._S*c*p42 - 2._S*p32 
   p33 = 15._S * sbc * s2 
   ! 
   !   calcul des polynomes de legendre / pole magnetique (79n,71w) 
   clmlmg = cos(xlon-xlmg) 
   sp = sbc*cpmg*clmlmg + c*spmg 
   if (sp   1._S) then 
     sp = 1.0_S 
   end if 
   if (sp < -1._S) then 
     sp = -1._S 
   end if 
   ! 
   cmg = sp ! pole magnetique 
   cmg2 = cmg * cmg 
   cmg4 = cmg2 * cmg2 
   p10mg = cmg 
   p20mg = 1.5_S*cmg2 - 0.5_S 
   p40mg = 4.375_S*cmg4 - 3.75_S*cmg2 + 0.375_S 
   ! 
   !   heure locale 
   hl0 = hl + pi 
   ch = cos(hl0) 
   sh = sin(hl0) 
   c2h = ch*ch - sh*sh 
   s2h = 2._S * ch * sh 
   c3h = c2h*ch - s2h*sh 
   s3h = s2h*ch + c2h*sh 
   ! 
   !   flux 
   fmfb = f - fbar 
   fbm150 = fbar - 150._S 
   ! 
   !   kp 
   apt = ap(4) + 2._S*abs(asin(cmg))/pi*(ap(3)-ap(4)) 
   akp = bint(apt) 
   ! 
   !   calcul de la fonction g(l) / tinf, t120, tp120 
   ikp = 1 
   call gldtm94(fmfb,fbm150,akp,day,tt,gdelt,1._S,ikp) 
   dtt(1) = 1._S + gdelt 
   tinf = tt(1) * dtt(1) 
   t120 = t0(1) 
   tp120 = tp(1) 
   !  Attention BUG: pour certaines configurations tinf est inferieur 
   !  a t120 ce qui est impossible. 
   !  on impose donc tinf =t120+1 si tinf <ou= t120 
   if (tinf <= t120) then 
     tinf = t120 + 1._S 
   end if 
   ! 
   !   calcul des concentrations n(z): h, he, o, n2, o2 
   sigma = tp120 / (tinf-t120) 
   dzeta = (re+zlb) / (re+alti) 
   zeta = (alti-zlb) * dzeta 
   sigzeta = sigma * zeta 
   expsz = exp(-sigzeta) 
   tz = tinf - (tinf-t120)*expsz 
   ! 
   !   calcul de la fonction g(l) / h, he, o, n2, o2 
   ikp = 2 
   call gldtm94(fmfb,fbm150,akp,day,h,gdelh,0._S,ikp) 
   dh(1) = exp(gdelh) 
   dbase(1) = h(1) * dh(1) 
   call gldtm94(fmfb,fbm150,akp,day,he,gdelhe,0._S,ikp) 
   dhe(1) = exp(gdelhe) 
   dbase(2) = he(1) * dhe(1) 
   call gldtm94(fmfb,fbm150,akp,day,o,gdelo,1._S,ikp) 
   xdo(1) = exp(gdelo) 
   dbase(3) = o(1) * xdo(1) 
   call gldtm94(fmfb,fbm150,akp,day,az2,gdelaz2,1._S,ikp) 
   daz2(1) = exp(gdelaz2) 
   dbase(4) = az2(1) * daz2(1) 
   call gldtm94(fmfb,fbm150,akp,day,o2,gdelo2,1._S,ikp) 
   do2(1) = exp(gdelo2) 
   dbase(5) = o2(1) * do2(1) 
   ! 
   glb = gsurf / (1._S+zlb/re)**2 
   glb = glb / (sigma*rgas*tinf) 
   t120tz = t120 / tz 
   tinftz = tinf / tz 
   do i = 1,5 
     gamma = REAL(ma(i),SINGLE) * glb 
     upapg = 1._S + alefa(i) + gamma 
     fz(i) = t120tz**upapg * exp(-sigzeta*gamma) 
     !   concentrations en h, he, o, n2, o2 
     cc(i) = dbase(i) * fz(i) 
     !   densites en h, he, o, n2, o2 
     dbc(i) = cc(i) * vma(i) 
     ! 
     !   densite totale 
     ro = ro + dbc(i) 
     !   derivee partielle / altitude (km) 
     ddh = ddh + dbc(i)*sigma*dzeta*(1._S+alefa(i)-upapg*tinftz) 
   end do 
   ! 
   !   masse moleculaire moyenne 
   ! 
   wmm = ro / (vma(1)*(cc(1)+cc(2)+cc(3)+cc(4)+cc(5))) 
 end subroutine dtm94 

   
 subroutine gldtm(f,fbar,akp,day,a,da,gdel,ff0,xlon,kle_eq) 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*aut R.Biancale / modif S. Bruinsma 
   !*ver dec 2000 
   !*rol calcul de la fonction g(l) evoluee 
   !     a = tableau des coefficients pour calcul de g(l) pour la 
   !         temperature ou chaque constituant 
   !     da= tableau des derivees partielles dg(l)/da 
   !     ff0=1 pour l'oxygene , l'azote , l'helium , la temperature 
   !     ff0=0 pour l'hydrogene 
   !     gdel=resultat du calcul de g(l) 
   ! 
   ! ARGUMENTS 
   ! --------- 
   !*ARG  f      E     flux instantane a t - 1j et derive 
   !*ARG  fbar   E     flux moyen a t et derive 
   !*ARG  akp    E     moyennes des ap 
   !*ARG  day    E     jour de l annee 
   !*ARG  a      E     tableau des coefficients pour calcul de g(l) 
   !*ARG  da     E/S   tableau des derivees partielles dg(l)/da 
   !*ARG  gdel   S     resultat du calcul de g(l) 
   !*ARG  ff0    E     =1 pour O,N,He,T  =0 pour H 
   !*ARG  xlon   E     longitude 
   !*ARG  kle_eq E     cle indiquant saison (-1=hiver, 0=equinox, 1=ete hem. nord) 
   ! 
   ! MODULES 
   ! ------- 
   !*MOD  datmo    E   npara 
   !               NU  libatm 
   !*MOD  hlocal   E   c2h, c3h, ch, s2h, s3h, sh 
   !               NU  hl0 
   !*MOD  plgdtm   E   p10, p10mg, p11, p20, p20mg, p21, p22, p30, p31, p32, 
   !               E   p33, p40, p40mg, p41, p42, p50, p51 
   !               NU  p52 
   !*MOD  pourzpf  E   natmo 
   !               NU  coeff_corr, nasp, nb_lib_stoc, nbaro, ncdga, 
   !               NU  ncdgc, ncdgd, ncdgl, nfroto, ngrav, ngravd, nmar, 
   !               NU  nnbias, npolh, npolk, npolt, npolx, npoly, 
   !               NU  nporto, npreso, nrbias, nredif, nst, nst_tra, 
   !               NU  ntbias, ntherm, ntopo, nxbias, nybias, nzbias, 
   !               NU  tfrot, tnbias, tpres, trbias, ttbias, txbias, 
   !               NU  tybias, tzbias 
   ! 
   ! APPELANTS 
   ! --------- 
   ! dtm 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use common_pourzpf 
   use common_plgdtm 
   use common_hlocal 
   use common_datmo 
   use parameters_divers 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
   ! 
   !.. Formal Arguments .. 
   real(SINGLE), dimension(2), intent(in) :: f 
   real(SINGLE), dimension(2), intent(in) :: fbar 
   real(SINGLE), dimension(4), intent(in) :: akp 
   real(SINGLE), intent(in) :: day 
   real(SINGLE), dimension(nlatm), intent(in) :: a 
   real(SINGLE), dimension(nlatm), intent(inout) :: da 
   real(SINGLE), intent(out) :: gdel 
   real(SINGLE), intent(in) :: ff0 
   real(SINGLE), intent(in) :: xlon 
   ! 
   !.. Local Scalars .. 
   integer :: i,ikp,ikpm,kle_eq 
   real(SINGLE) :: rot = .017214206_S 
   real(SINGLE) :: rot2 = .034428412_S 
   real(SINGLE) :: c2fi,clfl,cos2te,coste,dakp,dakpm,dkp,dkpm,f0,f1f,fp,fp1, & 
                   rsin2te,rsinte,slfl,a74,a77,a78,a88,a89,a90,a91 
   ! 
   !.. Local Arrays .. 
   real(SINGLE), dimension(2) :: fbm150,fmfb 
   ! 
   !.. Intrinsic Functions .. 
   intrinsic cos, sin 
   ! 
   ! ... Executable Statements ... 
   ! 
   !                2*pi/24   ,      2*pi/86400 
   ! 
   !   termes de latitude 
   da(2)  = p20 
   da(3)  = p40 
   da(74) = p10 
   da(77) = p30 
   da(78) = p50 
   da(79) = p60 
   a74    = a(74)  !seule les coeff. legendre impairs changent de signe 
   a77    = a(77) 
   a78    = a(78) 
   if(kle_eq.eq. -1) then      !hiver hemisphere nord 
     a74 = -a74 
     a77 = -a77 
     a78 = -a78 
   endif 
   if(kle_eq .eq. 0) then       !equinox printemps ou automne 
     call coefintp(day,a74) 
     call coefintp(day,a77) 
     call coefintp(day,a78) 
   endif 
   !   termes de flux 
   fmfb(1) = f(1) - fbar(1) 
   fmfb(2) = f(2) - fbar(2) 
   fbm150(1) = fbar(1) - 150._S 
   fbm150(2) = fbar(2) 
   da(4) = fmfb(1) 
   da(6) = fbm150(1) 
   da(4) = da(4) + a(70)*fmfb(2) 
   da(6) = da(6) + a(71)*fbm150(2) 
   da(70) = fmfb(2) * (a(4)+2._S*a(5)*da(4)+a(82)*p10+a(83)*p20+a(84)*p30) 
   da(71) = fbm150(2) * (a(6)+2._S*a(69)*da(6)+a(85)*p10+a(86)*p20+a(87)*p30) 
   da(5) = da(4) * da(4) 
   da(69) = da(6) * da(6) 
   da(82) = da(4) * p10 
   da(83) = da(4) * p20 
   da(84) = da(4) * p30 
   da(85) = da(6) * p10 
   da(86) = da(6) * p20 
   da(87) = da(6) * p30 
   !   termes de kp 
   ikp = 62 
   ikpm = 67 
   c2fi = 1._S - p10mg*p10mg 
   dkp = akp(1) + (a(ikp)+c2fi*a(ikp+1))*akp(2) 
   dakp = a(7) + a(8)*p20mg + a(68)*p40mg + & 
          2._S*dkp*(a(60)+a(61)*p20mg+a(75)*2._S*dkp*dkp) 
   da(ikp) = dakp * akp(2) 
   da(ikp+1) = da(ikp) * c2fi 
   dkpm = akp(3) + a(ikpm)*akp(4) 
   dakpm = a(64) + a(65)*p20mg + a(72)*p40mg + & 
           2._S*dkpm*(a(66)+a(73)*p20mg+a(76)*2._S*dkpm*dkpm) 
   da(ikpm) = dakpm * akp(4) 
   da(7) = dkp 
   da(8) = p20mg * dkp 
   da(68) = p40mg * dkp 
   da(60) = dkp * dkp 
   da(61) = p20mg * da(60) 
   da(75) = da(60) * da(60) 
   da(64) = dkpm 
   da(65) = p20mg * dkpm 
   da(72) = p40mg * dkpm 
   da(66) = dkpm * dkpm 
   da(73) = p20mg * da(66) 
   da(76) = da(66) * da(66) 
   !   fonction g(l) non periodique 
   f0 = a(4)*da(4) + a(5)*da(5) + a(6)*da(6) + a(69)*da(69) + a(82)*da(82) + & 
        a(83)*da(83) + a(84)*da(84) + a(85)*da(85) + a(86)*da(86) + & 
        a(87)*da(87) 
   f1f = 1._S + f0*ff0 
   f0 = f0 + a(2)*da(2) + a(3)*da(3) + a74*da(74) + a77*da(77) + a(7)*da(7) & 
        + a(8)*da(8) + a(60)*da(60) + a(61)*da(61) + a(68)*da(68) + & 
        a(64)*da(64) + a(65)*da(65) + a(66)*da(66) + a(72)*da(72) + & 
        a(73)*da(73) + a(75)*da(75) + a(76)*da(76) + a78*da(78) + a(79)*da(79) 
   !   termes annuels symetriques en latitude 
   da(9) = cos(rot*(day-a(11))) 
   da(10) = p20 * da(9) 
   !   termes semi-annuels symetriques en latitude 
   da(12) = cos(rot2*(day-a(14))) 
   da(13) = p20 * da(12) 
   !   termes annuels non symetriques en latitude 
   coste = cos(rot*(day-a(18))) 
   da(15) = p10 * coste 
   da(16) = p30 * coste 
   da(17) = p50 * coste 
   !   terme  semi-annuel  non symetrique  en latitude 
   cos2te = cos(rot2*(day-a(20))) 
   da(19) = p10 * cos2te 
   da(39) = p30 * cos2te 
   da(59) = p50 * cos2te 
   !   termes diurnes (et couples annuel) 
   da(21) = p11 * ch 
   da(22) = p31 * ch 
   da(23) = p51 * ch 
   da(24) = da(21) * coste 
   da(25) = p21 * ch * coste 
   da(26) = p11 * sh 
   da(27) = p31 * sh 
   da(28) = p51 * sh 
   da(29) = da(26) * coste 
   da(30) = p21 * sh * coste 
   !   termes semi-diurnes (et couples annuel) 
   da(31) = p22 * c2h 
   da(37) = p42 * c2h 
   da(32) = p32 * c2h * coste 
   da(33) = p22 * s2h 
   da(38) = p42 * s2h 
   da(34) = p32 * s2h * coste 
   da(88) = p32 * c2h 
   da(89) = p32 * s2h 
   da(90) = p52 * c2h 
   da(91) = p52 * s2h 
   a88    = a(88) 
   a89    = a(89) 
   a90    = a(90) 
   a91    = a(91) 
   if(kle_eq .eq. -1) then 
     a88 = -a88 
     a89 = -a89 
     a90 = -a90 
     a91 = -a91 
   endif 
   if(kle_eq .eq. 0) then 
     call coefintp(day,a88) 
     call coefintp(day,a89) 
     call coefintp(day,a90) 
     call coefintp(day,a91) 
   endif 
   da(92) = p62 * c2h 
   da(93) = p62 * s2h 
   !   termes ter-diurnes 
   da(35) = p33 * c3h 
   da(36) = p33 * s3h 
   !   fonction g(l) periodique 
   fp = a(9)*da(9) + a(10)*da(10) + a(12)*da(12) + a(13)*da(13) + a(15)*da(15) & 
        + a(16)*da(16) + a(17)*da(17) + a(19)*da(19) + a(21)*da(21) + & 
        a(22)*da(22) + a(23)*da(23) + a(24)*da(24) + a(25)*da(25) + & 
        a(26)*da(26) + a(27)*da(27) + a(28)*da(28) + a(29)*da(29) + & 
        a(30)*da(30) + a(31)*da(31) + a(32)*da(32) + a(33)*da(33) + & 
        a(34)*da(34) + a(35)*da(35) + a(36)*da(36) + a(37)*da(37) + & 
        a(38)*da(38) + a(39)*da(39) + a(59)*da(59) + & 
        a88 * da(88) + a89 * da(89) + a90 * da(90) + a91 * da(91) + & 
        a(92)*da(92) + a(93)*da(93) 
   ! 
   !   termes d'activite magnetique 
   da(40) = p10 * coste * dkp 
   da(41) = p30 * coste * dkp 
   da(42) = p50 * coste * dkp 
   da(43) = p11 * ch * dkp 
   da(44) = p31 * ch * dkp 
   da(45) = p51 * ch * dkp 
   da(46) = p11 * sh * dkp 
   da(47) = p31 * sh * dkp 
   da(48) = p51 * sh * dkp 
   ! 
   !   fonction g(l) periodique supplementaire 
   fp = fp + a(40)*da(40) + a(41)*da(41) + a(42)*da(42) + a(43)*da(43) + & 
        a(44)*da(44) + a(45)*da(45) + a(46)*da(46) + a(47)*da(47) + & 
        a(48)*da(48) 
   dakp = (a(40)*p10+a(41)*p30+a(42)*p50)*coste + & 
          (a(43)*p11+a(44)*p31+a(45)*p51)*ch + & 
          (a(46)*p11+a(47)*p31+a(48)*p51)*sh 
   da(ikp) = da(ikp) + dakp*akp(2) 
   da(ikp+1) = da(ikp) + dakp*c2fi*akp(2) 
   !   termes de longitude 
   clfl = cos(xlon) 
   da(49) = p11 * clfl 
   da(50) = p21 * clfl 
   da(51) = p31 * clfl 
   da(52) = p41 * clfl 
   da(53) = p51 * clfl 
   slfl = sin(xlon) 
   da(54) = p11 * slfl 
   da(55) = p21 * slfl 
   da(56) = p31 * slfl 
   da(57) = p41 * slfl 
   da(58) = p51 * slfl 
   ! 
   !   fonction g(l) periodique supplementaire 
   fp = fp + a(49)*da(49) + a(50)*da(50) + a(51)*da(51) + a(52)*da(52) + & 
        a(53)*da(53) + a(54)*da(54) + a(55)*da(55) + a(56)*da(56) + & 
        a(57)*da(57) + a(58)*da(58) 
   ! 
   !   fonction g(l) totale (couplage avec le flux) 
   gdel = f0 + fp*f1f 
   ! 
   !   derivees partielles 
   if (natmo/=0 .and. npara/=0) then 
     !   derivees des phases des termes annuels et semi-annuels 
     da(11) = (a(9)+a(10)*p20) * rot * sin(rot*(day-a(11))) 
     da(14) = (a(12)+a(13)*p20) * rot2 * sin(rot2*(day-a(14))) 
     rsinte = rot * sin(rot*(day-a(18))) 
     da(18) = & 
       rsinte * & 
       (a(15)*p10+a(16)*p30+a(17)*p50+(a(24)*p11+a(25)*p21)*ch+ & 
        (a(29)*p11+a(30)*p21)*sh+a(32)*p32*c2h+a(34)*p32*s2h) 
     rsin2te = rot2 * sin(rot2*(day-a(20))) 
     da(20) = rsin2te * (a(19)*p10+a(39)*p30+a(59)*p50) 
     !   derivees des phases dues aux termes supplementaires 
     da(18) = da(18) + rsinte*dkp*(a(40)*p10+a(41)*p30+a(42)*p50) 
   end if 
   !   couplages en flux 
   fp1 = 1._S + fp*ff0 
   da(4) = da(4) * fp1 
   da(5) = da(5) * fp1 
   da(6) = da(6) * fp1 
   da(69) = da(69) * fp1 
   da(70) = da(70) * fp1 
   da(71) = da(71) * fp1 
   da(82) = da(82) * fp1 
   da(83) = da(83) * fp1 
   da(84) = da(84) * fp1 
   da(85) = da(85) * fp1 
   da(86) = da(86) * fp1 
   da(87) = da(87) * fp1 
   do i = 9,59 
     da(i) = da(i) * f1f 
   end do 
 end subroutine gldtm 

  
 subroutine gldtm92(fmfb,fbm150,akp,day,a,da,gdel,ff0) 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*aut r.biancale 
   !*ver sept 1991 
   !*rol calcul de la fonction g(l) completee voir hedin(73/86) 
   !     a = tableau des coefficients pour calcul de g(l) pour la 
   !         temperature ou chaque constituant 
   !     da= tableau des derivees partielles dg(l)/da 
   !     ff0=1 pour oxygene atomique, azote et temperature 
   !     ff0=0 pour l helium 
   !     gdel=resultat du calcul de g(l) 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use common_pourzpf 
   use common_plgdtm 
   use common_hlocal 
   use common_datmo 
   use parameters_divers 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
   ! 
   !.. Formal Arguments .. 
   real(SINGLE), intent(in) :: fmfb 
   real(SINGLE), intent(in) :: fbm150 
   real(SINGLE), intent(in) :: akp 
   real(SINGLE), intent(in) :: day 
   real(SINGLE), dimension(nlatm), intent(in) :: a 
   real(SINGLE), dimension(nlatm), intent(inout) :: da 
   real(SINGLE), intent(out) :: gdel 
   real(SINGLE), intent(in) :: ff0 
   ! 
   !.. Local Scalars .. 
   integer :: i 
   real(SINGLE) :: rot = .017214206_S 
   real(SINGLE) :: rot2 = .034428412_S 
   real(SINGLE) :: coste,f0,f1f,fp,fp1,rsinte 
   ! 
   !.. Intrinsic Functions .. 
   intrinsic cos, sin 
   ! 
   ! ... Executable Statements ... 
   ! 
   !                2*pi/24   ,      2*pi/86400 
   ! 
   !   termes de latitude 
   da(2) = p20 
   da(3) = p40 
   !   termes de flux 
   da(4) = fmfb 
   da(5) = da(4) * da(4) 
   da(6) = fbm150 
   !   termes de kp 
   da(7) = akp 
   da(8) = p20 * akp 
   da(39) = p40 * akp 
   !   fonction g(l) non periodique 
   f0 = a(4)*da(4) + a(5)*da(5) + a(6)*da(6) 
   f1f = 1._S + f0*ff0 
   f0 = f0 + a(2)*da(2) + a(3)*da(3) + a(7)*da(7) + a(8)*da(8) + a(39)*da(39) 
   !   termes annuels symetriques en latitude 
   da(9) = cos(rot*(day-a(11))) 
   da(10) = p20 * da(9) 
   !   termes semi-annuels symetriques en latitude 
   da(12) = cos(rot2*(day-a(14))) 
   da(13) = p20 * da(12) 
   !   termes annuels non symetriques en latitude 
   coste = cos(rot*(day-a(18))) 
   da(15) = p10 * coste 
   da(16) = p30 * coste 
   da(17) = p50 * coste 
   !   et couples en activite magnetique 
   da(40) = p10 * coste * akp 
   da(41) = p30 * coste * akp 
   da(42) = p50 * coste * akp 
   !   terme  semi-annuel  non symetrique  en latitude 
   da(19) = p10 * cos(rot2*(day-a(20))) 
   !   termes diurnes (et couples annuel) 
   da(21) = p11 * ch 
   da(22) = p31 * ch 
   da(23) = p51 * ch 
   da(24) = da(21) * coste 
   da(25) = p21 * ch * coste 
   da(26) = p11 * sh 
   da(27) = p31 * sh 
   da(28) = p51 * sh 
   da(29) = da(26) * coste 
   da(30) = p21 * sh * coste 
   !   termes semi-diurnes (et couples annuel) 
   da(31) = p22 * c2h 
   da(37) = p42 * c2h 
   da(32) = p32 * c2h * coste 
   da(33) = p22 * s2h 
   da(38) = p42 * s2h 
   da(34) = p32 * s2h * coste 
   !   termes ter-diurnes 
   da(35) = p33 * c3h 
   da(36) = p33 * s3h 
   !   fonction g(l) periodique 
   fp = a(9)*da(9) + a(10)*da(10) + a(12)*da(12) + a(13)*da(13) + a(15)*da(15) & 
        + a(16)*da(16) + a(17)*da(17) + a(19)*da(19) + a(21)*da(21) + & 
        a(22)*da(22) + a(23)*da(23) + a(24)*da(24) + a(25)*da(25) + & 
        a(26)*da(26) + a(27)*da(27) + a(28)*da(28) + a(29)*da(29) + & 
        a(30)*da(30) + a(31)*da(31) + a(32)*da(32) + a(33)*da(33) + & 
        a(34)*da(34) + a(35)*da(35) + a(36)*da(36) + a(37)*da(37) + & 
        a(38)*da(38) + a(40)*da(40) + a(41)*da(41) + a(42)*da(42) 
   ! 
   !   fonction g(l) totale (couplage avec le flux) 
   gdel = f0 + fp*f1f 
   ! 
   !   derivees partielles 
   if (natmo/=0 .and. npara/=0) then 
     !   derivees des phases des termes annuels et semi-annuels 
     da(11) = (a(9)+a(10)*p20) * rot * sin(rot*(day-a(11))) 
     da(14) = (a(12)+a(13)*p20) * rot2 * sin(rot2*(day-a(14))) 
     rsinte = rot * sin(rot*(day-a(18))) 
     da(18) = & 
       rsinte* & 
       (a(15)*p10+a(16)*p30+a(17)*p50+(a(24)*p11+a(25)*p21)*ch+ & 
        (a(29)*p11+a(30)*p21)*sh+a(32)*p32*c2h+a(34)*p32*s2h) + & 
       rsinte*akp*(a(40)*p10+a(41)*p30+a(42)*p50) 
     da(20) = & 
       a(19)*p10*rot2*sin(rot2*(day-a(20))) + & 
       (a(40)*p10+a(41)*p30+a(42)*p50)*coste 
     !   couplages en ap 
     !      fac=a(7)+a(8)*p20+a(39)*p40 
     !   couplages en flux 
     fp1 = 1._S + fp*ff0 
     da(4) = da(4) * fp1 
     da(5) = da(5) * fp1 
     da(6) = da(6) * fp1 
     do i = 9,42 
       da(i) = da(i) * f1f 
     end do 
   end if 
 end subroutine gldtm92 

  
 subroutine gldtm94(fmfb,fbm150,akp,day,a,gdel,ff0,ikp) 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*auteur 
   !*version decembre 94 
   !*role calcul de la fonction g(l) evoluee 
   !     a = tableau des coefficients pour calcul de g(l) pour la 
   !         temperature ou chaque constituant 
   !     ff0=1 pour l'oxygene , l'azote , la temperature 
   !     ff0=0 pour l'hydrogene , l'helium' 
   !     ikp=1 pour la temperature 
   !     ikp=2 pour chaque constituant 
   !     gdel=resultat du calcul de g(l) 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use common_plgdtm 
   use common_hlocal 
   use parameters_divers 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
   ! 
   !.. Formal Arguments .. 
   real(SINGLE), intent(in) :: fmfb 
   real(SINGLE), intent(in) :: fbm150 
   real(SINGLE), intent(in) :: akp 
   real(SINGLE), intent(in) :: day 
   real(SINGLE), dimension(nlatm94), intent(in) :: a 
   real(SINGLE), intent(out) :: gdel 
   real(SINGLE), intent(in) :: ff0 
   integer, intent(in) :: ikp 
   ! 
   !.. Local Scalars .. 
   real(SINGLE) :: rot = .017214206_S 
   real(SINGLE) :: rot2 = .034428412_S 
   real(SINGLE) :: cos11,cos14,cos18,cos20,f0,f1f,fp,tansl,tasl,tdca,tflu,tkp, & 
                   tphi,tsansl,tsasl,tsdca,tterd 
   ! 
   !.. Intrinsic Functions .. 
   intrinsic cos, exp 
   ! 
   ! ... Executable Statements ... 
   ! 
   !                2*pi/24   ,      2*pi/86400 
   ! 
   !   termes de latitude 
   tphi = a(2)*p20 + a(3)*p40 + a(37)*p10 
   !   termes de flux 
   tflu = a(4)*fmfb + a(5)*fmfb*fmfb + a(6)*fbm150 + a(38)*fbm150*fbm150 
   !   termes de kp 
   if (ikp == 1) then 
     tkp = a(7)*akp + a(8)*p20mg*akp + a(39)*exp(akp) 
   end if 
   if (ikp == 2) then 
     tkp = a(7)*akp + a(8)*p20mg*akp + a(39)*akp*akp 
   end if 
   !   fonction g(l) non periodique 
   f1f = 1._S + tflu*ff0 
   f0 = tphi + tflu + tkp 
   !   termes annuels symetriques en latitude 
   cos11 = cos(rot*(day-a(11))) 
   tasl = (a(9)+a(10)*p20) * cos11 
   !   termes semi-annuels symetriques en latitude 
   cos14 = cos(rot2*(day-a(14))) 
   tsasl = (a(12)+a(13)*p20) * cos14 
   !   termes annuels non symetriques en latitude 
   cos18 = cos(rot*(day-a(18))) 
   tansl = (a(15)*p10+a(16)*p30+a(17)*p50) * cos18 
   !   terme  semi-annuel  non symetrique  en latitude 
   cos20 = cos(rot2*(day-a(20))) 
   tsansl = (a(19)*p10) * cos20 
   !   termes diurnes (et couples annuel) 
   tdca = a(21)*p11*ch + a(22)*p31*ch + a(23)*p51*ch + a(24)*p11*ch*cos18 + & 
          a(25)*p21*ch*cos18 + a(26)*p11*sh + a(27)*p31*sh + a(28)*p51*sh + & 
          a(29)*p11*sh*cos18 + a(30)*p21*sh*cos18 
   !   termes semi-diurnes (et couples annuel) 
   tsdca = a(31)*p22*c2h + a(32)*p32*c2h*cos18 + a(33)*p22*s2h + & 
           a(34)*p32*s2h*cos18 
   !   termes ter-diurnes 
   tterd = a(35)*p33*c3h + a(36)*p33*s3h 
   !   fonction g(l) periodique 
   fp = tasl + tsasl + tansl + tsansl + tdca + tsdca + tterd 
   ! 
   !   fonction g(l) totale (couplage avec le flux) 
   ! 
   gdel = f0 + fp*f1f 
 end subroutine gldtm94 

  

 subroutine lecdtm 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*aut r.biancale 
   !*ver fev 2001 - modification pour mars (annaig 08/99; s. bruinsma) 
   !*rol lecture des donnees relatives au modele d'atmospherique dtm 
   !     (si ikle(4)=2(dtm) ou 7(dtmars)) 
   ! 
   ! MODULES 
   ! ------- 
   !*MOD  central  E?  corps_central 
   !*MOD  fic      E   nfic, nlpr 
   !               NU  iu10, iu11, iu12, iu13, iu15, iu3, iu5, iu6, iu8, 
   !               NU  iu9, iuhorlgps 
   !*MOD  pardtm   S?  ar, az, az2, co, co2, h, h2, he, o, o2, t0, t0_dtm, tp, 
   !               S?  tt 
   ! 
   ! APPELANTS 
   ! --------- 
   ! entree 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use parameters_divers 
   use common_pardtm 
   use common_fic 
   use common_central 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
   ! 
   !.. Local Scalars .. 
   character(LEN=8) :: ct120,ctp120 
   character(LEN=80) :: titre 
   integer :: i,j,ncol,ncol1,ni,npdtm 
   ! 
   !.. Local Arrays .. 
   real(SINGLE), dimension(12) :: p 
   ! 
   ! ... Executable Statements ... 
   ! 
   ! 
   write (nlpr,10004) 
   ! 
   read (nfic,10000) titre 
   if (corps_central(1:5)=='terre') then 
   write (nlpr,10003) titre 
   ct120 = titre(57:64) 
   ctp120 = titre(73:80) 
   read (nfic,10001) npdtm, ncol 
   ncol1 = ncol + 1 
   do i = 1,npdtm 
     read (nfic,10002,end = 1000) ni, (p(j), j = 1,ncol) 
     do j = ncol1,9 
       p(j) = 0._S 
     end do 
     tt(i) = p(1) 
     h(i) = p(2) 
     he(i) = p(3) 
     o(i) = p(4) 
     az2(i) = p(5) 
     o2(i) = p(6) 
     az(i) = p(7) 
     t0_dtm(i) = p(8) 
     tp(i) = p(9) 
   end do 
   if (ncol < 8) then 
     read (ct120,'(f8.3)') t0_dtm(1) 
   end if 
   if (ncol < 9) then 
     read (ctp120,'(f8.4)') tp(1) 
   end if 
   do i = npdtm+1,nlatm 
     tt(i) = 0._S 
     h(i) = 0._S 
     he(i) = 0._S 
     o(i) = 0._S 
     az2(i) = 0._S 
     az(i) = 0._S 
     t0_dtm(i) = 0._S 
     tp(i) = 0._S 
   end do 
   else if (corps_central(1:4)=='mars') then 
     write (nlpr,10003) titre 
     read (nfic,10001) npdtm, ncol 
     ncol1 = ncol + 1 
     do i = 1,npdtm 
       read (nfic,10005,end = 1000) ni, (p(j), j = 1,ncol) 
       do j = ncol1,12 
         p(j) = 0._S 
       end do 
       tt(i) = p(1) 
       t0(i) = p(2) 
       tp(i) = p(3) 
       co2(i)= p(4) 
       o(i)  = p(5) 
       az2(i)= p(6) 
       ar(i) = p(7) 
       co(i) = p(8) 
       o2(i) = p(9) 
       he(i) = p(10) 
       h(i)  = p(11) 
       h2(i) = p(12) 
     end do 
     do i = npdtm+1,nlatm 
       tt(i) = 0._S 
       t0_dtm(i) = 0._S 
       tp(i) = 0._S 
       co2(i) = 0._S 
       o(i) = 0._S 
       az2(i) = 0._S 
       ar(i) = 0._S 
       co(i) = 0._S 
       o2(i) = 0._S 
       he(i) = 0._S 
       h(i) = 0._S 
       h2(i) = 0._S 
     end do 
 ! tempete de poussiere: parametrer DANS fichier dtm-mars 
     if(npdtm  = 74) then 
       tau   = tt(npdtm-3)  !opacite (IR) max. de l'atmosphere martienne pendant tempete 
       xlsdeb= tt(npdtm-2)  !Ls du debut de la tempete  (deg) 
       xlsfin= tt(npdtm-1)  !Ls de la fin de la tempete (deg) 
       taubg = tt(npdtm)    !opacite de l'atmosphere martienne moyenne (background dust opacity) 
       if(taubg.lt.0.15_S) taubg=0._S !average opacity during MGS Mapping, and thus in model 
       write(nlpr,*) ' ATTENTION, tempete de poussiere: tau=',tau,' debut=',xlsdeb 
     else 
       tau   = 0._S 
       xlsdeb= 0._S 
       xlsfin= 0._S 
       taubg = 0._S 
     endif 
   end if 
   return 
  
   1000 write (nlpr,*) ' e.o.f. sur tape', nfic, ' dans lecdtm' 
   stop 
   ! 
   ! ... Format Declarations ... 
   ! 
   !*fon les formats 
   !     ----------- 
   10000 format (a80) 
   10001 format (2i4) 
   10002 format (i4,9(e13.6,9x)) 
   10003 format (5x,a80) 
   10004 format (/,' ****  modele d atmosphere dtm  ****') 
   10005 format (i4,12(e13.6,9x)) 
 !610   format(5x,'natmo :',i5) 
 end subroutine lecdtm 

  
 subroutine libdtm(incpos) 
   ! 
   !*********************************************************************** 
   ! 
   ! HISTORIQUE 
   ! ---------- 
   ! 
   ! ROLE 
   ! ---- 
   !*aut j.c marty 
   !*ver dec 1992 
   !*rol lecture des donnees relatives a la liberation du modele 
   !     d'atmosphere dtm (si ikle(4)=2) 
   ! 
   ! ARGUMENTS 
   ! --------- 
   !*ARG  incpos  E? indice de position dans le tableau des inconnues 
   ! 
   ! MODULES 
   ! ------- 
   !*MOD  central  E?     corps_central 
   !*MOD  datmo    S      npara 
   !               E?/S?  libatm 
   !*MOD  dopdop   S?     llib, zpf 
   !               NU     icode, lks1, poi, qat, zdj, zfemi, zfrec, zqm 
   !*MOD  fic      E      nfic 
   !               E?     nlpr 
   !               NU     iu10, iu11, iu12, iu13, iu15, iu3, iu5, iu6, iu8, 
   !               NU     iu9, iuhorlgps 
   !*MOD  itersup  S?     nds 
   !               NU     inuti1, ipoli1, libs, nis, npar1, nparb1, 
   !               NU     npard1, nparm1, npart1, qlibs 
   !*MOD  pardtm   E?     ar, az, az2, co, co2, h, h2, he, o, o2, t0, t0_dtm, 
   !               E?     tp, tt 
   !*MOD  poursup  E?/S?  natmo1 
   !               NU     ks, nasp1, nbaro1, ncdga_bas1, ncdgd_bas1, 
   !               NU     ncdgg_bas1, ncdgl_bas1, ncdgpd_bas1, 
   !               NU     ncdgpr_bas1, ngrav1, ngravd1, nlovh1, 
   !               NU     nlovk1, nmar1, nnute1, nnutp1, npolh1, 
   !               NU     npolk1, npolt1, npolx1, npoly1, nqua1, nst1, 
   !               NU     nst1_tra, ntopo1 
   !*MOD  pourzpf  E?/S?  natmo 
   !               NU     coeff_corr, nasp, nb_lib_stoc, nbaro, ncdga, 
   !               NU     ncdgc, ncdgd, ncdgl, nfroto, ngrav, ngravd, 
   !               NU     nlovh, nlovk, nmar, nnbias, nnute, nnutp, 
   !               NU     npolh, npolk, npolt, npolx, npoly, nporto, 
   !               NU     npreso, nqua, nrbias, nredif, nst, nst_tra, 
   !               NU     ntbias, ntherm, ntopo, nxbias, nybias, 
   !               NU     nzbias, tfrot, tnbias, tpres, trbias, ttbias, 
   !               NU     txbias, tybias, tzbias 
   ! 
   ! APPELANTS 
   ! --------- 
   ! entree 
   ! 
   !*********************************************************************** 
   ! 
   use f90_kind 
   use parameters_divers 
   use common_pourzpf 
   use common_poursup 
   use common_pardtm 
   use common_itersup 
   use common_fic 
   use common_dopdop 
   use common_datmo 
   use common_central 
   ! 
   !.. Implicit Declarations .. 
   implicit none 
   ! 
   !.. Parameters .. 
     ! 
   !.. Formal Arguments .. 
   integer, intent(in) :: incpos 
   ! 
   !.. Local Scalars .. 
   integer :: i,ind,j 
   ! 
   ! ... Executable Statements ... 
   ! 
   read (nfic,10000,end = 1000) npara 
   if (npara == 0) then 
     write (nlpr,10004) 
   else if (npara /= 0) then 
     if (corps_central(1:5)=='terre') then 
       write (nlpr,10002) 
       ind = incpos 
       read (nfic,10001,end = 1000) ((libatm(i,j), i = 1,nlatm), j = 1,npara) 
       do j = 1,npara 
         do i = 1,nlatm 
           if (libatm(i,j) == 1) then 
             natmo = natmo + 1 
             ind = ind + 1 
             if (j == 1) then 
               zpf(ind) = tt(i) 
             end if 
             if (j == 2) then 
               zpf(ind) = h(i) 
             end if 
             if (j == 3) then 
               zpf(ind) = he(i) 
             end if 
             if (j == 4) then 
               zpf(ind) = o(i) 
             end if 
             if (j == 5) then 
               zpf(ind) = az2(i) 
             end if 
             if (j == 6) then 
               zpf(ind) = o2(i) 
             end if 
             if (j == 7) then 
               zpf(ind) = az(i) 
             end if 
             if (j == 8) then 
               zpf(ind) = t0_dtm(i) 
             end if 
             if (j == 9) then 
               zpf(ind) = tp(i) 
             end if 
             llib(ind) = 1 
           elseif (libatm(i,j) == 2) then ! pour iter. sup. 
             natmo1 = natmo1 + 1 
             ind = ind + 1 
             if (j == 1) then 
               zpf(ind) = tt(i) 
             end if 
             if (j == 2) then 
               zpf(ind) = h(i) 
             end if 
             if (j == 3) then 
               zpf(ind) = he(i) 
             end if 
             if (j == 4) then 
               zpf(ind) = o(i) 
             end if 
             if (j == 5) then 
               zpf(ind) = az2(i) 
             end if 
             if (j == 6) then 
               zpf(ind) = o2(i) 
             end if 
             if (j == 7) then 
               zpf(ind) = az(i) 
             end if 
             if (j == 8) then 
               zpf(ind) = t0_dtm(i) 
             end if 
             if (j == 9) then 
               zpf(ind) = tp(i) 
             end if 
             llib(ind) = 0 
             nds(ind) = 1 
           end if 
         end do 
       end do 
       write (nlpr,10003) natmo 
       if (natmo1 /= 0) then 
         natmo = natmo1 
       end if 
     else if (corps_central(1:4)=='mars') then !on travaille sur mars 
       write (nlpr,10002) 
       ind = incpos 
       read (nfic,10001,end = 1001) ((libatm(i,j), i = 1,nlatm), j = 1,npara) 
       do j = 1,npara 
         do i = 1,nlatm 
           if (libatm(i,j) == 1) then 
             natmo = natmo + 1 
             ind = ind + 1 
             if (j == 1) then 
               zpf(ind) = tt(i) 
             end if 
             if (j == 2) then 
               zpf(ind) = t0(i) 
             end if 
             if (j == 3) then 
               zpf(ind) = tp(i) 
             end if 
             if (j == 4) then 
               zpf(ind) = co2(i) 
             end if 
             if (j == 5) then 
               zpf(ind) = o(i) 
             end if 
             if (j == 6) then 
               zpf(ind) = az2(i) 
             end if 
             if (j == 7) then 
               zpf(ind) = ar(i) 
             end if 
             if (j == 8) then 
               zpf(ind) = co(i) 
             end if 
             if (j == 9) then 
               zpf(ind) = o2(i) 
             end if 
             if (j == 10) then 
               zpf(ind) = he(i) 
             end if 
             if (j == 11) then 
               zpf(ind) = h(i) 
             end if 
             if (j == 12) then 
               zpf(ind) = h2(i) 
             end if 
             llib(ind) = 1 
             elseif (libatm(i,j) == 2) then ! pour iter. sup. 
             natmo1 = natmo1 + 1 
             ind = ind + 1 
             if (j == 1) then 
               zpf(ind) = tt(i) 
             end if 
             if (j == 2) then 
               zpf(ind) = t0(i) 
             end if 
             if (j == 3) then 
               zpf(ind) = tp(i) 
             end if 
             if (j == 4) then 
               zpf(ind) = co2(i) 
             end if 
             if (j == 5) then 
               zpf(ind) = o(i) 
             end if 
             if (j == 6) then 
               zpf(ind) = az2(i) 
             end if 
             if (j == 7) then 
               zpf(ind) = ar(i) 
             end if 
             if (j == 8) then 
               zpf(ind) = co(i) 
             end if 
             if (j == 9) then 
               zpf(ind) = o2(i) 
             end if 
             if (j == 10) then 
               zpf(ind) = he(i) 
             end if 
             if (j == 11) then 
               zpf(ind) = h(i) 
             end if 
             if (j == 12) then 
               zpf(ind) = h2(i) 
             end if 
             llib(ind) = 0 
             nds(ind) = 1 
           end if 
         end do 
       end do 
       write (nlpr,10003) natmo 
       if (natmo1 /= 0) then 
         natmo = natmo1 
       end if 
     end if 
  end if 
 1000 return 
 1001 return 
     ! 
     ! ... Format Declarations ... 
     ! 
 10000 format (i5) 
 10001 format (87i1) 
 10002 format (/,' ****  liberation des coefficients du modele dtm ****') 
 10003 format (5x,'****  natmo :',i5) 
 10004 format (/,' ****  coefficients du modele dtm non liberes ****') 
 end subroutine libdtm