Helplist

From krc
Jump to: navigation, search

Contents

INTRODUCTION

This document is intended to help the expert user of KRC set up an input file that addresses their goals and will generate the kind of output they desire. It is assumed that the user is familiar with the KRC journal article: H.H. Kieffer, Thermal model for analysis of Mars infrared mapping, J. Geophys. Res. Planets, (2012)

The evolution of KRC code is contained in: evolve.txt

A crude diagram of the call architecture is in: flow.txt

METHOD

The program is designed to compute surface and subsurface temperatures for a global set of latitudes at a full set of seasons, with enough depth to capture the annual thermal wave, and to compute seasonal condensation mass. For historic reasons, the code has substantial optimization. There are generalities that allow this code set to be used for any solid body with any spin vector, in any orbit (around any star); this is also the source of some of the complexity.

Method is explicit forward finite differences with exponentially increasing layer thickness and binary time increase with depths where allowed by stability. Depth parameter is scaled to the diurnal thermal skin depth. Initially starts at 18 hours with the mean temperature of a perfect conductor. Second degree perturbation is applied at the end (midnight) of the (third) day; this jumps the mean temperature of all layers and the lower boundary to equal the mean surface temperature.

Boundary condition treatment
Perturbation solution of quartic equation at surface for each iteration; temperature gradient assumed uniform in top interval.
Lower boundary may be insulating or constant-temperature.
Atmospheric Radiation
KRC uses a one-layer atmosphere that is grey in both the solar and infra-red regions. parametric atmosphere. The default atmospheric parameters are based on estimates of Mars’ gas and aerosol properties.
Delta-Eddington model for insolation; direct onto sloped surface and diffuse, with possible twilight extension. Atmosphere temperature based on Delta-Eddington solar absorption and IR opacity.

Keplerian orbital motion; seasons are at uniform increments of time. Mean orbital elements are pre-calculated for any epoch (all planets and several comets) by the PORB code set.

Units are SI, except for use of days for orbital motion and rotation period.

Options
Different Physical properties below a set layer (IC).
Regional slope.
Three ways to handle seasonal global pressure variation.
Atmosphere condensation
Global integral of CO2 frost-gas budget can control surface pressure.
Allows different surface elevation for each latitude zone.
Zonal frost saturation temperature tracks local surface pressure.
Option for cap albedo to depend upon mean daily insolation.

Convergence Notes

Convergence prediction routine can’t jump more than one time constant (TAU=X**2/2) <math> \tau =x^2/2 </math> for the total thickness. Therefore, if X(N1) is small, make DDT smaller than usual. If DELJUL is much smaller than (X(N1))**2/2 <math> X_{N1}^2/2 </math>, then DDT can be as large as 0.3. Otherwise DDT must be about 0 for the prediction routine to work well (it assumes the 3rd derivative to be 0).

INPUT FILE

All parameters for KRC are set by a formatted text file. An example is master.inp, which has default values for a 19 latitude set for a run of three martian years, with the last output to disk. Parameter values are listed below their titles, which are in many cases identical to the code name, and last charater of the title is above the last location in the field. Thus, integer values MUST be aligned. Titles with a leading "[" indicate that the value is not used. The recommended procedure is to copy master.inp and edit only the values you wish to change. The number of lines of Latitudes and Elevations must match the value of N4, e.g., 2 lines for N4=11:20, entries beyond the N4 position may be left blank or contain the end of the line. The 13 lines following Elevations are a geometry matrix for Mars orientation and orbit in 2005, and should not be touched; they can be replaced by running PORBMN carefully.

The title lines are skipped, so that you may put comments there carefully.

The first input line is always KOLD,KEEP (I*), which sets file usage; these are described in Disk Binary Files section. If KOLD=0, then a full set of input values is read. If and only if there is a third non-zero integer, then will read next card as 6 debug flags, IDB1 to IBD6, which are normally zero. See Debug Options. The next (normally second) line is free text where you can outline the purpose of your run

Change lines may follow immediately after the geometry matrix (see Parameter Changes). The end of definition of a "case" is indicated by a "0/" line. Two successive "0/" lines ends the run.

Items with numbers inset 2 spaces below are computed, not input. The source code for ’krccom.inc’ indicates which subroutine sets many of the parameters; as the routine name in lowercase just below the parameter name.

 - - - - - - - - - - - - - - - - - - - -
Type 4  Title (20A4) 80 characters of anything to appear at top of each page.

Type 1  Real parameters     (8F10.2)    ================================
   Surface Properties
