![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 22732
- Timestamp:
- 05/01/18 11:13:06 (7 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 11 edited
-
src/c/classes/Elements/Element.h (modified) (1 diff)
-
src/c/classes/Elements/Penta.h (modified) (1 diff)
-
src/c/classes/Elements/Seg.h (modified) (1 diff)
-
src/c/classes/Elements/Tetra.h (modified) (1 diff)
-
src/c/classes/Elements/Tria.cpp (modified) (2 diffs)
-
src/c/classes/Elements/Tria.h (modified) (1 diff)
-
src/c/classes/Params/TransientArrayParam.cpp (modified) (4 diffs)
-
src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (modified) (12 diffs)
-
src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h (modified) (1 diff)
-
src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
-
test/Archives/Archive470.arch (modified) ( previous)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Element.h
r22625 r22732 278 278 virtual int NumberofNodesVelocity(void)=0; 279 279 virtual void PicoUpdateBoxid(int* pmax_boxid_basin)=0; 280 virtual void PicoUpdateFirstBox()=0; 281 virtual void PicoUpdateNextBox(int loopboxid)=0; 280 virtual void PicoUpdateBox(int loopboxid)=0; 282 281 virtual void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0; 283 282 virtual int PressureInterpolation()=0; -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r22647 r22732 140 140 int NumberofNodesVelocity(void); 141 141 void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; 142 void PicoUpdateFirstBox(){_error_("not implemented yet");}; 143 void PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");}; 142 void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; 144 143 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding); 145 144 int PressureInterpolation(); -
issm/trunk-jpl/src/c/classes/Elements/Seg.h
r22647 r22732 129 129 int NumberofNodesVelocity(void){_error_("not implemented yet");}; 130 130 void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; 131 void PicoUpdateFirstBox(){_error_("not implemented yet");}; 132 void PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");}; 131 void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; 133 132 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");}; 134 133 int PressureInterpolation(void){_error_("not implemented yet");}; -
issm/trunk-jpl/src/c/classes/Elements/Tetra.h
r22647 r22732 137 137 int NumberofNodesVelocity(void); 138 138 void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; 139 void PicoUpdateFirstBox(){_error_("not implemented yet");}; 140 void PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");}; 139 void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; 141 140 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");}; 142 141 int PressureInterpolation(void); -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r22647 r22732 2817 2817 this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid)); 2818 2818 }/*}}}*/ 2819 void Tria::PicoUpdate FirstBox(){/*{{{*/2819 void Tria::PicoUpdateBox(int loopboxid){/*{{{*/ 2820 2820 2821 2821 if(!this->IsIceInElement() || !this->IsFloating()) return; 2822 2822 2823 int boxid;2824 this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);2825 if(boxid!=0) return;2826 2827 int basinid, maxbox, num_basins, numnodes, M;2828 IssmDouble time, gamma_T, overturning_coeff,thickness;2829 IssmDouble pressure, T_star,p_coeff, q_coeff,toc_farocean,soc_farocean;2830 IssmDouble* boxareas = NULL;2831 2832 /*Get variables*/2833 IssmDouble rhoi = this->GetMaterialParameter(MaterialsRhoIceEnum);2834 IssmDouble rhow = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);2835 IssmDouble earth_grav = this->GetMaterialParameter(ConstantsGEnum);2836 IssmDouble rho_star = 1033.; // kg/m^32837 IssmDouble nu = rhoi/rhow;2838 IssmDouble latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum);2839 IssmDouble c_p_ocean = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);2840 IssmDouble lambda = latentheat/c_p_ocean;2841 IssmDouble a = -0.0572; // K/psu2842 IssmDouble b = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum); //K2843 IssmDouble c = 7.77e-4;2844 IssmDouble alpha = 7.5e-5; // 1/K2845 IssmDouble Beta = 7.7e-4; // K2846 2847 /* Get parameters and inputs */2848 this->parameters->FindParam(&time,TimeEnum);2849 this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);2850 this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);2851 this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);2852 this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);2853 this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);2854 this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);2855 this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);2856 this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);2857 Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);2858 IssmDouble area_box1 = boxareas[basinid*maxbox+boxid];2859 IssmDouble g1 = area_box1*gamma_T;2860 IssmDouble s1 = soc_farocean/(nu*lambda);2861 2862 /*Define new inputs*/2863 IssmDouble basalmeltrates_shelf[NUMVERTICES];2864 IssmDouble potential_pressure_melting_point[NUMVERTICES];2865 IssmDouble Tocs[NUMVERTICES];2866 IssmDouble Socs[NUMVERTICES];2867 IssmDouble overturnings[NUMVERTICES];2868 2869 /* Start looping on the number of nodes and calculate ocean vars */2870 Gauss* gauss=this->NewGauss();2871 for(int i=0;i<NUMVERTICES;i++){2872 gauss->GaussVertex(i);2873 thickness_input->GetInputValue(&thickness,gauss);2874 pressure = (rhoi*earth_grav*1e-4)*thickness;2875 T_star = a*soc_farocean+b-c*pressure-toc_farocean;2876 p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));2877 q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));2878 2879 /* To avoid negatives under the square root */2880 if((0.25*pow(p_coeff,2)-q_coeff)<0){2881 q_coeff = 0.25*p_coeff*p_coeff;2882 }2883 2884 Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));2885 Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);2886 potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;2887 basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);2888 overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));2889 }2890 2891 this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);2892 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);2893 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);2894 this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum);2895 2896 /*Cleanup and return */2897 delete gauss;2898 xDelete<IssmDouble>(boxareas);2899 }2900 /*}}}*/2901 void Tria::PicoUpdateNextBox(int loopboxid){/*{{{*/2902 2903 if(!this->IsIceInElement() || !this->IsFloating()) return;2904 2905 2823 int boxid; 2906 2824 this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); 2907 2825 if(loopboxid!=boxid) return; 2908 2826 2909 int basinid, maxbox, num_basins, numnodes,M; 2910 IssmDouble gamma_T, overturning_coeff; 2911 IssmDouble thickness, toc_farocean, soc_farocean; 2827 int basinid, maxbox, num_basins, numnodes, M; 2828 IssmDouble gamma_T, overturning_coeff, thickness; 2912 2829 IssmDouble pressure, T_star,p_coeff, q_coeff; 2913 IssmDouble* toc_weighted_avg = NULL; 2914 IssmDouble* soc_weighted_avg = NULL; 2915 IssmDouble* overturning_weighted_avg = NULL; 2916 IssmDouble* boxareas = NULL; 2830 IssmDouble* boxareas = NULL; 2917 2831 2918 2832 /*Get variables*/ … … 2931 2845 IssmDouble Beta = 7.7e-4; // K 2932 2846 2933 /* Get parameters and inputs */2847 /* Get non-box-specific parameters and inputs */ 2934 2848 this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum); 2935 2849 this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum); 2936 2850 this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum); 2937 2851 this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); 2938 this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);2939 this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);2940 this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);2941 2852 this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); 2942 2853 this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); 2943 2944 2945 _assert_(basinid<=num_basins); 2854 Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); 2855 _assert_(basinid<=num_basins); 2856 2946 2857 IssmDouble area_boxi = boxareas[basinid*maxbox+boxid]; 2947 IssmDouble mean_toc = toc_weighted_avg[basinid];2948 IssmDouble mean_soc = soc_weighted_avg[basinid];2949 IssmDouble mean_overturning = overturning_weighted_avg[basinid];2950 2858 IssmDouble g1 = area_boxi*gamma_T; 2951 IssmDouble g2 = g2/(nu*lambda);2952 2859 2953 2860 IssmDouble basalmeltrates_shelf[NUMVERTICES]; 2954 IssmDouble potential_pressure_melting_point[NUMVERTICES]; 2955 IssmDouble Tocs[NUMVERTICES]; 2956 IssmDouble Socs[NUMVERTICES]; 2957 2958 /* Start looping on the number of nodes and calculate ocean vars */ 2959 Gauss* gauss=this->NewGauss(); 2960 for(int i=0;i<NUMVERTICES;i++){ 2961 gauss->GaussVertex(i); 2962 thickness_input->GetInputValue(&thickness,gauss); 2963 pressure = (rhoi*earth_grav*1.e-4)*thickness; 2964 T_star = a*mean_soc+b-c*pressure-mean_toc; 2965 Tocs[i] = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc)); 2966 Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda)); 2967 potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; 2968 basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); 2969 } 2970 2971 this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); 2972 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); 2973 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); 2974 2975 /*Cleanup and return */ 2976 delete gauss; 2861 IssmDouble potential_pressure_melting_point[NUMVERTICES]; 2862 IssmDouble Tocs[NUMVERTICES]; 2863 IssmDouble Socs[NUMVERTICES]; 2864 2865 /* First box calculations */ 2866 if(boxid==0){ 2867 /* Get box1 parameters and inputs */ 2868 IssmDouble time, toc_farocean, soc_farocean; 2869 this->parameters->FindParam(&time,TimeEnum); 2870 this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum); 2871 this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum); 2872 IssmDouble s1 = soc_farocean/(nu*lambda); 2873 IssmDouble overturnings[NUMVERTICES]; 2874 2875 /* Start looping on the number of verticies and calculate ocean vars */ 2876 Gauss* gauss=this->NewGauss(); 2877 for(int i=0;i<NUMVERTICES;i++){ 2878 gauss->GaussVertex(i); 2879 thickness_input->GetInputValue(&thickness,gauss); 2880 pressure = (rhoi*earth_grav*1e-4)*thickness; 2881 T_star = a*soc_farocean+b-c*pressure-toc_farocean; 2882 p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha)); 2883 q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha))); 2884 2885 /* To avoid negatives under the square root */ 2886 if((0.25*pow(p_coeff,2)-q_coeff)<0){ 2887 q_coeff = 0.25*p_coeff*p_coeff; 2888 } 2889 2890 Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff)); 2891 Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]); 2892 potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; 2893 basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); 2894 overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i])); 2895 } 2896 2897 this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); 2898 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); 2899 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); 2900 this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum); 2901 2902 /*Cleanup and return*/ 2903 delete gauss; 2904 } 2905 2906 /* Subsequent box calculations */ 2907 else { 2908 /* Get subsequent box parameters and inputs */ 2909 IssmDouble* toc_weighted_avg = NULL; 2910 IssmDouble* soc_weighted_avg = NULL; 2911 IssmDouble* overturning_weighted_avg = NULL; 2912 this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum); 2913 this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum); 2914 this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum); 2915 IssmDouble mean_toc = toc_weighted_avg[basinid]; 2916 IssmDouble mean_soc = soc_weighted_avg[basinid]; 2917 IssmDouble mean_overturning = overturning_weighted_avg[basinid]; 2918 IssmDouble g2 = g1/(nu*lambda); 2919 2920 /* Start looping on the number of verticies and calculate ocean vars */ 2921 Gauss* gauss=this->NewGauss(); 2922 for(int i=0;i<NUMVERTICES;i++){ 2923 gauss->GaussVertex(i); 2924 thickness_input->GetInputValue(&thickness,gauss); 2925 pressure = (rhoi*earth_grav*1.e-4)*thickness; 2926 T_star = a*mean_soc+b-c*pressure-mean_toc; 2927 Tocs[i] = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc)); 2928 Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda)); 2929 potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; 2930 basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); 2931 } 2932 2933 this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); 2934 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); 2935 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); 2936 2937 /*Cleanup and return*/ 2938 xDelete<IssmDouble>(toc_weighted_avg); 2939 xDelete<IssmDouble>(soc_weighted_avg); 2940 xDelete<IssmDouble>(overturning_weighted_avg); 2941 delete gauss; 2942 } 2943 2944 /*Cleanup and return*/ 2977 2945 xDelete<IssmDouble>(boxareas); 2978 xDelete<IssmDouble>(toc_weighted_avg);2979 xDelete<IssmDouble>(soc_weighted_avg);2980 xDelete<IssmDouble>(overturning_weighted_avg);2981 2946 } 2982 2947 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r22647 r22732 111 111 int NumberofNodesVelocity(void); 112 112 void PicoUpdateBoxid(int* pmax_boxid_basin); 113 void PicoUpdateFirstBox(); 114 void PicoUpdateNextBox(int loopboxid); 113 void PicoUpdateBox(int loopboxid); 115 114 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding); 116 115 int PressureInterpolation(); -
issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp
r22628 r22732 24 24 _assert_(in_values && in_time); 25 25 26 enum_type=in_enum_type;27 M=in_M; //Number of time steps28 N=in_N; //Number of rows29 interpolation=interpolation_on;26 this->enum_type=in_enum_type; 27 this->M=in_M; //Number of rows 28 this->N=in_N; //Number of timesteps 29 this->interpolation=interpolation_on; 30 30 31 values=xNew<IssmDouble>(M*N);31 this->values=xNew<IssmDouble>(M*N); 32 32 xMemCpy<IssmDouble>(values,in_values,M*N); 33 33 34 t imesteps=xNew<IssmDouble>(M);35 xMemCpy<IssmDouble>(timesteps,in_time, M);34 this->timesteps=xNew<IssmDouble>(N); 35 xMemCpy<IssmDouble>(timesteps,in_time,N); 36 36 } 37 37 /*}}}*/ … … 100 100 IssmDouble output; 101 101 bool found; 102 _assert_(row>=0 && row<this->M); 102 103 103 /*Ok, we have the time , go through the timesteps, and figure out which interval we104 /*Ok, we have the time and row, go through the timesteps, and figure out which interval we 104 105 *fall within. Then interpolate the values on this interval: */ 105 106 if(time<this->timesteps[0]){ … … 118 119 if(time==this->timesteps[i]){ 119 120 /*We are right on one step time: */ 120 output =this->values[row*this->N+i];121 output = this->values[row*this->N+i]; 121 122 found=true; 122 123 break; //we are done with the time interpolation. … … 124 125 else{ 125 126 if(this->timesteps[i]<time && time<this->timesteps[i+1]){ 126 /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */127 IssmDouble deltat =this->timesteps[i+1]-this->timesteps[i];128 IssmDouble alpha =(time-this->timesteps[i])/deltat;127 /*ok, we have the interval [i:i+1]. Interpolate linearly for now: */ 128 IssmDouble deltat = this->timesteps[i+1]-this->timesteps[i]; 129 IssmDouble alpha = (time-this->timesteps[i])/deltat; 129 130 if(interpolation==true) output=(1.0-alpha)*this->values[row*this->N+i] + alpha*this->values[row*this->N+i+1]; 130 131 else output=this->values[row*this->N+i]; -
issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp
r22625 r22732 11 11 int maxbox; 12 12 13 /*First, reset all melt to 0 */ 14 for(int i=0;i<femmodel->elements->Size();i++){ 15 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 16 int numvertices = element->GetNumberOfVertices(); 17 IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices); 18 element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); 19 xDelete<IssmDouble>(values); 20 } 21 22 femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); 23 UpdateBoxIdsPico(femmodel); 24 ComputeBoxAreasPico(femmodel); 25 UpdateFirstBoxPico(femmodel); 26 for(int i=1;i<maxbox;i++){ 27 ComputeAverageOceanvarsPico(femmodel, i-1); 28 UpdateNextBoxPico(femmodel, i); 13 /*First, reset all melt to 0 */ 14 for(int i=0;i<femmodel->elements->Size();i++){ 15 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 16 int numvertices = element->GetNumberOfVertices(); 17 IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices); 18 element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); 19 xDelete<IssmDouble>(values); 20 } 21 22 femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); 23 UpdateBoxIdsPico(femmodel); 24 ComputeBoxAreasPico(femmodel); 25 for(int i=0;i<maxbox;i++){ 26 UpdateBoxPico(femmodel,i); 27 ComputeAverageOceanvarsPico(femmodel,i); 29 28 } 30 29 }/*}}}*/ … … 32 31 void UpdateBoxIdsPico(FemModel* femmodel){/*{{{*/ 33 32 34 int i,k,numvertices,num_basins,maxbox,basinid; 35 IssmDouble dist_max,val; 36 int* nd=NULL; 33 int numvertices,num_basins,maxbox,basinid; 37 34 IssmDouble* dmax_basin=NULL; 38 35 IssmDouble* distances=NULL; … … 49 46 50 47 /*find maximum distance to grounding line per domain and per basin*/ 51 dist_max=-1.;52 for(i =0;i<num_basins;i++){dmax_basin[i]=-1;}53 for(i =0;i<femmodel->elements->Size();i++){48 IssmDouble dist_max=-1.; 49 for(int i=0;i<num_basins;i++){dmax_basin[i]=-1;} 50 for(int i=0;i<femmodel->elements->Size();i++){ 54 51 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 55 52 if(!element->IsIceInElement() || !element->IsFloating()) continue; … … 58 55 element->GetInputListOnVertices(&distances[0],DistanceToGroundinglineEnum); 59 56 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); 60 for( k=0; k<numvertices; k++){57 for(int k=0; k<numvertices; k++){ 61 58 if(fabs(distances[k])>dist_max){dist_max=fabs(distances[k]);} 62 59 if(fabs(distances[k])>dmax_basin[basinid]){dmax_basin[basinid]=fabs(distances[k]);} … … 66 63 67 64 /*Define maximum number of boxes per basin*/ 68 nd=xNew<int>(num_basins); 69 for(i=0; i<num_basins;i++){ 70 val=sqrt(dmax_basin[i]/dist_max)*(maxbox-1); 71 k=0; 72 while(k<val+.5){k++;} 65 int* nd=xNew<int>(num_basins); 66 for(int i=0; i<num_basins;i++){ 67 IssmDouble val=sqrt(dmax_basin[i]/dist_max)*(maxbox-1); 68 69 #ifdef _HAVE_ADOLC_ 70 /*Do not use floor when AD is on*/ 71 int k=0; while(k<val+.5){k++;} 73 72 nd[i]=k; 73 74 #else 75 nd[i]= reCast<int>(floor(val)); 76 #endif 74 77 } 75 78 76 79 /*Assign box numbers*/ 77 for(i =0;i<femmodel->elements->Size();i++){80 for(int i=0;i<femmodel->elements->Size();i++){ 78 81 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 79 82 element->PicoUpdateBoxid(nd); … … 114 117 115 118 }/*}}}*/ 116 void UpdateFirstBoxPico(FemModel* femmodel){/*{{{*/ 117 118 for(int i=0;i<femmodel->elements->Size();i++){ 119 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 120 element->PicoUpdateFirstBox(); 121 } 122 119 void UpdateBoxPico(FemModel* femmodel, int loopboxid){/*{{{*/ 120 for(int i=0;i<femmodel->elements->Size();i++){ 121 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 122 element->PicoUpdateBox(loopboxid); 123 } 123 124 }/*}}}*/ 124 125 void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid){/*{{{*/ 125 126 126 int i,k,p,num_basins, basinid, maxbox, M;127 int num_basins, basinid, maxbox, M; 127 128 IssmDouble area, toc, soc, overturning; 128 129 IssmDouble* boxareas=NULL; 129 130 IssmDouble* overturning_weighted_avg=NULL; 130 Element* element;131 Gauss* gauss;132 131 133 132 femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum); 134 133 femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); 135 M=maxbox*num_basins;136 134 femmodel->parameters->FindParam(&boxareas,&M, BasalforcingsPicoBoxAreaEnum); 137 IssmDouble* toc_weighted_avg =xNewZeroInit<IssmDouble>(num_basins);138 IssmDouble* soc_weighted_avg =xNewZeroInit<IssmDouble>(num_basins);139 IssmDouble* toc_sumweightedavg =xNewZeroInit<IssmDouble>(num_basins);140 IssmDouble* soc_sumweightedavg =xNewZeroInit<IssmDouble>(num_basins);141 IssmDouble* overturning_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins);135 IssmDouble* toc_weighted_avg = xNewZeroInit<IssmDouble>(num_basins); 136 IssmDouble* soc_weighted_avg = xNewZeroInit<IssmDouble>(num_basins); 137 IssmDouble* toc_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins); 138 IssmDouble* soc_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins); 139 IssmDouble* overturning_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins); 142 140 143 141 /* Compute Toc and Soc weighted avg (boxes 0 to n-1) */ 144 for(i =0;i<femmodel->elements->Size();i++){145 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));142 for(int i=0;i<femmodel->elements->Size();i++){ 143 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 146 144 if(!element->IsIceInElement() || !element->IsFloating()) continue; 147 145 int el_boxid; … … 153 151 154 152 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); 155 gauss=element->NewGauss(1); gauss->GaussPoint(0);153 Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); 156 154 tocs_input->GetInputValue(&toc,gauss); 157 155 socs_input->GetInputValue(&soc,gauss); … … 166 164 ISSM_MPI_Allreduce(soc_weighted_avg,soc_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); 167 165 168 for( k=0;k<num_basins;k++){169 p=k*maxbox+boxid;166 for(int k=0;k<num_basins;k++){ 167 int p=k*maxbox+boxid; 170 168 if(boxareas[p]==0) continue; 171 169 toc_sumweightedavg[k] = toc_sumweightedavg[k]/boxareas[p]; … … 179 177 if(boxid==0){ 180 178 overturning_weighted_avg=xNewZeroInit<IssmDouble>(num_basins); 181 for(i =0;i<femmodel->elements->Size();i++){182 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));179 for(int i=0;i<femmodel->elements->Size();i++){ 180 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 183 181 if(!element->IsIceInElement() || !element->IsFloating()) continue; 184 182 int el_boxid; … … 189 187 190 188 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); 191 gauss=element->NewGauss(1); gauss->GaussPoint(0);189 Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); 192 190 overturnings_input->GetInputValue(&overturning,gauss); 193 191 delete gauss; … … 199 197 ISSM_MPI_Allreduce(overturning_weighted_avg,overturning_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); 200 198 201 for(k=0;k<num_basins;k++){ 202 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); 203 p=k*maxbox+boxid; 204 if(boxareas[p]==0) continue; 199 for(int k=0;k<num_basins;k++){ 200 int p=k*maxbox+boxid; 201 if(boxareas[p]==0.) continue; 205 202 overturning_sumweightedavg[k] = overturning_sumweightedavg[k]/boxareas[p]; 206 203 } … … 216 213 xDelete<IssmDouble>(soc_weighted_avg); 217 214 xDelete<IssmDouble>(boxareas); 218 219 }/*}}}*/ 220 void UpdateNextBoxPico(FemModel* femmodel, int loopboxid){/*{{{*/ 221 for(int i=0;i<femmodel->elements->Size();i++){ 222 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 223 element->PicoUpdateNextBox(loopboxid); 224 } 225 226 }/*}}}*/ 215 }/*}}}*/ -
issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h
r22625 r22732 13 13 void UpdateBoxIdsPico(FemModel* femmodel); 14 14 void ComputeBoxAreasPico(FemModel* femmodel); 15 void Update FirstBoxPico(FemModel* femmodel);15 void UpdateBoxPico(FemModel* femmodel, int loopboxid); 16 16 void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid); 17 void UpdateNextBoxPico(FemModel* femmodel, int loopboxid);18 17 19 18 #endif /* _FloatingiceMeltingRatePicox_H*/ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
r22646 r22732 235 235 iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature"); 236 236 _assert_(M>1 && N>1); 237 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp, M,N));237 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,N,M)); 238 238 xDelete<IssmDouble>(transparam); 239 239 iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity"); 240 240 _assert_(M>1 && N>1); 241 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp, M,N));241 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,N,M)); 242 242 xDelete<IssmDouble>(transparam); 243 243 break;
Note:
See TracChangeset
for help on using the changeset viewer.
