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