       PROGRAM CoupledSUA

CZ UPDATES BY MEZ/JFK SEPT 2011 to PRESENT

CZ hydro64: modified TA,TB,TC to incorp. fix of derivation re DR(x); need to fix for B.C.s; use RHO_M to calc YT
CZ hydro63: cs130: turned off DRX(J), replaced w/ DR(J); JFK also found sign error on flux diagnostic calc, fixing this in Output.f gets mass conservation to within a few percent.
CZ      changed index on TA,TB,TC back to MZ-1 from MZ-2
CZ hydro62dr: cs129dr: JFK found missing factor of 1/DR in his diffusion derivation; introduced this using DRX(J) which increases by factor of 1.003 per timestep; also made restart.edit5;
CZ     this combination allowed code to converge in <10k steps. 
CZ hydro62: Rev 129: updates to Output.f conservation debug
CZ Rev 128: Skipped short-lived sp NORM_MMR(0) (fixed primary problem); still have discontinuity in e-, NO2+,O2+ at same level; updates to Output.f for conservation debug
CZ Rev 127: Fixed RXTOT(NZ)=0; added debug; set incoming and outgoing UNe(NZ)=0 & UNUM(I,NZ)=0; changing EDDYA; 
CZ (Rev 126 skipped - used for debugging)
CZ Rev 125: completed updates of heating w/ moving UNUM and other variables to boundary
CZ hydro61: Rev 124: change tf_heating and subroutines to calc q(J+1), add COLUNUMB,TCOLUNUMB,put x__ (UNUM on bdrys) in fmodules 
CZ Rev 123: Mod Difco to calc PRES(J) if not available from subsonic; limit delT to 50K, added many warnings looking for T issue; 
CZ             works OK for 300 runs but q are not done right, eventually temps go out of bounds
CZ Rev 122: hydro60: Mods to NORM_MMR(0) to separate DEN, RHO, WT from midpoint variables USOL, UNUM, etc.; changed Output.f FLUXKZ calc to use J, J-1; 
CZ              in heating subroutines tf_heating and tf_cond, added qx(NZ) = qx(NZ-1); fixed NZ in adjust_exobase
CZ Rev 121: hydro59: Continuing to separate bdry variables from midpoint, esp. DEN, RHO in NORM_MMR;
CZ             removed old FLUX, replaced w/ FLUXKZ in fmodules.f; this version runs 10K steps, but DT = 10^-10
CZ Rev 120: hydro58: Added (midpoint) R_M, DEN_M, RHO_M to fmodules.f (atmosphere); fix calcs of WT, DEN; replaced WTM with WTP;
CZ          Fixed Ti NaN problem by setting clei, clin at NZ equal to value at NZ-1 in tf_heating;
CZ Rev 119: Now runs with exohold=0 (turned off); Fixed and updated Output.f; turned off debug statements; mod Difco and Output so only print out after last time step
CZ Rev 118: changed JS=NZ-1 (was NZ); added prints to EXO_DBG and IZDEBUGEXO; NORM_MMR no longer changes UNUM when USOL(I,NZ-1)<=SMALL;
CZ Rev 117: hydro57: changed Output.f to print fluxes using new calcs and on grid boundaries
CZ Rev 116: first version w/ new JK diffusion that runs; needs exohold=1; fixed typo in Difco, changed cp(NZ) to cp(NZ-1)
CZ         Need to fix:  (1)fluxes should be similar to reaction rates in interior grid points (diagnosic?); (2)BCs at top and bottom; (3)deal w/ changing exobase
CZ Rev 115: Added Error notification to Rev 114 chg; mod calcs of GZ, r2 now go to 600; removed bad recalc of zk; 
CZ         set drp=DR(J), drm=DR(J-1); removed tfm lines;
CZ Rev 114: Set TA, TB, TC NaN to 0
CZ Rev 113: Removed unused subroutine Densty (left as stub to make);
CZ Rev 112: Mod Rates and Photo to use NZ-1
CZr Rev 111: Rerwrote Difco, modified main and Output to fix flux calcs and grid
CZ Rev 110: Stripped out: global mean code, variations on grids, conversion of Earth to Mars, other commented out code (total over 1200 lines);
CZ     hydro56: mod Difco, subsonic, and fmodules to calculate pressure, calc WT and DEN at the midpoints; define DR in fmodules;
CZ Rev 105: IZLOGGRID = -1 now uses JK grid as per his subsonic code where R(J) are grid boundaries; mod Difco to fix calc of HD, HDA,
CZ Rev 104: changed IZLOGGRID =1: now uses JK grid but sets boundaries at log 0.5 above grid centers; moved rootpi to denom. in main code calc of ujeans 
CZ Rev 103: modified calc of grid at bottom end (intermediate rev, do not use)
CZ Rev 102: Implemented JK chemistry mods Difco and main; re-enabled U_ADV_save in second instance; added IZHUSOLDBG to turn off USOL<0 warnings;
CZ    set R0 value only when IZLOGGRID=2 (otherwise use read-in value);
CZ Rev 101: moved some minor sp. MMR into planet definitions (complete later)
CZ Rev 100: renumbered, no change
CZ Rev r82a15: added DEBUGUINF;
CZ Rev r82a14: added option for u_inf=-2; modified BC output; step output incl. u_inf and ujeans; 
CZ      uncommented U_ADV(J) = U_ADV_save(J)
CZ Rev r82a13: fixed i_exo initialization; exobase i_exo allowed to vary w/o limit if <11; IZPRNTFLG8 controls display of u_inf calc from sum of VEFF(I); 
CZ ***end of incorporating changes from Rev 82b through 86
CZ Rev r82a12: cleaned up MMR control section; added boundary condition printout at end; added IZPRNTFLG7
CZ Rev r82a11: changed calc of TB(1)=TB(2); set TA(1)=TC(1)=0 for LBC=1; added IZBCDEBUG;
CZ Rev r82a10: LBOUND renamed LBOUND1 in main and subroutines; added HO2 and OH to calc of flux in subsonic_endoutput;
CZ   added more output to tridiagonal solver failed message.
CZ Rev r82a09: Boundary condition flags added -  IZUBC, IZLBC;
CZ Rev r82a08: changed EDD(J) back to 1.0E6; fixed line numbers; added IZFLUXDEBUG and table of chem in/out; 
CZ ***above incorporate changes after Rev r82a
CZ Rev r82a07: set iexohold >0 to hold exobase at level determined at time step N=iexohold (no calls to find_exobase
CZ     or adjust_exobase); set OMMR=0.01 and other minor constituents to TF levels; mod Output.f to print fluxes of 3 species
CZ Rev r82a06: additional control on i_exo limits to +/-1 for all cases
CZ Rev r82a05: re-enabled calc of u_inf from weighted sum of VEFF(I); fixed subsonic so U_ADV(NZ)=u_inf for all cases;
CZ     fixed find_exobase and adjust_exobase for jeans mode; step output now replaces u_inf with ujeans for u_inf_ini<0
CZ Rev r82a03: fixed ujeans mode when u_inf_ini<0 
CZ Rev r82a02: completed changes from a01; set N2 = -1 to calculate by subtraction
CZ Rev r82a01: changes in chemistry for modern Earth (He,O,H,etc)
CZ Rev r82a: fixed problem with r82; always opens Inputs10.in
CZ Rev r82: PLANETNO determines which Inputs10xxx.in file is called - Do Not Use
CZ Rev r81: set up PLANETNO to fix values of minor species, select correct lower bdry conditions in fmodules.f, set psihi;
CZ      added conditions on psihi to work for psi<0.01
CZ Rev r80a: fix printout of lower boundary conditions
CZ Rev r79: hydro50: JK fixed FUP, FDOWN in Output.f; O2MMR and planet now read in from Inputs10.in; COMMR set by planet;
CZ       commented out subsonic mach error and warning; print out lower boundary conditions
CZ Rev r78: hydro49: N2MMR now read in from Inputs09.in
CZ Rev r77: set Np1=NZ+3; set exobase=ipsimaxgrid+3 (was+2)
CZ Rev r76: control exobase to always be at 2 grid points above psimax; tightened lower limit on psimax, esp. after N=100; set Np1=NZ+2
CZ Rev r75: calculate crossover mass at homopause (was exobase); reduced EDD to 7.5E5 to move homopause from 2 to 1
CZ Rev r74: limit exobase to no higher than 2 grid points above critical point (ipsimaxgrid)
CZ Rev r73: added utden column showning MMR(Htot)
CZ Rev r72c: fixed flux = sum(MMR) calcs
CZ Rev r72b: replaced MMR(O),O+ with H flux 
CZ Rev r72a: fixed utden.out, added MMR(O),MMR(O+)
CZ Rev r72: finished change in call to Output.f listed in r69a
CZ Rev r71: calculates H mass flux at every altitude for utden
CZ Rev r70: added three-point averaging to USOL(LH) and LH2 (SMOOTH3=1)
CZ Rev r69a: send N and TIME when calling Output.f;
CZ Rev r69: print out grid point J with chemistry;
CZ Rev r68: turn psimax lower limit back on; drop Tmin to 30K from 100K; moved psi>0.98 error after loop
CZ     mod utden to print protons/cm2-s; removed dt=0.6dt when exobase moved, replaced with max timestep 1.E2
CZ     parameters.f modified to increase upper limit on printout
CZ Rev r67: continued cleaning up grid setup, fixed bug in higher grid points; limit exobase to max J=401; 
CZ Rev r66: cleaning up grid setup; removed Grid.f (left stub for make)
CZ Rev r65a: removed feedback to increase psimax if below psihi - 0.2; fixed bug in Output.f (did not print DR(J))
CZ Rev r65: time step decreases if exobase changed position; changed Output.f to output H+ flux (was He flux);
CZ     cleaned up outputs including MSKIP replaces remaining ISKIP
CZ Rev r64: separate sub-loop to control psi using u_inf; calc wind equation, check psimax, and repeat
CZ    until psimax is under control (up to a max number of times ipsimaxcount), before running rho and chemistry
CZ    moved u_inf control of psimax from end of the time-stepping loop to right after call_subsonic(
CZ    removed some commented out code
CZ Rev r63: added ipsimax grid to output; reduced mulitpliers on u_inf to control psimax
CZ Rev r62: commented psihi back in; increased NEWTMAX to 1E5; increased NZchemMXN=30;
CZ     added CROSDBG, EXODEBUG
CZ Rev r61: removed u_inf recalc after call find_exobase (leaving in after adjust_exobase);
CZ     fixed bugs related to u_inf compensation
CZ      added print statements
CZ Rev r60: added u_inf to printout;
CZ Rev r59a: added recalculation after call find_exobase (was only after call adjust_exobase)
CZ Rev r59: recalculate u_inf if exobase moves
CZ Rev r58: u_inf_ini set = initial u_inf, use to decide if U_ADV(NZ)=u_inf;
CZ      comment out control on u_infinity;
CZ Rev r57c: commented out initial check of psi
CZ Rev 57b: allow reverse control on u_infinity
CZ Rev r57: insert control on u_infinity to hold psi steady
CZ Rev r56: calc psimax, report instead of U; comment out psi warnings
CZ Rev r55a: removed restrictions on species for TSOL warning (now warns for any sp < -0.001)
CZ Rev r55: still more options for CONVEM
CZ Rev r54: more options for CONVEM and EARTHMARS
CZ Rev r53: for u_inf = -1, uses ujeans
CZ Rev r52:  added EARTHMARS
CZ Rev r51: bigM now set in this program, passed to subsonic.f and adjust_exobase (was a parameter)
CZ      added bigM to CONVEM
CZ Rev r50: parameters.f change bigM for Mars; recalc altitudes for each grid point for each iteration of R_PLT; 
CZ Rev 49: hydr43: got deposition vel. working VDEP=VDEP0; set O2, CO to constant 10e-3 MMR;
CZ       turned off neg. conc. warning for sp. 2, 8;
CZ       changed 5 back to O2 decrease by 2%
CZ Rev 48: turned off NZ = 57 options; added N2MMR; changed O,CO,H to 1cm/s dep vel in fmodules.f;
CZ       Changed 5 to O2 decrease by 10%
CZ Rev 47: option to choose Earth or Mars, or convert Earth to Mars
CZ Rev 46: read in H2MMR, CO2MMR; fixed grid point 1 printout;
CZ Rev 45: do not use
CZ Rev 44a: now calc mass flux at NZ-3
CZ Rev 43: added control of CO2 MMR; H2MMR now a variable
CZ Rev 42: removed J=5 dif_limited flux
CZ Rev 41: mods dif_limited_flux - do not use
CZ Rev 40b: change target H conc under IZH2X1P01 = 9
CZ Rev 40: addl. controls on U_ADV
CZ Rev 39: added IZH2X1P01 = 10 to decrease H2 to specified level at homopause; increase temp max to 1000K
CZ     for error 6.2; changed i_exo_save limit from 300 to 600; chemkluge=0; 
CZ Rev 38: if chemkluge = 1: when chemistry tridiag solver fails, add a small number to each value, XL min SMALL 
CZ Rev 37: under chemkluge = 1, limits USOL(LH,J) and USOL(LH2,J) to SMALL or greater; remove recalc of U_ADV;
CZ     under tempkluge = 1, moved limit on delT to before error 5.2
CZ Rev 36: hydr40: added more limits on chemistry changes, moved Rev 35 changes to chemkluge = 1;
CZ    moved Rev 31 delT change to tempkluge = 1; relaxed delUNUM limits to +/- 50%; Tmin now J = 8 - NZ;
CZ    NEWTMAX = 100 (was 1E3)
CZ Rev 35: tightened limits on delUNUM
CZ Rev 34: fixed UNUMHmin and H2min;
CZ Rev 33: added dif_limit_flux, UNUM(LH2,NZ), UNUM(LH,NZ) to idlpass.misc
CZ Rev 32a: error 9 limited to species 1-14 (was 1-17)
CZ Rev 32: limits changes on USOL, UNUM; limits TSOL > -0.001; EMAX and EMAXK errors now variables
CZ      EMAXMX and EMAXKMX; utden now includes dif_limited_flux; maximum time step DTMX
CZ Rev 31: hydr39: replaced 3-point smoothing with limit on magnitude of delT
CZ Rev 30: hydr38: added 3-point smoothing of T(i) - causes much lower temperatures
CZ Rev 29: Added J to printout after error 4
CZ Rev 28: for EMAX, replaced max grid point with MMR .GE. 1E-4 and UNUM .GE. 1E10; 
CZ        max DT = 1E11 (was 1E15)
CZ Rev 27: EMAX counts only if species MMR at least 1E-4
CZ Rev 26: fixed error handling typos
CZ Rev 25a: fixed idlpass_status readin
CZ Rev 25: hydr37: T(1)=190, A_EDDY=7E26, B_EDDY=0.3
CZ Rev r24: added error for T>500, J<30; additional error EMAXK after T(i) calc 
CZ Rev r23a: utden now includes MMR and N for H and H2
CZ Rev r22a: Error for EMAX, EMAXK > 10 (was 100)
CZ Rev r22:  now T(i) < 100; when IERRDEST=1, errors call Output.f and subsonic_end_output.f (utden.out)
CZ Rev r21,21a: intermediate revs with improved error handling
CZ Rev r20a: added IF/THEN to check for T(i)<50, changed EMAX/EMAXK limits to 100
CZ Rev r19: fixed calc of comp H
CZ Rev r18: hydr35: direct control of U_ADV=u_inf; removed u_fac and other calc of u_inf, Inputs07.in;
CZ     add parameters to idlpass_misc; re-set U_ADV(NZ)=u_inf after save/return;
CZ Rev r17: hydr34: Production ver 2 - add idlpass_misc.out to create table
CZ Rev r16: hydr33: Production ver1
CZ Rev r15: Exits when psi GE 0.98, EMAX or EMAXK > 1E15 (and N >= 2), NaN,
CZ Rev r14: use Inputs06.in; mod order of idlpass_sp.out; added info on vol fraction of H
CZ Rev r13: hydr32: grid spacing and wind speed now included in Inputs05.f; wind speed now based on ujeans, not u_inf
CZ Rev r12: grid spacing now a variable ZGRIDFAC
CZ Rev r11: added u_inf_fac factor multiplying wind at exobase
CZ Rev r10: experimenting with control of DEN(1); IZMMRDEBUG now controls printout of H2 conc. 
CZ Rev r09: cleaned up r08
CZ Rev r08: mod Densty.f to fix DEN(1) = 2.5E13; was fixing RHO(1) = 1.107e-9 g/cm3
CZ Rev r07: attempt to control density at 100 km w/ temp T(1)
CZ Rev r06: mod density correction, move back into time-stepping loop
CZ Rev r05: removed density correction from time-stepping loop to after read in from restart.in
CZ Rev r04: hydr31: if statements in time-stepping loop set grid point 1 (100 km) to constant number density
CZ Number density of modern Earth atmosphere at 100 km = 1.194E19 particles/m3 = 1.194E13 part/cm3
CZ    from U.S. std atmosphere1976 http://tpsx.arc.nasa.gov/cgi-perl/alt.pl
CZ Rev r03: EDD set = 1E6
CZ Rev r02: profiles of EDD and DI(LH2) added as columns in utden.out; 
CZ       corrected Number Fractions to Mass Mixing Fractions in output
CZ Rev r01: hydr30: print out EDD(1) and DI(LH2,1)
CZ Rev q04: add IZLOGGRID = 3 option for JK grid that starts at 95 km (later removed)
CZ Rev q03: increase T(2) through T(10) as needed so T(1) is never > T(J)
CZ Rev q02: increase T(2) by half the increase in T(1)
CZ Rev q01: Inputs04.in, includes T(1) set
CZ Rev p7: mods to printout, print flags
CZ Rev p6: changed eddy diff calc: kappam=kappam(1+A_EDDY*exp(-Zcm/(B_EDDY*HSCALE(1))))
CZ Rev p5: fixed dif_limited_flux calc, prints out ftot @ 100 km
CZ Rev p4: intermediate rev
CZ Rev p3: hydr28: change calc of dif_limited_flux to normalize to homopause, not R(NZ)
CZ Rev p2: use Inputs03.in to load eddy diff factor A_EDDY, modify inputs #1-3 to make space
CZ Rev p1: hydr27: remove restriction on U-ADV(NZ) again;
CZ      first cut change heat transfer by adding coefficient A
CZ Rev n11c: calc diff. limited flux at both J=1 and J=5;
CZ Rev n11b: intermediate rev
CZ Rev n11a: put restriction back in U-ADV(NZ)<0.9, mod utden.out to include WTM; 
CZ    calc. diff limited flux at J=1 (was J=5); fix typos in n11
CZ Rev n11: hydr26:
CZ    In subsonic.f, removed restriction U_ADV(NZ)<0.8, calculate b/H at J=1 (was J=5);
CZ        calculate mass flux, fluxM
CZ    In cs, mod utden.out to include fluxM
CZ Rev n10a: hydr25: moved f107 to Inputs02.in
CZ Rev n09: fixed f107 EUV solar min/mean/max (ITF_SFLX_FLG) after emails w/ TF 2March2012
CZ Rev n08: updates hydr23 and hydr24
CZ hydr24: changed O2 back to fixed, O2 MMR=1E-3, He num density = 3E8
CZ hydr23: mod fmodules.f, made O,O2,H,H+,O+ constant dep (was const mix ratio); mod params, set damping=0.5 (was 0.9)
CZ Rev n07c: emakx>4e-1 increases DT
CZ Rev n07b: emax>0.3 and .0.2 conditions increase DT
CZ Rev n07a: reduced speed of increase in DT based on emaxk; added emax>0.3 condition, later also emax>0.1
CZ Rev n07: modified code that changes atmospheric constituents
CZ Rev n06 more debug features
CZ Rev n05: improved debugging
CZ Rev n04: mod timestep DT, cleaned up output
CZ Rev n03: modified H2 controls
CZ Rev n02: turned off code setting EMAXK=100  modified time step control
CZ Rev n01: moved NSTEPS step counter to Inputs.in
CZ Rev m8: additional Inputs.in switches to incr. CO2, decr. O2, O3, O
CZ Rev m7: made Newton limit 1E4 a variable
CZ Rev m6: test for OH, 03, HO2 all within 1E4 of each other before doing Newton
CZ Rev m5: do CHEMPLONE as input to Newton; above newton limit NH, do CHEMPLONE(OH), then O3, HO2; 
CZ Rev m4: renamed to cs_, separated chemMXN for Newton species from chemMXLL for long lived species, chemMXE for EMAX
CZ Rev m3: improved output
CZ Rev m2: added Newton species printout, MSKIP replaces ISKIP, extra print flags
CZ Rev m1: limit Newton method to max altitude of NZchemMX
CZ Rev m: fixed grid bug where program reverts to linear grid up high; mod Grid.f
CZ Rev k5: makes max exobase height a variable, i_exo_max (was 600); fixed nexodelt bug
CZ Rev k4: Set max grid point for temp calcs NZtempMX, in parameters.f
CZ Rev k3: Set max grid point to perform chemistry calcs, NZchemMX, in parameters.f
CZ Rev k2: makes damping (amount of old solution to include in new solution) a variable (was 0.9)
CZ Rev k1: makes exobase movement +/- NEXODELT
CZ Rev j: limits how much exobase can move i_exo +/-4, also Rev g mod disabling lowering of ceiling now only if IZ300.NE.0
CZ Rev h: forces ceiling up to 287 grid points; prevent this from forcing up exobase
CZ Rev g: disables lowering atm ceiling NZ; limits DT>1E-10; increases acceptable emaxk
CZ Rev f: attempt to force code to go to 300 grid points by setting IZ300 flag (not yet working)
CZ Rev e: fixed utden.out to correctly show Z in km using zk array
CZ
CZ Rev a,b,c,d: 
CZ 1) Replaced Feng's linear grid w/ JK log grid,  
CZ 2) Added capabilities to multiply H2 by x3, x10/3, or by 1/2.76... before time stepping, or 
CZ     by x3 or x10 in 1% increments at each time step
CZ 3) Added Inputs.in input file which turns features as in 2) on and off, sets values
CZ 4) Added output files idlpass.out and idlpass_sp.out to pass values to IDL for display and printing
CZ 4) Commented out tetn > 2700 warning in tf_heating .f
CZ 5) Moved close (710) statement to end, removed contains statements to eliminate errors
CZ 5) Added IZDEBUG flag and multiple print statements for debugging
CZ 6) Added IZPRNTxx flags which can be used to turn off selected print statements
CZ 7) Modified print 100 and Output.f to print out EMAXK and other useful variables
CZ 8) Added write(80, statements to provide headers in tf_profiles output
CZ 9) Modified makefile to turn off warnings
CZ


C-PK The hydrogen budget calculations have been revised.
C 100km grid, molecular diffusion,
C Escape of H2 and H
C
C         THIS PROGRAM IS DESIGNED FOR EXTREMELY LOW (PRE-PHOTOSYNTHETIC
C     O2 LEVELS.  THE PHOTOLYSIS OF O2 IN THE SCHUMANN-RUNGE BANDS CAN
C     BE CALCULATED EITHER BY THE METHOD OF ALLEN AND FREDERICK OR BY
C     EXPONENTIAL SUMS.  THE LATTER METHOD SHOULD BE USED FOR LOW-O2
C     ATMOSPHERES (I.E. USE IO2 = 1).  LIKEWISE, NO PHOTOLYSIS SHOULD
C     BE CALCULATED USING THE CIESLIK AND NICOLET METHOD (INO = 1), 
C     WHICH HAS BEEN MODIFIED TO AGREE WITH THE BAND INTENSITIES OF 
C     FREDERICK AND HUDSON (1979). 
C          THIS VERSION OF THE PROGRAM HAS A NEW OPTION FOR INCLUDING
C     TRANSPORT OF SPECIES THAT ONE DOES NOT WISH TO INCLUDE IN THE BIG
C     REVERSE EULER MATRIX.  THEY CAN BE SOLVED SEPARATELY USING A TRI-
C     DIAGONAL INVERSION METHOD.  S8 PARTICLES ARE TREATED THIS WAY IN
C     THIS PROGRAM BECAUSE THEY ARE VIRTUALLY NON-EXISTENT UP HIGH.  
C     THE TOP OF THE TRIDIAGONAL MATRIX IS AT GRID POINT MZ, WHICH CAN
C     BE SET EQUAL TO OR LESS THAN NZ.
C
C       THIS PROGRAM IS A ONE-DIMENSIONAL MODEL OF THE PRIMORDIAL
C     ATMOSPHERE.  THE MIXING RATIOS OF THE LONG-LIVED SPECIES  
C     ARE CALCULATED FROM THE EQUATION
C
C         DF/DT  =  (1/N)*D/DZ(KN*DF/DZ + WNF) + P/N - LF 
C-TF         DCi/DT  =  (1/rho*r^2)*D/DZ(KN*DF/DZ + WNF) + P/RHO*mi - LCi
C
C     WHERE
C     F = MASS MIXING RATIO (USOL)
C     mi = molecular weight (WGT)
C     K = EDDY DIFFUSION COEFFICIENT (EDD)
C     N = TOTAL NUMBER DENSITY (DEN)
C     RHO = TOTAL MASS DENSITY
C     L = CHEMICAL LOSS FREQUENCY (XL)
C     P = CHEMICAL PRODUCTION RATE (XP)
C     W = FALL VELOCITY (WFALL) FOR PARTICLES, POSITIVE DOWNWARD (ALL
C         FLUXES, BY CONTRAST, ARE POSITIVE UPWARD)
C
C          THIS PROGRAM HAS TWO SWITCHES FOR ALLOWING VARIABLE BOUNDARY
C     CONDITIONS.  ISULF ALLOWS THE SULFUR (SO2 AND/OR H2S) INFLUX TO 
C     VARY, AND IH2 ALLOWS THE H2 MIXING RATIO TO FLOAT IN SUCH A MAN-
C     NER AS TO BALANCE THE HYDROGEN BUDGET.  THESE FLAGS SHOULD BE  
C     TURNED OFF (SET = 0) FOR MOST CALCULATIONS; THEY ARE NOT GUARAN-
C     TEED TO WORK VERY WELL.
C
C          THE SYSTEM OF PDES IS SOLVED USING THE REVERSE EULER      
C     METHOD.  LETTING THE TIME STEP GO TO INFINITY GIVES YOU NEWTONS
C     METHOD, E.G. IT REVERTS TO AN EFFICIENT STEADY-STATE SOLVER.  
C
C          THE LIST OF SUBROUTINES IS AS FOLLOWS:
C     (1) GRID   -  SETS UP THE ALTITUDE GRID [NOT USED, DONE IN MAIN PROGRAM]
C     (2) RATES  -  DEFINES CHEMICAL REACTION RATES AND RAINOUT RATE
C     (3.2) PHOTO   - COMPUTES PHOTOLYSIS RATES (CALLS MSCAT)        
C     (3.3) O3PHOT  - COMPUTES COEFFICIENTS USED TO FIND O(1D) QUANTUM
C                     YIELDS IN O3 PHOTOLYSIS
C     (4)   DENSTY  - COMPUTES ATMOSPHERIC NUMBER DENSITIES FROM HYDRO-
C                     STATIC EQUILIBRIUM
C     (5)   DIFCO   - COMPUTES DK = K*N BETWEEN GRID POINTS; ALSO FINDS
C                     DIFFUSION LIFETIMES H*H/K
C     (5.2) SATRAT  - COMPUTES SATURATION H2O MIXING RATIOS AT ALL HEIGH
C                     FINDS MANABE/WETHERALD RH DISTRIBUTION IN TROPOSPH  
C     (6) OUTPUT -  PRINTS OUT RESULTS
C     (7) DOCHEM - DOES CHEMISTRY FOR ALL SPECIES AT ALL GRID
C                  POINTS BY CALLING CHEMPL
C     (8) CHEMPL - COMPUTES CHEMICAL PRODUCTION AND LOSS RATES
C                  FOR ONE SPECIES AT ALL GRID POINTS
C     (9) LTNING -  COMPUTES LIGHTNING PRODUCTION RATES FOR O2 AND
C                   N2 BASED ON CHAMEIDES' RESULTS
C
C          OTHER DEFINED FUNCTIONS INCLUDE:
C     (1) TBDY   -  COMPUTES 3-BODY REACTION RATES
C     (2) E1     - EXPONENTIAL INTEGRAL OF ORDER ONE
C
C ***** REACTION LIST *****
C ***** NEUTRAL REACTION LIST *****
C-TF  Nitrogen Chemistry
C     1)  N4S + N2D + M  =  N2 + M + 12.18eV
C     2)  N4S + O2 = NO  + O   + 1.4eV
C     3)  N4S + NO = N2  + O   + 2.68eV
C     4)  N2D + O  = N4S + O   + 2.38eV
C     5)  N2D + O2 = NO  + O1D + 1.84eV
C     6)  N2D + O2 = NO  + O3P + 3.76eV
C     7)  N2D + NO = N2  + O   + 5.63eV
C     8)  N2D      = N4S + HV(5200A)
C
C-TF Oxygen Chemistry  ignore O2(1Sg), O2(1Dg), etc.
C     9)  O1D + M  = O3P + M   + 1.96eV
C    10)  O1D      = O   + hv
C    11)  O1D + H2O= OH  + OH  + 1.23eV
C    12)  O   + O  + M   = O2  + M  + 5.1eV
C    13)  O   + O2 + M   = O3  + M  + 1.10eV
C    14)  O   + O3 = O2  + O2  + 4.06eV
C
C-TF Carbon Chemistry
C    15)  CO  + O  + M   = CO2 + M + 5.51eV
C    16)  C   + O2 = CO  + O   + 5.98eV
C
C-TF Hydrogen Chemistry
C    17)  H2  + O1D= H   + OH  + 1.88eV
C    18)  H2  + O  = H   + OH  + 0.08eV
C    19)  H2  + M  = H   + H   + M + -4.52eV
C    20)  H   + O2 = O   + OH  + -0.72eV
C    21)  H   + O3 = OH  + O2  + 3.34eV
C    22)  H   + H  + M   = H2  + M + 4.52eV
C    23)  H   + H2O= H2  + OH  + -0.65eV
C
C-TF OH Chemistry
C    24)  OH  + N4S= NO  + H   + 2.1eV
C    25)  OH  + O  = H   + O2  + 0.72eV
C    26)  OH  + CO = CO2 + H   + 1.07eV
C    27)  OH  + C  = CO  + H   + 6.70eV
C    28)  OH  + H2 = H2O + H   + 0.65eV
C    29)  OH  + OH = H2O + O   + 0.73eV
C    30)  OH  + H  + M   = H2O + M + 5.17eV
C    31)  OH  + H  = H2  + O   + 0.08eV
C
C ***** ION REACTION LIST *****
C-TF  Nitrogen Ion Chemistry
C    32)  N2I + O2 = O2I + N2  + 3.52eV
C    33)  N2I + O  = NOI + N2D + 0.7eV
C    34)  N2I + O  = OI  + N2  + 1.96eV
C    35)  N2I + NO = NOI + N2  + 6.25eV
C    36)  N2I + CO2= CO2I+ N2  + 1.81eV
C    37)  N2I + CO = COI + N2  + 1.57eV
C    38)  NI  + O2 = OI  + NO  + 1.28eV
C    39)  NI  + O2 = O2I + N2D + 0.1eV
C    40)  NI  + O2 = O2I + N4S + 2.486eV
C    41)  NI  + O2 = NOI + O   + 6.699eV
C    42)  NI  + O  = OI  + N4S + 0.98eV
C    43)  NI  + NO = NOI + N4S + 5.29eV
C    44)  NI  + CO2= CO2I+ N4S + 0.78eV
C    45)  NI  + CO2= COI + NO  + 1.57eV
C    46)  NI  + CO = COI + N4S + 0.54eV
C    47)  NI  + CO = NOI + C   + 0.63eV
C    48)  NI  + H  = HI  + N4S + 0.90eV
C
C-TF  Oxygen Ion Chemistry
C    49)  O2I + N2 = NOI + NO  + 0.93eV
C    50)  O2I + N4S= NOI + O   + 4.21eV
C    51)  O2I + NO = NOI + O2  + 2.813eV
C    52)  OI  + NO = NOI + O   + 4.36eV
C    53)  OI  + CO2= O2I + CO  + 1.2eV
C    54)  OI  + H2 = OHI + H   + 0.36eV
C    55)  OI  + H  = HI  + O   + 0.02eV
C    56)  OI  + N2 = NOI + N4S + 1.0888eV
C    57)  OI  + O2 = O2I + O   + 1.556eV
C    58)  OI  + N2D= NI  + O3P + 1.45eV
C    59)  OI2P+ N2 = N2I + O   + 3.02eV
C    60)  OI2P+ N2 = NI  + NO  + 0.70eV
C    61)  OI2P+ O  = OI  + O   + 5.2eV
C    62)  OI2P     = OI  + hv(2470A)
C    63)  OI2P     = OI2D+ hv(7320A)
C    64)  OI2D+ N2 = OI  + N2  + 3.31eV
C    65)  OI2D+ N2 = N2I + O   + 1.33eV
C    66)  OI2D+ O  = OI  + O   + 3.31eV
C    67)  OI2D+ O2 = O2I + O   + 4.865eV
C    68)  OI2D     = OI  + HV    (3726/3729A)
C
C-TF Hydrogen Ion Chemistry
C    69)  H2I + O  = OHI + H   + 2.17eV
C    70)  H2I + H2 = H3I + H   + 1.70eV
C    71)  H2I + H  = HI  + H2  + 1.83eV
C    72)  HI  + O  = OI  + H   + -0.02eV
C    73)  HI  + NO = NOI + H   + 4.34eV
C-TF HI+H2(ground) -> H2I + H loses 1.83eV energy. Do not know how
C-TF much energy is released from H2(v>=4) state.
C    74)  HI  + H2(v>=4) = H2I + H
C    75)  H3I + H  = H2I + H2  + -1.70eV
C
C-TF some CO2/CO related reactions from Fox 1979
C    76)  CO2I+ O  = O2I + CO  + 1.33eV		Fehsenfeld 1970
C    77)  CO2I+ O  = OI  + CO2 + 0.13eV		Fehsenfeld 1970
C    78)  CO2I+ NO = NOI + CO2 + 4.51eV		Fehsenfeld 1970
C    79)  CO2I+ H  = HI  + CO2 + 0.17eV		Nagy 1983 in Venus book
C    80)  COI + O  = OI  + CO  + 0.39eV		Nagy 1983 in Venus book
C    81)  COI + NO = NOI + CO  + 4.75eV		Nagy 1983 in Venus book
C    82)  COI + CO2= CO2I+ CO  + 0.24eV		Fehsenfeld 1966
C
C-TF Recombination reactions
C    83)  N2I + e  = N4S + N4S + 5.82eV  (10%)	Smithtro 2005
C    84)  N2I + e  = N4S + N2D + 3.44eV  (90%)	Smithtro 2005
C    85)  NI  + e  = N   + HV
C    86)  O2I + e  = O3P + O3P + 6.99eV  (22%)
C    87)  O2I + e  = O3P + O1D + 5.02eV  (42%)
C    88)  O2I + e  = O1D + O1D + 3.06eV  (36%)
C    89)  OI  + e  = O   + HV
C    90)  NOI + e  = N4S + O   + 2.75eV  (20%)
C    91)  NOI + e  = N2D + O   + 0.38eV  (80%)
C    92)  HI  + e  = H   + HV
C    93)  H2I + e  = H   + H   + 10.91eV
C    94)  H3I + e  = H2  + H   + 9.21eV
C    95)  H3I + e  = H   + H   + H + 4.69eV
C    96)  OHI + e  = O   + H   + 8.74eV
C-TF the energy associated with A98 is strange
C    97)  CO2I+ e  = CO  + O   + 4.56eV		
C    98)  COI + e  = C   + O   + 2.87eV
C
C-TF de-excitation
C    99)  N2D + e  = N4S + e   + 2.38eV
C   100)  OI2P+ e  = OI  + e   + 5.0eV
C   101)  OI2P+ e  = OI2D+ e   + 1.69eV
C   102)  OI2D+ e  = OI  + e   + 3.31eV
C
C-TF HO2, H2O2 related reactions
C   103)  H   + O2 + M = HO2 + M + 2.11eV
C   104)  HO2 + H  = H2O + O   + 2.34eV   (2%)   from JPL94     8.1e-11
C   105)  HO2 + H  = H2  + O2  + 2.41eV   (8%)   from JPL94
C   106)  HO2 + H  = OH  + OH  + 1.61eV   (90%)  from JPL94
C   107)  HO2 + OH = H2O + O2  + 3.06eV
C   108)  OH  + O3 = HO2 + O2  + 1.73eV
C   109)  HO2 + O  = OH  + O2  + 2.33eV
C   110)  HO2 + O3 = OH  + 2O2 + 1.23eV
C   111)  HO2 + HO2= H2O2+ O2  + 1.71eV
C   112)  H2O2+ OH = HO2 + H2O + 1.35eV
C   113)  H2O2+ O  = HO2 + OH  + 3.44eV
C
C-TF Photodissociation 
C   201)  N2  + HV = N4S + N2D        N2 excited to predissociation level
C   202)  O2  + HV = O3P + O1D        <1749A 
C   203)  O2  + HV = O3P + O3P        <2422A
C   204)  O3  + HV = O2  + O1D
C   205)  O3  + HV = O2  + O3P
C   206)  NO  + HV = N4S + O          <1910A
C   207)  CO2 + HV = CO  + O3P        <2247A
C   208)  CO2 + HV = CO  + O1D	      
C   209)  CO  + HV = C   + O
C   210)  H2O + HV = H   + OH
C   211)  H2O2+ HV = OH  + OH
C
C-TF Photoionization
C   301)  N2  + HV = N2I + e          <796A, include electron impact
C   302)  N2  + HV = NI  + N4S + e    dissociative ionization <510A
C   303)  N2  + HV = NI  + N2D + e
C   304)  N4S + HV = NI  + e
C   305)  O2  + HV = O2I + e          <1026A, include electron impact
C   306)  O2  + HV = OI  + O   + e    <662A
C   307)  O2  + HV = OI2P+ O   + e    
C   308)  O2  + HV = OI2D+ O   + e    
C   309)  O   + HV = OI  + e          <911A
C   310)  O   + HV = OI2P+ e
C   311)  O   + HV = OI2D+ e
C   312)  NO  + HV = NOI + e          (LyA or <1340A) 
C   313)  H2  + HV = H2I + e
C   314)  H2  + HV = HI  + H   + e
C   315)  H   + HV = HI  + e
C
C-TF ignore CO2 and CO ionization for now
C   316)  CO2 + HV = CO2I+ e
C   317)  CO  + HV = COI + e
C
C        THIS PROGRAM DOES THE CHEMISTRY AUTOMATICALLY.  THE CHEMICAL
C     REACTIONS ARE ENTERED ON DATA CARDS IN FIVE 10-DIGIT COLUMNS
C     STARTING IN COLUMN 11, I.E.                         
C
C         REAC1     REAC2     PROD1     PROD2     PROD3
C
C     THE IMPORTANT PARAMETERS DESCRIBING THE CHEMISTRY ARE
C        NR   = NUMBER OF REACTIONS
C        NSP  = NUMBER OF CHEMICAL SPECIES
C        NSP1 = NSP + 1 (INCLUDES HV)
C        NSP2 = NSP + 2 (INCLUDES "M")
C-TF     NSP3 = NSP + 3 (INCLUDES "NUL")
C-TF     NSP4 = NSP + 4 (INCLUDES "e")
C        NQ   = NUMBER OF SPECIES FOR WHICH A DIFFUSION EQUATION
C               IS SOLVED AND WHICH ARE IN THE BIG, REVERSE EULER MATR
C        NQT = TOTAL NUMBER OF SPECIES FOR WHICH TRANSPORT IS INCLUDED
C              (THOSE NOT IN THE BIG MATRIX ARE SOLVED WITH A STEADY
C              STATE TRIDIAGONAL INVERSION METHOD)
C        NMAX = MAXIMUM NUMBER OF REACTIONS IN WHICH AN INDIVIDUAL
C               SPECIES PARTICIPATES
C
C        PHOTOLYSIS REACTIONS ARE IDENTIFIED BY THE SYMBOL HV (NOT
C     COUNTED IN EVALUATING NSP).  THREE-BODY REACTIONS ARE WRITTEN
C     IN TWO-BODY FORM, SO THE DENSITY FACTOR MUST BE INCLUDED IN
C     THE RATE CONSTANT.
C        THE CHEMICAL REACTION SCHEME IS STORED IN THE FOLLOWING MATRICE
C
C     ISPEC(NSP4) = VECTOR CONTAINING THE HOLLERITH NAMES OF THE
C                  CHEMICAL SPECIES.  THE LAST 4 ENTRIES ARE HV, M, NUL and e.
C     JCHEM(5,NR) = MATRIX OF CHEMICAL REACTIONS.  THE FIRST TWO ARE
C                   REACTANTS, THE LAST THREE ARE PRODUCTS.
C     ILOSS(2,NSP,NMAX) = MATRIX OF LOSS PROCESSES.  ILOSS(1,I,L)
C                         HOLDS REACTION NUMBER J, ILOSS(2,I,L) HOLDS
C                         REACTANT NUMBER.
C     IPROD(NSP,NMAX) = MATRIX OF PRODUCTION PROCESSES.  IPROD(I,L)
C                       HOLDS REACTION NUMBER J.
C     NUML(NSP) = NUMBER OF NON-ZERO ELEMENTS FOR EACH ROW OF ILOSS
C     NUMP(NSP) = NUMBER OF NON-ZERO ELEMENTS FOR EACH ROW OF IPROD
C
      use parameters
      use flds_atmosphere
      use flds_chemistry
      use flds_transport
      use flds_species
      use flds_heating
      
