docs: def gts7(Input, flags, output):
      '''
  /*     Thermospheric portion of NRLMSISE-00
   *     See GTD7 for more extensive comments
   *     alt > 72.5 km! 
   */
   '''
      zn1 = [120.0, 110.0, 100.0, 90.0, 72.5]
      mn1 = 5
      dgtr=1.74533E-2;
      dr=1.72142E-2;
      alpha = [-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0]
      altl = [200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0]
      za = pdl[1][15];
      zn1[0] = za;
  
      for j in range(9):
          output.d[j]=0;
  
      #/* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
      if (Input.alt>zn1[0]):
          tinf = ptm[0]*pt[0] * \
                      (1.0+flags.sw[16]*globe7(pt,Input,flags));
      else:
          tinf = ptm[0]*pt[0];
      output.t[0]=tinf;
  
      #/*  GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
      if (Input.alt>zn1[4]):
          g0 = ptm[3]*ps[0] * \
              (1.0+flags.sw[19]*globe7(ps,Input,flags));
      else:
          g0 = ptm[3]*ps[0];
      tlb = ptm[1] * (1.0 + flags.sw[17]*globe7(pd[3],Input,flags))*pd[3][0];
      s = g0 / (tinf - tlb);
  
  #/*      Lower thermosphere temp variations not significant for
  # *       density above 300 km */
      if (Input.alt<300.0):
          meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags.sw[18]*glob7s(ptl[0], Input, flags));
          meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags.sw[18]*glob7s(ptl[1], Input, flags));
          meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags.sw[18]*glob7s(ptl[2], Input, flags));
          meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags.sw[18]*flags.sw[20]*glob7s(ptl[3], Input, flags));
          meso_tgn1[1]=ptm[8]*pma[8][0]*(1.0+flags.sw[18]*flags.sw[20]*glob7s(pma[8], Input, flags))*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
      else:
          meso_tn1[1]=ptm[6]*ptl[0][0];
          meso_tn1[2]=ptm[2]*ptl[1][0];
          meso_tn1[3]=ptm[7]*ptl[2][0];
          meso_tn1[4]=ptm[4]*ptl[3][0];
          meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
  	
  
      z0 = zn1[3];
      t0 = meso_tn1[3];
      tr12 = 1.0;
  
      #/* N2 variation factor at Zlb */
      g28=flags.sw[21]*globe7(pd[2], Input, flags);
  
      #/* VARIATION OF TURBOPAUSE HEIGHT */
      zhf=pdl[1][24]*(1.0+flags.sw[5]*pdl[0][24]*sin(dgtr*Input.g_lat)*cos(dr*(Input.doy-pt[13])));
      output.t[0]=tinf;
      xmm = pdm[2][4];
      z = Input.alt;
  
  
      #/**** N2 DENSITY ****/
  
      #/* Diffusive density at Zlb */
      db28 = pdm[2][0]*exp(g28)*pd[2][0];
      #/* Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[2];
      #/* Turbopause */
      zh28=pdm[2][2]*zhf;
      zhm28=pdm[2][3]*pdl[1][5]; 
      xmd=28.0-xmm;
      #/* Mixed density at Zlb */
      tz = [0]
      b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
      if ((flags.sw[15]) and (z<=altl[2])):
          #/*  Mixed density at Alt */
          global dm28
          dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          #/*  Net density at Alt */
          output.d[2]=dnet(output.d[2],dm28,zhm28,xmm,28.0);
      
  
  
      #/**** HE DENSITY ****/
  
      #/*   Density variation factor at Zlb */
      g4 = flags.sw[21]*globe7(pd[0], Input, flags);
      #/*  Diffusive density at Zlb */
      db04 = pdm[0][0]*exp(g4)*pd[0][0];
      #/*  Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[0];
      if ((flags.sw[15]) and (z<altl[0])):
          #/*  Turbopause */
          zh04=pdm[0][2];
          #/*  Mixed density at Zlb */
          RandomVariable = [output.t[1]]
          b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          #/*  Mixed density at Alt */
          RandomVariable = [output.t[1]]
          global dm04
          dm04=densu(z,b04,tinf,tlb,xmm,0.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          zhm04=zhm28;
          #/*  Net density at Alt */
          output.d[0]=dnet(output.d[0],dm04,zhm04,xmm,4.);
          #/*  Correction to specified mixing ratio at ground */
          rl=log(b28*pdm[0][1]/b04);
          zc04=pdm[0][4]*pdl[1][0];
          hc04=pdm[0][5]*pdl[1][1];
          #/*  Net density corrected at Alt */
          output.d[0]=output.d[0]*ccor(z,rl,hc04,zc04);
      
  
  
      #/**** O DENSITY ****/
  
      #/*  Density variation factor at Zlb */
      g16= flags.sw[21]*globe7(pd[1],Input,flags);
      #/*  Diffusive density at Zlb */
      db16 =  pdm[1][0]*exp(g16)*pd[1][0];
      #/*   Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],RandomVariable,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[1];
      if ((flags.sw[15]) and (z<=altl[1])):
          #/*   Turbopause */
          zh16=pdm[1][2];
          #/*  Mixed density at Zlb */
          RandomVariable = [output.t[1]]
          b16=densu(zh16,db16,tinf,tlb,16.0-xmm,(alpha[1]-1.0), RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          #/*  Mixed density at Alt */
          RandomVariable = [output.t[1]]
          global dm16
          dm16=densu(z,b16,tinf,tlb,xmm,0.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          zhm16=zhm28;
          #/*  Net density at Alt */
          output.d[1]=dnet(output.d[1],dm16,zhm16,xmm,16.);
          rl=pdm[1][1]*pdl[1][16]*(1.0+flags.sw[1]*pdl[0][23]*(Input.f107A-150.0));
          hc16=pdm[1][5]*pdl[1][3];
          zc16=pdm[1][4]*pdl[1][2];
          hc216=pdm[1][5]*pdl[1][4];
          output.d[1]=output.d[1]*ccor2(z,rl,hc16,zc16,hc216);
          #/*   Chemistry correction */
          hcc16=pdm[1][7]*pdl[1][13];
          zcc16=pdm[1][6]*pdl[1][12];
          rc16=pdm[1][3]*pdl[1][14];
          #/*  Net density corrected at Alt */
          output.d[1]=output.d[1]*ccor(z,rc16,hcc16,zcc16);
      
  
  
      #/**** O2 DENSITY ****/
  
      #/*   Density variation factor at Zlb */
      g32= flags.sw[21]*globe7(pd[4], Input, flags);
      #/*  Diffusive density at Zlb */
      db32 = pdm[3][0]*exp(g32)*pd[4][0];
      #/*   Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],RandomVariable,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[3];
      if (flags.sw[15]):
          if (z<=altl[3]):
              #/*   Turbopause */
              zh32=pdm[3][2];
              #/*  Mixed density at Zlb */
              RandomVariable = [output.t[1]]
              b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
              output.t[1] = RandomVariable[0]
              #/*  Mixed density at Alt */
              RandomVariable = [output.t[1]]
              global dm32
              dm32=densu(z,b32,tinf,tlb,xmm,0.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
              output.t[1] = RandomVariable[0]
              zhm32=zhm28;
              #/*  Net density at Alt */
              output.d[3]=dnet(output.d[3],dm32,zhm32,xmm,32.);
              #/*   Correction to specified mixing ratio at ground */
              rl=log(b28*pdm[3][1]/b32);
              hc32=pdm[3][5]*pdl[1][7];
              zc32=pdm[3][4]*pdl[1][6];
              output.d[3]=output.d[3]*ccor(z,rl,hc32,zc32);
          
          #/*  Correction for general departure from diffusive equilibrium above Zlb */
          hcc32=pdm[3][7]*pdl[1][22];
          hcc232=pdm[3][7]*pdl[0][22];
          zcc32=pdm[3][6]*pdl[1][21];
          rc32=pdm[3][3]*pdl[1][23]*(1.+flags.sw[1]*pdl[0][23]*(Input.f107A-150.));
          #/*  Net density corrected at Alt */
          output.d[3]=output.d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
      
  
  
      #/**** AR DENSITY ****/
  
      #/*   Density variation factor at Zlb */
      g40= flags.sw[21]*globe7(pd[5],Input,flags);
      #/*  Diffusive density at Zlb */
      db40 = pdm[4][0]*exp(g40)*pd[5][0];
      #/*   Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[4];
      if ((flags.sw[15]) and (z<=altl[4])):
          #/*   Turbopause */
          zh40=pdm[4][2];
          #/*  Mixed density at Zlb */
          RandomVariable = [output.t[1]]
          b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          #/*  Mixed density at Alt */
          RandomVariable = [output.t[1]]
          global dm40
          dm40=densu(z,b40,tinf,tlb,xmm,0.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          zhm40=zhm28;
          #/*  Net density at Alt */
          output.d[4]=dnet(output.d[4],dm40,zhm40,xmm,40.);
          #/*   Correction to specified mixing ratio at ground */
          rl=log(b28*pdm[4][1]/b40);
          hc40=pdm[4][5]*pdl[1][9];
          zc40=pdm[4][4]*pdl[1][8];
          #/*  Net density corrected at Alt */
          output.d[4]=output.d[4]*ccor(z,rl,hc40,zc40);
        
  
  
      #/**** HYDROGEN DENSITY ****/
  
      #/*   Density variation factor at Zlb */
      g1 = flags.sw[21]*globe7(pd[6], Input, flags);
      #/*  Diffusive density at Zlb */
      db01 = pdm[5][0]*exp(g1)*pd[6][0];
      #/*   Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[6];
      if ((flags.sw[15]) and (z<=altl[6])):
          #/*   Turbopause */
          zh01=pdm[5][2];
          #/*  Mixed density at Zlb */
          RandomVariable = [output.t[1]]
          b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          #/*  Mixed density at Alt */
          RandomVariable = [output.t[1]]
          global dm01
          dm01=densu(z,b01,tinf,tlb,xmm,0.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          zhm01=zhm28;
          #/*  Net density at Alt */
          output.d[6]=dnet(output.d[6],dm01,zhm01,xmm,1.);
          #/*   Correction to specified mixing ratio at ground */
          rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01);
          hc01=pdm[5][5]*pdl[1][11];
          zc01=pdm[5][4]*pdl[1][10];
          output.d[6]=output.d[6]*ccor(z,rl,hc01,zc01);
          #/*   Chemistry correction */
          hcc01=pdm[5][7]*pdl[1][19];
          zcc01=pdm[5][6]*pdl[1][18];
          rc01=pdm[5][3]*pdl[1][20];
          #/*  Net density corrected at Alt */
          output.d[6]=output.d[6]*ccor(z,rc01,hcc01,zcc01);
  
  
  
      #/**** ATOMIC NITROGEN DENSITY ****/
  
      #/*   Density variation factor at Zlb */
      g14 = flags.sw[21]*globe7(pd[7],Input,flags);
      #/*  Diffusive density at Zlb */
      db14 = pdm[6][0]*exp(g14)*pd[7][0];
      #/*   Diffusive density at Alt */
      RandomVariable = [output.t[1]]
      output.d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      dd=output.d[7];
      if ((flags.sw[15]) and (z<=altl[7])): 
          #/*   Turbopause */
          zh14=pdm[6][2];
          #/*  Mixed density at Zlb */
          RandomVariable = [output.t[1]]
          b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          #/*  Mixed density at Alt */
          RandomVariable = [output.t[1]]
          global dm14
          dm14=densu(z,b14,tinf,tlb,xmm,0.,RandomVariable,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
          output.t[1] = RandomVariable[0]
          zhm14=zhm28;
          #/*  Net density at Alt */
          output.d[7]=dnet(output.d[7],dm14,zhm14,xmm,14.);
          #/*   Correction to specified mixing ratio at ground */
          rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14);
          hc14=pdm[6][5]*pdl[0][1];
          zc14=pdm[6][4]*pdl[0][0];
          output.d[7]=output.d[7]*ccor(z,rl,hc14,zc14);
          #/*   Chemistry correction */
          hcc14=pdm[6][7]*pdl[0][4];
          zcc14=pdm[6][6]*pdl[0][3];
          rc14=pdm[6][3]*pdl[0][5];
          #/*  Net density corrected at Alt */
          output.d[7]=output.d[7]*ccor(z,rc14,hcc14,zcc14);
      
  
  
      #/**** Anomalous OXYGEN DENSITY ****/
  
      g16h = flags.sw[21]*globe7(pd[8],Input,flags);
      db16h = pdm[7][0]*exp(g16h)*pd[8][0];
      tho = pdm[7][9]*pdl[0][6];
      RandomVariable = [output.t[1]]
      dd=densu(z,db16h,tho,tho,16.,alpha[8],RandomVariable,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
      output.t[1] = RandomVariable[0]
      zsht=pdm[7][5];
      zmho=pdm[7][4];
      zsho=scalh(zmho,16.0,tho);
      output.d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));
  
  
      #/* 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]);
      db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);
  
  
  
      #/* temperature */
      z = sqrt(Input.alt*Input.alt);
      RandomVariable = [output.t[1]]
      ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, RandomVariable, ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
      output.t[1] = RandomVariable[0]
      if (flags.sw[0]): # pragma: no cover
          for i in range(9):
              output.d[i]=output.d[i]*1.0E6;
          output.d[5]=output.d[5]/1000;
      return