プログラミング学習の記録 #23（C）：温度構造
ボックスモデルにおいて各ボックスで拡散係数が異なる場合の拡散現象を計算するプログラムコードを用いて、境界条件を与えた場合の温度構造を計算するプログラムのプログラムコードを以下に記録しておこうと思う。プログラミング言語は、C言語を用いている。
境界条件を与えた場合の温度構造
海洋の鉛直方向の温度構造を想定している。海洋表層温度$${T_{0} = 293.15 ~ \mathrm{K}}$$、第100層（深度600m程度）の温度$${T_{100} = 276.15 ~ \mathrm{K}}$$として、温度構造を計算する。以下の計算においては、約100年未満で終端状態に達した。
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// Input the numbers of layer and depth.
#include <time.h>
#define layer 101
#define depth 606.0
// Input the nunmbers of time steps and end time.
#define year_rate 0.0027397300/4.0000
#define year_max 100.0
// the function of diffusion.
double get_diffusion(int m,double y0,double y,double y1,double M0,double M,double M1,double z0,double z,double z1)
{
double diff;
diff = ((y0 + y) / 2.0) * (M0 - M) / pow((z - z0),2) + ((y1 + y) / 2.0) * (M1 - M) / pow((z1 - z),2) ;
return diff;
}
int main()
{
FILE *fp;
int m;
double t,dt;
double T[layer],diff_T[layer];
double T_surface,T_bottom;
double y[layer],yT[layer];
double z[layer],z_bottom[layer],dz[layer],z_surface;
// reseting time.
t = 0.0 ;
dt = year_rate ;
// depth.
z_surface = 0.0 ;
z_bottom[0] = depth / (1.0 * layer) ;
dz[0] = z_bottom[0] - z_surface ;
z[0] = (z_bottom[0] + z_surface) / 2.0 ;
for(m = 1 ; m < layer ; m = m + 1)
{
z_bottom[m] = z_bottom[m-1] + depth / (1.0 * layer) ;
dz[m] = z_bottom[m] - z_bottom[m-1] ;
z[m] = (z_bottom[m] + z_bottom[m-1]) / 2.0 ;
}
// boundary conditions.
T_surface = 273.15 + 20.0 ;
T_bottom = 273.15 + 3.0 ;
// initial numbers of T.
for(m = 0 ; m < layer ; m = m + 1)
{
T[m] = 273.15 + 10.0 ;
}
// Input initial conditions.
T[0] = T_surface ;
T[layer-1] = T_bottom ;
// Input conditions of diffusion coefficient.
for(m = 0 ; m < layer ; m = m + 1)
{
if(z[m] < 50.0) // surface mixed layer.
{
y[m] = 20000.0 ;
}else if(z[m] < 300.0) // thermal layer.
{
y[m] = 5000.0 ;
}else
{
y[m] = 2500.0 ;
}
}
// OK.
for(m = 0 ; m < layer ; m = m + 1)
{
yT[m] = y[m] * 1.0 ;
}
// Calculating time variation.
while(t < year_max)
{
for(m = 0 ; m < layer ; m = m + 1)
{
if(m == 0)
{
diff_T[m] = get_diffusion(m,yT[m],yT[m],yT[m+1],T_surface,T[m],T[m+1],2*z[m]-z[m+1],z[m],z[m+1]) ;
}else if(m == layer - 1)
{
diff_T[m] = get_diffusion(m,yT[m-1],yT[m],yT[m],T[m-1],T[m],T_bottom,z[m-1],z[m],2*z[m]-z[m-1]) ;
}else
{
diff_T[m] = get_diffusion(m,yT[m-1],yT[m],yT[m+1],T[m-1],T[m],T[m+1],z[m-1],z[m],z[m+1]) ;
}
}
for(m = 0 ; m < layer ; m = m + 1)
{
T[m] = T[m] + diff_T[m] * dt ;
}
t = t + dt ;
}
// outputting into CSV.
fp = fopen("test1_distribution.csv","w") ;
if(fp == NULL)
{
printf("File error.\n");
exit(1);
}
fprintf(fp,"layer,depth,T\n");
for(m = 0 ; m < layer ; m = m + 1)
{
fprintf(fp,"%2d,%5.1f,%10.5f\n",m,z[m],T[m]);
}
fclose(fp);
return 0;
}温暖化する場合の温度構造
上記のプログラムに、海洋温暖化の効果を追加する。海洋温暖化率（warming rate）は、$${0.61~\mathrm{K ~ yr^{-1}}}$$とした。
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// Input the numbers of layer and depth.
#define layer 101
#define depth 505.0
// Input the nunmbers of time steps and end time.
#define year_rate 0.0027397300/8.0000
#define year_max 200.0
// the function of diffusion.
double get_diffusion(int m,double y0,double y,double y1,double M0,double M,double M1,double z0,double z,double z1)
{
double diff;
diff = ((y0 + y) / 2.0) * (M0 - M) / pow((z - z0),2) + ((y1 + y) / 2.0) * (M1 - M) / pow((z1 - z),2) ;
return diff;
}
int main()
{
FILE *fp;
int m,i;
double t,dt;
double T[layer],diff_T[layer];
double T_surface,T_bottom;
double y[layer],yT[layer];
double z[layer],z_bottom[layer],dz[layer],z_surface;
double warming_rate;
// reseting time.
t = 0.0 ;
dt = year_rate ;
// depth.
z_surface = 0.0 ;
z_bottom[0] = depth / (1.0 * layer) ;
dz[0] = z_bottom[0] - z_surface ;
z[0] = (z_bottom[0] + z_surface) / 2.0 ;
for(m = 1 ; m < layer ; m = m + 1)
{
z_bottom[m] = z_bottom[m-1] + depth / (1.0 * layer) ;
dz[m] = z_bottom[m] - z_bottom[m-1] ;
z[m] = (z_bottom[m] + z_bottom[m-1]) / 2.0 ;
}
// Input boundary conditions.
T_surface = 273.15 + 20.0 ;
T_bottom = 273.15 + 3.0 ;
// OK.
// Input initial conditions.
for(m = 0 ; m < layer ; m = m + 1)
{
T[m] = 273.15 + 10.0 ;
}
T[0] = T_surface ;
T[layer-1] = T_bottom ;
// OK.
// Input conditions of diffusion coefficient.
for(m = 0 ; m < layer ; m = m + 1)
{
if(z[m] < 50.0) // surface mixed layer.
{
y[m] = 10000.0 ;
}else if(z[m] < 300.0) // thermal layer.
{
y[m] = 5000.0 ;
}else
{
y[m] = 2500.0 ;
}
}
for(m = 0 ; m < layer ; m = m + 1)
{
yT[m] = y[m] * 3.0 ;
}
// OK.
fp = fopen("test_process.csv","w") ;
if(fp == NULL)
{
printf("File error.\n");
exit(1);
}
fprintf(fp,"t,m=0,m=10,m=20,m=50,m=80,m=100\n");
i = 0 ;
// Calculating time variation.
while(t < year_max)
{
for(m = 0 ; m < layer ; m = m + 1)
{
if(m == 0)
{
diff_T[m] = get_diffusion(m,yT[m],yT[m],yT[m+1],T_surface,T[m],T[m+1],2*z[m]-z[m+1],z[m],z[m+1]) ;
}else if(m == layer - 1)
{
diff_T[m] = get_diffusion(m,yT[m-1],yT[m],yT[m],T[m-1],T[m],T_bottom,z[m-1],z[m],2*z[m]-z[m-1]) ;
}else
{
diff_T[m] = get_diffusion(m,yT[m-1],yT[m],yT[m+1],T[m-1],T[m],T[m+1],z[m-1],z[m],z[m+1]) ;
}
}
for(m = 0 ; m < layer ; m = m + 1)
{
T[m] = T[m] + diff_T[m] * dt ;
}
// outputting into CSV.
i = i + 1 ;
if(i * year_rate > year_max / 100.0)
{
fprintf(fp,"%10.4f, %10.5f, %10.5f, %10.5f, %10.5f, %10.5f, %10.5f\n",t,T[0],T[10],T[20],T[50],T[80],T[100]);
i = 0 ;
}
// grobal warming.
warming_rate = 0.00610 ; // K yr-1 [Japan Meteorological Agency,2024]
T_surface = T_surface + warming_rate * dt ;
t = t + dt ;
}
// outputting into CSV.
fp = fopen("test2_distribution.csv","w") ;
if(fp == NULL)
{
printf("File error.\n");
exit(1);
}
fprintf(fp,"layer,depth,T\n");
for(m = 0 ; m < layer ; m = m + 1)
{
fprintf(fp,"%2d,%5.1f,%10.5f\n",m,z[m],T[m]);
}
fclose(fp);
return 0;
}-----
動け！タイムライン
いいなと思ったら応援しよう！