CZ hydro61 added
cz      USE flds_debugging
C      
      implicit none
C      
      integer, dimension(JJ) :: lwav
      integer :: NTMP, NTMP1, NTF_STOP_FLG, NZZ1
      integer :: JCHG

CZ Added integer flags to turn printing on (1) or off (0) to screen

      INTEGER :: NQ3
      INTEGER :: IERRDEST = 1  ! 1 = when error detected, print out utden
      REAL :: PROD

      INTEGER :: IZPRNTGRID = 0   ! set = 1 to print altitudes of grid points
      integer :: IZPRNTFLG1 = 0    ! printout: PRINT 197,J,(JCHAR(M,J),M=1,5),TF_EV(J)
      INTEGER :: IZPRNTFLG2 = 0    ! controls printout of species names
      INTEGER :: IZPRNTFLG3 = 0    ! controls printout of data read in from global mean model
      Integer :: IZPRNTFLG4 = 1   ! DO NOT USE; SET = 1; chem vs. dif statements, tf_chem, tf_transport, DUSOL, etc.
      INTEGER :: IZPRNTFLG5 = 0   ! "after SUBSONIC" printout fluxes and MMRs of H-species
      INTEGER :: IZPRNTFLG6 = 0    ! dif @ J = statements - fluxes and fractions of H-containing species
      INTEGER :: IZPRNTFLG7 = 0      ! controls "WARNING: TSOL <" printouts
      INTEGER :: IZPRNTFLG8 = 0      ! controls u_inf summation printouts     
      INTEGER :: IZPRNTFLG10 = 0   ! prints J, R(J), UNUM(LO,J), DEN(J), WT(J), RHO(J), H2VOLC value
      INTEGER :: IZPRNTFLG11 = 1    ! controls tf_heating.f
      INTEGER :: IZPRNTFLG12 = 1   ! controls Output.f outputs  
      
      INTEGER :: NEWTON_DEBUG = 0

CZ Grid debugger
      INTEGER :: GRD_DBG =0
CZ Exobase debuggers
      INTEGER :: EXO_DBG = 0
      INTEGER :: EXODEBUG = 0
CZ Debug flag for exobase i_exo and i_exo_save      
      INTEGER :: IZDEBUGEXO = 0 ! Also set in subroutines adjust_exobase and find_exobase
      INTEGER :: DEBUGUINF = 0
      
CZ NaN Debugger      
      INTEGER :: NAN_DBG = 0
CZ  MEZ newton debugger    
      INTEGER :: NEWT_DBG = 0
CZ Debug flag turns on additional print statements
      INTEGER :: IZDEBUG = 0
CZ Debug flag for 
      INTEGER :: IZMMRDEBUG = 2 !=1 prints MMR, VMR of H2 every step; =2 prints at 2nd last step
CZ 
      INTEGER :: CROSDBG = 0 !debug crossover mass
CZ
      integer :: IZFLUXDEBUG = 0 ! debug flux calculation and updated chemistry; =100 prints all species; other +integer prints only that species; 0 = no print
      INTEGER :: IZHUSOLDBG = 0  ! gives warnings of USOL(H or H2) < small 

CZ Flags to control boundary conditions
CZ: If IZUBC=0, MB(I)=1, and VEFF(I) = 0; If IZUBC=1, VEFF(I) = VEFF0(I) (input in fmodules.f)
      INTEGER :: IZUBC = 0  ! Normally =1
CZ: IZLBC=0, no change; If IZLBC=1, set all lower BC to 1, constant mixing ratio
      INTEGER :: IZLBC = 1  ! Normally = 0
CZ Flag to debug boundary conditions; set = 1 to print extra info on BCs      
      integer :: IZBCDEBUG = 0 

      integer :: TEDEBUG = 0

      INTEGER :: EDEBUG = 0   ! electron number UNe

CZ Other variables added by MEZ      
      INTEGER :: IZH2X1P01,IZGRIDSTR,IZLOGGRID
     | ,IZ300,J1,J2,cretry,NZREAD,ipsimaxgrid,ipsimaxcount
      integer :: i_exo_chg = 0 !was temp. used to note that exobase had moved up or down
      INTEGER :: iexohold = 0 ! normal operation =0; set >0 to hold at level determined at step N=iexohold
      REAL :: A_EDDY,B_EDDY,TGRID1
C      
      double precision, dimension(JJ) :: wav,phflux,sigh2,sigh
     | ,effh2,effh,eflux
C      
      double precision, dimension(NZM) :: u02
     | ,fluxN,psi,alam,alamp,r0k,ujeans,RHO_ini
     | ,UNe_ini
     | ,rho_hydro,unum_h_hydro,unum_n2_hydro
     | ,vmro2,vmrn2,vmro,vmrh2,vmrh,vmrn
     | ,T_save,RHO_save,U_ADV_save,WT_save,fluxN_save,DEN_save,q_save
     | ,Te_save,Ti_save,delTe,delT
     | ,PRES_MKS
     | ,TA,TB,TC,TY,TSOL
     | ,Te_INI,Ti_INI,T_INI
CZ Added by MEZ     
     | ,zl,R1,zk1,fluxM,fluxM_save,UOLDH,UOLDH2,zkb,zlb
     | , RratioL,RratioU,RHOR2,RHOR2_M
CZ hydro58     
cz     | ,R_M,RHO_M,DEN_M  ! defined in fmodules
    
     
      double precision :: DJACN(NAQ,NAQ),RHSN(NAQ)
     | ,IPVTN(NAQ),F(NAQ),FP(NAQ),UNUM1(NSP4),
     | UNUM1_SAVE(NSP4)
     
      double precision :: FVAL(NQ,NZM),FV(NQ,NZM)
     1 ,DJAC(LDA,NEQ),RHS(NEQ),IPVT(NEQ),PTBR(NZM)

      double precision :: USOL_OLD(13,NZM),USOL_TRI_OLD(1,NZM)
     | ,USH_OLD(16,NZM),UNUM_OLD(30,NZM)

CZ Added by MEZ
      DOUBLE PRECISION :: fvH,fvH2,fvH2O,fvtot,fvtot1,HfluxNcomp,
     | delUSOL, delUNUM,UNUMH2min,UNUMHmin,Tmin3,psimax,
     | u_inf_ini, u_inf_adj, delU_ADV, delZold, delZnew,zk0,zl0,RL0
     | ,RBND0,zlb0,zkb0
     
      double precision, dimension(NZM) :: Tz
      DOUBLE PRECISION, DIMENSION(NSP,NZM) :: USAVEUN,USAVEUN1,USAVE1
      logical, dimension(NZM) :: MASK1
      double precision, dimension(NQ) :: SGFLUX
