プログラミング学習の記録 #24（C）：水温の鉛直構造
前回の記事では、ボックスモデルにおいて、境界条件を与えた場合の温度構造を計算するプログラムのプログラムコードを記録した。しかし、海洋の水温構造を再現するためには、より実際的には水塊の鉛直運動や水平運動を考慮する必要があり、ボックスモデルによる単純な拡散現象の計算のみでは困難である。そこで、現実的な流体現象を扱うことを諦めつつ、水温の鉛直構造の観測データをボックスモデルで導入する場合に生データを参照するのではなく、仮想的な数値計算を用いることによって導入する場合のプログラムコードを以下に記録しておこうと思う。プログラミング言語は、C言語を用いている。
水温の鉛直構造
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// Input the numbers of layer and depth.
#define layer 101
#define depth 1010.0
// Input the nunmbers of time steps and end time.
#define year_rate 0.0027397300/2.0000
#define year_max 200.0
// function.
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,l;
double t,dt;
double T[layer],diff_T[layer];
double T_surface,T_bottom;
double yT[layer];
double gravity_T[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 ;
}
// Input boundary conditions.
T_surface = 273.0 + 20.0 ;
T_bottom = 273.0 + 2.0 ;
// OK.
// Input initial conditions.
for(m = 0 ; m < layer ; m = m + 1)
{
T[m] = 273.0 + 15.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] < 30.0) // surface mixed layer.
{
yT[m] = 30000.0 ;
}else if(z[m] < 50.0) // surface mixed layer.
{
yT[m] = 20000.0 ;
}else if(z[m] < 200.0) // thermal layer.
{
yT[m] = 15000.0 ;
}else if(z[m] < 400.0) // thermal layer.
{
yT[m] = 10000.0 ;
}else if(z[m] < 600.0) // thermal layer.
{
yT[m] = 3000.0 ;
}else if(z[m] < 800.0) // thermal layer.
{
yT[m] = 4000.0 ;
}else // deep layer.
{
yT[m] = 20000.0 ;
}
}
// OK.
// Calculating time variation.
while(t < year_max)
{
// temperature.
l = 0 ;
for(m = 0 ; m < layer ; m = m + 1)
{
if(m = 0)
{
T[m] = T_surface ;
}else if(m = layer - 1)
{
T[m] = T_bottom ;
}
if(l = 0)
{
if(z[m] > 300.0)
{
T[m] = 273.15 + 14.0 ;
l = 1 ;
}
}else if(l = 1)
{
if(z[m] > 700.0)
{
T[m] = 273.15 + 4.0 ;
l = 2 ;
}
}
}
// diffusion.
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("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;
}このプログラムを実行して得られた計算結果は、以下の図の通りである。
そして、気象庁Webサイトに掲載されている「日本近海の海面水温に見られる十年規模変動」より、「水温の鉛直プロファイルの例と主躍層」を以下に引用する。この図は、気象庁の海洋気象観測船で観測した北緯31度・東経144度における冬季の水温プロファイルと主躍層を表した図である。
これらの図を重ねると、以下のようになる。
おおよそ似たような形の図になった。この調整は、各ボックスにおける拡散係数や境界条件としての温度を調整することで行った。
-----
動け！タイムライン
いいなと思ったら応援しよう！
