00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <math.h>
00004 #include "globalt.h"
00005 #include "globalc.h"
00006 #include "glasgowmatc.h"
00007
00008
00009
00010
00011
00012
00013
00014 glasgowmatc::glasgowmatc ()
00015 {
00016 model=2;
00017 ra=287.0;
00018 rv=461.5;
00019 stef=5.67e-8;
00020
00021 rhoc = 0.0;
00022 rhos = 0.0;
00023 por0 = 0.0;
00024 k0 = 0.0;
00025 rhol0 = 0.0;
00026 sssp = 0.0;
00027 t0 = 0.0;
00028 rhov0 = 0.0;
00029 pginf = 0.0;
00030 emmi = 0.0;
00031 alph = 0.0;
00032 hq = 0.0;
00033 crhoair = 0.0;
00034 tfirestart = 0.0;
00035 }
00036
00037 glasgowmatc::~glasgowmatc ()
00038 {}
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void glasgowmatc::matcond_u (matrix &d,long ri,long ci,long ipp)
00058 {
00059 long m;
00060
00061 m = d.m;
00062
00063 switch (m){
00064 case 1:{
00065 if (ci == 2){
00066 matrix s(d.m,d.m),sd(d.m,d.n);
00067
00068 matcond1d_u (sd,ri,ci,ipp);
00069
00070 Cmu->matstiff (s,ipp);
00071 mxm(s,sd,d);
00072
00073 destrm (s); destrm (sd);
00074 }
00075 else
00076 matcond1d_u (d,ri,ci,ipp);
00077 break;
00078 }
00079 case 3:{
00080 if (ci == 2){
00081 matrix s(d.m,d.m),sd(d.m,d.n);
00082
00083 matcond2d_u (sd,ri,ci,ipp);
00084
00085 Cmu->matstiff (s,ipp);
00086 mxm(s,sd,d);
00087
00088 destrm (s); destrm (sd);
00089 }
00090 else
00091 matcond2d_u (d,ri,ci,ipp);
00092 break;
00093 }
00094 case 6:{
00095 if (ci == 2){
00096 matrix s(d.m,d.m),sd(d.m,d.n);
00097
00098 matcond3d_u (sd,ri,ci,ipp);
00099
00100 Cmu->matstiff (s,ipp);
00101 mxm(s,sd,d);
00102
00103 destrm (s); destrm (sd);
00104 }
00105 else
00106 matcond3d_u (d,ri,ci,ipp);
00107 break;
00108 }
00109 default:{
00110 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00111 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00112 }
00113 }
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 void glasgowmatc::matcap_u (matrix &d,long ri,long ci,long ipp)
00128 {
00129 long m;
00130
00131 m = d.m;
00132
00133 switch (m){
00134 case 1:{
00135 matcap1d_u (d,ri,ci,ipp);
00136 break;
00137 }
00138 case 3:{
00139 matcap2d_u (d,ri,ci,ipp);
00140 break;
00141 }
00142 case 6:{
00143 matcap3d_u (d,ri,ci,ipp);
00144 break;
00145 }
00146 default:{
00147 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00148 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00149 }
00150 }
00151 }
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 void glasgowmatc::matcond1d_u (matrix &d,long ri,long ci,long ipp)
00166 {
00167 double k;
00168 double rhov,pg,t;
00169 k = 0.0;
00170
00171 rhov = Cmu->ip[ipp].av[0];
00172 pg = Cmu->ip[ipp].av[1];
00173 t = Cmu->ip[ipp].av[2];
00174
00175 if((ri == 0) && (ci == 0))
00176 k = kuv(t,pg,rhov);
00177 if((ri == 0) && (ci == 1))
00178 k = kup(t,pg,rhov);
00179 if((ri == 0) && (ci == 2))
00180 k = kut(t,pg,rhov);
00181
00182 d[0][0] = k;
00183 }
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 void glasgowmatc::matcond2d_u (matrix &d,long ri,long ci,long ipp)
00197 {
00198 double k;
00199 double rhov,pg,t;
00200 k = 0.0;
00201
00202 rhov = Cmu->ip[ipp].av[0];
00203 pg = Cmu->ip[ipp].av[1];
00204 t = Cmu->ip[ipp].av[2];
00205
00206 if((ri == 0) && (ci == 0))
00207 k = kuv(t,pg,rhov);
00208 if((ri == 0) && (ci == 1))
00209 k = kup(t,pg,rhov);
00210 if((ri == 0) && (ci == 2))
00211 k = kut(t,pg,rhov);
00212
00213 fillm(0.0,d);
00214
00215 d[0][0] = k;
00216 d[1][0] = k;
00217 d[2][0] = 0.0;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 void glasgowmatc::matcond3d_u (matrix &d,long ri,long ci,long ipp)
00232 {
00233 double k;
00234 double rhov,pg,t;
00235 k = 0.0;
00236
00237 rhov = Cmu->ip[ipp].av[0];
00238 pg = Cmu->ip[ipp].av[1];
00239 t = Cmu->ip[ipp].av[2];
00240
00241 if((ri == 0) && (ci == 0))
00242 k = kuv(t,pg,rhov);
00243 if((ri == 0) && (ci == 1))
00244 k = kup(t,pg,rhov);
00245 if((ri == 0) && (ci == 2))
00246 k = kut(t,pg,rhov);
00247
00248 fillm(0.0,d);
00249
00250 d[0][0]=k;
00251 d[1][0]=k;
00252 d[2][0]=k;
00253 d[3][0]=0.0;
00254 d[4][0]=0.0;
00255 d[5][0]=0.0;
00256 }
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270 void glasgowmatc::matcap1d_u (matrix &d,long ri,long ci,long ipp)
00271 {
00272 double c;
00273 double rhov,pg,t;
00274 c = 0.0;
00275
00276 rhov = Cmu->ip[ipp].av[0];
00277 pg = Cmu->ip[ipp].av[1];
00278 t = Cmu->ip[ipp].av[2];
00279
00280 if((ri == 0) && (ci == 0))
00281 c = cuv();
00282 if((ri == 0) && (ci == 1))
00283 c = cup();
00284 if((ri == 0) && (ci == 2))
00285 c = cut();
00286
00287 fillm(0.0,d);
00288
00289 d[0][0] = c;
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 void glasgowmatc::matcap2d_u (matrix &d,long ri,long ci,long ipp)
00305 {
00306 double c;
00307 double rhov,pg,t;
00308 c = 0.0;
00309
00310 rhov = Cmu->ip[ipp].av[0];
00311 pg = Cmu->ip[ipp].av[1];
00312 t = Cmu->ip[ipp].av[2];
00313
00314 if((ri == 0) && (ci == 0))
00315 c = cuv();
00316 if((ri == 0) && (ci == 1))
00317 c = cup();
00318 if((ri == 0) && (ci == 2))
00319 c = cut();
00320
00321 fillm(0.0,d);
00322
00323 d[0][0] = c;
00324 d[1][0] = c;
00325 d[2][0] = 0.0;
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 void glasgowmatc::matcap3d_u (matrix &d,long ri,long ci,long ipp)
00341 {
00342 double c;
00343 double rhov,pg,t;
00344 c = 0.0;
00345
00346 rhov = Cmu->ip[ipp].av[0];
00347 pg = Cmu->ip[ipp].av[1];
00348 t = Cmu->ip[ipp].av[2];
00349
00350 if((ri == 0) && (ci == 0))
00351 c = cuv();
00352 if((ri == 0) && (ci == 1))
00353 c = cup();
00354 if((ri == 0) && (ci == 2))
00355 c = cut();
00356
00357 fillm(0.0,d);
00358
00359 d[0][0] = c;
00360 d[1][0] = c;
00361 d[2][0] = c;
00362 d[3][0] = 0.0;
00363 d[4][0] = 0.0;
00364 d[5][0] = 0.0;
00365 }
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 void glasgowmatc::matcond_l (matrix &d,long ri,long ci,long ipp)
00384 {
00385 long m;
00386
00387 m = d.n;
00388
00389 switch (m){
00390 case 1:{
00391 matcond1d_l (d,ri,ci,ipp);
00392 break;
00393 }
00394 case 3:{
00395 matcond2d_l (d,ri,ci,ipp);
00396 break;
00397 }
00398 case 6:{
00399 matcond3d_l (d,ri,ci,ipp);
00400 break;
00401 }
00402 default:{
00403 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00404 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00405 }
00406 }
00407 }
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 void glasgowmatc::matcap_l (matrix &d,long ri,long ci,long ipp)
00421 {
00422 long m;
00423
00424 m = d.n;
00425
00426 switch (m){
00427 case 1:{
00428 matcap1d_l (d,ri,ci,ipp);
00429 break;
00430 }
00431 case 3:{
00432 matcap2d_l (d,ri,ci,ipp);
00433 break;
00434 }
00435 case 6:{
00436 matcap3d_l (d,ri,ci,ipp);
00437 break;
00438 }
00439 default:{
00440 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00441 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00442 }
00443 }
00444 }
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 void glasgowmatc::matcond1d_l (matrix &d,long ri,long ci,long ipp)
00459 {
00460 double k;
00461 double rhov,pg,t;
00462 k = 0.0;
00463
00464 rhov = Cml->ip[ipp].av[0];
00465 pg = Cml->ip[ipp].av[1];
00466 t = Cml->ip[ipp].av[2];
00467
00468 if((ri == 0) && (ci == 0))
00469 k = kvu();
00470 if((ri == 1) && (ci == 0))
00471 k = kpu();
00472 if((ri == 2) && (ci == 0))
00473 k = ktu();
00474
00475 d[0][0] = k;
00476 }
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489 void glasgowmatc::matcond2d_l (matrix &d,long ri,long ci,long ipp)
00490 {
00491 double k;
00492 double rhov,pg,t;
00493 k = 0.0;
00494
00495 rhov = Cml->ip[ipp].av[0];
00496 pg = Cml->ip[ipp].av[1];
00497 t = Cml->ip[ipp].av[2];
00498
00499 if((ri == 0) && (ci == 0))
00500 k = kvu();
00501 if((ri == 1) && (ci == 0))
00502 k = kpu();
00503 if((ri == 2) && (ci == 0))
00504 k = ktu();
00505
00506 fillm(0.0,d);
00507
00508
00509 d[0][0] = k;
00510 d[0][1] = k;
00511 d[0][2] = 0.0;
00512 }
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 void glasgowmatc::matcond3d_l (matrix &d,long ri,long ci,long ipp)
00526 {
00527 double k;
00528 double rhov,pg,t;
00529 k = 0.0;
00530
00531 rhov = Cml->ip[ipp].av[0];
00532 pg = Cml->ip[ipp].av[1];
00533 t = Cml->ip[ipp].av[2];
00534
00535 if((ri == 0) && (ci == 0))
00536 k = kvu();
00537 if((ri == 1) && (ci == 0))
00538 k = kpu();
00539 if((ri == 2) && (ci == 0))
00540 k = ktu();
00541
00542 fillm(0.0,d);
00543
00544 d[0][0]=k;
00545 d[0][1]=k;
00546 d[0][2]=k;
00547 d[0][3]=0.0;
00548 d[0][4]=0.0;
00549 d[0][5]=0.0;
00550 }
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 void glasgowmatc::matcap1d_l (matrix &d,long ri,long ci,long ipp)
00565 {
00566 double c;
00567 double rhov,pg,t;
00568 c = 0.0;
00569
00570 rhov = Cml->ip[ipp].av[0];
00571 pg = Cml->ip[ipp].av[1];
00572 t = Cml->ip[ipp].av[2];
00573
00574 if((ri == 0) && (ci == 0))
00575 c = cvu();
00576 if((ri == 1) && (ci == 0))
00577 c = cpu();
00578 if((ri == 2) && (ci == 0))
00579 c = ctu();
00580
00581 fillm(0.0,d);
00582
00583 d[0][0] = c;
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 void glasgowmatc::matcap2d_l (matrix &d,long ri,long ci,long ipp)
00599 {
00600 double c;
00601 double rhov,pg,t;
00602 c = 0.0;
00603
00604 rhov = Cml->ip[ipp].av[0];
00605 pg = Cml->ip[ipp].av[1];
00606 t = Cml->ip[ipp].av[2];
00607
00608 if((ri == 0) && (ci == 0))
00609 c = cvu();
00610 if((ri == 1) && (ci == 0))
00611 c = cpu();
00612 if((ri == 2) && (ci == 0))
00613 c = ctu();
00614
00615 fillm(0.0,d);
00616
00617 d[0][0] = c;
00618 d[0][1] = c;
00619 d[0][2] = 0.0;
00620 }
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634 void glasgowmatc::matcap3d_l (matrix &d,long ri,long ci,long ipp)
00635 {
00636 double c;
00637 double rhov,pg,t;
00638 c = 0.0;
00639
00640 rhov = Cml->ip[ipp].av[0];
00641 pg = Cml->ip[ipp].av[1];
00642 t = Cml->ip[ipp].av[2];
00643
00644 if((ri == 0) && (ci == 0))
00645 c = cvu();
00646 if((ri == 1) && (ci == 0))
00647 c = cpu();
00648 if((ri == 2) && (ci == 0))
00649 c = ctu();
00650
00651 fillm(0.0,d);
00652
00653 d[0][0] = c;
00654 d[0][1] = c;
00655 d[0][2] = c;
00656 d[0][3] = 0.0;
00657 d[0][4] = 0.0;
00658 d[0][5] = 0.0;
00659 }
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677 void glasgowmatc::rhs_u1 (matrix &d,long ri,long ci,long ipp)
00678 {
00679 long m;
00680
00681 m = d.m;
00682
00683 switch (m){
00684 case 1:{
00685 if (ci == 2){
00686 matrix s(d.m,d.m),sd(d.m,d.n);
00687
00688 rhs1d_u1 (sd,ri,ci,ipp);
00689
00690 Cmu->matstiff (s,ipp);
00691 mxm(s,sd,d);
00692
00693 destrm (s); destrm (sd);
00694 }
00695 else
00696 rhs1d_u1 (d,ri,ci,ipp);
00697 break;
00698 }
00699 case 3:{
00700 if (ci == 2){
00701 matrix s(d.m,d.m),sd(d.m,d.n);
00702
00703 rhs2d_u1 (sd,ri,ci,ipp);
00704
00705 Cmu->matstiff (s,ipp);
00706 mxm(s,sd,d);
00707
00708 destrm (s); destrm (sd);
00709 }
00710 else
00711 rhs2d_u1 (d,ri,ci,ipp);
00712 break;
00713 }
00714 case 6:{
00715 if (ci == 2){
00716 matrix s(d.m,d.m),sd(d.m,d.n);
00717
00718 rhs3d_u1 (sd,ri,ci,ipp);
00719
00720 Cmu->matstiff (s,ipp);
00721 mxm(s,sd,d);
00722
00723 destrm (s); destrm (sd);
00724 }
00725 else
00726 rhs3d_u1 (d,ri,ci,ipp);
00727 break;
00728 }
00729 default:{
00730 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00731 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00732 }
00733 }
00734 }
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748 void glasgowmatc::rhs1d_u1 (matrix &d,long ri,long ci,long ipp)
00749 {
00750 double k;
00751 double rhov,pg,t;
00752 k = 0.0;
00753
00754 rhov = Cmu->ip[ipp].av[0];
00755 pg = Cmu->ip[ipp].av[1];
00756 t = Cmu->ip[ipp].av[2];
00757
00758 if((ri == 0) && (ci == 0))
00759 k = fuv1(t,pg,rhov);
00760 if((ri == 0) && (ci == 1))
00761 k = fup1(t,pg,rhov);
00762 if((ri == 0) && (ci == 2))
00763 k = fut1(t,pg,rhov);
00764
00765 d[0][0] = k;
00766 }
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779 void glasgowmatc::rhs2d_u1 (matrix &d,long ri,long ci,long ipp)
00780 {
00781 double k;
00782 double rhov,pg,t;
00783 k = 0.0;
00784
00785 rhov = Cmu->ip[ipp].av[0];
00786 pg = Cmu->ip[ipp].av[1];
00787 t = Cmu->ip[ipp].av[2];
00788
00789 if((ri == 0) && (ci == 0))
00790 k = fuv1(t,pg,rhov);
00791 if((ri == 0) && (ci == 1))
00792 k = fup1(t,pg,rhov);
00793 if((ri == 0) && (ci == 2))
00794 k = fut1(t,pg,rhov);
00795
00796 fillm(0.0,d);
00797
00798 d[0][0] = k;
00799 d[1][0] = k;
00800 d[2][0] = 0.0;
00801 }
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814 void glasgowmatc::rhs3d_u1 (matrix &d,long ri,long ci,long ipp)
00815 {
00816 double k;
00817 double rhov,pg,t;
00818 k = 0.0;
00819
00820 rhov = Cmu->ip[ipp].av[0];
00821 pg = Cmu->ip[ipp].av[1];
00822 t = Cmu->ip[ipp].av[2];
00823
00824 if((ri == 0) && (ci == 0))
00825 k = fuv1(t,pg,rhov);
00826 if((ri == 0) && (ci == 1))
00827 k = fup1(t,pg,rhov);
00828 if((ri == 0) && (ci == 2))
00829 k = fut1(t,pg,rhov);
00830
00831 fillm(0.0,d);
00832
00833 d[0][0]=k;
00834 d[1][0]=k;
00835 d[2][0]=k;
00836 d[3][0]=0.0;
00837 d[4][0]=0.0;
00838 d[5][0]=0.0;
00839 }
00840
00841
00842
00843
00844
00845 void glasgowmatc::read (XFILE *in)
00846 {
00847 xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf",&rhoc,&rhos,&por0,&k0,&rhol0,&sssp,&t0,&rhov0,&pginf);
00848 xfscanf (in,"%lf %lf %lf %lf %lf",&emmi,&alph,&hq,&crhoair,&tfirestart);
00849 }
00850
00851
00852
00853
00854
00855 double glasgowmatc::f_rhovinf ()
00856 {
00857 double rhovinf;
00858
00859 rhovinf=rhov0*0.8;
00860
00861 return(rhovinf);
00862 }
00863
00864
00865
00866
00867
00868
00869 double glasgowmatc::f_tinf (double time)
00870 {
00871 double tinf;
00872
00873
00874
00875
00876
00877 tinf=((t0-273.15)+345.0*log10(8.0*(time-tfirestart)/60.0+1.0))+273.15;
00878
00879 return(tinf);
00880 }
00881
00882
00883
00884
00885
00886 double glasgowmatc::f_pv (double rhov, double t)
00887 {
00888 double pv;
00889
00890 pv=rhov*rv*t;
00891
00892 return(pv);
00893 }
00894
00895
00896
00897
00898
00899 double glasgowmatc::f_psat(double t)
00900 {
00901 double psat;
00902
00903 psat = -0.000000001437422219446870*t*t*t*t*t*t +
00904 0.000004424390583021230*t*t*t*t*t -
00905 0.003928080821257910*t*t*t*t +
00906 1.591032529443030*t*t*t -
00907 325.8874385048470*t*t +
00908 32147.79527519750*t -
00909 1154663.603250870;
00910
00911 return(psat);
00912 }
00913
00914
00915
00916
00917
00918 double glasgowmatc::f_dpsatdt (double t)
00919 {
00920 double dpsatdt;
00921
00922 dpsatdt=32147.7952751975 -
00923 651.774877009694*t +
00924 4.77309758832909*t*t -
00925 0.01571232328503164*t*t*t +
00926 0.00002212195291510615*t*t*t*t -
00927 8.62453331668122e-9*t*t*t*t*t;
00928
00929 return(dpsatdt);
00930 }
00931
00932
00933
00934
00935
00936 double glasgowmatc::f_rhow (double t)
00937 {
00938 double rhow;
00939
00940
00941
00942 rhow=1000.0;
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981 return(rhow);
00982 }
00983
00984
00985
00986
00987
00988
00989 double glasgowmatc::f_rhow0 (double t)
00990 {
00991 double rhow0;
00992
00993
00994
00995 rhow0=1000.0;
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034 return(rhow0);
01035 }
01036
01037
01038
01039
01040
01041 double glasgowmatc::f_drhowdt (double t)
01042 {
01043 double drhowdt;
01044
01045
01046 drhowdt=0.0;
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080 return(drhowdt);
01081 }
01082
01083
01084
01085
01086
01087 double glasgowmatc::f_por (double t)
01088 {
01089 double tc,por,pora,porb,porc,pord;
01090
01091 tc=t-273.15;
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101 pora= (-1.0/85750000.0);
01102 porb= (27.0/1715000.0);
01103 porc= (-24.0/8575.0);
01104 pord= (389.0/343.0);
01105
01106 if(tc < 100.0){
01107 por=por0;
01108 }
01109 if((tc >= 100.0) && (tc <= 800.0)){
01110 por=por0*((pora*tc*tc*tc)+(porb*tc*tc)+(porc*tc)+pord);
01111 }
01112 if(tc > 800.0){
01113 por=por0*3.0;
01114 }
01115
01116 return(por);
01117 }
01118
01119
01120
01121
01122
01123
01124
01125 double glasgowmatc::f_dpordt (double t)
01126 {
01127 double dpordt,tc,pora,porb,porc,pord;
01128
01129 tc=t-273.15;
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139 pora=-(1.0/85750000.0);
01140 porb=(27.0/1715000.0);
01141 porc=-(24.0/8575.0);
01142 pord=(389.0/343.0);
01143
01144 if(tc < 100.0){
01145 dpordt=0.0;
01146 }
01147 if((tc >= 100.0) && (tc <= 800.0)){
01148 dpordt=por0*((3*pora*tc*tc)+(2*porb*tc)+porc);
01149 }
01150 if(tc > 800.0){
01151 dpordt=0.0;
01152 }
01153
01154 return(dpordt);
01155 }
01156
01157
01158
01159
01160
01161
01162 double glasgowmatc::f_fracl0 (double t)
01163 {
01164 double rhow0,fracl0;
01165
01166 rhow0 = f_rhow0(t);
01167
01168 fracl0=rhol0/rhow0;
01169 if(fracl0 > por0)
01170 fracl0=por0;
01171
01172 return(fracl0);
01173 }
01174
01175
01176
01177
01178
01179
01180 double glasgowmatc::f_h (double rhov, double t)
01181 {
01182 double pv,psat,h;
01183
01184 pv = f_pv(rhov,t);
01185 psat = f_psat(t);
01186
01187 h=pv/psat;
01188
01189 return(h);
01190 }
01191
01192
01193
01194
01195
01196
01197 double glasgowmatc::f_m (double t)
01198 {
01199 double m;
01200 double tc;
01201
01202 tc=t-273.15;
01203
01204 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
01205
01206 if (m > 1.0){
01207 m = 1.0;
01208 }
01209
01210 return(m);
01211 }
01212
01213
01214
01215
01216
01217
01218
01219 double glasgowmatc::f_mi (double t)
01220 {
01221 double m,mi;
01222 double tc;
01223
01224 tc=t-273.15;
01225
01226 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
01227 mi=1.0/m;
01228
01229 if (m > 1.0){
01230 m = 1.0;
01231 mi = 1.0;
01232 }
01233
01234 return(mi);
01235 }
01236
01237
01238
01239
01240
01241
01242 double glasgowmatc::f_dmdt (double t)
01243 {
01244 double m,dmdt;
01245 double tc;
01246
01247 tc=t-273.15;
01248
01249 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
01250 dmdt= -1.0*((54635.0*(10.0+tc)))/(27417.5 + 20.0*tc + tc*tc)/(27417.5 + 20.0*tc + tc*tc);
01251
01252
01253 if (m > 1.0){
01254 dmdt = 0.0;
01255 }
01256
01257 return(dmdt);
01258 }
01259
01260
01261
01262
01263
01264
01265 double glasgowmatc::f_dmidt (double t)
01266 {
01267 double m,dmdt,dmidt;
01268 double tc;
01269
01270 tc=t-273.15;
01271
01272 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
01273 dmdt= -1.0*((54635.0*(10.0+tc)))/(27417.5 + 20.0*tc + tc*tc)/(27417.5 + 20.0*tc + tc*tc);
01274
01275 dmidt = -1.0/m/m*dmdt;
01276
01277
01278 if (m > 1.0){
01279 dmidt = 0.0;
01280 }
01281
01282 return(dmidt);
01283 }
01284
01285
01286
01287
01288
01289 double glasgowmatc::f_fracl (double rhov, double t)
01290 {
01291 double fracl0,rhow0,fracl,h;
01292 double aaa,bbb,ccc,ddd;
01293 double mi,rhow,m,pv,psat,por;
01294
01295 h = f_h(rhov,t);
01296 fracl0 = f_fracl0(t);
01297 rhow0 = f_rhow0(t);
01298 mi = f_mi(t);
01299 rhow = f_rhow(t);
01300 m = f_m(t);
01301 pv = f_pv(rhov,t);
01302 psat = f_psat(t);
01303 por = f_por(t);
01304
01305
01306
01307 if(h <= 0.96){
01308
01309 fracl=exp(mi*log((fracl0*rhow0/rhoc)*h))*rhoc/rhow;
01310
01311 }
01312 if((h > 0.96) && (h < 1.04)){
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322 */
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339 aaa=(-3887.499999985135*fracl0*rhow0*m + 3906.249999985064*exp(mi*log(0.96))*rhoc*
01340 exp(mi*log((fracl0*rhow0)/rhoc))*(0.04166666666666362 + m))/(m*rhow);
01341
01342 bbb=(11663.249999955411*fracl0*rhow0*m - 11718.749999955196*exp(mi*log(0.96))*rhoc*
01343 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042222222222221294 + m))/(m*rhow);
01344
01345 ccc=(-11645.279999955485*fracl0*rhow0*m + 11699.999999955271*exp(mi*log(0.96))*rhoc*
01346 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042824074074075444 + m))/(m*rhow);
01347
01348 ddd=(3870.02879998521*fracl0*rhow0*m - 3886.9999999851375*exp(mi*log(0.96))*rhoc*
01349 exp(mi*log((fracl0*rhow0)/rhoc))*(0.043478260869569095 + m))/(m*rhow);
01350
01351 fracl = aaa*(pv/psat)*(pv/psat)*(pv/psat) + bbb*(pv/psat)*(pv/psat) + ccc*(pv/psat) + ddd;
01352
01353
01354 }
01355 if(h >= 1.04){
01356 fracl=(fracl0*rhow0*(1.0+(0.12*(h-1.04))))/rhow;
01357 }
01358
01359
01360 if(fracl > por)
01361 fracl=por;
01362
01363 return(fracl);
01364 }
01365
01366
01367
01368
01369
01370 double glasgowmatc::f_dfracldt (double rhov, double t)
01371 {
01372 double dfracldt;
01373 double fracl0,rhow0,ff,ffc,h;
01374 double aaa,bbb,ccc,daaadt,dbbbdt,dcccdt,dddddt;
01375 double a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3;
01376 double rhow,m,mi,pv,psat;
01377 double drhowdt,dmdt,dmidt,dpsatdt;
01378 double derf,derff,help;
01379
01380 h = f_h(rhov,t);
01381 fracl0 = f_fracl0(t);
01382 rhow0 = f_rhow0(t);
01383 rhow = f_rhow(t);
01384 m = f_m(t);
01385 mi = f_mi(t);
01386 pv = f_pv(rhov,t);
01387 psat = f_psat(t);
01388
01389 drhowdt = f_drhowdt(t);
01390 dmdt = f_dmdt(t);
01391 dmidt = f_dmidt(t);
01392 dpsatdt = f_dpsatdt(t);
01393
01394 a1 = -3887.499999985135;
01395 a2 = 3906.249999985064;
01396 a3 = 0.04166666666666362;
01397
01398 b1 = 11663.249999955411;
01399 b2 = -11718.749999955196;
01400 b3 = 0.042222222222221294;
01401
01402 c1 = -11645.279999955485;
01403 c2 = 11699.999999955271;
01404 c3 = 0.042824074074075444;
01405
01406 d1 = 3870.02879998521;
01407 d2 = -3886.9999999851375;
01408 d3 = 0.043478260869569095;
01409
01410 ff = fracl0*rhow0;
01411 ffc= fracl0*rhow0/rhoc;
01412
01413
01414
01415 if(h <= 0.96){
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431 derf = ffc*rhov*rv/psat - ffc*pv/psat/psat*dpsatdt;
01432 derff = (exp(mi*log(ffc*pv/psat)))*(dmidt*log(ffc*pv/psat) +
01433 mi*derf/(ffc*pv/psat));
01434 dfracldt = derff*rhoc/rhow - exp(mi*log(ffc*pv/psat))*rhoc/rhow/rhow*drhowdt;
01435 }
01436 if((h > 0.96) && (h < 1.04)){
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450 */
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464 aaa=(-3887.499999985135*fracl0*rhow0*m + 3906.249999985064*exp(mi*log(0.96))*rhoc*
01465 exp(mi*log((fracl0*rhow0)/rhoc))*(0.04166666666666362 + m))/(m*rhow);
01466
01467 bbb=(11663.249999955411*fracl0*rhow0*m - 11718.749999955196*exp(mi*log(0.96))*rhoc*
01468 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042222222222221294 + m))/(m*rhow);
01469
01470 ccc=(-11645.279999955485*fracl0*rhow0*m + 11699.999999955271*exp(mi*log(0.96))*rhoc*
01471 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042824074074075444 + m))/(m*rhow);
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507 daaadt = -1.0*a1*ff/rhow/rhow*drhowdt + a2*a3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
01508 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
01509 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
01510 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
01511 a2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
01512 *dmidt/rhow +
01513 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
01514
01515
01516 dbbbdt = -1.0*b1*ff/rhow/rhow*drhowdt + b2*b3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
01517 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
01518 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
01519 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
01520 b2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
01521 *dmidt/rhow +
01522 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
01523
01524
01525 dcccdt = -1.0*c1*ff/rhow/rhow*drhowdt + c2*c3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
01526 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
01527 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
01528 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
01529 c2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
01530 *dmidt/rhow +
01531 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
01532
01533
01534 dddddt = -1.0*d1*ff/rhow/rhow*drhowdt + d2*d3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
01535 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
01536 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
01537 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
01538 d2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
01539 *dmidt/rhow +
01540 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574 help = (rhov*rv*psat - pv*dpsatdt)/psat/psat;
01575 dfracldt = daaadt*(pv/psat)*(pv/psat)*(pv/psat) + dbbbdt*(pv/psat)*(pv/psat) + dcccdt*(pv/psat) + dddddt;
01576 dfracldt = dfracldt + 3.0*aaa*(pv/psat)*(pv/psat)*help + 2.0*bbb*(pv/psat)*help + ccc*help;
01577
01578
01579 }
01580 if(h >= 1.04){
01581
01582
01583
01584
01585 help = 0.12*(rhov*rv*psat - pv*dpsatdt)/psat/psat;
01586 dfracldt = fracl0*rhow0*(help*rhow - (1.0+(0.12*(h+1.04)))*drhowdt)/rhow/rhow;
01587 }
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599 return(dfracldt);
01600 }
01601
01602
01603
01604
01605
01606 double glasgowmatc::f_dfracldrhov (double rhov, double t)
01607 {
01608 double dfracldrhov;
01609 double fracl0,rhow0,h;
01610 double rhow,m,pv,psat;
01611 double aaa,bbb,ccc;
01612
01613 h = f_h(rhov,t);
01614 fracl0 = f_fracl0(t);
01615 rhow0 = f_rhow0(t);
01616 rhow = f_rhow(t);
01617 m = f_m(t);
01618 pv = f_pv(rhov,t);
01619 psat = f_psat(t);
01620
01621 if(h <= 0.96){
01622
01623
01624
01625 dfracldrhov=(fracl0*rhow0*exp((-1.0 + 1.0/m)*log((fracl0*rhow0*pv)/(rhoc*psat)))*rv*t)/(m*psat*rhow);
01626 }
01627 if((h > 0.96) && (h < 1.04)){
01628
01629
01630
01631
01632
01633
01634
01635 */
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649 aaa=(-3887.499999985135*fracl0*rhow0*m + 3906.249999985064*exp((1.0/m)*log(0.96))*rhoc*
01650 exp((1.0/m)*log((fracl0*rhow0)/rhoc))*(0.04166666666666362 + m))/(m*rhow);
01651
01652 bbb=(11663.249999955411*fracl0*rhow0*m - 11718.749999955196*exp((1.0/m)*log(0.96))*rhoc*
01653 exp((1.0/m)*log((fracl0*rhow0)/rhoc))*(0.042222222222221294 + m))/(m*rhow);
01654
01655 ccc=(-11645.279999955485*fracl0*rhow0*m + 11699.999999955271*exp((1.0/m)*log(0.96))*rhoc*
01656 exp((1.0/m)*log((fracl0*rhow0)/rhoc))*(0.042824074074075444 + m))/(m*rhow);
01657
01658
01659
01660 dfracldrhov = 3.0*aaa*(pv/psat)*(pv/psat/psat)*(rv*t) + 2.0*bbb*(pv/psat/psat)*(rv*t) + ccc*(rv*t/psat);
01661
01662
01663 }
01664 if(h >= 1.04){
01665
01666 dfracldrhov=(0.12*fracl0*rhow0*rv*t)/(psat*rhow);
01667 }
01668 return(dfracldrhov);
01669 }
01670
01671
01672
01673
01674
01675
01676 double glasgowmatc::f_mcbwrel (double t)
01677 {
01678 double tc,mcbwrel,mcbw0;
01679
01680 tc=t-273.15;
01681
01682
01683
01684
01685
01686
01687
01688 mcbw0=0.09*rhoc;
01689 if(tc <= 200.0){
01690 mcbwrel=0.0;
01691 }
01692 if((tc > 200.0) && (tc <= 300.0)){
01693 mcbwrel=rhoc*(7.0e-4*(tc-200.0));
01694 }
01695 if((tc > 300.0) && (tc <= 800.0)){
01696 mcbwrel=rhoc*(0.4e-4*(tc-300.0)+0.07);
01697 }
01698 if(tc > 800.0){
01699 mcbwrel=rhoc*0.09;
01700 }
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733 return(mcbwrel);
01734 }
01735
01736
01737
01738
01739
01740 double glasgowmatc::f_dmcbwreldt (double t)
01741 {
01742 double tc,dmcbwreldt,mcbw0;
01743
01744 tc=t-273.15;
01745
01746
01747
01748
01749
01750
01751 mcbw0=0.09*rhoc;
01752 if(tc <= 200.0){
01753 dmcbwreldt=0.0;
01754 }
01755 if((tc > 200.0) && (tc <= 300.0)){
01756 dmcbwreldt=rhoc*7.0e-4;
01757 }
01758 if((tc > 300.0) && (tc <= 800.0)){
01759 dmcbwreldt=rhoc*0.4e-4;
01760 }
01761 if(tc > 800.0){
01762 dmcbwreldt=0.0;
01763 }
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796 return(dmcbwreldt);
01797 }
01798
01799
01800
01801
01802
01803 double glasgowmatc::f_fracd (double t)
01804 {
01805 double fracd,mcbwrel,rhow;
01806
01807 mcbwrel = f_mcbwrel(t);
01808 rhow = f_rhow(t);
01809 fracd=mcbwrel/rhow;
01810
01811 return(fracd);
01812 }
01813
01814
01815
01816
01817
01818 double glasgowmatc::f_dfracddt(double t)
01819 {
01820 double dfracddt,dmcbwreldt,rhow;
01821
01822 dmcbwreldt = f_dmcbwreldt(t);
01823 rhow = f_rhow(t);
01824
01825 dfracddt=dmcbwreldt/rhow;
01826
01827 return(dfracddt);
01828 }
01829
01830
01831
01832
01833
01834 double glasgowmatc::f_fracg(double rhov, double t)
01835 {
01836 double fracg,por,fracl;
01837
01838 por = f_por(t);
01839 fracl = f_fracl(rhov,t);
01840
01841 fracg=por-fracl;
01842
01843 return(fracg);
01844 }
01845
01846
01847
01848
01849
01850 double glasgowmatc::f_s(double rhov, double t)
01851 {
01852 double s,fracl,por;
01853
01854 por = f_por(t);
01855 fracl = f_fracl(rhov,t);
01856
01857 s=fracl/por;
01858
01859 return(s);
01860 }
01861
01862
01863
01864
01865
01866 double glasgowmatc::f_pl(double pc, double pg, double t)
01867 {
01868 double pl,s,rhov;
01869
01870 s = f_s(rhov,t);
01871 pc = f_pc(rhov,t);
01872
01873 switch (model){
01874 case 1:{
01875 pl=pg;
01876 break;
01877 }
01878 case 2:{
01879 if(s <= sssp){
01880 pl = 0.0;
01881 }
01882 else{
01883 pl = pg - pc;
01884 }
01885 break;
01886 }
01887 default:{
01888 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
01889 }
01890 }
01891 return(pl);
01892 }
01893
01894
01895
01896
01897
01898 double glasgowmatc::f_pc (double rhov, double t)
01899 {
01900 double pc,s,h,rhow;
01901
01902 s = f_s(rhov,t);
01903 h = f_h(rhov,t);
01904 rhow = f_rhow(t);
01905
01906 switch (model){
01907 case 1:{
01908 pc=0.0;
01909 break;
01910 }
01911 case 2:{
01912 if(s <= sssp){
01913
01914 pc=-1.0*rv*t*rhow*log(h);
01915 }
01916 else{
01917 pc=-1.0*rv*t*rhow*log(h);
01918 if(pc < 0.0) {
01919 pc=0.0;
01920 }
01921 }
01922 break;
01923 }
01924 default:{
01925 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
01926 }
01927 }
01928 return(pc);
01929 }
01930
01931
01932
01933
01934
01935 double glasgowmatc::f_dpcdt (double rhov, double t)
01936 {
01937 double dpcdt,pc,s,h,rhow,drhowdt,psat,dpsatdt;
01938
01939 s = f_s(rhov,t);
01940 h = f_h(rhov,t);
01941 rhow = f_rhow(t);
01942 drhowdt = f_drhowdt(t);
01943 psat = f_psat(t);
01944 dpsatdt = f_dpsatdt(t);
01945
01946 switch (model){
01947 case 1:{
01948 dpcdt=0.0;
01949 break;
01950 }
01951 case 2:{
01952 if(s <= sssp){
01953 dpcdt=0.0;
01954 }
01955 else{
01956 pc=-1.0*rv*t*rhow*log(h);
01957 if(pc < 0.0) {
01958 dpcdt=0.0;
01959 }
01960 else{
01961 dpcdt=rv*rhow*((t*dpsatdt/psat)-((1.0+t*drhowdt/rhow)*log(rhov*rv*t/psat))-1.0);
01962 }
01963 }
01964 break;
01965 }
01966 default:{
01967 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
01968 }
01969 }
01970
01971 return(dpcdt);
01972 }
01973
01974
01975
01976
01977
01978 double glasgowmatc::f_dpcdrhov (double rhov, double t)
01979 {
01980 double dpcdrhov,pc,s,h,rhow;
01981
01982 s = f_s(rhov,t);
01983 h = f_h(rhov,t);
01984 rhow = f_rhow(t);
01985
01986 switch (model){
01987 case 1:{
01988 dpcdrhov=0.0;
01989 break;
01990 }
01991 case 2:{
01992 if(s <= sssp){
01993 dpcdrhov=0.0;
01994 }
01995 else{
01996 pc=-1.0*rv*t*rhow*log(h);
01997 if(pc < 0.0) {
01998 dpcdrhov=0.0;
01999 }
02000 else{
02001 dpcdrhov=-rv*t*rhow/rhov;
02002 }
02003 }
02004 break;
02005 }
02006 default:{
02007 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02008 }
02009 }
02010
02011 return(dpcdrhov);
02012 }
02013
02014
02015
02016
02017
02018 double glasgowmatc::f_sln (double rhov, double t)
02019 {
02020 double sln,s;
02021
02022 s = f_s(rhov,t);
02023
02024 switch (model){
02025 case 1:{
02026 sln=0.0;
02027 break;
02028 }
02029 case 2:{
02030 if(s <= sssp){
02031 sln=0.0;
02032 }
02033 else{
02034 sln=(s-sssp)/(1.0-sssp);
02035 }
02036 break;
02037 }
02038 default:{
02039 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02040 }
02041 }
02042
02043 return(sln);
02044 }
02045
02046
02047
02048
02049
02050 double glasgowmatc::f_ppore (double rhov, double pg, double t, double pginf)
02051 {
02052 double ppore;
02053 double sln,pc,s,h,rhow;
02054
02055 s = f_s(rhov,t);
02056 h = f_h(rhov,t);
02057 rhow = f_rhow(t);
02058
02059 switch (model){
02060 case 1:{
02061 ppore=pg-pginf;
02062 pc=0.0;
02063 sln=0.0;
02064 break;
02065 }
02066 case 2:{
02067 if(s <= sssp){
02068 pc=0.0;
02069 sln=0.0;
02070 ppore=pg-pginf;
02071 }
02072 else{
02073 pc=-1.0*rv*t*rhow*log(h);
02074 if(pc < 0.0) {
02075 pc=0.0;
02076 }
02077 else{}
02078 sln=(s-sssp)/(1.0-sssp);
02079 ppore=pg-(sln*pc)-pginf;
02080 }
02081 break;
02082 }
02083 default:{
02084 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02085 }
02086 }
02087
02088 return(ppore);
02089 }
02090
02091
02092
02093
02094
02095 double glasgowmatc::f_sb (double rhov, double t)
02096 {
02097 double sb,s;
02098
02099 s = f_s(rhov,t);
02100
02101 switch (model){
02102 case 1:{
02103 sb=0.0;
02104 break;
02105 }
02106 case 2:{
02107 if(s <= sssp){
02108 sb=s;
02109 }
02110 else{
02111 sb=sssp;
02112 }
02113 break;
02114 }
02115 default:{
02116 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02117 }
02118 }
02119 return(sb);
02120 }
02121
02122
02123
02124
02125
02126 double glasgowmatc::f_pa (double rhov, double pg, double t)
02127 {
02128 double pa,pv;
02129
02130 pv = f_pv(rhov,t);
02131
02132 pa=pg-pv;
02133
02134 return(pa);
02135 }
02136
02137
02138
02139
02140
02141 double glasgowmatc::f_rhoa (double rhov, double pg, double t)
02142 {
02143 double rhoa,pa;
02144
02145 pa = f_pa(rhov,pg,t);
02146
02147 rhoa=pa/ra/t;
02148
02149 return(rhoa);
02150 }
02151
02152
02153
02154
02155
02156 double glasgowmatc::f_rhog (double rhov, double pg, double t)
02157 {
02158 double rhog,rhoa;
02159
02160 rhoa = f_rhoa(rhov,pg,t);
02161
02162 rhog=rhoa+rhov;
02163
02164
02165
02166 return(rhog);
02167 }
02168
02169
02170
02171
02172
02173 double glasgowmatc::f_dav (double rhov, double pg, double t)
02174 {
02175 double dif,delta,tor,dav;
02176
02177
02178
02179
02180
02181
02182
02183 dif=1.87*exp(2.072*log(t))/pg*1.0e-5;
02184
02185 delta=0.5;
02186 tor=3.0;
02187 dav=dif*delta/tor/tor;
02188
02189
02190
02191
02192
02193 return(dav);
02194 }
02195
02196
02197
02198
02199
02200 double glasgowmatc::f_davex (double pg, double tinf)
02201 {
02202 double davex,tor,delta,difex;
02203
02204
02205
02206
02207
02208
02209
02210 delta=0.5;
02211 tor=3.0;
02212
02213
02214
02215 difex=1.87*exp(2.072*log(tinf))/pg*1.0e-5;
02216
02217 davex=difex*delta/tor/tor;
02218
02219 return(davex);
02220 }
02221
02222
02223
02224
02225
02226 double glasgowmatc::f_db(double rhov, double t)
02227 {
02228 double db,db0,dbtref,s;
02229
02230 s = f_s(rhov,t);
02231
02232
02233 db0=1.57e-11;
02234 dbtref=295.0;
02235 if(s <= sssp)
02236 db=db0*exp(-2.08*(s/sssp)*(t/dbtref));
02237 else
02238 db=0.0;
02239
02240 return(db);
02241 }
02242
02243
02244
02245
02246
02247
02248 double glasgowmatc::f_kk (double t)
02249 {
02250 double kk,por;
02251
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265 por = f_por(t);
02266
02267
02268 kk=exp((2.0/3.0)*log(por/por0))*k0;
02269
02270
02271 kk = k0;
02272
02273
02274 return(kk);
02275 }
02276
02277
02278
02279
02280
02281 double glasgowmatc::f_kg (double rhov, double t)
02282 {
02283 double kg,kb,km,s;
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305 s = f_s(rhov,t);
02306
02307
02308 kb=2.2748;
02309 km=1.0/kb;
02310
02311 kg=(sqrt(1.0-s))*exp((2.0*km)*log(1.0-exp((1.0/km)*log(s))));
02312
02313
02314
02315
02316
02317
02318 return(kg);
02319 }
02320
02321
02322
02323
02324
02325 double glasgowmatc::f_kl (double rhov, double t)
02326 {
02327 double kl,km,kb,s;
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349 s = f_s(rhov,t);
02350
02351
02352 kb=2.2748;
02353 km=1.0/kb;
02354
02355 kl=sqrt(s)*exp(2.0*log(1.0-exp(km*log(1.0-exp((1.0/km)*log(s))))));
02356
02357
02358
02359
02360
02361
02362 return(kl);
02363 }
02364
02365
02366
02367
02368
02369 double glasgowmatc::f_mul (double t)
02370 {
02371 double mul;
02372
02373
02374
02375
02376
02377 if(t <= 647.3)
02378
02379
02380 mul=0.6612*exp((-1.562)*log(t-229.0));
02381 else
02382 mul=0.000043;
02383
02384 return(mul);
02385 }
02386
02387
02388
02389
02390
02391
02392 double glasgowmatc::f_muv (double t)
02393 {
02394 double muv,t0muv,amuv,muv0;
02395
02396
02397
02398
02399
02400 if(t <= 647.3){
02401 t0muv=273.15;
02402 muv0=8.85e-6;
02403 amuv=3.53e-8;
02404 muv=muv0+amuv*(t-t0muv);
02405 }
02406 else
02407 muv=0.00004313;
02408
02409 return(muv);
02410 }
02411
02412
02413
02414
02415
02416
02417 double glasgowmatc::f_mua (double t)
02418 {
02419 double mua,t0mua,amua,bmua,mua0;
02420
02421
02422
02423
02424
02425 t0mua=273.15;
02426 mua0=17.17e-6;
02427 amua=4.73e-8;
02428 bmua=2.22e-11;
02429
02430
02431 mua=mua0+(amua*(t-t0mua))-(bmua*exp(2.0*log(t-t0mua)));
02432
02433 return(mua);
02434 }
02435
02436
02437
02438
02439
02440
02441 double glasgowmatc::f_mug (double rhov, double pg, double t)
02442 {
02443 double mug,mua,muv,rhoa;
02444
02445 mua = f_mua(t);
02446 muv = f_muv(t);
02447 rhoa = f_rhoa(rhov,pg,t);
02448
02449
02450 if((rhoa == 0.0) && (rhov == 0.0))
02451 mug=0.0;
02452 else
02453 mug=((rhoa*mua)+(rhov*muv))/(rhoa+rhov);
02454
02455
02456
02457
02458 return(mug);
02459 }
02460
02461
02462
02463
02464
02465
02466 double glasgowmatc::f_cl (double t)
02467 {
02468 double cl;
02469
02470
02471
02472
02473
02474 if(t <= 647.3)
02475
02476
02477
02478 cl=((2.4768*t)+3368.2)+exp(31.4447657616636*log((1.08542631988638*t)/513.15));
02479 else
02480 cl=24515.0;
02481
02482 return(cl);
02483 }
02484
02485
02486
02487
02488
02489
02490 double glasgowmatc::f_cv (double t)
02491 {
02492 double cv;
02493
02494
02495
02496
02497
02498 if(t <= 647.3)
02499
02500
02501
02502 cv=((7.1399*t)-443.0)+exp(29.4435287521143*log((1.13771502228162*t)/513.15));
02503 else
02504 cv=45821.04;
02505
02506 return(cv);
02507 }
02508
02509
02510
02511
02512
02513
02514 double glasgowmatc::f_ca (double t)
02515 {
02516 double ca;
02517
02518
02519
02520
02521 ca = -0.00000009849367018147350*t*t*t + 0.0003564362577698610*t*t - 0.1216179239877570*t +1012.502552163240;
02522
02523 return(ca);
02524 }
02525
02526
02527
02528
02529
02530
02531 double glasgowmatc::f_cs (double t)
02532 {
02533 double cs,tc;
02534
02535 tc=t-273.15;
02536
02537
02538
02539
02540
02541 cs=900.0+80.0*(tc/120.0)-(4.0*tc*tc/120.0/120.0);
02542
02543
02544
02545 return(cs);
02546 }
02547
02548
02549
02550
02551
02552
02553 double glasgowmatc::f_keff (double t)
02554 {
02555 double keff,tc;
02556
02557 tc=t-273.15;
02558
02559
02560
02561
02562
02563 keff=2.0-(0.24*(tc/120.0))+(0.012*tc*tc/120.0/120.0);
02564
02565
02566
02567
02568 return(keff);
02569 }
02570
02571
02572
02573
02574
02575
02576 double glasgowmatc::f_fracs (double rhov, double t)
02577 {
02578 double fracs,fracg,fracl;
02579
02580 fracg = f_fracg(rhov,t);
02581 fracl = f_fracl(rhov,t);
02582
02583 fracs=1.0-fracg-fracl;
02584
02585 return(fracs);
02586 }
02587
02588
02589
02590
02591
02592
02593 double glasgowmatc::f_crho (double rhov, double pg, double t)
02594 {
02595 double crho,fracs,cs,fracl,rhow,cl,cv,fracg,rhoa,ca;
02596
02597 fracs = f_fracs(rhov,t);
02598 cs = f_cs(t);
02599 fracl = f_fracl(rhov,t);
02600 rhow = f_rhow(t);
02601 cl = f_cl(t);
02602 fracg = f_fracg(rhov,t);
02603 cv = f_cv(t);
02604 rhoa = f_rhoa(rhov,pg,t);
02605 ca = f_ca(t);
02606
02607 crho=(rhos*fracs*cs)+(fracl*rhow*cl)+(fracg*rhov*cv)+(fracg*rhoa*ca);
02608
02609 return(crho);
02610 }
02611
02612
02613
02614
02615
02616
02617 double glasgowmatc::f_le (double t)
02618 {
02619 double le;
02620
02621
02622
02623
02624
02625
02626
02627
02628
02629 if(t <= 647.3)
02630
02631
02632 le=2.672e5*exp(0.38*log(647.3-t));
02633 else
02634 le=0.0;
02635
02636 return(le);
02637 }
02638
02639
02640
02641
02642
02643
02644 double glasgowmatc::f_ld ()
02645 {
02646 double ld;
02647
02648
02649 ld=2400.0e3;
02650
02651 return(ld);
02652 }
02653
02654
02655
02656
02657
02658
02659 double glasgowmatc::f_hrad(double t, double tinf)
02660 {
02661 double hrad;
02662
02663
02664
02665
02666
02667 hrad=emmi*stef*((t*t)+(tinf*tinf))*(t+tinf);
02668
02669 return(hrad);
02670 }
02671
02672
02673
02674
02675
02676 double glasgowmatc::f_hqr(double t, double tinf)
02677 {
02678 double hrad,hqr;
02679
02680 hrad=f_hrad(t,tinf);
02681
02682
02683
02684
02685 hqr=hq+hrad;
02686
02687 return(hqr);
02688 }
02689
02690
02691
02692
02693
02694
02695 double glasgowmatc::f_beta (double pg, double tinf)
02696 {
02697 double beta,davex;
02698
02699 davex = f_davex(pg,tinf);
02700
02701
02702
02703
02704 beta=(hq/crhoair)*exp((2.0/3.0)*log(davex/alph));
02705
02706 return(beta);
02707 }
02708
02709
02710
02711
02712
02713
02714 double glasgowmatc::f_alpha (double t)
02715 {
02716
02717
02718 double tt,alpha;
02719
02720 tt = t - 273.15;
02721 tt = (tt - 20.0)/100.0;
02722 if(tt >= 0.0 && tt <= 6.0)
02723 alpha = 6.0e-5/(7.0 - tt);
02724 else
02725 alpha = 0.0;
02726
02727 alpha = 10.0e-6;
02728
02729 return(alpha);
02730 }
02731
02732
02733
02734
02735
02736 double glasgowmatc::ktt (double t,double pg,double rhov)
02737 {
02738 double ktt,keff;
02739
02740 keff = f_keff(t);
02741
02742 switch (model){
02743 case 1:{
02744 ktt = keff;
02745 break;
02746 }
02747 case 2:{
02748 double klt,fracl,rhow,sb,db,s,por,dfracldt,dpordt,kk,kl,dpcdt,mul,le;
02749
02750 fracl = f_fracl(rhov,t);
02751 rhow = f_rhow(t);
02752 sb = f_sb(rhov,t);
02753 db = f_db(rhov,t);
02754 s = f_s(rhov,t);
02755 por = f_por(t);
02756 dfracldt = f_dfracldt(rhov,t);
02757 dpordt = f_dpordt(t);
02758 kk = f_kk(t);
02759 kl = f_kl(rhov,t);
02760 dpcdt = f_dpcdt(rhov,t);
02761 mul = f_mul(t);
02762 le = f_le(t);
02763
02764 klt=fracl*rhow*(((sb*db/s/por)*(dfracldt-(fracl*dpordt/por)))-(1.0-(sb/s))*(kk*kl*dpcdt/mul));
02765 ktt = keff-(le*klt);
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775 break;
02776 }
02777 default:{
02778 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02779 }
02780 }
02781
02782 return (ktt);
02783 }
02784
02785 double glasgowmatc::ktp (double t,double pg,double rhov)
02786 {
02787 double ktp,klp;
02788 double sb,s,fracl,rhow,kk,kl,mul,le;
02789
02790 sb = f_sb(rhov,t);
02791 s = f_s(rhov,t);
02792 fracl = f_fracl(rhov,t);
02793 rhow = f_rhow(t);
02794 kk = f_kk(t);
02795 kl = f_kl(rhov,t);
02796 mul = f_mul(t);
02797 le = f_le(t);
02798
02799 switch (model){
02800 case 1:{
02801 ktp= -le*fracl*rhow*kk*kl/mul;
02802 break;
02803 }
02804 case 2:{
02805 klp=(1.0-(sb/s))*fracl*rhow*kk*kl/mul;
02806 ktp= -le*klp;
02807 break;
02808 }
02809 default:{
02810 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02811 }
02812 }
02813
02814 return (ktp);
02815 }
02816
02817 double glasgowmatc::ktv (double t,double pg,double rhov)
02818 {
02819 double ktv;
02820 double klv,db,sb,s,fracl,rhow,kk,kl,mul,le,dfracldrhov,dpcdrhov,por;
02821
02822 sb = f_sb(rhov,t);
02823 db = f_db(rhov,t);
02824 s = f_s(rhov,t);
02825 fracl = f_fracl(rhov,t);
02826 dfracldrhov = f_dfracldrhov(rhov,t);
02827 dpcdrhov = f_dpcdrhov(rhov,t);
02828 rhow = f_rhow(t);
02829 kk = f_kk(t);
02830 por = f_por(t);
02831 kl = f_kl(rhov,t);
02832 mul = f_mul(t);
02833 le = f_le(t);
02834
02835 switch (model){
02836 case 1:{
02837 ktv=0.0;
02838 break;
02839 }
02840 case 2:{
02841 klv=fracl*rhow*((sb*db*dfracldrhov/s/por)-((1.0-(sb/s))*kk*kl*dpcdrhov/mul));
02842 ktv= -le*klv;
02843 break;
02844 }
02845 default:{
02846 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02847 }
02848 }
02849
02850 return (ktv);
02851 }
02852
02853 double glasgowmatc::kat (double t,double pg,double rhov)
02854 {
02855 double kat,dav,fracg,rhog;
02856
02857 dav = f_dav(rhov,pg,t);
02858 fracg = f_fracg(rhov,t);
02859 rhog = f_rhog(rhov,pg,t);
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869 kat= -1.0*dav*fracg*rhov*pg/rhog/ra/t/t;
02870
02871 return (kat);
02872 }
02873
02874 double glasgowmatc::kap (double t,double pg,double rhov)
02875 {
02876 double kap;
02877 double dav,fracg,rhoa,rhog,kk,kg,mug;
02878
02879 dav = f_dav(rhov,pg,t);
02880 fracg = f_fracg(rhov,t);
02881 rhoa = f_rhoa(rhov,pg,t);
02882 rhog = f_rhog(rhov,pg,t);
02883 kk = f_kk(t);
02884 kg = f_kg(rhov,t);
02885 mug = f_mug(rhov,pg,t);
02886
02887 kap= (kk*kg*fracg*rhoa/mug)+(dav*fracg*rhov/rhog/ra/t);
02888
02889 return (kap);
02890 }
02891
02892 double glasgowmatc::kav (double t,double pg,double rhov)
02893 {
02894 double kav;
02895 double dav,fracg,rhoa,rhog;
02896
02897 dav = f_dav(rhov,pg,t);
02898 fracg = f_fracg(rhov,t);
02899 rhoa = f_rhoa(rhov,pg,t);
02900 rhog = f_rhog(rhov,pg,t);
02901
02902 kav= -1.0*(rhoa+(rv*rhov/ra))*dav*fracg/rhog;
02903
02904 return (kav);
02905 }
02906
02907 double glasgowmatc::kmt (double t,double pg,double rhov)
02908 {
02909 double kmt,klt,kvt;
02910 double dav,fracg,rhoa,rhog;
02911
02912 dav = f_dav(rhov,pg,t);
02913 fracg = f_fracg(rhov,t);
02914 rhoa = f_rhoa(rhov,pg,t);
02915 rhog = f_rhog(rhov,pg,t);
02916
02917 switch (model){
02918 case 1:{
02919 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
02920 klt=0.0;
02921 kmt=kvt+klt;
02922 break;
02923 }
02924 case 2:{
02925 double db,sb,s,fracl,rhow,kk,kl,mul,le,dfracldt,dpcdt,por,dpordt;
02926
02927 sb = f_sb(rhov,t);
02928 db = f_db(rhov,t);
02929 s = f_s(rhov,t);
02930 fracl = f_fracl(rhov,t);
02931 dfracldt = f_dfracldt(rhov,t);
02932 dpcdt = f_dpcdt(rhov,t);
02933 rhow = f_rhow(t);
02934 kk = f_kk(t);
02935 por = f_por(t);
02936 dpordt = f_dpordt(t);
02937 kl = f_kl(rhov,t);
02938 mul = f_mul(t);
02939 le = f_le(t);
02940
02941 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
02942 klt=fracl*rhow*(((sb*db/s/por)*(dfracldt-(fracl*dpordt/por)))-(1.0-(sb/s))*(kk*kl*dpcdt/mul));
02943 kmt=kvt+klt;
02944 break;
02945 }
02946 default:{
02947 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02948 }
02949 }
02950
02951 return (kmt);
02952 }
02953
02954 double glasgowmatc::kmp (double t,double pg,double rhov)
02955 {
02956 double kmp,kvp,klp;
02957 double dav,fracg,rhoa,rhog;
02958 double sb,s,fracl,rhow,kk,kl,mul,le;
02959 double kg,mug;
02960
02961 dav = f_dav(rhov,pg,t);
02962 fracg = f_fracg(rhov,t);
02963 rhoa = f_rhoa(rhov,pg,t);
02964 rhog = f_rhog(rhov,pg,t);
02965 sb = f_sb(rhov,t);
02966 s = f_s(rhov,t);
02967 fracl = f_fracl(rhov,t);
02968 rhow = f_rhow(t);
02969 kk = f_kk(t);
02970 kl = f_kl(rhov,t);
02971 kg = f_kg(rhov,t);
02972 mul = f_mul(t);
02973 mug = f_mug(rhov,pg,t);
02974 le = f_le(t);
02975
02976 switch (model){
02977 case 1:{
02978 kvp=fracg*rhov*((kk*kg/mug)-(dav/rhog/ra/t));
02979 klp=fracl*rhow*(kk*kl/mul);
02980 kmp=kvp+klp;
02981 break;
02982 }
02983 case 2:{
02984 kvp=fracg*rhov*((kk*kg/mug)-(dav/rhog/ra/t));
02985 klp=(1.0-(sb/s))*fracl*rhow*kk*kl/mul;
02986 kmp=kvp+klp;
02987 break;
02988 }
02989 default:{
02990 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
02991 }
02992 }
02993 return(kmp);
02994 }
02995
02996 double glasgowmatc::kmv (double t,double pg,double rhov)
02997 {
02998 double kmv,kvv,klv;
02999 double dav,fracg,fracl,rhoa,rhog,sb,db,s,dfracldrhov,rhow,kk,kl,mul,por,dpcdrhov;
03000
03001 dav = f_dav(rhov,pg,t);
03002 fracg = f_fracg(rhov,t);
03003 rhoa = f_rhoa(rhov,pg,t);
03004 rhog = f_rhog(rhov,pg,t);
03005 sb = f_sb(rhov,t);
03006 db = f_db(rhov,t);
03007 s = f_s(rhov,t);
03008 fracl = f_fracl(rhov,t);
03009 dfracldrhov = f_dfracldrhov(rhov,t);
03010 rhow = f_rhow(t);
03011 kk = f_kk(t);
03012 kl = f_kl(rhov,t);
03013 mul = f_mul(t);
03014 por = f_por(t);
03015 dpcdrhov = f_dpcdrhov(rhov,t);
03016
03017 switch (model){
03018 case 1:{
03019 kvv=(fracg*dav/rhog)*(rhoa+(rv*rhov/ra));
03020 klv=0.0;
03021 kmv=kvv+klv;
03022 break;
03023 }
03024 case 2:{
03025 kvv=(fracg*dav/rhog)*(rhoa+(rv*rhov/ra));
03026 klv=fracl*rhow*((sb*db*dfracldrhov/s/por)-((1.0-(sb/s))*kk*kl*dpcdrhov/mul));
03027 kmv=kvv+klv;
03028 break;
03029 }
03030 default:{
03031 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
03032 }
03033 }
03034 return (kmv);
03035 }
03036
03037 double glasgowmatc::ctt (double t,double pg,double rhov)
03038 {
03039 double ctt;
03040 double crho,le,ld,rhow,drhowdt,fracl,fracd,dfracldt,dfracddt;
03041
03042 crho = f_crho(rhov,pg,t);
03043 rhow = f_rhow(t);
03044 drhowdt = f_drhowdt(t);
03045 fracl = f_fracl(rhov,t);
03046 fracd = f_fracd(t);
03047 dfracldt = f_dfracldt(rhov,t);
03048 dfracddt = f_dfracddt(t);
03049 ld = f_ld();
03050 le = f_le(t);
03051
03052 ctt = (crho+((ld+le)*(fracd*drhowdt+rhow*dfracddt))-(le*(fracl*drhowdt+rhow*dfracldt)));
03053
03054 return(ctt);
03055 }
03056
03057 double glasgowmatc::ctp ()
03058 {
03059 return 0.0;
03060 }
03061
03062 double glasgowmatc::ctv (double t,double pg,double rhov)
03063 {
03064 double ctv;
03065 double le,rhow,dfracldrhov;
03066
03067 le = f_le(t);
03068 rhow = f_rhow(t);
03069 dfracldrhov = f_dfracldrhov(rhov,t);
03070
03071 ctv = (-le*rhow*dfracldrhov);
03072
03073 return(ctv);
03074 }
03075
03076 double glasgowmatc::cat (double t,double pg,double rhov)
03077 {
03078 double cat;
03079 double dpordt,rhoa,fracg,dfracldt;
03080
03081 dpordt = f_dpordt(t);
03082 rhoa = f_rhoa(rhov,pg,t);
03083 fracg = f_fracg(rhov,t);
03084 dfracldt = f_dfracldt(rhov,t);
03085
03086 cat = ((rhoa*(dpordt-dfracldt))-(fracg*pg/ra/t/t));
03087
03088 return(cat);
03089 }
03090
03091 double glasgowmatc::cap (double t,double pg,double rhov)
03092 {
03093 double cap;
03094 double fracg;
03095
03096 fracg = f_fracg(rhov,t);
03097
03098 cap = (fracg/ra/t);
03099
03100 return(cap);
03101 }
03102
03103 double glasgowmatc::cav (double t,double pg,double rhov)
03104 {
03105 double cav;
03106 double fracg,rhoa,dfracldrhov;
03107
03108 fracg = f_fracg(rhov,t);
03109 rhoa = f_rhoa(rhov,pg,t);
03110 dfracldrhov = f_dfracldrhov(rhov,t);
03111
03112 cav = (-1.0*(fracg*rv/ra)-(rhoa*dfracldrhov));
03113
03114 return(cav);
03115 }
03116
03117 double glasgowmatc::cmt (double t,double pg,double rhov)
03118 {
03119 double cmt;
03120 double dpordt,rhow,drhowdt,fracl,fracd,dfracldt,dfracddt;
03121
03122 dpordt = f_dpordt(t);
03123 rhow = f_rhow(t);
03124 drhowdt = f_drhowdt(t);
03125 fracl = f_fracl(rhov,t);
03126 fracd = f_fracd(t);
03127 dfracldt = f_dfracldt(rhov,t);
03128 dfracddt = f_dfracddt(t);
03129
03130 cmt = ((rhov*dpordt)+((fracl-fracd)*drhowdt)+((rhow-rhov)*dfracldt)-(rhow*dfracddt));
03131
03132 return(cmt);
03133 }
03134
03135 double glasgowmatc::cmp ()
03136 {
03137 return 0.0;
03138 }
03139
03140 double glasgowmatc::cmv (double t,double pg,double rhov)
03141 {
03142 double cmv;
03143 double fracg,rhow,dfracldrhov;
03144
03145 fracg = f_fracg(rhov,t);
03146 rhow = f_rhow(t);
03147 dfracldrhov = f_dfracldrhov(rhov,t);
03148
03149 cmv = (fracg+((rhow-rhov)*dfracldrhov));
03150
03151 return(cmv);
03152 }
03153
03154
03155
03156
03157
03158
03159
03160
03161
03162
03163
03164
03165
03166
03167
03168 double glasgowmatc::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
03169 {
03170 long k;
03171 double new_trc,rhov,pg,t;
03172 new_trc = 0.0;
03173
03174 k=Gtt->give_dof(nn,0);
03175 if (k>0) {rhov = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03176 if (k==0) {rhov = 0.0;}
03177 if (k<0) {rhov = Tb->lc[0].pv[0-k-1].getval();}
03178 k=Gtt->give_dof(nn,1);
03179 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03180 if (k==0) {pg = 0.0;}
03181 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
03182 k=Gtt->give_dof(nn,2);
03183 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03184 if (k==0) {t = 0.0;}
03185 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
03186
03187 if((ri == 0) && (ci == 0))
03188 new_trc = get_transmission_transcoeff_mv(t,pg,rhov,bc,ipp);
03189 if((ri == 0) && (ci == 1))
03190 new_trc = 0.0;
03191 if((ri == 0) && (ci == 2))
03192 new_trc = get_transmission_transcoeff_mt(t,pg,rhov,bc,ipp);
03193
03194 if((ri == 1) && (ci == 0))
03195 new_trc = 0.0;
03196 if((ri == 1) && (ci == 1))
03197 new_trc = 0.0;
03198 if((ri == 1) && (ci == 2))
03199 new_trc = 0.0;
03200
03201 if((ri == 2) && (ci == 0))
03202 new_trc = 0.0;
03203 if((ri == 2) && (ci == 1))
03204 new_trc = 0.0;
03205 if((ri == 2) && (ci == 2))
03206 new_trc = get_transmission_transcoeff_tt(t,pg,rhov,bc,ipp);
03207
03208 new_trc = new_trc*trc;
03209
03210 return (new_trc);
03211 }
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226 double glasgowmatc::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
03227 {
03228 long k;
03229 double new_nodval,rhov,pg,t;
03230 new_nodval = 0.0;
03231
03232 k=Gtt->give_dof(nn,0);
03233 if (k>0) {rhov = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03234 if (k==0) {rhov = 0.0;}
03235 if (k<0) {rhov = Tb->lc[0].pv[0-k-1].getval();}
03236 k=Gtt->give_dof(nn,1);
03237 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03238 if (k==0) {pg = 0.0;}
03239 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
03240 k=Gtt->give_dof(nn,2);
03241 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03242 if (k==0) {t = 0.0;}
03243 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
03244
03245 if((ri == 0) && (ci == 0))
03246 new_nodval = get_transmission_nodval_mv(nodval,t,pg,rhov,bc,ipp);
03247 if((ri == 0) && (ci == 1))
03248 new_nodval = 0.0;
03249 if((ri == 0) && (ci == 2))
03250 new_nodval = get_transmission_nodval_mt(nodval,t,pg,rhov,bc,ipp);
03251
03252 if((ri == 1) && (ci == 0))
03253 new_nodval = 0.0;
03254 if((ri == 1) && (ci == 1))
03255 new_nodval = 0.0;
03256 if((ri == 1) && (ci == 2))
03257 new_nodval = 0.0;
03258
03259 if((ri == 2) && (ci == 0))
03260 new_nodval = 0.0;
03261 if((ri == 2) && (ci == 1))
03262 new_nodval = 0.0;
03263 if((ri == 2) && (ci == 2))
03264 new_nodval = get_transmission_nodval_tt(nodval,t,pg,rhov,bc,ipp);
03265
03266 return (new_nodval);
03267 }
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282
03283 double glasgowmatc::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
03284 {
03285 long k;
03286 double flux,rhov,pg,t;
03287 flux = 0.0;
03288
03289 k=Gtt->give_dof(nn,0);
03290 if (k>0) {rhov = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03291 if (k==0) {rhov = 0.0;}
03292 if (k<0) {rhov = Tb->lc[0].pv[0-k-1].getval();}
03293 k=Gtt->give_dof(nn,1);
03294 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03295 if (k==0) {pg = 0.0;}
03296 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
03297 k=Gtt->give_dof(nn,2);
03298 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
03299 if (k==0) {t = 0.0;}
03300 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
03301
03302 if((ri == 0) && (ci == 0))
03303 flux = get_transmission_flux_mv(nodval,t,pg,rhov,bc,ipp);
03304 if((ri == 0) && (ci == 1))
03305 flux = 0.0;
03306 if((ri == 0) && (ci == 2))
03307 flux = get_transmission_flux_mt(nodval,t,pg,rhov,bc,ipp);
03308
03309 if((ri == 1) && (ci == 0))
03310 flux = 0.0;
03311 if((ri == 1) && (ci == 1))
03312 flux = 0.0;
03313 if((ri == 1) && (ci == 2))
03314 flux = 0.0;
03315
03316 if((ri == 2) && (ci == 0))
03317 flux = 0.0;
03318 if((ri == 2) && (ci == 1))
03319 flux = 0.0;
03320 if((ri == 2) && (ci == 2))
03321 flux = get_transmission_flux_tt(nodval,t,pg,rhov,bc,ipp);
03322
03323 return (flux);
03324 }
03325
03326
03327
03328
03329
03330
03331
03332
03333
03334
03335
03336 double glasgowmatc::get_transmission_transcoeff_mv(double t,double pg,double rhov,long bc,long ipp)
03337 {
03338 double trc;
03339
03340 switch (bc){
03341 case 30:{
03342 trc=1.0;
03343 break;
03344 }
03345 case 31:{
03346 double beta;
03347
03348 beta = f_beta(pg,t);
03349
03350 trc = beta;
03351 break;
03352 }
03353 case 32:{
03354 double beta,tinf;
03355
03356 tinf = f_tinf(Tp->time);
03357 beta = f_beta(pg,tinf);
03358
03359 trc=beta;
03360 break;
03361 }
03362 case 36:{
03363 trc=0.0;
03364 break;
03365 }
03366 default:{
03367 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03368 exit(0);
03369 }
03370 }
03371
03372 return(trc);
03373 }
03374
03375
03376
03377
03378
03379
03380
03381
03382
03383
03384
03385
03386 double glasgowmatc::get_transmission_nodval_mv(double bv,double t,double pg,double rhov,long bc,long ipp)
03387 {
03388 double nodval;
03389
03390 switch (bc){
03391 case 30:{
03392 nodval = bv;
03393 break;
03394 }
03395 case 31:{
03396 double beta;
03397
03398 beta = f_beta(pg,t);
03399
03400 nodval = beta*bv;
03401 break;
03402 }
03403 case 32:{
03404 double beta,tinf;
03405
03406 tinf = f_tinf(Tp->time);
03407 beta = f_beta(pg,tinf);
03408
03409 nodval = beta*bv;
03410 break;
03411 }
03412 case 36:{
03413 double psat;
03414
03415 psat = f_psat(t);
03416
03417 nodval = bv*psat/rv/t;
03418
03419 nodval = nodval - rhov;
03420 break;
03421 }
03422 default:{
03423 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03424 exit(0);
03425 }
03426 }
03427 return(nodval);
03428 }
03429
03430
03431
03432
03433
03434
03435
03436
03437
03438
03439
03440
03441 double glasgowmatc::get_transmission_flux_mv(double bv,double t,double pg,double rhov,long bc,long ipp)
03442 {
03443 double flux;
03444
03445 switch (bc){
03446 case 30:{
03447 flux = (bv - rhov);
03448 break;
03449 }
03450 case 31:{
03451 double beta;
03452
03453 beta = f_beta(pg,t);
03454
03455 flux = beta*(bv - rhov);
03456 break;
03457 }
03458 case 32:{
03459 double beta,tinf;
03460
03461 tinf = f_tinf(Tp->time);
03462 beta = f_beta(pg,tinf);
03463
03464 flux = beta*(bv - rhov);
03465 break;
03466 }
03467 case 36:{
03468 double psat;
03469
03470 psat = f_psat(t);
03471
03472 flux = bv*psat/rv/t;
03473
03474 flux = flux - rhov;
03475 break;
03476 }
03477 default:{
03478 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03479 exit(0);
03480 }
03481 }
03482 return(flux);
03483 }
03484
03485
03486
03487
03488
03489
03490
03491
03492
03493
03494 double glasgowmatc::get_transmission_transcoeff_tt(double t,double pg,double rhov,long bc,long ipp)
03495 {
03496 double trc;
03497
03498 switch (bc){
03499 case 30:{
03500 trc=1.0;
03501 break;
03502 }
03503 case 31:{
03504 double hqr;
03505
03506 hqr = f_hqr(t,t);
03507
03508 trc=hqr;
03509 break;
03510 }
03511 case 32:{
03512 double hqr,tinf;
03513
03514 tinf = f_tinf(Tp->time);
03515 hqr = f_hqr(t,tinf);
03516
03517 trc=hqr;
03518 break;
03519 }
03520 default:{
03521 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03522 exit(0);
03523 }
03524 }
03525 return(trc);
03526 }
03527
03528
03529
03530
03531
03532
03533
03534
03535
03536
03537
03538
03539
03540 double glasgowmatc::get_transmission_nodval_tt(double bv,double t,double pg,double rhov,long bc,long ipp)
03541 {
03542 double nodval;
03543
03544 switch (bc){
03545 case 30:{
03546 nodval = bv;
03547 break;
03548 }
03549 case 31:{
03550 double hqr;
03551
03552 hqr = f_hqr(t,bv);
03553
03554 nodval = hqr*bv;
03555 break;
03556 }
03557 case 32:{
03558 double hqr,tinf;
03559
03560 tinf = f_tinf(Tp->time);
03561 hqr = f_hqr(t,tinf);
03562
03563 nodval = hqr*tinf;
03564 break;
03565 }
03566 default:{
03567 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03568 exit(0);
03569 }
03570 }
03571 return(nodval);
03572 }
03573
03574
03575
03576
03577
03578
03579
03580
03581
03582
03583
03584
03585
03586 double glasgowmatc::get_transmission_flux_tt(double bv,double t,double pg,double rhov,long bc,long ipp)
03587 {
03588 double flux;
03589
03590 switch (bc){
03591 case 30:{
03592 flux = (bv - t);
03593 break;
03594 }
03595 case 31:{
03596 double hqr;
03597
03598 hqr = f_hqr(t,bv);
03599
03600 flux = hqr*(bv - t);
03601 break;
03602 }
03603 case 32:{
03604 double hqr,tinf;
03605
03606 tinf = f_tinf(Tp->time);
03607 hqr = f_hqr(t,tinf);
03608
03609 flux = hqr*(tinf - t);
03610 break;
03611 }
03612 default:{
03613 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03614 exit(0);
03615 }
03616 }
03617 return(flux);
03618 }
03619
03620
03621
03622
03623
03624
03625
03626
03627
03628
03629
03630 double glasgowmatc::get_transmission_transcoeff_mt(double t,double pg,double rhov,long bc,long ipp)
03631 {
03632 double trc;
03633 double kvt,keff,hqr,tinf;
03634 double dav,fracg,rhoa,rhog;
03635
03636 dav = f_dav(rhov,pg,t);
03637 fracg = f_fracg(rhov,t);
03638 rhoa = f_rhoa(rhov,pg,t);
03639 rhog = f_rhog(rhov,pg,t);
03640
03641 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
03642 keff = f_keff (t);
03643
03644 switch (bc){
03645 case 30:{
03646
03647 trc=-1.0*kvt/keff;
03648 break;
03649 }
03650 case 31:{
03651 hqr = f_hqr(t,t);
03652
03653 trc=-1.0*kvt*hqr/keff;
03654 break;
03655 }
03656 case 32:{
03657 tinf = f_tinf(Tp->time);
03658 hqr = f_hqr(t,tinf);
03659
03660 trc=-1.0*kvt*hqr/keff;
03661 break;
03662 }
03663 case 36:{
03664 trc=0.0;
03665 break;
03666 }
03667 default:{
03668 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03669 exit(0);
03670 }
03671 }
03672
03673 return(trc);
03674 }
03675
03676
03677
03678
03679
03680
03681
03682
03683
03684
03685
03686
03687 double glasgowmatc::get_transmission_nodval_mt(double bv,double t,double pg,double rhov,long bc,long ipp)
03688 {
03689 double nodval;
03690 double kvt,keff,hqr,tinf;
03691 double dav,fracg,rhoa,rhog;
03692
03693 dav = f_dav(rhov,pg,t);
03694 fracg = f_fracg(rhov,t);
03695 rhoa = f_rhoa(rhov,pg,t);
03696 rhog = f_rhog(rhov,pg,t);
03697
03698 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
03699
03700 keff = f_keff (t);
03701
03702 switch (bc){
03703 case 30:{
03704
03705 nodval = -1.0*kvt*bv/keff;
03706 break;
03707 }
03708 case 31:{
03709 hqr = f_hqr(t,bv);
03710
03711 nodval = -1.0*kvt*hqr*bv/keff;
03712 break;
03713 }
03714 case 32:{
03715 tinf = f_tinf(Tp->time);
03716 hqr = f_hqr(t,tinf);
03717
03718 nodval = -1.0*kvt*hqr*tinf/keff;
03719 break;
03720 }
03721 case 36:{
03722 nodval=0.0;
03723 break;
03724 }
03725 default:{
03726 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03727 exit(0);
03728 }
03729 }
03730 return(nodval);
03731 }
03732
03733
03734
03735
03736
03737
03738
03739
03740
03741
03742
03743 double glasgowmatc::get_transmission_flux_mt(double bv,double t,double pg,double rhov,long bc,long ipp)
03744 {
03745 double flux;
03746 double kvt,keff,hqr,tinf;
03747 double dav,fracg,rhoa,rhog;
03748
03749 dav = f_dav(rhov,pg,t);
03750 fracg = f_fracg(rhov,t);
03751 rhoa = f_rhoa(rhov,pg,t);
03752 rhog = f_rhog(rhov,pg,t);
03753
03754 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
03755
03756 keff = f_keff (t);
03757
03758 switch (bc){
03759 case 30:{
03760
03761 flux = -1.0*kvt*(bv - t)/keff;
03762 break;
03763 }
03764 case 31:{
03765 hqr = f_hqr(t,bv);
03766
03767 flux = -1.0*kvt*hqr*(bv - t)/keff;
03768 break;
03769 }
03770 case 32:{
03771 tinf = f_tinf(Tp->time);
03772 hqr = f_hqr(t,tinf);
03773
03774 flux = -1.0*kvt*hqr*(tinf - t)/keff;
03775 break;
03776 }
03777 case 36:{
03778 flux=0.0;
03779 break;
03780 }
03781 default:{
03782 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03783 exit(0);
03784 }
03785 }
03786 return(flux);
03787 }
03788
03789
03790
03791
03792
03793
03794
03795
03796
03797
03798
03799 double glasgowmatc::get_othervalue(long compother,long ipp, double *r)
03800 {
03801 double other;
03802
03803 switch (compother){
03804 case 0:{
03805 other = r[0];
03806 break;
03807 }
03808 case 1:{
03809 other = r[1];
03810 break;
03811 }
03812 case 2:{
03813 other = r[2];
03814 break;
03815 }
03816 case 3:{
03817 double h;
03818 h = f_h(r[0],r[2]);
03819
03820 other = h;
03821 break;
03822 }
03823 case 4:{
03824 double s;
03825 s = f_s(r[0],r[2]);
03826
03827 other = s;
03828 break;
03829 }
03830 case 5:{
03831 double pv;
03832 pv = f_pv(r[0],r[2]);
03833
03834 other = pv;
03835 break;
03836 }
03837 case 6:{
03838 double pl;
03839 pl = f_pl(r[0],r[1],r[2]);
03840
03841 other = pl;
03842 break;
03843 }
03844 case 7:{
03845 double pc;
03846 pc = f_pc(r[0],r[2]);
03847
03848 other = pc;
03849 break;
03850 }
03851 default:{
03852 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
03853 }
03854 }
03855 return (other);
03856 }
03857
03858
03859
03860
03861
03862
03863 void glasgowmatc::print_othervalue_name(FILE *out,long compother)
03864 {
03865 switch (compother){
03866 case 0:{
03867 fprintf (out,"Vapour content (kg/m3) ");
03868 break;
03869 }
03870 case 1:{
03871 fprintf (out,"Gas pressure (Pa) ");
03872 break;
03873 }
03874 case 2:{
03875 fprintf (out,"Temperature (K) ");
03876 break;
03877 }
03878 case 3:{
03879 fprintf (out,"Relative humidity () ");
03880 break;
03881 }
03882 case 4:{
03883 fprintf (out,"Degree of saturation () ");
03884 break;
03885 }
03886 case 5:{
03887 fprintf (out,"Pore water vapor pressure (Pa)");
03888 break;
03889 }
03890 case 6:{
03891 fprintf (out,"Pore water pressure (Pa) ");
03892 break;
03893 }
03894 case 7:{
03895 fprintf (out,"Capilary pressure (Pa) ");
03896 break;
03897 }
03898 default:{
03899 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
03900 }
03901 }
03902 }
03903
03904
03905
03906
03907
03908
03909
03910
03911
03912 double glasgowmatc::kuu (double t,double pg,double rhov)
03913 {
03914 double kuu;
03915
03916 switch (model){
03917 case 1:{
03918 kuu = 0.0;
03919 break;
03920 }
03921 case 2:{
03922 kuu = 0.0;
03923 break;
03924 }
03925 default:{
03926 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
03927 }
03928 }
03929
03930 return (kuu);
03931 }
03932
03933
03934 double glasgowmatc::kut (double t,double pg,double rhov)
03935 {
03936 double kut;
03937
03938 switch (model){
03939 case 1:{
03940 kut = -1.0*f_alpha (t);
03941 break;
03942 }
03943 case 2:{
03944 kut = -1.0*f_alpha (t);
03945 break;
03946 }
03947 default:{
03948 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
03949 }
03950 }
03951
03952
03953
03954 return (kut);
03955 }
03956
03957
03958 double glasgowmatc::kut2 (double t,double pg,double rhov)
03959 {
03960 double kut,sln,dpcdt;
03961
03962 sln = f_sln(rhov,t);
03963 dpcdt = f_dpcdt(rhov,t);
03964
03965 switch (model){
03966 case 1:{
03967 kut = sln*dpcdt;
03968 break;
03969 }
03970 case 2:{
03971 kut = sln*dpcdt;
03972 break;
03973 }
03974 default:{
03975 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
03976 }
03977 }
03978
03979 kut = 0.0;
03980
03981 return (kut);
03982 }
03983
03984 double glasgowmatc::kup (double t,double pg,double rhov)
03985 {
03986 double kup;
03987
03988 switch (model){
03989 case 1:{
03990 kup = -1.0;
03991 break;
03992 }
03993 case 2:{
03994 kup = -1.0;
03995 break;
03996 }
03997 default:{
03998 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
03999 }
04000 }
04001
04002 kup = 0.0;
04003
04004 return (kup);
04005 }
04006
04007
04008 double glasgowmatc::kuv (double t,double pg,double rhov)
04009 {
04010 double kuv,sln,dpcdrhov;
04011
04012 sln = f_sln(rhov,t);
04013 dpcdrhov = f_dpcdrhov(rhov,t);
04014
04015 switch (model){
04016 case 1:{
04017 kuv = sln*dpcdrhov;
04018 break;
04019 }
04020 case 2:{
04021 kuv = sln*dpcdrhov;
04022 break;
04023 }
04024 default:{
04025 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
04026 }
04027 }
04028
04029 kuv = 0.0;
04030
04031 return (kuv);
04032 }
04033
04034
04035 double glasgowmatc::ktu ()
04036 {
04037 double ktu;
04038
04039 ktu = 0.0;
04040
04041 return (ktu);
04042 }
04043
04044 double glasgowmatc::kpu ()
04045 {
04046 double kpu;
04047
04048 kpu = 0.0;
04049
04050 return (kpu);
04051 }
04052
04053
04054 double glasgowmatc::kvu ()
04055 {
04056 double kvu;
04057
04058 kvu = 0.0;
04059
04060 return (kvu);
04061 }
04062
04063
04064 double glasgowmatc::cvu ()
04065 {
04066 double cvu;
04067
04068 cvu = 0.0;
04069
04070 return (cvu);
04071 }
04072
04073 double glasgowmatc::cpu ()
04074 {
04075 double cpu;
04076
04077 cpu = 0.0;
04078
04079 return (cpu);
04080 }
04081
04082
04083 double glasgowmatc::ctu ()
04084 {
04085 double ctu;
04086
04087 ctu = 0.0;
04088
04089 return (ctu);
04090 }
04091
04092
04093 double glasgowmatc::cuv ()
04094 {
04095 double cuv;
04096
04097 cuv = 0.0;
04098
04099 return (cuv);
04100 }
04101
04102 double glasgowmatc::cup ()
04103 {
04104 double cup;
04105
04106 cup = 0.0;
04107
04108 return (cup);
04109 }
04110
04111
04112 double glasgowmatc::cut ()
04113 {
04114 double cut;
04115
04116 cut = 0.0;
04117
04118 return (cut);
04119 }
04120
04121
04122
04123
04124
04125
04126 double glasgowmatc::emod()
04127 {
04128 double emod0;
04129
04130 emod0 = 3.0e+10;
04131
04132 return(emod0);
04133 }
04134
04135
04136
04137
04138
04139
04140
04141
04142
04143
04144 double glasgowmatc::fuv1(double t,double pg,double rhov)
04145 {
04146 return(0.0);
04147 }
04148
04149
04150
04151
04152
04153
04154
04155
04156
04157
04158
04159 double glasgowmatc::fup1(double t,double pg,double rhov)
04160 {
04161 return(0.0);
04162 }
04163
04164
04165
04166
04167
04168
04169
04170
04171
04172
04173 double glasgowmatc::fut1(double t,double pg,double rhov)
04174 {
04175 double fut1;
04176
04177
04178 switch (model){
04179 case 1:{
04180 fut1 = -1.0*f_alpha(t);
04181 break;
04182 }
04183 case 2:{
04184 fut1 = -1.0*f_alpha(t);
04185 break;
04186 }
04187 default:{
04188 fprintf (stderr,"\n\n unknown model type in glasgowmatc.cpp (%s, line %d).\n",__FILE__,__LINE__);
04189 }
04190 }
04191
04192
04193
04194 return(fut1);
04195 }