C
C ***** SET MODEL PARAMETERS *****
C     ZY = SOLAR ZENITH ANGLE (IN DEGREES)
C     AGL = DIURNAL AVERAGING FACTOR FOR PHOTORATES
C     LTIMES = COUNTER FOR PHOTORATE SUBROUTINE
C     ISEASON = TELLS WHETHER P AND T VARY WITH TIME (THEY DON'T FOR
C               ISEASON < 3)
C     IZYO2 = TELLS WHETHER SOLAR ZENITH ANGLE VARIES WITH TIME (0 SAYS
C             IT DOESN'T; 1 SAYS IT DOES)
C     IO2 = 0 FOR ALLEN AND FREDERICK O2 SCHUMANN-RUNGE COEFFICIENTS
C         = 1 FOR EXPONENTIAL SUM FITS (FOR LOW-O2 ATMOSPHERES)
C     INO = 0 FOR ALLEN AND FREDERICK NO PREDISSOCIATION COEFFICIENTS
C         = 1 FOR MODIFIED CIESLIK AND NICOLET FORMULATION
C     EPSJ = AMOUNT BY WHICH TO PERTURB THINGS FOR JACOBIAN CALCULATION
C     ZTROP = TROPOPAUSE HEIGHT (ABOVE WHICH H2O BEHAVES AS A NONCONDENS
C             ABLE GAS)
C     PRONO  = COLUMN-INTEGRATED NO PRODUCTION RATE FROM LIGHTNING IN
C              PRESENT ATMOSPHERE
C     VOLFLX = ESTIMATED VOLCANIC OUTGASSING RATE OF SO2 PLUS H2S
C     H2VOLC = ESTIMATED VOLCANIC OUTGASSING RATE OF H2
C     ISULF  = 0 FOR CONSTANT MIXING RATIOS OR CONSTANT FLUXES OF
C              SO2 AND H2S, 1 FOR VARIABLE FLUXES (OPTION 1 MAKES
C              THE TOTAL SULFUR LOSS EQUAL TO VOLFLX, INCLUDING
C              SURFACE DEPOSITION OF SO2 AND H2S)
C     IH2    = 0 FOR CONSTANT H2 MIXING RATIO, 1 FOR MIXING RATIO
C              THAT VARIES SO AS TO BALANCE THE H2 BUDGET
C     DT = INITIAL TIME STEP
C     TSTOP = TIME AT WHICH CALCULATION IS TO STOP
C     NSTEPS = NUMBER OF TIME STEPS TO RUN (IF TSTOP IS NOT REACHED)
C
      double precision :: ZY = 60.
      double precision :: AGL = 0.5
      double precision :: LTIMES = 0
      integer :: ISEASON = 1
      integer :: IZYO2 = 0
      integer :: IO2 = 1
      integer :: INO = 1
      double precision :: EPSJ = 1.E-4
      double precision :: PRONO = 1.E9
      double precision :: VOLFLX = 3.E9
      integer :: ISULF = 0
      integer :: IH2 = 0
CZ Added alternative IH2 = 1, this program probably doesn't use the variable IH2 anyway
C      integer :: IH2 = 1
C	
      double precision :: rootpi, xx
      double precision :: DT = 1.E-8
      double precision :: DTINV
      double precision :: TIME = 0.
      double precision :: TSTOP = 1.E14

CZ  MOVED value to Inputs.in
      integer :: NSTEPS

CZ Variables added by MEZ
      INTEGER :: IZ300A = 0
      INTEGER :: NZMAX = 600
      INTEGER :: IZDEF01 = 1  ! If = 1, eliminates divide by 10 on timestep
      INTEGER :: i_exo_max = 401
      INTEGER :: IF107  ! integer version
      INTEGER :: tempkluge = 0   ! set to 1 to limit temp excursions to +/-10% per step
      INTEGER :: chemkluge = 0   ! set to 1 to force chemistry code to converge!
      integer :: SMOOTH3 = 1    ! set to 1 to do 3-point smoothing of H and H2
CZ2
      INTEGER :: i_exo_savez 
      
      REAL :: NEWTMAX = 1.e5
      REAL :: psihi     ! upper limit on psi; set negative to turn off
      REAL :: ZGRIDFAC  ! divisor which sets grid spacing; now read in thru Inputs05.in and later
CZ Limits on change in EMAX, EMAXK      
      real :: EMAXMX = 100
      REAL :: EMAXKMX = 100      
      REAL :: DTMX = 1.E2   ! Max size of timesteps (was 1E15, later 1E12)
      
      REAL :: CO2MMR  ! Inputs
      REAL :: H2MMR    ! Inputs
      REAL :: N2MMR    ! Inputs (or calc from remainer if set = -1)
      REAL :: O2MMR    ! Inputs
      REAL :: COMMR  ! set based on planet
CZ2    
      REAL :: OMMR = 0.01    ! TF: 0.01;    Walker: 0.02 @ 100km
      REAL :: HMMR = 2.46E-8   ! TF: 2.46E-8; Walker: 2.0E-7 @ 120km
      REAL :: HEMMR = 1.29E-6 ! TF: 1.29E-6; Walker: 6.6E-7 @ 120km
     
      REAL :: PLANETNO
      INTEGER :: MARS
      INTEGER :: EARTHMARS = 0 ! 1 =Mars except Earth mass & radius; 2= Mars exc Earth mass, rad, G
      INTEGER :: EARTH
      INTEGER :: CONVEM = 0 ! CONVERT EARTH TO MARS (Set MARS=1)

CZ CONVEM=1: G, bigM, and R_PLT all start at Earth mass and radius, change to Mars
CZ CONVEM=2: bigM and R_PLT start at Earth values
CZ CONVEM=3: only R_PLT starts with Earth value
CZ CONVEM=4: only bigM starts with Earth value

       double precision :: bigM
     
C
C   PRINT OUT RESULTS EVERY NPR TIME STEPS
      integer :: N, NPR, MP=1, NN=0, NZTF_FLG=0, NQTF_FLG=0
     | ,MTF
C      
      real :: PRN,PM,u_inf,x,glbmeanread,restart_flg,read_jmatrix,x_tf
     | ,uni_grid,DZ_TF,read_sub_flg,fh2,phidif,tinf,emax,SM,BOVERH2
     | ,BOVERH2O,fh_dif,fh2_dif,fh2o_dif,dif_limited_flux,crossoverm
     | ,XP,XL,xs,dx,EREL,UMAX,glbmean_heating_flg,VDEP(NQ),u02i,ve
     | ,veu,veu2,alampi,ujeansi,flux_jeans,phie,dusol,dusol2,zmax
     | ,usol_tot,usol_tot2,dt_save,emaxk,chg,dt0,RL,dif,tmp_tf
     | ,tmp_tf1,phi_counterg_h,phi_counterg_h2,phi_counterg_h2o
     | ,tf_damping,TFM,TFM_STEPS=1000.
C
      double precision, dimension(NZM) :: qntot_glb, qnc_glb, qsrc_glb
     | , qsrb_glb, qn1d_glb, xlen_glb,qjoul_glb, qnlya_glb,qnspe_glb
     | , xnoc_glb, xco2c_glb, xo3pc_glb, xircool_glb
C      
C
      INTEGER :: I,jfc,J,IT,JTROP,M,K,IDO,II,IERR,IS,JS,MS
     | ,INEWT,NH,NQ1,KK,INN,IFIRST,L,INFO,NRAIN,LTEST
     | ,i_exo,MZ,MZ1,LB,MB,JJJ,nflag,J_TF,imaxk,i_exo_flg,IROW,LR
     | ,K1,K2,ISKIP,ITF_SFLX_FLG,iread_sub_flg

CZ hydro64 added to check flux calc

      double precision, dimension(NSP,NZM) :: FT1,FT3A,FT3B 
     | ,FT3C,FT3,FLUXKZ1


      NTF_STOP_FLG = 0
C      

      NZ1=NZ-1
      NZ2=NZ1-1
CZ 
      UMAX = 0.0
CZ      
      rootpi=sqrt(pi)
      DTINV = 1./DT
C      
C   DO PHOTORATES EVERY MP TIME STEPS
      NPHOT = 0

      PM = MP
C
      H2VOLC = 2.e10
      IF (IZPRNTFLG10.eq. 1) PRINT 515, H2VOLC
 515  FORMAT(6X,'H2VOLC = ',1PE10.3)
C
      DO I=1,NSP
      WGT(I) = WGTM(I)*protonM
      ENDDO
C
C      u_inf = 2.2e4
       u_inf = SMALL
C			      
C ****************************************************************
C Input files
      OPEN(UNIT=310,FILE='InputData/PHOTO_DATA.CH4')
CZ      
      OPEN(UNIT=710,FILE='restart.in')
      OPEN(UNIT=910,FILE='InputData/CHEM.DAT.SUA')
      OPEN(UNIT=100,FILE='InputData/EUV.DAT')
C-TF  INPUT FILES NEEDED FOR THE SUBSONIC MODEL
      open(unit=104,file='InputData/kasting2.dat',status='old')
C      
      OPEN(UNIT=110,FILE='InputData/FARUV.DAT')
      OPEN(UNIT=111,FILE='InputData/co2cool_x_t_fi_col.dat')
      OPEN(UNIT=112,FILE='InputData/uars_solaruv_flx_fsig.dat')
      OPEN(UNIT=113,FILE='InputData/uars_solaruv_flx_fsigt_o3.dat')
      OPEN(UNIT=114,FILE='InputData/uars_solaruv_flx_fsigt_o3-1d.dat')
C
C================================================================================
C
C Output files
      OPEN(UNIT=810,FILE='restart.out')
      OPEN(UNIT=140,FILE='rxtot.out.dat')
      OPEN(UNIT=150,FILE='int.rates.out.dat')
      OPEN(UNIT=220,FILE='tautot_v_wav.out.dat')
      OPEN(UNIT=240,FILE='freq.out.dat')
      OPEN(UNIT=79,FILE='tpf_absor_spect.out.dat')
      OPEN(UNIT=80,FILE='tf_profiles.dat')
      OPEN(UNIT=82,FILE='tf_transport.dat')
      OPEN(UNIT=83,FILE='tf_chem.dat')
      OPEN(UNIT=85,FILE='single_species.dat')
C      
      OPEN(UNIT=850,FILE='time_evol.dat')
C      
C-TF  OUTPUT FILES FOR SUBSONIC MODEL      
      open(unit=101,file='subsonic.out')
      open(unit=103,file='utden.out')

C ****************************************************************
CZ Additional files

CZ Output file for column densities of species
      OPEN(UNIT=6904,FILE='idlpass_sp.out')
CZ Output file for status of run
      OPEN(UNIT=6942,FILE='idlpass_status.out')
CZ Outputs
      OPEN(UNIT=6954,FILE='idlpass_misc.out')

      OPEN(UNIT=6999,FILE='Inputs10.in')
CZ Read in Inputs and Function Flags to control coupled_sua functionality
      READ(6999,*) H2MMR,CO2MMR,N2MMR,O2MMR,A_EDDY,B_EDDY,IZH2X1P01,
     | IF107,IZLOGGRID,IZ300,TGRID1,NSTEPS,ZGRIDFAC,PLANETNO,u_inf_ini
CZ
       f107 = IF107
       f107a = f107
CZ2
      if (u_inf_ini .GE. 0) then
         u_inf = u_inf_ini
         print*,"u_inf_ini>0, using u_inf to control to constant psi"
      endif	 

      if (u_inf_ini .LT. 0.) THEN
         PRINT*,"u_inf_ini<0, U_ADV(NZ) calculated using Jeans escape"
      ENDIF
CZ2

CZ hydro62dr
cz      do J=1,NZM
cz       DRX(J)=1.0
cz      ENDDO
      

      IF(DEBUGUINF .EQ. 1) PRINT*,"after Read, u_inf =",u_inf
       
       IF (PLANETNO.EQ.1.0 .OR. PLANETNO.EQ.1.1) THEN ! Earth
          EARTH = 1
          MARS = 0
       ENDIF
       
       IF (PLANETNO .EQ. 2.0) THEN ! Early Mars
         EARTH = 0
         MARS = 1
       ENDIF
       
       IF (PLANETNO.EQ.1.0 .OR. PLANETNO.EQ.2.0) THEN ! Early H-rich planet
         COMMR = 1.E-3
         psihi = 0.6
CZ
        IF (IZLBC .EQ. 1) THEN
           DO I = 1,NQ
             LBOUNDERLY(I) = 1
           ENDDO
        ENDIF   
CZ	 
         DO 6990 I = 1,NQ
           VDEP(I)=VDEPERLY(I)
           SGFLUX(I)=SGFLUXERLY(I)
          LBOUND1(I)=LBOUNDERLY(I)
 6990 CONTINUE	  
      ENDIF

       IF (PLANETNO.EQ.1.1) THEN ! Modern Earth
         COMMR = 1.36E-5  ! TF: 1.36E-5;   JK: 3.E-5
         OMMR = 0.01    ! TF: 0.01;    Walker: 0.02 @ 100km
         HMMR = 2.46E-8   ! TF: 2.46E-8; Walker: 2.0E-7 @ 120km	 
CZ2	 
         psihi = 0.001
CZ
        IF (IZLBC .EQ. 1) THEN
           DO I = 1,NQ
             LBOUNDMOD(I) = 1
           ENDDO
        ENDIF   
CZ 
         DO 6991 I = 1,NQ
           VDEP(I)=VDEPMOD(I)
           SGFLUX(I)=SGFLUXMOD(I)
          LBOUND1(I)=LBOUNDMOD(I)
 6991 CONTINUE
       ENDIF
       
       if (u_inf_ini .ge. 0) PRINT*,"psihi set to ",psihi
       if (u_inf_ini .lt. 0) print*,"Jeans mode on, psihi not used"

CZ**************
CZ More input files (moved after reading Inputs10.in)

      IF (EARTH .EQ. 1) THEN
      OPEN(UNIT=410,FILE='InputData/PLANET.EARTH')

       bigM =5.975E27    ! mass of the Earth (gram)
      ENDIF
      
cz      IF (MARS .EQ. 1) OPEN(UNIT=410,FILE='InputData/PLANET.MARS')      
      IF (MARS .EQ. 1) THEN
        G = 373.0 ! Earth G = 980.7
        bigM =6.419E26    ! mass of Mars (gram)
        R_PLT = 3.396E8 !radius of Mars=3396km
        FSCALE = 0.43 !Solar Flux Scaling - Earth=1.0, Mars=0.43
        ALB = 0.215 !Surface Albedo - Earth=0.25, Mars=0.215
cz        DELZ = 2.0E5 !Height of each atm layer - Earth=1.0E5, Mars=2.0E5 (cm)
        DELZ = 1.0E5 !Height of each atm layer - Earth=1.0E5, Mars=2.0E5 (cm)	
        ZTROP = 1.5E6 !Height of the tropopause - typically Earth=11km, Mars=15km
        JTROP = 11
      endif
CZ      
      IF (N2MMR .LT. 0.) THEN
        N2MMR=1.-H2MMR-HMMR-CO2MMR-COMMR-O2MMR-OMMR-HEMMR-0.005
        print*,"N2MMR set to",N2MMR
      ENDIF	
CZ2

CZ*************
       
      PRINT*,"H2MMR =",H2MMR
      PRINT*,"CO2MMR =",CO2MMR
      PRINT*,"N2MMR =",N2MMR
      PRINT*,"O2MMR =",O2MMR
      PRINT*,"COMMR =",COMMR
CZ2
      print*,"OMMR =",OMMR
      PRINT*,"HMMR =",HMMR
      PRINT*,"HEMMR =",HEMMR
      
      PRINT*,"LWR BDRY:"
      print 1119, LBOUND1
      PRINT*,"  O  O2  N2  He   H  H2  CO2  CO N4S NO  H2O  O+  N+  H+"
      PRINT*,"0=Const Vdep;1=Const MR; 2=Const Flux"
      PRINT*,""      
      PRINT*,"A_EDDY =",A_EDDY
      PRINT*,"B_EDDY =",B_EDDY
      PRINT*,"IZH2X1P01 =",IZH2X1P01
      PRINT*,"f107 =",f107
      PRINT*,"IZLOGGRID =",IZLOGGRID
      PRINT*,"IZ300 = ",IZ300
      PRINT*,"TGRID1 = ",TGRID1
      PRINT*,"NSTEPS = ",NSTEPS
      PRINT*, "ZGRIDFAC =",ZGRIDFAC
      print*,"PLANETNO = ",PLANETNO
      IF (PLANETNO .EQ. 1.0)PRINT*,"Early Earth"
      IF (PLANETNO .EQ. 1.1)PRINT*,"Modern Earth"
      IF (PLANETNO .EQ. 2.0)PRINT*,"Early Mars"      
CZ2      
      PRINT*, "u_inf_ini =",u_inf_ini
      PRINT*,""

 1119 format(14I4)
 
      NPR = NSTEPS

      IF (PLANETNO .EQ. 1.0)PRINT*,""

CZ
      WRITE(850,3213) NSTEPS,NZM,NSP,TFM_STEPS
 3213 FORMAT(3I5,1PE10.2)
C
      print*, "in sua.f:0"
C-TF  READ IN DATA FOR Formichev CO2 cooling COMPUTATION
      jfc=0
      do j=1,50
      read(111,1112) jfc,atm_formichev(j,1),atm_formichev(j,2),
     | atm_formichev(j,3)
      atm_formichev(j,4)=0.
      atm_formichev(j,5)=0.
      atm_formichev(j,6)=0.
      atm_formichev(j,7)=0.
      atm_formichev(j,8)=0.
      atm_formichev(j,9)=0.
      atm_formichev(j,10)=0.
 1112 format(4x,I2,8x,F5.2,6x,F5.1,5x,1PE8.2)     
C      print 1112, jfc,atm_formichev(j,1),atm_formichev(j,2),
C     | atm_formichev(j,3) 
      enddo
      do j=51,67
      read(111,1113) jfc,atm_formichev(j,1),atm_formichev(j,2),
     | atm_formichev(j,3),atm_formichev(j,4),atm_formichev(j,5),
     | atm_formichev(j,6),atm_formichev(j,7),atm_formichev(j,8),
     | atm_formichev(j,9),atm_formichev(j,10)
 1113 format(4x,I2,8x,F5.2,6x,F5.1,5x,1PE8.2
     | ,4x,1PE8.2,3x,1PE8.2,4x,1PE8.2,4x,1PE8.2,4x,1PE8.2,5x,F5.2
     | ,4x,1PE8.2)
C      print 1113, jfc,atm_formichev(j,1),atm_formichev(j,2),
C     | atm_formichev(j,3),atm_formichev(j,4),atm_formichev(j,5),
C     | atm_formichev(j,6),atm_formichev(j,7),atm_formichev(j,8),
C     | atm_formichev(j,9),atm_formichev(j,10)
      enddo
      close(111)
C-TF  END READING DATA FOR Formichev CO2 cooling

C-TF  READ IN DATA FOR UARS_SOLARUV (Guy Brassier data)
C      print*, "reading in data for UARS_SOLARUV"
      do j=1,170
C-TF  fsig(:,1:3) are the absorption cross section for O2, H2O, and CO2      
      read(112,1114), jfc,wv_low(j),wv_high(j),sminflx(j),smedflx(j)
     | ,smaxflx(j),x,fsig(j,1),fsig(j,2),fsig(j,3)
 1114 format(4x,I3,4x,1PE10.4,8(1x,1PE10.4))
C      print 1114, jfc,wv_low(j),wv_high(j),sminflx(j),smedflx(j)
C     | ,smaxflx(j),x,fsig(j,1),fsig(j,2),fsig(j,3)
      enddo
C
 1115 format(4x,I3,5x,1PE8.2,5(3x,1PE8.2))
      read(113,1115), jfc, (sigtn(it),it=1,6)
C      print 1115, jfc, (sigtn(it),it=1,6)
      do j=1,170
C-TF  sigtn      
      read(113,1115), jfc, (fsigt(j,1,it),it=1,6)
C      print 1115, jfc, (fsigt(j,1,it),it=1,6)
      enddo
      close(113)
C      
      read(114,1115), jfc, (sigtn(it),it=1,6)
      do j=1,96
      read(114,1115), jfc, (fsigt(j,2,it),it=1,6)
C      print 1115, jfc, (fsigt(j,2,it),it=1,6)
      enddo
      do j=97,170
      do it=1,6
      fsigt(j,2,it)=0.
      enddo
      enddo
      close(114)
C-TF  END READING DATA FOR UARS_SOLARUV
C
C   Read the EUV data file
      read(104,4001)
 4001 format(/)
      do 4000 j=1,JJ
      read(104,4002) lwav(j),phflux(j),sigh2(j),sigh(j),effh2(j),
     |   effh(j),eflux(j)
C      eflux(J) = eflux(J)*10
 4002 format(1x,i3,6(1x,e9.2))
      wav(j) = lwav(j)
C      write(101,4002) lwav(j),phflux(j),sigh2(j),sigh(j),effh2(j),
C     2   effh(j),eflux(j)
C      PRINT 4002, lwav(j),phflux(j),sigh2(j),sigh(j),effh2(j),
C     2   effh(j),eflux(j)
 4000 continue
      close(104)
C-TF  END READING EUV DATA FILE (OLD FOR LOWER ATMOSPHERE) 
C
C
C ***** READ THE PLANET PARAMETER DATAFILE *****
      if (EARTH .EQ. 1) THEN
      READ(410,102) G,R_PLT,FSCALE,ALB,DELZ,ZTROP,JTROP
      ENDIF
 102  FORMAT(F5.1/,E7.3/,F4.2/,F5.3/,E5.1/,E5.1/,I2)
      IF (MARS .EQ. 1 ) THEN
         PRINT*, "Opening input data for Mars"
      endif
      IF (EARTH .EQ. 1 ) THEN
         PRINT*, "Opening input data for Earth"
      endif      
        PRINT*, "G = ",G
        PRINT*, "R_PLT = ",R_PLT
        PRINT*,"bigM = ",bigM
        PRINT*, "FSCALE = ",FSCALE
        PRINT*, "ALB = ",ALB
        PRINT*, "DELZ = ",DELZ
        PRINT*, "ZTROP = ",ZTROP
        PRINT*, "JTROP = ",JTROP

C ***** READ THE CHEMISTRY DATA CARDS *****
C
C-TF  IF RUNNING ON NCAR MACHINES, USE THE FOLLOWING
C      READ(910,200) JCHAR
C      PRINT 201,(J,(JCHAR(M,J),M=1,5),J=1,NR)
C
C-TF  IF RUNNING ON URANUS, USE THE FOLLOWING
C      READ(910,200) JCHEM
C      PRINT 201,(J,(JCHEM(M,J),M=1,5),J=1,NR)

 199  FORMAT(10X,A4,2X,A3,2X,A4,2X,A4,1X,A1,f5.2)
 197  FORMAT(1X,I3,1H,5X,A4,4H +  ,A3,7H  =    ,A4,4H +  ,A4,4X,A1,f5.2)
      DO 198 J=1,NR
C-TF  IF RUNNING ON NCAR MACHINES, USE THE FOLLOWING
      READ(910,199) (JCHAR(M,J),M=1,5),TF_EV(J)
      
CZ  Changed print to be conditional on flag=1    
      IF (IZPRNTFLG1.eq. 1) 
     | PRINT 197,J,(JCHAR(M,J),M=1,5),TF_EV(J)
C
C-TF  IF RUNNING ON URANUS, USE THE FOLLOWING
C      READ(910,199) (JCHEM(M,J),M=1,5),TF_EV(J)
C      PRINT 197,J,(JCHEM(M,J),M=1,5),TF_EV(J)
C
 198  CONTINUE
C
C-TF 200  FORMAT(10X,A8,2X,A8,2X,A8,2X,A8,2X,A8)
C-TF 201  FORMAT(1X,I3,1H),5X,A8,4H +  ,A8,7H  =    ,A8,4H +  ,A8,4X,A8)
 200  FORMAT(10X,A4,2X,A3,2X,A4,2X,A4,1X,A1)
 201  FORMAT(1X,I3,1H),5X,A4,4H +  ,A3,7H  =    ,A4,4H +  ,A4,4X,A1)
 
         IF (IZPRNTFLG10.eq. 1) PRINT 202,NQ,NZ,NAQ
 202  FORMAT(//1X,"NQ=",I2,5X,"NZ=",I3,5X,"NAQ=",I7)
C ***** REPLACE HOLLERITH LABELS WITH SPECIES NUMBERS IN JCHEM *****
      DO 5 J=1,NR
      DO 5 M=1,5
C      
C-TF  IF RUNNING ON NCAR MACHINES, USE THE FOLLOWING
        IF(JCHAR(M,J).EQ.' ') GO TO 5
        DO 6 I=1,NSP+4
        IF(JCHAR(M,J).NE.ISPEC(I)) GO TO 6
C	
C-TF  IF RUNNING ON URANUS, USE THE FOLLOWING
C        IF(JCHEM(M,J).EQ.1H ) GO TO 5
C        DO 6 I=1,NSP+4
C        IF(JCHEM(M,J).NE.ISPEC(I)) GO TO 6
C
        JCHEM(M,J) = I
        GO TO 5
   6    CONTINUE
      IERR = J
C      
C-TF  IF RUNNING ON NCAR MACHINES, USE THE FOLLOWING
      IF (IZPRNTFLG2.eq. 1) 
     | print*, M,J,JCHAR(M,J),I,ISPEC(I)
C-TF  IF RUNNING ON URANUS, USE THE FOLLOWING
C      print*, M,J,JCHEM(M,J),I,ISPEC(I)
C
      GO TO 25
   5  CONTINUE
C
C Read character array for P&L tables, "int.rates.out.dat"
      REWIND 910
      READ(910,200) CHEM
C     WRITE(150,201)(J,(CHEM(M,J),M=1,5),J=1,NR)
C
C ***** FILL UP CHEMICAL PRODUCTION AND LOSS MATRICES *****
      DO 7 M=1,2
      N = 3-M
      DO 7 J=1,NR
      I = JCHEM(M,J)
      IF(I.LT.1.OR.I.GT.NSP) GO TO 7
      NUML(I) = NUML(I) + 1
      IF(NUML(I).GT.NMAX) GO TO 20
      K = NUML(I)
      ILOSS(1,I,K) = J
      ILOSS(2,I,K) = JCHEM(N,J)
   7  CONTINUE
C
      DO 8 M=3,5
      DO 8 J=1,NR
      I = JCHEM(M,J)
      IF(I.LT.1.OR.I.GT.NSP) GO TO 8
      NUMP(I) = NUMP(I) + 1
      IF(NUMP(I).GT.NMAX) GO TO 20
      K = NUMP(I)
      IPROD(I,K) = J
   8  CONTINUE
C
CZ MEZ modified boundary conditions (note this is inside global mean loop)
      if (IZUBC .EQ. 1) THEN
      DO K=1,NQ
         VEFF(K) = VEFF0(K)
      enddo
      ENDIF

      if (IZUBC .EQ. 0) THEN
      DO K=1,NQ
         MBOUND(K)=0 !0=constant effusion velocity      
         VEFF(K) = 0.
      enddo
      ENDIF

      DO 40 K=1,NQ
  40    VEFF(K) = VEFF0(K)
  
      DO 402 J=1,NZ
      DO 411 I=1,NSP
      USOL(I,J)=SMALL
CZ      
      UNUM1_SAVE(I)=1.
CZ      
  411 CONTINUE
  402 CONTINUE

C
C ***** READ THE INPUT DATAFILE *****
C
      glbmean_heating_flg = 1.
      glbmeanread  = 0.
      restart_flg  = 1.
      read_sub_flg = 0.
CZ Added new integer flag to start from restart (UNIT=710) rather than NCAR results (UNIT=58)      
      iread_sub_flg = 0
      read_Jmatrix = 0.
C-TF

C-TF  READ FROM RESTART FILE****************************
      IF (restart_flg .eq. 1) THEN
CZ modified
C-TF  READ from restart FILE: restart.in
      READ(710,*) NZ,NZ1,NZ2,NZ_INI,R0
CZ
       NZREAD = NZ
       
      NZ_INI=57  ! SET  = 57  for Feng's mode       
CZ
      PRINT*, "READING DATA FROM RESTART.IN FILE"

CZ ?? Need to drop many of these parameters from restart.in, and replace UNUM with USOL
      DO J=1,NZ
      READ(710,*) R(J),xx,RHO(J),T(J),Ti(J),Te(J),EDD(J)
     | ,U_ADV(J),RBND(J),xx,DRC(J),drp(J),drm(J),WT(J),UNe(J)
     | ,(UNUM(I,J),I=1,NSP4)
      ENDDO ! End reading data from restart.in
CZ 
      if(EDEBUG.GT.0)PRINT*,"UNe=",(UNe(J),J=1,NZ)

      IF (IZDEBUG.EQ.1) PRINT*, "UNUM WAS READ FROM RESTART FILE"
CZ MEZ Added per JK, replaces varying EDD w/ 1.0e6 or 7.5e5
      DO 5773 J=1,NZM
         EDD(J) = 1.0E6
cz             EDD(J) = 7.5E5
 5773 CONTINUE

CZ hydro61 added to remove extra layer from earlier restart.in files
      J=NZ
      UNe(J)=0
      DO I=1,NSP4
      UNUM(I,J)=0
      ENDDO

      if(EDEBUG.GT.0)PRINT*,"UNe=",(UNe(J),J=1,NZ)

CZ MEZ added input temperature at first grid point from Inputsxx.in file
      T(1) = TGRID1
      JCHG = 0
      do J=1,10
            IF (T(J).LT.T(1)) THEN
            T(J) = T(1)
            JCHG = J
            ENDIF
      ENDDO
      
      PRINT*," "
      PRINT*,"T(1) SET TO ",T(1)
      IF (JCHG.GT.0) PRINT 7971,JCHG,T(1)
 7971 FORMAT("TEMPS UP TO GRID PT. ",I3," SET TO ",1P1E10.1)      
      PRINT*," "

CZ  use JK grid as boundaries
         DO J=1,600
           zl(J) = 2. + (J-1)/ZGRIDFAC 
           zk(J) = 10.**zl(J)
           R(J) = 1.e5*(zk(J)) + R_PLT
           GZ(J) = G * (R_PLT/R(J))*(R_PLT/R(J))
           r2(J) = R(J)*R(J)
         ENDDO  ! end loop to set radii of grid boundaries

CZ hydro56, hydro58
        drm(1) = 0
        DO J=1,599
         DR(J) = (R(J+1)-R(J))  ! size of layer, boundary to boundary
         drp(J) = DR(J)
         R_M(J) = (R(J+1)+R(J))/2	 
        ENDDO
	
CZ hydro56
        DO J=2,598
         drm(J) = DR(J-1) ! hydro56
        ENDDO

       DO J=1,598
         DRC(J) = (DR(J+1)+DR(J))/2.      ! distance between centers of grids
       ENDDO
	
           if (IZPRNTGRID.EQ.1) THEN
            DO J=1,99

            PRINT 6701,J,R(J)/1E5,(R(J)-R_PLT)/1E5,J,DR(J)
            ENDDO
           ENDIF	


 6701 FORMAT("   R(",I2,") = ",1P1E13.6," km,  Alt = ",
     | E10.3," km, DR(",I2,") =",1P1E10.3)

CZ       
      IF (NAN_DBG.EQ.1) THEN
         PROD = UNUM(17,1)*UNUM(18,1)*UNUM(19,1)
         PRINT*,"1126 PROD =",PROD
      ENDIF
CZ      
	 
      IF (EXO_DBG.EQ.1) THEN
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF

CZ hydro58: fixed the boundary/midpoint issues; WT is at bdry, WT_M is at midpoint

CZ Calc WTP and DEN from RHO, WT which are read in from restart.in (RHO later calc in subsonic)
      DO J=1,NZ
        WTP(J) = WT(J)/protonM      
        DEN(J) = RHO(J)/WT(J)
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0090, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
      ENDDO

cz      PRINT*,"at 1314 DEN(NZ-1)=",DEN(NZ-1)
cz      PRINT*,"at 1314 DEN(NZ)=",DEN(NZ)

CZ calculate midpoint RHO, DEN, and WT
      DO J=1,NZ-1 
        WT_M(J) = SQRT(WT(J+1)*WT(J))  
        WTP_M(J) = WT_M(J)/protonM
        RHO_M(J) = SQRT(RHO(J)*RHO(J+1))
        DEN_M(J) = RHO_M(J)/WT_M(J)
       enddo


CZ end hydro58 boundary/midpoint fixes       

      DO J=1,NZ
      
      HSCALE(J) = bol_k/(WT(J)*GZ(J))*T(J)
C      
      Te_INI(J) = Te(J)
      Ti_INI(J) = Ti(J)
      T_INI(J) = T(J)
      RHO_INI(J) = RHO(J)
      UNUM_INI(1:NSP4,J)=UNUM(1:NSP4,J)
      UNe_INI(J)=UNe(J)
C
      IF (IZPRNTFLG10.eq.1) THEN 
         PRINT 3000, J, R(J),
     |   UNUM(LO,J), DEN_M(J), WT(J), RHO(J)
          PRINT*,"J,zk(J)",J,zk(J)
          PRINT*,"DEN_M(J),RHO_INI(J)",DEN_M(J),RHO_INI(J)       
      endif	  
C      
      ENDDO
      
C
C-TF   Converting to mass mixing ratios
Cz                1 -- CONVERT UNUM INTO USOL
       CALL NORM_MMR(1)
C       
      ENDIF      ! END READING FROM RESTART FILE**********************
C
CZ hydro56: removed these, redefined drp, drm above
cz       DO J=NZ+1,NZM 
cz       DRC(J) = DRC(NZ)
cz       DR(J) = DR(NZ)
cz       drp(J) = drp(NZ) 
cz       drm(J) = drm(NZ)
CZ
cz       r2(J) = R(J)**2
cz       GZ(J) = G * (R_PLT/R(J))**2
cz       ENDDO
       
C-TF  CHECK TO SEE IF THE INPUT IS IN HYDROSTATIC EQUILIBRIUM
      PRINT*, ''
      PRINT*, "CHECK TO SEE IF THE INPUT IS IN HYDROSTATIC EQUILIBRIUM:"
      PRINT*, ''
C
      VEFF(LH) = 1E3
CZ      
      if (IZUBC .EQ. 0) VEFF(LH) = 0
CZ        
cz      u_inf = VEFF(LH)*USOL(LH,NZ)  ! set upper boundary bulk motion V
      u_inf = SMALL
      N = 0
CZ added fluxM to call statement and subsonic (cs_n11.f)
CZ replaced u_inf w/ u_fac in call statement; subsonic will now use wind based on u_fac and ujeans (cs_r13.f)
CZ back to u_inf but now under direct control

      CALL SUBSONIC(N,u_inf,rootpi,ujeans,fluxN,fluxM,psi,bigM)
CZ
      IF(DEBUGUINF .EQ. 1) 
     |   PRINT*,"returning from first call to subsonic, u_inf =",u_inf

CZ hydro64: recalculate midpoint RHO, DEN, and WT
      DO J=1,NZ-1 
        WT_M(J) = SQRT(WT(J+1)*WT(J))  
        WTP_M(J) = WT_M(J)/protonM
        RHO_M(J) = SQRT(RHO(J)*RHO(J+1))
        DEN_M(J) = RHO_M(J)/WT_M(J)
       enddo
CZ
      DO J=1,NZ
      if (IZPRNTFLG5.eq.1) 
     | PRINT 5998, J, RHO(J), RHO_INI(J),WT(J),GZ(J)
C     | , USOL(LOI,J)*RHO(J)/WGT(LOI), UNUM(LOI,J), UNUM_INI(LOI,J)
C     | , USOL(LO,J)*RHO(J)/WGT(LO), UNUM(LO,J), UNUM_INI(LO,J)
     | , U_ADV(J)
      ENDDO
 5998 FORMAT(1x,"After SUBSONIC:",I5,1P14E18.8)
C
      LINERT = 0
C      
 2009 FORMAT(3x,"J", 6x,"RHO",7x,"U",8x,"T",8x,"DEN",8x,"WT"
     | ,8x,"C(H)",5x,"C(N2)"
     | ,5x,"n(H)",5x,"n(N2)")
      WRITE(101,2009)

CZ  this includes both boundary and midpoint variables
      DO J=1,NZ
      write(101,3000), J, RHO(J), U_ADV(J), T(J), DEN_M(J), WTP(J)
     | , USOL(LH,J), USOL(LN2,J), UNUM(LH,J), UNUM(LN2,J)
     | , rho_hydro(J)

C
      ENDDO
      WRITE(101,2009)
      
      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ       
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF
      
C 
C      PRINT 5008, J, R(J)-R_PLT, DEN(J), UNe(J)
C     | , (UNUM_OLD(I,J),I=NQT,NQT+10)
 5008 FORMAT (I5,1P20E10.2)       
C
C      PRINT 5104,USOL(1,1)*DEN(1)/WGT(1),DEN(1),RHO(1),UNe(NZ),RHO(NZ)
5104  FORMAT (1x,1P5e11.3)
C  
C   Calculate diffusion-limited flux, assuming a fixed total number
C   density at the lower boundary
      fh2 = UNUM(LH,1)/DEN_M(1)
      phidif = 2.5e13 * fh2
      Tinf=T(NZ)  !  changed back to NZ 
C    
C ***** PRINT OUT INITIAL DATA *****
cz      PRINT*,"1278 EMAX =",EMAX
      
      CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
      
cz      PRINT*,"1282 EMAX =",EMAX      
CZ       
      IF (NAN_DBG.EQ.1) THEN
         PROD = UNUM(17,1)*UNUM(18,1)*UNUM(19,1)
         PRINT*,"1306 PROD =",PROD
      ENDIF
CZ initialize i_exo
      i_exo = NZ
      
CZ-0100
CSSS ***** START THE TIME-STEPPING LOOP *****
      DO 1 N=1,NSTEPS
C
      NN_STEPS = N

            IF (IZDEBUG.EQ.1) PRINT*, "BEGINNING OF TIME-STEPPING LOOP"
      IF(IZDEBUGEXO.EQ.1) PRINT*,"Start Time loop, i_exo =", i_exo	    
CZ       
      IF (NAN_DBG.EQ.1) THEN
         PROD = UNUM(17,1)*UNUM(18,1)*UNUM(19,1)
         PRINT*,"1422 PROD =",PROD
      ENDIF
CZ  
      NN = NN + 1
      IS = 36
      JS = NZ-1  ! hydro57 Was NZ
      EMAX = 0
      MS = (NN-1)/MP
      SM = (NN-1)/PM
      IF(NN.EQ.NSTEPS) SM = MS

CZ Save current values of all species to compare against new
      IF (chemkluge .EQ. 1) THEN
      DO 5275 J=1,NZ-1 ! hydro57 was NZ
      DO 5274 I=1,NSP
      USAVE1(I,J)=USOL(I,J)
      USAVEUN1(I,J)=UNUM(I,J)
 5274  CONTINUE
 5275  CONTINUE
      ENDIF

CZ Section added for planet conversion
      IF (CONVEM .GT. 0) THEN
        IF  (G.GT.373)   THEN
           G = 0.999 * G
           PRINT*, "Reducing G, now ",G
        ENDIF
        IF  (G.LT.373)   THEN
           G = 373
           PRINT*, "Setting G to",G
        ENDIF	

        IF  (R_PLT.GT.3.396E8)   THEN
           R_PLT = 0.999 * R_PLT
           PRINT*, "Reducing R_PLT, now ",R_PLT
           IF (IZLOGGRID.ge.1) THEN
            DO 5234 J=1,NZ
             zl(J) = 2. + (J-1)/ZGRIDFAC 
             zk(J) = 10.**zl(J)
             R(J) = 1.e5*(zk(J)) + R_PLT
 5234    CONTINUE
	     PRINT*,"Recalculating Altitudes of Grid Points"
	     ENDIF
        ENDIF
	 
        IF  (R_PLT.LT.3.396E8)   THEN
           R_PLT = 3.396E8
           PRINT*, "Setting R_PLT to",R_PLT
	
           IF (IZLOGGRID.ge.1) THEN
            DO 5236 J=1,NZ
             zl(J) = 2. + (J-1)/ZGRIDFAC 
             zk(J) = 10.**zl(J)
             R(J) = 1.e5*(zk(J)) + R_PLT
 5236    CONTINUE
 	     PRINT*,"Recalculating Altitudes of Grid Points"
           ENDIF
	
	     if (IZPRNTGRID.EQ.1) THEN
              PRINT 6736,R(J)/1E5, (R(J)-R_PLT)/1E5,J
 6736     FORMAT("R(J) = ",E10.3," km,  Alt = ", 
     |        E10.3," km,  J = ",I4)     
            ENDIF
	    
	    
         ENDIF

        IF  (bigM.GT.6.419E26 .and. (NZ .LE. 272))   THEN
           PRINT*,"NZ = ",NZ
           bigM = 0.999 * bigM
           PRINT*, "Reducing bigM, now ",bigM
        ENDIF
        IF  (bigM.LT.6.419E26)    THEN
           bigM = 6.419E26
           PRINT*, "Setting bigM to",bigM
        ENDIF
	
        IF  (FSCALE.GT.0.43)   THEN
           FSCALE = 0.999 * FSCALE
           PRINT*, "Reducing FSCALE, now ",FSCALE
        ENDIF
        IF  (FSCALE.LT.0.43)   THEN
           FSCALE = 0.43
           PRINT*, "Setting FSCALE to",FSCALE
        ENDIF
	
        IF  (ALB.GT. 0.215)   THEN
           ALB = 0.999 * ALB
           PRINT*, "Decreasing ALB, now ",ALB
        ENDIF
        IF  (ALB.LT. 0.215)   THEN
           ALB = 0.215
           PRINT*, "Setting ALB to",0.215
        ENDIF
	
        IF  (ZTROP .LT. 1.5E6)   THEN
           ZTROP = 1.01 * ZTROP
           PRINT*, "Increasing ZTROP, now ",ZTROP
        ENDIF
        IF  (ZTROP .GT. 1.5E6)   THEN
           ZTROP = 1.5E6
           PRINT*, "Setting ZTROP to ",ZTROP
        ENDIF
	
         do 5231 J=1,NZ ! recalc if planet conversion in process
  	     GZ(J) = G * (R_PLT/R(J))**2
	     HSCALE(J) = bol_k/(WT(J)*GZ(J))*T(J)
cz	     RHO(J) = RHO(J-1)* exp(-1.E5*(zk(j) - zk(j-1))/HSCALE(J)) * 
cz     |         T(j-1)/T(j) * WT(J)/WT(J-1)
 5231  CONTINUE

C                1 -- CONVERT UNUM INTO USOL
       CALL NORM_MMR(1)
      ENDIF

        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0200, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-0200
CZ This section added by MEZ to control and stabilize chemical species MMR levels to those chosen by user
      NZ1 = NZ-1

CZ CONTROL N2

      IF  (USOL(LN2,1).LT.N2MMR)   THEN
         print*,"USOL(LN2,1)=",USOL(LN2,1)
         PRINT*, " INCREASING N2 BY 1%"
         DO J=1,NZ1 ! hydro56 was NZ
            USOL(LN2,J)=USOL(LN2,J)/(0.99)
            UNUM(LN2,J)=UNUM(LN2,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LN2,1).GT.(N2MMR*1.015))  THEN
         PRINT*, " DECREASING N2 BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LN2,J)=USOL(LN2,J)*(0.99)
            UNUM(LN2,J)=UNUM(LN2,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	 

CZ CONTROL CO2

      IF  (USOL(LCO2,1).LT.CO2MMR)   THEN
         PRINT*, " INCREASING CO2 BY 1%"
         DO J=1,NZ1 ! hydro56 was NZ
            USOL(LCO2,J)=USOL(LCO2,J)/(0.99)
            UNUM(LCO2,J)=UNUM(LCO2,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LCO2,1).GT.(CO2MMR*1.015))  THEN
         PRINT*, " DECREASING CO2 BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LCO2,J)=USOL(LCO2,J)*(0.99)
            UNUM(LCO2,J)=UNUM(LCO2,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	      

CZ CONTROL CO
      IF  (USOL(LCO,1).LT.COMMR)   THEN
         PRINT*, " INCREASING CO BY 1%"
         DO J=1,NZ1 ! hydro56 was NZ
            USOL(LCO,J)=USOL(LCO,J)/(0.99)
            UNUM(LCO,J)=UNUM(LCO,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LCO,1).GT.(COMMR*1.015))  THEN
         PRINT*, " DECREASING CO BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LCO,J)=USOL(LCO,J)*(0.99)
            UNUM(LCO,J)=UNUM(LCO,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF

CZ CONTROL O2
      IF  (USOL(LO2,1).LT.O2MMR)   THEN
         PRINT*, " INCREASING O2 BY 1%"
         DO J=1,NZ1 ! hydro56 was NZ
            USOL(LO2,J)=USOL(LO2,J)/(0.99)
            UNUM(LO2,J)=UNUM(LO2,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LO2,1).GT.(O2MMR*1.015))  THEN
         PRINT*, " DECREASING O2 BY 1%"
         DO J=1,NZ1 ! hydro56 was NZ
            USOL(LO2,J)=USOL(LO2,J)*(0.99)
            UNUM(LO2,J)=UNUM(LO2,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	 

CZ CONTROL O 

      IF  (USOL(LO,1).LT.OMMR)   THEN
         PRINT*, " INCREASING O BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LO,J)=USOL(LO,J)/(0.99)
            UNUM(LO,J)=UNUM(LO,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LO,1).GT.(OMMR*1.015))  THEN
         PRINT*, " DECREASING O BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ 
            USOL(LO,J)=USOL(LO,J)*(0.99)
            UNUM(LO,J)=UNUM(LO,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
 

CZ CONTROL H2      
      IF ( (USOL(LH2,1).LT.H2MMR)) THEN !.AND. (IZH2X1P01.EQ.9) ) THEN
         PRINT*, " INCREASING H2 BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ 
            USOL(LH2,J)=USOL(LH2,J)/(0.99)
            UNUM(LH2,J)=UNUM(LH2,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
cz      PRINT*,"IZH2X1P01 = ",IZH2X1P01
      IF ( (USOL(LH2,1).GT.H2MMR*1.015)) THEN !.AND. (IZH2X1P01.EQ.10) ) THEN
         PRINT*,"USOL(LH2,1) = ",USOL(LH2,1)      
         PRINT*, " DECREASING H2 BY 1%"
         DO J=1,NZ1 ! hydro56 was NZ
            USOL(LH2,J)=USOL(LH2,J)*(0.99)
            UNUM(LH2,J)=UNUM(LH2,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	

CZ CONTROL H

      IF  (USOL(LH,1).LT.HMMR)   THEN
         PRINT*, " INCREASING H BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LH,J)=USOL(LH,J)/(0.99)
            UNUM(LH,J)=UNUM(LH,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LH,1).GT.(HMMR*1.015))  THEN
         PRINT*, " DECREASING H BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ  
            USOL(LH,J)=USOL(LH,J)*(0.99)
            UNUM(LH,J)=UNUM(LH,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	


CZ CONTROL HE
      IF  (USOL(LHE,1).LT.HEMMR)   THEN
         PRINT*, " INCREASING HE BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LHE,J)=USOL(LHE,J)/(0.99)
            UNUM(LHE,J)=UNUM(LHE,J)/(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	
      
      IF  (USOL(LHE,1).GT.(HEMMR*1.015))  THEN
         PRINT*, " DECREASING HE BY 1%"
         DO J=1,NZ1  ! hydro56 was NZ
            USOL(LHE,J)=USOL(LHE,J)*(0.99)
            UNUM(LHE,J)=UNUM(LHE,J)*(0.99)
         ENDDO ! END OF LOOP 
      ENDIF	 

      IF (IZMMRDEBUG.EQ.1) THEN
         PRINT 7359,USOL(LH2,1),UNUM(LH2,1)/DEN_M(1),DEN_M(1)
 7359 FORMAT(1X,"At Grid Pt. 1, MMR of H2 = ",F7.5,2x,
     | "VMR of H2 = ",F6.4,2X,"DEN = ",1P1E10.3," /cm3")
          PRINT 7358,USOL(LCO2,1),UNUM(LCO2,1)/DEN_M(1),RHO(1)
 7358 FORMAT(1X,"At Grid Pt. 1, MMR of CO2 = ",F7.5,2x,
     | "VMR of CO2 = ",F6.4,2X,"RHO = ",1P1E10.3," g/cm3")
      ENDIF

      IF ((IZMMRDEBUG.EQ.2).and.(N.EQ.NSTEPS)) THEN
         PRINT 7359,USOL(LH2,1),UNUM(LH2,1)/DEN_M(1),DEN_M(1)
         PRINT 7358,USOL(LCO2,1),UNUM(LCO2,1)/DEN_M(1),RHO(1)
      ENDIF

        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0300, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-0300
C-TF chemical reaction rate coefficients are set up in Rates.f      
C
      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF

            IF (IZDEBUG.EQ.1) PRINT*, "CALLING RATES"
      CALL RATES
C      
C-TF set up the components in the Jacobian matrix in DIFCO.f
C
            IF (IZDEBUG.EQ.1) PRINT*, "CALLING DIFCO"
      CALL DIFCO(NN,NSTEPS)
C
C     ZY = SOLAR ZENITH ANGLE (IN DEGREES)
C     LTIMES = COUNTER FOR PHOTORATE SUBROUTINE
C     ISEASON = TELLS WHETHER P AND T VARY WITH TIME (THEY DON'T FOR
C               ISEASON < 3)
C     IZYO2 = TELLS WHETHER SOLAR ZENITH ANGLE VARIES WITH TIME (0 SAYS
C             IT DOESN'T; 1 SAYS IT DOES)
C     IO2 = 0 FOR ALLEN AND FREDERICK O2 SCHUMANN-RUNGE COEFFICIENTS
C         = 1 FOR EXPONENTIAL SUM FITS (FOR LOW-O2 ATMOSPHERES)
C     INO = 0 FOR ALLEN AND FREDERICK NO PREDISSOCIATION COEFFICIENTS
C         = 1 FOR MODIFIED CIESLIK AND NICOLET FORMULATION
      IDO = 0
C      IF (NN.EQ.NSTEPS) IDO = 1
            IF (IZDEBUG.EQ.1) PRINT*, "CALLING PHOTO"
      CALL PHOTO(ZY,AGL,LTIMES,ISEASON,IZYO2,IO2,INO,IDO)

            IF (IZDEBUG.EQ.1) PRINT*, "RETURNED FROM PHOTO"
	    
      IF (EXO_DBG.EQ.1) THEN
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF	    
C
C   TIME-DEPENDENT BOUNDARY CONDITIONS
C-TF  Escape of hydrogen: VEFF(H) = (Bi/Ha)/Nt, Nt=total number density
C-TF  These are diffusion-limited escape flux at the lower boundary  
C
      IF(NN.EQ.NSTEPS) THEN 

CZ MEZ changed code to normalize at homopause, not R(NZ)
cz      J=5  ! ~100 km altitude
      J=1   ! 100 km altitude w/ JK log grid
      BOVERH = DI(LH,J)*DEN(J)/HSCALE(J)
      BOVERH2 = DI(LH2,J)*DEN(J)/HSCALE(J)
      BOVERH2O = DI(LH2O,J)*DEN(J)/HSCALE(J)
CZ      
      fvH = UNUM(LH,J)/DEN_M(J)
      fvH2 = UNUM(LH2,J)/DEN_M(J)      
      fvH2O = UNUM(LH2O,J)/DEN_M(J)
      fvtot = fvH + 2*fvH2 + 2*fvH2O
      fvtot1 = fvtot + 1.
CZ
      FH_dif = (BOVERH*fvH/fvtot1)*(R(J)/R(1))**2
      FH2_dif = 2*(BOVERH2*fvH2/fvtot1)*(R(J)/R(1))**2
      FH2O_dif = 2*(BOVERH2O*fvH2O/fvtot1)*(R(J)/R(1))**2
CZ      
      dif_limited_flux = FH_dif + FH2_dif + FH2O_dif
C      
      PRINT 603, J, BOVERH, BOVERH2, BOVERH2O
     | ,FH_dif, FH2_dif, FH2O_dif
     | ,fvH
     | ,fvtot
     | ,dif_limited_flux
 603  FORMAT(/1x,"@J=",I5,"  BOVERH =",1PE10.3,"  BOVERH2 =",1PE10.3
     1 ,"  BOVERH2O =",1PE10.3,"  PHI(LH) = ",1PE10.3
     2 ,"  PHI(LH2) = ",1PE10.3, "  PHI(LH2O) = ",1PE10.3/
     3 ,"volume fraction of H at 100 km = ",1PE10.3/
     4 ,"total vol. fraction of H at 100 km = ",1PE10.3/     
     5 ,"diffusion limited flux at 100 km = ", 1PE10.3,/)


 7603  FORMAT(/1x,"@J=",I5,"  BOVERH =",1PE10.3,"  BOVERH2 =",1PE10.3
     1 ,"  BOVERH2O =",1PE10.3,"  PHI(LH) = ",1PE10.3
     2 ,"  PHI(LH2) = ",1PE10.3, "  PHI(LH2O) = ",1PE10.3/
     3 ,"volume fraction of H at 110 km = ",1PE10.3/
     4 ,"total vol. fraction of H at 110 km = ",1PE10.3/       
     5 ,"diffusion limited flux at 110 km = ", 1PE10.3,/)
      
C
cz       CrossoverM = WGTM(LH) +
cz     | bol_k*T(NZ)*VEFF(LH)/DI(LH,NZ)/GZ(NZ)/(UNUM(LH,NZ)/DEN(NZ))/
cz     | protonM
C
cz       PRINT*, "Crossover mass =", CrossoverM
cz     | , bol_k*T(NZ)*VEFF(LH)/DI(LH,NZ)/GZ(NZ)

CZ Replaced above crossover mass calc with calc at homopause per JFK
       CrossoverM = WGTM(LH) +
     | bol_k*T(1)*VEFF(LH)/DI(LH,1)/GZ(1)/(UNUM(LH,1)/DEN_M(1))/
     | protonM

       PRINT*, "Crossover mass =", CrossoverM
     | , bol_k*T(1)*VEFF(LH)/DI(LH,1)/GZ(1)

        IF(CROSDBG .EQ. 1) THEN
          PRINT*,"T(NZ) =",T(NZ)
          PRINT*,"VEFF(LH) =",VEFF(LH)
          PRINT*,"DI(LH,NZ) =",DI(LH,NZ)
          PRINT*,"GZ(NZ) =",GZ(NZ)
        ENDIF
C     
      ENDIF

      IF (NN.LT.NSTEPS) GO TO 18
C      PRINT 97
  97  FORMAT(//1X,"PHOTOLYSIS RATES")
C      PRINT 98
  98  FORMAT(/5X,"Z",7X,"PO2",6X,"PO2D",5X,"PCO2",5X,"PCO2D",4X,
     2  "PH2O",5X,"PO3",6X,"PO3D",7X,"PNO")
C      PRINT 99,(Z(I),PO2(I),PO2D(I),PCO2(I),PCO2D(I),PH2O(I),
C     2  PO3(I),PO3D(I),PNO(I),
C     3  I=1,NZ,3)
  99  FORMAT(1X,1P9E9.2)
C  
  18  CONTINUE

        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0400, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-0400    NEWTON SPECIES
      INEWT = 20
      
      IF(IZDEBUGEXO.EQ.1) PRINT*,"Setting up Jacobian, i_exo =", i_exo
 
      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF

C ***** SET UP THE JACOBIAN MATRIX AND RIGHT-HAND SIDE *****
C
      DO 17 J=1,NAQ
      DO 17 K=1,NAQ
  17  DJACN(J,K) = 0.
      DO 19 K=1,NAQ
  19  RHSN(K) = 0.
C
C   COMPUTE CHEMISTRY TERMS AT ALL GRID POINTS
C-TF FVAL=Pi/N-Li*fi
      IDO = 0
      IF (NN.EQ.NSTEPS) IDO = 1
C
C-TF  SAVE SPECIES FOR RECALCULATION OF THE STEPS
      DO 275 J=1,NZ1 ! hydro56 was NZ
      DO 274 I=1,NSP
      USAVE(I,J)=USOL(I,J)
CZ Added     
      USAVEUN(I,J)=UNUM(I,J)
CZ      
 274  CONTINUE
 275  CONTINUE
C
C-TF  GO TO 1616 TO JUMP OVER THE NEWTON SPECIES
C      GO TO 1616
C-TF  NH can be adjusted
C

      NH = min(NZ1,NZchemMXN) ! hydro56 was NZ

cz      print*,"NZ1=",NZ1
cz      PRINT*,"NH=",NH
 
CZ      NH = NZ !- 30
      NQ1 = NQ + 1
C
C   START NEWTON ITERATION
      DO 1010 J=1,NH   ! ALTITUDE LEVELS NEWTON METHOD IS TRIED
C      IF (J .EQ. 2) STOP
C
      DO 1515 I=1,NSP4
      UNUM1(I)=UNUM(I,J)
C
      IF (NEWTON_DEBUG .EQ. 1) THEN
      PRINT 9801, NN, J, I, UNUM1(I), ISPEC(I)
9801  FORMAT ("after assigning UNUM1:",1x,3I3,1PE10.2,5x,A8)
C      PRINT 9802, N, 1, UNUM(NQ1,1), UNUM(NQT,1)
9802  FORMAT ("UNUM at level 1:",1x,2I3,1P2E10.2)
      ENDIF
C      
1515  CONTINUE      
C
C
C-TF  USE CHEMPLONE TO PREDICT CHEMICAL EQUILIBRIUM NUMBER DENSITIES     
C-TF  FOR INITIALIZATION OF THE NEWTON SPECIES
      NQ3=NQ1+2
      IF (NEWTON_DEBUG .EQ. 1) THEN
      PRINT*,  ''
      PRINT*, "@ level :",J
      PRINT 9011, (UNUM1(K),K=NQ1,NQ3)
9011  FORMAT (1x,"initial unum:", 1P3E10.2)
      ENDIF
      
      PROD = UNUM1(NQ1)*UNUM1(NQ1+1)*UNUM1(NQ3)
      IF(PROD.EQ.PROD*0.5+1 .AND. PROD.NE.2.)THEN
         PRINT*, "Error 2 - NEWTON SPECIES NaN - coupledSUA terminating"
         PRINT*, "N = ",N
         WRITE(6942,6943) 2.
           IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,TIME,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WT,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
          ENDIF
         GOTO 2911
      ENDIF
      
      DO KK = NAQ,1,-1
      K = NQ + KK



C-TF  USE THE CHEMISTRY TO PREDICT THE STEADY STATE NUMBER DENSITY
C-TF  ONLY AT THE FIRST STEP
cz      IF (J .EQ. 1) THEN

CZ Per J. Kasting, change to use CHEMPLONE as the first guess for Newton for all altitudes J
        IF(chemkluge .eq. 1 ) THEN
	      if(UNUM1(K) .LT. 1) UNUM1(K) = 1
	      UNUM1_SAVE(15) = UNUM1(15)
	      UNUM1_SAVE(16) = UNUM1(16)
	      UNUM1_SAVE(17) = UNUM1(17)
        endif	      
	      
      CALL CHEMPLONE(UNUM1,J,K,XP,XL)
            IF (chemkluge .eq. 1 .and. XL .LT. SMALL) then
                 PRINT*, "WARNING, XL = ",XL
                 XL = SMALL
           endif
      UNUM1(K) = XP/XL
            PROD = UNUM1(15)*UNUM1(16)*UNUM1(17)

      IF (PROD.EQ.PROD*0.5+1 .AND. PROD.NE.2.) THEN 
         PRINT 7103,UNUM1(15),UNUM1(16),UNUM1(17),XP,XL,J    
 7103 FORMAT("ERROR: NEWTON SPEC. NaN, LINE 1935,  UNUM = ",1P3E10.2,
     |     "  XP = ",1P1E10.2,"  XL = ",1P1E10.2," J =",I5)      
         PRINT*, " "
cz         PRINT*,"reaction rates A=",A
          IF (chemkluge .eq. 1) THEN
	     UNUM1(15) = UNUM1_SAVE(15)
	     UNUM1(16) = UNUM1_SAVE(16)
	     UNUM1(17) = UNUM1_SAVE(17)	     
             goto 2051
          ENDIF


        PRINT*,"Mass Mixing Fractions" 
      PRINT*," of Long Lived Neutral Species:"
       PRINT 6921
      
      NZZ1 = MIN(NZ1,MGRIDMX)  ! hydro56 was NZ
      DO 7119 J1=1,NZZ1,MSKIP
         PRINT 6920,J1,(USOL(I,J1),I=1,11),(R(J1)-R_PLT)/1E5
 7119 CONTINUE
       PRINT 6921
CZ
      PRINT*, " "
      PRINT*,"Particle Number Densities per cubic centimeter " 
      PRINT*,"of Ions (O+,N+,H+) and Newton Method Species (O3,HO2,OH):"
       PRINT*,"Grid   O+          N+         H+        O3         HO2
     |       OH        Alt (km)"
      
      NZZ1 = MIN(NZ1,MGRIDMX) ! hydro56 was NZ
      DO 7116 J2=1,NZZ1,MSKIP
         PRINT 6918,J,(UNUM(I,J2),I=12,17),(R(J2)-R_PLT)/1E5  !hydro57 was 6.371E8
 7116 CONTINUE
 
       PRINT*,"Grid   O+          N+         H+        O3         HO2
     |       OH        Alt (km)"

        PRINT*,"Error 2 - NaN calculated, coupledSUA terminating"
        PRINT*,"N = ",N
        WRITE(6942,6943) 2.
	     IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,TIME,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
            ENDIF
        GOTO 2911
       
       ENDIF

 2051  CONTINUE

cz      ENDIF

cz      IF (J .NE. 1) THEN
cz      UNUM1(K) = UNUM(K,J-1)
cz      ENDIF

      ENDDO  ! KK loop

      IF (NEWTON_DEBUG .EQ. 1) THEN
      PRINT*,  ''
      PRINT*, "@ level :",J
      PRINT 9012, (UNUM1(K),K=NQ1,NQ3)
9012  FORMAT ("predicted unum:", 1P3E10.2)
      IF (PROD.EQ.PROD*0.5+1 .AND. PROD.NE.2.) THEN
        PRINT*,"**NEWTON SPECIES NaN**"
      ENDIF	
      ENDIF

      DO 1011 INN=1,INEWT   ! TRY THE NEWTON METHOD FOR INEWT TIMES
      
C-TF  COMPUTE THE CHEMICAL PRODUCTION/LOSS TERM FOR SPECIES
C-TF  AT ALTITUDE LEVEL J

CZ      PRINT*, "CALLING CHEMONE"
      
      CALL CHEMONE(UNUM1,NAQ,J,F,IFIRST)
C
CZ Check for NaN in Newton species
      PROD = UNUM1(NQ1)*UNUM1(NQ1+1)*UNUM1(NQ3)
      IF (PROD.EQ.PROD*0.5+1 .AND. PROD.NE.2.) THEN
        PRINT 7104,UNUM1(15),UNUM1(16),UNUM1(17)
      ENDIF	
 7104 FORMAT("**NEWTON SPEC. NaN, LINE 2062**,  UNUM = ",1P3E10.2)

CZ Test ratios between species to see if Newton method needed
      IF ((UNUM1(17)/UNUM1(15).GT.NEWTMAX)) THEN
        IF(NEWT_DBG.EQ.1) PRINT 7102,J,",  OH/O3 >",NEWTMAX,",  UNUM=",
     |     UNUM1(15),UNUM1(16),UNUM1(17)	
        GO TO 1010
      ENDIF

      IF ((UNUM1(17)/UNUM1(16).GT.NEWTMAX)) THEN
        IF(NEWT_DBG.EQ.1)PRINT 7102,J,",  OH/HO2 > ",NEWTMAX,",  UNUM=",
     |     UNUM1(15),UNUM1(16),UNUM1(17)	
        GO TO 1010
      ENDIF

      IF ((UNUM1(16)/UNUM1(17).GT.NEWTMAX)) THEN
        IF(NEWT_DBG.EQ.1)PRINT 7102,J,",  HO2/OH > ",NEWTMAX,",  UNUM=",
     |     UNUM1(15),UNUM1(16),UNUM1(17)	
        GO TO 1010
      ENDIF
      
      IF ((UNUM1(16)/UNUM1(15).GT.NEWTMAX)) THEN
        IF(NEWT_DBG.EQ.1)PRINT 7102,J,",  OH/O3 > ",NEWTMAX,",  UNUM=",
     |     UNUM1(15),UNUM1(16),UNUM1(17)	
        GO TO 1010
      ENDIF      

      IF ((UNUM1(15)/UNUM1(17).GT.NEWTMAX)) THEN
        IF(NEWT_DBG.EQ.1)PRINT 7102,J,",  O3/OH > ",NEWTMAX,",  UNUM=",
     |     UNUM1(15),UNUM1(16),UNUM1(17)	
        GO TO 1010
      ENDIF

      IF ((UNUM1(15)/UNUM1(16).GT.NEWTMAX)) THEN
        IF(NEWT_DBG.EQ.1)PRINT 7102,J,",  O3/HO2 > ",NEWTMAX,",  UNUM=",
     |     UNUM1(15),UNUM1(16),UNUM1(17)	
        GO TO 1010
      ENDIF
cz      GO TO 1010  ! Use to skip Newton
 7102 FORMAT("Newton Skipped, J =",I3,A,1P1E10.1,A,1P3E11.2)

CZ
      IF (NEWTON_DEBUG .EQ. 1) THEN
      PRINT*,  ''
      PRINT*, "@ level :",J
      PRINT 9897, (UNUM1(K),K=NQ1,NQ3)
      PRINT 9897, F
 9897 FORMAT(1x,"1st time calling chemone",1P3E10.2)
      PROD = UNUM1(NQ1)*UNUM1(NQ1+1)*UNUM1(NQ3)
      IF (PROD.EQ.PROD*0.5+1 .AND. PROD .NE. 2) THEN 
        PRINT*,"**NEWTON SPECIES NaN**" 
      ENDIF	
      ENDIF
C      
C-TF  COMPUTE THE CHEMICAL PRODUCTION/LOSS TERMS FOR
C-TF  PERTURBED SPECIES K AT ALTITUDE LEVEL J

CZ      PRINT*, "STARTING LOOP TO COMPUTE CHEM PROD/LOSS TERMS"
CZ In order to find dF/dx, where x = UNUM, need to perturb unum by adding dx = epsilon x UNUM
CZ Then call CHEMONE to calculate new F = FP.  Then fill the Jacobian matrix with dF/dx values
CZ where dF = F - FP

      DO 1313 K=1,NAQ
        KK = K + NQ
        XS = UNUM1(KK)
        DX = EPSJ*UNUM1(KK)
        UNUM1(KK)=UNUM1(KK)+DX
        CALL CHEMONE(UNUM1,NAQ,J,FP,IFIRST)
C
C-TF  FILL THE MATRIX WITH dFi/dCi
        DO 1414 L=1,NAQ
        DJACN(L,K) = (FP(L) - F(L))/DX

C      PRINT 9898, K, L, DJACN(L,K), DX, EPSJ, UNUM1(KK)
9898  FORMAT (1x,"Filling Matrix for SGEFA:",2I3,1P13E10.2)

1414    CONTINUE
        UNUM1(KK)=XS
1313  CONTINUE
C
C-TF comment out the following print statement
      IF (NEWTON_DEBUG .EQ. 1) THEN
      PRINT 9997, INN, J
9997  FORMAT (/1x,"NEWTON STEP:",I3,2x,"AT LEVEL:", I3)
      PRINT 9998, (UNUM1(K),K=NQ1,NQT)
9998  FORMAT ("UNUM INITIAL", 1P5E10.2)
      PRINT 9993,DX
9993  FORMAT ("DX follows...", 1P5E10.2)

      PRINT 9994,FP
9994  FORMAT ("FP follows...", 1P5E10.2)
      PRINT 9996,F
9996  FORMAT ("Fo follows...", 1P10E10.2)
      PRINT 9995,DJACN
9995  FORMAT ("DJACN follows...", 1P10E10.2)
      endif

      PROD = UNUM1(NQ1)*UNUM1(NQ1+1)*UNUM1(NQ3)
      IF (PROD.EQ.PROD*0.5+1) THEN
        PRINT*,"Error 2  - NEWTON SPECIES NaN, coupledSUA terminated"
        PRINT*,"N = ",N
        WRITE(6942,6943) 2.
	   IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)  
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
          ENDIF
          GOTO 2911
      ENDIF
9999  FORMAT (1P7E12.3)

CZ Factor Jacobian DJACN
      CALL SGEFA(DJACN,NAQ,NAQ,IPVTN,INFO)
      IF (INFO.EQ.0) GO TO 1020
      PRINT 1021, INFO,INN,NRAIN,J
C-TF the following format can only handle NRAIN<1000. be careful
1021  FORMAT(//1X,"NEWTON SOLVER FAILED SGEFA",5X,"INFO =",I3,
     2  "GRID STEP =",I3,2X,"NRAIN =",I3,2x,"AT LEVEL ",I3)
      PRINT *,"IPVTN follows..."
      PRINT 9999, DJACN
C
      NTF_STOP_FLG = 1
C  
CZ Solve Jabobian DJACN
1020  CALL SGESL(DJACN,NAQ,NAQ,IPVTN,F,0)

Cz Check to see how closely converged the solution is; if F/UNUM1 > 1E-6
      LTEST = 0
      DO 2222 K=1,NAQ
      KK = K + NQ
      UNUM1(KK)=UNUM1(KK)-F(K)
C
      IF (ABS(F(K)/UNUM1(KK)) .GT. 1.E-6) LTEST = 1
2222  CONTINUE

      IF (NEWTON_DEBUG .EQ. 1) THEN 
      PRINT 9991,-F
9991  FORMAT ("DX follows...", 1P10E10.2)
      PRINT 9992,INN,UNUM1(NQ1),UNUM1(NQ+2),UNUM1(NQT)
9992  FORMAT ("UNUM UPDATED follows...", I3,1x,1P5E10.2)
      ENDIF

      IF (LTEST.EQ.0) GO TO 1516
1011  CONTINUE
C
      PRINT 2500,I,NRAIN
C-TF the following format can only handle NRAIN<1000. be careful
2500  FORMAT(//1X,"NEWTON SOLVER FAILED TO CONVERGE "/,5X,
     2  "GRID STEP =",I3,2X,"NRAIN =",I3)
C
      NTF_STOP_FLG = 1
C      STOP
C
1516  CONTINUE

cz      PRINT*, "06 R(271) = ",R(271) 

      IF (NEWTON_DEBUG .EQ. 1) THEN
      PRINT*, "NEW VALUES TAKEN FOR LEVEL", J
      PRINT 9990
9990  FORMAT (/1x)
      ENDIF
C
CZ NQ=14             ! # of long-lived species
CZ NQ1 = NQ+1 
CZ NAQ=3             ! # of Newtonian species
CZ  NQT=NQ+NAQ        ! # of long-lived species + Newtonian species

CZ So NQ1,NQT cycles over the Newtonian species:

      DO 2223 K=NQ1,NQT
      
      UNUM(K,J)=UNUM1(K)
      IF (chemkluge .eq. 1 .AND. UNUM(K,J) .LE. 1.) THEN
         UNUM(K,J) = 1.
      ENDIF
      
      USOL(K,J)=UNUM(K,J)*WGT(K)/RHO(J)
CZ 
      IF (chemkluge .eq. 1) THEN

      IF (USOL(LH,J) .LT. SMALL) THEN
      print*,"WARNING: @2717 USOL(H) = ",USOL(LH,J)
      PRINT*,"J = ",J
      USOL(LH,J) = SMALL
      ENDIF
      IF (USOL(LH2,J) .LT. SMALL) THEN
      print*,"WARNING: @2722 USOL(H2) = ",USOL(LH2,J)
      USOL(LH2,J) = SMALL
      PRINT*,"J = ",J      
      ENDIF      
CZ
      IF (USAVE(K,J) .LT. SMALL) then
          USAVE(K,J) = SMALL
      ENDIF
      
      delUSOL = (USOL(K,J)-USAVE(K,J))/USAVE(K,J)
      EREL = ABS(delUSOL)
      
      delUNUM = (UNUM(K,J)-USAVEUN(K,J))/USAVEUN(K,J)
      
      IF (delUNUM .GT. 0.5) THEN
         UNUM(K,J) = 1.5*USAVEUN(K,J)
      ENDIF

      IF (delUNUM .LT. -0.5) THEN
         UNUM(K,J) = 0.5*USAVEUN(K,J)
      ENDIF
CZ      

CZ Added to limit UNUM > 0
CZ3 changed from 1 to 1.E-5
      if (UNUM(K,J) .LT. 1) UNUM(K,J) = SMALL
      USOL(K,J)=UNUM(K,J)*WGT(K)/RHO(J)
CZ  
      ENDIF ! chemkluge = 1

      IF (USOL(LH,J) .LT. SMALL) THEN
      print*,"WARNING: @2755 USOL(H) = ",USOL(LH,J)
      PRINT*,"J = ",J
      USOL(LH,J) = SMALL
      ENDIF
      IF (USOL(LH2,J) .LT. SMALL) THEN
      print*,"WARNING: @2760 USOL(H2) = ",USOL(LH2,J)
      USOL(LH2,J) = SMALL
      PRINT*,"J = ",J      
      ENDIF         
      
CZ Added to look for two altitudes set to UNUM = 1 for same species
      IF(UNUM(K,J) .EQ. 1 .AND. UNUM(K,J-1) .EQ. 1
     | .AND. (N .GT. 10) .AND. J .GE. 2 
     | .AND. K .LE. 14)     THEN
        PRINT*,"Error 9, UNUM(K,J) = UNUM(K,J-1) =1, J =" 
        PRINT*, "Alt. Grid J = ",J
	PRINT*, "Species K = ", K
        WRITE(6942,6943) 9.
            IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)  
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)     
            ENDIF
cz        WRITE(6942,6949) USOL(LH2,J)
        GOTO 2911      
      ENDIF
CZ

CZ Limit emax to levels below NZchemMXE  
CZ added USOL (MMR) at least 1 part in 10^4

CZ For rev cs_r28.f, replaced NZchemMXE with min UNUM(), USOL()
      IF ((EREL .GT. EMAX) .AND. (UNUM(K,J) .GE. 1.E10) 
cz     |    .AND. (J.LE.NZchemMXE)
     |   .AND. (USOL(K,J) .GE. 1.E-4))   THEN
        UMAX = USOL(K,J)
        IS = K
        JS = J
        EMAX = EREL
cz        PRINT*,"EMAX =", EMAX
cz        PRINT*,"UMAX =",UMAX
cz        PRINT*,"J= ",J
cz        PRINT*,"K= ",K	
      ENDIF

      if((EMAX .GT. EMAXMX) .AND. (N .GT. 1)
     |  .AND. (J .LE. NZchemMXE)) THEN
        PRINT*,"Error 4.1, EMAX > limit: coupledSUA terminated; EMAX = "
        PRINT 6949, EMAX
        PRINT*,"J = ",J
        PRINT*,"UNUM(K,J) = "
     	PRINT 6949, UNUM(K,J)
        PRINT*,"USOL(K,J) = "
     	PRINT 6949, USOL(K,J)
        WRITE(6942,6943) 4.1
	     IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
            ENDIF
        GOTO 2911
      ENDIF
      
      IF (UNUM(K,J) .le. 0.) THEN
      PRINT*, " "
cz      PRINT*, "AFTER NEWTON: NEGATIVE UNUM at (I,J):",K,ISPEC(K),J
      PRINT 7101, K,ISPEC(K),J
7101  FORMAT("After Newton, Negative UNUM at Species Number:",
     |I3,3X,"Species Name:",2X,A4,3X, "Grid Point:",I3)
      PRINT*, " "

C        IF ( (K .eq. LOH) .AND. (J .ge. 20) ) THEN
C           UNUM(K,J) = SMALL
C        ELSE
           STOP
C        ENDIF
      ENDIF
2223  CONTINUE
C
1010  CONTINUE   ! ALTITUDE LOOPS END************************
C 
CZ Set levels of Newton species above NH to SMALL
      DO J=NH+1,NZ1  ! hydro56 was NZ
       DO K=NQ1,NQT
        UNUM(K,J) = SMALL  !UNUM(K,NH)
        USOL(K,J) = SMALL  !USOL(K,NH)
       ENDDO
      ENDDO
C
CZ Above Newton limit NZchemMXN, do Chemplone for O3, then OH, HO2

      DO 6923 J=NH+1,NZ1   ! ALTITUDE LEVELS ABOVE NEWTON METHOD ! hydro56 was NZ
      
CZ DO OH first, then HO2 and O3  
      DO 6922 I=17,15,-1
         UNUM1(I) = UNUM(I,J-1)

         CALL CHEMPLONE(UNUM1,J,I,XP,XL)
           IF (chemkluge .eq. 1 .and. XL .LT. SMALL) then
                 PRINT*, "WARNING, XL = ",XL
                 XL = SMALL
           endif	 
	   
         UNUM1(I) = XP/XL

 6922  CONTINUE

 6923  CONTINUE

C-TF  END OF NEWTON STEPS
1616  CONTINUE
C
            IF (IZDEBUG.EQ.1) PRINT*, "CALLING DOCHEM"
      CALL DOCHEM(FVAL,N)
C
 1613 FORMAT(1x,"for ",A8,2x,"YP, YL:",I4, 1P10E10.2)
C      
C      PRINT 1614, A(55,NZ)*UNUM(LH,NZ)*UNUM(LOI,NZ)/
C     | A(72,NZ)/UNUM(LO,NZ), A(55,NZ), A(72,NZ)
 1614 FORMAT(1x,"xihp=",1P10E10.2)     
C
C      PRINT 9001, N,(TL(I),I=1,NQ)
9001  FORMAT (/1x,"After CALLING DOCHEM ",I3,1x,1P13E10.2)

        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0500, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-0500   LOOSELY-COUPLED SPECIES
      IF(IZDEBUGEXO.EQ.1) PRINT*,"Starting chemistry, i_exo =", i_exo
C ***** SOLVE FOR LOOSELY COUPLED SPECIES USING A 
C       TRIDIAGONAL INVERSION *****
C
      MZ=NZ
      MZ1 = MZ - 1

cz      PRINT*,"In main, MZ,MZ1=",MZ,MZ1
      DO 58 I=1,NQ    ! Tridiag species loop
C
C   COMPUTE ADVECTION TERMS FOR PARTICLES
C   TA = LOWER DIAGONAL, TB = DIAGONAL, TC = UPPER DIAGONAL, TY =
C   RIGHT-HAND SIDE
      DO 70 J=1,NZM ! hydro56; was MZ
      TA(J) = 0.
      TB(J) = 0.
      TC(J) = 0.
  70  TY(J) = 0.
C
      IF (IZDEBUG .GT. 0) THEN
      PRINT*,"YL(I,MZ1)=",YL(I,MZ1)
      ENDIF
      
      DO 44 J=1,MZ1 ! hydro56; was NZ
      TB(J) = YL(I,J)
C
C      IF (I .eq. LOI) PRINT 1019,I,J,Z(J),YP(LOI,J)
C     |  ,YL(I,J), DU(I,J), DL(I,J)
 1019 FORMAT ("YP(O+)=",2I5,1P10E10.2)      
C
CZ hydro64 use RHO_M instead of RHO
      TY(J) = YP(I,J)*WGT(I)/RHO_M(J)
C
  44  CONTINUE
            IF (IZDEBUG.EQ.1) PRINT*, "CONTINUING FROM 44"

cz      DO J = 1, MZ
cz        RHOR2(J) = RHO(J)*R(J)**2
cz        RHOR2_M(J) = RHO_M(J)*R_M(J)**2
cz      ENDDO

cz      DO J = 2, MZ-1
cz       RratioL(J) = RHOR2(J)/(RHOR2_M(J)*DR(J)*(DR(J-1)+DR(J)))
cz       RratioU(J) = RHOR2(J+1)/(RHOR2_M(J)*DR(J)*(DR(J)+DR(J+1)))
cz      ENDDO

CZ hydro62 added TA1,TB1,TC1
CZ hydro64 separated DR(J from D variables in Difco, multiplied here, to fix bug found in derivation
      DO 45 J=2,MZ-1 ! hydro56 changed to MZ-2, hydro63 back to MZ-1
      TA(J) = - DL(I,J) + (DHL(I,J) + DTL(I,J))*DR(J)
      TB(J) = TB(J) + DD(I,J)
     |      + (DHL(I,J) + DTL(I,J))*DR(J-1)
     |      - (DHU(I,J) + DTU(I,J))*DR(J+1)
      TC(J) = - DU(I,J) - (DHU(I,J) + DTU(I,J))*DR(J)

     IF(FLUXDBG.EQ.1)THEN
      
       TA1(I,J) = TA(J)
       TB1(I,J) = TB(J)
       TC1(I,J) = TC(J)

       K = I
       FT1(K,J) = 1/(DR(J)+DR(J-1))
       FT3A(K,J) = 2*(USOL(K,J) - USOL(K,J-1))*DK(K,J)
       FT3B(K,J) = DTK(K,J) + HD(K,J) 
       FT3C(K,J) = DR(J)*USOL(K,J-1)+DR(J-1)*USOL(K,J)
       FT3(K,J) = FT3A(K,J) + FT3B(K,J)*FT3C(K,J)

       FLUXKZ1(K,J)= -FT1(K,J)*FT3(K,J)

         IF(K .EQ. LHE) THEN
         PRINT*,"FLUX(K,J)=",FLUXKZ1(K,J)," J=",J,"species = ",ISPEC(K)
         ENDIF

      ENDIF

cz      TA(J) = RratioL(J)*(-2*DK(I,J)+(DTK(I,J)+HD(I,J)*DR(J)))
cz
cz      TC(J) = RratioU(J)*(-2*DK(I,J+1)-(DTK(I,J+1)+HD(I,J)*DR(J)))
cz
cz      TB(J) =  RratioL(J)*(2*DK(I,J)+(DTK(I,J)+HD(I,J))*DR(J-1))
cz     | + RratioU(J)*(2*DK(I,J+1)-(DTK(I,J+1)+HD(I,J+1))*DR(J+1))
cz     | + YL(I,J)


      TA2(I,J) = TA(J)
      TB2(I,J) = TB(J)
      TC2(I,J) = TC(J)


  45  CONTINUE
C      PRINT*, "BEFORE SETTING UP BCs, @N=", N, "I=",I
C      
C   BOUNDARY CONDITIONS
C ***** LOWER BOUNDARY CONDITIONS *****
      LB = LBOUND1(I)
C
C   CONSTANT DEPOSITION VELOCITY
      IF(LB.NE.0) GO TO 16 ! if LB=0 then have constant deposition vel.
cz  VDEP(I)=0   !this statement set all dep vel to zero; also data input is VDEP0, not VDEP
CZ Add conditional statements for different planets?

      TB(1) = TB(1) + DU(I,1) - VDEP(I)/DR(1) ! ?? Need to multiply some terms by correct DR(
     |  - DHU(I,1) - DTU(I,1)
      TC(1) = - DU(I,J) - (DHU(I,J) + DTU(I,J))*DR(J) 
      GO TO 15
C
C   CONSTANT MIXING RATIO   ?
  16  IF(LB.NE.1) GO TO 31
cz      TB(1) = DU(I,1) - DHU(I,1) - DTU(I,1) ! Replaced w/ below per JFK
       TB(1) = TB(2)  
CZ Added
CZ? Want to set TA(1)=TC(1)=0 for all boundary conditions??
       TA(1) = 0
       TC(1) = 0
       if (IZBCDEBUG .ne. 0) then
       print*," "
       PRINT*,"TA(1) and TC(1) set to zero for  ",ISPEC(I)
       endif
CZ  
C      TC(1) = 0.
      TY(1) = TB(1) * USOL(I,1)
      GO TO 15
C
C   CONSTANT UPWARD FLUX: SGFLUX in unit of cm-2 s-1
  31  CONTINUE
      TB(1) = TB(1) + DU(I,1) - DHU(I,1) - DTU(I,1)
      TC(1) = - DU(I,1) - DHU(I,1) - DTU(I,1)
      TY(1) = TY(1) + SGFLUX(I)*WGT(I)/DR(1)/RHO(1)
  15  CONTINUE
            IF (IZDEBUG.EQ.1) PRINT*, "CONTINUING FROM 15"

            IF (IZBCDEBUG .gt. 0) then	    
	     PRINT*,"TA(1)",TA(1)
	     PRINT*,"TB(1)",TB(1)
	     PRINT*,"TC(1)",TC(1)
	     PRINT*,"TY(1)",TY(1)	     
            ENDIF	     
C
      TA1(I,1) = TA(1)
      TB1(I,1) = TB(1)
      TC1(I,1) = TC(1)

C ***** UPPER BOUNDARY CONDITIONS *****
CZ IF IZUBC=0, Set all upper B.C. to const eff velocity, and all VEFF = 0
        IF (IZUBC .EQ. 0) then
          VEFF(I) = 0
          MBOUND(I)=0 !0=constant effusion velocity	  
        endif  
CZ	
      MB = MBOUND(I)
cz      PRINT*,"1938 EMAX =",EMAX
C   CONSTANT EFFUSION VELOCITY
      IF(MB.NE.0) GO TO 29
C-TF  COMPUTE THE JEANS ESCAPE VELOCITY HERE
C
      u02i = bol_k*T(NZ)/WGT(I)
      ve = sqrt(2.*bigG*bigM/R(NZ))
      veu = ve - U_ADV(NZ)
      veu2 = veu*veu
      alampi = veu2/(2.*u02i)
CZ moved rootpi to denom per JFK email Nov 2011      
cz      ujeansi = rootpi/2. * exp(-alampi)*sqrt(2.*u02i)*(1.+alampi)
      ujeansi = 1./(2.*rootpi) * exp(-alampi)*sqrt(2.*u02i)*(1.+alampi)

      flux_jeans=ujeansi*UNUM(I,NZ1) ! hydro56; was NZ
C      dif_limited_flux=VEFF(I)*UNUM(I,1)
C
      phie = 2.8e8*UNUM(LH,NZ1)*T(NZ)/1000.*flux_jeans/(5e7*2.75e4) ! hydro56; was NZ
C
C      PRINT 1047, I, flux_jeans, dif_limited_flux, phie+1.36*flux_jeans
C     | , T(NZ), UNUM(LH,NZ), (phie+1.36*flux_jeans)/UNUM(LH,NZ), ujeansi
 1047 FORMAT(1x,"Compare Jeans:", I5, 1P10E10.2)
C
      VEFF(I)=ujeansi  ! Jeans effusion velocity
CZ2
      if (N .GE. NSTEPS-1) THEN
      print*,"VEFF set to ",VEFF(I)," for ",ISPEC(I) !," species # ",I
      ENDIF
CZ      
      IF (NN .eq. NSTEPS) PRINT 1727, ISPEC(I), UNUM(I,NZ1)  ! hydro56 was NZ
     | , VEFF(I), VEFF(I)*UNUM(I,NZ1)
     | , ujeansi, flux_jeans
     | , phie/UNUM(I,NZ1), phie
 1727 FORMAT(/,1x,"for ",A8,1x,'n_i=',1P1E10.2,/
     | ,1x,"VEFF=",1P1E10.2,3x, "F_i=", 1P1E10.2,/
     | ,1x,'Jeans V=',1P1E10.2,3x, "Jeans F=",1P1E10.2,/
     | ,1x,'nonthermal V=',1P1E10.2,3x, 'nonthermal F=',1P1E10.2,/)    

      if(TB(MZ1).EQ.TB(MZ1)*0.5+1) then
      print*,"at 2578, TB(MZ1)=",TB(MZ1)
      endif

CZ ! hydro56; was MZ
cz      TA(MZ1) = - DL(I,MZ1)  + DHL(I,MZ1) + DTL(I,MZ1) ?? Put this back in??
      TB(MZ1) = TB(MZ1) + DL(I,MZ1) + VEFF(I)/DR(MZ1)
     |   + DHL(I,MZ1) + DTL(I,MZ1)

      IF (TB(J).EQ.TB(J)*0.5+1) THEN
        PRINT*,"WARNING, TB(J)=",TB(J)
        PRINT*,"J=",J
        print*,"MZ1=",MZ1
        PRINT*,"DL(I,J)=",DL(I,J)
        PRINT*,"DHL(I,J)=",DHL(I,J)
        PRINT*,"DTL(I,J)",DTL(I,J)	
        print*,"DR(J)=",DR(J)
        PRINT*,"RHO(J)=",RHO(J)
        PRINT*,"VEFF(J)=",VEFF(J)	
      ENDIF

      TB(MZ) = TB(MZ1)
      TY(MZ) = TY(MZ1)

      IF (IZDEBUG.EQ.1) then
      print*,"DR(MZ1)=",DR(MZ1)
      PRINT*,"VEFF(I)=",VEFF(I)
      PRINT*,"TB(MZ1)=",TB(MZ1)
      endif
      
      GO TO 30
C
C   CONSTANT UPWARD FLUX: SMFLUX in unit of cm-2 s-1
  29  CONTINUE
            IF (IZDEBUG.EQ.1) PRINT*,"CONTINUING FROM 29"

CZ ! hydro56; was MZ
cz      TA(MZ1) = - DL(I,MZ1)  + DHL(I,MZ1) + DTL(I,MZ1)

      print*,"at 2605, TB(MZ1)=",TB(MZ1)

      TB(MZ1) = TB(MZ1) + DL(I,MZ1) + DHL(I,MZ1) + DTL(I,MZ1)
      TY(MZ1) = TY(MZ1) - SMFLUX(I)*WGT(I)/DR(MZ1)/RHO(MZ1)
C
      TB(MZ) = TB(MZ1)
      TY(MZ) = TY(MZ1)

C      PRINT 1009, I, MZ, TY(MZ), SMFLUX(I)*WGT(I)/DR(MZ)/RHO(MZ)
C     | ,TY(MZ)+SMFLUX(I)*WGT(I)/DR(MZ)/RHO(MZ)
 1009 FORMAT ("in sua.f: SMFLUX", 2I5, 1P10E10.2)
  30  CONTINUE

        IF (IZBCDEBUG .gt. 0) then	    
          PRINT*,"TA(MZ1)=",TA(MZ1)
          PRINT*,"TB(MZ1)=",TB(MZ1)
          PRINT*,"TC(MZ1)=",TC(MZ1)
          PRINT*,"TY(MZ1)=",TY(MZ1)	     
        ENDIF	 

            IF (IZDEBUG.EQ.1) PRINT*, "CONTINUING FROM 30"  

CZ cs114  CHECK FOR NaNs in TA, TB, TC
            IF (IZDEBUG.EQ.1) PRINT*, "MZ=",MZ        
      DO J = 1, MZ-1
      IF (TA(J).EQ.TA(J)*0.5+1 .AND. PROD.NE.2.) THEN
        PRINT*,"TA(J)=",TA(J)
        PRINT*,"J=",J
        TA(J)=0
        PRINT*, "Reset to TA=",TA(J)
      ENDIF	
      IF (TB(J).EQ.TB(J)*0.5+1) THEN
        PRINT*,"WARNING, TB(J)=",TB(J)
        PRINT*,"J=",J
        PRINT*,"DL(I,J)=",DL(I,J)
        PRINT*,"DHL(I,J)=",DHL(I,J)
        PRINT*,"DTL(I,J)",DTL(I,J)	
        print*,"DR(J)=",DR(J)
        PRINT*,"RHO(J)=",RHO(J)
        PRINT*,"VEFF(J)=",VEFF(J)	
        TB(J)=0
        PRINT*, "Reset to TB=",TB(J)
      ENDIF	
      IF (TC(J).EQ.TC(J)*0.5+1) THEN
        PRINT*,"TC(J)",TC(J)
        PRINT*,"J=",J
        TC(J)=0
        PRINT*, "Reset to TC=",TC(J)
      ENDIF	
      enddo
      
      TB(MZ) = TB(MZ1)
      TY(MZ) = TY(MZ1)

CZ 1=O, 2=O2, 3=N2, 4=He, 5=H, 6=H2, 7=CO2
      IF(IZFLUXDEBUG.EQ.100 .OR. 
     |  IZFLUXDEBUG.EQ.I) THEN
        PRINT*," "
        PRINT*,"SPECIES NUMBER I=",I,"  SPECIES:  ",ISPEC(I)

        PRINT 1708
CZ ! hydro56; was NZ/MZ	
        DO J=1,NZ1       
        PRINT 1709,J,TA(J),TB(J),TC(J),TY(J),DL(I,J),DHL(I,J),DTL(I,J)
     |    ,USOL(I,J),UNUM(I,J)	
        ENDDO
	
 1708 FORMAT(4x,"J",6x,"TA",8x,"TB",8x,"TC",8x,"TY"
     | ,8X,"DL",8X,"DHL",7X,"DTL",6x,"USOL",6X,"UNUM")	
 1709 FORMAT(I4,1X,1P9E10.2)
 
      ENDIF      

C-TF  CHEMICAL EQUILIBRIUM: Ci=TY/TB*WGT(I)/RHO

 1214 FORMAT(1x,"@NN=",I5,2x,"Species",2x,A8,2x,"UNUM(",I3," )="
     | ,2x,1P10E10.2)
C
      IF ( (1.eq.1) .AND.
     |     ( ( NN .eq. NSTEPS ) .OR. ( NN .eq. 1) ) .AND.
C     |    (NN .ge. 400) .AND.
     |     (I .eq. LH) 
     |   ) THEN
C       IF ( I .eq. LO) THEN
C      IF ( ( I .eq. LINERT ) .OR. (I .eq. LH) ) THEN
      PRINT 1213, N, I, ISPEC(I), NZ1
 1213 FORMAT(/1x,"ELEMENTS IN THE TRIDIAGONAL MATRIX: @N=",I5
     | , "I=",I5,2x,A8,2x,"NZ1=",I4)
C      PRINT*, " J     TA      TB         TC        TY       CH       
C     |  USOL"
      PRINT 1215
 1215 FORMAT(20x,"CHEM",8x,"FUP",7x,"FLOW",7x,"DU",7x,"DHU",7x,"DTU"
     | ,7x,"DL",7x,"DHL",7x,"DTL",7x,"RHS",7x,"DUSOL1",4x,"DUSOL2")
C
      JJJ=1
      DUSOL = USOL(I,JJJ+1) - USOL(I,JJJ)
      IF (abs(DUSOL/USOL(I,JJJ)) .le. SMALL) DUSOL=0
      TF_transport1(JJJ) =
     | - DU(I,JJJ)*DUSOL
     | - DHU(I,JJJ)*(USOL(I,JJJ+1) + USOL(I,JJJ))
     | - DTU(I,JJJ)*(USOL(I,JJJ+1) + USOL(I,JJJ))
      TF_transport2(JJJ) = 0.
      TF_chem(JJJ) = YP(I,JJJ)*WGT(I)/RHO(JJJ)-YL(I,JJJ)*USOL(I,JJJ)

CZ Added flag
       IF (IZPRNTFLG4.eq.1) PRINT 999,JJJ,TF_chem(JJJ),
     |       TF_transport1(JJJ), TF_transport2(JJJ)
C      PRINT 999, JJJ, TF_chem, TA(JJJ), TB(JJJ), TC(JJJ), TY(JJJ)
C      
      DO JJJ=2,NZ1 ! hydro56; was NZ
      DUSOL = USOL(I,JJJ+1) - USOL(I,JJJ)
      IF (abs(DUSOL/USOL(I,JJJ)) .le. SMALL) DUSOL=0
      TF_transport1(JJJ) =
     | - DU(I,JJJ)*DUSOL
     | - DHU(I,JJJ)*(USOL(I,JJJ+1) + USOL(I,JJJ))
     | - DTU(I,JJJ)*(USOL(I,JJJ+1) + USOL(I,JJJ))
C     
      DUSOL2 = USOL(I,JJJ) - USOL(I,JJJ-1)
      IF (abs(DUSOL2/USOL(I,JJJ)) .le. SMALL) DUSOL2=0
      TF_transport2(JJJ) = 
     | - DL(I,JJJ)*DUSOL2
     | - DHL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
     | - DTL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
C
      TF_chem(JJJ) = YP(I,JJJ)*WGT(I)/RHO(JJJ)-YL(I,JJJ)*USOL(I,JJJ)
C
CZ Added flag
      IF (IZPRNTFLG4.eq.1) PRINT 999,JJJ,TF_chem(JJJ), 
     | TF_transport1(JJJ), TF_transport2(JJJ)
     | , -DU(I,JJJ)*DUSOL, -DHU(I,JJJ)*(USOL(I,JJJ+1) + USOL(I,JJJ))
     | , -DTU(I,JJJ)*(USOL(I,JJJ+1) + USOL(I,JJJ))
     | , -DL(I,JJJ)*DUSOL2, -DHL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
     | , -DTL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
     | , -TF_transport1(JJJ)+TF_transport2(JJJ)
     | , DUSOL, DUSOL2
C     | , YP(I,JJJ),WGT(I),RHO(JJJ),YL(I,JJJ),USOL(I,JJJ)
C      PRINT 999, JJJ, TF_chem, TA(JJJ), TB(JJJ), TC(JJJ), TY(JJJ)
      ENDDO
      JJJ=MZ1 ! hydro56; was MZ
C      TF_transport1(JJJ) = SMFLUX(I)*WGT(I)/DR(MZ)/RHO(MZ)
      TF_transport1(JJJ) = VEFF(I)*USOL(I,MZ1)/DR(MZ1) ! hydro56; was NZ
      DUSOL2 = USOL(I,JJJ) - USOL(I,JJJ-1)
      IF (abs(DUSOL2/USOL(I,JJJ)) .le. SMALL) DUSOL2=0
      TF_transport2(JJJ) =
     | - DL(I,JJJ)*DUSOL2
     | - DHL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
     | - DTL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
      TF_chem(JJJ) = YP(I,JJJ)*WGT(I)/RHO(JJJ)-YL(I,JJJ)*USOL(I,JJJ)
C
CZ Added flag
      IF (IZPRNTFLG4.eq.1) PRINT 999,JJJ,TF_chem(JJJ),
     | TF_transport1(JJJ), TF_transport2(JJJ)
     | , 0., 0., 0.
     | , -DL(I,JJJ)*DUSOL2, -DHL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
     | , -DTL(I,JJJ)*(USOL(I,JJJ-1) + USOL(I,JJJ))
     | , -TF_transport1(JJJ)+TF_transport2(JJJ)
     | , 0., DUSOL2
C     
C      PRINT 999, JJJ, TF_chem, TA(JJJ), TB(JJJ), TC(JJJ), TY(JJJ)
C     
      PRINT 1215
      PRINT*, ''
      ENDIF
 999  FORMAT (1x,"CHEM vs. DIF:",I3, 1P16E10.2)

      DO J=1,MZ1 ! hydro56; was NZ
      TSOL(J) = TY(J)   ! SAVE TY FOR DEBUG
      ENDDO
CZ
           cretry = 0
 7683 CONTINUE
CZ
      IF ( (I .ne. LOI) .AND. (I .ne. LNI) ) THEN
           IF (IZDEBUG.EQ.1) PRINT*, "CALLING SGTSL 1"
           IF ((i_exo.GE.599) .AND. (IZDEBUG.EQ.1))  PRINT*, "MZ = ", MZ 

CZ ! hydro56; was NZ, MZ
           CALL SGTSL(MZ1,TA,TB,TC,TSOL,NFLAG)
           IF (IZDEBUG.EQ.1) THEN
	     PRINT*, "RETURNING FROM SGTSL 1"
             PRINT 7195, MZ1,NFLAG 
 7195     FORMAT("MZ-1 = ",I,"       NFLAG = ",I)
             PRINT*," " 
	     IF (NFLAG.NE.0) THEN
                PRINT*,"TA = ",TA
                PRINT*," "   
	        PRINT*,"TB = ",TB
                PRINT*," " 
                PRINT*,"TC = ",TC
                PRINT*," " 
	        PRINT*,"TSOL = ",TSOL
              ENDIF
           ENDIF	     
      ENDIF
 
CZ ?? May want to add option to run these all the way down to level 1

      IF (I .eq. LOI) THEN  ! level 13 is the lower boundary for O+
        IF (NZ_INI .eq. 57) THEN
        NTMP=12
        NTMP1=NTMP+1
              IF (IZDEBUG.EQ.1) PRINT*, "CALLING SGTSL 2"	
CZ ! hydro56; was NZ/MZ	      
        CALL SGTSL(MZ1-NTMP,TA(NTMP1:NZ1),TB(NTMP1:NZ1),TC(NTMP1:NZ1)
     | , TSOL(NTMP1:NZ1),NFLAG)
C
C        CALL SGTSL(MZ-12,TA(13:NZ),TB(13:NZ),TC(13:NZ)
C     | , TSOL(13:NZ),NFLAG)
        ENDIF
C
        IF (NZ_INI .eq. 67) THEN
              IF (IZDEBUG.EQ.1) PRINT*, "CALLING SGTSL 3"	
CZ ! hydro56; was NZ/MZ	      
        CALL SGTSL(MZ-22,TA(23:NZ1),TB(23:NZ1),TC(23:NZ1)
     | , TSOL(23:NZ1),NFLAG)
        ENDIF
C			
        IF (NZ_INI .eq. 77) THEN
              IF (IZDEBUG.EQ.1) PRINT*, "CALLING SGTSL 4"	
CZ ! hydro56; was NZ/MZ	      
        CALL SGTSL(MZ-32,TA(33:NZ1),TB(33:NZ1),TC(33:NZ1)
     | , TSOL(33:NZ1),NFLAG)
        ENDIF
C			     
      ENDIF
C
      IF (I .eq. LNI) THEN  ! set level 20 as the lower boundary for N+ may have discontinuity problem
        IF (NZ_INI .eq. 57) THEN
        NTMP=0 !19
cz        IF (chemkluge .EQ. 1) NTMP = 19
        NTMP1=NTMP+1
            IF (IZDEBUG.EQ.1) PRINT*,"CALLING SGTSL 5"	
	    
CZ ! hydro56; was NZ/MZ	    
        CALL SGTSL(MZ1-NTMP,TA(NTMP1:NZ1),TB(NTMP1:NZ1),TC(NTMP1:NZ1)
     | , TSOL(NTMP1:NZ1),NFLAG)
C     
C       do J=1,NZ
C       print 4101, NN,J, TSOL(J), USAVE(I,J), TSOL(J)/USAVE(I,J)
 4101  format(1x,"after calling SGTSL for N+:",2I4,1P10E10.2)       
C       enddo
        ENDIF ! INI .eq. 57
C
        IF (NZ_INI .ne. 57) THEN
         PRINT*, "lower boundary of N+ not specified"
         STOP
        ENDIF ! ne. 57
      ENDIF ! I .eq. LNI 

C      CALL TRIDAG_LEI(TSOL,MZ,1,MZ,TA,TB,TC,TY)
C      
      IF (NFLAG.NE.0) THEN
      PRINT 400, N,NFLAG,I,ISPEC(I)
      NTF_STOP_FLG = 1
        PRINT*,"TIME STEP N = ",N
        WRITE(6942,6943) 3.
CZ Added to debug tridiag solver failure	
        PRINT*," "
        PRINT*,"SPECIES NUMBER I=",I,"  SPECIES:  ",ISPEC(I)

        PRINT 1708	
        DO J=1,NZ1  ! hydro56 was NZ      
        PRINT 1709,J,TA(J),TB(J),TC(J),TY(J),DL(I,J),DHL(I,J),DTL(I,J)
     |    ,USOL(I,J),UNUM(I,J)	
        ENDDO

      PRINT*," "
      print*, "Lower Boundary Conditions:"
      print 1710,LBOUND1

      print*, "Upper Boundary Conditions:"
      print 1710,MBOUND
      
      PRINT*,"VEFF="
      PRINT 1711,VEFF
      
 1710 FORMAT(15I3)      
 1711 FORMAT(1P15F7.2)
      

        IF (chemkluge .eq. 1 .and. cretry .le. 5) Then
	     DO J = 1, NZ1  ! hydro56 was NZ 
	       TA(J) = TA(J) + SMALL*10.**cretry
	       TB(J) = TB(J) + SMALL*10.**cretry
	       TC(J) = TC(J) + SMALL*10.**cretry
	       TSOL(J) = TSOL(J) + SMALL*10.**cretry	       
	     ENDDO
             PRINT*,"retry count = ", cretry
	     cretry = cretry + 1
	     
        GOTO 7683
        endif
        GOTO 2911
      ENDIF
 400  FORMAT(//1X,"TRIDIAGONAL SOLVER FAILED AT N =",I3,2X,
     2  "NFLAG =",I2,2X,"SPECIES #",I2,2x,"SPECIES: ",A6)
C
CZ Altitude Loop for Chemistry

      DO 59 J=1,NZ1 ! hydro56; was NZ
CZ      

      IF (USOL(LH,J) .LT. SMALL) THEN
      print*,"WARNING: @3274 USOL(H) = ",USOL(LH,J)
      PRINT*,"J = ",J
      USOL(LH,J) = SMALL
      ENDIF
      IF (USOL(LH2,J) .LT. SMALL) THEN
      print*,"WARNING: @3279 USOL(H2) = ",USOL(LH2,J)
      USOL(LH2,J) = SMALL
      PRINT*,"J = ",J      
      ENDIF          
CZ
C      IF (I .eq. LH) THEN

      IF ((TSOL(J) .LE. -0.001)) THEN !.AND.(I.NE.2).AND.(I.NE.8)) THEN
        if(IZPRNTFLG7 .EQ. 1) THEN
         PRINT*,"WARNING: TSOL(J) = ",TSOL(J)
         PRINT*,"Species I = ",I
         PRINT*,"Altitude J = ",J
        ENDIF
        TSOL(J) = -0.001
      ENDIF

CZ Damping
        USOL(I,J)=TSOL(J)*(1-DAMPING) + USAVE(I,J)*DAMPING
C      
      IF (NZ_INI .eq. 57) THEN
        IF ( (I .eq. LOI) .AND. (J .le. 13) ) THEN  ! IGNORE THE VARIATION OF O+
        USOL(I,J)=USAVE(I,J)
        GO TO 1615  
        ENDIF
      ENDIF 

 998  FORMAT(1x, I3,1P14E10.2)
C 
CZ Added to control changes in USOL and UNUM
      IF (USAVE(K,J) .LT. SMALL) then
          USAVE(K,J) = SMALL
      ENDIF
      
      delUSOL = (USOL(K,J)-USAVE(K,J))/USAVE(K,J)
      EREL = ABS(delUSOL)
      
      delUNUM = (UNUM(K,J)-USAVEUN(K,J))/USAVEUN(K,J)
      
      IF (delUNUM .GT. 1.0) THEN
         UNUM(K,J) = 2.0*USAVEUN(K,J)
      ENDIF

      IF (delUNUM .LT. -0.5) THEN
         UNUM(K,J) = 0.5*USAVEUN(K,J)
      ENDIF

CZ Added to limit UNUM > 0
      if (UNUM(K,J) .LT. 1) UNUM(K,J) = 1
      USOL(K,J)=UNUM(K,J)*WGT(K)/RHO(J)

      IF (USOL(LH,J) .LT. SMALL) THEN
      print*,"WARNING: @3395 USOL(H) = ",USOL(LH,J)
      PRINT*,"J = ",J
      USOL(LH,J) = SMALL
      ENDIF
      IF (USOL(LH2,J) .LT. SMALL) THEN
      print*,"WARNING: @3400 USOL(H2) = ",USOL(LH2,J)
      USOL(LH2,J) = SMALL
      PRINT*,"J = ",J      
      ENDIF        

CZ Added to look for two altitudes set to UNUM = 1 for same species
      IF(UNUM(K,J) .EQ. 1 .AND. UNUM(K,J-1) .EQ. 1
     | .AND. (N .GT. 10) .AND. J .GE. 2
     | .AND. K .LE. 14) THEN
        PRINT*,"Error 9, UNUM(K,J) = UNUM(K,J-1) =1, J =" 
        PRINT*, "Alt. Grid J = ",J
	PRINT*, "Species K = ", K
        WRITE(6942,6943) 9.
            IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)  
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)     
            ENDIF
