diff --git a/csswm_config.txt b/csswm_config.txt index beb79c8..952d977 100644 --- a/csswm_config.txt +++ b/csswm_config.txt @@ -1,9 +1,10 @@ -CSSWM_OUTPUTPATH=/data/Aaron/TMIF_CSSWM/EquatorialWaveTest +CSSWM_OUTPUTPATH=/data/Aaron/TMIF_CSSWM/Gill/csswm CSSWM_DT=200 -CSSWM_TIMEEND=1200000 -CSSWM_OUTPUTSTEP=1 +CSSWM_TIMEEND=4320000 +CSSWM_OUTPUTSTEP=500 CSSWM_GRAVITY=0.2391 CSSWM_DIFFUSION_KX=200000. CSSWM_DIFFUSION_KY=200000. CSSWM_DIFFUSION_TS=0.06 -CSSWM_ADDFORCING_TIME=60000 \ No newline at end of file +CSSWM_ADDFORCING_TIME=4320000 +CSSWM_H_NUDGE_TIME=86400 \ No newline at end of file diff --git a/src/bp_wind.cpp b/src/bp_wind.cpp index b853846..5e61b26 100644 --- a/src/bp_wind.cpp +++ b/src/bp_wind.cpp @@ -368,7 +368,7 @@ void CSSWM::BP_wind_interpolation(CSSWM &model) { #endif } - +/* void CSSWM::BP_wind_interpolation2(CSSWM &model) { double B, A1, A2, V1, V2, V3, V4; @@ -507,3 +507,165 @@ void CSSWM::BP_wind_interpolation2(CSSWM &model) { } #endif } +*/ + + +void CSSWM::BP_wind_interpolation2(CSSWM &model) { + double B, A1, A2, V1, V2, V3, V4; + + int p1, p2, i1, j1, i2, j2, reversed, lonlat; + double uIP, vIP; + int I1, I2_1, I2_2, J1, J2_1, J2_2; + for (int pp = 0; pp < 24; pp++) { + #if defined(SecondOrderSpace) + p1 = match[pp][0], p2 = match[pp][1], i1 = match[pp][2], j1 = match[pp][3], i2 = match[pp][4], j2 = match[pp][5], reversed = match[pp][6], lonlat = match[pp][7]; + #elif defined(FourthOrderSpace) + p1 = match_ouTTer[pp][0], p2 = match_ouTTer[pp][1], i1 = match_ouTTer[pp][2], j1 = match_ouTTer[pp][3], i2 = match_ouTTer[pp][4], j2 = match_ouTTer[pp][5], reversed = match_ouTTer[pp][6], lonlat = match_ouTTer[pp][7]; + #endif + for (int idx = 0; idx < NX; idx++) { + I1 = i1 == -1 ? idx : i1, J1 = j1 == -1 ? idx : j1; + #if defined(SecondOrderSpace) + I2_1 = i2 == -1 ? reversed ? checkIP[NX-1-idx][0] : checkIP[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP[NY-1-idx][0] : checkIP[idx][0] : j2; + I2_2 = i2 == -1 ? reversed ? checkIP[NX-1-idx][1] : checkIP[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP[NY-1-idx][1] : checkIP[idx][1] : j2; + #elif defined(FourthOrderSpace) + I2_1 = i2 == -1 ? reversed ? checkIP_ouTTer[NX-1-idx][0] : checkIP_ouTTer[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP_ouTTer[NY-1-idx][0] : checkIP_ouTTer[idx][0] : j2; + I2_2 = i2 == -1 ? reversed ? checkIP_ouTTer[NX-1-idx][1] : checkIP_ouTTer[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP_ouTTer[NY-1-idx][1] : checkIP_ouTTer[idx][1] : j2; + #endif + if (lonlat == 0) { + B = model.csswm[p1].lat[I1][J1]; + A1 = model.csswm[p2].lat[I2_1][J2_1], A2 = model.csswm[p2].lat[I2_2][J2_2]; + // V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2]; + // V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2]; + + V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1]; + V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1]; + V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2]; + V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2]; + } + else { + B = model.csswm[p1].lon[I1][J1]; + A1 = model.csswm[p2].lon[I2_1][J2_1], A2 = model.csswm[p2].lon[I2_2][J2_2]; + // V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2]; + // V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2]; + + V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1]; + V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1]; + V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2]; + V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2]; + + if (A1 > A2 && (p1 == 0 || p2 == 0)) A2 += 2 * M_PI; + if (A1 > B && B < A2) B += 2 * M_PI; + } + + uIP = model.interpolate(A1, A2, V1, V2, B); + vIP = model.interpolate(A1, A2, V3, V4, B); + + // TODO:transfer to from contravariant A to covariant B + if (i1 == -1) { + if (j1 == 0) { + #if defined(SecondOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_D[I1][0] * uIP + model.csswm[p1].IP1_D[I1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_D[I1][2] * uIP + model.csswm[p1].IP1_D[I1][3] * vIP; + #elif defined(FourthOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_D[I1][0] * uIP + model.csswm[p1].IP_ouTTer_D[I1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_D[I1][2] * uIP + model.csswm[p1].IP_ouTTer_D[I1][3] * vIP; + #endif + } + else if (j1 == NY-1) { + #if defined(SecondOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_U[I1][0] * uIP + model.csswm[p1].IP1_U[I1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_U[I1][2] * uIP + model.csswm[p1].IP1_U[I1][3] * vIP; + #elif defined(FourthOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_U[I1][0] * uIP + model.csswm[p1].IP_ouTTer_U[I1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_U[I1][2] * uIP + model.csswm[p1].IP_ouTTer_U[I1][3] * vIP; + #endif + } + } + else if (j1 == -1) { + if (i1 == 0) { + #if defined(SecondOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_L[J1][0] * uIP + model.csswm[p1].IP1_L[J1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_L[J1][2] * uIP + model.csswm[p1].IP1_L[J1][3] * vIP; + #elif defined(FourthOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_L[J1][0] * uIP + model.csswm[p1].IP_ouTTer_L[J1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_L[J1][2] * uIP + model.csswm[p1].IP_ouTTer_L[J1][3] * vIP; + #endif + } + else if (i1 == NX-1) { + #if defined(SecondOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP1_R[J1][0] * uIP + model.csswm[p1].IP1_R[J1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP1_R[J1][2] * uIP + model.csswm[p1].IP1_R[J1][3] * vIP; + #elif defined(FourthOrderSpace) + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTTer_R[J1][0] * uIP + model.csswm[p1].IP_ouTTer_R[J1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTTer_R[J1][2] * uIP + model.csswm[p1].IP_ouTTer_R[J1][3] * vIP; + #endif + } + } + } + } + + #if defined(FourthOrderSpace) + for (int pp = 0; pp < 24; pp++) { + p1 = match_ouTer[pp][0], p2 = match_ouTer[pp][1], i1 = match_ouTer[pp][2], j1 = match_ouTer[pp][3], i2 = match_ouTer[pp][4], j2 = match_ouTer[pp][5], reversed = match_ouTer[pp][6], lonlat = match_ouTer[pp][7]; + for (int idx = 0; idx < NX; idx++) { + if (lonlat == 0) { + I1 = i1 == -1 ? idx : i1, J1 = j1 == -1 ? idx : j1; + I2_1 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][0] : checkIP_ouTer[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][0] : checkIP_ouTer[idx][0] : j2; + I2_2 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][1] : checkIP_ouTer[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][1] : checkIP_ouTer[idx][1] : j2; + + B = model.csswm[p1].lat[I1][J1]; + A1 = model.csswm[p2].lat[I2_1][J2_1], A2 = model.csswm[p2].lat[I2_2][J2_2]; + // V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2]; + // V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2]; + + V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1]; + V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1]; + V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2]; + V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2]; + } + else { + I1 = i1 == -1 ? idx : i1, J1 = j1 == -1 ? idx : j1; + I2_1 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][0] : checkIP_ouTer[idx][0] : i2, J2_1 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][0] : checkIP_ouTer[idx][0] : j2; + I2_2 = i2 == -1 ? reversed ? checkIP_ouTer[NX-1-idx][1] : checkIP_ouTer[idx][1] : i2, J2_2 = j2 == -1 ? reversed ? checkIP_ouTer[NY-1-idx][1] : checkIP_ouTer[idx][1] : j2; + + B = model.csswm[p1].lon[I1][J1]; + A1 = model.csswm[p2].lon[I2_1][J2_1], A2 = model.csswm[p2].lon[I2_2][J2_2]; + // V1 = model.csswm[p2].up[I2_1][J2_1], V2 = model.csswm[p2].up[I2_2][J2_2]; + // V3 = model.csswm[p2].vp[I2_1][J2_1], V4 = model.csswm[p2].vp[I2_2][J2_2]; + + V1 = model.gUpper[I2_1][J2_1][0] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][1] * model.csswm[p2].vp[I2_1][J2_1]; + V3 = model.gUpper[I2_1][J2_1][2] * model.csswm[p2].up[I2_1][J2_1] + model.gUpper[I2_1][J2_1][3] * model.csswm[p2].vp[I2_1][J2_1]; + V2 = model.gUpper[I2_2][J2_2][0] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][1] * model.csswm[p2].vp[I2_2][J2_2]; + V4 = model.gUpper[I2_2][J2_2][2] * model.csswm[p2].up[I2_2][J2_2] + model.gUpper[I2_2][J2_2][3] * model.csswm[p2].vp[I2_2][J2_2]; + + if (A1 > A2 && (p1 == 0 || p2 == 0)) A2 += 2 * M_PI; + if (A1 > B && B < A2) B += 2 * M_PI; + } + + uIP = model.interpolate(A1, A2, V1, V2, B); + vIP = model.interpolate(A1, A2, V3, V4, B); + + if (i1 == -1) { + if (j1 == 1) { + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_D[I1][0] * uIP + model.csswm[p1].IP_ouTer_D[I1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_D[I1][2] * uIP + model.csswm[p1].IP_ouTer_D[I1][3] * vIP; + } + else if (j1 == NY-2) { + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_U[I1][0] * uIP + model.csswm[p1].IP_ouTer_U[I1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_U[I1][2] * uIP + model.csswm[p1].IP_ouTer_U[I1][3] * vIP; + } + } + else if (j1 == -1) { + if (i1 == 1) { + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_L[J1][0] * uIP + model.csswm[p1].IP_ouTer_L[J1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_L[J1][2] * uIP + model.csswm[p1].IP_ouTer_L[J1][3] * vIP; + } + else if (i1 == NX-2) { + model.csswm[p1].up[I1][J1] = model.csswm[p1].IP_ouTer_R[J1][0] * uIP + model.csswm[p1].IP_ouTer_R[J1][1] * vIP; + model.csswm[p1].vp[I1][J1] = model.csswm[p1].IP_ouTer_R[J1][2] * uIP + model.csswm[p1].IP_ouTer_R[J1][3] * vIP; + } + } + } + } + #endif +} diff --git a/src/construction.hpp b/src/construction.hpp index 4183207..5d699f3 100644 --- a/src/construction.hpp +++ b/src/construction.hpp @@ -41,6 +41,8 @@ class CSSWM { double du[NX][NY][2]; double dv[NX][NY][2]; #endif + + double f[NX][NY]; }; CSSWM(); @@ -79,6 +81,7 @@ class CSSWM { double diffusion_ky = KY; double diffusion_ts = TIMETS; double addforcingtime = ADDFORCINGTIME; + double csswm_h_nudge_time = CSSWM_H_NUDGE_TIME; // *********************************************************************************** @@ -195,6 +198,7 @@ class CSSWM { public: static void DiffusionAll(CSSWM &); static void timeFilterAll(CSSWM &); + static void NudgeH(CSSWM &); }; diff --git a/src/define.hpp b/src/define.hpp index 6db39f7..fa2c925 100644 --- a/src/define.hpp +++ b/src/define.hpp @@ -34,6 +34,7 @@ #define ALPHA0 (0) // #define Advection // #define GravityWave +// #define DeformationalFlow // #define SteadyGeostrophy // #define Barotropic // #define Mountain @@ -49,5 +50,6 @@ #define TIMETS (0.06) #define ADDFORCINGTIME (100. * 60.) +#define CSSWM_H_NUDGE_TIME (0.) #endif diff --git a/src/initialField.cpp b/src/initialField.cpp index 0e001b2..29324c6 100644 --- a/src/initialField.cpp +++ b/src/initialField.cpp @@ -34,6 +34,8 @@ void CSSWM::Init::Init2d(CSSWM & model) { #ifdef DeformationalFlow model.csswm[p].hp[i][j] = DeformationalFlowH(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]); + model.csswm[p].up[i][j] = 0.; + model.csswm[p].vp[i][j] = 0.; #endif #ifdef GravityWave @@ -48,6 +50,8 @@ void CSSWM::Init::Init2d(CSSWM & model) { (model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * SteadyGeostrophyV(model.csswm[p].lon[i][j]); model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * SteadyGeostrophyU(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]) + (model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * SteadyGeostrophyV(model.csswm[p].lon[i][j]); + + model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); #endif #ifdef Barotropic @@ -57,6 +61,7 @@ void CSSWM::Init::Init2d(CSSWM & model) { (model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * 0; model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * BarotropicU(model.csswm[p].lat[i][j]) + (model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * 0; + model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); #endif #ifdef Mountain @@ -65,6 +70,7 @@ void CSSWM::Init::Init2d(CSSWM & model) { (model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * MountainV(model.csswm[p].lon[i][j]); model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * MountainU(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]) + (model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * MountainV(model.csswm[p].lon[i][j]); + model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); #endif #ifdef RossbyHaurwitz @@ -73,10 +79,11 @@ void CSSWM::Init::Init2d(CSSWM & model) { (model.gLower[i][j][0] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][1] * model.csswm[p].IA[i][j][3]) * RossbyHaurwitzV(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]); model.csswm[p].vp[i][j] = (model.gLower[i][j][2] * model.csswm[p].IA[i][j][0] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][2]) * RossbyHaurwitzU(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]) + (model.gLower[i][j][2] * model.csswm[p].IA[i][j][1] + model.gLower[i][j][3] * model.csswm[p].IA[i][j][3]) * RossbyHaurwitzV(model.csswm[p].lon[i][j], model.csswm[p].lat[i][j]); + model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); #endif #ifdef EquatorialWave - if (p == 0) { + if (p == 1) { model.csswm[p].h_forcing[i][j] = EquatorialWaveH(model.csswm[p].x[i][j], model.csswm[p].y[i][j]); } else { @@ -86,6 +93,11 @@ void CSSWM::Init::Init2d(CSSWM & model) { model.csswm[p].hp[i][j] = 10454.608791605699 + model.csswm[p].h_forcing[i][j]; model.csswm[p].up[i][j] = 0.; model.csswm[p].vp[i][j] = 0.; + + // double f0 = 0; + // double beta = 2.5 * 10E-11; + // model.csswm[p].f[i][j] = f0 + beta * model.csswm[p].lat[i][j] * 180. / M_PI * (111000.) / 11.5; + model.csswm[p].f[i][j] = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); #endif #if defined(Uniform) @@ -302,7 +314,7 @@ double CSSWM::Init::RossbyHaurwitzV(double lon, double lat) { } double CSSWM::Init::EquatorialWaveH(double x, double y) { - double h0 = 0.08; - double a = 100000., b = 100000.; - return -h0 * exp(-(x*x / (2*(a*a)) + y*y / (2*(b*b)))); + double h0 = 5E-2; + double a = 400. * 1E3, b = 200. * 1E3; + return -h0 * exp(-(x*x / ((a*a)) + y*y / ((b*b)))); } \ No newline at end of file diff --git a/src/iteration.cpp b/src/iteration.cpp index 5027d3a..d146baa 100644 --- a/src/iteration.cpp +++ b/src/iteration.cpp @@ -206,7 +206,6 @@ void CSSWM::Iteration::ph_pt_4(CSSWM &model) { void CSSWM::Iteration::pu_pt_4(CSSWM &model) { double dx_for_u = 0, dy_for_u = 0; double pgH_px = 0, pU2_px = 0, pUV_px = 0, pV2_px = 0, rotationU = 0; - double f; #ifdef Mountain double pgHs_px = 0.; #endif @@ -216,18 +215,6 @@ void CSSWM::Iteration::pu_pt_4(CSSWM &model) { for (int j = 2; j < NY-2; j++) { dx_for_u = 0.5 * (model.csswm[p].x[i+1][j] - model.csswm[p].x[i-1][j]); dy_for_u = 0.5 * (model.csswm[p].y[i][j+1] - model.csswm[p].y[i][j-1]); - - #if defined(SteadyGeostrophy) || defined(Mountain) - f = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); - #elif defined(Barotropic) || defined(RossbyHaurwitz) - f = 2 * OMEGA * sin(model.csswm[p].lat[i][j]); - #elif defined(EquatorialWave) - double f0 = 0; - double beta = 2.5 * 10E-11; - f = f0 + beta * model.csswm[p].lat[i][j] * 180. / M_PI * (111000.); - #else - f = 0; - #endif pgH_px = model.gravity / (12.*dx_for_u) * (-1*model.csswm[p].h[i+2][j] + 8*model.csswm[p].h[i+1][j] - 8*model.csswm[p].h[i-1][j] + 1*model.csswm[p].h[i-2][j]); @@ -251,7 +238,7 @@ void CSSWM::Iteration::pu_pt_4(CSSWM &model) { rotationU = (((-1.*model.csswm[p].v[i+2][j] + 8.*model.csswm[p].v[i+1][j] - 8.*model.csswm[p].v[i-1][j] + 1.*model.csswm[p].v[i-2][j]) / (12.*dx_for_u)) - ((-1.*model.csswm[p].u[i][j+2] + 8.*model.csswm[p].u[i][j+1] - 8.*model.csswm[p].u[i][j-1] + 1.*model.csswm[p].u[i][j-2]) / (12.*dy_for_u)) - + model.sqrtG[i][j] * f) + + model.sqrtG[i][j] * model.csswm[p].f[i][j]) * (model.gUpper[i][j][2] * model.csswm[p].u[i][j] + model.gUpper[i][j][3] * model.csswm[p].v[i][j]); @@ -280,7 +267,6 @@ void CSSWM::Iteration::pu_pt_4(CSSWM &model) { void CSSWM::Iteration::pv_pt_4(CSSWM &model) { double dx_for_v = 0, dy_for_v = 0; double pgH_py = 0, pU2_py = 0, pUV_py = 0, pV2_py = 0, rotationV = 0; - double f; #ifdef Mountain double pgHs_py = 0.; #endif @@ -291,18 +277,6 @@ void CSSWM::Iteration::pv_pt_4(CSSWM &model) { dx_for_v = 0.5 * (model.csswm[p].x[i+1][j] - model.csswm[p].x[i-1][j]); dy_for_v = 0.5 * (model.csswm[p].y[i][j+1] - model.csswm[p].y[i][j-1]); - #if defined(SteadyGeostrophy) || defined(Mountain) - f = 2 * OMEGA * (-cos(model.csswm[p].lon[i][j]) * cos(model.csswm[p].lat[i][j]) * sin(ALPHA0) + sin(model.csswm[p].lat[i][j]) * cos(ALPHA0)); - #elif defined(Barotropic) || defined(RossbyHaurwitz) - f = 2 * OMEGA * sin(model.csswm[p].lat[i][j]); - #elif defined(EquatorialWave) - double f0 = 0; - double beta = 2.5 * 10E-11; - f = f0 + beta * model.csswm[p].lat[i][j] * 180. / M_PI * (111000.); - #else - f = 0; - #endif - pgH_py = model.gravity / (12.*dy_for_v) * (-1.*model.csswm[p].h[i][j+2] + 8.*model.csswm[p].h[i][j+1] - 8.*model.csswm[p].h[i][j-1] + 1.*model.csswm[p].h[i][j-2]); pU2_py = 0.5 /(12.*dy_for_v) * @@ -325,7 +299,7 @@ void CSSWM::Iteration::pv_pt_4(CSSWM &model) { rotationV = (((-1.*model.csswm[p].v[i+2][j] + 8.*model.csswm[p].v[i+1][j] - 8.*model.csswm[p].v[i-1][j] + 1.*model.csswm[p].v[i-2][j]) / (12.*dx_for_v)) - ((-1.*model.csswm[p].u[i][j+2] + 8.*model.csswm[p].u[i][j+1] - 8.*model.csswm[p].u[i][j-1] + 1.*model.csswm[p].u[i][j-2]) / (12.*dy_for_v)) - + model.sqrtG[i][j] * f) * + + model.sqrtG[i][j] * model.csswm[p].f[i][j]) * (model.gUpper[i][j][0] * model.csswm[p].u[i][j] + model.gUpper[i][j][1] * model.csswm[p].v[i][j]); @@ -441,15 +415,15 @@ void CSSWM::Iteration::TimeMarching(CSSWM &model) { CSSWM::NumericalProcess::DiffusionAll(model); #endif - // // Time filter - // #if defined(TIMEFILTER) && !defined(AB2Time) - // CSSWM::NumericalProcess::timeFilterAll(model); - // #endif - - #if defined(TIMEFILTER) + // Time filter + #if defined(TIMEFILTER) && !defined(AB2Time) CSSWM::NumericalProcess::timeFilterAll(model); #endif + if (model.csswm_h_nudge_time != 0) { + CSSWM::NumericalProcess::NudgeH(model); + } + // Add forcing (1 speed) // Coupling time: 600 diff --git a/src/main.cpp b/src/main.cpp index 3f31eec..0df3b59 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -20,6 +20,7 @@ int main(void) { double csswm_diffusion_ky = std::stod(configs["CSSWM_DIFFUSION_KY"]); double csswm_diffusion_ts = std::stod(configs["CSSWM_DIFFUSION_TS"]); double csswm_addforcing_time = std::stod(configs["CSSWM_ADDFORCING_TIME"]); + double csswm_h_nudge_time = std::stod(configs["CSSWM_H_NUDGE_TIME"]); // int csswm_nx = (int) (90 / csswm_dx + 2); // int csswm_ny = (int) (90 / csswm_dy + 2); @@ -40,6 +41,7 @@ int main(void) { model.diffusion_ky = csswm_diffusion_ky; model.diffusion_ts = csswm_diffusion_ts; model.addforcingtime = csswm_addforcing_time; + model.csswm_h_nudge_time = csswm_h_nudge_time; CSSWM::Outputs::create_all_directory(model); CSSWM::Init::Init2d(model); diff --git a/src/numericalProcess.cpp b/src/numericalProcess.cpp index 53d6eb3..3e01b94 100644 --- a/src/numericalProcess.cpp +++ b/src/numericalProcess.cpp @@ -30,4 +30,15 @@ void CSSWM::NumericalProcess::timeFilterAll(CSSWM &model) { } } return; +} + +void CSSWM::NumericalProcess::NudgeH(CSSWM &model) { + for (int p = 0; p < 6; p++) { + for (int i = 0; i < NX; i++) { + for (int j = 0; j < NY; j++) { + model.csswm[p].hp[i][j] -= model.dt / (model.csswm_h_nudge_time) * (model.csswm[p].hp[i][j]-10454.608791605699); + } + } + } + return; } \ No newline at end of file diff --git a/src/transform.cpp b/src/transform.cpp index ae68f7e..6e73418 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -261,18 +261,29 @@ void CSSWM::Cube2Cube_matrix() { get_A(A_tmp, p2, alpha_B_tmp, beta_B_tmp); get_gUpper(gUpper_tmp, alpha_B_tmp, beta_B_tmp); + // matrixMul(gLower_tmp, IA_tmp, tmp1); + // matrixMul(A_tmp, gUpper_tmp, tmp2); + // double tmp1_4[4], tmp2_4[4]; + // count = 0; + // for (int i = 0; i < 2; i++) { + // for (int j = 0; j < 2; j++) { + // tmp1_4[count] = tmp1[i][j]; + // tmp2_4[count] = tmp2[i][j]; + // count++; + // } + // } + // matrixMul(tmp1_4, tmp2_4, mult); + matrixMul(gLower_tmp, IA_tmp, tmp1); - matrixMul(A_tmp, gUpper_tmp, tmp2); - double tmp1_4[4], tmp2_4[4]; + double tmp1_4[4]; count = 0; for (int i = 0; i < 2; i++) { for (int j = 0; j < 2; j++) { tmp1_4[count] = tmp1[i][j]; - tmp2_4[count] = tmp2[i][j]; count++; } } - matrixMul(tmp1_4, tmp2_4, mult); + matrixMul(tmp1_4, A_tmp, mult); if (i1 == -1) { if (j1 == 0) { @@ -366,18 +377,29 @@ void CSSWM::Cube2Cube_matrix() { get_A(A_tmp, p2, alpha_B_tmp, beta_B_tmp); get_gUpper(gUpper_tmp, alpha_B_tmp, beta_B_tmp); + // matrixMul(gLower_tmp, IA_tmp, tmp1); + // matrixMul(A_tmp, gUpper_tmp, tmp2); + // double tmp1_4[4], tmp2_4[4]; + // count = 0; + // for (int i = 0; i < 2; i++) { + // for (int j = 0; j < 2; j++) { + // tmp1_4[count] = tmp1[i][j]; + // tmp2_4[count] = tmp2[i][j]; + // count++; + // } + // } + // matrixMul(tmp1_4, tmp2_4, mult); + matrixMul(gLower_tmp, IA_tmp, tmp1); - matrixMul(A_tmp, gUpper_tmp, tmp2); - double tmp1_4[4], tmp2_4[4]; + double tmp1_4[4]; count = 0; for (int i = 0; i < 2; i++) { for (int j = 0; j < 2; j++) { tmp1_4[count] = tmp1[i][j]; - tmp2_4[count] = tmp2[i][j]; count++; } } - matrixMul(tmp1_4, tmp2_4, mult); + matrixMul(tmp1_4, A_tmp, mult); if (i1 == -1) { if (j1 == 1) {