いや、最終確認というよりきちんとモデルの中（コード）のパラメーターを確認するのは初めてかもしれない。
確認してわかったことは、PISCESは改良されつづけているようだが改良の変更点を明確にしたマニュアルが現時点では残念ながらないようだ。僕はAumontが２００６年に出した論文の補足資料としてウェブ上にあるマニュアル（ftp://ftp.agu.org/apend/gb/2005gb002591/）を参考に論文を書いていただけに残念だ。これ以後PISCES式やパラメーターは微調整されていろんなエンドユーザーに使われているんだがモデルを紹介するときはみな上記論文を引用している。これじゃblack boxを使っているのと同じじゃないか…。
っということで上記論文からNEMO3.4.1版まででの相違点・変更をパラメーターを周辺にまとめてみることに。
&nampisbio ! biological parameters
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
nrdttrc = 1 ! time step frequency for biology
wsbio = 2. ! POC sinking speed
xkmort = 2.E-7 ! half saturation constant for mortality
ferat3 = 10.E-6 ! Fe/C in zooplankton
wsbio2 = 30. ! Big particles sinking speed
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampislim ! parameters for nutrient limitations
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
conc0 = 1.E-6 ! Phosphate half saturation
conc1 = 8.E-6 ! Phosphate half saturation for diatoms
conc2 = 1.E-9 ! Iron half saturation for phyto
conc2m = 3.E-9 ! Max iron half saturation for phyto
conc3 = 3.E-9 ! Iron half saturation for diatoms
conc3m = 8.E-9 ! Maxi iron half saturation for diatoms
xsizedia = 1.E-6 ! Minimum size criteria for diatoms
xsizephy = 1.E-6 ! Minimum size criteria for phyto
concnnh4 = 1.E-7 ! NH4 half saturation for phyto
concdnh4 = 8.E-7 ! NH4 half saturation for diatoms
xksi1 = 2.E-6 ! half saturation constant for Si uptake
xksi2 = 333.E-8 ! half saturation constant for Si/C
xkdoc = 417.E-6 ! half-saturation constant of DOC remineralization
concfebac = 1.E-11 ! Half-saturation for Fe limitation of Bacteria
qnfelim = 7.E-6 ! Optimal quota of phyto
qdfelim = 7.E-6 ! Optimal quota of diatoms
caco3r = 16.E-2 ! mean rain ratio
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampisprod ! parameters for phytoplankton growth
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
pislope = 2. ! P-I slope
pislope2 = 2. ! P-I slope for diatoms
excret = 5.E-2 ! excretion ratio of phytoplankton
excret2 = 5.E-2 ! excretion ratio of diatoms
ln_newprod = .true. ! Enable new parame. of production (T/F)
bresp = 333.E-5 ! Basal respiration rate
chlcnm = 33.E-3 ! Minimum Chl/C in nanophytoplankton
chlcdm = 5.E-2 ! Minimum Chl/C in diatoms
chlcmin = 33.E-4 ! Maximum Chl/c in phytoplankton
！後のミニマムとマキシマム逆。
fecnm = 40.E-6 ! Maximum Fe/C in nanophytoplankton
fecdm = 40.E-6 ! Minimum Fe/C in diatoms
grosip = 151.E-3 ! mean Si/C ratio
これはマニュアルでは0.14の定数として記されている。
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampismort ! parameters for phytoplankton sinks
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
wchl = 1.E-3 ! quadratic mortality of phytoplankton
wchld = 2.E-2 ! maximum quadratic mortality of diatoms
mprat = 1.E-2 ! phytoplankton mortality rate
mprat2 = 1.E-2 ! Diatoms mortality rate
mpratm = 1.E-2 ! Phytoplankton minimum mortality rate
こいつは使われていないのになぜか存在する。
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampismes ! parameters for mesozooplankton
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
part2 = 75.E-2 ! part of calcite not dissolved in mesozoo guts
grazrat2 = 7.E-1 ! maximal mesozoo grazing rate
resrat2 = 5.E-3 ! exsudation rate of mesozooplankton
mzrat2 = 3.E-2 ! mesozooplankton mortality rate
xprefc = 1. ! zoo preference for phyto
xprefp = 3.E-1 ! zoo preference for POC
xprefz = 1. ! zoo preference for zoo
xprefpoc = 3.E-1 ! zoo preference for poc
このへんのインデックス要注意。今まで「2」はmesozooplanktonを表すために使われていたのにこれらのパラメーターのときだけ逆になっとる（この場合xpref2cはmicrozooplanktonに関連している）。
xthresh2zoo = 1.E-8 ! zoo feeding threshold for mesozooplankton
xthresh2dia = 1.E-8 ! diatoms feeding threshold for mesozooplankton
xthresh2phy = 1.E-8 ! nanophyto feeding threshold for mesozooplankton
xthresh2poc = 1.E-8 ! poc feeding threshold for mesozooplankton
xthresh2 = 2.E-7 ! Food threshold for grazing
xkgraz2 = 20.E-6 ! half sturation constant for meso grazing
epsher2 = 3.E-1 ! Efficicency of Mesozoo growth
sigma2 = 6.E-1 ! Fraction of mesozoo excretion as DOM
unass2 = 3.E-1 ! non assimilated fraction of P by mesozoo
grazflux = 2.E3 ! flux-feeding rate
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampiszoo ! parameters for microzooplankton
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
part = 5.E-1 ! part of calcite not dissolved in microzoo gutsa
grazrat = 3.0 ! maximal zoo grazing rate
resrat = 3.E-2 ! exsudation rate of zooplankton
mzrat = 1.E-3 ! zooplankton mortality rate
捕食されるのはmesozooplanktonのみとマニュアルには書いてあったのにコードではこれも計算されてるっぽい…。そして変動はzooplankton量に対しsquare dependencyのはずなのにコードを見る限りでは二乗されてないと思うんだが。。。
xpref2c = 1.E-1 ! Microzoo preference for POM
xpref2p = 1. ! Microzoo preference for Nanophyto
xpref2d = 5.E-1 ! Microzoo preference for Diatoms
xthreshdia = 1.E-8 ! Diatoms feeding threshold for microzooplankton
xthreshphy = 1.E-8 ! Nanophyto feeding threshold for microzooplankton
xthreshpoc = 1.E-8 ! POC feeding threshold for microzooplankton
xthresh = 2.E-7 ! Food threshold for feeding
xkgraz = 20.E-6 ! half sturation constant for grazing
epsher = 3.E-1 ! Efficiency of microzoo growth
sigma1 = 6.E-1 ! Fraction of microzoo excretion as DOM
unass = 3.E-1 ! non assimilated fraction of phyto by zoo
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampisrem ! parameters for remineralization
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
xremik = 25.E-2 ! remineralization rate of DOC
xremip = 25.E-3 ! remineralisation rate of POC
これと一致するパラメータはマニュアルの式では記述されてない。
nitrif = 5.E-2 ! NH4 nitrification rate
xsirem = 3.E-3 ! remineralization rate of Si
xsiremlab = 25.E-3 ! fast remineralization rate of Si
xsilab = 31.E-2 ! Fraction of labile biogenic silica
xlam1 = 5.E-3 ! scavenging rate of Iron
oxymin = 1.E-6 ! Half-saturation constant for anoxia
ligand = 6.E-10 ! Ligands concentration
/
!'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
&nampiscal ! parameters for Calcite chemistry
!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
kdca = 6. ! calcite dissolution rate constant (1/time)
nca = 1. ! order of dissolution reaction (dimensionless)
2013年6月27日木曜日
2013年6月24日月曜日
Sensitivity Analysisだん
I ended up rewriting the script (sensitivityanalysis.sh) to resolve the precision errors involving floats with ACEnet clusters' Bash. It is probably better to avoid using decimals in Bash. I have adjusted the script and namelist_pisces accordingly (see the notes on the script for details).
Here is my script and it should be the final version:
# Sensitivity Analysis of PISCES model parameters
# INSTRUCTIONS
#
# The value before the exponent (E) in namelist_pisces for all PISCES parameters has to be an integer. If it involves a number after the decimal place, you have to rewrite it using E. This is to avoid floating precision errors involved in Bash.
# For Example, If a1 = 0.034E-5, rewrite it as a1 = 34.E-7. If a2 = 99.9, rewrite it as a2 = 999.E-1. etc...
#
#END of INSTRUCTIONS
# Start of Input. Modify below #
# set of names of parameters to analyze.
set_parnam="wsbio xkmort ferat3 wsbio2 conc0 conc1 conc2 conc2m conc3 conc3m xsizedia xsizephy concnnh4 concdnh4 xksi1 xksi2 xkdoc concfebac qnfelim qdfelim caco3r pislope pislope2 excret excret2 bresp chlcnm chlcdm chlcmin fecnm fecdm grosip wchl wchld mprat mprat2 mpratm part2 grazrat2 resrat2 mzrat2 xprefc xprefp xprefz xprefpoc xthresh2zoo xthresh2dia xthresh2phy xthresh2poc xthresh2 xkgraz2 epsher2 sigma2 unass2 grazflux part grazrat resrat mzrat xpref2c xpref2p xpref2d xthreshdia xthreshphy xthreshpoc xthresh xkgraz epsher sigma1 unass xremik xremip nitrif xsirem xsiremlab xsilab xlam1 oxymin ligand kdca nca"
# set of the corresponding values to the above set (excluding exponents (E) if there are any).
set_stdval="2. 2. 10. 30. 1. 8. 1. 3. 3. 8. 1. 1. 1. 8. 2. 333. 417. 1. 7. 7. 16. 2. 2. 5. 5. 333. 33. 5. 33. 40. 40. 151. 1. 2. 1. 1. 1. 75. 7. 5. 3. 1. 3. 1. 3. 1. 1. 1. 1. 2. 20. 3. 6. 3. 2. 5. 3. 3. 1. 1. 1. 5. 1. 1. 1. 2. 20. 3. 6. 3. 25. 25. 5. 3. 25. 31. 5. 1. 6. 6. 1."
# End of Input. Do NOT modify below #
parnam=($set_parnam)
stdval=($set_stdval)
count=${#parnam[@]}
floats=1 # Number of float precision
for it in `seq 1 $count`
do
para=0
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
for newnum in $(seq `echo "scale=${floats};${stdval[$it-1]}/2" | bc` `echo "scales=${floats};${stdval[$it-1]}/10" | bc` `echo "scales=${floats};${stdval[$it-1]}/2+${stdval[$it-1]}" | bc`)
do
sed -i -r "s/${parnam[$it-1]}\s+=\s+[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${newnum}/" namelist_pisces
./opa
para=$((${para}+1))
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
mv LABSEA_1m*_ptrc_T.nc ../MY_SRC/SA/${parnam[$it-1]}${para}.nc
rm LABSEA_*restart*.nc
done
echo "end of run for ${parnam[$it-1]}" >> SA_stat.txt
sed -i -r "s/${parnam[$it-1]} += +[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${stdval[$it-1]}/" namelist_pisces
このスクリプトはquequeで走らせるのではなくhead nodeで走らせないと10倍の時間がかかる。１Dの場合head nodeで走らせたとしても２０％ぐらいのCPUしか使用されないので
nohup sensitivityanalysis.sh > SA.out
というふうにすることに。
またいつか使うときが来るといいけど。
Here is my script and it should be the final version:
# Sensitivity Analysis of PISCES model parameters
# INSTRUCTIONS
#
# The value before the exponent (E) in namelist_pisces for all PISCES parameters has to be an integer. If it involves a number after the decimal place, you have to rewrite it using E. This is to avoid floating precision errors involved in Bash.
# For Example, If a1 = 0.034E-5, rewrite it as a1 = 34.E-7. If a2 = 99.9, rewrite it as a2 = 999.E-1. etc...
#
#END of INSTRUCTIONS
# Start of Input. Modify below #
# set of names of parameters to analyze.
set_parnam="wsbio xkmort ferat3 wsbio2 conc0 conc1 conc2 conc2m conc3 conc3m xsizedia xsizephy concnnh4 concdnh4 xksi1 xksi2 xkdoc concfebac qnfelim qdfelim caco3r pislope pislope2 excret excret2 bresp chlcnm chlcdm chlcmin fecnm fecdm grosip wchl wchld mprat mprat2 mpratm part2 grazrat2 resrat2 mzrat2 xprefc xprefp xprefz xprefpoc xthresh2zoo xthresh2dia xthresh2phy xthresh2poc xthresh2 xkgraz2 epsher2 sigma2 unass2 grazflux part grazrat resrat mzrat xpref2c xpref2p xpref2d xthreshdia xthreshphy xthreshpoc xthresh xkgraz epsher sigma1 unass xremik xremip nitrif xsirem xsiremlab xsilab xlam1 oxymin ligand kdca nca"
# set of the corresponding values to the above set (excluding exponents (E) if there are any).
set_stdval="2. 2. 10. 30. 1. 8. 1. 3. 3. 8. 1. 1. 1. 8. 2. 333. 417. 1. 7. 7. 16. 2. 2. 5. 5. 333. 33. 5. 33. 40. 40. 151. 1. 2. 1. 1. 1. 75. 7. 5. 3. 1. 3. 1. 3. 1. 1. 1. 1. 2. 20. 3. 6. 3. 2. 5. 3. 3. 1. 1. 1. 5. 1. 1. 1. 2. 20. 3. 6. 3. 25. 25. 5. 3. 25. 31. 5. 1. 6. 6. 1."
# End of Input. Do NOT modify below #
parnam=($set_parnam)
stdval=($set_stdval)
count=${#parnam[@]}
floats=1 # Number of float precision
for it in `seq 1 $count`
do
para=0
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
for newnum in $(seq `echo "scale=${floats};${stdval[$it-1]}/2" | bc` `echo "scales=${floats};${stdval[$it-1]}/10" | bc` `echo "scales=${floats};${stdval[$it-1]}/2+${stdval[$it-1]}" | bc`)
do
sed -i -r "s/${parnam[$it-1]}\s+=\s+[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${newnum}/" namelist_pisces
./opa
para=$((${para}+1))
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
mv LABSEA_1m*_ptrc_T.nc ../MY_SRC/SA/${parnam[$it-1]}${para}.nc
rm LABSEA_*restart*.nc
done
echo "end of run for ${parnam[$it-1]}" >> SA_stat.txt
sed -i -r "s/${parnam[$it-1]} += +[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${stdval[$it-1]}/" namelist_pisces
このスクリプトはquequeで走らせるのではなくhead nodeで走らせないと10倍の時間がかかる。１Dの場合head nodeで走らせたとしても２０％ぐらいのCPUしか使用されないので
nohup sensitivityanalysis.sh > SA.out
というふうにすることに。
またいつか使うときが来るといいけど。
2013年6月23日日曜日
同じshellなのに？
sensitivityanalysis.shが一晩走らせたら終わっていたのでMatlabで結果を解析してみると変なエラーがあることに気づいた。それは、二つのパラメータ(excret, excret2)のときだけ１１回されるはずのシミュレーションが１０回しかされていなかった。この二つが共通して他のパラメータと違う点はデフォルト値が0.05ということ。
そしていろいろ調べた結果：
for it in `seq 0.025 0.005 0.075`; do echo "$it"; done
すると不思議なことに双方で結果に違いがあった。
Macbook:
0.025
0.03
0.035
0.04
0.045
0.05
0.055
0.06
0.065
0.07
0.075
Placentia:
0.025
0.03
0.035
0.04
0.045
0.05
0.055
0.06
0.065
0.07
もちろん正しいのは前者。なぜか後者では最後のiterationがされていない。これを他の数字で検証してみる。
for it in `seq 0.025 0.005 0.085`; do echo "$it"; done
すると、
0.025
0.03
0.035
0.04
0.045
0.05
0.055
0.06
0.065
0.07
0.075
0.08
0.085
という解答が両者から得られた。
0.075はマジックナンバーなのか？？？
他のユーザーに同じことしてもらってもし同じ解答が得られたらACEnetに聞いてみよう。
さらに調べてみると他にもxsiremlab, nitrif, xremip, xpref2c, chlcdmも１１回目のイタレーションは行われていなかった。これらのデフォルト値はそれぞれ、0.025, 0.05, 0.025, 0.1, 0.05である。これらの数字の規則性は何か。
もしACEnetが何もできないというのであればifを使ってスクリプトをtweakしないといけないな。
さらに調べてみると他にもxsiremlab, nitrif, xremip, xpref2c, chlcdmも１１回目のイタレーションは行われていなかった。これらのデフォルト値はそれぞれ、0.025, 0.05, 0.025, 0.1, 0.05である。これらの数字の規則性は何か。
もしACEnetが何もできないというのであればifを使ってスクリプトをtweakしないといけないな。
2013年6月22日土曜日
for do done (shell scripting)
「seq: zero increment」再び！
少数の計算でまた不都合が発生してしていた。今回は「0.05×0.1=0」となってしまう現象。scale=6だから0.005は表示されるはずだと思ったんだがなぁ。
「bc -l」で解決！
ネットで検索して解決策が見つかった。そしてそれはあまりにもシンプルすぎた。
bcに「-l」のオプションを付け加えるだけでかけ算も割り算も支障なくできるようになった。このオプションデフォルトでいるだろって感じだが、bcを使うときはつけるように心がけたい。
「sensitivityanalysis.sh」最終的なスクリプト！
またいつか同じような計算をするかもしれないという将来的なことも考えてPISCESモデルパラメータの各名称＆デフォルト値をリストダウンしてスクリプトを作成しておくことにする。これはバージョン3.4.1なので他ではirrelevantな場合多いにあり。
# Start of Input. Modify below #
# set of names of parameters to analyze.
set_parnam="wsbio xkmort ferat3 wsbio2 conc0 conc1 conc2 conc2m conc3 conc3m xsizedia xsizephy concnnh4 concdnh4 xksi1 xksi2 xkdoc concfebac qnfelim qdfelim caco3r pislope pislope2 excret excret2 bresp chlcnm chlcdm chlcmin fecnm fecdm grosip wchl wchld mprat mprat2 mpratm part2 grazrat2 resrat2 mzrat2 xprefc xprefp xprefz xprefpoc xthresh2zoo xthresh2dia xthresh2phy xthresh2poc xthresh2 xkgraz2 epsher2 sigma2 unass2 grazflux part grazrat resrat mzrat xpref2c xpref2p xpref2d xthreshdia xthreshphy xthreshpoc xthresh xkgraz epsher sigma1 unass xremik xremip nitrif xsirem xsiremlab xsilab xlam1 oxymin ligand kdca nca"
# set of standard PISCES values of the above parameters (omitting Exponents if there are any).
floats=6 # number of decimal places. ⇦⇦⇦これは別に必要なかったが簡素化するため。
# End of Input. Do NOT modify below #
parnam=($set_parnam)
stdval=($set_stdval)
count=${#parnam[@]}
for it in `seq 1 $count`
do
para=0
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
for newnum in $(seq `echo "scale=${floats};${stdval[$it-1]}/2" | bc -l` `echo "scales=${floats};${stdval[$it-1]}/10" | bc -l` `echo "scales=${floats};${stdval[$it-1]}/2+${stdval[$it-1]}" | bc -l`)
do
sed -i -r "s/${parnam[$it-1]}\s+=\s+[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${newnum}/" namelist_pisces
./opa
para=$((${para}+1))
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
mv LABSEA_1m*_ptrc_T.nc ../MY_SRC/SA/${parnam[$it-1]}${para}.nc
rm LABSEA_*restart*.nc
done
echo "end of run for ${parnam[$it-1]}" >> SA_stat.txt
sed -i -r "s/${parnam[$it-1]} += +[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${stdval[$it-1]}/" namelist_pisces
done
2013年6月21日金曜日
sensitivity analysis続き
ACEnet: Placentia
我がメモリアル大のスパコンクラスターPlacentiaのhead nodeとshort.qで計算時間に違いが。前者だと５分で終わる計算が後者だと１時間もかかった。違いは何か。メモリ(
RAM)量？とりあえずhead nodeのときは容量がわからない。short.qではデフォルト値が２GB。
short.qで８GBメモリリクエストすると５０分で計算できた。大して変わらない。Placentiaはユーザー数が多いせいか本来より長い時間かかるような気がする。
PISCESのパラメーターが８０ちかくあることに気づいたのでループを新しくスクリプトに追加することにした。以下がコード：
# Start of Input. Modify below #
set_parnam="conc3 conc3m" # set of names of parameters to analyze.
set_stdval="3.0 8.0" # set of standard PISCES values of the above parameters.
floats=6 # number of decimal places.
# End of Input. Do NOT modify below #
parnam=($set_parnam)
stdval=($set_stdval)
count=${#parnam[@]}
for it in `seq 1 $count`
do
para=0
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
for newnum in $(seq `echo "scale=${floats};${stdval[$it-1]}*0.5" | bc` `echo "scales=${floats};${stdval[$it-1]}*0.1" | bc` `echo "scales=${floats};${stdval[$it-1]}*0.5+${stdval[$it-1]}" | bc`)
do
sed -i -r "s/${parnam[$it-1]}\s+=\s+[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${newnum}/" namelist_pisces
./opa
para=$((${para}+1))
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
mv LABSEA_1m*_ptrc_T.nc ../MY_SRC/SA/${parnam[$it-1]}${para}.nc
rm LABSEA_*restart*.nc
done
echo "end of run for ${parnam[$it-1]}" >> SA_stat.txt
sed -i -r "s/${parnam[$it-1]} += +[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${stdval[$it-1]}/" namelist_pisces
done
ちょっと手こずったのは青色でハイライトされているところ。少数計算はechoとbcコマンドを使ってできることは知ったが割り算で答えが１以下になる場合は０となってしまう問題があることに気づいた。例えば:
echo "1/2" | bc
と入力しても0.5という答えにはならず「０」となってしまう。
どうしたら良いか考えた結果、割り算をかけ算に直すことにした。例えばさっきの場合だと１÷２ではなくて１× 0.5と置き換えることにした。割り算は極力避けた方がいいかもな。たぶん解決策はあるんだろうけど。
このスクリプトをhead nodeで試してみると２つのパラメータで合計２０分で計算できた。ひとつあたり１０分計算。早い。しかもcpuもそんなに使わない（２０％程度）なのでたぶん走らせといても何も言われないだろう。nohupで夜行列車モードだな。
short.qで８GBメモリリクエストすると５０分で計算できた。大して変わらない。Placentiaはユーザー数が多いせいか本来より長い時間かかるような気がする。
PISCESのパラメーターが８０ちかくあることに気づいたのでループを新しくスクリプトに追加することにした。以下がコード：
# Start of Input. Modify below #
set_parnam="conc3 conc3m" # set of names of parameters to analyze.
set_stdval="3.0 8.0" # set of standard PISCES values of the above parameters.
floats=6 # number of decimal places.
# End of Input. Do NOT modify below #
parnam=($set_parnam)
stdval=($set_stdval)
count=${#parnam[@]}
for it in `seq 1 $count`
do
para=0
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
for newnum in $(seq `echo "scale=${floats};${stdval[$it-1]}*0.5" | bc` `echo "scales=${floats};${stdval[$it-1]}*0.1" | bc` `echo "scales=${floats};${stdval[$it-1]}*0.5+${stdval[$it-1]}" | bc`)
do
sed -i -r "s/${parnam[$it-1]}\s+=\s+[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${newnum}/" namelist_pisces
./opa
para=$((${para}+1))
echo $para >> SA_stat.txt
grep "${parnam[$it-1]} " namelist_pisces >> SA_stat.txt
mv LABSEA_1m*_ptrc_T.nc ../MY_SRC/SA/${parnam[$it-1]}${para}.nc
rm LABSEA_*restart*.nc
done
echo "end of run for ${parnam[$it-1]}" >> SA_stat.txt
sed -i -r "s/${parnam[$it-1]} += +[0-9]+\.*[0-9]*/${parnam[$it-1]} = ${stdval[$it-1]}/" namelist_pisces
done
ちょっと手こずったのは青色でハイライトされているところ。少数計算はechoとbcコマンドを使ってできることは知ったが割り算で答えが１以下になる場合は０となってしまう問題があることに気づいた。例えば:
echo "1/2" | bc
と入力しても0.5という答えにはならず「０」となってしまう。
どうしたら良いか考えた結果、割り算をかけ算に直すことにした。例えばさっきの場合だと１÷２ではなくて１× 0.5と置き換えることにした。割り算は極力避けた方がいいかもな。たぶん解決策はあるんだろうけど。
このスクリプトをhead nodeで試してみると２つのパラメータで合計２０分で計算できた。ひとつあたり１０分計算。早い。しかもcpuもそんなに使わない（２０％程度）なのでたぶん走らせといても何も言われないだろう。nohupで夜行列車モードだな。
登録:
投稿 (Atom)