cz        WRITE(6942,6949) USOL(LH2,J)
        GOTO 2911      
      ENDIF
CZ

CZ Limit emax to grid points at or below NZchemMXLL
CZ For rev cs_r28.f, replaced NZchemMXE with min UNUM()
      IF ((EREL .GT. EMAX) .AND. (UNUM(I,J) .GE. 1.E10) 
cz     |    .AND. (J.LE.NZchemMXE)
     |   .AND. (USOL(I,J) .GE. 1.E-4))   THEN
        UMAX = USOL(I,J)
        IS = I
        JS = J
        UMAX = USOL(I,J)
        EMAX = EREL
      ENDIF

       if((EMAX .GT. EMAXMX) .AND. (N .GT. 1)
     |      .AND. (J .LE. NZchemMXE)) THEN
        PRINT*,"Error 4.2, EMAX > limit: coupledSUA terminated; EMAX = "
                PRINT 6949, EMAX
        PRINT*,"J = ",J
        PRINT*,"UNUM(I,J) = "
	        PRINT 6949, UNUM(I,J)
        PRINT*,"USOL(I,J) = "
	        PRINT 6949, USOL(I,J)
        WRITE(6942,6943) 4.2
	     IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
            ENDIF
        GOTO 2911
      ENDIF
C     
 1615 CONTINUE
