Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 27239) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 27240) @@ -160,5 +160,5 @@ T = T - CtoK; // convert dT from degC/m to degC/cm - dT = dT/100; + dT = dT/100.0; // temperature coefficient F @@ -410,5 +410,5 @@ // determine new grain size (mm) - gsz[i] = 1e-1*(gdn[i]/.99+(1-1*gdn[i]/.99)*(gsp[i]/.99*3+(1-gsp[i]/.99)*4)); + gsz[i] = max(1e-1*(gdn[i]/.99+(1.0-1.0*gdn[i]/.99)*(gsp[i]/.99*3.0+(1.0-gsp[i]/.99)*4.0)),Gdntol*2.0); } @@ -756,5 +756,5 @@ IssmDouble zT=0.0; IssmDouble zQ=0.0; - IssmDouble zratio=1; + IssmDouble zratio=1.0; IssmDouble dt=0.0; IssmDouble max_fdt=0.0; @@ -815,5 +815,5 @@ //zT and zQ are percentage of z0 (Foken 2008) - zratio=10; + zratio=10.0; zT=z0/zratio; zQ=z0/zratio; @@ -886,5 +886,5 @@ // NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability dt=1e12; - for(int i=0;i 1/2 @@ -977,5 +977,5 @@ // calculate the Bulk Richardson Number (Ri) - Ri = pow(100000/pAir,0.286)*(2.0*9.81*(Ta - Ts)) / (Tz*(Ta + Ts)* pow(V/(Vz),2.0)); + Ri = pow(100000./pAir,0.286)*(2.0*9.81*(Ta - Ts)) / (Tz*(Ta + Ts)* pow(V/(Vz),2.0)); IssmDouble PhiM; @@ -993,5 +993,5 @@ } else { - coefM =pow (1.0-18*Ri,-0.25); + coefM =pow (1.0-18.0*Ri,-0.25); } @@ -1013,11 +1013,11 @@ // if stable //coefM = pow(1.0-5.0*Ri,2.0); //Fitzpatrick et al., 2017, from Brock et al., 2010 - coefM=1+5.3*min((Ri/(1-5*Ri)),0.5); - coefH=1+8*min((Ri/(1-5*Ri)),0.5); + coefM=1.0+5.3*min((Ri/(1.0-5.0*Ri)),0.5); + coefH=1.0+8.0*min((Ri/(1.0-5.0*Ri)),0.5); } else { //coefM =pow(1.0-16.0*max(Ri,-1.0),0.75); //Fitzpatrick et al., 2017, from Brock et al., 2010 - coefM=pow(1-19*max(Ri/1.5,-2.0),-0.25); - coefH=0.95*pow(1-11.6*max(Ri/1.5,-2.0),-0.5); + coefM=pow(1.0-19.0*max(Ri/1.5,-2.0),-0.25); + coefH=0.95*pow(1.0-11.6*max(Ri/1.5,-2.0),-0.5); } @@ -1032,10 +1032,10 @@ if (Ri > 0.0+Ttol){ // if stable - coefM=1+15*Ri*pow(1+Ri,1/2); - coefH=1; + coefM=1.0+15.0*Ri*pow(1.0+Ri,1./2.); + coefH=1.0; } else { - coefM=pow(1-15*Ri/(1+75*pow(0.4/log(Tz/zT),2)*pow(Tz/zT*fabs(Ri),1/2)),-1); - coefH=1; + coefM=pow(1.0-15.0*Ri/(1.0+75.0*pow(0.4/log(Tz/zT),2)*pow(Tz/zT*fabs(Ri),1./2.)),-1); + coefH=1.0; } @@ -1046,15 +1046,15 @@ else { IssmDouble a1=1.0; - IssmDouble b1=2/3; + IssmDouble b1=2.0/3.0; IssmDouble c1=5.0; IssmDouble d1=0.35; - IssmDouble PhiMz; - IssmDouble PhiHz; + IssmDouble PhiMz=0.0; + IssmDouble PhiHz=0.0; IssmDouble PhiMz0=0.0; IssmDouble PhiHzT=0.0; IssmDouble PhiHzQ=0.0; - IssmDouble zL; - IssmDouble zLT; - IssmDouble zLM; + IssmDouble zL=0.0; + IssmDouble zLT=0.0; + IssmDouble zLM=0.0; if (Ri > 0.0+Ttol){ @@ -1062,5 +1062,5 @@ if(Ri < 0.2-Ttol){ - zL = Ri/(1-5*Ri); + zL = Ri/(1.0-5.0*Ri); } else{ @@ -1068,12 +1068,12 @@ } //zL = min(zL, 0.5); //Sjoblom, 2014 - zLM=max(zL/Vz*z0,1e-3); - zLT=max(zL/Tz*zT,1e-3); + zLM=max(zL/Vz*z0,1.e-3); + zLT=max(zL/Tz*zT,1.e-3); //Ding et al. 2020, from Beljaars and Holtslag (1991) - PhiMz=-1*(a1*zL + b1*(zL-c1/d1)*exp(-1*d1*zL) + b1*c1/d1); - PhiHz=-1*(pow(1+2*a1*zL/3,1.5) + b1*(zL-c1/d1)*exp(-1*d1*zL) + b1*c1/d1 - 1); - PhiMz0=-1*(a1*zLM + b1*(zLM-c1/d1)*exp(-1*d1*zLM) + b1*c1/d1); - PhiHzT=-1*(pow(1+2*a1*zLT/3,1.5) + b1*(zLT-c1/d1)*exp(-1*d1*zLT) + b1*c1/d1 - 1); + PhiMz=-1.*(a1*zL + b1*(zL-c1/d1)*exp(-1.*d1*zL) + b1*c1/d1); + PhiHz=-1.*(pow(1.+2.*a1*zL/3.,1.5) + b1*(zL-c1/d1)*exp(-1.*d1*zL) + b1*c1/d1 - 1.0); + PhiMz0=-1.*(a1*zLM + b1*(zLM-c1/d1)*exp(-1.*d1*zLM) + b1*c1/d1); + PhiHzT=-1.*(pow(1.+2.*a1*zLT/3.,1.5) + b1*(zLT-c1/d1)*exp(-1.*d1*zLT) + b1*c1/d1 - 1.0); PhiHzQ=PhiHzT; @@ -1091,9 +1091,9 @@ if (true){ //Sjoblom, 2014 - xm=pow(1.0-19*zL,-0.25); - PhiMz=2*log((1+xm)/2.0) + log((1+pow(xm,2))/2.0) - 2*atan(xm) + Pi/2; + xm=pow(1.0-19.0*zL,-0.25); + PhiMz=2.0*log((1.+xm)/2.0) + log((1.+pow(xm,2))/2.0) - 2.*atan(xm) + Pi/2.; xh=0.95*pow(1.0-11.6*zL,-0.5); - PhiHz=2*log((1+pow(xh,2))/2.0); + PhiHz=2.0*log((1.0+pow(xh,2))/2.0); } else{ //Ding et al., 2020 @@ -1101,9 +1101,9 @@ xmM=pow(1.0-16*zLM,0.25); xmT=pow(1.0-16*zLT,0.25); - PhiMz=2*log((1+xm)/2.0) + log((1+pow(xm,2))/2.0) - 2*atan(xm) + Pi/2; - PhiMz0=2*log((1+xmM)/2.0) + log((1+pow(xmM,2))/2.0) - 2*atan(xmM) + Pi/2; - - PhiHz=2*log((1+pow(xm,2))/2.0); - PhiHzT=2*log((1+pow(xmT,2))/2.0); + PhiMz=2.0*log((1.+xm)/2.0) + log((1.+pow(xm,2))/2.0) - 2.0*atan(xm) + Pi/2.0; + PhiMz0=2.0*log((1.+xmM)/2.0) + log((1.+pow(xmM,2))/2.0) - 2.0*atan(xmM) + Pi/2.0; + + PhiHz=2.0*log((1.+pow(xm,2))/2.0); + PhiHzT=2.0*log((1.+pow(xmT,2))/2.0); PhiHzQ=PhiHzT; @@ -1121,5 +1121,5 @@ //// Sensible Heat // calculate the sensible heat flux [W m-2](Patterson, 1998) - shf = dAir * C * CA * (Ta - Ts) * pow(100000/pAir,0.286); + shf = dAir * C * CA * (Ta - Ts) * pow(100000./pAir,0.286); // adjust using Monin-Obukhov stability theory @@ -1147,5 +1147,5 @@ // Latent heat flux [W m-2] - lhf = C * L * (eAir - eS) / (461.9*(Ta+Ts)/2); + lhf = C * L * (eAir - eS) / (461.9*(Ta+Ts)/2.0); // adjust using Monin-Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface, @@ -1521,5 +1521,5 @@ gdnNew = min(max(1.29 - 0.17*V,0.20),1.0); gspNew = min(max(0.08*V + 0.38,0.5),0.9); - reNew=max(1e-1*(gdnNew/.99+(1-1*gdnNew/.99)*(gspNew/.99*3+(1-gspNew/.99)*4))/2,Gdntol); + reNew=max(1e-1*(gdnNew/.99+(1.0-1.0*gdnNew/.99)*(gspNew/.99*3.0+(1.0-gspNew/.99)*4.0))/2.0,Gdntol); break; @@ -1579,5 +1579,5 @@ gdn[0] = dfall; gsp[0] = sfall; - re[0] = max(1e-1*(gdn[0]/.99+(1-1*gdn[0]/.99)*(gsp[0]/.99*3+(1-gsp[0]/.99)*4))/2,Gdntol); + re[0] = max(1e-1*(gdn[0]/.99+(1.0-1.0*gdn[0]/.99)*(gsp[0]/.99*3.0+(1.0-gsp[0]/.99)*4.0))/2.0,Gdntol); } } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27239) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27240) @@ -34,5 +34,5 @@ void albedo(IssmDouble** a, IssmDouble** adiff, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid); void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble dswdiff, IssmDouble as, IssmDouble asdiff, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); -void thermo(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, int tcIdx, int eIdx, IssmDouble teValue, IssmDouble dulwrfValue, IssmDouble teThresh, IssmDouble Ws, IssmDouble dt0, IssmDouble dzMin, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT, bool isdeltaLWup); +void thermo(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, int tcIdx, int eIdx, IssmDouble teValue, IssmDouble dulwrfValue, IssmDouble teThresh, IssmDouble Ws, IssmDouble dt0, IssmDouble dzMin, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT, bool isdeltaLWup); void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* pRa, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid); void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pF, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble Ra, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid);