Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26056) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26057) @@ -972,4 +972,5 @@ xDelete(newthickness); xDelete(cumdeltathickness); + xDelete(deltathickness); xDelete(newbase); xDelete(newsurface); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26056) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26057) @@ -2310,5 +2310,5 @@ } else{ - values[i]=meltratefactor[i]*tanh((base[i]-bed[i])/thresholdthickness)*(upperdepthmelt-base[i]); + values[i]=meltratefactor[i]*tanh((base[i]-bed[i])/thresholdthickness);//*(upperdepthmelt-base[i]); } } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26056) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26057) @@ -2961,5 +2961,5 @@ else{ // ice is close to calvinghaf: smooth transition between limitation and free calving. //heaviside: 0 for floating, 1 for grounded - heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI); + heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(M_PI*(groundedice-calvinghaf)/haf_eps)/(2.*M_PI); calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate; meltingrate=heaviside*meltingrate+0.; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26056) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26057) @@ -1512,12 +1512,12 @@ /*retrieve coordinates: lat,long,radius */ ::GetVerticesCoordinates(&llr_list[0][0],vertices,NUMVERTICES,spherical); - x1=llr_list[0][0]/180*PI; y1=llr_list[0][1]/180*PI; z1=llr_list[0][2]; - x2=llr_list[1][0]/180*PI; y2=llr_list[1][1]/180*PI; z2=llr_list[1][2]; - x3=llr_list[2][0]/180*PI; y3=llr_list[2][1]/180*PI; z3=llr_list[2][2]; + x1=llr_list[0][0]/180.*M_PI; y1=llr_list[0][1]/180.*M_PI; z1=llr_list[0][2]; + x2=llr_list[1][0]/180.*M_PI; y2=llr_list[1][1]/180.*M_PI; z2=llr_list[1][2]; + x3=llr_list[2][0]/180.*M_PI; y3=llr_list[2][1]/180.*M_PI; z3=llr_list[2][2]; /*compute great circle distance between vertices */ - arc12=2.*asin(sqrt(pow(sin((x2-x1)/2),2.0)+cos(x1)*cos(x2)*pow(sin((y2-y1)/2),2))); - arc23=2.*asin(sqrt(pow(sin((x3-x2)/2),2.0)+cos(x2)*cos(x3)*pow(sin((y3-y2)/2),2))); - arc31=2.*asin(sqrt(pow(sin((x1-x3)/2),2.0)+cos(x3)*cos(x1)*pow(sin((y1-y3)/2),2))); + arc12=2.*asin(sqrt(pow(sin(0.5*(x2-x1)),2)+cos(x1)*cos(x2)*pow(sin(0.5*(y2-y1)),2))); + arc23=2.*asin(sqrt(pow(sin(0.5*(x3-x2)),2)+cos(x2)*cos(x3)*pow(sin(0.5*(y3-y2)),2))); + arc31=2.*asin(sqrt(pow(sin(0.5*(x1-x3)),2)+cos(x3)*cos(x1)*pow(sin(0.5*(y1-y3)),2))); /*semi parameter */ @@ -2119,4 +2119,5 @@ /*Assign output pointers*/ + xDelete(timesteps); *pvalues=values; *ptimes=times; @@ -3641,5 +3642,5 @@ else{ // ice is close to calvinghaf: smooth transition between limitation and free calving. //heaviside: 0 for floating, 1 for grounded - heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI); + heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*M_PI); calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate; meltingrate=heaviside*meltingrate+0.; @@ -5173,13 +5174,13 @@ dx = x_element - xx[i]; dy = y_element - yy[i]; dist = sqrt(pow(dx,2)+pow(dy,2)); // distance between vertex and elemental centroid [m] - alpha = dist*360.0/(2*PI*earth_radius) * PI/180.0; // [in radians] 360 degree = 2*pi*earth_radius + alpha = dist*360.0/(2*M_PI*earth_radius) * M_PI/180.0; // [in radians] 360 degree = 2*pi*earth_radius /*Compute azimuths, both north and east components: */ - ang = PI/2 - atan2(dy,dx); // this is bearing angle! + ang = M_PI/2 - atan2(dy,dx); // this is bearing angle! Y_azim = cos(ang); X_azim = sin(ang); /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ - int index=reCast(alpha/PI*(M-1)); + int index=reCast(alpha/M_PI*(M-1)); U_elastic[i] += U_elastic_precomputed[index]; Y_elastic[i] += H_elastic_precomputed[index]*Y_azim; @@ -5187,14 +5188,14 @@ /*Add all components to the pUp solution vectors:*/ - U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i]; - Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i]; - X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i]; + U_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*U_elastic[i]; + Y_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*Y_elastic[i]; + X_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*X_elastic[i]; /*North-south, East-west components */ if (hemi == -1) { - ang2 = PI/2 - atan2(yy[i],xx[i]); + ang2 = M_PI/2 - atan2(yy[i],xx[i]); } else if (hemi == 1) { - ang2 = PI/2 - atan2(-yy[i],-xx[i]); + ang2 = M_PI/2 - atan2(-yy[i],-xx[i]); } if (hemi != 0){ @@ -5203,6 +5204,6 @@ N_elastic[i] += H_elastic_precomputed[index]*N_azim; E_elastic[i] += H_elastic_precomputed[index]*E_azim; - N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i]; - E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i]; + N_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/(4*M_PI*pow(earth_radius,2))*I*E_elastic[i]; } } @@ -5308,6 +5309,6 @@ if(longe>180)longe=longe-360; - late=late/180*PI; - longe=longe/180*PI; + late=late/180.*M_PI; + longe=longe/180.*M_PI; /*}}}*/ @@ -5344,5 +5345,5 @@ /*Compute alpha angle between centroid and current vertex: */ - lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; + lati=latitude[i]/180.*M_PI; longi=longitude[i]/180.*M_PI; delPhi=fabs(lati-late); delLambda=fabs(longi-longe); @@ -5362,5 +5363,5 @@ /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ - int index=reCast(alpha/PI*(M-1)); + int index=reCast(alpha/M_PI*(M-1)); U_elastic[i] += U_elastic_precomputed[index]; N_elastic[i] += H_elastic_precomputed[index]*N_azim; @@ -5455,9 +5456,9 @@ longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; - late=90-late; - if(longe>180)longe=(longe-180)-180; - - late=late/180*PI; - longe=longe/180*PI; + late=90.-late; + if(longe>180.)longe=(longe-180.)-180.; + + late=late/180.*M_PI; + longe=longe/180.*M_PI; /*}}}*/ re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0; @@ -5477,7 +5478,7 @@ * ALL in geographic coordinates * */ - dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; - dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; - dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; + dI_list[0] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; + dI_list[1] = -4*M_PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; + dI_list[2] = +4*M_PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; } else if(masks->isiceonly[this->lid]){ @@ -5492,7 +5493,7 @@ deltathickness_input->GetInputAverage(&I); - dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; - dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; - dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; + dI_list[0] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; + dI_list[1] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; + dI_list[2] = +4*M_PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; } @@ -5621,8 +5622,8 @@ /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ - lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; - delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda; + lati=latitude[i]/180.*M_PI; longi=longitude[i]/180.*M_PI; + delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - indices[i]=alpha/PI*reCast(M-1); + indices[i]=alpha/M_PI*reCast(M-1); index=reCast(indices[i]); @@ -6145,5 +6146,5 @@ /*element radius: */ - IssmDouble re=sqrt(area/PI); + IssmDouble re=sqrt(area/M_PI); /*figure out gravity center of our element: */ Index: /issm/trunk-jpl/src/c/main/esmfbinders.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 26056) +++ /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 26057) @@ -153,4 +153,8 @@ } /*}}}*/ + /*TODO: we need 2 initialize routines, the second one will be empty for now + * In: ESMF config, ESMF Field bundle + */ + /*FISOC*/ void InitializeISSM_FISOC(int argc, char** argv,ISSM_MPI_Comm comm_init){ /*{{{*/ @@ -165,5 +169,5 @@ /*Initialize ESMC Mesh*/ int pdim; /*parametric dimension is the same as the domain dimensions */ - int sdim = 3; /*coordinates of each vertex is always 2 (just x,y for now) */ + int sdim = 2; /*coordinates of each vertex is always 2 (just x,y for now) */ ESMC_CoordSys_Flag coordsys = ESMC_COORDSYS_CART; /*Cartesian coordinate system by default */ femmodel->parameters->FindParam(&pdim,DomainDimensionEnum); @@ -171,9 +175,7 @@ if(rc!=ESMF_SUCCESS) _error_("could not create EMSC_Mesh"); - /*How to install ESMG with lib/libesmf.so - * do indices need to be 1 based? - * ESMF installation, how can we get everything in /lib + /* * What do we do with vertices at the boundary, declare twice? - * Do we need a restart file or save femmodel somewhere?*/ + * */ /*Add nodes (which are ISSM Vertices)*/ @@ -182,5 +184,5 @@ int *nodeOwner = xNew(numnodes); IssmDouble *nodeCoord = xNew(sdim*numnodes); - for (int i=0;ielements->Size();i++){ + for (int i=0;ivertices->Size();i++){ Vertex* vertex = xDynamicCast(femmodel->vertices->GetObjectByOffset(i)); nodeId[i] = vertex->Sid()+1; @@ -213,11 +215,4 @@ xDelete(elemConn); - /*DO we need to create fields here? https://earthsystemmodeling.org/docs/nightly/develop/ESMC_crefdoc/node5.html#SECTION05024400000000000000*/ - //ESMC_InterArrayInt *gridToFieldMap, // in - //ESMC_InterArrayInt *ungriddedLBound, // in - //ESMC_InterArrayInt *ungriddedUBound, // in - //ESMC_Field esmf_shelfbase = ESMC_FieldCreateMeshTypeKind(mesh,ESMC_TYPEKIND_R8,???,gridToFieldMap,ungriddedLBound,ungriddedUBound,"ShelfTopography",&rc); - if(rc!=ESMF_SUCCESS) _error_("could not create EMSC_Field"); - /*Create restart file for later */ femmodel->Restart();