CZ       PRINT*, "CONTINUING FROM 1615" 

      IF (USOL(I,J) .LE. 0) THEN
CZ Removed stop for USOL<0, replaced with previous TF fix
        USOL(I,J) = SMALL  
CZ
      IF(USOL(LH2,J) .le. 0) THEN
        PRINT*,"Error 7, H2 CONC. < 0, USOL(LH2,J) = ",USOL(LH2,J)
        PRINT*, "J = ",J
        WRITE(6942,6943) 7.
            IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)  
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)
            ENDIF
        WRITE(6942,6949) USOL(LH2,J)
        GOTO 2911      
      ENDIF
CZ
      IF(USOL(LH,J) .le. 0) THEN
        PRINT*,"Error 8, H CONC. < 0, USOL(LH,J) = ",USOL(LH,J)
        PRINT*, "J = ",J
        WRITE(6942,6943) 8.
        WRITE(6942,6949) USOL(LH,J)
             IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)  
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
            ENDIF
        GOTO 2911      
      ENDIF
CZ
CZ      PRINT 1930,ISPEC(I),USOL(I,J),J,USAVE(I,J),TSOL(J)
 1930 FORMAT(/1x,"NEGATIVE USOL(",A8," )=",1P1E10.2,2x," at J=",I5
     | ,2x,"USAVE=",1PE10.2,2x,"TSOL=",1PE10.2)     
C
 1900 FORMAT(1x,I5, 1P10E10.2)      
CZ      STOP
C      
 1991 CONTINUE      
C      
      ENDIF   ! USOL < 0!!!
C      
  59  CONTINUE ! Tridiag altitude loop J
C
   58 CONTINUE  ! TRI-DIAGONAL SPECIES LOOP I

CZ Added per JK to provide 3-point smoothing of H and H2
      IF (SMOOTH3 .EQ. 1) THEN
      DO J=1, NZ1  ! hydro56 was NZ   
      UOLDH(J) = USOL(LH,J)
      UOLDH2(J) = USOL(LH2,J)
      ENDDO
      
      DO J=2,NZ-2  ! hydro56 was NZ-1
      USOL(LH,J) = 0.25*UOLDH(J-1) + 0.5*UOLDH(J) +
     |  0.25*UOLDH(J+1)      
      USOL(LH2,J) = 0.25*UOLDH2(J-1) + 0.5*UOLDH2(J) +
     |  0.25*UOLDH2(J+1)       
       ENDDO
       ENDIF
CZ      
 992  FORMAT(5x, "J",5x,"TSOL",8x,"Te",8x,"Ti",8x,
     | "T",8x,"RHO",8x,"WT",7x,"USOL")
C
      ISP = ISPEC(IS)
      ZMAX = zk(JS)
CZ2
      if (u_inf_ini .ge. 0) then
        PRINT 100, N,EMAX,ISP,ZMAX,JS,psimax,ipsimaxgrid,DT,TIME,EMAXK,
     |   imaxk,i_exo,u_inf
      endif
CZ2      
      if (u_inf_ini .lt. 0) then
        PRINT 1733, N,EMAX,ISP,ZMAX,JS,psimax,ipsimaxgrid,DT,TIME,EMAXK,
     |   imaxk,i_exo,u_inf,ujeans(NZ)
      endif      
CZ2
CZ TEST
      IF (IZDEBUG.EQ.1) CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
     
 100   FORMAT(1X,"N=",I5,1X,"EMAX=",1PE9.2," FOR ",A4,
     2  1x,"AT Z=",E8.2,1X,"J=",1x,I3,2x,"psimax=",0P,F6.3,2X,
     3  "@J= ",I3,2X,"DT=",1PE9.2,2X,"TIME=",1PE9.2
     4  ,2x,'EMAXK=',1P1E10.2,"  @J =",I4,"  exobase @J =",I4
     5  ,2x,"u_inf=",1P1E10.2)

CZ2     
 1733  FORMAT(1X,"N=",I5,1X,"EMAX=",1PE9.2," FOR ",A4,
     2  1x,"AT Z=",E8.2,1X,"J=",1x,I3,2x,"psimax=",0P,1P1E10.2,2X,
     3  "@J= ",I3,2X,"DT=",1PE9.2,2X,"TIME=",1PE9.2
     4  ,2x,'EMAXK=',1P1E10.2,"  @J=",I4," exobase @J =",I4      
     5  ,2x,"u_inf=",1P1E10.2," ujeans(NZ)=",1P1E10.2)
CZ2
C
 317  CONTINUE

        DO J = NZ-1,NZ
        if(DEN(J).EQ.DEN(J)*0.5+1) THEN
          print*,"WARNING, before NORM_MMR, DEN(J) =",DEN(J)
          print*,"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO 

C
C-TF  BECAUSE USOL OF LONG LIVED SPECIES AND NEWTON SPECIES 
C-TF  MAY HAVE CHANGED DURING THE CURRENT TIME STEP, WE NEED 
C-TF  TO NORMALIZE THE MMR (USOL), RECALCULATE NUMBER DENSITY 
C-TF  (UNUM & DEN), RECOMPUTE MEAN MOLECULAR WEIGHT (WT). ALL 
C-TF  THESE ARE DONE BY ASSUMING CONSTANT TOTAL MASS DENSITY.
C
Cz 0 -- NORMALIZE USOL AND COMPUTE DEN,UNUM,WT,WTP FROM THE NEW USOL
      CALL NORM_MMR(0)   ! NORMALIZE USOL
C 
      if(EDEBUG.GT.0)PRINT*,"UNe=",(UNe(J),J=1,NZ)


      DO J=1,NZ-1 !hydro56 was NZ; hydro61: changed to J+1, added scale height to scale up to bdry
      xno(J+1)=UNUM(LO,J) *exp(-DR(J)/(2*HSCALE(J))) 
      xnno(J+1)=UNUM(LNO,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnn2(J+1)=UNUM(LN2,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xno2(J+1)=UNUM(LO2,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnhe(J+1)=UNUM(LHE,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnh(J+1)=UNUM(LH,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnh2(J+1)=UNUM(LH2,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnh2o(J+1)=UNUM(LH2O,J)*exp(-DR(J)/(2*HSCALE(J)))  
      xnoh(J+1)=UNUM(LOH,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnco2(J+1)=UNUM(LCO2,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xnco(J+1)=UNUM(LCO,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xne(J+1)=UNe(J)*exp(-DR(J)/(2*HSCALE(J))) 
      xno3(J+1)=UNUM(LO3,J)*exp(-DR(J)/(2*HSCALE(J)))  
      xno21d(J+1)=UNUM(LO21d,J)*exp(-DR(J)/(2*HSCALE(J)))  
      xno21s(J+1)=UNUM(LO21s,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xio2p(J+1)=UNUM(LO2I,J)*exp(-DR(J)/(2*HSCALE(J)))  
      xinop(J+1)=UNUM(LNOI,J)*exp(-DR(J)/(2*HSCALE(J)))  
      xiop(J+1)=UNUM(LOI,J)*exp(-DR(J)/(2*HSCALE(J))) 
      xihp(J+1)=UNUM(LHI,J)*exp(-DR(J)/(2*HSCALE(J)))
      xnn4s(J)=UNUM(LN4S,J)*exp(DR(J)/(2*HSCALE(J)))       
      ENDDO

      J=1 
      xno(J)=UNUM(LO,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnno(J)=UNUM(LNO,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnn2(J)=UNUM(LN2,J)*exp(DR(J)/(2*HSCALE(J))) 
      xno2(J)=UNUM(LO2,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnhe(J)=UNUM(LHE,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnh(J)=UNUM(LH,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnh2(J)=UNUM(LH2,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnh2o(J)=UNUM(LH2O,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnoh(J)=UNUM(LOH,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnco2(J)=UNUM(LCO2,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnco(J)=UNUM(LCO,J)*exp(DR(J)/(2*HSCALE(J))) 
      xne(J)=UNe(J)*exp(DR(J)/(2*HSCALE(J))) 
      xno3(J)=UNUM(LO3,J)*exp(DR(J)/(2*HSCALE(J))) 
      xno21d(J)=UNUM(LO21d,J)*exp(DR(J)/(2*HSCALE(J))) 
      xno21s(J)=UNUM(LO21s,J)*exp(DR(J)/(2*HSCALE(J))) 
      xio2p(J)=UNUM(LO2I,J)*exp(DR(J)/(2*HSCALE(J))) 
      xinop(J)=UNUM(LNOI,J)*exp(DR(J)/(2*HSCALE(J))) 
      xiop(J)=UNUM(LOI,J)*exp(DR(J)/(2*HSCALE(J))) 
      xihp(J)=UNUM(LHI,J)*exp(DR(J)/(2*HSCALE(J))) 
      xnn4s(J)=UNUM(LN4S,J)*exp(DR(J)/(2*HSCALE(J))) 

      if(EDEBUG.gt.0)PRINT*,"xne=",(xne(J),J=1,NZ)

      IF (NTF_STOP_FLG .eq. 1) GO TO 22
C
      DO 2000 J=1,NZ  ! hydro56 was NZ; hydro61 back to NZ, added xn__
C
C      UNUM(LHE,J) = UNUM_INI(LHE,J)
C      
      vmrn2(J)=xnn2(J)/DEN(J)
      vmro2(J)=xno2(J)/DEN(J)
      vmro(J)=xno(J)/DEN(J)
      vmrh2(J)=xnh2(J)/DEN(J)
      vmrh(J)=xnh(J)/DEN(J)
      vmrn(J)=xnn4s(J)/DEN(J)
 2000 CONTINUE

C    
C-TF  ADD SUBSONIC MODEL HERE TO SOLVE FOR T, RHO, and U_ADV
C
C-TF  maybe this is where the recalculation can be added (What does Feng mean??)
      DO 1718 J=1,NZ ! hydro59 back to NZ
      T_save(J) = T(J)
      Te_save(J) = Te(J)
      Ti_save(J) = Ti(J)
      RHO_save(J) = RHO(J)
      DEN_save(J) = DEN(J)
      U_ADV_save(J) = U_ADV(J)
      WT_save(J) = WT(J)
      fluxN_save(J) = fluxN(J)
CZ Added fluxM
      fluxM_save(J) = fluxM(J)

      if(DEN(J).EQ.DEN(J)*0.5+1)THEN
      PRINT*,"WARNING: before call to tf_cond, DEN(J)=",DEN(J),
     |  "J =",J      
      ENDIF


 1718 continue
 
C 
C   RETRY TIME STEP IF EMAX EXCEEDS 30 PERCENT
 1717 continue
C
      dt_save = DT
      TIME = TIME + DT

        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0600, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-0600  THERMAL CONDUCTIVITIES AND HEATING/COOLING RATES
C-TF  Compute neutral and electron thermal conductivity 
C
            IF (IZDEBUG.EQ.1) PRINT*,"CALLING tf_cond"
CZ hydro61 shortened call, use xn__ variables in fmodules      
      CALL tf_cond(T, Te, DEN, xne, cond, conde)
C
      CALL tf_cp(NZ,T,vmrn2,vmro,vmro2,vmrh2,vmrh,vmrn,cp)    ! hydro59 changed to NZ-1; hydro61 back to NZ
      
CZ hydro59 added
cz      cp(NZ) = cp(NZ-1) !*RHO(NZ-1)/RHO(NZ-2)  ! hydro60 mod; use ratio of mass densities, but not at NZ in case not yet calc

C      
C   Calculate heating and cooling rates. Note the geometrical heating
C   effect: more energy is available at large distances.
C
            IF (IZDEBUG.EQ.1) PRINT*,"CALLING tf_heating"

      CALL tf_heating(glbmean_heating_flg,NN,i_exo
     | ,eflux,sigh2,nsteps,IZPRNTFLG11)
C
      DTINV = 1./DT
      
            IF (IZDEBUG.EQ.1) PRINT*,"CALLING setup_tridag"      
      CALL setup_tridag(DTINV,NN,NSTEPS,A_EDDY,B_EDDY,R_PLT)
C
 3013 FORMAT(1x,"before solving T:",I3,1P10E10.2)      
 
       IF (IZDEBUG.EQ.1) then
         print*,"NZ =",NZ
         print*,"T(1) =",T(1)	       
         PRINT*,"AK =",AK
         PRINT*,"BK =",BK
         PRINT*,"CK =",CK
         PRINT*,"rhsk =",rhsk
     	 PRINT*,"CALLING TRIDAG"
       endif
	    
      call TRIDAG(ak,bk,ck,rhsk,delT,NZ)  ! CZ changed back to NZ
      
C      
      CALL setup_tridag_lei(DTINV)
C
            IF (IZDEBUG.EQ.1) PRINT*, "CALLING TRIDAG_LEI, NZ = ",NZ ! CZ changed back to NZ
CZ      IF (NZ.GE.599) PRINT*, "C = ",cke

      if(TEDEBUG.GE.1) PRINT*,"before TRIDAG_LEI, Te(NZ)=",Te(NZ),
     |    "Te(NZ-1)=",Te(NZ-1)

      CALL TRIDAG_LEI(delTe,NZ,1,NZ,ake,bke,cke,rhske) ! CZ changed back to NZ
C
      if(TEDEBUG.GE.1) PRINT*,"after TRIDAG_LEI, delTe(NZ)=",delTe(NZ),
     |    "  delTe(NZ-1)=",delTe(NZ-1),"  NZ=",NZ

C   Compute new temperatures and maximum change in T at any height
      emaxk = 0.
      do 2556 i=1,NZ  ! CZ changed back to NZ

       Te(i) = delTe(i)**(1./3.5)
       
      if(TEDEBUG.GE.1) PRINT*,"before DAMPING, Te(NZ)=",Te(NZ),
     |    "  Te(NZ-1)=",Te(NZ-1),"  NZ=",NZ       
       
CZ Damping       
       Te(i) = Te(i)*(1-DAMPING) + Te_save(i)*DAMPING
CZ
CZ limit temperature change delTe to no more than +/- 10% of Te
CZ if tempkluge = 1

      IF ((Te(i) - Te_save(i))/Te_save(i) .LT. -0.1 
     | .AND. tempkluge .EQ. 1) THEN
        Te(i) = 0.9 * Te_save(i)
      ENDIF	

      IF ((Te(i) - Te_save(i))/Te_save(i) .GT. 0.1 
     | .AND. tempkluge .EQ. 1) THEN
        Te(i) = 1.1 * Te_save(i)
      ENDIF	
    
       
C       GO TO 2557   ! Ignore restriction on DT by the variation of Te       
       chg = abs(Te(i)-Te_save(i))/Te_save(i)
CZ Add NZtempMX       
        if ((chg.gt.emaxk).AND.(i.LE.NZtempMX)) then
       emaxk = chg
       imaxk = i
       endif
 2557 CONTINUE       
C       
C
CZ limit temperature change delT to no more than +/- 10% of T
CZ delT/abs(delT) gives the correct sign
      if (abs(delT(i)/T(i)) .GT. 0.1 .AND. abs(delT(i)) .GT. 1E-20
     |   .AND. tempkluge .EQ. 1)     then
            delT(i) = 0.1*T(i)*(delT(i)/abs(delT(i)))
      endif

      chg = abs(delT(i)/T(i))
CZ Add NZtempMX      
      if ((chg.gt.emaxk).AND.(i.LE.NZtempMX)) then
      imaxk = i
      emaxk = chg
      endif
CZ
      if((EMAXK .GT. EMAXKMX) .AND. (N .GT. 1)) THEN
      PRINT*,"Error 5.2, EMAXK > limit: coupledSUA terminated; EMAXK = "
     |   ,EMAXK
        PRINT*,"Level = ",i
        PRINT*,"delT(i) = ",delT(i)
        PRINT*,"T(i) = ", T(i)	
        WRITE(6942,6943) 5.2
              IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)  
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
            ENDIF
        GOTO 2911
      ENDIF
CZ

      if(delT(i).ge.50.) then
       print*,"WARNING, @3327, delT(i)=",delT(i)
       print*,"T(i)=",T(i)
       print*,"ak(i)=",ak(i)
       print*,"bk(i)=",bk(i)
       print*,"ck(i)=",ck(i)
       print*,"rhsk(i)=",rhsk(i)
       print*,"i=",i,"NZ=",NZ
       delT(i) = 50.
       print*,"setting delT(i) to",delT(i)
      endif

      T(i) = T(i) + delT(i)

CZ Damping      
      T(i) = T(i)*(1-DAMPING) + T_save(i)*DAMPING
CZ
      IF((T(i) .LT. 30) .and. (N .GT. 2)) THEN
         PRINT*,"Error 6.1:  T(i) < 30 K, T(i) = ",T(i)
         PRINT*,"i = ",i
         WRITE(6942,6943) 6.1
        WRITE(6942,6949) T(i)
          IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
          ENDIF      
      GOTO 2911 
      ENDIF
CZ
cz temp increased for cs_r39; later disabled
      IF((T(i) .GT. 1000) .and. (i .LT. 30) .and. 0 .eq. 1) THEN
         PRINT*,"Error 6.2:  T(i) > 1000 K, T(i) = ",T(i)
         PRINT*,"i = ",i
         WRITE(6942,6943) 6.2
        WRITE(6942,6949) T(i)
          IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
          ENDIF      
      GOTO 2911 
      ENDIF
    
C-TF COMPUTE Ti first assuming equilibrium
C-TF need to add in ion-neutral chemical reaction heating later (may be important)
C
      IF ( Te(i) .lt. T(i) ) Te(i) = T(i)
      Ti(i) = ( clei(i)*Te(i) + clin(i)*T(i)
C     |                        + qjoul(i)*RHO(i)
     |        ) / ( clei(i) + clin(i) )
C
      if (Ti(i) .eq. T(i)*0.5+1) then
cz       print*,"Error, Ti(i)=",Ti(i),"i=",i
cz       print*,"clei(i)=",clei(i)
cz       print*,"clei(i)=",clin(i)       
       if (Ti(i-1).lt. 1500) Ti(i) = Ti(i-1)
       if (Ti(i-1).ge. 1500) then
         Ti(i) = Ti(i-2)
         Ti(i-1) = Ti(i-2)
        endif
       print*,"Ti(i) NaN, set to",Ti(i)
      endif


 3000 format(1x,i3,2x,1p16e10.2)
C       
      IF (T(i) .lt. 100) THEN
        PRINT*," T(i) = ",T(i)
        PRINT*,"At grid point i = ", i
      ENDIF
 2556 continue   ! update T, Te, and Ti
C
 3111 format(/1x,'delT    delTe')

 3112 FORMAT(1x, 1P2E10.2)
C
CZ Added if statement - used only when in Jeans escape mode (u_inf_ini < 0)
CZ2 equiv to line number 1947 in TF original code
CZ2 Added u_inf_ini<-2 condition, sets minimum u_inf 
      IF (u_inf_ini .lt. 0) THEN
        u_inf=0.
        DO I=1,NQ
          IF (MBOUND(I) .eq. 0) THEN 
            u_inf = u_inf + VEFF(I)*USOL(I,NZ1) ! hydro56 was NZ
            IF (IZPRNTFLG8 .EQ. 1) 
     |      print*,ISPEC(I)," added to u_inf; u_inf= ",u_inf
          ENDIF
        ENDDO
	
        if ((u_inf_ini .eq. -2).and.(u_inf .lt. 0.1)) THEN
        	u_inf=0.1
            PRINT*,"u_inf below minimum, set to 0.1"
       endif
	
         if (u_inf_ini .le. -10. .and. u_inf .lt. -u_inf_ini) then
          u_inf = -u_inf_ini
          print*,"u_inf below minimum, set to ",u_inf
        endif
      ENDIF
      
      IF(DEBUGUINF .EQ. 1) PRINT*,"Line 4069 u_inf =",u_inf      
CZ
CZ Begin loop to control psimax using u_inf
      ipsimaxcount = 0

 4701 continue
 
         DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0700, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO
 
CZ-0700  CALL SUBSONIC TO CALC WIND AND DENSITY
CZ Added fluxM
      CALL SUBSONIC(N,u_inf,rootpi,ujeans,fluxN,fluxM,psi,bigM)
CZ Above line equiv to original line 1965

CZ Find the maximum psi and grid point of max psi
      psimax = 0.
      ipsimaxgrid = 0
      
      DO J=1,NZ
      IF(psi(J).gt.psimax) THEN
         psimax = psi(J)
         ipsimaxgrid = J
      ENDIF	 
      ENDDO
      
CZ Check temperatures at levels from 1 to NZ, make sure minimum is greater than 30 K
      Tz = T(1:NZ) ! CZ changed back to NZ
      MASK1(1:NZ) = .TRUE.  ! CZ changed back to NZ
      MASK1(NZ+1:NZM) = .FALSE.  ! CZ changed back to NZ+1    
      
      IF((MINVAL(Tz,MASK1) .LT. 30) .and. (N .GT. 2)) THEN

         PRINT*,"Error 6.1  T < 30 K, Tmin = ",MINVAL(Tz,MASK1)
         PRINT*,"Location of Tmin = ",MINLOC(Tz,MASK1)
cz         PRINT*," T = ",Tz	 
        PRINT*,"MASK1 = ",MASK1 	 
        WRITE(6942,6943) 6.1
        WRITE(6942,6949) MINVAL(Tz,MASK1)
            IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
            ENDIF
        GOTO 2911
      ENDIF
      
CZ Count and limit the number of times through the control loop controlling psimax using u_inf
        ipsimaxcount = ipsimaxcount + 1
CZ ipsimaxcount was 100, changed to 10 for r76	
        if (ipsimaxcount .gt. 40) THEN
         PRINT*,"Exceeded max loop count for psi, count =",ipsimaxcount
         GOTO 4702
        ENDIF
CZ Control psi using u_inf
      IF (u_inf_ini .GT. 0) THEN
        IF(psimax.gt.psihi + 0.05) then
           u_inf = 0.98 * u_inf 
           GOTO 4701
        endif   
       IF(psimax.gt.psihi) then	
          u_inf = 0.99 * u_inf
           GOTO 4701
        endif
        IF(psimax.lt.psihi-0.02 .and. N .GT. 100) then	
          u_inf = 1.01 * u_inf
          GOTO 4701	  
        endif
        IF(psimax.lt.psihi-0.01 .and. N .GT. 300) then	
          u_inf = 1.005 * u_inf
          GOTO 4701	  
        endif
        IF(psimax.lt.psihi-0.003 .and. N .GT. 500) then	
          u_inf = 1.003 * u_inf
          GOTO 4701	  
        endif

        IF(psimax.lt.psihi-0.1) then
          u_inf =  1.02 * u_inf
          GOTO 4701	  
        endif

        IF(psimax.lt.psihi*0.99) then
          u_inf =  1.02 * u_inf
          GOTO 4701	  
        endif

        IF(psimax.gt.psihi*1.01) then
          u_inf =  0.98 * u_inf
          GOTO 4701	  
        endif

       ENDIF

      IF(DEBUGUINF .EQ. 1) PRINT*,"after psimax control, u_inf =",u_inf

 4702 CONTINUE

      do J = 1, NZ
      IF (psi(J) .ge. 0.98) THEN
        PRINT*,"Error 1 - Psi > .98, term cs_xx; Psi = ",psi(J)
        print*, "J =",J
        PRINT*, "N = ",N
        WRITE(6942,6943) 1.
        WRITE(6942,6949) psi(NZ)
        IF (IERRDEST .EQ. 1) THEN
               PRINT*, "Calling Output"
              CALL OUTPUT(N,NSTEPS,0.,IZPRNTFLG12,EMAX,EMAXK)
              PRINT*, "Calling SUBSONIC_END_OUTPUT"
             CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |       delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)	      
        ENDIF
        GOTO 2911
      ENDIF
      enddo

      chg = 0.
      DO J=1,NZ  ! hydro59 back to NZ
C      
      RHO(J) = RHO(J)*(1-DAMPING) + RHO_save(J)*DAMPING
      chg = abs(RHO(J)-RHO_save(J))/RHO_save(J)

CZ hydro58; calc new DEN from new RHO calc in subsonic
      DEN(J) = RHO(J)/WT(J)
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0780, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"RHO_SAVE =",RHO_SAVE
          print*,"WT(J) =",WT(J)
        ENDIF



CZ Added condition to limit emaxk to below max grid point
      if((chg.gt.emaxk).AND.(J.le.NZtempMX)) then
      imaxk = J
      emaxk = chg
      endif
      ENDDO
C      
 2013 CONTINUE
CZ This is the subsonic.out output (unit 101) hydro56 changed some to NZ1
      write(101,2033) NN,DT,NZ,zk(NZ),U_ADV(NZ),T(NZ),RHO(NZ1),
     | fluxN(NZ),
     | emaxk,imaxk,T(imaxk),T_save(imaxk),Te(imaxk),Te_save(imaxk)
 2033 format(1x,'L=',i5,2x,'dt=',1pe8.1,2x,'N=',i3,2x,'z=',1pe8.1,
     |  2x,'u=',e8.1,2x,'T=',e9.2,2x,'rho=',e9.2,2x,'flux=',e8.1,2x,
     |  'emaxk=',e8.1, ' at i=',i3, 1P4E10.2)
CZ2
      IF(EXODEBUG.EQ.1) THEN
       print*,"current exobase =",i_exo
       PRINT*,"current u_inf =",u_inf
       PRINT*,"iexohold =",iexohold
      ENDIF 
CZ2
       IF (iexohold .eq. 0 .or. N .lt. iexohold) then ! call find_exobase & adjust_exobase only if 
CZ2       exobase hold is off, iexohold=0 or code has run steps up to iexohold

        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-0800, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-0800 FIND NEW EXOBASE
C-TF  COMPUTE exobase level in order to readjust the top of the grid
CZ added this statement, moved back after find_exobase:  IF (i_exo.lt.600) THEN
      
      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
         PRINT*,"3448 i_exo = ", i_exo      
         PRINT*,"CALLING find_exobase"
      ENDIF
       
       i_exo_chg = 0
CZ2 
       i_exo_savez = i_exo
CZ2       
      CALL find_exobase(i_exo,i_exo_flg,ipsimaxgrid,u_inf_ini)

      IF(EXODEBUG.EQ.1) THEN
        Print*,"Returning from find_exobase, i_exo =",i_exo
      ENDIF	
cz      PRINT*,"U_ADV(NZ) =",U_ADV(NZ)

      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
      PRINT 6956,i_exo,R(i_exo)/1E5,
     |     (R(i_exo)-R_PLT)/1E5
         PRINT*,"Exobase at NZ=",NZ     
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF

 6956 FORMAT( "Exobase found at level i_exo = ",I4,2X,
     | "R(J) = ",F6.0," km,  Alt = ",F6.0," km")      
CZ
      IF (i_exo.lt.i_exo_max) THEN
CZ
      IF ( (i_exo .ne. NZ) .OR. (i_exo_flg .ne. 1) ) THEN   !change in exobase
       NZ_OLD = NZ
       i_exo_chg = 1
cz         PRINT 6952,i_exo,NZ,i_exo_flg,NN
6952  FORMAT(8X,"Before Exobase Adjusted: i_exo = ",I4,2X,
     | "NZ = ",I4,2X,"i_exo_flg = ",I4,2X,"NN = ",I4)
CZ
CZ Added fluxM
CZ uslope used to extrapolate to point above old exobase
cz        uslope =  (zl(NZ)-zl(NZ-1))/(U_ADV(NZ)-U_ADV(NZ-1))

       CALL adjust_exobase(i_exo,rootpi,fluxN,fluxM,IZ300,IZ300A,bigM
     |  ,ipsimaxgrid)
CZ removed control for i_exo < 10  
      IF (i_exo .ge. 10) then ! allow i_exo to vary if <10
        IF (i_exo .gt. i_exo_savez+1) i_exo = i_exo_savez+1
      endif
      IF (i_exo .lt. i_exo_savez-1) i_exo = i_exo_savez-1      
CZ2     
       NZ1 = NZ-1
       NZ2 = NZ1-1

      IF (u_inf_ini .ge. 0) THEN  ! u_inf_ini > 0, u_inf is controlled to achieve psimax = 0.6
CZ This was our original method to impose u_infinity on the top of the atmosphere
cz         U_ADV(NZ) = u_inf
CZ Original exobase level read in from restart file is NZREAD
CZ Previous exobase level is NZ_OLD
CZ Original u_infinity is u_inf_ini
CZ Adjusted u_infinity is u_inf_adj
CZ These routines assume that NEXODELT = 1 (exobase can only change +/-1)

         IF(EXODEBUG.EQ.1) PRINT*,"Returning from adJust_exobase" 
	 
         IF (NZ .EQ. NZ_OLD+1) THEN  !exobase moved up 1
	     IF(EXODEBUG.EQ.1) PRINT*,"Exobase moved up 1, NZ =",NZ
	      delU_ADV = (U_ADV(NZ_OLD-1)-U_ADV(NZ_OLD))
	      delZnew = (zl(NZ_OLD)-zl(NZ_OLD+1)) 
	      delZold = (zl(NZ_OLD-1)-zl(NZ_OLD))
	      u_inf = u_inf + delU_ADV*delZnew/delZold
	      IF(EXODEBUG.EQ.1)THEN
	        print*,"delU_ADV =",delU_ADV
	        print*,"delZnew =",delZnew	      
	        PRINT*,"delZold =",delZold
                print*,"new u_inf =",u_inf
              ENDIF		
         ENDIF     
	 
         IF (NZ .EQ. NZ_OLD-1) THEN  !exobase moved down 1
           IF(EXODEBUG.EQ.1) then
	      PRINT*,"Exobase moved down 1, NZ =",NZ
	      print*,"NZ_OLD =",NZ_OLD
            ENDIF  
              u_inf = U_ADV(NZ_OLD-1)
cz            u_inf = u_inf * 0.97 
           IF(EXODEBUG.EQ.1) print*,"new u_inf =",u_inf	      
         ENDIF
	 
        U_ADV(NZ) = u_inf
	
      ENDIF ! u_inf_ini > 0, control u_inf to achieve psimax = 0.6

6953  FORMAT(8X,"After Exobase Adjusted:  i_exo = ",I4,2X,
     | "NZ = ",I4,2X,"i_exo_flg = ",I4,2X,"NN = ",I4)
CZ
      ENDIF   !change in exobase
CZ      
      ENDIF  !exobase has not reached max level
CZ2
      ENDIF ! iexohold is off (exobase free to move)
CZ2

      IF(DEBUGUINF .EQ. 1) PRINT*,"Line 4356, u_inf =",u_inf

C      IF (NZ .gt. NZ_INI) GO TO 22
C
 2001 CONTINUE !seems to not be used except in commented out code

CZ-0900  ADJUST TIME STEP AND CHECK TIME AND # STEPS
C   Adjust the time step based on how fast the solution is changing
      dt0 = DT
CZ     if(emaxk.gt.3.e-2) DT = dt0*0.5
CZ      if(emaxk.gt.2.e-2) DT = dt0*0.9

      IF (NAN_DBG.EQ.1) PRINT*, "EMAXK = ", EMAXK
      
CZ Control DT based on emaxk      
      IF(emaxk.gt.1.0) then
          DT = dt0*0.1
         PRINT 7198, EMAXK, DT
         GO TO 7197
      ENDIF

      IF(emaxk.gt.2.e-1) then
          DT = dt0*0.5
         PRINT 7198, EMAXK, DT
         GO TO 7197
      ENDIF

      IF(emaxk.gt.1.e-1) then
          DT = dt0*0.9
         PRINT 7198, EMAXK, DT
         GO TO 7197
      ENDIF
      
CZ Control DT based on EMAX
      IF(emax.gt.3.e-1) then
          DT = dt0*0.2
         PRINT 7190, EMAX, DT
         GO TO 7197
      ENDIF

      IF(emax.gt.2.e-1) then
          DT = dt0*0.8
         PRINT 7190, EMAX, DT
         GO TO 7197
      ENDIF

      if(emaxk.lt.1.e-4) THEN
cz         DT = dt0*10.
         DT = dt0*2.0
      	 GO TO 7197
      ENDIF
      
       if(emaxk.lt.1.e-3) THEN
         DT = dt0*1.3
      	 GO TO 7197
      ENDIF
 
       if(emaxk.lt.2.5e-2) THEN
         DT = dt0*1.02
      	 GO TO 7197
      ENDIF
 
       if(emaxk.lt.4e-2) THEN
         DT = dt0*1.01
      	 GO TO 7197
      ENDIF
      
       if(emaxk.lt.6e-2) THEN
cz         DT = dt0*1.003
      	 GO TO 7197
      ENDIF

CZ
C      
 7197  IF(DT.GT.3E15) PRINT*,"DT = ", DT   
        if(DT .GT. DTMX) DT=DTMX
      
 7198 FORMAT("EMAXK = ",1P1E10.2,2X,"DT = ",1P1E10.2)      
 7190 FORMAT("emax = ",1P1E10.2,2X,"DT = ",1P1E10.2)   
 
 7150  CONTINUE

      IF (NAN_DBG.EQ.1) PRINT 7199 , DT
 7199 FORMAT("After dt0 adjustment, DT = ",1P1E10.2)      

CZ Don't change DT if DT <= 1E-10, instead increase grid size to 300
      if(DT .LE. 1E-10) THEN
        DT=1E-10
	     IZ300A = 0
      IF((IZ300.EQ.1).AND.(NN.GT.500).AND.(NZ.LT.NZMAX)
CZ2
     |    .and. (iexohold .eq. 0)) THEN
          IZ300A = 1
CZ Added IZ300,IZ300A,fluxM	  
         CALL adjust_exobase(i_exo,rootpi,fluxN,
     |   fluxM,IZ300,IZ300A,bigM,ipsimaxgrid)
        PRINT*,"Adjusting atmosphere ceiling, NZ =",NZ
       ENDIF	
        GO TO 281 !moves program forward to complete the run through the loop
      ENDIF
CZ

CZ IZDEF01=1 cancels timestep short
CZ 281 moves program forward to complete the run through the loop
      IF((emaxk.LT.0.33).OR.(IZDEF01.EQ.1)) GO TO 281
CZ if emaxk > 0.33, save old values and drop DT by factor of 10
      NZ = NZ_OLD
      NZ1 = NZ-1
      NZ2 = NZ1-1
C
      PRINT*,"emaxk>0.33, dropping timestep and reducing DT"
      DT = 0.1*dt_save
      TIME = TIME - dt_save
      
      do J=1,NZ  ! hydro59 back to NZ
      T(J) = T_save(J)
      Te(J) = Te_save(J)
      Ti(J) = Ti_save(J)
      RHO(J) = rho_save(J)
      DEN(J) = den_save(J)
      
      U_ADV(J) = U_ADV_save(J) !MEZ had commented out
      
      WT(J) = WT_save(J)
      fluxN(J) = fluxN_save(J)
      fluxM(J) = fluxM_save(J)      
C      DO I=1,NSP
C      USOL(I,J) = USAVE(I,J)
C      ENDDO

      WRITE(101, 2035), DT
      PRINT 2035, DT

 2035 FORMAT(1x,"Maintaining NZ, DT reduced to ", 1E9.2)

      enddo
 
      
CZ 
      IF (u_inf_ini .GE. 0) U_ADV(NZ) = u_inf
      
CZ 1717 jumps program back and drops DT by 10 if emaxk > 0.33
      GO TO 1717
C      
 281  CONTINUE
C      
       IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
       PRINT*,"Converting new USOL to UNUM"
      ENDIF
      
CZ hydro59 fixed and moved before call to NORM_MMR

      DO J=1,NZ-1 
        WT_M(J) = SQRT(WT(J+1)*WT(J))  
        WTP_M(J) = WT_M(J)/protonM
        RHO_M(J) = SQRT(RHO(J)*RHO(J+1))
        DEN_M(J) = RHO_M(J)/WT_M(J)
       enddo

C-TF  END OF SUBSONIC MODEL

cz      if(USOL(1,NZ-1).EQ.USOL(1,NZ-1)*0.5)THEN
cz       PRINT*,"WARNING, USOL(1,NZ-1) = ",USOL(1,NZ-1)
cz       PRINT*,"NZ-1=",NZ-1       
cz      ENDIF

C                2 -- CONVERT USOL INTO UNUM
       CALL NORM_MMR(2)

      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF

C-TF  Jump out of the time loop when maximum timesteps reached
      IF(NN.EQ.NSTEPS) THEN
      print*, "NSTEPS exceeded!"
      GO TO 22 
      ENDIF
C
      IF ( (TIME.GT.TSTOP) .AND. (NN .GT. 1000) ) THEN
C JFK Let's let it stop at TSTOP
CZ Changed back to previous way, requires at least 1000 time steps
CZ      IF (TIME.GT.TSTOP) THEN
      print*, "tstop exceeded!"
      NN = NSTEPS - 1
      ENDIF
C
      WRITE(850,3212) NN,TIME,NZ
      
      DO J=1,NZ1 ! hydro56 was NZ
      WRITE(850,3113) zk(J),RHO(J),WTP(J),T(J),Te(J),Ti(J)
     | ,q(J),cond(J),cp(J),UNe(J)
     | ,qntot_glb(J)
     | ,(UNUM(I,J),I=1,NSP)
      ENDDO
C
 3212 FORMAT(I5,1P1E10.2)      
 3113 FORMAT(1P11E16.6)      
 
        DO J = NZ-1,NZ
        IF(DEN(J).EQ.DEN(J)*0.5+1)THEN
          PRINT*,"WARNING CZ-1000, DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
        ENDDO

CZ-1000 CONTROLS ON CHANGES IN CHEMISTRY 
      IF (chemkluge .eq. 1) THEN

      IF (USOL(LH,J) .LT. SMALL) THEN
      print*,"WARNING: @4470 USOL(H) = ",USOL(LH,J)
      PRINT*,"J = ",J
      USOL(LH,J) = SMALL
      ENDIF
      IF (USOL(LH2,J) .LT. SMALL) THEN
      print*,"WARNING: @4475 USOL(H2) = ",USOL(LH2,J)
      USOL(LH2,J) = SMALL
      PRINT*,"J = ",J      
      ENDIF          
CZ
      IF (USAVE1(K,J) .LT. SMALL) then
          USAVE1(K,J) = SMALL
      ENDIF
      
      delUSOL = (USOL(K,J)-USAVE1(K,J))/USAVE1(K,J)
      EREL = ABS(delUSOL)
      
      delUNUM = (UNUM(K,J)-USAVEUN1(K,J))/USAVEUN1(K,J)
      
      IF (delUNUM .GT. 0.5) THEN
         UNUM(K,J) = 1.5*USAVEUN1(K,J)
      ENDIF

      IF (delUNUM .LT. -0.5) THEN
         UNUM(K,J) = 0.5*USAVEUN1(K,J)
      ENDIF
CZ      

CZ Added to limit UNUM > 0
      if (UNUM(K,J) .LT. 1) UNUM(K,J) = 1
      USOL(K,J)=UNUM(K,J)*WGT(K)/RHO(J)
      
CZ  
      ENDIF ! chemkluge = 1

      IF (EXO_DBG.EQ.1) THEN ! Added to debug NaN problems when exobase moves
         PRINT*,"Exobase at NZ=",NZ      
       DO I=1,12
         PRINT*,"UNUM(I,NZ-1)=",UNUM(I,NZ-1),"I=",ISPEC(I)
       enddo
      ENDIF

            IF (IZDEBUG.EQ.1) PRINT*,"END OF TIME-STEPPING LOOP"
      IF(DEBUGUINF .EQ. 1) PRINT*,"End of time loop, u_inf =",u_inf
C
   1  CONTINUE
C

CSSS ***** END THE TIME-STEPPING LOOP *****
C
  22  CONTINUE
 
      WRITE(850,3212) NSTEPS,TIME,NZ1 ! hydro56 was NZ
      DO J=1,NZ1
      WRITE(850,3113) zk(J),RHO(J),WTP(J),T(J),Te(J),Ti(J)
     | ,q(J),cond(J),cp(J),UNe(J)
     | ,qntot_glb(J)
     | ,(UNUM(I,J),I=1,NSP)
      ENDDO
      CLOSE(850)
C  
      PRINT*, "Calling Output"
      CALL OUTPUT(N,NSTEPS,TIME,IZPRNTFLG12,EMAX,EMAXK)
C
      write(101,2008) TIME
 2008 format(/1x,'TIME =',1pe10.3)
C
CZ Added fluxM, WTP
      PRINT*, "Calling SUBSONIC_END_OUTPUT"
      CALL SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |  delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     |       dif_limited_flux)
C  
C-PK Write to file used for TPF spectra
      WRITE(79,230)
 230    FORMAT(/1X,"Alt [km]",3X,"Pres [Pa]",3X,"Temp [K]",4X,"f(O3)",
     2    5X,"f(CO2)",5X,"f(H2O)"/)
      DO J=NZ1,1,-1  ! hydro56 was NZ
        PRES_MKS(J) = DEN(J)*bol_k*T(J)*0.1
        WRITE(79,231) J,PRES_MKS(J),T(J)
      END DO
 231    FORMAT(3X,I3,4X,6(1PE11.3))
      WRITE(79,230)
C
C print out P&L tables with integrated rxn rates, "int.rates.out.dat"
      DO 702 I=1,NSP
         ISP = ISPEC(I)
         WRITE(150,703) ISP,TP(I)
 703     FORMAT(/A8,12X,"PRODUCTION RXS",14X,"INT RX RATE",4X,
     2      "TP = ",1PE9.2)
       DO 704 K=1,NR 
          IF(JCHEM(3,K).EQ.I .OR. JCHEM(4,K).EQ.I .OR. 
     2       JCHEM(5,K).EQ.I)THEN
           IF(RAT(K).NE.0.) WRITE(150,705) I,(CHEM(J,K),J=1,5),RAT(K)
 705       FORMAT(1X,I3,1H),1X,A7,3H + ,A7,3H = ,A7,3H + ,A6,2X,A4,
     2      1PE10.3)
          ENDIF
 704   CONTINUE
C
         WRITE(150,706) ISP,TL(I)
 706     FORMAT(/A8,15X,"LOSS RXS",16X,"INT RX RATE",4X,"TL = ",1PE9.2)
       DO 707 K=1,NR 
          IF(JCHEM(1,K).EQ.I .OR. JCHEM(2,K).EQ.I)THEN
             IF(RAT(K).NE.0.) WRITE(150,705) K,(CHEM(J,K),J=1,5),RAT(K)
          ENDIF
 707   CONTINUE
 702  CONTINUE
C
C-TF  print out profiles: tf_profile.dat
C      PRINT*, "R_PLT=", R_PLT, "USH(1,1)=", USH(1,1)

cz      PRINT*, "DU,DL",DU(11,61),DL(11,61)
cz      PRINT*, "N2(NZ1):", NZ1, USOL(5,NZ1), NN

CZ  UNIT 80 == tf_profiles.dat
      WRITE(80,6900)
 6900 FORMAT(1x,"NN",5x,"NZ1",3x,"NQ",2x,"NQT",2x,"NR"
     1 ,3x,"NSP",1x,"lmax",1x,"nsigbin",3x,"R0",9x,"R_PLT",9x,
     2 "TIME")
CZ
      WRITE(80,900), NN,NZ,NQ,NQT,NR,NSP,lmax,nsigbin,R0,R_PLT,TIME
CZ
      WRITE(80,*)" "
      WRITE(80,*)"THE FOLLOWING ROW IS CURRENTLY NOT USED"
CZ
      WRITE(80,899), LO,   LO2,   LN2,   LHE,    LH,   LH2,  LCO2,   LCO
     1 ,           LN4S,   LNO,  LH2O,   LOI,   LNI,   LHI,   LO3,  LHO2
     2 ,            LOH,  LH2I,  LH3I, LOI2P, LOI2D,  LN2I,  LCOI, LCO2I
     3 ,           LO2I,  LNOI,  LOHI,  LN2D,  LO1D,    LC, LH2O2, LO21s
     4 ,          LO21d
 899  FORMAT (33I3)     
 900  FORMAT (8I5, 1P3E13.5)      
C
CZ
      WRITE(80,*)" "
      WRITE(80,6901)
 6901 FORMAT(5x,"R",8x,"DEN",7x,"RHO",8x,"T",9x,"Ti"
     1 ,8x,"Te",8x,"EDD",5x,"U_ADV",7x,"DR",7x,"DRC",7x,
     2 "WTP")
CZ 
      DO 901 J=1,NZ1  ! hydro56 was NZ
      WRITE(80,902),R(J),DEN(J),RHO(J),T(J),Ti(J),Te(J),EDD(J)
     | ,U_ADV(J),DR(J),DRC(J),WTP(J)
 901  CONTINUE
 902  FORMAT (1P11E10.3)
CZ
      WRITE(80,6901)
      WRITE(80,*)" "
CZ
      WRITE(80,*)" "
      WRITE(80,6902)

 6902 FORMAT(2x,"J",4x," O ",7x," O2",7x," N2",7x," HE"
     1 ,7x," H ",7x," H2",7x,"CO2",7x," CO",7x,"N4s",7x," NO",7x,
     2 "H2O",7x," OI",7x," NI",7x," HI",7x," O3",7x,"HO2",7x," OH",
     3 7x,"H2I",7x,"H3I",6x,"OI2p",6x,"OI2d",7x,"N2I",7x,"COI",6x,
     4 "CO2I",7x,"O2I",7x,"NOI",7x,"OHI",7x,"N2d",7x,"O1d",7x," C ",
     5 6x,"H2O2",6x,"O21s",6x,"O21d")
CZ 
      DO 903 J=1,NZ1  ! hydro56 was NZ
      WRITE(80,904),J,(USOL(I,J),I=1,NSP)
 903  CONTINUE
CZ
      WRITE(80,6902)
      WRITE(80,*)" "
CZ
      WRITE(80,*)" "
      WRITE(80,6902)

      DO 913 J=1,NZ1 ! hydro56 was NZ
      WRITE(80,904),J,(UNUM(I,J),I=1,NSP)
C      PRINT 904,J,(UNUM(I,J),I=1,NSP)
 913  CONTINUE      
CZ
      WRITE(80,6902)
      WRITE(80,*)" "
CZ
      WRITE(80,6903) 
      DO 933 J=1,NZ1  ! hydro56 was NZ
      WRITE(80,934),UNe(J) 
 933  CONTINUE

      WRITE(80,6903)
      WRITE(80,*)" "
 6903 FORMAT(3x,"UNe(J)")
CZ 
      WRITE(80,*)"WRITING TCOLUNUM"

      DO 943 J=1,NZ1  ! hydro56 was NZ
      WRITE(80,944),(TCOLUNUM(I,J),I=1,NQT)
 943  CONTINUE      
C 
      WRITE(80,*)" "
      WRITE(80,*)"WRITING COLUNUM"
CZ
      DO J=1,NZ1  ! hydro56 was NZ
      WRITE(80,944),(COLUNUM(I,J),I=1,NQT)
C      PRINT 904, J, COLUNUM(1,J), TCOLUNUM(1,J), UNe(J), UNUM(1,J)
      ENDDO
CZ
      WRITE(80,*)" "      
      WRITE(80,6934), sflux_lya
      WRITE(80,6944), f107
      WRITE(80,*)" "
CZ

 904  FORMAT (I3,1P33E10.3)
 934  FORMAT (1PE10.3)
 944  FORMAT (1P15E10.3)
CZ
 6934  FORMAT ("sflux_lya = ",1PE10.3)
 6944  FORMAT ("f107 = ",1P15E10.3)
CZ
C-TF  Energy of chemical reactions: tf_profile.dat
CZ
      WRITE(80,*)" "
      WRITE(80,6905)
 6905 FORMAT(2x,"R",2x,"TF_EV")
CZ 
      DO 912 I=1,NR
      WRITE(80,915), I, TF_EV(I)
C      print 915, I, TF_EV(I)
 912  CONTINUE

      WRITE(80,6905)
      WRITE(80,*)" "
CZ

 915  FORMAT (I3, f7.2)
C      PRINT*, "EV printed"
C
  888 FORMAT(1PE10.3)
C  
C-TF  heating rate profiles: tf_profile.dat
CZ
      WRITE(80,*)" "
      WRITE(80,6908)

 6908 FORMAT(2x,"J",5x,"o3pc",8x,"noc",9x,"metac",7x,"co2c"
     1 ,8x,"ircool",7x,"chem",7x,"len",9x,"npe",9x,"chap",9x,"har",8x,
     2 "herz",8x,"hug",9x,"srb",9x,"lya",10x,"src",9x,"joul",8x,"night",
     3 7x,"adv",9x,"cond",9x,"q",7x,"q-qircool")
CZ
      DO J=1,NZ
      WRITE(80,898), J, qo3pc(J),qnoc(J),qmetac(J),qco2c(J),qircool(J)
     |, qchem(J), qlen(J), qnpe(J), qchap(J), qhar(J), qherz(J), qhug(J)
     |, qsrb(J), qlya(J), qsrc(J), qjoul(J), qnight(J)
     |, q_adv(J), q_cond(J), q(J), q(J)-qircool(J)
C      PRINT 898, J, qo3pc(J),qjoul(J),qnight(J),q_cond(J)
      ENDDO
  898 FORMAT(I3, 1P21E12.3)  
C      PRINT*, "q printed"
      WRITE(80,6908)
      WRITE(80,*)" "
CZ
C-TF  solar radiation flux and wavelengths
CZ
      WRITE(80,*)" "
      WRITE(80,6906)
 6906 FORMAT(3x,"wave1",5x,"wave2",5x,"sflux")
CZ
      DO L=1,lmax
      WRITE(80,897), wave1(L), wave2(L), sflux(L)
C      PRINT 897, wave1(L), wave2(L), sflux(L)
      ENDDO
CZ
      WRITE(80,6906)
      WRITE(80,*)" "
CZ
CZ
      WRITE(80,*)" "
      WRITE(80,6907)
 6907 FORMAT(3x,"wv_low",4x,"wv_high",3x,"usflx")
CZ
      DO L=1,nsigbin
      WRITE(80,897), wv_low(L), wv_high(L), usflx(L)
      ENDDO
 897  FORMAT(1P3E10.3)
      WRITE(80,6907)      
CZ      
      CLOSE(80)
C 
C     TRANSPORT TERMS: tf_transport.dat
      DO 905 J=1,NZ1  ! hydro56 was NZ
      DO 905 I=1,NQ
      WRITE(82,906),USOL(I,J),DU(I,J),DL(I,J),DHU(I,J),DHL(I,J)
     | ,DTU(I,J),DTL(I,J),DI(I,J),DT
 905  CONTINUE
 906  FORMAT (1P8E10.2)
C
C     Chemical terms: tf_chem.dat
      DO 907 J=1,NZ-1  !   hydro56 was NZ
      DO 907 I=1,NQ
      WRITE(83,908) FVAL(I,J)*WGTM(I)/WTP(J)
 907  CONTINUE
 908  FORMAT (1PE10.3)
      DO 909 I=1,NQ
      WRITE(83,910) WGTM(I)
 909  CONTINUE
 910  FORMAT (1PE10.3)
C
CZ hydro61 added
      J=NZ
      UNe(J)=0
      DO I=1,NSP4
      UNUM(I,J)=0
      ENDDO

CZ hydro61 added zk(J) (Alt)
C-TF  WRITE RESULTS TO restart FILE: restart.out
      WRITE(810,502) NZ,NZ1,NZ2,NZ_INI,R0
      xx=0.
      DO J=1,NZ
      WRITE(810,501) R(J),zk(J),RHO(J),T(J),Ti(J),Te(J),EDD(J)
     | ,U_ADV(J),RBND(J),DR(J),DRC(J),drp(J),drm(J),WT(J),UNe(J)
     | ,(UNUM(I,J),I=1,NSP4)
      ENDDO

      write(810,504)
      write(810,503)
      
  503 FORMAT(5X,"R(J)",7X,"Alt",7X,"RHO(J)",6X,"T(J)",7X,"Ti(J)",6X
     |,"Te(J)",6X,"EDD(J)",3x,"U_ADV(J)",4X,"RBND(J)",5X,"DR(J)",5X
     | ,"DRC(J)",5X,"drp(J)",5X,"drm(J)",5X,"WT(J)",6X,"UNe(J)"
     | ,8x,"O",9x,"O2",9x,"N2",9x,"HE",10x,"H",9x,"H2",8x,"CO2"
     | ,9x,"CO",8x,"N4S"9x,"NO",8x,"H2O"
     | ,9x,"O+",9x,"N+",9x,"H+",9X,"O3",8X,"HO2",9X,"OH"
     | ,8x,"H2+",8x,"H3+",7x,"OI2P",7x,"OI2D",8x,"N2+",8x,"CO+"
     | ,7x,"CO2+",8x,"O2+",8x,"NO+"
     | ,8x,"OH+",8x,"N2D",8x,"O1D",10x,"C",7x,"H2O2",7x,"O21s"
     | ,7x,"O21d",9x,"HV",10x,"M",8x,"NUL",10x,"e")



 504  FORMAT(" ")

 501  FORMAT(1P52E11.4)
 502  FORMAT(4I4,1P1E10.3)
C 
      GO TO 21
  20  PRINT 300,I
 300  FORMAT(//1X,"NMAX EXCEEDED FOR SPECIES ",I3)
      GO TO 21
  25  PRINT 301,IERR
 301  FORMAT(//1X,"ERROR IN REACTION ",I3)
C
  21  CONTINUE
      PRINT *, 'TAUO2 = ', TAUO2, "WT(1)=", WT(1)


C     TRANSPORT/CHEMICAL/RHS TERMS OF ONE SPECIES: single_species.dat
      I = LCO2
      DO 918 J=1,NZ
      K = I + (J-1)*NQ

      WRITE(85,919),R(J)-R_PLT,DU(I,J),DL(I,J),DHU(I,J),DHL(I,J)
     | ,DTU(I,J),DTL(I,J),DI(I,J),FVAL(I,J)*WGT(I)/WT(J),RHS(K)
 918  CONTINUE
 919  FORMAT (1P12E10.2)
C
      GO TO 1330
      PRINT 1331
1331  FORMAT(/1X,"PERCENTAGE CHANGES OF ALL SPECIES",/)
      IROW = 10
      LR = NQ/IROW + 1
      RL = FLOAT(NQ)/IROW + 1
      DIF = RL - LR
      IF (DIF.LT.0.001) LR = LR - 1
C
1332  FORMAT(1X,1P13E9.2)
1333  FORMAT(/5X,"Z",6X,13(A8,1X))
      DO 1339 L=1,LR                 
      K1 = 1 + (L-1)*IROW                                 
      K2 = K1 + IROW - 1                                  
      IF (L.EQ.LR) K2 = NQ
      PRINT 1333, (ISPEC(K),K=K1,K2)                      
      DO 1342 I=1,NZ1,MSKIP  ! hydro56 was NZ                                  
1342  PRINT 1332, zk(I),((USOL(K,I)-USAVE(K,I))/USAVE(K,I),K=K1,K2)
      PRINT 1333, (ISPEC(K),K=K1,K2)
1339  CONTINUE

1330  CONTINUE
C
 2007 FORMAT(3x,"J", 6x,"RHO",7x,"U",8x,"T",8x,"DEN",8x,"WT"
     | ,8x,"C(H)",5x,"C(N2)"
     | ,5x,"n(H)",5x,"n_hydro(H)",3x,"n(N2)",4x,"n_hydro(N2)")
      WRITE(101,2007)
      do J=1,NZ1  ! hydro56 was NZ
      write(101,3000), J, RHO(J), U_ADV(J), T(J), DEN(J), WTP(J)
     | , USOL(LH,J), USOL(LN2,J)
     | , UNUM(LH,J), unum_h_hydro(J), UNUM(LN2,J), unum_n2_hydro(J)
     | , DI(LH,J), DI(LN2,J), EDD(J)
      enddo
      WRITE(101,2007)
C
C-TF  Escape of hydrogen: VEFF(H) = (Bi/Ha)/Nt, Nt=total number density
C-TF  These are diffusion-limited escape flux at the lower boundary
C
       J=5  ! ~100 km altitude
       BOVERH = DI(LH,J)*DEN_M(J)/HSCALE(J)
       BOVERH2 = DI(LH2,J)*DEN_M(J)/HSCALE(J)
       BOVERH2O = DI(LH2O,J)*DEN_M(J)/HSCALE(J)
       FH_dif = UNUM(LH,J)*BOVERH/DEN_M(J)*(R(J)/R(1))**2
       FH2_dif = UNUM(LH2,J)*2*BOVERH2/DEN_M(J)*(R(J)/R(1))**2
       FH2O_dif = UNUM(LH2O,J)*2*BOVERH2O/DEN_M(J)*(R(J)/R(1))**2

   63  FORMAT(/1x,"@J=",I5,"  BOVERH =",1PE10.3,"  BOVERH2 =",1PE10.3
     1 ,"  BOVERH2O =",1PE10.3,"  PHI(LH) = ",1PE10.3
     2 ,"  PHI(LH2) = ",1PE10.3, "  PHI(LH2O) = ",1PE10.3/
     3 ,"diffusion limited flux = ", 1PE10.3,/)
C
        CrossoverM = WGTM(LH) +
     | bol_k*T(NZ1)*VEFF(LH)/DI(LH,NZ1)/GZ(NZ1)
     | /(UNUM(LH,NZ1)/DEN_M(NZ1))/
     | protonM  ! hydro56 was NZ
        PRINT*, "Crossover mass =", CrossoverM
     | , bol_k*T(NZ)*VEFF(LH)/DI(LH,NZ1)/GZ(NZ1) ! hydro56 was NZ ??
C
      IF (IZPRNTFLG6.eq.1) PRINT 1729
 1729 FORMAT(/,13x,'J',4x,'f(H)',7x,'f(H2)',5x,'f(H2O)'
     | ,4x,'phi_l(H)',3x,'phi_l(H2)',2x,'phi_l(H2O)'
     | ,2x,'phi_cg(H)',2x,'phi_cg(H2)',1x,'phi_cg(H2O)'
     | ,1x,'phi_lim')
     
CZ ?? What is all this tmp_tf stuff?  
      DO J=2,NZ2  ! hydro56 was NZ1
      tmp_tf1 = DRC(J)+DRC(J-1)
      tmp_tf = WTP(J)/WGTM(LH)*(USOL(LH,J+1)-USOL(LH,J-1))/tmp_tf1 + 
     | UNUM(LH,J)/RHO(J)*(WT(J+1)-WT(J-1))/tmp_tf1
      phi_counterG_H = -(EDD(J)+DI(LH,J))*DEN_M(J)*tmp_tf
      FH_dif = DI(LH,J)*UNUM(LH,J)/HSCALE(J)*(R(J)/R(1))**2
     | *(1-WGTM(LH)/WTP(J))
C
      tmp_tf = WTP(J)/WGTM(LH2)*(USOL(LH2,J+1)-USOL(LH2,J-1))/tmp_tf1 +
     | UNUM(LH2,J)/RHO(J)*(WT(J+1)-WT(J-1))/tmp_tf1
      phi_counterG_H2 = -(EDD(J)+DI(LH2,J))*DEN_M(J)*tmp_tf*2
      FH2_dif = DI(LH2,J)*UNUM(LH2,J)/HSCALE(J)*2*(R(J)/R(1))**2
     | *(1-WGTM(LH2)/WTP(J))
C
      tmp_tf = WTP(J)/WGTM(LH2O)*(USOL(LH2O,J+1)-USOL(LH2O,J-1))
     | /tmp_tf1 + UNUM(LH2O,J)/RHO(J)*(WT(J+1)-WT(J-1))/tmp_tf1
      phi_counterG_H2O = -(EDD(J)+DI(LH2O,J))*DEN_M(J)*tmp_tf*2
      FH2O_dif = DI(LH2O,J)*UNUM(LH2O,J)/HSCALE(J)*2*(R(J)/R(1))**2
     | *(1-WGTM(LH2O)/WTP(J)) 
C      
      IF (IZPRNTFLG6.eq.1) PRINT 1728, J
     | , UNUM(LH,J)/DEN_M(J), UNUM(LH2,J)/DEN_M(J)*2, 
     |  UNUM(LH2O,J)/DEN_M(J)*2
     | , FH_dif, FH2_dif, FH2O_dif
     | , phi_counterG_H, phi_counterG_H2, phi_counterG_H2O
     | , FH_dif+FH2_dif+FH2O_dif
     | +phi_counterG_H+phi_counterG_H2+phi_counterG_H2O
     | , 1-WGTM(LH2)/WTP(J)
      ENDDO
 1728 FORMAT(1x,"dif @ J=",I5,1P14E11.2)
      IF (IZPRNTFLG6.EQ.1) PRINT 1729
C 
      x_tf = 0.
      DO I=NQ+1,NSP
      x_tf = x_tf + UNUM(I,NZ1)  ! hydro56 was NZ
      ENDDO
      PRINT*, "@ J=NZ, all short lived species VMR=", x_tf/DEN_M(NZ1)  ! hydro56 was NZ
C

      DO 6915 J=1,NZ1  ! hydro56 was NZ
CZ Write species levels to idlpass_sp.out, along with grid point J and altitude above surface in km     
      WRITE(6904,3904),J,(R(J)-R_PLT)/1E5,(USOL(I,J),I=1,NSP)  !hydro57 was 6.371E8
 3904  FORMAT (I3,1P34E10.3)

 6915  CONTINUE
CZ
      PRINT*, " "
      PRINT*,"Mass Mixing Fractions" 
      PRINT*," of Long Lived Neutral Species:"
       PRINT 6921
 6921 FORMAT(1X,"Grid",4X,"O",10X,"O2",9X,"N2",9X,"He",9X,"H",10X,
     | "H2",9X,"CO2",8X,"CO",9X,"N4S",8X,"NO",9X,"H2O",5X,"Alt (km)")
      
      NZZ1 = MIN(NZ1,MGRIDMX)  ! hydro56 was NZ
      DO 6919 J=1,NZZ1,MSKIP
         PRINT 6920,J,(USOL(I,J),I=1,11),(R(J)-R_PLT)/1E5 !hydro57 was 6.371E8
 6920  FORMAT (I3,1P12E11.2)	 
 6919 CONTINUE
       PRINT 6921
CZ
      PRINT*, " "
      PRINT*,"Particle Number Densities per cubic centimeter " 
      PRINT*,"of Ions (O+,N+,H+) and Newton Method Species (O3,HO2,OH):"
       PRINT*,"Grid   O+          N+         H+        O3         HO2
     |       OH        Alt (km)"
      
      NZZ1 = MIN(NZ1,MGRIDMX)  ! hydro56 was NZ
      DO 6916 J=1,NZZ1,MSKIP
         PRINT 6918,J,(UNUM(I,J),I=12,17),(R(J)-R_PLT)/1E5 !hydro57 was 6.371E8
 6918  FORMAT (I3,1P7E11.2)	 
 6916 CONTINUE
 
       PRINT*,"Grid   O+          N+         H+        O3         HO2
     |       OH        Alt (km)"
cz      ENDIF 
CZ
      PRINT*,""

           PRINT*, "final G = ",G
           PRINT*, "final R_PLT = ",R_PLT
           PRINT*, "final bigM = ",bigM
           PRINT*, "final FSCALE = ",FSCALE
           PRINT*, "final ALB = ",ALB
           PRINT*, "final DELZ = ",DELZ	   
           PRINT*, "final ZTROP = ",ZTROP	
           PRINT*, "final JTROP = ",JTROP	

      PRINT*," "
      print*, "Lower Boundary Conditions:"
      PRINT*,"O  O2 N2 He H  H2 CO2 CO N4S NO H2O O+ N+ H+"
      print 1710, LBOUND1

      print*, "Upper Boundary Conditions:"
      PRINT*,"O  O2 N2 He H  H2 CO2 CO N4S NO H2O O+ N+ H+"      
      print 1710, MBOUND

      print*,""
      PRINT*,"VEFF="
      PRINT*,"   O     O2     N2     He     H      H2     CO2     CO    N
     | 4S   NO     H2O    O+     N+     H+"  
      PRINT 1711,VEFF

cz      HfluxNcomp = fluxN(NZ)*WTP(NZ)*((zk(NZ)+6371)/6471)**2
      HfluxNcomp = fluxM(NZ-3)/protonM  
      UNUMH2min = 1.E20
      UNUMHmin = 1.E20
      do J = 1, NZ1  ! hydro56 was NZ
         IF (UNUM(LH2,J).LT.UNUMH2min) THEN
            UNUMH2min = UNUM(LH2,J)
         ENDIF	    
         IF (UNUM(LH,J).LT.UNUMHmin) THEN
            UNUMHmin = UNUM(LH,J)
         ENDIF      
      ENDDO

      Tmin3 = 1E4
      do J = 8, NZ1 ! hydro56 was NZ
         IF (T(J) .LT. Tmin3) THEN
            Tmin3 = T(J)
         ENDIF	    
      ENDDO

CZ Write to idlpass_misc
      write(6954,6949)U_ADV(NZ)
      write(6954,6949)PSI(NZ)
      write(6954,6949)T(NZ)  ! CZ changed back to NZ
      write(6954,6949)ujeans(NZ)
      write(6954,6949)WTP(NZ)
      write(6954,6949)fluxN(NZ)      
      write(6954,6949)HfluxNcomp
      write(6954,6949)dif_limited_flux
      write(6954,6949)UNUM(LH2,NZ1)  ! hydro56 was NZ
      write(6954,6949)UNUMH2min
      write(6954,6949)UNUM(LH,NZ1)      ! hydro56 was NZ
      write(6954,6949)UNUMHmin
      write(6954,6949)Tmin3

CZ If No Errors, set csstatus = 0 
        WRITE(6942,6943) 0.
 6943 FORMAT(f5.1)
 6949 FORMAT(1P1E9.3)

 2911 CONTINUE

CZ Moved close statement to end
      CLOSE(710)
      END PROGRAM CoupledSUA

CZZZ
C End of main program
C *******************************************************************
C
      Function thercon(T,den_o,den_o2,den_n2,den_h,den_h2,den_co2,
     | den_tot)
C   This subroutine calculates the thermal conductivity of H2 from the
C   formula given by Saxena and Saxena, J.Phys. A.: Gern. Phys. 3, 309-
C   320, 1970.
C
C      appakH2 = 20.37 + 8.2e-2*T + 3.56e-6*T*T
C      thercon = appakH2 * 4.184 * 100.
       appakH2 = 4.184 * 100.* (20.37 + 8.2e-2*T + 3.56e-6*T*T)
C
C-TF  Approximation for the thermal conductivity is given in Rees 1989 Book
C-TF  pp126:
C
      tt=T**0.69
      thercon=( ( (den_o2+den_n2)*56.+den_o*75.9+den_h*379 )*tt + 
     | den_h2*appakH2 ) / den_tot
C
      thercon=( (den_o2+den_n2)*56.+den_o*75.9)*tt/ den_tot

      GOTO 5000
C-TF  Schunk and Nagy 2000 book contains more specific constants but with 
C-TF  different indices:
C      
      thercon=1./den_tot*( den_o2 *36. *T**0.75 + 
     |                     den_o  *54. *T**0.75 +
     |                     den_co2*0.82*T**1.28 + 
     |                     den_h  *235.*T**0.75 ) 
      IF (T .gt. 150) THEN
        thercon=thercon + 1./den_tot*(den_n2*36. *T**0.75 + 
     |                                den_h2*223.*T**0.77)
      ENDIF
      IF (T .lt. 150) THEN
        thercon=thercon + 1./den_tot*(den_n2*5.63 *T**1.12 +
     |                                den_h2*103. *T**0.92)
      ENDIF
C
5000  CONTINUE
C      print*, "glbmean thercon=", thercon

      return
      end
C   End of Function thercon(T)
C *******************************************************************
C
      SUBROUTINE TRIDAG(A,B,C,R,U,M)
C      PARAMETER (M=401)
      implicit none
C
      integer :: M,J,JJ
      double precision, dimension(M) :: A,B,C,R,U
C
      double precision :: GAM(M),BET
C      
C      print*, "in TRIDAG:0"

      IF(B(1).EQ.0.) THEN
      PRINT*, "in TRIDAG.f: B(1)=ZERO!"
      PAUSE
      ENDIF
C      print*, "in TRIDAG:1"

      BET=B(1)
      U(1)=R(1)/BET
CZ            IF (IZDEBUG.EQ.1) PRINT*, "About to calc C(J-1)"      
      DO 11 J=2,M
       GAM(J)=C(J-1)/BET
       BET=B(J)-A(J)*GAM(J)

CZ      PRINT*, "About to calc C(JJ)"   
       IF(BET.EQ.0.) THEN
       PRINT*, "in TRIDAG.f: BET=ZERO! at J=",J
       DO JJ=1,J
       PRINT 1000, JJ,A(JJ),B(JJ),C(JJ),R(JJ),U(JJ)
     | ,GAM(JJ),B(J)-A(J)*GAM(J)
       ENDDO
 1000  FORMAT(1x,"A,B,C,R,U,GAM:",I3,1P10E10.2)       
C 
       STOP
C
       PAUSE
       ENDIF
C       
       U(J)=(R(J)-A(J)*U(J-1))/BET
11    CONTINUE
C
      DO 12 J=M-1,1,-1
        U(J)=U(J)-GAM(J+1)*U(J+1)
12    CONTINUE
      RETURN
      END
C   End of SUBROUTINE TRIDAG
C *******************************************************************
C
      SUBROUTINE SUBSONIC_END_OUTPUT(ujeans,psi,alam,
     |  delT,fluxN,fluxM,WTP,Tinf,fh2,phidif,i_exo,TIME,zk,
     | dif_limited_flux)
CZ This subroutine writes output file utden.out (unit 103)
CZ Added fluxM
      use parameters,only: NZM,protonM,bol_k,NSP
      use flds_atmosphere,only: Z,RHO,DEN,U_ADV,T,NZ,DI,USOL,UNUM
CZ Added EDD for utden printout
     |  ,EDD

CZ Added Flux; hydro59: changed to FLUXKZ
      USE flds_transport,ONLY:FLUXKZ

C      
      implicit none
C			
      integer :: i_exo,i
      double precision, dimension(NZM) :: ujeans,psi,alam,delT
      DOUBLE PRECISION, DIMENSION(NZM) :: fluxN,fluxM,WTP,zk,fluxH
      double precision, dimension(NZM) :: FLUXHtot,USOLHtot
      double precision :: Tinf,fh2,phidif,TIME
      REAL :: dif_limited_flux
cz      DOUBLE PRECISION, DIMENSION(NSP,NZM) :: USOL,UNUM
CZ
CZ Added calcsof H flux from Feng()
      DO i=1,NZ-1  ! hydro56 was NZ
CZ sum of MMRs times baryon particle flux p+n      
          USOLHtot(i) = (USOL(6,i) + USOL(5,i) + USOL(14,i))
          fluxH(i)=fluxN(i)*WTP(i)*USOLHtot(i)
CZ sum of H species fluxes calc by Feng; hydro59: replaced w/ FLUXKZ 	  
           FLUXHtot(i)=FLUXKZ(5,i)+2*FLUXKZ(6,i)+FLUXKZ(14,i)
CZ Added HO2 and OH	   
     |	   +FLUXKZ(16,i)+FLUXKZ(17,i)	   
      ENDDO
CZ
      write(103,224), NZ-1 ! hydro56 was NZ
 224  format(I5)
CZ Added mass flux, molecular weight in protons WT; added extra rows of data beyond exobase
CZ Added DI(LH2,i)
CZ Added H_flux
      write(103,204)
 204  format(3x,'i',7x,'Z',7x,'rho',9x,'den',8x,'u',9x,'ujeans',
     *  6x,'psi',9x,'WTP',8x,'T',9x,'delT',6x,'p+n_flux',4x,'H_flux',
     *  3x,'mass_flux',4x,'DI(H2)',6x,'EDD',6x,'MMR(H2)',4x,
     * 'MMR(H)',2x,'MMR(Htot)',2x,'HtotF')
cz     'MMR(O)',4x,'MMR(O+)')
      write(103,214)
 214  format(1x,'grid#',3x,'(km)',5x,'(g/cm3)',5x,'(cm-3)',5x,
     |  '(cm/s)',5x,'(cm/s)',5x,'mach^2',5x,'protons',6x,'(K)',8x,'(K)',
     | 5x,'(p/cm2-s)',2x,'(p/cm2-s)',2x,'(g/cm2-s)',2x,'(/cm2-s)',
     | 4x,'(/cm2-s)',4x,'ratio',6x,'ratio',5x,'ratio')
cz     'ratio',6x,'ratio')
      do 11 i=1,NZ  ! MEZ had NZ+15
      write(103,205) i,zk(i),rho(i),den(i),U_ADV(i),
     |  ujeans(i),psi(i),WTP(i),T(i),delT(i),WTP(i)*fluxN(i),fluxH(i),
     |   fluxM(i),
     |  DI(6,i),EDD(i),USOL(6,i),USOL(5,i),USOLHtot(i),FLUXHtot(i)
cz     ,USOL(1,i),USOL(12,i)
 205  format(1x,i4,1p18e11.3)
CZ USOL(5,i) = H, USOL(6,i) = H2 
CZ 
  11  continue
      write(103,204)
      write(103,207)
 207  format(1x,i2)
      write(103,218)
 218  format(3x,'Tinf',4x,'den(1)',5x,'fh2',7x,'phidif',4x,'max(psi)',
     *  2x,'time',5x,'i(max_psi)',4x,'i(exobase)')
      write(103,209), Tinf, den(1), fh2, phidif, maxval(psi),TIME,
     *  maxloc(psi),i_exo
 209  format(1x,1p6e10.3,2x,i3,10x,i3,/)
CZ
      write(103,287), dif_limited_flux
 287  FORMAT("Diffusion limited flux at 100 km = ", 1PE10.3,/)
CZ
C      
      RETURN
      END
C   End of SUBROUTINE SUBSONIC_END_OUTPUT
C *******************************************************************
C
      SUBROUTINE find_exobase(i_exo,i_exo_flg,ipsimaxgrid,u_inf_ini)
CZ2 add u_inf_ini to call statement      
      use parameters,only: siginv,NEXODELT
      use flds_atmosphere,only: DEN,NZ,HSCALE,zk
      implicit none
C
      integer :: i,i_exo,i_exo_flg,i_exo_save,ipsimaxgrid
      double precision :: prod,u_inf_ini
cz      
      INTEGER :: IZDEBUGEXO = 0 
cz
      IF(IZDEBUGEXO.EQ.1) PRINT*,"FindExobase1 i_exo =", i_exo
      i_exo_save=0
CZ Removed condition i_exo > 10      
      IF ((i_exo .LE. 600)) i_exo_save=i_exo
cz      
      IF(IZDEBUGEXO.EQ.1) PRINT*,"FindExobase1 i_exo_save =", i_exo_save
       
      i_exo_flg=0
      do 12 i=1,NZ ! hydro60 back to NZ
      prod = DEN(i)*HSCALE(i)
      if (prod .lt. 0.5*siginv) then
C      if (prod .lt. siginv) then
        i_exo_flg=1   ! found the exobase
        go to 13
      endif
C      PRINT 10, i, zk(i), DEN(i), HSCALE(i), prod, siginv
  10  FORMAT(1x,"in find_exobase:",I3,1P5E10.2)      
  12  continue
  13  continue
      i_exo=i
      
CZ Add control to hold exobase no higher than 2 grid points above critical point       
cz      IF(i_exo .gt. ipsimaxgrid+2) i_exo = ipsimaxgrid+2
CZ Add control to hold exobase at exactly 3 points above critical point
CZ2 but only for mode where we are controlling psimax 
      if (u_inf_ini .ge. 0) i_exo = ipsimaxgrid+3
      
CZ Add control on how much i_exo can change      
       IF (i_exo_save .ge. 10 .and. i_exo .ge. 10)  THEN !CZ allow large change in i_exo if starts <10
        IF(i_exo .GT. i_exo_save+NEXODELT) i_exo=i_exo_save+NEXODELT
        IF(i_exo .LT. i_exo_save-NEXODELT) i_exo=i_exo_save-NEXODELT
      ENDIF

      IF(IZDEBUGEXO.EQ.1)PRINT*,"End FindExobase i_exo =", i_exo   
C
      RETURN
      END
C   End of SUBROUTINE find_exobase
C *******************************************************************
C-TF  adjust the top boundary location of the model
C
      SUBROUTINE adjust_exobase(i_exo,rootpi,fluxN,fluxM,IZ300,IZ300A,
     |   bigM,ipsimaxgrid)
C      
      use parameters,only: NZM,bol_k,bigG,siginv,NSP,NSP4,NEXODELT
      use flds_atmosphere
      use flds_species,only: WGT, LOI, LNOI, LO2I, LHI, LNI
C     
      implicit none
C
C-TF  input parameters
      integer :: i_exo,IZ300,IZ300A,ipsimaxgrid
      double precision :: rootpi, fluxN(NZM),
     |      fluxM(NZM),bigM
CZ2     
     | ,u_inf_ini
      
C-TF  internal variables
      integer :: Nsave, i, Nm5, ii, Np1, i_exo_save
      double precision :: ve, veu, veu2, r0av, rhsu, ui2, rhsu1
     | , u02av, rhsr, prod
      double precision, dimension(NZM) :: u02,r0_cri,alam,alamp,ujeans
cz      
      INTEGER :: IZDEBUGEXO = 0 
cz
      Nsave = NZ
      IF (IZ300A.EQ.1) NZ = NZ+1 ! cz added for option to force exobase up to 300
      IF (NZ .lt. NZM) Np1 = NZ+3
      IF (NZ .eq. NZM) NZ = NZM

C-TF  If i_exo<NZ, the exobase is below the current top boundary
C-TF  move it down.
CZ Force code to raise NZ if IZ300 is not zero
CZ
      i_exo_save = i_exo
      IF(IZDEBUGEXO.EQ.1) PRINT*,"AdjustExobase1 i_exo =", i_exo   

      if((i_exo .lt. NZ).AND.(IZ300.EQ.0)) go to 24
C      
C-TF  Now the exobase is at or above the top of the grid, so need 
C-TF  to expand the grid to see if it should move up. 
C-TF  To do this, first calculate T, U_ADV, and RHO outside the 
C-TF  old exobase
C
C      write(101,1000)
      IF(IZDEBUGEXO.EQ.1) print 1000
 1000 format(/1x,'Grid moves up')
C
C   Assume that T is constant above the old exobase
      do 16 i=NZ+1,Np1  ! CZ changed back to NZ+1
C      
      T(i) = T(NZ) ! CZ changed back to NZ
      Te(i) = Te(NZ) ! CZ changed back to NZ
      Ti(i) = Ti(NZ) ! CZ changed back to NZ
      EDD(i) = EDD(NZ)  ! hydro60 back to NZ
      WT(i) = WT(NZ)   ! hydro60 back to NZ
      WTP(i) = WTP(NZ)   ! hydro60 back to NZ      

      USOL(:,i-1) = USOL(:,NZ1)**2/USOL(:,NZ-2)  ! hydro60 changed to i-1, NZ1, NZ-2
      USAVE(:,i-1) = USOL(:,i-1)

      HSCALE(i) = bol_k*T(i)/(WT(i)*GZ(i))
C      
      u02(i) = bol_k*T(i)/WT(i)
      r0_cri(i) = bigG*bigM/(2.*u02(i))
      alam(i) = R(i)/HSCALE(i)
      ve = sqrt(2.*bigG*bigM/R(i))
      veu = ve - U_ADV(i)
      veu2 = veu*veu
      alamp(i) = veu2/(2.*u02(i))
      ujeans(i) = rootpi/2. * exp(-alamp(i))*sqrt(2.*u02(i))*
     *  (1.+alamp(i))

      if(ujeans(i).ge.1.E3)THEN
       PRINT*,"In adjust_exobase, ujeans(i)=",ujeans(i)
       print*," exp(-alamp(i))=", exp(-alamp(i))
       print*,"sqrt(2.*u02(i))=",sqrt(2.*u02(i))
       print*,"(1.+alamp(i))",(1.+alamp(i))
      ENDIF

  16  continue
C
C      write(101,1001) T(NZM)
 1001 format(5x,'T(NZM) =',1pe10.3)
C
C-TF  Recalculate U_ADV  
CZ 
      do 17 i=NZ+1,Np1  ! changed back to NZ+1
      r0av = 0.5*(r0_cri(i)+r0_cri(i-1))
      rhsu = alog(T(i)/T(i-1)) - alog(WT(i)/WT(i-1)) + 2.*r0av*
     *  (1./R(i-1)-1./R(i)) -2.*alog(R(i)/R(i-1))
      ui2 = U_ADV(i)*U_ADV(i)  !   hydro56 was i-1
      rhsu1 = rhsu/(1.-ui2/u02(i-1))
cz      U_ADV(i) = U_ADV(i-1)*exp(rhsu1)
  17  continue
C
C      write(101,1002) U_ADV(NZM)
 1002 format(5x,'U_ADV(NZM) =',1pe10.3)
C
cz       PRINT*,"Np1 =",Np1
C-TF  Now solve the density equation upward from the previous upper boundary
      do 18 i=NZ+1,Np1 ! hydro59 back to NZ+1
      u02av = 0.5*(u02(i)+u02(i-1))
      r0av = 0.5*(r0_cri(i)+r0_cri(i-1))
      rhsr = - (U_ADV(i)**2-U_ADV(i-1)**2)/(2.*u02av) 
     |       + 2.*r0av*(1./R(i)-1./R(i-1)) 
     |       - alog(T(i)/T(i-1)) 
     |       + alog(WT(i-1)/WT(i-2)) ! hydro56 was i, i-1
      RHO(i) = RHO(i-1)*exp(rhsr)  ! hydro58 back to i, i-1
      DEN(i) = RHO(i)/WT(i)        ! hydro58 changed to i
      
CZ hydro59 added
      RHO_M(i-1) = SQRT(RHO(i)*RHO(i-1))
      DEN_M(i-1) = SQRT(DEN(i)*DEN(i-1))  
      
CZ       
C
C      IF (NZ .ge. 73) PRINT*, "in adjust_exobase:", i, DEN(i)
CZ hydro60 changed all i below to  i-1      
      DO ii=1,NQ
      UNUM(ii,i-1) = RHO_M(i-1)*USOL(ii,i-1)/WGT(ii)
      ENDDO
      DO ii=NQ+1,NSP
      UNUM(ii,i-1) = UNUM(ii,NZ1)*r2(NZ1)/r2(i-1) ! hydro56 was NZ
      ENDDO
      UNUM(NSP+1,i-1) = 1.
      UNUM(NSP+2,i-1) = DEN(i)
      UNUM(NSP+3,i-1) = 1.
      UNUM(NSP+4,i-1) = UNUM(LOI,i-1)+UNUM(LNOI,i-1)
     |             +UNUM(LO2I,i-1)+UNUM(LHI,i-1)+UNUM(LNI,i-1)
      UNe(i-1) = UNUM(NSP+4,i-1)

C
CZ Added fluxM 
      fluxN(i) = DEN(i)*U_ADV(i)*r2(i)/r2(1)
      fluxM(i) = RHO(i)*U_ADV(i)*r2(i)/r2(1)      
  18  continue
        
cz      print*,"in adj_exo, UNe="
cz      do i=1,NZ
cz      print*,UNe(i)
cz      ENDDO  
C
C      write(101,1003) rho(NZM)
 1003 format(5x,'rho(NZM) =',1pe10.3)
C
C   Finally, recalculate the exobase height starting from the old value
CZ Disable this for IZ300 mode
cz      IF (IZ300.EQ.0) THEN
      i_exo_save = i_exo
      
      IF(IZDEBUGEXO.EQ.1)PRINT*,"AdjustExobase2 i_exo_save =",i_exo_save      
      
      do 22 i=NZ,Np1 ! CZ changed back to NZ
      i_exo = i
      prod = DEN(i)*HSCALE(i)
      if (prod .lt. 0.5*siginv) go to 24
C      if (prod .lt. siginv) go to 24
  22  continue
cz        ENDIF
C
CZ Add control to keep exobase no higher than 2 grid points above critical point
cz      IF(i_exo .gt. ipsimaxgrid+2) i_exo = ipsimaxgrid+2
CZ Add control to hold exobase at exactly 3 points above critical point
CZ2 But only for case where we are controlling psimax 
       if (u_inf_ini .ge. 0.) i_exo = ipsimaxgrid+3

CZ Control how much exobase can move per time step
  24  CONTINUE

      IF(IZDEBUGEXO.EQ.1) then
        PRINT*,"AdjustExobase2a before NEXODELT limit, i_exo =",i_exo  
        print*,"NEXODELT =",NEXODELT  
      endif
       
CZ removed control for i_exo < 10  
       IF (i_exo_save .ge. 10 .and. i_exo .ge. 10)  THEN   !CZ allow large change in i_exo if starts <10
        IF(i_exo .GT. (i_exo_save+NEXODELT)) i_exo=i_exo_save+NEXODELT
        IF(i_exo .LT. (i_exo_save-NEXODELT)) i_exo = i_exo_save-NEXODELT
      endif
      
      IF(IZDEBUGEXO.EQ.1)PRINT*,"AdjustExobase3 i_exo_save=",i_exo_save        
      IF(IZDEBUGEXO.EQ.1) PRINT*,"AdjustExobase3 i_exo =", i_exo       
       
CZ Disable lowering atmosphere ceiling if IZ300 not equal 0
CZ Added debug message
       IF ((IZ300.EQ.0).OR.(NZ.LT.i_exo)) then
         if(IZDEBUGEXO.eq.1) print*,"NZ=",NZ,"i_exo=",i_exo       
         NZ = i_exo
         if(IZDEBUGEXO.eq.1) print*,"updated NZ=",NZ	 
       endif
      Nm5 = Nsave - 1.
C      
C-TF  set a lower limit for the atmosphere to collapse in one iteration
      if(NZ .lt. Nm5) NZ = Nm5

 2000 FORMAT(/1x,"in adjust_exobase:",3I3)
C
         if(IZDEBUGEXO.eq.1) print*,"Returning from Adjust Exobase"

      RETURN
      END
C     END of SUBROUTINE adjust_exobase      
C *******************************************************************
C
      SUBROUTINE tesolv(delTe)
C
      use flds_atmosphere,only: NZ
      implicit none
C      
      integer :: i
      double precision :: delTe(NZ) ! 
C
      DO i=1,NZ !
      delTe(i) = 0.
      ENDDO

      RETURN
      END
C     End of SUBROUTINE tesolv
C *******************************************************************
C *******************************************************************
C     
      SUBROUTINE setup_tridag(DTINV,NN,NSTEPS,A_EDDY,B_EDDY,R_PLT)
CZ MEZ modified for cs_p1 to include A_EDDY, eddy diffusion coeff.
      use parameters,only: protonM,bol_k,M_EVAPO_FLG
      use flds_atmosphere,only: R,r2,T,GZ,RHO,U_ADV,drp,drm,NZ,NZ1
CZ MEZ added HSCALE      
     | ,HSCALE 
      use flds_heating,only: ak,bk,ck,rhsk,cond,q,cp
     | ,q_cond1,q_adv1,q_cond2,q_adv2,q_cond,q_adv,q_v,q_grav
      implicit none
C      
      double precision :: DTINV
C
      integer :: i,NN,NSTEPS
      double precision :: rip,rip2,rim,rim2,appakm,appakp,drp2,drm2,dr2
     | ,denom1,denom2,d
CZ
      REAL :: A_EDDY,B_EDDY,EDDYFAC,Zcm,R_PLT

CZ      
      do i=2,NZ1 ! CZ changed back to NZ-1
      Zcm = R(i) - (R_PLT)
      rip = (R(i)+R(i+1))/2.
      rip2 = rip*rip
      rim = (R(i)+R(i-1))/2.
      rim2 = rim*rim
      appakm = cond(i-1)
      appakp = cond(i)
CZ
      EDDYFAC = A_EDDY*exp(-Zcm/(B_EDDY*HSCALE(1))) 
      appakm = appakm*(1+EDDYFAC)
      appakp = appakp*(1+EDDYFAC)
      
CZ
      drp2 = drp(i)*drp(i)
      drm2 = drm(i)*drm(i)
      dr2 = drp(i) + drm(i)
      denom1 = drp(i)*drm(i)*dr2
      denom2 = RHO(i)*cp(i)*r2(i)*dr2
C
C   Set things up in tridiagonal form for neutral temperature
C   (using inverse Euler method)
C
      q_cond1(i) = 2.*rim2*appakm/drm(i)/denom2
      q_cond2(i) = 2.*rip2*appakp/drp(i)/denom2
     
       IF(q_cond1(i).gt.1.E8)THEN
        print*," "
        PRINT*,"WARNING, q_cond1(i)=",q_cond1(i)
        print*,"i =",i
        PRINT*,"appakm=",appakm
        PRINT*,"EDDYFAC=",EDDYFAC
        print*,"1/denom2=",1/denom2
        print*,"RHO(i)=",RHO(i)
        PRINT*,"cp(i)=",cp(i)       
       ENDIF

C      
      IF (M_EVAPO_FLG .eq. 0) THEN
C
      q_adv1(i) = U_ADV(i)*drp2/denom1
      q_adv2(i) = U_ADV(i)*drm2/denom1
C
      q_grav(i) = U_ADV(i)*GZ(i)      
      ENDIF
C
      IF (M_EVAPO_FLG .eq. 1) THEN   ! Ignore adiabatic cooling related to U_ADV
      q_adv1(i) = 0.
      q_adv2(i) = 0.
      q_grav(i) = 0.
      ENDIF
C
      ak(i) =  q_adv1(i) + q_cond1(i)
      bk(i) = - ( q_adv1(i)-q_adv2(i) ) - (q_cond1(i)+q_cond2(i))
      ck(i) = -q_adv2(i) + q_cond2(i)
      d = q(i)/cp(i) - q_grav(i)/cp(i)
      rhsk(i) = ak(i)*T(i-1) + bk(i)*T(i) + ck(i)*T(i+1) + d
C
      q_cond(i) = cp(i) * ( q_cond1(i)*(T(i-1)-T(i))
     |                    + q_cond2(i)*(T(i+1)-T(i)) )
      q_adv(i)  = cp(i) * ( q_adv1(i) *(T(i-1)-T(i))
     |                    - q_adv2(i) *(T(i+1)-T(i)) )
     |          - q_grav(i)
C     
 3000 format(1x,"in setup_tridag:",i3,2x,1p26e10.2)
C
      if(rhsk(i).ge.1.e8)then
       print*," "
       print*,"WARNING, in setup_tridag, rhsk(i)=",rhsk(i)
       print*,"i=",i         
       print*,"ak(i)=",ak(i)
       print*,"bk(i)=",bk(i)
       print*,"ck(i)=",ck(i)
       print*,"d=",d
       print*,"cp(i)=",cp(i)
       PRINT*,"q(i)=",q(i)
       PRINT*,"q_grav(i)=",q_grav(i) 

       PRINT*,"q_adv1(i)=",q_adv1(i)
       PRINT*,"q_adv2(i)=",q_adv2(i)       
       PRINT*,"q_cond1(i)=",q_cond1(i)
       PRINT*,"q_cond2(i)=",q_cond2(i)
     
      endif

      enddo
       
C   Set boundary conditions
      ak(1) = 0.
      bk(1) = bk(2)
      ck(1) = 0.
      rhsk(1) = 0.
C
      ak(NZ) = - bk(NZ1)  ! CZ changed back to NZ, NZ1
      bk(NZ) = bk(NZ1)  ! CZ changed back to NZ, NZ1
      ck(NZ) = 0.  ! CZ changed back to NZ
C     
C-TF  Whether use evaporative cooling or adiabatic cooling is the same at the top boundary
      rhsk(NZ) = bk(NZ)*(T(NZ)-T(NZ1)) + q(NZ)/cp(NZ)  ! CZ changed back to NZ, NZ-1; cp changed to NZ-1
     

CZ hydro56: Check for NaN or Inf      
      if(rhsk(NZ).EQ.rhsk(NZ)*0.5+1)THEN
      PRINT*,"WARNING: rhsk(NZ) = ",rhsk(NZ)
      PRINT*,"T(NZ) =",T(NZ)
      print*,"q(NZ) =",q(NZ)
      PRINT*,"cp(NZ) =",cp(NZ-1)      
      ENDIF
      
      if(bk(NZ).EQ.bk(NZ)*0.5+1)then
      print*," "
      print*,"WARNING: bk(NZ) =",bk(NZ)
      print*,"i =",i       
      print*,"q_adv1(NZ1) =",q_adv1(NZ1)
      PRINT*,"q_adv2(NZ1) =",q_adv2(NZ1)
      PRINT*,"q_cond1(NZ1) =",q_cond1(NZ1)
      print*,"q_cond2(NZ1) =",q_cond2(NZ1)
      print*,"drm(NZ1) =",drm(NZ1)
      print*,"drp(NZ1) =",drp(NZ1)
      print*,"cp(NZ1) =",cp(NZ1)
      print*,"cond(NZ1) =",cond(NZ1)
      print*,"cond(NZ-2) =",cond(NZ-2)
      PRINT*,"rim2/drm(i)/denom2",rim2/drm(i)/denom2
      print*,"drp2/denom1=",drp2/denom1
      print*,"drm2/denom1=",drm2/denom1
     
      endif

CZ This flag was set to 0 in Feng's code
      IF (M_EVAPO_FLG .eq. 0) THEN   ! Ignore adiabatic cooling related to U_ADV
      rhsk(NZ) = rhsk(NZ) - U_ADV(NZ)*GZ(NZ)/cp(NZ)  !hydro60 back to NZ
      ENDIF
C
C   Compute the matrix (I/dt-J) and store as vectors a, b, and c
      do 2555 i=1,NZ  ! changed back to NZ
      ak(i) = - ak(i)  
      bk(i) = DTINV - bk(i)
      ck(i) = -ck(i)
 2555 continue

      RETURN
      END
C     END of SUBROUTINE setup_tridag
C *******************************************************************
C *******************************************************************
C----------------------------------------------------------------------
C
       SUBROUTINE TRIDAG_LEI(DELTA,MAX,IF,L,A,B,C,D)
C       
C-TF  FROM Jiuhou Lei in Jan. 2007
C
C Solving A SYSTEM OF LINEAR SIMULTANEOUS EQUATIONS WITH A
C TRIDIAGONAL COEFF MATRIX. THE EQNS ARE NUMBERED FROM IF TO L,
C & THEIR  SUB-DIAG. , DIAG. ,& SUPER-DIAG COEFFS. ARE STORED
C IN THE ARRAYS A, B, C. THE RIGHT HAND SIDE OF THE VECTOR IS
C STORED IN D. THE COMPUTED SOLUTION VECTOR IS STORED
C IN THE ARRAY DELTA. THIS ROUTINE COMES FROM CARNAHAN, LUTHER,
C AND WILKES, APPLIED NUMERICAL METHODS, WILEY, 1969, PAGE 446;
C
C Input:
C      MAX: The size of A,B,C,D,and DELTA;
C      IF: The index of the first data in A,B,C,D,and DELTA;
C      L: The index of the last data in A,B,C,D,and DELTA;
C      A,B,C,D: Coefficients of LINEAR SIMULTANEOUS EQUATIONS;
C
C OutPut:
C      DELTA:   SOLUTION VECTOR;
C
       IMPLICIT REAL*8 (A-H,O-Z)
       INTEGER MAX,IF,L
       DIMENSION A(1:MAX),B(1:MAX),C(1:MAX),D(1:MAX),DELTA(1:MAX)
       DIMENSION ALPHA(1:MAX),GAMMA(1:MAX),Te(1:MAX)
       INTEGER IFP1,I,K,LAST
C.....  COMPUTE INTERMEDIATE ARRAYS ALPHA & GAMMA
       ALPHA(IF)=B(IF)
       GAMMA(IF)=D(IF)/ALPHA(IF)
       IFP1=IF+1
            IF (IZDEBUG.EQ.1) PRINT*, "About to calc C(I-1)"          
       DO 1 I=IFP1,L
          ALPHA(I)=B(I)-A(I)*C(I-1)/ALPHA(I-1)
          GAMMA(I)=(D(I)-A(I)*GAMMA(I-1))/ ALPHA(I)
1      CONTINUE
C.....  COMPUTE FINAL SOLUTION VECTOR V
       DELTA(L)=GAMMA(L)
       LAST=L-IF
            IF (IZDEBUG.EQ.1) PRINT*, "About to calc C(L-K)"       
       DO 2 K=1,LAST
          I=L-K
          IF ((K.EQ.1) .AND. (IZDEBUG.EQ.1))  PRINT*,"L-K = I = ",I
          DELTA(I)=GAMMA(I)-C(I)*DELTA(I+1)/ALPHA(I)
          Te(I) = DELTA(I)**(1./3.5)
C       print 9999, I, DELTA(I), Te(I)
9999   FORMAT(1x, "in tridag_lei:", I5, 1P14E10.2)

2      CONTINUE
       RETURN
       END
C *******************************************************************
C *******************************************************************
C
      SUBROUTINE setup_tridag_lei(DTINV)
C
      use parameters,only: erg,bol_k
      use flds_atmosphere,only: R,r2,drp,drm,T,Ti,Te,UNe,DR,DRC,RHO
     | ,NZ,NZ1
      use flds_heating,only: ake,bke,cke,rhske,conde
     | ,clei,clin,clen,qecool,qse
      implicit none
C
      double precision :: DTINV
C 
      INTEGER k
      double precision :: rip, rim, temp1, qtophe, ce 
C
C-TF  EQ TO SOLVE HERE IS: 
C
! Start coefficients and rhs for trsolv:
C      do k=lev0,lev1-1
      do k=2,NZ  ! CZ   changed back to NZ
C      ce = 1.5*UNe(k)*bol_k 
      rip = (R(k)+R(k+1))/2.
      rim = (R(k)+R(k-1))/2.
      
C      p_coef(k,i) = 2./7.*sindipmag(i)/(h_mid(k,i)*dz**2) ! s1
C      temp1 = 2./7./r2(k)/dr

CZ ?? This electron/ion stuff is the only place DRC is used; also clei clin are calculated using USOL/UNUM, but go up to NZ, 
CZ      so had to add lines setting value at NZ equal to value at NZ-1 ??

C      r_coef(k,i) = p_coef(k,i)*tek0(k+1,i)/h_int(k+1,i)  ! s3
C      cke(k) = temp1*conde(k)*rip**2/drp(k)
      cke(k) = 2./7*conde(k)*rip**2/DR(k)/DRC(k)
      
C      p_coef(k,i) = p_coef(k,i)*tek0(k  ,i)/h_int(k  ,i)  ! s1
C      ake(k) = temp1*conde(k-1)*rim**2/drm(k)
      ake(k) = 2./7*conde(k-1)*rim**2/DR(k)/DRC(k-1)
      
C      q_coef(k,i) = -(p_coef(k,i)+r_coef(k,i))            ! s2
      bke(k) = -(ake(k)+cke(k)) - (clei(k)+clen(k))/(Te(k)**2.5)*r2(k)
      
C      rhs(k,i) = 0.                                       ! s4
      rhske(k) = (- qse(k) - clei(k)*Ti(k) - clen(k)*T(k) )*r2(k)
C
C-TF  by adding the following line, we proved that the solver works
C-TF  when rhske(k) has hte same sign as bke(k). Because of this problem,
C-TF  the treatment used in TIEGCM is needed here. Namely, to seperate
C-TF  te from ti and tn in qecool(k).
C
C      if (rhske(k) .gt. 0.) rhske(k) = -rhske(k)
C
C      print 9999, k, ake(k), bke(k), cke(k), rhske(k), qse(k), 
C     | clei(k)*Ti(k), clen(k)*T(k)      
C     | temp1, R(k) ,conde(k), conde(k+1), drm(k), drp(k)
9999  FORMAT(1x, "in setup_tridag_lei:", I5, 1P14E10.2)

C      enddo ! k=lev0,lev1-1
      enddo
C
      rip = (R(1)+R(2))/2.
C      temp1 = 2./7./r2(1)/dr
C
C      ake(1) = temp1*conde(1)*rip**2/drp(1)
      ake(1) = 0.
C      
C      cke(1) = temp1*conde(1)*rip**2/drp(1)
      cke(1) = 2./7.*conde(1)*rip**2/DR(1)/DRC(1)
      cke(1) = 0.
C      
C      bke(1) = -2*ake(1) - cke(1) - (clei(1)+clen(1))/(Te(1)**2.5)
      bke(1) = -(ake(1)+cke(1)) - (clei(1)+clen(1))/(Te(1)**2.5)
C      
C      rhske(1) = - 2*ake(1)*T(1)**3.5 - qse(1) 
C     |           - clei(1)*Ti(1)-clen(1)*T(1)
      rhske(1) = bke(1)*Te(1)**3.5
C      
C      ake(1) = 0.
C      
C-TF  Note that the ratio between rhske(1) and bke(1) gives about the
C-TF  right Te.
C      
      rim = (R(NZ)+R(NZ1))/2.  !   CZ changed back to NZ, NZ1
      qtophe = 3e9*erg/DR(NZ)  !  
C      qtophe = 1e9*erg
C      qtophe = 0.
C      temp1 = r2(NZ)*dr
      
C      ake(NZ)=2./7./temp1*conde(NZ1)*rim**2/drm(NZ)
      ake(NZ)=2./7.*conde(NZ1)*rim**2/DR(NZ)/DRC(NZ1)
      cke(NZ)=0.
C      
      bke(NZ)=-(ake(NZ)+cke(NZ))
     |        -( clei(NZ)+clen(NZ) )/( Te(NZ)**2.5 )*r2(NZ)
      rhske(NZ)= ( - qse(NZ)
     |             - clei(NZ)*Ti(NZ)
     |             - clen(NZ)*T(NZ)-qtophe) * r2(NZ)
C
      RETURN
      END
C *******************************************************************
C *******************************************************************
C
C-TF  NUMBER DENSITIES (UNUM)
C-TF  MASS MIXING RATIOS (USOL)

CZ hydro59: more _M fixes
CZ hydro58: fixed RHO, RHO_M, DEN_M, etc midpoint/bdry confusion

      SUBROUTINE NORM_MMR(NORM_TYPE)
C      
      use parameters,only: NSP,NQT,protonM,SMALL,NQ
     | ,MTF_COUNT_NSP,MTF_COUNT_GLBMEAN,bol_k
      use flds_atmosphere,only: RHO,DEN,RHO_M,DEN_M,USOL,USAVE,UNUM
     | ,UNe,WT,WTP,WT_M,WTP_M,R,DR
     | ,UNUM_INI,HSCALE,T,GZ,NZ,NZ_INI
      use flds_species,only: WGT,LO,LO2,LH,LOH,LH2,LCO,LNO,LN4S,LO3
     | ,LO1D,LCO2,LN2,LH2O, LOI,LNOI,LO2I,LHI,LNI
      implicit none
C      
      INTEGER :: NORM_TYPE, I,J, NTF_DEBUG=1
      double precision :: USOL_TOT,USOL_TOT2
C
C-TF  NORM_TYPE: 
C                0 -- NORMALIZE USOL AND COMPUTE DEN,UNUM,WT,WTP FROM THE NEW USOL
C                1 -- CONVERT UNUM INTO USOL
C                2 -- CONVERT USOL INTO UNUM
C
      IF (NORM_TYPE .EQ. 1) THEN 
      DO J=1,NZ-1  ! hydro56 was NZ
      DO I=1,NSP
      USOL(I,J)=UNUM(I,J)*WGT(I)/RHO_M(J)
C      
      IF (USOL(I,J) .lt. SMALL) THEN
       USOL(I,J)=SMALL
       UNUM(I,J)=SMALL*RHO_M(J)/WGT(I)
      ENDIF
C        
C      if (I .eq. LOI) then
C      PRINT 999, J, UNUM(I,J), WGT(I)/1.66e-24, RHO(J), USOL(I,J)
C      endif
 999  FORMAT(1x,"in NORM_MMR:", I4, 1P10E10.2)      
C
      ENDDO
      ENDDO
      ENDIF   ! END NORM_TYPE .EQ. 1 *********************
C
C
      IF (NORM_TYPE .EQ. 2) THEN ! CONVERT USOL INTO UNUM
      DO J=1,NZ-1  !  hydro56 was NZ
      DO I=1,NQ
C        IF (I .eq. LH) THEN
C         USOL(LO,J) = USOL(LO,J) - (USAVE(I,J) - USOL(I,J))
C         USOL(I,J) = USAVE(I,J)
C        ENDIF

CZ hydro57 if USOL<SMALL, use USOL from layer below
      IF(USOL(I,J).LE.SMALL .AND. J.GE.5) THEN
cz       print*,"WARNING, IN NORM_MMR(2), USOL(I,J) =",USOL(I,J) ! Warning added hydro59
cz       PRINT*,"I=",I," J=",J
        IF(J.GE.2.) USOL(I,J)=USOL(I,J-1)
      ENDIF
CZ hydro57
        UNUM(I,J) = RHO_M(J)*USOL(I,J)/WGT(I)

CZ	
        if(UNUM(I,J).EQ.UNUM(I,J)*0.5) THEN
          PRINT*,"Error, in NORM_MMR, UNUM(I,J)=",UNUM(I,J)
          print*,"RHO_M(J)=",RHO_M(J),"USOL(I,J)=",USOL(I,J)
          print*,"I=",I,"J=",J
        ENDIF
	
C
      ENDDO
C-TF  RESET TOTAL NUMBER DENSITY AND BALANCE ELECTRON DENSITY
       UNUM(NSP+2,J) = DEN_M(J)
       UNUM(NSP+4,J) = UNUM(LOI,J)+UNUM(LNOI,J)
     |               + UNUM(LO2I,J)+UNUM(LHI,J)+UNUM(LNI,J)
       UNe(J)=UNUM(NSP+4,J)
C
C       UNe(J)=UNe_ini(J)
      ENDDO
C
      ENDIF   ! END NORM_TYPE .EQ. 2  ****************************
C
C     NORMALIZE USOL OF LONG-LIVED NEUTRAL SPECIES - BEGIN NORM TYPE 0

CZ 1. Add up all MMRs (USOLs) and divide by total

      IF (NORM_TYPE .EQ. 0) THEN 
CZ hydro60      
      DO J=1,NZ-1
        RHO_M(J) = SQRT(RHO(J)*RHO(J+1))
      ENDDO
      
      
      DO 2000 J=1,NZ-1   !  hydro56 was NZ
      DEN_M(J)=0.
      USOL_TOT=0.
      USOL_TOT2=1-( USOL(LOI,J)+USOL(LNI,J)+USOL(LHI,J) )
C
      DO 3001 I=1,NQ
      USOL_TOT=USOL_TOT+USOL(I,J)
 3001 CONTINUE
C
CZ hydro61 cs128 put this back in
      IF (MTF_COUNT_NSP .eq. 0) GO TO 1991  ! Skip short lived species in normalization
C       
      DO I=NQ+1,NQT
      USOL_TOT=USOL_TOT+USOL(I,J)
      ENDDO
C      
      DO 2002 I=NQT+1,NSP
      USOL(I,J)=UNUM(I,J)*WGT(I)/RHO_M(J) ! hydro58 _M 
      USOL_TOT=USOL_TOT+USOL(I,J)
 2002 CONTINUE
C
 1991 CONTINUE
      DO 2003 I=1,NQ
C
CZ hydro61 cs128 put this back in
      IF (NZ_INI .eq. 57) THEN
        IF ( (J .le. 13) .AND. (I .eq. LOI) ) GO TO 1700
      ENDIF
C
      USOL(I,J)=USOL(I,J)/USOL_TOT       ! NORMALIZE ALL MASS MIXING RATIOS
      UNUM(I,J)=USOL(I,J)*RHO_M(J)/WGT(I)  !cz Find new UNUMs from USOLs and RHO_M
C      
C-TF  DO NOT CHANGE USOL and UNUM OF O+ and N+ at lower thermosphere where 
C-TF  their concentrations are low.  CZ removed this
 1700 CONTINUE   
C
      DEN_M(J)=DEN_M(J)+UNUM(I,J)              ! COUNT ALL SPECIES INTO DEN hydro58: DEN_M

      IF (UNUM(I,J) .lt. 0) GO TO 1980
C      
C      PRINT 4112, I, J, DEN(J), UNUM(I,J)
4112  FORMAT ("RENORM:",2I4,1P10E10.2)
C
      USOL_TOT2=USOL_TOT2+USOL(I,J)
 2003 CONTINUE
C
      IF ( (MTF_COUNT_GLBMEAN .eq. 1) .OR. (MTF_COUNT_NSP .eq. 0) ) 
     | GO TO 1992  ! skip short lived species in normalization
C
      DO I=NQ+1,NSP
      USOL(I,J)=USOL(I,J)/USOL_TOT
      UNUM(I,J)=USOL(I,J)*RHO_M(J)/WGT(I)
      DEN_M(J)=DEN_M(J)+UNUM(I,J)
      IF (UNUM(I,J) .lt. 0) GO TO 1980
      USOL_TOT2=USOL_TOT2+USOL(I,J)
      ENDDO
C
 1992 CONTINUE
C
CZ hydro60 now convert to WT_M per JFK b/c WT varies much slower than DEN
      WT_M(J) = RHO_M(J)/DEN_M(J)
      WTP_M(J) = WT_M(J)/protonM   

      IF(J .GT. 1) THEN
       WT(J) = (DR(J)*WT_M(J-1)+DR(J-1)*WT_M(J))/(DR(J)+DR(J+1))
      ELSE
       WT(1) = WT_M(1)
      ENDIF

      WTP(J) = WT(J)/protonM

      DEN(J) = RHO(J)/WT(J)

        IF(DEN(J).EQ.DEN(J)*0.5)THEN
          PRINT*,"WARNING in NORM_MMR(0), DEN(J) =",DEN(J),"J=",J
          PRINT*,"RHO(J) =",RHO(J),"WT(J) =",WT(J)
        ENDIF
C
      HSCALE(J) = bol_k/(WT(J)*GZ(J))*T(J)
C
C-TF  RESET TOTAL NUMBER DENSITY AND BALANCE ELECTRON DENSITY
      UNUM(NSP+2,J) = DEN_M(J)
      UNUM(NSP+4,J) = UNUM(LOI,J)+UNUM(LNOI,J)
     |              + UNUM(LO2I,J)+UNUM(LHI,J)+UNUM(LNI,J)
C
      UNe(J)=UNUM(NSP+4,J)
C
 2000 CONTINUE
C
cz      DEN(1)=RHO(1)/WT(1) ! RHO(1) and WT(1) should be constants under most B.C.

      WT(NZ) = WT_M(NZ-1)
      DEN(NZ) = RHO(NZ)/WT(NZ-1)


      ENDIF   ! END NORM_TYPE .EQ. 0

      RETURN
C
 1980 CONTINUE  ! Error handling
      PRINT*, "NEGATIVE UNUM at I, J", I, J
      PRINT*, "Stop in NORM_MMR"
      STOP
C			      
      END
C   End of SUBROUTINE NORM_MMR
C *******************************************************************
 