1   ALB Surface albedo
2   EMIS    Surface emissivity
3   SKRC    Surface thermal inertia [J m^-2 s^-1/2 K^-1] { cal cm * 4.184e4}
4   COND2   Lower material conductivity (IC>0)
5   DENS2   Lower material density (IC>0)
6   PERIOD  Length of solar day in days (of 86400 seconds)
7   SPHT    Surface specific heat [J/(kg K)]  {cal/(g K) * 4184.}
8   DENS    Surface density [kg/m^3] {g/cubic cm. *10}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   Atmospheric Properties
9   CABR    IR opacity of dust-free atmosphere of PTOTAL surface pressure
10  AMW Molecular weight of the atmosphere
11  [ABRPHA UNUSED [Phase of ABRAMP, degrees relative to midnight] 
12  PTOTAL  Global annual mean surface pressure at 0 elev., Pascal[=.01mb]
                If KPREF=2, global average of atmosphere plus cap system.
13  FANON   Mass-fraction of mean atmosphere that is non-condensing
14  TATM    Atm temp for scale-height calculations
-   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   
15  TDEEP   Fixed bottom temperature. Used if IB>=1.
16  SPHT2   Lower material specific heat (IC>0)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  Dust  & Slope Properties
17  TAUD    Mean visible opacity of dust, solar wavelengths
18  DUSTA   Single scattering albedo of dust
19  TAURAT  Ratio of thermal to visible opacity of dust
20  TWILI   Twilight extension angle [deg]
21  ARC2    Henyey-Greenstein asymmetry factor
  moon        = eclipse start time in local Hours
22  [ARC3   NOT USED  coeff. for planetary heating 
  moon        = eclipse duration in seconds 0=no eclipse
23  SLOPE   Ground slope, degrees dip. Only pit may slope beyond pole.
24  SLOAZI  Slope azimuth, degrees east from north. <-360 is a pit
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    Frost Properties
25  TFROST  Minimum Frost saturation temperature
          may be overridden by local saturation temperature (LVFT)
26  CFROST  Frost latent heat [J/kg] {cal/gm*4184.  [ Not used if
27  AFROST  Frost albedo, may be overridden (LVFA)  [ TFROST never
28  FEMIS   Frost emissivity            [ reached
29  AF1 constant term in linear relation of albedo to solar flux
30  AF2 linear term in relation of albedo to solar flux units=1/flux
           Afrost = AF1 + AF2 * <cos incidence> SOLCON / DAU^2
31  FROEXT  Frost required for unity scattering attenuation coeff. [kg/m^2]
                the greater of this and 0.01 is always used.
32  fd32    UNUSED
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   Thermal Solution Parameters
33  RLAY    Layer thickness ratio
34  FLAY    First layer thickness (in skin depths)
35  CONVF   Safety factor for classical numerical convergence
         0 for no binary time division of lower layers
         >0.8 for binary time division. Larger is more conservative
36  DEPTH   Total model depth (scaled) (overrides FLAY if not 0.)
37  DRSET   Perturbation factor in jump convergence. If = 0., then
          all layers reset to same average as surface layer. Else,
                  does quadratic curve between surface and bottom averages
38  DDT Convergence limit of temperature RMS 2nd differences
39  GGT Surface boundary condition iteration test on temperature
40  DTMAX   Convergence test: RMS layer T changes in a day
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   Orbit Geometry & Constants
41  DJUL    Starting Julian date of run -2440000(N5>0)
42  DELJUL  Increment between seasons in Julian days (if N5>1)
43  SDEC    Solar declination in degrees. (if Not LPROB)
44  DAU Distance from Sun in astronomical units (if Not LPROB)
45  SUBS    Aerocentric longitude of Sun, in degrees. For printout 
          only. Computed from date unless N5=0(for printout only)
46  SOLCON  Solar constant Applied Optics 1977 v.16, p.2693: 1367.9 W/m^2
                    1366.2 Based on figure in Frohlich, Observations of 
            irradiance variations, Space Sci. Rev.,94,15-24,2000
47  GRAV    Surface gravity.  MKS-units
48  AtmCp   Specific heat at constant pressure of the atmosphere [J/kg/K]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   Temperature dependent conductivity.  Ignored unless LKOFT set.
49      ConUp0  Constant coef for upper material 
50      ConUp1  Linear   in k=c0+c1x+c2x^2+c3x^3 where x=(T-220)*0.01
51      ConUp2  Quadratic    " 
52      ConUp3  Cubic coeff. "
53      ConLo0  Constant coef for lower material 
54      ConLo1  Linear      as for ConUp above
55      ConLo2  Quadratic    "
56      ConLo3  Cubic coeff. "
   Temperature dependent specific heat.  Ignored unless LKOFT set.
57      SphUp0  Constant coef for upper material 
58      SphUp1  Linear   in k=c0+c1x+c2x^2+c3x^3 where x=(T-220)*0.01
59      SphUp2  Quadratic    " 
60      SphUp3  Cubic coeff. "
61      SphLo0  Constant coef for lower material 
62      SphLo1  Linear      as for   SphUp above
63      SphLo2  Quadratic    "
64      SphLo3  Cubic coeff. "
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
COMPUTED REAL*4 VALUE   
  65    HUGE  = 3.3E38   nearly largest  REAL*4 value
  66    TINY  = 2.0E-38  nearly smallest REAL*4 value
  67    EXPMIN = 86.80   neg exponent that would almost cause underflow
  68    FSPARE  Spare
  69    FLOST   Atm frost 'lost' in the atm. in last day at current lat./season
  70    RGAS  = 8.3145   ideal gas constant  (MKS=J/mol/K)
  71    TATMIN  Atmosphere saturation temperature
  72    PRES    Local surface pressure at current season
  73    OPACITY Solar opacity for current elevation and season
  74    TAUIR   current thermal opacity at the zenith
  75    TAUEFF  effective current thermal opacity 
  76    TATMJ   One-layer atmosphere temperature
  77    SKYFAC  fraction of upper hemisphere that is sky
  78    TFNOW   frost condensation temperature at current latitude
  79    AFNOW   frost albedo  at current latitude
  80    PZREF   Current surface pressure at 0 elevation, [Pascal]
  81    SUMF    Global average columnar mass of frost [MKS]
  82    TEQUIL  Equilibrium temperature (no diurnal variation)
  83    TBLOW   Numerical limit (Blowup) temperature
  84    HOURO   Output Hour requested for "one-point" model
  85    SCALEH  Atmospheric scale height
  86    BETA    Atmospheric IR absorption
  87    DJU5    Current Julian date (offset 2440000 ala PORB convention)
  88    DAM Half length of daylight in degrees
  89    EFROST  Frost on the ground at current latitude [kg/m^2] {g/cm^2 * 10.} 
  90    DLAT    Current latitude
  91    COND    Top material Thermal conductivity (for printout only)
  92    DIFFU   Top material Thermal diffusivity (for printout only)
  93    SCALE   Top material Diurnal skin depth (for printout only)
  94    PIVAL   pi
  95    SIGSB   Stephan-Boltzman constant (set in KRC)
  96    RADC    Degrees/radian

Type 2  Integer Parameters  (8I10) ====================================

1   N1  # layers (including fake first layer) (lim MAXN1)
2   N2  # 'times' per day (lim MAXN2). Must be an even number, 
          should be a multiple of N24 and NMHA.
3   N3  Maximum # days to iterate for solution (lim MAXN3)
          98sep03 This can be 1, but then must use  DELJUL ~= PERIOD
        If N3 lt 3, first day starts on midnight. else at 18H 
4   N4  # latitudes (lim MAXN4=19). Global integrations done for N4>8
5   N5  # 'seasons' total for this run. If 0, then DAU and SDEC will be 
        used as entered for a single season.
6   N24 # 'hours' per day stored, should be divisior of  N2 (lim MAXNH)
7   IB  Bottom control: 0=insulating, 1=constant temperature 
          2=start all layers =TDEEP & constant temperature 
8   IC  First layer (remember that 1 is air) of changed properties. 
          if 3 to N1-2.   > N1-2 (e.g., 999) =homogeneous
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
9   NRSET   # days before reset of lower layers first season;  >N3=no reset
10  NMHA    # 'hour angles' per day for printout (no limit)
11  NRUN    Run #; appears in some printout
12  JDISK   Season count that disk output is to begin. 0=none
13  IDOWN   Season at which to read change cards
14  I14 Index in FD of flexible print
15  I15  ""
16  KPREF   Mean global pressure control. 0=constant
          1= follows Viking Lander curve  2=reduced by global frost, but
              then N4 must be >8, and latitudes must be monotonic increasing
              and must include both polar regions (no warning for your failure)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
17  K4OUT   Disk output control: See details in DISK BINARY FILES section
        Three modes of direct access Fortran files;  one case per file.
              -=KRCCOM(once), then TSF & TPF;
              0=KRCCOM,LATCOM each season
           1:49=KRCCOM,DAYCOM for the last latitude; each season
         Modes of bin5 file for multiple cases
           51=(Hours, 2 min/max,  lat, seasons, cases)
           52=(hours,  7 items, lat, seasons, cases)
           54=[many seasons, 5 items,lats, cases]
           55=[many seasons,9 items, cases]
           56=[packed T hour and depth, latitude,season,case]
18  JBARE   J5 season count at end of which to set frost amount to 0. 0=never
19  NMOD    Spacing of season for notification. minimum of 1
20      IDISK2  Last season to disk for which to print notice
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
COMPUTED I*4 VALUES
  21    KOLD    Season index for reading starting conditions
  22    KVALB   Flag: to use seasonal surface albedo ALB
  23    KVTAU   Flag: 1:TAUD=SEASTAU(SUBS)  2:CLIMTAU opacities for dust and ice
  24    ID24(4) spare
  28    NFD Number of real items read in
  39    NID Number of integer items read in
  30    NLD Number of logical items read in
  31    N1M1    Temperature vrs depth printout limit (N1-1)
  32    NLW Temperature vrs depth printout increment
  33    JJO Index of starting time of first day
  34    KKK Total # separately timed layers
  35    N1PIB   N1+IB Used to control reset of lowest layer
  36    NCASE   Count of input parameter sets in one run
  37    J2  Index of current time of day
  38    J3  Index of current day of iteration
  39    J4  Index of current latitude
  40    J5  Index of current "season"

Type 3  Logical Parameters  (10L7) ====================================

1   LP1 Print program description. TPRINT(1) 
2   LP2 Print all parameters and change cards (2)
3   LP3 Print hourly conditions on last day (3)
4   LP4 Print daily convergence summary (4)
5   LP5 Print latitude summary (5)
6   LP6 Print TMIN and TMAX versus latitude and layer (6)
7   LPGLOB  Print global parameters each season
8   LVFA    Use variable frost albedo. Uses AF1 & AF2 (real # 29,30)
9   LVFT    Use variable frost temperatures
10  LKOFT   Use temperature-dependent conductivity and specific heat
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
11  LPORB   Call PORB1 just after full input set
12  LKEY    Read change item from terminal after main input set
13  LSC Read change cards from input file at start of each season
14  LNOTIF  spare
15  LOCAL   Use each layer for scaling depth
16  LD16    Print hourly table to fort.76 [TLATS] 
17  LPTAVE  Print <T>-<TSUR> at midnight for each layer [TDAY]
18  LD18    spare
19  LD19    Output to fort.79 [TLATS] insolation and atm.rad. arrays 
20  LONE    (Computed) Set TRUE if KRC is in the "one-point" mode
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
followed in 'krccom' by: 
[real*4] TITLE(20)  80-character title
[real*4] DAYTIM(5)  20-character run  date and time
================================================================

Latitude(s) (10F7.2) N4 latitudes in degrees, no internal separations.
Latitudes to be in order; south to north. [[If last latitude is
.LE. 0, will assume symmetric results for global integrations]]

Elevation(s) (10F7.2) N4 values in Km corresponding to latitudes

Orbital Parameters (LPORB=T) Format identical to that produced by PORB
program set ASCII file output. So these can be directly pasted with an
editor. see PORBCM.INC
 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

PARAMETER CHANGES [pc]

Fortran List Directed. Change the values in KRCCOM White-separated, a "/" terminates the read and leaves remaining values unchanged. The 4 items are:

Integer
Integer Numeric_value
’Text’
/ Comment
    1: Type (integer) see table below
    2: Index in array (integer), as listed in table above
    3: New value, numeric, will read as real and convert. 0.=false.
    4: Reason, text string within single quotes
   [after a / (forward slash) nothing is read, so you can use for comments]

The print file will list each change as read, followed by the title of the changed item. It is a good idea to look at this print to be sure you changed what you intended.

 Type     Meaning                                                        Valid Index

   0   End of Current Changes                                                 any
   1   Real Parameter                                                         1:NFDR
   2   Integer Parameter                                                      1:NIDR
   3   Logical Parameter                                                      1:NLDR
   4   New Latitude Card(s) Follow                                            any
   5   New Elevation Card(s) Follow                                           any
   6   New Orbital Parm Cards Follow (LPORB Must be True)                     any
   7   Text becomes new Title                                                 any
   8   Text becomes new disk or season-variation file name
         if index=22, call SEASALB to read variable ALBEDO
         if index=23, call SEASTAU to read variable TAUD
         if index=24, call CLIMTAU to read Mars climate
   9   Complete new set of input follows                                      any
  10   Text becomes new One-Point input file name
  11   This is a set of parameters for "one-point" model 
          For this type, 9 values must appear in a rigid format
  12   Set of 2*4 coefficents for T-dep. conductivity.  List-directed IO
  13   Set of 2*4 coefficents for T-dep. specific heat. List-directed IO 

For 12 and 13, 8 white-space-separated coefficients must follow after the type on the same line, with no intervening index or text.

To start variable albedo, use input card: 
  8 22 0 'AlbedoFileName' / Variable albedo text file name
Can revert to constant albedo by hokey technique of using a bad name. E.g.,
  8 22 0 'badName' / turn variable albedo off
Files of text table of value versus season will be read at the start of a
run. These will apply to ALL latitudes. See example  valb1.tab 
Variable Tau done the same way, with 22 being replaced with 23
 CLIMTAU files have dust and ice opacity over season and latitude. Uses BINF5 to 
read binary array (72 seasons, 36 latitudes, 2=dust/ice) of opacities

CONTENTS OF COMOMS

COMMON /KRCCOM/ Input and transfer variables. See krccom.inc
COMMON /DAYCOM/ Layer and time-of-day items. See daycom.inc
COMMON /FILCOM/ File names. See filcom.inc
COMMON /HATCOM/ Store post-2003 items. See hatcom.inc
COMMON /LATCOM/ Latitude-dependent items. See latcom.inc
COMMON /PORBCM/ PORB system geometry matrix. See probcm.inc
COMMON /UNITS/ Logical units for I/O and errors. See units.inc

Because the binding routines to IDL are intolerant of any errors, changes to KRCCOM, DAYCOM and LATCOM are avoid if possible. Rather, in 2004July HATCOM was added as a "catch-all" for any new items.

A listing of all Fortran commons can be generated by these Linux commands:
cd /home/hkieffer/krc/src [replace top part of path with local installation]
rm allinc.txt
cat krccom.inc latcom.inc daycom.inc hatcom.inc filcom.inc units.inc porbcm.inc <math>></math> allinc.txt

ERROR RETURNS

 "Parameter error in TDAY(1)" : Convergence factor < .8 classic. 
        Instability anticipated.  
 "UNSTABLE; Layer..... TDAY(1): 

DRSET: 0=>     Reset by delta_average_T for each layer:
                 else: reset by {linear + DRSET*quadratic}*{<surf>-<botm>}
TDAY: LRESET   Reset midnight T's for all but top layer.
      LDAY     Last day computations

DISK BINARY FILES [dbf]

The routine TDISK is used to read or write direct-access binary files or bin5 files. The first season to write is specified by JDISK, all following seasons will go to the same file. For direct-access files, each file record consists of KRCCOM plus LATCOM or KRCCOM plus DAYCOM.

Disk output is largely controlled by the KRC and TSEAS routines.

Items which control file I/O

KOLD & KEEP on first input line
  KOLD: 0= input card set follows;  else=disk record number to start from,
         then will read any change cards.
        If LPORB in old file was True, then there must be a PORB card set 
          as the set of lines following the KEEP,KOLD line
  KEEP: 0= close disk file after reading seasonal record KOLD;
       >0= value of JJJJJ at which to start saving seasons in same disk 
        file [overrides JDISK].
  To start from a prior seasonal run, need to determine the record 
  corresponding to the desired season;
         KOLD=J5_target - JDISK(old) ; >0
         set KEEP=1, change card J5=number of new seasons, set K4OUT.

JDISK sets the first season to save results

N5    sets the last season to run

K4OUT sets the record content:
 -      Will output first record of KRCCOM,ALAT,ELEV, then records of TSF & TPF
 0      Will output records of KRCCOM+LATCOM. Usual for large data-base.
 +n<=50 Will output records of KRCCOM+DAYCOM for the last computed latitude.

<math>></math>50 Will write custom bin5 file at the end of a run, with dimensionality from 3 to 5 (more possible). All 5x outputs allow multiple cases, each with a "prefix" for each case consisting of 4 size integers (converted to Float) followed by KRCCOM; after this may come vectors of parameters versus season. The next-to-last dimension is increased to allow room for the prefix to be embedded in the bin5 array. KRC input items that would change any of the bin5 dimensions are not allowed to change between cases. Each dimension is adjusted to the necessary size. Each case has the same structure; this simplifies coding although some items are then present redundantly. The number of cases allowed is set by the size of case one, and printed as MASE at the end of the first case in the print output. Cases beyond the maximum that can be stored will be executed, but not saved.

The first 4 words of the prefix, and of thus of the bin5 array, are:

(1)=FLOAT(NWKRC)   ! Number of words in KRCCOM
(2)=FLOAT(IDX)     ! 1-based index of dimension with extra values
(3)=FLOAT(NDX)     ! Number of those extra
(4)=FLOAT(NSOUT)   ! [Available of other use]

    51=(N24 hours, 2: TSF TPF, N4 lats, NDX+ seasons, cases)
The prefix section contains: sub_array(seasons,5)(0-based index)
  0)=DJU5  1)=SUBS  2)=PZREF  3)=TAUD  4)=SUMF

    52=(N24 hours, 7 items, N4 lats, NDX+ seasons, cases)
The 7 items are:  1)=TSF  2)=TPF  3)=TAF  4)=DOWNVIS  5)=DOWNIR
6) packed with [NDJ4,DTM4,TTA4,      followed by TIN(2+
7) packed with [FROST4,AFRO4,HEATMM, followed by TAX(2+
  The number of layers for TIN and TAX is the smaller of: the number computed 
and that fit here.
The prefix is identical to Type 51

    54= (seasons, 5 items, NDX +nlat, cases)
        Items are (0-based index): 
        0= TSF=surface temperature at 1 am, 1= TSF at 13 hours,
        2= HEATMM=heat flow, 3= FROST4=frost amount, 
        4= TTB4 = predicted mean bottom temperature
        The prefix contains DJU5 

    55= (seasons,NDX+ items,cases).  For seasonal studies at one latitude
        ITEMS intended to be recoded as needed. Initial version is 9 items:
        [Tsur@ 1am,3am,1pm, spare, Tplan @1am,1pm, Surface heat flow,
        frost budget, T_bottom]
        The prefix contains DJU5        
         Can hold very large number of seasons and cases. 
        THIS MODE DOES NOT SUPPORT CONTINUATION RUNS

    56= [vectors&items, latitudes, NDX+ seasons, cases]
The first dimension is: TSF for all hours, TPF at all hours, 
 T4 for all layers at midnight, then FROST4,HEATMM,TTA4
The prefix is identical to Type 51

Once a disk file is opened, any records written will go into that file until a new filename is specified (Type 8 Change line), which closes the current file. It is best to ensure that output file does not already exist. If the file already exists, new output may be written in same area, even if larger than needed.

HANDY THINGS

The first "hour" in printout and output arrays is 1/24 (strictly, 1/N24) of a sol after midnight. E.g., the last time is midnight, not the first.

Atmospheric scale height, SCALEH, depends upon physical constants GRAV [input] and TATMAVE which (2007nov) is TATM [input] for the first season and thereafter the diurnal average of the prior season.

To run and save various cases for a single season, set N5 and JDISK to 1.

To extract a detailed day by saving DAYCOM to disk, set JDISK=N5, set a new file name, and set K4OUT to desired latitude index (normally 1).

To run continuously with output every K ((1-3) days, set DELJUL=K*PERIOD. This will force prediction terms to near 0.

setting N3=1 will turn off all prediction.
set GGT large (to avoid iteration for convergence)
set NRSET=999 (to avoid reset of layers)

To continue run with new parameters (e.g., DELJUL)

3 21 1 ’flag set to continue’

Note: changing DELJUL will cause reset of DJUL
Must increase the value of N5: e.g., 2 5 [bigger] ’Increase stopping season’ Reset will not occur because J5 continues incrementing

ASCII Output Files

krc.prt  General results. Stuff output is controlled by LP1:6 & LPGLOB

fort.76
tlats.f: mimic Mike Mellon ASCII files
        if (ld16) then
          write(76,761)subs,dlat,alb,skrc,taud,pres
 761      format(/,'      Ls      Lt       A       I    TauD       P'
 762      format(f7.2,f9.3,f8.3,f9.3)
            write(76,762)qh,tsfh(i),adgr(j),qs

          do i=1,n24
            j=(i*n2)/n24
            qh=i*qhs
            qs=(1.-alb)*asol(j) ! absorbed insolation
            write(76,762)qh,tsfh(i),adgr(j),qs
          enddo

fort.78
tlats.f:  for average and maximum:
        if (ld18) write(78,*)cosi_(i), t_(i),ADGs(i),ADGP(i)
        if (ld18) write(78,*)j5,j4,sol,ave_a,adgir,c52,beta

fort.79
tlats.f:   for each time-step
       if (ld19) write(79,*)adgr(jj),qa,direct,diffuse
   col 1 = downgoing thermal radiation
   col 2 = total insolation reaching surface
   col 3 = direct  fraction of insolation
   col 4 = diffuse fraction of insolation

To run two material types

Set IC to the first layer to have the lower material properties ( >= 3)
Set COND2 to the lower material conductivity
Set DENS2 to the lower material density
Set SPHT2 to the lower material specific heat
If LOCAL is False, then initial setting of all layer thicknesses is based upon the scale of the upper material; if it is set True, the thickness of the lower layers is set by their scale.
TDAY no longer allows unstable (thin) layers, and will increase the thickness of the layer IC to satisfy the convergence safety factor FCONV if needed. However, the code to check on convergence was retained.

Setting temperature-dependant properties

Basic Flag is L10=LKOFT . If this is true, then the 8 input parmeters ConUp0 to ConLo3 must be set to yield thermal conductivity as a function of temperature for the upper and lower materials. <math> k=c0 +c_1x + c_2x^2 +c_3x^3 </math> where <math>x=(T-200.)*0.01</math>

Correspondingly, the 8 input parmeters SphUp0 to SphLo3 must be set for specific heat

One way to generate the coefficients is to run for each of the upper and lower materials the IDL procedure KOFTOP, which can call all of the temperature-dependant routines. KOFTOP allows change of its parameters, including grain radius and pressure, and will print the required parameters ready for input to KRC.

Below are sample coefficients for thermal conductivity based on Sylvain Piqueux’s numerical model for un-cemented soils; the fit error is <math><0.1</math>% over 120-320K. Left column is grain radius in micrometers, then the four normalized coefficients ready for inclusion in a KRC input file, followed by the thermal inertia at 220K for nominal density and specific heat.

 R(mu)         c0        c1        c2        c3      Iner  
    10.     0.008274  0.000735 -0.000376  0.000148    89.8 
    20.     0.012379  0.001280 -0.000629  0.000250   109.9 
    50.     0.021485  0.002647 -0.001201  0.000483   144.7 
   100.     0.032051  0.004528 -0.001874  0.000761   176.8 
   200.     0.046023  0.007569 -0.002743  0.001129   211.8 
   500.     0.068387  0.014075 -0.003874  0.001687   258.2 
  1000.     0.086303  0.021288 -0.004146  0.002099   290.1 
  2000.     0.103743  0.030909 -0.003141  0.002535   318.0 
  5000.     0.127172  0.049907  0.002019  0.003469   352.1 
 10000.     0.149810  0.074734  0.011546  0.004939   382.2 
 20000.     0.185706  0.119913  0.030938  0.007877   425.5 
 50000.     0.283361  0.250283  0.089327  0.016714   525.6 

RUNNING THE "ONE-POINT" MODE

A parameter initialization file Mone.inp is provided. It sets the KRC system into a reasonable mode for one-point calculations. Do not change that file unless you have read this entire file.

A line near the end of that file points to a file ’one.inp’ which can contain any number of one-point conditions. ’one.inp’ is intended to be edited to contain the cases you want; however, it must maintain the input format of the sample file.

First Line is any title you wish. It must be present.
The second line is an alignment guide for the location lines. It must be there.

Each following line must start with an ’11 ’; this is a code that tells the full-up KRC that is a one-point line. The next 9 fields are read with a fixed format, and each item should be aligned with the last character of the Column title. All items must be present, each line must extend at least to the m in Azim; comments may extend beyond that, but they will not appear in the output file. Be sure to have a <math><</math>CR<math>></math> at the end of the last input line.

The fields (after the 11) in the one-point input are:
    Ls  L_sub_S season, in degrees
    Lat Aerographic latitude in degrees
    Hour    Local time, in 1/24'ths of a Martian Day
    Elev    Surface elevation (relative to a mean surface Geoid), in Km
    Alb Bolometric Albedo, dimensionless
    Inerti  Thermal Inertia, in SI units
    Opac    Atmospheric dust opacity in the Solar wavelength region
    Slop_   Regional slope, in degrees from horizontal
    Azim    Azimuth of the down-slope direction, Degrees East of North.

The two additional columns in the output file are:
    TkSur   Surface kinetic temperature
    TbPla   Planetary bolometric brightness temperature

Try running the binary file first. If that fails, a Makefile is provided to complile and link the program; simply enter "make krc" and pray. If this fails, have your local guru look over the Makefile for local dependancies. Suggestions of making the Makefile more universal are welcome.

To run the program, change to the directory where the program was built, and enter "krc". You should get a prompt: ?* Input file name or / for default = Mone.inp
If the initialization file still has this name and is in the same directory, enter a single "/" and <math><</math>CR<math>></math>. Otherwise, enter the full pathname to the initialization file, with no quotes and no blanks.

A second prompt is for the name of the output file: ?* Print file name or / for default = krc.prt
Again, if this is satisfactory, simply enter / <math><</math>CR<math>></math> , else enter the desired file path-name.

Comments on the One-point model

The initialization file of 2002mar08 is set to compute the temperatures at the season requested without seasonal memory. It uses layers that extend to 5 diurnal skin depths. It does not treat the seasonal frost properly, so don’t believe the results near the edge of the polar cap. Execution time on a circa 2001 PC may be the order of 0.01 seconds per case.

The underlying model is the full version of KRC. By modifying the initialization file, you can compute almost anything you might want. If you choose to try this, best to read all of this document.

DEBUG OPTIONS [debug]

tcard.f:75:     IF (IDB2.GE.5) WRITE(IOSP,*) 'TCARD-A',IQ
tcard.f:123:    IF (IDB1.GE.1) PRINT *,'Before PORB0'
tcard.f:125:    IF (IDB1.GE.1) PRINT *,'AFTER PORB0'
tcard.f:349:    IF (IDB1.NE.0) WRITE(IOSP,*)'TCARD Exit: IRET=',IRET,NFD,ID(1) 
tday.f:63:      IF (IDB2.GE.5) WRITE(IOSP,*) 'TDAY IQ,J4=',IQ,J4,jjo
tday.f:544: 9   IF (IDB2.GE.6) WRITE(IOSP,*) 'TDAYx'
tdisk.f:90:     IF (IDB3.NE.0) WRITE(IOSP,*)'TDISKa ',KODE,KREC,NCASE,J5,K4OUT
tdisk.f:424:       IF (IDB3.GE.3) WRITE(IOSP,*)'TDISKc  KREC=',KREC,LOPN2,IOD2,I
tdisk.f:431:    IF (IDB3.GE.3) WRITE(IOSP,*)'TDISKx  KREC=',KREC
tlats.f:56:     LQ1=IDB2.GE.3           ! once per season or latitude
tlats.f:+       LQ2=IDB2.GE.6           ! each day
tlats.f:+       IF (IDB2.NE.0) WRITE(IOSP,*)'TLATSa',N3,N4,J5,LATM,LQ1,LQ2
tlats.f:422: 9      IF (IDB2.GE.3) WRITE(IOSP,*)'TLATSx',N1,N1PIB,N2,N24,J3
tseas.f:41:     IF (IDB1.NE.0) WRITE(IOSP,*)'MSEASa',IQ,IR,J5,LSC,N5,LONE
in tlats.f:
 98:    IF (LQ1) THEN
       WRITE(75,*) 'J5+',J5,SUBS,SLOPE,SLOAZI,SKYFAC
       WRITE(75,*) 'MXX+',MXX,DIP,SAZ
       WRITE(75,*) 'PXX+',PXX
top of Lat loop
130:    IF (LQ1) PRINT *,'TLAT1 J5,TBLOW=',J5,TBLOW
170:    IF (LQ1) THEN
       WRITE(75,*)'FXX+',FXX,J4,DLAT
       WRITE(75,*)'RXX=',RXX  ! R should be 90 deg from F
       WRITE(75,*)'TXX=',TXX
221:    IF (LQ1) print *,'TLATS: J4,SOLR...',J4,SOLR,ACOSLIM,COSIAM(1)
 top of time loop
299:       IF (LQ1.AND.(MOD(JJ,24).EQ.1)) THEN
              WRITE(75,*)'HXX+',HXX,JJ
              WRITE(75,*)'ANG:',ANGLE,COSI,COS2,DIRECT,QI
303:    IF (LQ2) WRITE(IOSP,*),'TLatc',JJ,COSI,COS3,DIRECT,DIFFUSE 
309:    IF (LD19) WRITE(79,777) QA,QI,DIRECT,DIFFUSE,BOUNCE
 end time loop
341:    IF (LQ1) then 
           PRINT *,'AVEA ...',AVEA,AVEE,AVEI,AVEH
           PRINT *,'CABR...',CABR,TAUD,TAUIR,FACTOR,TAUEFF
           PRINT *,'BETA...',BETA,QS,SIGSB 
           PRINT *,'TAEQ4,TSEQ4,TEQUIL',TAEQ4,TSEQ4,TEQUIL
356:       IF (LQ1) PRINT *,'TSUR,TBOT',TEQUIL,TSUR,TBOT
           IF (LQ1) PRINT *,'XCEN',XCEN 
379:    IF (LQ1) PRINT *,'TTJ',TTJ
453:    IF (LD16) THEN
       WRITE(76,761)SUBS,DLAT,ALB,SKRC,TAUD,PRES
           loop on hour: WRITE(76,762)QH,TSFH(I),ADGR(J),QS,TPFH(I)
end lat loop

Set LD19 to write bottom-of-atmosphere downgoing fluxes to separate file for every time-step for every latitude, every time tlats is called.

READING TYPE 5x FILES

IDL routines do not access files directly unless specifically listed.

DEFINEKRC Define structures in IDL that correspond for Fortran commons
Calls: None == None other than IDL library
Firm code of common definitions. Must be recoded if a Fortran *.inc changes.

READKRCCOM Read a KRCCOM structure from a bin5 file
Uses 3-element HOLD array. Returns a structure of krccom
Options to open or close bin5 file or read one case
Calls: DEFINEKRC
Files: bin5
HOLD is: 0]=logical unit 1]=number of words in a case 2]=# cases in the file

MAKEKRCVAL
Make string of selected KRC inputs: Key=val
Calls: DEFINEKRC

KRCHANGE Find changes in KRC input values in common KRCCOM
Calls: READKRCCOM MAKEKRCVAL
Reads and stores krccom for first case. For each additional case, makes a list of any changes in the flaot, integer or logical input values.

KRCCOMLAB Print KRC common input items
all items via arguments
Calls: None

KRCLAYER Compute center depth of KRC layers
all items via arguments
Calls: None

KRCCOMLAB Print KRC common input items
all items via arguments
Calls: None

KRCSIZES Compute array and common sizes for KRC Fortran
Test procedure to compute array sizes or hours.
Must recode if any size in *.inc changes
Calls: None

NOTES ON HOW SOME ASPECTS OF THE CODE WORK

New file name

TCARD reads a card of Type 8, (and index is not 22 or 23) it calls TDISK(4,0), which closes current file and sets LOPN2=.FALSE. TCARD then moves new file name into common.
KRC checks if current (new) values of N5 and JDISK call for file output; with LOPN2=.FALSE., KRC calls CALL TDISK (1,0) to open new file.

End of a case and end of a run

TCARD sets KOUNT=0 at entry; this is incremented for every card except those of type 0 (or less) or type 11 (one-point mode). When type 0 is encountered, if KOUNT is positive, does normal check of changes before return with IR=1 to indicate start of a new case; if KOUNT is zero, returns with IR=5 and prints ’END OF DATA ON INPUT UNIT.'

Setting one-point mode

This can be done only in the first case, and there is no way to leave the one-point mode except to end the run.

TCARD encounters: "10 * filename" as change card in the initial case. sets this as new input file name, then returns with IRET=4.
KRC closes prior input file, opens the new one, and reads past first two lines then calls TCARD to read first one-point line and sets LONE=true and drops into the top of the "case" loop.
The master one-point should have a single latitude, no binary output file.
The small number of layers, days to converge, and seasons ignores the seasonal effect.

One-point request values are read by TCARD @ 310, which computes starting DJUL.

TPRINT does linear interpolation of TOUT, which has N2 points be sol. To get Tp, does interpolation of Tp-Ts at the hour points, and adds to interpolated Ts.

Starting conditions and date

Initial N5-JDISK sets the size of output files. There could be any number of interior seasons where parameter changes are made; based on successive values of IDOWN.

KRC initially calls TCARD(1
For each case loop, sets IQ=TCARD_return. If one-point mode, sets IQ=1.

TSEAS uses IQ as key. If this is 1, then sets J5=0 and sets DJU5 to season -1., else, increments J5 and increments DJU5 with current DELDUL. This allows use of variable resolution dates. (so J5 never 0 when TCARD(2 called)
If J5 equals IDISK2, then TSEAS calls TCARD(2 to read changes, and proceeds to next season.

TLATS uses J5 as the key; if it is <math><= 1</math>, then starts from equilibrium conditions, else uses predictions from prior season.

The default is that change cards cause a fresh calculation of starting conditions. Exceptions are when J5=IDOWN<math>>0</math> at TCARD entry.

Changing parameters within a seasonal run = Continue from memory.

When J5 reaches IDOWN, TSEAS calls TCARD, which will set IRET=3 before reading the new parameters. May change DELJUL to get finer seasonal resolution, but must NOT change N5

Use: Normal restrictions for what may not change for Type 5x files apply. E.g., type 56 must NOT change number of latitudes nor total number of seasons.

Set N5 to be the total number of seasons desired, including those after any number of parameter changes; it must NOT be changed later.

Set IDOWN to the season at the beginning of which wish to (first) change parameters. The next set of changes could include a revised (larger) IDOWN.

Use of common PORBCM

Contents are described in porbcm.inc
PORBCM is filled by TCARD calling PORB0, which reads the first 60 items in 5G15.7 from the input file and sets the value of PI. KRC references porbcm.inc but does not use it.

TSEAS uses a few items to calculate LsubS; 
  CALL PORB (DJU5,DAU,PEA,PEB,HFA,HFB)  ! FIND CURRENT PLANET POSITION 
  TANOM = ATAN2 (PHFXX(2),PHFXX(1))     ! TRUE ANOMOLY
  SUBS = MOD(((TANOM+SLP)*RAD)+360.,360.)       ! L-SUB-S IN DEGREES
    SLP        L-SUB-S AT PERIAPSIS
    PHFXX(3)   VECTOR FROM FOCUS (SUN) TO PLANET: ORBITAL PLANE COORDINATES
[ TYEAR uses the value for length of year. ]

Lower boundary condition and resetting (jumping) layer temperatures.

At the start of a case, TLATS sets the temperature profile linear with depth in one of three ways:

IB=0: top and bottom at equilbrium temperature
IIB=1: top at equilbrium temperature, the bottom at TDEEP
IIB=2: top and bottom at TDEEP

The kind of resetting is controlled by IB. In TCARD, if IB>0, then N1PIB=N1+1, else N1PIB=N1. T(N1+1) is not reset in the time calculations. In TDAY, for each time step, the temperature of the lower boundary is set equal to T(N1PIB), which results in either zero heat flow (IB=0) or a constant temperature.

Seasonal variation of albedo or opacity or “climate”

When TCARD encounters a type of 8 and an index of 23(tau) [or 22(albedo)], it transfers the text item into FVTAU (which is in COMMON /FILCOM/) and then calls SEASTAU with an Ls of -999 . SEASTAU when called with LSUB LT -90 calls (providing IOD3) READTXT360, which reads file. Maximum number of rows is 360, more will be ignored. First and last entry read are wrapped with +,- 360 to Ls to ensure no interpolation faults later. TCARD sets the variable Tau flag, KVTAU, true if table-read was successful, else it is set false. If KVTAU is set, TSEAS calls SEASTAU at start of each season, resetting TAUD.

If type 8 and index 23, the same as above except names are -ALB rather than -TAU.

If type 8 and index 24, then TCARD calls CLIMTAU to open and read a .bin5 climate file, and sets KVTAU=2.

CLIMTAU expects to read a .bin5 file (season,latitude,2) with dust and ice infrared opacities; this file is normally made by the IDL routine mopacity.pro. Season is assumed to the uniform in Ls from 0 to 360-delta and latitude assumes to be uniform from -90+delta/2 to +90-delta/2. CLIMTAU has firm-coded sizes, 72 seasons and 36 latitudes, and stores the file. Upon later calls from TLATS, it returns the two opacities at a requested Ls and latitude, using bi-linear interpolation.

Cap-dependent pressure

TSEAS:   BUF(1)=0.       ! flag for  TINT to compute areas
IF (N4.GT.8) CALL TINT (FROST4, BUF, SUMF)

Tlats
        IF (N4.GT.8) THEN       ! use global integrations
          PCAP = SUMF*GRAV      ! cap_frost equivalent surface pressure

 KPREF.EQ.2
 PZREF = PTOTAL - PCAP
          PCO2G = PCO2M -PCAP ! all changes are pure CO2


Personal tools