docs: def gtd7(Input, flags, output):
      '''The standard model subroutine (GTD7) always computes the 
      ‘‘thermospheric’’ mass density by explicitly summing the masses of
      the species in equilibrium at the thermospheric temperature T(z).
      '''
      mn3 = 5
      zn3 = [32.5,20.0,15.0,10.0,0.0]
      mn2 = 4
      zn2 = [72.5,55.0,45.0,32.5]
      zmix = 62.5
      soutput = nrlmsise_output()
  
      tselec(flags);
  
      #/* Latitude variation of gravity (none for sw[2]=0) */
      xlat=Input.g_lat;
      if (flags.sw[2]==0): # pragma: no cover
          xlat=45.0;
      glatf(xlat, gsurf, re);
  
      xmm = pdm[2][4];
  
      #/* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
      if (Input.alt>zn2[0]):
          altt=Input.alt;
      else:
          altt=zn2[0];
  
      tmp=Input.alt;
      Input.alt=altt;
  
      gts7(Input, flags, soutput);
      altt=Input.alt;
      Input.alt=tmp;
      if (flags.sw[0]): # pragma: no cover  #/* metric adjustment */
          dm28m= dm28*1.0E6;
      else:
          dm28m = dm28;
      output.t[0]=soutput.t[0];
      output.t[1]=soutput.t[1];
      if (Input.alt>=zn2[0]): 
          for i in range(9):
              output.d[i]=soutput.d[i];
          return
  
  
  #/*       LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
  #*         Temperature at nodes and gradients at end nodes
  #*         Inverse temperature a linear function of spherical harmonics
  #*/
      meso_tgn2[0]=meso_tgn1[1];
      meso_tn2[0]=meso_tn1[4];
      meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags.sw[20]*glob7s(pma[0], Input, flags));
      meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags.sw[20]*glob7s(pma[1], Input, flags));
      meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags.sw[20]*flags.sw[22]*glob7s(pma[2], Input, flags));
      meso_tgn2[1]=pavgm[8]*pma[9][0]*(1.0+flags.sw[20]*flags.sw[22]*glob7s(pma[9], Input, flags))*meso_tn2[3]*meso_tn2[3]/(pow((pma[2][0]*pavgm[2]),2.0));
      meso_tn3[0]=meso_tn2[3];
  
      if (Input.alt<zn3[0]):
  #/*       LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
  #*         Temperature at nodes and gradients at end nodes
  #*         Inverse temperature a linear function of spherical harmonics
  #*/
          meso_tgn3[0]=meso_tgn2[1];
          meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags.sw[22]*glob7s(pma[3], Input, flags));
          meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags.sw[22]*glob7s(pma[4], Input, flags));
          meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags.sw[22]*glob7s(pma[5], Input, flags));
          meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags.sw[22]*glob7s(pma[6], Input, flags));
          meso_tgn3[1]=pma[7][0]*pavgm[7]*(1.0+flags.sw[22]*glob7s(pma[7], Input, flags)) *meso_tn3[4]*meso_tn3[4]/(pow((pma[6][0]*pavgm[6]),2.0));
      
  
      #/* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */
  
      dmc=0;
      if (Input.alt>zmix):
          dmc = 1.0 - (zn2[0]-Input.alt)/(zn2[0] - zmix);
      dz28=soutput.d[2];
      
      #/**** N2 density ****/
      dmr=soutput.d[2] / dm28m - 1.0;
      tz = [0.0]
      output.d[2]=densm(Input.alt,dm28m,xmm, tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
      output.d[2]=output.d[2] * (1.0 + dmr*dmc);
  
      #/**** HE density ****/
      dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
      output.d[0] = output.d[2] * pdm[0][1] * (1.0 + dmr*dmc);
  
      #/**** O density ****/
      output.d[1] = 0;
      output.d[8] = 0;
  
      #/**** O2 density ****/
      dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
      output.d[3] = output.d[2] * pdm[3][1] * (1.0 + dmr*dmc);
  
      #/**** AR density ***/
      dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
      output.d[4] = output.d[2] * pdm[4][1] * (1.0 + dmr*dmc);
  
      #/**** Hydrogen density ****/
      output.d[6] = 0;
  
      #/**** Atomic nitrogen density ****/
      output.d[7] = 0;
  
      #/**** Total mass density */
      output.d[5] = 1.66E-24 * (4.0 * output.d[0] + 16.0 * output.d[1] + 28.0 * output.d[2] + 32.0 * output.d[3] + 40.0 * output.d[4] + output.d[6] + 14.0 * output.d[7]);
  
      if (flags.sw[0]): # pragma: no cover
          output.d[5]=output.d[5]/1000;
  
      #/**** temperature at altitude ****/
      global dd
      dd = densm(Input.alt, 1.0, 0, tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
      output.t[1]=tz[0];
      return