diff --git a/SRC/FEMAIN/FastEddy.c b/SRC/FEMAIN/FastEddy.c index 56a1204..e48be5d 100644 --- a/SRC/FEMAIN/FastEddy.c +++ b/SRC/FEMAIN/FastEddy.c @@ -352,9 +352,9 @@ int main(int argc, char **argv){ printf("FastEddy MAin timestepping loop: Reading new BdyPlanes at it=%d...\n",it); fflush(stdout); errorCode = timeIntBdyPlaneUpdates(); - //if((cellpertSelector==1)&&(cellpert_tvcp==1)){ // DME update CP parameters with dynamic LBCs - // errorCode = hydroCoreTVCP(dt,mpi_rank_world); - //} // end if((cellpertSelector==1)&&(cellpert_tvcp==1)) + if((cellpertSelector==1)&&(cellpert_tvcp==1)){ // update CP parameters with dynamic LBCs + errorCode = hydro_coreTVCP(dt); + } // end if((cellpertSelector==1)&&(cellpert_tvcp==1)) }//end if hydroBCs == 1 } MPI_Barrier(MPI_COMM_WORLD); diff --git a/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu index 282e1fb..67ddced 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu @@ -99,6 +99,19 @@ extern "C" int cuda_hydroCoreDeviceBuildCPmethod(int simTime_it){ return(errorCode); }//end cuda_hydroCoreDeviceBuildCPmethod() +/*----->>>>> extern "C" int cuda_hydroCoreTVCP(); ----------------------------------------------------------- +* Updates device-sided parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions +*/ +extern "C" int cuda_hydroCoreTVCP(){ + int errorCode = CUDA_CELLPERT_SUCCESS; + + cudaMemcpyToSymbol(cellpert_amp_d, &cellpert_amp, sizeof(float)); + cudaMemcpyToSymbol(cellpert_ktop_d, &cellpert_ktop, sizeof(int)); + cudaMemcpyToSymbol(cellpert_nts_d, &cellpert_nts, sizeof(int)); + + return(errorCode); +} //end cuda_hydroCoreTVCP() + __global__ void cudaDevice_hydroCoreUnitTestCompleteCellPerturbation(float* hydroFlds, float* randcp_d, int my_mpi, int numpx, int numpy){ int i,j,k,ijk; diff --git a/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice_cu.h index 7d28c92..7216e5d 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice_cu.h +++ b/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice_cu.h @@ -31,9 +31,14 @@ extern "C" int cuda_cellpertDeviceSetup(); extern "C" int cuda_cellpertDeviceCleanup(); /*----->>>>> extern "C" int cuda_hydroCoreDeviceBuildCPmethod(); -------------------------------------------------- - * * This routine provides the externally callable cuda-kernel call to perform a call to cell perturbation method */ +* This routine provides the externally callable cuda-kernel call to perform a call to cell perturbation method */ extern "C" int cuda_hydroCoreDeviceBuildCPmethod(int simTime_it); +/*----->>>>> extern "C" int cuda_hydroCoreTVCP(); ----------------------------------------------------------- +* Updates device-sided parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions +*/ +extern "C" int cuda_hydroCoreTVCP(); + __global__ void cudaDevice_hydroCoreUnitTestCompleteCellPerturbation(float* hydroFlds, float* randcp_d, int my_mpi, int numpx, int numpy); /*----->>>>> __device__ void cudaDevice_CellPerturbation(); -------------------------------------------------- diff --git a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu index e6fd987..7f38505 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu @@ -121,7 +121,7 @@ extern "C" int cuda_hydroCoreDeviceSetup(){ cudaMemcpyToSymbol(corioConstHorz_d, &corioConstHorz, sizeof(float)); cudaMemcpyToSymbol(corioConstVert_d, &corioConstVert, sizeof(float)); cudaMemcpyToSymbol(corioLS_fact_d, &corioLS_fact, sizeof(float)); - cudaMemcpyToSymbol(kappa_d, &kappa, sizeof(float)); // DME + cudaMemcpyToSymbol(kappa_d, &kappa, sizeof(float)); cudaMemcpyToSymbol(L_v_d, &L_v, sizeof(float)); gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ @@ -195,7 +195,11 @@ extern "C" int cuda_hydroCoreDeviceSetup(){ if (surflayerSelector > 0) { errorCode = cuda_surfaceLayerDeviceSetup(); } - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMalloc calls*/ + + /* CELL PERTURBATION METHOD */ + if (cellpertSelector > 0) { + errorCode = cuda_cellpertDeviceSetup(); + } /* CANOPY */ if (canopySelector > 0){ @@ -271,6 +275,9 @@ extern "C" int cuda_hydroCoreDeviceCleanup(){ if (surflayerSelector > 0) { errorCode = cuda_surfaceLayerDeviceCleanup(); } + if (cellpertSelector > 0) { + errorCode = cuda_cellpertDeviceCleanup(); + } if (canopySelector > 0) { errorCode = cuda_canopyDeviceCleanup(); } diff --git a/SRC/HYDRO_CORE/CUDA/cuda_surfaceLayerDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_surfaceLayerDevice.cu index 6d69d55..eb2e68a 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_surfaceLayerDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_surfaceLayerDevice.cu @@ -598,7 +598,7 @@ __device__ void cudaDevice_offshoreRoughness(float* z0m, float* z0t, float* fric float alpha_charnock = 0.018; float alpha_charnock_mod; float wspd_1; - float air_vis = 1.5e-5; // kinematic air viscosity (DME: make it T dependent ...) + float air_vis = 1.5e-5; // kinematic air viscosity (make it T dependent ...) float z0_m2t_fact = 0.1; // ratio of z0t/z0m int z0_m2t_opt = 1; // 0; // ==0 (constant), ==1 (roughness Re dependent) float Ren; diff --git a/SRC/HYDRO_CORE/hydro_core.c b/SRC/HYDRO_CORE/hydro_core.c index 7351aa4..f43ee70 100644 --- a/SRC/HYDRO_CORE/hydro_core.c +++ b/SRC/HYDRO_CORE/hydro_core.c @@ -225,6 +225,9 @@ int cellpert_ndbc; /* number of cells normal to domain lateral boundaries int cellpert_kbottom; /* z-grid point where the perturbations start */ int cellpert_ktop; /* z-grid point where the perturbations end */ int cellpert_ktop_prev[4];/* z-grid point where the perturbations end array previous time step */ +int cellpert_tvcp; /* time-varying CP method selector: 0= off, 1= on (when hydroBCs == 1) */ +float cellpert_eckert; /* Eckert number for the potential temperature perturbations (when hydroBCs == 1) */ +float cellpert_tsfact; /* factor on the refreshing perturbation time scale (when hydroBCs == 1) */ /*--- Rayleigh Damping Layer ---*/ int dampingLayerSelector; // Rayleigh Damping Layer selector @@ -420,10 +423,12 @@ int hydro_coreGetParams(){ errorCode = queryIntegerParameter("cellpert_kbottom", &cellpert_kbottom, 1, 10, PARAM_MANDATORY); cellpert_ktop = 20; // Default to 20th grid point above surface errorCode = queryIntegerParameter("cellpert_ktop", &cellpert_ktop, 0, 200, PARAM_MANDATORY); - cellpert_ktop_prev[0] = cellpert_ktop; // use params in as previous cellpert_ktop value - cellpert_ktop_prev[1] = cellpert_ktop; // use params in as previous cellpert_ktop value - cellpert_ktop_prev[2] = cellpert_ktop; // use params in as previous cellpert_ktop value - cellpert_ktop_prev[3] = cellpert_ktop; // use params in as previous cellpert_ktop value + cellpert_tvcp = 0; // Default to off + errorCode = queryIntegerParameter("cellpert_tvcp", &cellpert_tvcp, 0, 1, PARAM_MANDATORY); + cellpert_eckert = 0.2; // Default to Ec = 0.2 + errorCode = queryFloatParameter("cellpert_eckert", &cellpert_eckert, 0.0, 10.0, PARAM_MANDATORY); + cellpert_tsfact = 1.0; // Default to cellpert_tsfact = 1.0 + errorCode = queryFloatParameter("cellpert_tsfact", &cellpert_tsfact, 0.0, 10.0, PARAM_MANDATORY); } // lsfSelector = 0; // Default to off @@ -757,6 +762,9 @@ int hydro_coreInit(){ printParameter("cellpert_ndbc", "number of cells normal to domain lateral boundaries"); printParameter("cellpert_kbottom", "z-grid point where the perturbations start"); printParameter("cellpert_ktop", "z-grid point where the perturbations end"); + printParameter("cellpert_tvcp", "time-varying CP method selector: 0= off, 1= on"); + printParameter("cellpert_eckert", "Eckert number for the potential temperature perturbations (when hydroBCs == 1)"); + printParameter("cellpert_tsfact", "factor on the refreshing perturbation time scale (when hydroBCs == 1)"); } printComment("----------: RAYLEIGH DAMPING LAYER ---"); printParameter("dampingLayerSelector", "Rayleigh damping layer selector: 0= off, 1= on."); @@ -938,6 +946,14 @@ int hydro_coreInit(){ MPI_Bcast(&cellpert_ndbc, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&cellpert_kbottom, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&cellpert_ktop, 1, MPI_INT, 0, MPI_COMM_WORLD); + //Initialize the ktop_prev vector on all ranks + cellpert_ktop_prev[0] = cellpert_ktop; // use params value as previous cellpert_ktop value + cellpert_ktop_prev[1] = cellpert_ktop; // use params value as previous cellpert_ktop value + cellpert_ktop_prev[2] = cellpert_ktop; // use params value as previous cellpert_ktop value + cellpert_ktop_prev[3] = cellpert_ktop; // use params value as previous cellpert_ktop value + MPI_Bcast(&cellpert_tvcp, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&cellpert_eckert, 1, MPI_FLOAT, 0, MPI_COMM_WORLD); + MPI_Bcast(&cellpert_tsfact, 1, MPI_FLOAT, 0, MPI_COMM_WORLD); } MPI_Bcast(&dampingLayerSelector, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&dampingLayerDepth, 1, MPI_FLOAT, 0, MPI_COMM_WORLD); @@ -1275,6 +1291,24 @@ int hydro_coreInit(){ fflush(stdout); } + if(cellpertSelector>0){ //Cell Perturbation parameters (time-varying when cellpert_tvcp == 1 and hydroBCs==1) + errorCode = sprintf(&fldName[0],"cellpert_amp"); + errorCode = ioRegisterVar(&fldName[0], "float", 1, dims1dTD, &cellpert_amp); + printf("cellpert:Variable = %s stored at %p, has been registered with IO.\n", + &fldName[0],&cellpert_amp); + fflush(stdout); + errorCode = sprintf(&fldName[0],"cellpert_nts"); + errorCode = ioRegisterVar(&fldName[0], "int", 1, dims1dTD, &cellpert_nts); + printf("cellpert:Variable = %s stored at %p, has been registered with IO.\n", + &fldName[0],&cellpert_nts); + fflush(stdout); + errorCode = sprintf(&fldName[0],"cellpert_ktop"); + errorCode = ioRegisterVar(&fldName[0], "int", 1, dims1dTD, &cellpert_ktop); + printf("cellpert:Variable = %s stored at %p, has been registered with IO.\n", + &fldName[0],&cellpert_ktop); + fflush(stdout); + } // end of cellpertSelector > 0 + if(canopySelector>0){ canopy_lad = memAllocateFloat3DField(Nxp, Nyp, Nzp, Nh, "canopy_lad"); errorCode = sprintf(&fldName[0],"CanopyLAD"); @@ -1363,7 +1397,7 @@ int hydro_coreInit(){ }else{ nBndyVars = Nhydro; nSurfBndyVars = 1; //Only allows tskin - } //end if moisture is on else not //JAS NOTE: Doesn't handle any AuxScalars or TKE-related Prog. variables. + } //end if moisture is on else not //NOTE: Doesn't handle any AuxScalars or TKE-related Prog. variables. XZBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Nx*Nz*sizeof(float) ); YZBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Ny*Nz*sizeof(float) ); XYBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Nx*Ny*sizeof(float) ); @@ -2321,6 +2355,209 @@ int hydro_coreScatterFieldBndyPlanes(int Nfields){ return(errorCode); }//end hydro_coreScatterFieldBndyPlanes() +/*----->>>>> int hydro_coreTVCP(); ----------------------------------------------------------- +* Updates model parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions. +*/ +int hydro_coreTVCP(float dt){ + int errorCode = HYDRO_CORE_SUCCESS; + + int bdy_id; // 0==North; 1==East; 2==South; 3==West; + float cellpert_amp_array[4]; + int cellpert_ktop_array[4]; + int cellpert_nts_array[4]; + float cellpert_amp_aver; + int cellpert_ktop_aver; + int cellpert_nts_aver; + + /*--- North boundary ---*/ + bdy_id = 0; + hydro_coreTVCP_LBCparams(bdy_id, XZBdyPlanesGlobal, cellpert_amp_array, cellpert_ktop_array, cellpert_ktop_prev, cellpert_nts_array, dt); + /*--- East boundary ---*/ + bdy_id = 1; + hydro_coreTVCP_LBCparams(bdy_id, YZBdyPlanesGlobal, cellpert_amp_array, cellpert_ktop_array, cellpert_ktop_prev, cellpert_nts_array, dt); + /*--- South boundary ---*/ + bdy_id = 2; + hydro_coreTVCP_LBCparams(bdy_id, XZBdyPlanesGlobal, cellpert_amp_array, cellpert_ktop_array, cellpert_ktop_prev, cellpert_nts_array, dt); + /*--- West boundary ---*/ + bdy_id = 3; + hydro_coreTVCP_LBCparams(bdy_id, YZBdyPlanesGlobal, cellpert_amp_array, cellpert_ktop_array, cellpert_ktop_prev, cellpert_nts_array, dt); + + cellpert_amp_aver = 0.25*(cellpert_amp_array[0]+cellpert_amp_array[1]+cellpert_amp_array[2]+cellpert_amp_array[3]); + cellpert_ktop_aver = 0.25*(cellpert_ktop_array[0]+cellpert_ktop_array[1]+cellpert_ktop_array[2]+cellpert_ktop_array[3]); + cellpert_nts_aver = 0.25*(cellpert_nts_array[0]+cellpert_nts_array[1]+cellpert_nts_array[2]+cellpert_nts_array[3]); + cellpert_amp = cellpert_amp_aver; + cellpert_ktop = cellpert_ktop_aver; + cellpert_nts = cellpert_nts_aver; + + cellpert_ktop_prev[0] = cellpert_ktop; + cellpert_ktop_prev[1] = cellpert_ktop; + cellpert_ktop_prev[2] = cellpert_ktop; + cellpert_ktop_prev[3] = cellpert_ktop; + +//#define TVCP_DEBUG +#ifdef TVCP_DEBUG + printf("cellpert_amp_array[%d] = %f \n",0,cellpert_amp_array[0]); + printf("cellpert_amp_array[%d] = %f \n",1,cellpert_amp_array[1]); + printf("cellpert_amp_array[%d] = %f \n",2,cellpert_amp_array[2]); + printf("cellpert_amp_array[%d] = %f \n",3,cellpert_amp_array[3]); + printf("cellpert_ktop_array[%d] = %d \n",0,cellpert_ktop_array[0]); + printf("cellpert_ktop_array[%d] = %d \n",1,cellpert_ktop_array[1]); + printf("cellpert_ktop_array[%d] = %d \n",2,cellpert_ktop_array[2]); + printf("cellpert_ktop_array[%d] = %d \n",3,cellpert_ktop_array[3]); + printf("cellpert_nts_array[%d] = %d \n",0,cellpert_nts_array[0]); + printf("cellpert_nts_array[%d] = %d \n",1,cellpert_nts_array[1]); + printf("cellpert_nts_array[%d] = %d \n",2,cellpert_nts_array[2]); + printf("cellpert_nts_array[%d] = %d \n",3,cellpert_nts_array[3]); +#endif +if(mpi_rank_world == 0){ + printf("---------------------- \n"); + printf("Updating CP parameters \n"); + printf("---------------------- \n"); + printf("cellpert_amp = %f K \n",cellpert_amp); + printf("cellpert_ktop = %d grid points \n",cellpert_ktop); + printf("cellpert_nts = %d time steps \n",cellpert_nts); + fflush(stdout); +} + + return(errorCode); +}//end hydro_coreTVCP() + +/*----->>>>> int hydro_coreTVCP_LBCparams(); ----------------------------------------------------------- +* Computes model parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions. +*/ +int hydro_coreTVCP_LBCparams(int bdy_id, float* var_LBCplane, float* cellpert_amp_array, + int* cellpert_ktop_array, int* cellpert_ktop_prev, int* cellpert_nts_array, float dt){ + int errorCode = HYDRO_CORE_SUCCESS; + int bndyfldStride_LZ; + int Nl, k, l, ind_lk; + int bdy_lh; + int ku_s, kv_s, krho_s, kth_s; + float u_tmp, v_tmp, rho_tmp, th_tmp, ws_tmp; + float uz_tmp[Nz]; + float thz_tmp[Nz]; + float ws_tmp_2, ws_diff, ws_abl; + float th_min, th_diff, th_diff_2; + float ntsf_tmp; + int ku_max; + int kth_zi; + int kzi_tmp, kzi_tmp_old; + int k_low = 3; + float dzidt_tend = 0.05; + int dzidt_max; + float dzi_tend_sign; + + if(bdy_id==0 || bdy_id==2){ // N or S boundary planes + Nl = Nx; + }else{ // E or W boundary planes + Nl = Ny; + } + bndyfldStride_LZ = Nl*Nz; + + bdy_lh = 0; //Initialize for safety + if(bdy_id == 0){ // N + bdy_lh = 1; + } else if(bdy_id == 1){ // E + bdy_lh = 1; + } else if(bdy_id == 2){ // S + bdy_lh = 0; + } else if(bdy_id == 3){ // W + bdy_lh = 0; + } + + + //printf("bdy_id = %d Nl = %d \n",bdy_id, Nl); + // //fflush(stdout); + + ku_s = (U_INDX*2+bdy_lh)*bndyfldStride_LZ; + kv_s = (V_INDX*2+bdy_lh)*bndyfldStride_LZ; + krho_s = (RHO_INDX*2+bdy_lh)*bndyfldStride_LZ; + kth_s = (THETA_INDX*2+bdy_lh)*bndyfldStride_LZ; + + for(k=0; k < Nz; k++){ + uz_tmp[k] = 0.0; + thz_tmp[k] = 0.0; + for(l=0; l < Nl; l++){ + ind_lk = l*Nz+k; + u_tmp = var_LBCplane[ku_s+ind_lk]; + v_tmp = var_LBCplane[kv_s+ind_lk]; + rho_tmp = var_LBCplane[krho_s+ind_lk]; + th_tmp = var_LBCplane[kth_s+ind_lk]; + + u_tmp = u_tmp/rho_tmp; + v_tmp = v_tmp/rho_tmp; + th_tmp = th_tmp/rho_tmp; + ws_tmp = sqrt(pow(u_tmp,2.0)+pow(v_tmp,2.0)); + + uz_tmp[k] = uz_tmp[k] + ws_tmp; + thz_tmp[k] = thz_tmp[k] + th_tmp; + } + uz_tmp[k] = uz_tmp[k]/(float)(Nl); + thz_tmp[k] = thz_tmp[k]/(float)(Nl); + } + + ws_tmp = 0.0; + th_min = 500.0; + for(k=0; k < Nz; k++){ + // wind speed maximum // + ws_tmp_2 = uz_tmp[k]; + ws_diff = ws_tmp_2 - ws_tmp; + if(ws_diff > 0.0){ + ku_max = k; + } + ws_tmp = ws_tmp_2; + // min theta // + th_tmp = thz_tmp[k]; + th_diff = th_tmp - th_min; + if (th_diff < 0.0){ + th_min = th_tmp; + } + } + + // closest level to min theta + 1.5 K // + th_diff_2 = 500.0; + for(k=0; k < Nz; k++){ + th_diff = fabs(thz_tmp[k] - (th_min + 1.5)); + if(th_diff < th_diff_2){ + th_diff_2 = th_diff; + kth_zi = k; + } + } + kzi_tmp = floor(0.5*(ku_max+kth_zi)); + + // limit ABL rate of change // + kzi_tmp_old = cellpert_ktop_prev[bdy_id]; + dzidt_max = dzidt_tend*kzi_tmp_old; + dzi_tend_sign = copysignf(1.0,kzi_tmp - kzi_tmp_old); + + if(dzi_tend_sign >= 0.0){ + if(kzi_tmp_old+dzidt_max < kzi_tmp){ + kzi_tmp = kzi_tmp_old+dzidt_max; + } + } else{ + if(kzi_tmp_old+dzidt_max > kzi_tmp){ + kzi_tmp = kzi_tmp_old-dzidt_max; + } + } + + // ABL wind speed // + ws_abl = 0.0; + for(k=k_low; k < kzi_tmp; k++){ + ws_abl = ws_abl + uz_tmp[k]; + } + ws_abl = ws_abl/(float)(kzi_tmp-k_low); + + //printf("k_low = %d kzi_tmp = %d, ws_abl = %f\n",k_low,kzi_tmp,ws_abl); + // //fflush(stdout); + + cellpert_ktop_array[bdy_id] = kzi_tmp; + cellpert_amp_array[bdy_id] = pow(uz_tmp[kzi_tmp],2.0)/(cellpert_eckert*cp_gas); + ntsf_tmp = (float)(cellpert_gppc*cellpert_ndbc*cellpert_tsfact)*d_xi/ws_abl; + cellpert_nts_array[bdy_id] = floor(ntsf_tmp/dt); + + return(errorCode); +}//end hydro_coreTVCP_LBCparams() + + /*----->>>>> int hydro_coreStateLogDump(); ------------------------------------------------------- * Utility function to produce field state summaries for a desired set of hydro_core fields. */ diff --git a/SRC/HYDRO_CORE/hydro_core.h b/SRC/HYDRO_CORE/hydro_core.h index de6f880..c9446a6 100644 --- a/SRC/HYDRO_CORE/hydro_core.h +++ b/SRC/HYDRO_CORE/hydro_core.h @@ -33,7 +33,7 @@ /*#################------------------- HYDRO_CORE module variable declarations ---------------------#################*/ /* Parameters */ extern int Nhydro; /*Number of prognostic variable fields under hydro_core */ -extern int hydroBCs; /*selector for hydro BC set. 1= Dirichlet lateral, ceiling and surface boundary conditions (LAD), +extern int hydroBCs; /*selector for hydro BC set. 1= Dirichlet lateral, ceiling and surface boundary conditions (Limited Area Domain -- LAD), 2= periodicHorizVerticalAbl */ extern int hydroForcingWrite; /*switch for dumping forcing fields of prognostic variables. 0-off (default), 1= on*/ @@ -91,8 +91,10 @@ extern float *SURFBdyPlanesNext; /*Base Adress of memory block for surfaceVar /*---PRESSURE_GRADIENT_FORCE*/ extern int pgfSelector; /*Pressure Gradient Force (pgf) selector: 0=off, 1=on*/ + /*---BUOYANCY*/ extern int buoyancySelector; /*buoyancy Force selector: 0=off, 1=on*/ + /*---CORIOLIS*/ extern int coriolisSelector; /*coriolis Force selector: 0= none, 1= Horiz.-only, 2=Horz. & Vert.*/ extern float coriolisLatitude; /*Charactersitc latitude in degrees from equator of the LES domain*/ @@ -100,6 +102,7 @@ extern float corioConstHorz; /*Latitude dependent horizontal Coriolis term c extern float corioConstVert; /*Latitude dependent Vertical Coriolis term constant */ extern int coriolis_LAD; /*Coriolis force selector for LAD BC cases (hydroBCs==1): 0=off, 1=on*/ extern float corioLS_fact; /*large-scale factor on Coriolis term*/ + /*---TURBULENCE*/ extern int turbulenceSelector; /*turbulence scheme selector: 0= none, 1= Lilly/Smagorinsky */ extern int TKESelector; /* Prognostic TKE selector: 0= none, 1= Prognostic */ @@ -109,6 +112,7 @@ extern float c_s; /* Smagorinsky turbulence model constant used for turbulen extern float c_k; /* Lilly turbulence model constant used for turbulenceSelector = 1 with TKESelector > 0 */ extern float *sgstkeScalars; /* Base Adress of memory containing all prognostic "sgstke" variable fields */ extern float *sgstkeScalarsFrhs; /* Base Adress of memory containing all prognostic "sgstke" RHS forcing fields */ + /*---DIFFUSION*/ extern int diffusionSelector; /*diffusion Term-type selector: 0= none, 1= constant, 2= scalar turbulent-diffusivity*/ extern float nu_0; /* constant diffusivity used when diffusionSelector = 1 */ @@ -116,10 +120,12 @@ extern float *hydroTauFlds; /*Base address for scratch/work Tau tensor arra extern float *hydroDiffTauXFlds; /*Base address for diffusion TauX arrays for all prognostic fields*/ extern float *hydroDiffTauYFlds; /*Base address for diffusion TauY arrays for all prognostic fields*/ extern float *hydroDiffTauZFlds; /*Base address for diffusion TauZ arrays for all prognostic fields*/ + /*---ADVECTION*/ extern int advectionSelector; /*advection scheme selector: 0= 1st-order upwind, 2= 3rd-order QUICK */ extern float b_hyb; /*hybrid advection scheme parameter: 0.0= lower-order upwind, 1.0=higher-order cetered, 0.0 < b_hyb < 1.0 = hybrid */ + /*---SURFACE LAYER*/ extern int surflayerSelector; /*Monin-Obukhov surface layer selector: 0= off, 1= on */ extern float surflayer_z0; /* roughness length (momentum) */ @@ -146,6 +152,7 @@ extern float surflayer_ideal_amp; /*maximum amplitude of the idealized sinusoida extern float surflayer_ideal_qts; /*start time (seconds) for idealized sinusoidal surf. forcing of latent heat flux*/ extern float surflayer_ideal_qte; /*end time (seconds) for idealized sinusoidal surf. forcing of latent heat flux*/ extern float surflayer_ideal_qamp; /*maximum amplitude of idealized sinusoidal surface forcing of latent heat flux*/ + /*OFFSHORE ROUGNESS PARAMETERS*/ extern int surflayer_offshore; /* offshore selector: 0=off, 1=on */ extern int surflayer_offshore_opt; /* offshore roughness parameterization: ==0 (Charnock), ==1 (Charnock with variable alpha), ==2 (Taylor & Yelland), ==3 (Donelan), ==4 (Drennan), ==5 (Porchetta) */ @@ -156,12 +163,14 @@ extern float surflayer_offshore_cp; /* wave phase speed */ extern float surflayer_offshore_theta; /* wave/wind angle */ extern int surflayer_offshore_visc; /* viscous term on z0m: 0=off, 1=on (default) */ extern float* sea_mask; /* Base Address of memory containing sea mask 0,1 field */ + /*---CANOPY*/ extern int canopySelector; /* canopy selector: 0=off, 1=on */ extern int canopySkinOpt; /* canopy selector to use additional skin friction effect on drag coefficient: 0=off, 1=on */ extern float canopy_cd; /* non-dimensional canopy drag coefficient */ extern float canopy_lf; /* representative canopy element length scale */ extern float *canopy_lad; /* Base Address of memory containing leaf area density (LAD) field [m^{-1}] */ + /*---LARGE SCALE FORCING*/ extern int lsfSelector; /* large-scale forcings selector: 0=off, 1=on */ extern float lsf_w_surf; /* lsf to w at the surface */ @@ -182,6 +191,7 @@ extern float lsf_qv_zlev2; /* lsf to qv height 2 */ extern int lsf_horMnSubTerms; /* Switch 0=off, 1=on */ extern int lsf_numPhiVars; /* number of variables in the slabMeanPhiProfiles set (e.g. rho,u,v,theta,qv=5) */ extern float lsf_freq; /* large-scale forcing frequency (seconds) */ + /*---MOISTURE*/ extern int moistureSelector; /* moisture selector: 0=off, 1=on */ extern int moistureNvars; /* number of moisture species */ @@ -197,6 +207,7 @@ extern int moistureAdvSelectorQi; /* moisture advection scheme selector for non- extern float moistureCondTscale; /* relaxation time in seconds */ extern int moistureCondBasePres; /* selector to use base pressure for microphysics */ extern float moistureMPcallTscale; /* time scale for microphysics to be called */ + /*---FILTERS*/ extern int filterSelector; /* explicit filter selector: 0=off, 1=on */ extern int filter_6thdiff_vert; /* vertical 6th-order filter on w selector: 0=off, 1=on */ @@ -204,6 +215,7 @@ extern float filter_6thdiff_vert_coeff; /* vertical 6th-order filter w factor: extern int filter_6thdiff_hori; /* horizontal 6th-order filter on rho,theta,qv selector: 0=off, 1=on */ extern float filter_6thdiff_hori_coeff; /* horizontal 6th-order filter factor: 0.0=off, 1.0=full */ extern int filter_divdamp; /* divergence damping selector: 0=off, 1=on */ + /*---CELL PERTURBATION METHOD*/ extern int cellpertSelector; /* CP method selector: 0= off, 1= on */ extern int cellpert_sw2b; /* switch to do: 0= all four lateral boundaries, 1= only south & west boundaries, 2= only south boundary */ @@ -214,6 +226,10 @@ extern int cellpert_ndbc; /* number of cells normal to domain lateral bou extern int cellpert_kbottom; /* z-grid point where the perturbations start */ extern int cellpert_ktop; /* z-grid point where the perturbations end */ extern int cellpert_ktop_prev[4];/* z-grid point where the perturbations end array previous time step */ +extern int cellpert_tvcp; /* time-varying CP method selector: 0= off, 1= on (when hydroBCs == 1) */ +extern float cellpert_eckert; /* Eckert number for the potential temperature perturbations (when hydroBCs == 1) */ +extern float cellpert_tsfact; /* factor on the refreshing perturbation time scale (when hydroBCs == 1) */ + /*---RAYLEIGH DAMPING LAYER*/ extern int dampingLayerSelector; // Rayleigh Damping Layer selector extern float dampingLayerDepth; // Rayleigh Damping Layer Depth @@ -308,6 +324,17 @@ int hydro_coreScatterFieldBndyPlanes(int Nfields); */ int hydro_coreReadNextBndyPlanesFile(); +/*----->>>>> int hydroi_coreTVCP(); ----------------------------------------------------------- + * Updates model parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions. + */ +int hydro_coreTVCP(float dt); + +/*----->>>>> int hydro_coreTVCP_LBCparams(); ----------------------------------------------------------- + * Computes model parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions. + */ +int hydro_coreTVCP_LBCparams(int bdy_id, float* var_LBCplane, float* cellpert_amp_array, + int* cellpert_ktop_array, int* cellpert_ktop_prev, int* cellpert_nts_array, float dt); + /*----->>>>> int hydro_coreFldStateLogDump(float * Fld); -------------------------------------------------------- * Utility function to carry out a log-dump summary of the hydro_core state */ diff --git a/SRC/MEM_UTILS/mem_utils.c b/SRC/MEM_UTILS/mem_utils.c index 9f4cbee..2276b99 100644 --- a/SRC/MEM_UTILS/mem_utils.c +++ b/SRC/MEM_UTILS/mem_utils.c @@ -51,7 +51,7 @@ int mem_utilsInit(){ * Used to allocate memory for a 2-Dmensional field through the MEM_UTILS module. This field-memory space will * be aligned on ALIGN_SIZE bytes using posix_memalign. * */ -float * memAllocateFloat2DField(int iN, int jN, int halo_extent, char *fieldName){ // DME +float * memAllocateFloat2DField(int iN, int jN, int halo_extent, char *fieldName){ float *field; void *m_field; void *memsetReturnVal; @@ -78,7 +78,7 @@ float * memAllocateFloat2DField(int iN, int jN, int halo_extent, char *fieldName * Used to allocate memory for a 2-Dmensional field made of n vector of one spatial dimention through the MEM_UTILS module. This field-memory space will * be aligned on ALIGN_SIZE bytes using posix_memalign. * */ -float * memAllocateFloat2DFieldN1D(int nN, int kN, int halo_extent, char *fieldName){ // DME +float * memAllocateFloat2DFieldN1D(int nN, int kN, int halo_extent, char *fieldName){ float *field; void *m_field; void *memsetReturnVal; diff --git a/SRC/MEM_UTILS/mem_utils.h b/SRC/MEM_UTILS/mem_utils.h index df4a0b1..b3637a5 100644 --- a/SRC/MEM_UTILS/mem_utils.h +++ b/SRC/MEM_UTILS/mem_utils.h @@ -40,10 +40,10 @@ int mem_utilsInit(); * Used to allocate memory for a 2-Dmensional field through the MEM_UTILS module. This field-memory space will * be aligned on ALIGN_SIZE bytes using posix_memalign. * */ -float * memAllocateFloat2DField(int iN, int jN, int halo_extent, char *fieldName); // DME +float * memAllocateFloat2DField(int iN, int jN, int halo_extent, char *fieldName); /*----->>>>> float * memAllocateFloat2DFieldN1D(nN, kN, halo_extent, fieldName); ------------------------------- - * Used to allocate memory for a 2-Dmensional field made of n vector of one spatial dimention through the MEM_UTILS module. This field-memory space will + * Used to allocate memory for a 2-Dmensional field made of n vector of one spatial dimension through the MEM_UTILS module. This field-memory space will * be aligned on ALIGN_SIZE bytes using posix_memalign. * */ float * memAllocateFloat2DFieldN1D(int nN, int kN, int halo_extent, char *fieldName); diff --git a/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu b/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu index c8a1e7d..0ae6f08 100644 --- a/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu +++ b/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu @@ -120,9 +120,9 @@ extern "C" int cuda_timeIntDeviceCommence(int it){ printf("cuda_timeIntDeviceCommence:Updating device-side BdyPlanes at it=%d...\n",it); fflush(stdout); errorCode = cuda_hydroCoreDeviceBdyPlanesUpdate(); - // if((cellpertSelector==1)&&(cellpert_tvcp==1)){ // DME update CP parameters with dynamic LBCs - // errorCode = cuda_hydroCoreTVCP(dt,mpi_rank_world); - // } // end if((cellpertSelector==1)&&(cellpert_tvcp==1)) + if((cellpertSelector==1)&&(cellpert_tvcp==1)){ // update CP parameters with dynamic LBCs + errorCode = cuda_hydroCoreTVCP(); + } // end if((cellpertSelector==1)&&(cellpert_tvcp==1)) } }//end if hydroBCs==1