#Original Function Location: /themis/lib/dav_lib/public_library/latest/thermal_model.dvrc
/themis/lib/dav_lib/library/thermal_model.dvrc
Source Code for Function: "krc_one()" in "thermal_model.dvrc" (Public)
thermal_model_version=1.01
# load_krc - taken from r.fergason
# jd2lsubs - taken from r.fergason
# krc_themis - taken from r.fergason
# krc_one
# krc_bin
# process_bin52
# tmodel
# deddi
#mfp - taken from s. piqueux
#pore - taken from s. piqueux
#Kn - taken from s. piqueux
#k0 - taken from s. piqueux
#sigmoid - taken from s. piqueux
#k_gas - taken from s. piqueux
#k_radsol - taken from s. piqueux
#k_SCP - taken from s. piqueux
#k_BCCP - taken from s. piqueux
#k_CCP - taken from s. piqueux
#k_cell - taken from s. piqueux
#k_tot - taken from s. piqueux
#Cp - taken from s. piqueux
#I2d - taken from s. piqueux
#1.01 - 1/30/12 - added thermal model functions provided by s.piqueux
define load_krc(file_name, seas, head_byte, type) {
if ($ARGC == 0) {
printf(" \n")
printf("This script will read in binary krc K4OUT -1 \n")
printf("files or TES production tables into davinci. \n")
printf("This script will only work on a Linux machine. \n")
printf(" \n")
printf("Usage: a=load_it(file_name, seasons, max_lat, type) \n")
printf(" where file name = name of file within quotes \n")
printf(" seasons = number of seasons (687/DELJUL) \n")
printf(" latitudes = maximum number of latitudes (MAXNH) \n")
printf(' type = "krc" or "tes" \n')
printf(" \n")
printf("Assuming krc tables were produced on a Linux machine, \n")
printf("and TES production tables were produced on East. \n")
printf(" \n")
printf("2004 Oct 15 - r.fergason, initial development \n")
printf(" \n")
return(NULL)
}
# Leins:
# This script will only work on a Linux machine.
# This script will only work if MAXN4=24 and K4OUT is -1.
# Assuming krc is produced on a Linux machine (Little Endian)
# Assuming TES Production Tables were producted on east (Big Endian)
# Typical krc (K4OUT -1) header size for MAXNH=19 = 3648 bytes (912*4)
# Typical krc (K4OUT -1) header size for MAXNH=37 = 7104 bytes (1776*4)
# Typical TES production header size = 7104 bytes (1776*4)
file_name = $1
seas = $2
lats = $3
type = $4
if (type == "krc") {
# This is for typical krc K4OUT -1 output
maxnh=(lats)
maxn4=24
step=(maxnh*maxn4)*4
head_bytes=(maxnh*maxn4)*2*4
# Assuming krc were producted on a Linux machine. Format is lsb (Least Significant Byte First = Little Endian)
fmt="lsb_float"
} else if (type == "tes") {
# This is for TES Production Table output
maxnh=37
maxn4=24
step=(maxnh*maxn4)*4
head_bytes=(maxnh*maxn4)*2*4
# TES Production tables were calculated on east. Format is msb (Most Significant Byte First = Big Endian)
fmt="msb_float"
fmt="msb_float"
}
# create a matrix to store the temperature data; tsf-surface temperature, tpf-planetary temperature
tsf=create(maxnh, maxn4, seas, start=0, step=0, format=float)
tpf=create(maxnh, maxn4, seas, start=0, step=0, format=float)
# For loop designed to split the file, by season, into a structure
skip=head_bytes
for (i=1; i<=seas; i+=1) {
# Translated because fortran is column-major, davinci is row-major
tsf[,,i]=translate(load_raw(file_name,x=maxn4,y=maxnh,z=1,org=bsq,format=fmt,header=skip),x,y)
skip=skip+step
tpf[,,i]=translate(load_raw(file_name,x=maxn4,y=maxnh,z=1,org=bsq,format=fmt,header=skip),x,y)
skip=skip+step
}
# convert the information into a structure
output=struct(surf_temp=1, planet_temp=2, type=3)
output.surf_temp=tsf
output.planet_temp=tpf
output.type=type
return(output)
}
define jd2lsubs() {
###########################
#
# Hugh Kieffer's idl function l_sub_s.pro converted
# to a davinci script.
#
# r.fergason
# 2005 Aug 26
#
###########################
if ($ARGC == 0) {
printf(" \n")
printf("usage: a=lsubs(julian date) \n")
printf("returns aerographic longitude of the sun, in degrees \n")
printf(" $1 = Julian date, no offset - recomend double precision \n")
printf(" \n")
return(NULL)
}
date = $1
pi=acos(-1)
# fractional year
x=(date-2450322.34)/686.98
# insure in range 0 to 1
x=x-floor(x)
# radians from start of martian year
x=x*2.0*pi
a= create(8,1,1, start=0, step=0, format=double)
a[1,1,1]= -10.328371
a[2,1,1]= 57.296001
a[3,1,1]= -10.688416
a[4,1,1]= -3.2931221
a[5,1,1]= -0.62748339
a[6,1,1]= -5.0168313
a[7,1,1]= -0.050566287
a[8,1,1]= -0.41753547
f = a[1,1,1] + a[2,1,1]*x + a[3,1,1]*cos(x-a[4,1,1]) + a[5,1,1]*cos(2.*x-a[6,1,1]) + a[7,1,1]*cos(3.*x-a[8,1,1])
return(f)
###############
#
# Hugh's notes:
#
# 1. coefficents derived in fit_lsubs 97dec17 to porb.prt
# x= 0.00910607* ( date- 11010.0)
# a= [ -9.86388 , 57.5329 , -10.7144 , -3.31933 , -0.636553 , -5.10982]
#
# 2. following are fits to MICA run, every 5 days
# Inital data file has L_s to 2 decimal places., this dominates error!
# dblarr
# chisq= 2.1078586e-05
# sigma= 0.44947921 0.13704680 0.17420014 0.029852505 0.21318231 0.27769774 0.17337688 3.7167788
# Maximum errors= -0.013322718 0.012496469
# average absolute error= 0.0037515840
# form is a0 +a1 x + a2 cos(x-a3)+a4 cos ( 2x-a5) +... a_even cos(nx-a_odd)
# a= -10.328371 57.296001 -10.688416 -3.2931221 -0.62748339 -5.0168313 -0.050566287 -0.41753547
# Before applying formula, MUST scale date by:
# radian = 0.00914610* [ Julian date- 2450322.34 ]
#
##############
}
define krc_themis(check) {
# Added help and modified tail command to work on other linux systems
# J Bandfield 5/2008
if ($ARGC == 0) {
printf(" \n")
printf("Runs one point krc thermal model\n")
printf("$1 - ls\n")
printf("$2 - local time\n")
printf("$3 - local slope (degrees from horizontal)\n")
printf("$4 - azimuth of slope (degrees clockwise from north)\n")
printf("$5 - visible wavelength dust opacity\n")
printf("$6 - latitude\n")
printf("$7 - elevation (km)\n")
printf("$8 - albedo\n")
printf("$9 - number of inertia steps (increments of 20)\n")
printf("$10 - starting inertia value\n")
printf("check = 1 - Check table for constantly increasing values\n")
printf("This should be selected for daytime measurements (default = 1)\n\n")
printf("This script will only work on a Linux machine. \n")
printf("\n")
return(NULL)
}
if (HasValue(check)==0) {
check = 1
}
ls = $1
hour = $2
slope = $3
azim = $4
opac = $5
lat = $6
elev = $7
alb = $8
#num = $9
length = $9
value = $10
num=1
# create data portion of krc input file
#length = 99
#value=20
table=sprintf("%2d %5.1f %5.1f %4.1f %5.2f %4.2f %6.1f %4.2f %4.1f %4.0f", 11,ls,lat[,num,],hour,elev[,num,],alb,value,opac,slope,azim)
for (i=1; i<=length; i+=1) {
value=value+20
line=sprintf("%2d %5.1f %5.1f %4.1f %5.2f %4.2f %6.1f %4.2f %4.1f %4.0f", 11,ls,lat[,num,],hour,elev[,num,],alb,value,opac,slope,azim)
table=cat(table, line, axis=Y)
}
mytmpdir=$TMPDIR//"/krc"
system(sprintf("mkdir %s", mytmpdir))
write(table, filename=mytmpdir//"/table", type=ascii, force=1)
system(sprintf("cp %s/krc %s/krc",$DV_KRC,$TMPDIR))
system(sprintf("cp %s/krc.inp %s/krc", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/fake_out_krc %s/krc", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/one.header %s/krc", $DV_KRC, $TMPDIR))
system(sprintf("cd %s/krc; cat one.header table > one.inp", $TMPDIR))
system(sprintf("cd %s/krc; krc < fake_out_krc", $TMPDIR))
system(sprintf("cd %s/krc; tail -n +10 < krc.prt > output", $TMPDIR))
# This line will only work on some systems -- not Solaris!
system(sprintf("cd %s/krc; grep '11' output > table", $TMPDIR))
table=read_ascii(sprintf("%s/table", mytmpdir), format=float)
system(sprintf("rm -rf %s", mytmpdir))
# check that the values are always increasing and that there are no repeated temperature values at the top of the table
val=1
if (check==1) {
flag=1
maxval=dim(table)[2]
while (flag!=0 && val
if (table[11,val+1,] <= table[11,val,]) {
val=val+1
} else {
if (table[11,val+1,] > table[11,val,]) {
flag=0
}}
}
}
table=table[,val:length+1,]
return(table)
}
define krc_one(struct,alb,lat,lon,elev,tau,slope,azim,hour,bol,ls) {
if ($ARGC == 0) {
printf("\n")
printf("This runs krc in one point mode and derives TI from input temperature.\n")
printf("\n")
printf("$1 = surface kinetic or bolometric temperature\n")
printf("lat = latitude (default = 0)\n")
printf("lon = east longitude (default = 0)\n")
printf("elev = elevation km (default from lat,lon 2ppd MOLA)\n")
printf("alb = albedo (default from lat,lon 2ppd TES)\n")
printf("tau = dust visible opacity (default = 0.30)\n")
printf("slope = surface slope in degrees (default = 0)\n")
printf("azim = slope azimuth in degrees east from north (default = 0)\n")
printf("hour = hour of day for data (default=2)\n")
printf("bol = match to bolometric temps (default=0)\n")
printf("ls = Lsubs (default=0)\n")
printf("struct = Output ancillary info (default=0)\n")
printf("\n")
printf("\n")
return(0)
}
temp=$1
if (HasValue(bol)==0) {
bol = 0
}
if (HasValue(struct)==0) {
struct = 0
}
if (HasValue(lat)==0) {
lat = 0
}
if (HasValue(lon)==0) {
lon = 0
}
# Get elevation and albedo and ti from maps
tin=read(sprintf("%s/krc_bin/ti_map2ppd_v3.vicar", $DV_KRC))
albedo=read(sprintf("%s/krc_bin/albedo_2ppd.vicar", $DV_KRC))
elevation=read(sprintf("%s/krc_bin/mola_2ppd.vicar", $DV_KRC))
image_xy=geo_trans(lat,360-lon,2)
tin=tin[image_xy[1,,],image_xy[2,,],]
albedo=albedo[image_xy[1,,],image_xy[2,,],]
elevation=elevation[image_xy[1,,],image_xy[2,,],]/1000
if (HasValue(alb)==0) {
alb = albedo
}
if (HasValue(elev)==0) {
elev=elevation
}
printf("\nTES Inertia: %f\nTES Albedo: %f\nElevation (km): %f\n\n" , tin, albedo, elev)
if (HasValue(tau) == 0) {
tau=0.30
}
if (HasValue(ls) == 0) {
ls = 0
}
if (HasValue(slope) == 0) {
slope = 0
}
if (HasValue(azim) == 0) {
azim = 0
}
if (HasValue(hour) == 0) {
hour = 2
}
table=krc_themis(ls,hour,slope,azim,tau,lat,elev,alb,99,20,check=0)
# Check to see whether values are increasing or decreasing or are static at co2 temps.
stop=100
noco2=0
reverse=0
for (i=1; i<=99; i+=1) {
if (table[11,i,]!=table[11,i+1,] && HasValue(start)==0) {
noco2=1
start=i
if ((table[11,i,]>table[11,i+1,])) {
reverse=1
}
}
if (noco2==1 && reverse==1 && (table[11,i,]
stop=i
i=101
}
if (noco2==1 && reverse==0 && (table[11,i,]>table[11,i+1,])) {
stop=i
i=101
}
}
if (HasValue(start)==0) {
start=1
}
printf(""+start+"\n")
printf(""+stop+"\n")
printf(""+reverse+"\n")
printf(""+noco2+"\n")
ti=table[7,start:stop]
tk=table[11,start:stop]
tbol=table[12,start:stop]
if (reverse==1) {
ti=translate(ti,y,y,flip=1)
tk=translate(tk,y,y,flip=1)
tbol=translate(tbol,y,y,flip=1)
}
if (bol==1) {
tscale=tbol
} else {
tscale=tk
}
if (temp
if (reverse==0) {
out=0
} else {
out=2000
}
}
if (temp>max(tscale)) {
if (reverse==0) {
out=2000
} else {
out=0
}
}
if (temp>=min(tscale) && temp<=max(tscale)) {
out=interp(ti,tscale,temp)
}
if (struct==1) {
tout=out
out=struct(ti_scale,tk,tbol,tes_ti,alb,elev)
out.ti=tout
out.tes_ti=tin
out.alb=alb
out.elev=elev
out.ti_scale=ti
out.tk=tk
out.tbol=tbol
}
return(out)
}
define krc_bin(lat,lon,elev,alb,tau,ti,jbare,tdeep,slope,azimuth,layer,out,type,hour,info,mult,usage) {
if(HasValue(usage)==0) {
usage=0
}
if (usage == 1) {
printf("\n")
printf("This runs krc in binary mode with the following inputs.\n")
printf("\n")
printf("lat = latitude (default = 0)\n")
printf("lon = east longitude (default = 0)\n")
printf("elev = elevation km (default from lat,lon 2ppd MOLA)\n")
printf("alb = albedo (default from lat,lon 2ppd TES)\n")
printf("tau = dust visible opacity (default = 0.30)\n")
printf("ti = thermal inertia of upper layers (default from 2ppd TES)\n")
printf("jbare = force frost free at specified ls (default = no forcing)\n")
printf("tdeep = set bottom layer temperature (default = insulating)\n")
printf("slope = surface slope in degrees (default = 0)\n")
printf("azimuth = slope azimuth in degrees east from north (default = 0)\n")
printf("layer = layer number of changed properties (default = no layering)\n")
printf("out = bin5 output type, 52 or 53 (default = 52)\n")
printf("type = temperature type, 1 - kinetic, 2 - bolometer (default = 1)\n")
printf("hour = hour of day for data when out = 52, 0 = all 24 hours (default=0)\n")
printf("mult = run layers 6,7,8,9,10,12,14,16,18,20 0 - no, 1 - yes (default=0)\n")
printf("\n")
printf("The model has 29 layers. First layer is set at 0.2 skin depths.\n")
printf("Each successive layer is increased by a factor of 1.15.\n")
printf("\n")
printf("The model iterates for 2 years before outputting a full year of data\n")
printf("at 1/80 year intervals, starting at 0 with 24 hours from 1 to 24. Albedo and tau can be dynamic if input as \n")
printf("2 x n x 1 arrays. Ls = column 1 and alb/tau = column 2\n")
printf("\n")
return(0)
}
# Copy files over to temp directories
system(sprintf("mkdir %s/krc_bin", $TMPDIR))
system(sprintf("cp %s/krc %s/krc_bin/", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/krc_bin/krc_part1.inp %s/krc_bin/", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/krc_bin/krc_part2.inp %s/krc_bin/", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/krc_bin/krc_part3.inp %s/krc_bin/", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/krc_bin/var_header.ascii %s/krc_bin/", $DV_KRC, $TMPDIR))
system(sprintf("cp %s/krc_bin/fake_out_krc_bin %s/krc_bin/", $DV_KRC, $TMPDIR))
if (HasValue(lat) == 0) {
lat = 0
}
if (HasValue(lon) == 0) {
lon = 0
}
# Get elevation and albedo and ti from maps
tin=read(sprintf("%s/krc_bin/ti_map2ppd_v3.vicar", $DV_KRC))
albedo=read(sprintf("%s/krc_bin/albedo_2ppd.vicar", $DV_KRC))
elevation=read(sprintf("%s/krc_bin/mola_2ppd.vicar", $DV_KRC))
image_xy=geo_trans(lat,360-lon,2)
tin=tin[image_xy[1,,],image_xy[2,,],]
albedo=albedo[image_xy[1,,],image_xy[2,,],]
elevation=elevation[image_xy[1,,],image_xy[2,,],]/1000
if (HasValue(ti)==0) {
ti = tin
}
if (HasValue(alb)==0) {
alb = albedo
}
if (HasValue(elev)==0) {
elev=elevation
}
printf("\nTES Inertia: %f\nTES Albedo: %f\nElevation (km): %f\n\n" , tin, albedo, elev)
if (HasValue(tau) == 0) {
tau=0.30
}
if (HasValue(jbare) == 0) {
jbare = 0
}
if (HasValue(tdeep) == 0) {
tdeep = 0
}
if (HasValue(slope) == 0) {
slope = 0
}
if (HasValue(azimuth) == 0) {
azimuth = 0
}
if (HasValue(layer) == 0) {
layer = 0
}
if (HasValue(mult) ==0) {
mult = 0
}
if (HasValue(out) == 0) {
out = 52
}
if (HasValue(hour) == 0) {
hour = 0
}
if (HasValue(info) == 0) {
info = 0
}
# Set up switch for layering case
if (mult == 1) {
layer= 6
}
# Construct latitude and elevation lines
lat=printf("%7.2f\n",lat)
elev=printf("%7.2f\n",elev)
# Start construction of parameter change lines
params="2 5 240 'N5' / seasons"
params=cat(params,sprintf("2 12 161 'JDISK' /"),y)
params=cat(params,sprintf("1 42 8.5875 'DELJUL 1/80 year' /"),y)
params=cat(params,sprintf("1 3 %.2f 'TI' /",ti),y)
if (tdeep!=0) {
params=cat(params,sprintf("1 15 %.2f 'TDEEP' /",tdeep),y)
} else {
params=cat(params,sprintf("2 7 0 'IB' /"),y)
}
if (jbare !=0) {
params=cat(params,sprintf("2 18 %i 'JBARE'/",jbare),y)
}
params=cat(params,sprintf("1 23 %.2f 'Slope' /",slope),y)
params=cat(params,sprintf("1 24 %.2f 'Azimuth' /",azimuth),y)
params=cat(params,sprintf("2 8 %i 'IC' /",layer),y)
params=cat(params,sprintf("2 17 %i 'K4out' /",out),y)
params=cat(params,sprintf("8 0 0 '%s/krc_bin/outdata.bin5' /",$TMPDIR),y)
# Construct parameters for changing albedo or dust opacity values
# Start with formatting the albedo and dust arrays into the correct spacing and formats
if (max(dim(tau)) ==1) {
params=cat(params,sprintf("1 17 %.2f 'TauD' /",tau),y)
} else {
write(tau,sprintf("%s/krc_bin/taufile.tab", $TMPDIR),type=ascii,force=1)
system(sprintf("cd %s/krc_bin; cat var_header.ascii taufile.tab > taufile", $TMPDIR))
# params=cat(params,sprintf("0/"),y)
params=cat(params,sprintf("8 23 0 '%s/krc_bin/taufile' /", $TMPDIR),y)
}
if (max(dim(alb)) ==1) {
params=cat(params,sprintf("1 1 %.2f 'Albedo' /",alb),y)
} else {
write(alb,sprintf("%s/krc_bin/albfile.tab", $TMPDIR),type=ascii,force=1)
system(sprintf("cd %s/krc_bin; cat var_header.ascii albfile.tab > albfile", $TMPDIR))
# params=cat(params,sprintf("0/"),y)
params=cat(params,sprintf("8 22 0 '%s/krc_bin/albfile' /", $TMPDIR),y)
}
if (mult==1) {
for (i=6; i<=10; i+=1) {
params=cat(params,"2 5 240 'N5' / seasons",y)
params=cat(params,sprintf("1 42 8.5875 'DELJUL 1/80 year' /"),y)
params=cat(params,sprintf("2 8 %i 'IC' /",i),y)
params=cat(params,sprintf("0/"),y)
}
for (i=12; i<=20; i+=2) {
params=cat(params,"2 5 240 'N5' / seasons",y)
params=cat(params,sprintf("1 42 8.5875 'DELJUL 1/80 year' /"),y)
params=cat(params,sprintf("2 8 %i 'IC' /",i),y)
params=cat(params,sprintf("0/"),y)
}
}
params=cat(params,sprintf("0/"),y)
params=cat(params,sprintf("0/"),y)
params=cat(params,sprintf("0/"),y)
params=cat(params,sprintf("0/"),y)
# Start constructing krc_inp file using system commands
write(sprintf("%7.2f\n",lat),$TMPDIR+"/krc_bin/lat.inp",type=ascii,force=1)
write(sprintf("%7.2f\n",elev),$TMPDIR+"/krc_bin/elev.inp",type=ascii,force=1)
write(params,$TMPDIR+"/krc_bin/params.inp",type=ascii,force=1)
system(sprintf("cd %s/krc_bin; cat krc_part1.inp lat.inp > out1", $TMPDIR))
system(sprintf("cd %s/krc_bin; cat out1 krc_part2.inp > out2", $TMPDIR))
system(sprintf("cd %s/krc_bin; cat out2 elev.inp > out1", $TMPDIR))
system(sprintf("cd %s/krc_bin; cat out1 krc_part3.inp > out2", $TMPDIR))
system(sprintf("cd %s/krc_bin; cat out2 params.inp > krc.inp", $TMPDIR))
system(sprintf("cd %s/krc_bin; %s/krc_bin/krc < fake_out_krc_bin", $TMPDIR, $TMPDIR))
#system(sprintf("cd %s/krc_bin; krc < fake_out_krc_bin", $TMPDIR))
if (out==52) {
data=process_bin52(sprintf("%s/krc_bin/outdata.bin5", $TMPDIR),hour=hour,type=type,latitude=1,raw=0)
}
if (out==53) {
# Create ls axis for resampling
et=13046146
for (i=1; i<=79; i+=1) {
et=cat(double(et),double(et[,,1]+(3600*24*8.58725*i)),z)
}
oldls=lsubs(et)
oldls[,,1]=0
newls=create(1,1,80,start=0,step=4.5,format=float)
if (mult ==0) {
data=struct(hour,layermin,layermax,insolation)
data.hour=load_bin5(sprintf("%s/krc_bin/outdata.bin5", $TMPDIR))
data.hour=resample(translate(data.hour[,3:],y,z),oldls,newls).data
data.insolation=data.hour[110:,1]
data.layermin=data.hour[50:78,1]
data.layermax=data.hour[80:108,1]
data.tbol=data.hour[26:49,1]
data.tk=data.hour[2:25,1]
data.seasonmax=max(data.layermax,axis=y)
data.seasonmin=min(data.layermin,axis=y)
} else {
data=struct(hour,layermin,layermax,insolation)
data.hour=load_bin5(sprintf("%s/krc_bin/outdata.bin5", $TMPDIR))
data.hour=resample(translate(data.hour[,3:],z,y),oldls,newls).data
data.insolation=data.hour[110:,]
data.layermin=data.hour[50:78,]
data.layermax=data.hour[80:108,]
data.tbol=data.hour[26:49,]
data.tk=data.hour[2:25,]
data.seasonmax=max(data.layermax,axis=y)
data.seasonmin=min(data.layermin,axis=y)
}
remove_struct(data,"hour")
}
# Get Layer Info
if (info=1) {
line=atoi(syscall(sprintf("nl -b a %s/krc_bin/krc.prt | grep LAYER", $TMPDIR)))
linestart=line+1
nlines=dim(linestart)[2]
layer_info=create(7,29,1,start=0,step=0,format=float)
for (i=1; i<=nlines; i+=1) {
syscall(sprintf("tail -n +%i %s/krc_bin/krc.prt | head -29 > %s/krc_bin/layerinfo.ascii", linestart[,i], $TMPDIR, $TMPDIR))
layer_temp=ascii(sprintf("%s/krc_bin/layerinfo.ascii", $TMPDIR),x=7,y=29,z=1,format=float)
layer_info=cat(layer_info,layer_temp,z)
}
layer_info=layer_info[,,2:]
temp=data
data=struct(data,layer_thickness,center_depth)
data.data=temp
data.layer_thickness=translate(layer_info[2:3],x,y)
data.center_depth=translate(layer_info[4:6],x,y)
data.layer_top=data.center_depth[,1:2]-(data.layer_thickness[,1:2]*0.5)
data.alb=alb
}
# Clean up the mess that we've created
system(sprintf("rm -rf %s/krc_bin", $TMPDIR))
return(data)
}
define process_bin52(hour,type,latitude,raw) {
if ($ARGC == 0) {
printf("\n")
printf("Read in bin5 krc file and process to something more manageable\n")
printf("$1 = filename\n")
printf("hour = local_time (default=2)\n")
printf("An hour of 0 returns full diurnal info with atm info\n")
printf("type = temperature_type (1 = kinetic, 2 = bolometer default)\n")
printf("latitude = latitude_bin (1-10, default is 1)\n")
printf("\n")
return(0)
}
if (HasValue(hour)==0) {
hour=2
}
if (HasValue(type)==0) {
type=2
}
if (HasValue(latitude)==0) {
latitude=1
}
if (HasValue(raw)==0) {
raw=0
}
# Read in bin52 file
data=load_bin5($1)
# Get its dimensions 1, 2, and 4 are always 24, 3, and 81 respectively
dims=dim(data[1][1])
# Make lsubs resampling (assuming 0 start and 4.5 ls steps)
et=13046146
for (i=1; i<=79; i+=1) {
et=cat(double(et),double(et[,,1]+(3600*24*8.58725*i)),z)
}
datals=lsubs(et)
datals[,,1]=0
ls=create(1,1,80,start=0,step=4.5,format=float)
if (hour>0) {
newdata=resample(translate(data[hour][type][latitude,2:,],y,z),datals,ls).data
newdata=cat(newdata,resample(translate(data[hour][2][latitude,2:,],y,z),datals,ls).data,y)
newdata=cat(newdata,resample(translate(data[hour][3][latitude,2:,],y,z),datals,ls).data,y)
newdata=cat(newdata,resample(translate(data[hour][4][latitude,2:,],y,z),datals,ls).data,y)
newdata=cat(newdata,resample(translate(data[hour][5][latitude,2:,],y,z),datals,ls).data,y)
}
if (hour==0) {
newdata=resample(translate(data[1][1][latitude,2:,],y,z),datals,ls).data
newdata=cat(newdata,resample(translate(data[1][2][latitude,2:,],y,z),datals,ls).data,y)
newdata=cat(newdata,resample(translate(data[1][3][latitude,2:,],y,z),datals,ls).data,y)
newdata=cat(newdata,resample(translate(data[1][4][latitude,2:,],y,z),datals,ls).data,y)
newdata=cat(newdata,resample(translate(data[1][5][latitude,2:,],y,z),datals,ls).data,y)
for (i=2; i<=24; i+=1) {
temp=resample(translate(data[i][1][latitude,2:,],y,z),datals,ls).data
temp=cat(temp,resample(translate(data[i][2][latitude,2:,],y,z),datals,ls).data,y)
temp=cat(temp,resample(translate(data[i][3][latitude,2:,],y,z),datals,ls).data,y)
temp=cat(temp,resample(translate(data[i][4][latitude,2:,],y,z),datals,ls).data,y)
temp=cat(temp,resample(translate(data[i][5][latitude,2:,],y,z),datals,ls).data,y)
newdata=newdata//temp
}}
out=struct(tk,tbol,tatm,down_vis,down_ir,ls)
out.tk=newdata[,1]
out.tbol=newdata[,2]
out.tatm=newdata[,3]
out.down_vis=newdata[,4]
out.down_ir=newdata[,5]
out.ls=ls
if (raw==1) {
out.raw=data
}
return(out)
}
define tmodel(lat,ti,lrat,lthick,tau,alb,dense,hcap,tstart,em,tfrost,ls,iter,time_step,nlayers,outsteps,planet,usage) {
if(HasValue(usage)==0) {
usage=0
}
if (usage==1) {
printf("\n")
printf("Davinci based ghetto thermal model\n")
printf("Inspired by Kieffer's krc code\n")
printf("\n")
printf("lat = latitude (default = 0)\n")
printf("ti = thermal inertia (default = 150) or 2x1x1 conductivity factors for variable cond.\n")
printf("lrat = layer scaling ratio in skin depths (default = 1.2)\n")
printf("lthick = intital layer thickness in skin depths (default = 0.15)\n")
printf("tau = visible atmospheric opacity (default = 0.40) 0 turns off atm.\n")
printf("alb = surface albedo (default = 0.20)\n")
printf("hcap = heat capacity (default = 830) or 3x1x1 heat capacity factors for variable hc\n")
printf("dense = bulk density (default = 1500)\n")
printf("tstart = starting temperature or 0 = start with day average equillibrium (default)\n")
printf("em = surface emissivity (default = 0.98)\n")
printf("tfrost = CO2 condensation temperature (default = 148)\n")
printf("ls = ls of predicted temperature (default = 90)\n")
printf("iter = length of time in mars/lunar years to iterate (default = 0.10)\n")
printf("time_step = number of iterations/day (default = 360)\n")
printf("outsteps = number of output steps for the final day (default = 24)\n")
printf("nlayers = number of layers to use (default = 30)\n")
printf("planet = mars (1) or the moon (2) (default = 1 (mars)\n")
printf("usage = print the help seeing as there are no arguments to check (default = 0)\n")
printf("\n")
printf("No frost control!\n")
printf("Returns 24 hour kinetic temperature curve at specified Ls\n")
printf("\n")
return(0)
}
##############################
#
# Set up variables
#
##############################
if (HasValue(lat)==0) {
lat=0.
}
if (HasValue(ti)==0) {
ti=150.
}
if (HasValue(lrat)==0) {
lrat=1.2
}
if (HasValue(lthick)==0) {
lthick=0.15
}
if (HasValue(tau)==0) {
tau=0.40
}
if (HasValue(alb)==0) {
alb=0.20
}
if (HasValue(hcap)==0) {
hcap=830.
}
if (HasValue(dense)==0) {
dense=1500.
}
if (HasValue(tstart)==0) {
tstart=0
}
if (HasValue(em)==0) {
em=0.98
}
if (HasValue(tfrost)==0) {
tfrost=148.
}
if (HasValue(ls)==0) {
ls=90.
}
if (HasValue(iter)==0) {
iter=0.1
}
if (HasValue(time_step)==0) {
time_step=360
}
if (HasValue(nlayers)==0) {
nlayers=30
}
if (HasValue(planet)==0) {
planet=1
}
if (HasValue(outsteps)==0) {
outsteps=24
}
time_step=int(time_step)
pi=3.1415927
e=2.7182818
sb=5.67051*10.^-8
# Set solar constant for 1 AU in W/m-2
S0=1368.
##############################
#
# Set up orbital parameters
#
##############################
# Read in orbital parameters file
if (planet==1) {
orb=read(sprintf("%s/orbital_parameters/mars.hdf",$DV_SCRIPT_FILES))
} else {
orb=read(sprintf("%s/orbital_parameters/lunar.hdf",$DV_SCRIPT_FILES))
}
# Set up iterative steps (time_step steps/day) ~4 minutes
# Truncate to the nearest full day for housekeeping reasons
delta_t=orb.day_length/time_step
tot_t=orb.year_length*iter
days=int((tot_t/delta_t)/time_step)
steps=int(days*time_step)
# Find start time dealing with wrap around issues
final_t=interp(orb.time_continuous,orb.ls_continuous,ls)
times=create(1,steps,1,start=(final_t-(steps-1)*delta_t),step=delta_t,format=float)
# if times are less than 0 find the equivalent time of year positive value
times [where (times<0)] = (times+orb.year_length)
# Find solar distance (AU), declination and hour angle (from midnight) for each time step (radians!)
solar_dec=interp(orb.planetocentric_sun_lat,orb.time,times)
solar_dist=interp(orb.solar_distance,orb.time,times)
hour_angle=create(1,steps,1,start=0.,step=(360./time_step),format=float)
hour_angle=(hour_angle%360.)/180*pi
# Get solar_incidence angle for each time step
cos_solar_inc = sin(solar_dec)*sind(lat)-cos(solar_dec)*cosd(lat)*cos(hour_angle)
cos_solar_inc [where (cos_solar_inc<=0)]=0
##############################
#
# Set up atm and energy inputs
#
##############################
# Get solar energy input at top of atmosphere (normal and scaled by incidence)
Sm = S0/(solar_dist)^2
Sm_inc = Sm*cos_solar_inc
# We're going to do a very crude atmospheric model because we are lazy right now
# We will transform tau, g, and w a la Delta Eddington
g = 0.70
w = 0.94
f=g^2
gp = g/(1+g)
wp = ((1-f)*w)/(1-w*f)
taup = (1-w*f)*tau
# First start with direct light
#Sm_inc_direct=Sm_inc*e^(-taup/cos_solar_inc)
#Sm_inc_direct [where (cos_solar_inc==0)]=0
#Sm_inc_direct [where (Sm_inc_direct>=1368.)]=0
#Sm_inc_direct_absorbed=Sm_inc_direct*(1-alb)
# Second we'll try to estimate scattered light
# Because of Delta Eddington transformation, things are somewhat more isotropic
#scattered_up=(Sm_inc-Sm_inc_direct)*wp*(1-gp)
#scattered_down=(Sm_inc-Sm_inc_direct)*wp*(gp)
# Transform based on solar incidence (no effect when sun is directly overhead)
# Transformed up = transformed down at inc = 90
#factor=acosd(cos_solar_inc)/180.
#scattered_up=((1.-factor)*scattered_up+factor*scattered_down)
#scattered_down=((1.-factor)*scattered_down+factor*scattered_up)
#Sm_inc_indirect_absorbed=(1-alb)*scattered_down
# Full up Delta-Eddington Approximation is now implemented using deddi
eddington=deddi(gp,wp,taup,acosd(cos_solar_inc),alb)
Sm_inc_direct_absorbed=Sm_inc*eddington.direct
Sm_inc_direct_absorbed [where (cos_solar_inc==0)]=0
Sm_inc_direct_absorbed [where (Sm_inc_direct_absorbed>=1368.)]=0
Sm_inc_indirect_absorbed=Sm*eddington.bottom*(1-alb)
# Now for IR stuff
# Using a scaling approximation from Jakosky and Haberle 1991 and
# mini-TES data rather than trying to model the atm directly
# First set up clear atm downwelling flux
IR_down_clear=create(1,9,1,start=0,step=(pi*0.25),format=float)
IR_down_clear=IR_down_clear//cat(20.,17.5,15.,17.,23.,26.5,27.,23.,20.,y)
#IR_down_clear=IR_down_clear//cat(26.01,24.73,23.57,23.93,26.66,28.99,28.61,27.30,26.01,y)
IR_down_clear=translate(resample(translate(IR_down_clear[2],y,z),translate(IR_down_clear[1],y,z),translate(hour_angle[,1:time_step],y,z)).data,z,y)
# Scale to dust opacity by empirical guestimate
IR_down=((1.+tau)*IR_down_clear) + (10.*tau)
IR_down_absorbed=IR_down*(em)
if (tau==0) {
IR_down_absorbed=IR_down_absorbed*0.
}
#IR_down_absorbed=IR_down_clear*(em)
# Create IR downwelling array to match other data
temp=IR_down_absorbed
for (i=2; i<=days; i+=1) {
temp=cat(temp,IR_down_absorbed,y)
}
IR_down_absorbed=temp
temp=0
sky_absorbed=IR_down_absorbed+Sm_inc_indirect_absorbed+Sm_inc_direct_absorbed
##############################
#
# Set up layers
#
##############################
# Get day average equillibrium temperature
if (tstart==0) {
energy_in=avg(sky_absorbed[,1:time_step])
tstart=(energy_in/(em*sb))^0.25
}
tsurf=create(1,(steps+1),1,start=0,step=0,format=float)
tsurf[,1]=tstart
# Get thermal conductivity, skin depth, layer thicknesses...
# Allow for variable thermal conductivity and heat capacity...
# Heat capacity is three terms 3 x 1 x 1 -- A*T^2 + B*T + C
# Thermal conductivity is two terms 2 x 1 x 1 -- A + B*T^3
if (dim(ti)[1]==1) {
cond = ti^2/(hcap*dense)
# tdiffuse = cond/(hcap*dense)
skind = ((cond*orb.day_length)/(pi*dense*hcap))^0.5
} else {
if (dim(hcap)[1]!=3) {
hcap=0.//0.//float(hcap[1])
}
hcap_factors = hcap
if (dim(hcap)[2]==1) {
hcap = clone(hcap_factors[1]*tstart^2+hcap_factors[2]*tstart+hcap_factors[3],y=nlayers)
} else {
hcap = hcap_factors[1]*tstart^2+hcap_factors[2]*tstart+hcap_factors[3]
}
cond_factors = ti
if (dim(ti)[2]==1) {
cond = clone(cond_factors[1]+cond_factors[2]*tstart^3,y=nlayers)
} else {
cond = cond_factors[1]+cond_factors[2]*tstart^3
}
skind = ((cond[1,]*orb.day_length)/(pi*dense*hcap[1,]))^0.5
}
layer=create(1,nlayers,1,start=0,step=1,format=float)
delta_T=create(1,nlayers,1,start=0,step=0,format=double)
layer_T=create(1,nlayers,1,start=tstart,step=0,format=double)
layer_T_max=layer_T
layer_T_min=layer_T
layer_thick_sd=lthick*(lrat^layer)
layer_thick=skind*layer_thick_sd
layer_center=layer_thick[,1]*0.5
layer_center_sd=layer_thick_sd[,1]*0.5
for (i=2; i<=nlayers; i+=1) {
layer_center=cat(layer_center,(float(sum(layer_thick[,1:(i-1)]))+layer_thick[,i]*0.5),y)
layer_center_sd=cat(layer_center_sd,(float(sum(layer_thick_sd[,1:(i-1)]))+layer_thick_sd[,i]*0.5),y)
}
##############################
#
# Do thermal diffusion stuff
#
##############################
# Start iterations to temperature
# Heat flow into top layer is sky_absorbed radiance-sb*em*T^4
# Heat flow into bottom layer from beneath is 0
for (i=1; i<=steps; i+=1) {
# Start with 2(B(i)rho(i)C(i))^-1
brc = double(2*(layer_thick*dense*hcap)^-1)
# Also do B(i)/k(i)
bk = double(layer_thick/cond)
delta_T[,2:(nlayers-1)] = (brc[,2:(nlayers-1)] * ((layer_T[,3:nlayers]-layer_T[,2:(nlayers-1)]) / (bk[,2:(nlayers-1)]+bk[,3:nlayers]) - (layer_T[,2:(nlayers-1)]-layer_T[,1:(nlayers-2)]) / (bk[,2:(nlayers-1)]+bk[,1:(nlayers-2)])))*delta_t
delta_T[,nlayers] = (brc[,nlayers]*-(layer_T[,nlayers]-layer_T[,(nlayers-1)]) / (bk[,nlayers]+bk[,(nlayers-1)]))*delta_t
delta_T[,1] = (brc[,1]*((layer_T[,2]-layer_T[,1]) / (bk[,2]+bk[,1]) - (0.5*cond[,1]/(0.5*layer_thick[,1]))*(layer_T[,1]-tsurf[,i])))*delta_t
layer_T = layer_T+delta_T
layer_T_max [where (layer_T>layer_T_max)]=layer_T
layer_T_min [where (layer_T
tsurf[,i+1] = tsurf[,i] + ((sky_absorbed[,i]+(cond[,1]/(0.5*layer_thick[,1]))*(layer_T[,1]-tsurf[,i])-em*sb*tsurf[,i]^4) / ((cond[,1]/(0.5*layer_thick[,1])) + 4*em*sb*(tsurf[,i]^3)))
if (dim(ti)[1]!=1) {
hcap=hcap_factors[1]*layer_T^2+hcap_factors[2]*layer_T+hcap_factors[3]
cond=cond_factors[1]+cond_factors[2]*layer_T^3
}
#xplot(layer_T,axis=y)
}
# Resample output surface temperatures to get final day by specified amounts
# Index 1 is first step after midnight!
tsurf=tsurf[,(steps-time_step+1):]
tsurf=tsurf[,2:]
IR_down_absorbed=IR_down_absorbed[,(steps-time_step+1):]
Sm_inc_indirect_absorbed=Sm_inc_indirect_absorbed[,(steps-time_step+1):]
Sm_inc_direct_absorbed=Sm_inc_direct_absorbed[,(steps-time_step+1):]
scale=create(1,time_step,1,start=0,step=(orb.day_length/time_step),format=float)
hour_scale=create(1,int(outsteps),1,start=0,step=(orb.day_length/outsteps),format=float)
tsurf=translate(resample(translate(tsurf,y,z),translate(scale,y,z),translate(hour_scale,y,z)).data,z,y)
Sm_inc_indirect_absorbed=translate(resample(translate(Sm_inc_indirect_absorbed,y,z),translate(scale,y,z),translate(hour_scale,y,z)).data,z,y)
Sm_inc_direct_absorbed=translate(resample(translate(Sm_inc_direct_absorbed,y,z),translate(scale,y,z),translate(hour_scale,y,z)).data,z,y)
IR_down_absorbed=translate(resample(translate(IR_down_absorbed,y,z),translate(scale,y,z),translate(hour_scale,y,z)).data,z,y)
tsurf=cat(tsurf[,2:],tsurf[,1],y)
out=struct(layer_thick,layer_center,delta_T,layer_T,tsurf)
out.layer_thick=layer_thick
out.layer_center=layer_center
out.skin_depth=skind
out.delta_T=delta_T
out.layer_T=layer_T
out.layer_T_min=layer_T_min
out.layer_T_max=layer_T_max
out.tsurf=tsurf
out.IR_flux=cat(IR_down_absorbed[,2:],IR_down_absorbed[,1],y)
out.scattered_flux=cat(Sm_inc_indirect_absorbed[,2:],Sm_inc_indirect_absorbed[,1],y)
out.direct_flux=cat(Sm_inc_direct_absorbed[,2:],Sm_inc_direct_absorbed[,1],y)
out.hour_scale=create(1,int(outsteps),1,start=24./outsteps,step=24./outsteps,format=float)
out.scale=scale
return(out)
}
define deddi() {
if ($ARGC==0) {
printf("\n")
printf("Delta Eddington two stream approximation.\n")
printf("$1 = Assymetry parameter (g)\n")
printf("$2 = Single scattering albedo (w)\n")
printf("$3 = Opacity (tau)\n")
printf("$4 = Incidence angle (degrees)\n")
printf("$5 = Surface albedo\n")
printf("\n")
return(0)
}
# First off lets transform the input parameters
# The Delta-Eddington approximation basically
# throws some of the scattered radiance into the
# direct sunlight component since the angle at which
# it was scattered is small enough as to make little difference
f=$1^2
g=$1/(1.+$1)
tau=(1.-$2*f)*$3
w=((1.-f)*$2)/(1.-$2*f)
A=$5
mu=cosd($4)
e=2.7182818
pi=3.1415927
# Set up some of the basic functions to keep things looking relatively simple
# some notation is inspired by Kieffer's dedding.f program
k = (3.*(1.-w)*(1.-w*g))^0.5
p = (3.*(1.-w)/(1.-w*g))^0.5
alph = 3.*(1./pi)*w*(mu^2)*(1.+g*(1.-w))/(4.*(1.-(k^2)*(mu^2)))
beta = 3.*(1./pi)*w*mu*(1.+3.*g*(1.-w)*mu^2)/(4.*(1.-(k^2)*(mu^2)))
v1 = 1.+(2./3.)*p
v2 = 1.-(2./3.)*p
v3 = alph+(2./3.)*beta
v4 = (1.-A-(1.+A)*p*(2./3.))*e^(-k*tau)
v5 = (1.-A+(1.+A)*p*(2./3.))*e^(k*tau)
v6 = ((1.-A)*alph-(1.+A)*beta*(2./3.)+A*mu*(1./pi))*e^(-tau/mu)
# Calculate determinant and solve for c2, c2
det=v1*v5-v4*v2
c1 = (v3*v5-v6*v2)/det
c2 = (v1*v6-v4*v3)/det
I0top = c1+c2-alph
I1top = p*(c1-c2)-beta
I0bot = c1*e^(-k*tau) + c2*e^(k*tau) - alph*e^(-tau/mu)
I1bot = p*(c1*e^(-k*tau) - c2*e^(k*tau)) - beta*e^(-tau/mu)
out=struct(top,bottom,direct)
out.top = pi*(I0top-(2./3.)*I1top)
out.bottom = pi*(I0bot+(2./3.)*I1bot)
out.direct = e^(-tau/mu)*(1-A)
return(out)
}
define mfp() {
#Calculates the Mean Free Path of Gas molecules in a soil
if($ARGC == 0) {
printf (" Returns a scalar \n")
printf (" Calculates the mean free path in m of various gas molecules \n")
printf (" $1 = a string indicating the nature of the gas, between double quotes \n")
printf (" Supported gases are CO2, N2, O2, CH4, He, H2, SO2, and H2O \n")
printf (" $2 = Temperature in K \n")
printf (" $3 = Pressure in Pa \n")
printf (" The output is the mean free path in m \n")
return(null)
}
gas = $1
T = float($2)
P = float($3)
#Loads the constants associated with each gas
if (gas == "CO2") {
diam = 4.357
} else {
if (gas == "N2") {
diam = 4.048
} else {
if (gas == "O2") {
diam = 3.906
} else {
if (gas == "CH4") {
diam = 4.286
} else {
if (gas == "He") {
diam = 2.651
} else {
if (gas == "H2") {
diam = 3.260
} else {
if (gas == "SO2") {
diam = 4.604
} else {
if (gas == "H2O") {
diam = 3.042
}}}}}}}}
#Mean free path calculation
mfp = 1.380650424E-23*T/(P*3.14159265*(diam*1.E-10)^2*sqrt(2.))
return(mfp)
}
define pore() {
#Calculate Pore sizes in soils
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the typical pore size in homogeneous samples \n")
printf (" $1 = the grain diameter in m \n")
printf (" $2 = the porosity, <1 \n")
printf (" MANY MORE PORE SIZE FUNCTIONS EXIST AND SHOULD BE INCLUDED \n")
return(null)
}
diam = $1
phi=$2
#Assumes that the pores are homogeneous spheres
#Add more pore size formulas here and references
#Grain size mixtures ?
pore = 0.5*diam*(phi/(1-phi))^0.33
return (pore)
}
define Kn(){
#Calculates a Knudsen number
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the Knudsen number in the pores of soils \n")
printf (" $1 = the gas composition, between double quotes \n")
printf (" Supported gases are CO2, N2, O2, CH4, He, H2, SO2, and H2O \n")
printf (" $2 = Temperature in K \n")
printf (" $3 = Pressure in Pa \n")
printf (" $4 = the grain diameter in m \n")
printf (" $5 = the porosity, <1 \n")
printf (" The pore size is calculated from pore() => spherical pores \n")
return(null)
}
mfp = mfp($1,float($2),float($3))
pore = pore($4,$5)
Knudsen = mfp/pore
return (Knudsen)
}
define k0() {
#Calculates a gas thermal conductivity in the continuum (high P) domain
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates a gas thermal conductivity in the continuum (high P) domain \n")
printf (" $1 = the gas composition, between double quotes \n")
printf (" Supported gases are CO2, N2, O2, CH4, He, H2, SO2, and H2O \n")
printf (" $2 = Temperature in K \n")
return(null)
}
gas = $1
T = $2
#Constants given by DIPPR 801
if (gas == "CO2") {
a = 3.69
b = -0.3838
c = 964.0
d = 1860000
} else {
if (gas == "N2") {
a = 0.00033143
b = 0.7722
c = 16.323
d = 373.72
} else {
if (gas == "O2") {
a = 0.00044994
b = 0.7456
c = 56.699
d = 0.0
} else {
if (gas == "CH4") {
a = 0.0000083983
b = 1.4268
c = -49.654
d = 0.0
} else {
if (gas == "He") {
a = 0.00226
b = 0.7305
c = -18.63
d = 440.0
} else {
if (gas == "H2") {
a = 0.002653
b = 0.7452
c = 12.0
d = 0.0
} else {
if (gas == "SO2") {
a = 0.00035686
b = 0.73812
c = 795.65
d = 0.0
} else {
if (gas == "H2O") {
a = 0.0000062041
b = 1.3973
c = 0.0
d = 0.0
}}}}}}}}
k0_1 = (a*T^b)/(1+c/T+d/T^2)
return(k0_1)
}
define sigmoid() {
# Calculates the T normalized gas conductivity in a pore
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the normalized gas conductivity from Kn \n")
printf (" $1 = the gas composition, between double quotes \n")
printf (" Supported gases are CO2, N2, O2, CH4, He, H2, SO2, and H2O \n")
printf (" $2 = Temperature in K \n")
printf (" $3 = Pressure in Pa \n")
printf (" $4 = grain diameter in m \n")
printf (" $5 = the porosity, <1 \n")
printf (" The pore size is calculated from pore() => spherical pores \n")
return(null)
}
Kn = Kn($1,$2,$3,$4,$5)
buffer = -(log10(1/Kn)-2.15)/0.55
buffer_2 = 1/(1+exp(buffer))
return (buffer_2)
}
define k_gas() {
# Calculates the gas conductivity in a pore
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the true gas conductivity from Kn \n")
printf (" $1 = the gas composition, between double quotes \n")
printf (" Supported gases are CO2, N2, O2, CH4, He, H2, SO2, and H2O \n")
printf (" $2 = Temperature in K \n")
printf (" $3 = Pressure in Pa \n")
printf (" $4 = grain diameter in m \n")
printf (" $5 = the porosity, <1 \n")
printf (" The pore size is calculated from pore() => spherical pores \n")
return(null)
}
k0 = k0($1,$2)
norm = sigmoid($1,$2,$3,$4,$5)
k_gas_1 = k0 * norm
return (k_gas_1)
}
define k_radsol() {
# Calculates the radiative conductivity between grains
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates A and B, solid and radiative contributions, from Watson 1964 \n")
printf (" $1 = the grain diameter in m \n")
printf (" $2 = Temperature in K \n")
return(null)
}
A = ($1/30 + 0.32E-5)*1E-5
B = 3E-8/($1)
k = A*($2)^3+B
return (k)
}
define k_SCP() {
# Calculates the conductivity of a SCP cell
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the bulk conductivity of a SCP cell \n")
printf (" SCP porosity is 47.6% \n")
printf (" The solid/solid contribution is not included \n")
printf (" The radiative contribution is not included \n")
printf (" $1 = gas phase conductivity, in SI \n")
printf (" $2 = solid phase conductivity, in SI \n")
return(null)
}
k_gas=float($1)
k_sol=float($2)
X = log10(k_gas/k_sol)
Y_buffer = 0.0008*X^5+0.0119*X^4+0.0501*X^3-0.0367*X^2+0.3799*X-0.3106
Delta = 0.4343*ln(k_sol)+0.3011
Y = Y_buffer+Delta
k = 10^Y
return(k)
}
define k_BCCP() {
# Calculates the conductivity of a BCCP cell
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the bulk conductivity of a SCP cell \n")
printf (" BCCP porosity is 25.9% \n")
printf (" The solid/solid contribution is not included \n")
printf (" The radiative contribution is not included \n")
printf (" $1 = gas phase conductivity, in SI \n")
printf (" $2 = solid phase conductivity, in SI \n")
return(null)
}
k_gas=float($1)
k_sol=float($2)
X = log10(k_gas/k_sol)
Y_buffer = 0.0008*X^5+0.0134*X^4+0.0687*X^3+0.022*X^2+0.2942*X-0.3024
Delta = 0.4343*ln(k_sol)+0.301
Y = Y_buffer+Delta
k = 10^Y
return(k)
}
define k_CCP() {
# Calculates the conductivity of a CCP cell
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the bulk conductivity of a SCP cell \n")
printf (" BCCP porosity is 32.0% \n")
printf (" The solid/solid contribution is not included \n")
printf (" The radiative contribution is not included \n")
printf (" $1 = gas phase conductivity, in SI \n")
printf (" $2 = solid phase conductivity, in SI \n")
return(null)
}
k_gas=float($1)
k_sol=float($2)
X = log10(k_gas/k_sol)
Y_buffer = 0.0012*X^5+0.02*X^4+0.0993*X^3+0.0514*X^2+0.2253*X-0.3061
Delta = 0.4343*ln(k_sol)+0.301
Y = Y_buffer+Delta
k = 10^Y
return(k)
}
define k_cell() {
# Calculates the bulk conductivity of a sample
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the bulk conductivity of a cell \n")
printf (" The solid/solid contribution is not included \n")
printf (" The radiative contribution is not included \n")
printf (" $1 = gas phase conductivity, in SI \n")
printf (" $2 = solid phase conductivity, in SI \n")
printf (" $3 = the porosity, <1 \n")
return(null)
}
k_gas = float($1)
k_sol = float($2)
phi = float($3)
k = k_SCP(k_gas,k_sol)+(phi-0.259)*(k_CCP(k_gas,k_sol)-k_SCP(k_gas,k_sol))/(0.476-0.259)
return(k)
}
define k_tot() {
# Calculates the bulk conductivity of a soil, all included
if($ARGC == 0){
printf (" Returns a scalar \n")
printf (" Calculates the bulk conductivity of a cell \n")
printf (" The solid/solid contribution is included \n")
printf (" The radiative contribution is included \n")
printf (" $1 = gas composition, between double quotes \n")
printf (" Supported gases are CO2, N2, O2, CH4, He, H2, SO2, and H2O \n")
printf (" $2 = Temperature in K \n")
printf (" $3 = Pressure in Pa \n")
printf (" $4 = grain diameter in m \n")
printf (" $5 = the porosity, <1 \n")
printf (" The pore size is calculated from pore() => spherical pores \n")
printf (" $6 = solid phase conductivity, in SI \n")
return(null)
}
diam = float($4)
phi = float($5)
gas = $1
T = float($2)
P = float($3)
k_sol = float($6)
pore = pore(diam,phi)
mfp = mfp(gas,T,P)
Kn = Kn(gas,T,P,diam,phi)
sigmoid = sigmoid(gas,T,P,diam,phi)
k0 = k0(gas,T)
k_gas = k_gas(gas,T,P,diam,phi)
k_radsol = k_radsol(diam,T)
k_SCP = k_SCP(k_gas,k_sol)
k_BCCP = k_BCCP(k_gas,k_sol)
k_CCP = k_CCP(k_gas,k_sol)
k_cell = k_cell(k_gas,k_sol,phi)
k = k_cell + k_radsol
return(k)
}
define Cp(){
# Calculates the specific heat of basalt from ~400K to 77K
if($ARGC == 0){
printf (" Calculates the specific heat of basaltic material, SI \n")
printf (" From Fujii and Osako, 1973 \n")
printf (" $1: Temperature in K \n")
return(null)
}
Cp = (0.0650+(0.394E-3)*$1-(0.452E3)/$1^2)*4184.
return(Cp)
}
define I2d(){
if($ARGC == 0){
printf (" Calculates the grain size in m from a Thermal inertia in SI \n")
printf (" From Piqueux and Christensen 2009 \n")
printf (" $1: Temperature in K \n")
printf (" $2: Pressure in Pa \n")
printf (" $3: Thermal Inertia in SI \n")
return(null)
}
T = $1
P = $2
I = $3
Poro_k = 0.45
Poro_soil = 0.30
Dens = (1-Poro_soil)*2650
Cp = Cp(T)
Table = create(2,1000,1,start=0.,step=0.,format=float)
Table[1,1,1] = 15E-6
Table[2,1,1] = sqrt(Dens*Cp*k_tot("CO2",T,P,Table[1,1,1],Poro_k,2.))
for(i=2;i<=1000;i+=1){
Table[1,i,1] = Table[1,i-1,1]+1.E-6
Table[2,i,1] = sqrt(Dens*Cp*k_tot("CO2",T,P,Table[1,i,1],Poro_k,2.))
if (Table[2,i,1] > I){
d = Table[1,i,1]
i = 2000
}
}
return(d)
}