1 static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \n"; 2 /* 3 Contributed by Steve Froehlich, Illinois Institute of Technology 4 5 Usage: 6 mpiexec -n <np> ./ex5 [options] 7 ./ex5 -help [view petsc options] 8 ./ex5 -ts_type sundials -ts_view 9 ./ex5 -da_grid_x 20 -da_grid_y 20 -log_view 10 ./ex5 -da_grid_x 20 -da_grid_y 20 -ts_type rosw -ts_atol 1.e-6 -ts_rtol 1.e-6 11 ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4 12 */ 13 14 /* 15 ----------------------------------------------------------------------- 16 17 Governing equations: 18 19 R = s*(Ea*Ta^4 - Es*Ts^4) 20 SH = p*Cp*Ch*wind*(Ta - Ts) 21 LH = p*L*Ch*wind*B(q(Ta) - q(Ts)) 22 G = k*(Tgnd - Ts)/dz 23 24 Fnet = R + SH + LH + G 25 26 du/dt = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx) 27 dv/dt = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy) 28 dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts) 29 = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy) 30 dp/dt = -Div([u*p,v*p]) 31 = - u*dp/dx - v*dp/dy 32 dTa/dt = Fnet/Cp 33 34 Equation of State: 35 36 P = p*R*Ts 37 38 ----------------------------------------------------------------------- 39 40 Program considers the evolution of a two dimensional atmosphere from 41 sunset to sunrise. There are two components: 42 1. Surface energy balance model to compute diabatic dT (Fnet) 43 2. Dynamical model using simplified primitive equations 44 45 Program is to be initiated at sunset and run to sunrise. 46 47 Inputs are: 48 Surface temperature 49 Dew point temperature 50 Air temperature 51 Temperature at cloud base (if clouds are present) 52 Fraction of sky covered by clouds 53 Wind speed 54 Precipitable water in centimeters 55 Wind direction 56 57 Inputs are are read in from the text file ex5_control.txt. To change an 58 input value use ex5_control.txt. 59 60 Solvers: 61 Backward Euler = default solver 62 Sundials = fastest and most accurate, requires Sundials libraries 63 64 This model is under development and should be used only as an example 65 and not as a predictive weather model. 66 */ 67 68 #include <petscts.h> 69 #include <petscdm.h> 70 #include <petscdmda.h> 71 72 /* stefan-boltzmann constant */ 73 #define SIG 0.000000056703 74 /* absorption-emission constant for surface */ 75 #define EMMSFC 1 76 /* amount of time (seconds) that passes before new flux is calculated */ 77 #define TIMESTEP 1 78 79 /* variables of interest to be solved at each grid point */ 80 typedef struct { 81 PetscScalar Ts,Ta; /* surface and air temperature */ 82 PetscScalar u,v; /* wind speed */ 83 PetscScalar p; /* density */ 84 } Field; 85 86 /* User defined variables. Used in solving for variables of interest */ 87 typedef struct { 88 DM da; /* grid */ 89 PetscScalar csoil; /* heat constant for layer */ 90 PetscScalar dzlay; /* thickness of top soil layer */ 91 PetscScalar emma; /* emission parameter */ 92 PetscScalar wind; /* wind speed */ 93 PetscScalar dewtemp; /* dew point temperature (moisture in air) */ 94 PetscScalar pressure1; /* sea level pressure */ 95 PetscScalar airtemp; /* temperature of air near boundary layer inversion */ 96 PetscScalar Ts; /* temperature at the surface */ 97 PetscScalar fract; /* fraction of sky covered by clouds */ 98 PetscScalar Tc; /* temperature at base of lowest cloud layer */ 99 PetscScalar lat; /* Latitude in degrees */ 100 PetscScalar init; /* initialization scenario */ 101 PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */ 102 } AppCtx; 103 104 /* Struct for visualization */ 105 typedef struct { 106 PetscBool drawcontours; /* flag - 1 indicates drawing contours */ 107 PetscViewer drawviewer; 108 PetscInt interval; 109 } MonitorCtx; 110 111 /* Inputs read in from text file */ 112 struct in { 113 PetscScalar Ts; /* surface temperature */ 114 PetscScalar Td; /* dewpoint temperature */ 115 PetscScalar Tc; /* temperature of cloud base */ 116 PetscScalar fr; /* fraction of sky covered by clouds */ 117 PetscScalar wnd; /* wind speed */ 118 PetscScalar Ta; /* air temperature */ 119 PetscScalar pwt; /* precipitable water */ 120 PetscScalar wndDir; /* wind direction */ 121 PetscScalar lat; /* latitude */ 122 PetscReal time; /* time in hours */ 123 PetscScalar init; 124 }; 125 126 /* functions */ 127 extern PetscScalar emission(PetscScalar); /* sets emission/absorption constant depending on water vapor content */ 128 extern PetscScalar calc_q(PetscScalar); /* calculates specific humidity */ 129 extern PetscScalar mph2mpers(PetscScalar); /* converts miles per hour to meters per second */ 130 extern PetscScalar Lconst(PetscScalar); /* calculates latent heat constant taken from Satellite estimates of wind speed and latent heat flux over the global oceans., Bentamy et al. */ 131 extern PetscScalar fahr_to_cel(PetscScalar); /* converts Fahrenheit to Celsius */ 132 extern PetscScalar cel_to_fahr(PetscScalar); /* converts Celsius to Fahrenheit */ 133 extern PetscScalar calcmixingr(PetscScalar, PetscScalar); /* calculates mixing ratio */ 134 extern PetscScalar cloud(PetscScalar); /* cloud radiative parameterization */ 135 extern PetscErrorCode FormInitialSolution(DM,Vec,void*); /* Specifies initial conditions for the system of equations (PETSc defined function) */ 136 extern PetscErrorCode RhsFunc(TS,PetscReal,Vec,Vec,void*); /* Specifies the user defined functions (PETSc defined function) */ 137 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); /* Specifies output and visualization tools (PETSc defined function) */ 138 extern PetscErrorCode readinput(struct in *put); /* reads input from text file */ 139 extern PetscErrorCode calcfluxs(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); /* calculates upward IR from surface */ 140 extern PetscErrorCode calcfluxa(PetscScalar, PetscScalar, PetscScalar, PetscScalar*); /* calculates downward IR from atmosphere */ 141 extern PetscErrorCode sensibleflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar*); /* calculates sensible heat flux */ 142 extern PetscErrorCode potential_temperature(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); /* calculates potential temperature */ 143 extern PetscErrorCode latentflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); /* calculates latent heat flux */ 144 extern PetscErrorCode calc_gflux(PetscScalar, PetscScalar, PetscScalar*); /* calculates flux between top soil layer and underlying earth */ 145 146 int main(int argc,char **argv) 147 { 148 PetscInt time; /* amount of loops */ 149 struct in put; 150 PetscScalar rh; /* relative humidity */ 151 PetscScalar x; /* memory varialbe for relative humidity calculation */ 152 PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */ 153 PetscScalar emma; /* absorption-emission constant for air */ 154 PetscScalar pressure1 = 101300; /* surface pressure */ 155 PetscScalar mixratio; /* mixing ratio */ 156 PetscScalar airtemp; /* temperature of air near boundary layer inversion */ 157 PetscScalar dewtemp; /* dew point temperature */ 158 PetscScalar sfctemp; /* temperature at surface */ 159 PetscScalar pwat; /* total column precipitable water */ 160 PetscScalar cloudTemp; /* temperature at base of cloud */ 161 AppCtx user; /* user-defined work context */ 162 MonitorCtx usermonitor; /* user-defined monitor context */ 163 TS ts; 164 SNES snes; 165 DM da; 166 Vec T,rhs; /* solution vector */ 167 Mat J; /* Jacobian matrix */ 168 PetscReal ftime,dt; 169 PetscInt steps,dof = 5; 170 PetscBool use_coloring = PETSC_TRUE; 171 MatFDColoring matfdcoloring = 0; 172 PetscBool monitor_off = PETSC_FALSE; 173 174 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 175 176 /* Inputs */ 177 PetscCall(readinput(&put)); 178 179 sfctemp = put.Ts; 180 dewtemp = put.Td; 181 cloudTemp = put.Tc; 182 airtemp = put.Ta; 183 pwat = put.pwt; 184 185 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial Temperature = %g\n",(double)sfctemp)); /* input surface temperature */ 186 187 deep_grnd_temp = sfctemp - 10; /* set underlying ground layer temperature */ 188 emma = emission(pwat); /* accounts for radiative effects of water vapor */ 189 190 /* Converts from Fahrenheit to Celsuis */ 191 sfctemp = fahr_to_cel(sfctemp); 192 airtemp = fahr_to_cel(airtemp); 193 dewtemp = fahr_to_cel(dewtemp); 194 cloudTemp = fahr_to_cel(cloudTemp); 195 deep_grnd_temp = fahr_to_cel(deep_grnd_temp); 196 197 /* Converts from Celsius to Kelvin */ 198 sfctemp += 273; 199 airtemp += 273; 200 dewtemp += 273; 201 cloudTemp += 273; 202 deep_grnd_temp += 273; 203 204 /* Calculates initial relative humidity */ 205 x = calcmixingr(dewtemp,pressure1); 206 mixratio = calcmixingr(sfctemp,pressure1); 207 rh = (x/mixratio)*100; 208 209 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial RH = %.1f percent\n\n",(double)rh)); /* prints initial relative humidity */ 210 211 time = 3600*put.time; /* sets amount of timesteps to run model */ 212 213 /* Configure PETSc TS solver */ 214 /*------------------------------------------*/ 215 216 /* Create grid */ 217 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,20,20,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&da)); 218 PetscCall(DMSetFromOptions(da)); 219 PetscCall(DMSetUp(da)); 220 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 221 222 /* Define output window for each variable of interest */ 223 PetscCall(DMDASetFieldName(da,0,"Ts")); 224 PetscCall(DMDASetFieldName(da,1,"Ta")); 225 PetscCall(DMDASetFieldName(da,2,"u")); 226 PetscCall(DMDASetFieldName(da,3,"v")); 227 PetscCall(DMDASetFieldName(da,4,"p")); 228 229 /* set values for appctx */ 230 user.da = da; 231 user.Ts = sfctemp; 232 user.fract = put.fr; /* fraction of sky covered by clouds */ 233 user.dewtemp = dewtemp; /* dew point temperature (mositure in air) */ 234 user.csoil = 2000000; /* heat constant for layer */ 235 user.dzlay = 0.08; /* thickness of top soil layer */ 236 user.emma = emma; /* emission parameter */ 237 user.wind = put.wnd; /* wind spped */ 238 user.pressure1 = pressure1; /* sea level pressure */ 239 user.airtemp = airtemp; /* temperature of air near boundar layer inversion */ 240 user.Tc = cloudTemp; /* temperature at base of lowest cloud layer */ 241 user.init = put.init; /* user chosen initiation scenario */ 242 user.lat = 70*0.0174532; /* converts latitude degrees to latitude in radians */ 243 user.deep_grnd_temp = deep_grnd_temp; /* temp in lowest ground layer */ 244 245 /* set values for MonitorCtx */ 246 usermonitor.drawcontours = PETSC_FALSE; 247 PetscCall(PetscOptionsHasName(NULL,NULL,"-drawcontours",&usermonitor.drawcontours)); 248 if (usermonitor.drawcontours) { 249 PetscReal bounds[] = {1000.0,-1000., -1000.,-1000., 1000.,-1000., 1000.,-1000., 1000,-1000, 100700,100800}; 250 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer)); 251 PetscCall(PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds)); 252 } 253 usermonitor.interval = 1; 254 PetscCall(PetscOptionsGetInt(NULL,NULL,"-monitor_interval",&usermonitor.interval,NULL)); 255 256 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 257 Extract global vectors from DA; 258 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 259 PetscCall(DMCreateGlobalVector(da,&T)); 260 PetscCall(VecDuplicate(T,&rhs)); /* r: vector to put the computed right hand side */ 261 262 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 263 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 264 PetscCall(TSSetType(ts,TSBEULER)); 265 PetscCall(TSSetRHSFunction(ts,rhs,RhsFunc,&user)); 266 267 /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */ 268 PetscCall(DMSetMatType(da,MATAIJ)); 269 PetscCall(DMCreateMatrix(da,&J)); 270 PetscCall(TSGetSNES(ts,&snes)); 271 if (use_coloring) { 272 ISColoring iscoloring; 273 PetscCall(DMCreateColoring(da,IS_COLORING_GLOBAL,&iscoloring)); 274 PetscCall(MatFDColoringCreate(J,iscoloring,&matfdcoloring)); 275 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 276 PetscCall(MatFDColoringSetUp(J,iscoloring,matfdcoloring)); 277 PetscCall(ISColoringDestroy(&iscoloring)); 278 PetscCall(MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts)); 279 PetscCall(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring)); 280 } else { 281 PetscCall(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL)); 282 } 283 284 /* Define what to print for ts_monitor option */ 285 PetscCall(PetscOptionsHasName(NULL,NULL,"-monitor_off",&monitor_off)); 286 if (!monitor_off) { 287 PetscCall(TSMonitorSet(ts,Monitor,&usermonitor,NULL)); 288 } 289 PetscCall(FormInitialSolution(da,T,&user)); 290 dt = TIMESTEP; /* initial time step */ 291 ftime = TIMESTEP*time; 292 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"time %" PetscInt_FMT ", ftime %g hour, TIMESTEP %g\n",time,(double)(ftime/3600),(double)dt)); 293 294 PetscCall(TSSetTimeStep(ts,dt)); 295 PetscCall(TSSetMaxSteps(ts,time)); 296 PetscCall(TSSetMaxTime(ts,ftime)); 297 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 298 PetscCall(TSSetSolution(ts,T)); 299 PetscCall(TSSetDM(ts,da)); 300 301 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 302 Set runtime options 303 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 304 PetscCall(TSSetFromOptions(ts)); 305 306 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 307 Solve nonlinear system 308 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 309 PetscCall(TSSolve(ts,T)); 310 PetscCall(TSGetSolveTime(ts,&ftime)); 311 PetscCall(TSGetStepNumber(ts,&steps)); 312 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %" PetscInt_FMT " steps\n",(double)(ftime/3600),steps)); 313 314 if (matfdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 315 if (usermonitor.drawcontours) { 316 PetscCall(PetscViewerDestroy(&usermonitor.drawviewer)); 317 } 318 PetscCall(MatDestroy(&J)); 319 PetscCall(VecDestroy(&T)); 320 PetscCall(VecDestroy(&rhs)); 321 PetscCall(TSDestroy(&ts)); 322 PetscCall(DMDestroy(&da)); 323 324 PetscCall(PetscFinalize()); 325 return 0; 326 } 327 /*****************************end main program********************************/ 328 /*****************************************************************************/ 329 /*****************************************************************************/ 330 /*****************************************************************************/ 331 PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux) 332 { 333 PetscFunctionBeginUser; 334 *flux = SIG*((EMMSFC*emma*PetscPowScalarInt(airtemp,4)) + (EMMSFC*fract*(1 - emma)*PetscPowScalarInt(cloudTemp,4)) - (EMMSFC*PetscPowScalarInt(sfctemp,4))); /* calculates flux using Stefan-Boltzmann relation */ 335 PetscFunctionReturn(0); 336 } 337 338 PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */ 339 { 340 PetscScalar emm = 0.001; 341 342 PetscFunctionBeginUser; 343 *flux = SIG*(-emm*(PetscPowScalarInt(airtemp,4))); /* calculates flux usinge Stefan-Boltzmann relation */ 344 PetscFunctionReturn(0); 345 } 346 PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar *sheat) 347 { 348 PetscScalar density = 1; /* air density */ 349 PetscScalar Cp = 1005; /* heat capicity for dry air */ 350 PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 351 352 PetscFunctionBeginUser; 353 wndmix = 0.0025 + 0.0042*wind; /* regression equation valid for neutral and stable BL */ 354 *sheat = density*Cp*wndmix*(airtemp - sfctemp); /* calculates sensible heat flux */ 355 PetscFunctionReturn(0); 356 } 357 358 PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar *latentheat) 359 { 360 PetscScalar density = 1; /* density of dry air */ 361 PetscScalar q; /* actual specific humitity */ 362 PetscScalar qs; /* saturation specific humidity */ 363 PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 364 PetscScalar beta = .4; /* moisture availability */ 365 PetscScalar mr; /* mixing ratio */ 366 PetscScalar lhcnst; /* latent heat of vaporization constant = 2501000 J/kg at 0c */ 367 /* latent heat of saturation const = 2834000 J/kg */ 368 /* latent heat of fusion const = 333700 J/kg */ 369 370 PetscFunctionBeginUser; 371 wind = mph2mpers(wind); /* converts wind from mph to meters per second */ 372 wndmix = 0.0025 + 0.0042*wind; /* regression equation valid for neutral BL */ 373 lhcnst = Lconst(sfctemp); /* calculates latent heat of evaporation */ 374 mr = calcmixingr(sfctemp,pressure1); /* calculates saturation mixing ratio */ 375 qs = calc_q(mr); /* calculates saturation specific humidty */ 376 mr = calcmixingr(dewtemp,pressure1); /* calculates mixing ratio */ 377 q = calc_q(mr); /* calculates specific humidty */ 378 379 *latentheat = density*wndmix*beta*lhcnst*(q - qs); /* calculates latent heat flux */ 380 PetscFunctionReturn(0); 381 } 382 383 PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar *pottemp) 384 { 385 PetscScalar kdry; /* poisson constant for dry atmosphere */ 386 PetscScalar pavg; /* average atmospheric pressure */ 387 /* PetscScalar mixratio; mixing ratio */ 388 /* PetscScalar kmoist; poisson constant for moist atmosphere */ 389 390 PetscFunctionBeginUser; 391 /* mixratio = calcmixingr(sfctemp,pressure1); */ 392 393 /* initialize poisson constant */ 394 kdry = 0.2854; 395 /* kmoist = 0.2854*(1 - 0.24*mixratio); */ 396 397 pavg = ((0.7*pressure1)+pressure2)/2; /* calculates simple average press */ 398 *pottemp = temp*(PetscPowScalar((pressure1/pavg),kdry)); /* calculates potential temperature */ 399 PetscFunctionReturn(0); 400 } 401 extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1) 402 { 403 PetscScalar e; /* vapor pressure */ 404 PetscScalar mixratio; /* mixing ratio */ 405 406 dtemp = dtemp - 273; /* converts from Kelvin to Celsuis */ 407 e = 6.11*(PetscPowScalar(10,((7.5*dtemp)/(237.7+dtemp)))); /* converts from dew point temp to vapor pressure */ 408 e = e*100; /* converts from hPa to Pa */ 409 mixratio = (0.622*e)/(pressure1 - e); /* computes mixing ratio */ 410 mixratio = mixratio*1; /* convert to g/Kg */ 411 412 return mixratio; 413 } 414 extern PetscScalar calc_q(PetscScalar rv) 415 { 416 PetscScalar specific_humidity; /* define specific humidity variable */ 417 specific_humidity = rv/(1 + rv); /* calculates specific humidity */ 418 return specific_humidity; 419 } 420 421 PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar* Gflux) 422 { 423 PetscScalar k; /* thermal conductivity parameter */ 424 PetscScalar n = 0.38; /* value of soil porosity */ 425 PetscScalar dz = 1; /* depth of layer between soil surface and deep soil layer */ 426 PetscScalar unit_soil_weight = 2700; /* unit soil weight in kg/m^3 */ 427 428 PetscFunctionBeginUser; 429 k = ((0.135*(1-n)*unit_soil_weight) + 64.7)/(unit_soil_weight - (0.947*(1-n)*unit_soil_weight)); /* dry soil conductivity */ 430 *Gflux = (k*(deep_grnd_temp - sfctemp)/dz); /* calculates flux from deep ground layer */ 431 PetscFunctionReturn(0); 432 } 433 extern PetscScalar emission(PetscScalar pwat) 434 { 435 PetscScalar emma; 436 437 emma = 0.725 + 0.17*PetscLog10Real(PetscRealPart(pwat)); 438 439 return emma; 440 } 441 extern PetscScalar cloud(PetscScalar fract) 442 { 443 PetscScalar emma = 0; 444 445 /* modifies radiative balance depending on cloud cover */ 446 if (fract >= 0.9) emma = 1; 447 else if (0.9 > fract && fract >= 0.8) emma = 0.9; 448 else if (0.8 > fract && fract >= 0.7) emma = 0.85; 449 else if (0.7 > fract && fract >= 0.6) emma = 0.75; 450 else if (0.6 > fract && fract >= 0.5) emma = 0.65; 451 else if (0.4 > fract && fract >= 0.3) emma = emma*1.086956; 452 return emma; 453 } 454 extern PetscScalar Lconst(PetscScalar sfctemp) 455 { 456 PetscScalar Lheat; 457 sfctemp -=273; /* converts from kelvin to celsius */ 458 Lheat = 4186.8*(597.31 - 0.5625*sfctemp); /* calculates latent heat constant */ 459 return Lheat; 460 } 461 extern PetscScalar mph2mpers(PetscScalar wind) 462 { 463 wind = ((wind*1.6*1000)/3600); /* converts wind from mph to meters per second */ 464 return wind; 465 } 466 extern PetscScalar fahr_to_cel(PetscScalar temp) 467 { 468 temp = (5*(temp-32))/9; /* converts from farhrenheit to celsuis */ 469 return temp; 470 } 471 extern PetscScalar cel_to_fahr(PetscScalar temp) 472 { 473 temp = ((temp*9)/5) + 32; /* converts from celsuis to farhrenheit */ 474 return temp; 475 } 476 PetscErrorCode readinput(struct in *put) 477 { 478 int i; 479 char x; 480 FILE *ifp; 481 double tmp; 482 483 PetscFunctionBegin; 484 ifp = fopen("ex5_control.txt", "r"); 485 PetscCheck(ifp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Unable to open input file"); 486 for (i=0; i<110; i++) { PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 487 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 488 put->Ts = tmp; 489 490 for (i=0; i<43; i++) { PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 491 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 492 put->Td = tmp; 493 494 for (i=0; i<43; i++) { PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 495 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 496 put->Ta = tmp; 497 498 for (i=0; i<43; i++) { PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 499 PetscCheck(fscanf(ifp, "%lf", &tmp)== 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 500 put->Tc = tmp; 501 502 for (i=0; i<43; i++) { PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 503 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 504 put->fr = tmp; 505 506 for (i=0; i<43; i++) {PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 507 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 508 put->wnd = tmp; 509 510 for (i=0; i<43; i++) {PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 511 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 512 put->pwt = tmp; 513 514 for (i=0; i<43; i++) {PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 515 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 516 put->wndDir = tmp; 517 518 for (i=0; i<43; i++) {PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 519 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 520 put->time = tmp; 521 522 for (i=0; i<63; i++) {PetscCheck(fscanf(ifp, "%c", &x) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");} 523 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file"); 524 put->init = tmp; 525 PetscFunctionReturn(0); 526 } 527 528 /* ------------------------------------------------------------------- */ 529 PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx) 530 { 531 AppCtx *user = (AppCtx*)ctx; /* user-defined application context */ 532 PetscInt i,j,xs,ys,xm,ym,Mx,My; 533 Field **X; 534 535 PetscFunctionBeginUser; 536 PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 537 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 538 539 /* Get pointers to vector data */ 540 PetscCall(DMDAVecGetArray(da,Xglobal,&X)); 541 542 /* Get local grid boundaries */ 543 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 544 545 /* Compute function over the locally owned part of the grid */ 546 547 if (user->init == 1) { 548 for (j=ys; j<ys+ym; j++) { 549 for (i=xs; i<xs+xm; i++) { 550 X[j][i].Ts = user->Ts - i*0.0001; 551 X[j][i].Ta = X[j][i].Ts - 5; 552 X[j][i].u = 0; 553 X[j][i].v = 0; 554 X[j][i].p = 1.25; 555 if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001; 556 if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001; 557 } 558 } 559 } else { 560 for (j=ys; j<ys+ym; j++) { 561 for (i=xs; i<xs+xm; i++) { 562 X[j][i].Ts = user->Ts; 563 X[j][i].Ta = X[j][i].Ts - 5; 564 X[j][i].u = 0; 565 X[j][i].v = 0; 566 X[j][i].p = 1.25; 567 } 568 } 569 } 570 571 /* Restore vectors */ 572 PetscCall(DMDAVecRestoreArray(da,Xglobal,&X)); 573 PetscFunctionReturn(0); 574 } 575 576 /* 577 RhsFunc - Evaluates nonlinear function F(u). 578 579 Input Parameters: 580 . ts - the TS context 581 . t - current time 582 . Xglobal - input vector 583 . F - output vector 584 . ptr - optional user-defined context, as set by SNESSetFunction() 585 586 Output Parameter: 587 . F - rhs function vector 588 */ 589 PetscErrorCode RhsFunc(TS ts,PetscReal t,Vec Xglobal,Vec F,void *ctx) 590 { 591 AppCtx *user = (AppCtx*)ctx; /* user-defined application context */ 592 DM da = user->da; 593 PetscInt i,j,Mx,My,xs,ys,xm,ym; 594 PetscReal dhx,dhy; 595 Vec localT; 596 Field **X,**Frhs; /* structures that contain variables of interest and left hand side of governing equations respectively */ 597 PetscScalar csoil = user->csoil; /* heat constant for layer */ 598 PetscScalar dzlay = user->dzlay; /* thickness of top soil layer */ 599 PetscScalar emma = user->emma; /* emission parameter */ 600 PetscScalar wind = user->wind; /* wind speed */ 601 PetscScalar dewtemp = user->dewtemp; /* dew point temperature (moisture in air) */ 602 PetscScalar pressure1 = user->pressure1; /* sea level pressure */ 603 PetscScalar airtemp = user->airtemp; /* temperature of air near boundary layer inversion */ 604 PetscScalar fract = user->fract; /* fraction of the sky covered by clouds */ 605 PetscScalar Tc = user->Tc; /* temperature at base of lowest cloud layer */ 606 PetscScalar lat = user->lat; /* latitude */ 607 PetscScalar Cp = 1005.7; /* specific heat of air at constant pressure */ 608 PetscScalar Rd = 287.058; /* gas constant for dry air */ 609 PetscScalar diffconst = 1000; /* diffusion coefficient */ 610 PetscScalar f = 2*0.0000727*PetscSinScalar(lat); /* coriolis force */ 611 PetscScalar deep_grnd_temp = user->deep_grnd_temp; /* temp in lowest ground layer */ 612 PetscScalar Ts,u,v,p; 613 PetscScalar u_abs,u_plus,u_minus,v_abs,v_plus,v_minus; 614 615 PetscScalar sfctemp1,fsfc1,Ra; 616 PetscScalar sheat; /* sensible heat flux */ 617 PetscScalar latentheat; /* latent heat flux */ 618 PetscScalar groundflux; /* flux from conduction of deep ground layer in contact with top soil */ 619 PetscInt xend,yend; 620 621 PetscFunctionBeginUser; 622 PetscCall(DMGetLocalVector(da,&localT)); 623 PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 624 625 dhx = (PetscReal)(Mx-1)/(5000*(Mx-1)); /* dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5] */ 626 dhy = (PetscReal)(My-1)/(5000*(Mx-1)); /* dhy = 1/dy; */ 627 628 /* 629 Scatter ghost points to local vector,using the 2-step process 630 DAGlobalToLocalBegin(),DAGlobalToLocalEnd(). 631 By placing code between these two statements, computations can be 632 done while messages are in transition. 633 */ 634 PetscCall(DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT)); 635 PetscCall(DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT)); 636 637 /* Get pointers to vector data */ 638 PetscCall(DMDAVecGetArrayRead(da,localT,&X)); 639 PetscCall(DMDAVecGetArray(da,F,&Frhs)); 640 641 /* Get local grid boundaries */ 642 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 643 644 /* Compute function over the locally owned part of the grid */ 645 /* the interior points */ 646 xend=xs+xm; yend=ys+ym; 647 for (j=ys; j<yend; j++) { 648 for (i=xs; i<xend; i++) { 649 Ts = X[j][i].Ts; u = X[j][i].u; v = X[j][i].v; p = X[j][i].p; /*P = X[j][i].P; */ 650 651 sfctemp1 = (double)Ts; 652 PetscCall(calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1)); /* calculates surface net radiative flux */ 653 PetscCall(sensibleflux(sfctemp1,airtemp,wind,&sheat)); /* calculate sensible heat flux */ 654 PetscCall(latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat)); /* calculates latent heat flux */ 655 PetscCall(calc_gflux(sfctemp1,deep_grnd_temp,&groundflux)); /* calculates flux from earth below surface soil layer by conduction */ 656 PetscCall(calcfluxa(sfctemp1,airtemp,emma,&Ra)); /* Calculates the change in downward radiative flux */ 657 fsfc1 = fsfc1 + latentheat + sheat + groundflux; /* adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux */ 658 659 /* convective coefficients for upwinding */ 660 u_abs = PetscAbsScalar(u); 661 u_plus = .5*(u + u_abs); /* u if u>0; 0 if u<0 */ 662 u_minus = .5*(u - u_abs); /* u if u <0; 0 if u>0 */ 663 664 v_abs = PetscAbsScalar(v); 665 v_plus = .5*(v + v_abs); /* v if v>0; 0 if v<0 */ 666 v_minus = .5*(v - v_abs); /* v if v <0; 0 if v>0 */ 667 668 /* Solve governing equations */ 669 /* P = p*Rd*Ts; */ 670 671 /* du/dt -> time change of east-west component of the wind */ 672 Frhs[j][i].u = - u_plus*(u - X[j][i-1].u)*dhx - u_minus*(X[j][i+1].u - u)*dhx /* - u(du/dx) */ 673 - v_plus*(u - X[j-1][i].u)*dhy - v_minus*(X[j+1][i].u - u)*dhy /* - v(du/dy) */ 674 -(Rd/p)*(Ts*(X[j][i+1].p - X[j][i-1].p)*0.5*dhx + p*0*(X[j][i+1].Ts - X[j][i-1].Ts)*0.5*dhx) /* -(R/p)[Ts(dp/dx)+ p(dTs/dx)] */ 675 /* -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx */ 676 + f*v; 677 678 /* dv/dt -> time change of north-south component of the wind */ 679 Frhs[j][i].v = - u_plus*(v - X[j][i-1].v)*dhx - u_minus*(X[j][i+1].v - v)*dhx /* - u(dv/dx) */ 680 - v_plus*(v - X[j-1][i].v)*dhy - v_minus*(X[j+1][i].v - v)*dhy /* - v(dv/dy) */ 681 -(Rd/p)*(Ts*(X[j+1][i].p - X[j-1][i].p)*0.5*dhy + p*0*(X[j+1][i].Ts - X[j-1][i].Ts)*0.5*dhy) /* -(R/p)[Ts(dp/dy)+ p(dTs/dy)] */ 682 /* -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy */ 683 -f*u; 684 685 /* dT/dt -> time change of temperature */ 686 Frhs[j][i].Ts = (fsfc1/(csoil*dzlay)) /* Fnet/(Cp*dz) diabatic change in T */ 687 -u_plus*(Ts - X[j][i-1].Ts)*dhx - u_minus*(X[j][i+1].Ts - Ts)*dhx /* - u*(dTs/dx) advection x */ 688 -v_plus*(Ts - X[j-1][i].Ts)*dhy - v_minus*(X[j+1][i].Ts - Ts)*dhy /* - v*(dTs/dy) advection y */ 689 + diffconst*((X[j][i+1].Ts - 2*Ts + X[j][i-1].Ts)*dhx*dhx /* + D(Ts_xx + Ts_yy) diffusion */ 690 + (X[j+1][i].Ts - 2*Ts + X[j-1][i].Ts)*dhy*dhy); 691 692 /* dp/dt -> time change of */ 693 Frhs[j][i].p = -u_plus*(p - X[j][i-1].p)*dhx - u_minus*(X[j][i+1].p - p)*dhx /* - u*(dp/dx) */ 694 -v_plus*(p - X[j-1][i].p)*dhy - v_minus*(X[j+1][i].p - p)*dhy; /* - v*(dp/dy) */ 695 696 Frhs[j][i].Ta = Ra/Cp; /* dTa/dt time change of air temperature */ 697 } 698 } 699 700 /* Restore vectors */ 701 PetscCall(DMDAVecRestoreArrayRead(da,localT,&X)); 702 PetscCall(DMDAVecRestoreArray(da,F,&Frhs)); 703 PetscCall(DMRestoreLocalVector(da,&localT)); 704 PetscFunctionReturn(0); 705 } 706 707 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec T,void *ctx) 708 { 709 const PetscScalar *array; 710 MonitorCtx *user = (MonitorCtx*)ctx; 711 PetscViewer viewer = user->drawviewer; 712 PetscReal norm; 713 714 PetscFunctionBeginUser; 715 PetscCall(VecNorm(T,NORM_INFINITY,&norm)); 716 717 if (step%user->interval == 0) { 718 PetscCall(VecGetArrayRead(T,&array)); 719 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"step %" PetscInt_FMT ", time %8.1f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n",step,(double)time,(double)(((array[0]-273)*9)/5 + 32),(double)(((array[1]-273)*9)/5 + 32),(double)array[2],(double)array[3],(double)array[4],(double)array[5])); 720 PetscCall(VecRestoreArrayRead(T,&array)); 721 } 722 723 if (user->drawcontours) { 724 PetscCall(VecView(T,viewer)); 725 } 726 PetscFunctionReturn(0); 727 } 728 729 /*TEST 730 731 build: 732 requires: !complex !single 733 734 test: 735 args: -ts_max_steps 130 -monitor_interval 60 736 output_file: output/ex5.out 737 requires: !complex !single 738 localrunfiles: ex5_control.txt 739 740 test: 741 suffix: 2 742 nsize: 4 743 args: -ts_max_steps 130 -monitor_interval 60 744 output_file: output/ex5.out 745 localrunfiles: ex5_control.txt 746 requires: !complex !single 747 748 TEST*/ 749