Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 12292) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 12293) @@ -484,4 +484,5 @@ PatersonEnum, ArrheniusEnum, + SurfaceforcingsIspddEnum, /*}}}*/ MaximumNumberOfEnums Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12292) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12293) @@ -452,4 +452,5 @@ case PatersonEnum : return "Paterson"; case ArrheniusEnum : return "Arrhenius"; + case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd"; default : return "unknown"; Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 12292) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 12293) @@ -89,4 +89,5 @@ parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum)); parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum)); /*some parameters that did not come with the iomodel: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 12292) +++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 12293) @@ -20,4 +20,5 @@ int stabilization; bool dakota_analysis; + bool ispdd; /*Fetch data needed: */ @@ -27,4 +28,5 @@ iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); iomodel->FetchData(1,MeshElementsEnum); + iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum); /*Update elements: */ @@ -44,4 +46,5 @@ iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum); iomodel->FetchDataToInput(elements,MaskElementonwaterEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); @@ -66,5 +69,9 @@ iomodel->FetchDataToInput(elements,TemperatureEnum); } - + if(ispdd){ + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); + } /*Free data: */ iomodel->DeleteData(1,MeshElementsEnum); Index: /issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp =================================================================== --- /issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (revision 12292) +++ /issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (revision 12293) @@ -20,5 +20,5 @@ // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME) // OUTPUT: mass-balance (m/yr ice): agd(NA) -// mean annual surface temperature (degrees C): Tsurf(NA) +// mean annual sssurface temperature (degrees C): Tsurf(NA) int i, it, jj, itm; @@ -26,6 +26,6 @@ double signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - double signormc; // sigma of the temperature distribution for cloudy day - double siglimc=0, siglim0, siglim0c; + double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day + double siglimc, siglim0, siglim0c; double tstep, tsint, tint, tstepc; int NPDMAX = 1504, NPDCMAX = 1454; @@ -43,8 +43,4 @@ pds=(double*)xmalloc(NPDCMAX*sizeof(double)+1); - - // PDD constant - siglim = 2.5*signorm; - // initialize PDD (creation of a lookup table) tstep = 0.1; @@ -53,4 +49,9 @@ snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0))); siglim = 2.5*signorm; + siglimc = 2.5*signormc + siglim0 = siglim/DT + 0.5 + siglim0c = siglimc/DT + 0.5 + PDup = siglimc+PDCUT + itm = (int)(2*siglim/DT + 1.5); Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12292) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12293) @@ -462,4 +462,5 @@ else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; + else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else stage=5; } Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 12292) +++ /issm/trunk/src/c/modules/modules.h (revision 12293) @@ -93,4 +93,5 @@ #include "./ConstraintsStatex/ConstraintsStatex.h" #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h" +#include "./PositiveDegreeDayx/PositiveDegreeDayx.h" #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h" #include "./Dakotax/Dakotax.h" @@ -122,3 +123,4 @@ #include "./VecMergex/VecMergex.h" + #endif Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 12292) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 12293) @@ -2434,5 +2434,6 @@ // PDD and PD constants and variables double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - double siglimc=0, siglim0, siglim0c; + double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day + double siglimc, siglim0, siglim0c; double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) double DT = 0.02; @@ -2480,11 +2481,11 @@ /*Get material parameters :*/ rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - density=rho_ice/rho_water; + //rho_freshwater=matpar->GetRhoWater(); - sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months + sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months /*PDD constant*/ siglim = 2.5*signorm; + siglimc = 2.5*signormc siglim0 = siglim/DT + 0.5; siglim0c = siglimc/DT + 0.5; Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 12292) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 12293) @@ -2238,5 +2238,5 @@ int i,iqj,imonth; - double agd[NUMVERTICES]; // surface and basal + double agd[NUMVERTICES]; // surface mass balance double saccu[NUMVERTICES] = {0}; // yearly surface accumulation double smelt[NUMVERTICES] = {0}; // yearly melt @@ -2247,9 +2247,9 @@ double sconv; //rhow_rain/rhoi / 12 months - double rho_water,rho_ice,density; - double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. - double desfac = 0.5; //desert elevation factor - double s0p[NUMVERTICES]={0}; //should be set to elevation from precip source - double s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source + double rho_water,rho_ice,density; + double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper + double desfac = 0.5; // desert elevation factor + double s0p[NUMVERTICES]={0}; // should be set to elevation from precip source + double s0t[NUMVERTICES]={0}; // should be set to elevation from temperature source double st; // elevation between altitude of the temp record and current altitude double sp; // elevation between altitude of the prec record and current altitude @@ -2258,5 +2258,6 @@ // PDD and PD constants and variables double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - double siglimc=0, siglim0, siglim0c; + double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day + double siglimc, siglim0, siglim0c; double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) double DT = 0.02; @@ -2275,5 +2276,5 @@ double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] - double t[NUMVERTICES][12],prec[NUMVERTICES][12]; + double t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12]; double deltm=1/12; int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11}; @@ -2304,11 +2305,11 @@ /*Get material parameters :*/ rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - density=rho_ice/rho_water; + //rho_freshwater=matpar->GetRhoWater(); - sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months + sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months /*PDD constant*/ siglim = 2.5*signorm; + siglimc = 2.5*signormc siglim0 = siglim/DT + 0.5; siglim0c = siglimc/DT + 0.5; @@ -2335,5 +2336,5 @@ else {q = 1.0;} - qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent + qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent per month qmpt= q*prec[i][imonth]*sconv; qmp[i]= qmp[i] + qmpt; Index: /issm/trunk/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 12292) +++ /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 12293) @@ -16,4 +16,5 @@ /*parameters: */ bool save_results; + bool ispdd; /*activate formulation: */ @@ -22,8 +23,14 @@ /*recover parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum); + + if(ispdd){ + _printf_(VerboseSolution()," call positive degree day module\n"); + PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } _printf_(VerboseSolution()," call computational core\n"); solver_linear(femmodel); - + if(save_results){ _printf_(VerboseSolution()," saving results\n"); Index: /issm/trunk/src/m/classes/clusters/acenet.m =================================================================== --- /issm/trunk/src/m/classes/clusters/acenet.m (revision 12293) +++ /issm/trunk/src/m/classes/clusters/acenet.m (revision 12293) @@ -0,0 +1,165 @@ +%ACENET class definition +% +% Usage: +% cluster=acenet(); +% cluster=acenet('np',3); +% cluster=acenet('np',3,'login','username'); + +classdef acenet + properties (SetAccess=public) + % {{{1 + name='glacdyn.ace-net.ca' + %name='placentia.ace-net.ca' + %name='brasdor.ace-net.ca' + login='klemorza'; + np=10; + port=0; + queue='shortq'; + time=10; + % codepath='/usr/local/issm-r11321/bin'; % this one is for issm on acenet global + codepath='/home/klemorza/issm/bin'; % this one is for issm on my acenet directory + executionpath='/home/klemorza/issm/execution'; + %}}} + end + methods + function cluster=acenet(varargin) % {{{1 + %use provided options to change fields + options=pairoptions(varargin{:}); + + %initialize cluster using user settings if provided + if (exist([cluster.name '_settings'])==2), eval([cluster.name '_settings']); end + + %OK get other fields + cluster=AssignObjectFields(pairoptions(varargin{:}),cluster); + end + %}}} + function disp(cluster) % {{{1 + % display the object + disp(sprintf('class ''%s'' object ''%s'' = ',class(cluster),inputname(1))); + disp(sprintf(' name: %s',cluster.name)); + disp(sprintf(' login: %s',cluster.login)); + disp(sprintf(' np: %i',cluster.np)); + disp(sprintf(' port: %i',cluster.port)); + disp(sprintf(' queue: %s',cluster.queue)); + disp(sprintf(' time: %i',cluster.time)); + disp(sprintf(' codepath: %s',cluster.codepath)); + disp(sprintf(' executionpath: %s',cluster.executionpath)); + end + %}}} + function checkconsistency(cluster,md,solution,analyses) % {{{1 + + available_queues={'debug','shortq','longq'}; + queue_requirements_time=[60*1 60*3 60*17]; + queue_requirements_np=[32 128 256]; + + QueueRequirements(available_queues,queue_requirements_time,queue_requirements_np,cluster.queue,cluster.np,cluster.time) + end + %}}} + function BuildQueueScript(cluster,md) % {{{1 + + %retrieve parameters + modelname=md.miscellaneous.name; + solution=md.private.solution; + + %open file for writing: + fid=fopen([modelname '.queue'],'w'); + + %write instructions for launching a job on the cluster + % fprintf(fid,'#!/bin/sh\n'); + fprintf(fid,'#!/bin/bash\n'); + %fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s 2> %s.errlog >%s.outlog ',... + % cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname,modelname,modelname); + fprintf(fid,'#$ -cwd\n'); + fprintf(fid,'#$ -N issm\n'); + %fprintf(fid,'#$ -l h_rt=%i\n',cluster.time); + fprintf(fid,'#$ -l h_rt=10:0:0\n'); + %fprintf(fid,'#$ -l h_rt=1:0:0,test=true\n'); + fprintf(fid,'#$ -pe ompi* %i\n',cluster.np); + fprintf(fid,'#$ -j y\n'); + fprintf(fid,'module purge\n'); + fprintf(fid,'module load gcc openmpi/gcc\n'); + fprintf(fid,'module load issm\n'); + fprintf(fid,'\n'); + %fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s',cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname); + fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s 2> %s.errlog >%s.outlog ',... + cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname,modelname,modelname); + + % fprintf(fid,'#PBS -l select=%i:ncpus=1\n',cluster.np); + % fprintf(fid,'#PBS -N %s\n',modelname); + % fprintf(fid,'#PBS -l walltime=%i\n',time*60); %walltime is in seconds. + % fprintf(fid,'#PBS -q %s\n',queue); + % fprintf(fid,'#PBS -o %s.outlog \n',modelname); + % fprintf(fid,'#PBS -e %s.errlog \n',modelname); + % fprintf(fid,'export PBS_O_WORKDIR=%s\n',cluster.executionpath); + % fprintf(fid,'cd $PBS_O_WORKDIR\n'); + % fprintf(fid,'export OMP_NUM_THREADS=1\n'); + % fprintf(fid,'ulimit -s unlimited\n'); + % fprintf(fid,'ulimit -c 0\n'); + % fprintf(fid,'/opt/mpich/gm/intel10.1/bin/mpiexec -np %i %s/issm.exe %s %s %s',cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname); + + %close file + fclose(fid); + + end + %}}} + function LaunchQueueJob(cluster,md,options)% {{{1 + %retrieve parameters + modelname=md.miscellaneous.name; + solution=md.private.solution; + + %lauch command, to be executed via ssh + launchcommand=['cd ' cluster.executionpath ' && rm -rf ./' md.private.runtimename ' && mkdir ' md.private.runtimename ... + ' && cd ' md.private.runtimename ' && mv ../' md.private.runtimename '.tar.gz ./ && tar -zxf ' md.private.runtimename '.tar.gz && qsub ' modelname '.queue ']; + + if ~strcmpi(options.batch,'yes'), + + %compress the files into one zip. + compressstring=['tar -zcf ' md.private.runtimename '.tar.gz ' md.miscellaneous.name '.bin ' md.miscellaneous.name '.queue ' md.miscellaneous.name '.petsc ']; + if md.qmu.isdakota, + compressstring=[compressstring md.miscellaneous.name '.qmu.in']; + end + system(compressstring); + + disp('uploading input file and queueing script'); + issmscpout(md.cluster.name,md.cluster.executionpath,md.cluster.login,md.cluster.port,{[md.private.runtimename '.tar.gz']}); + + disp('launching solution sequence on remote cluster'); + issmssh(md.cluster.name,md.cluster.login,md.cluster.port,launchcommand); + + else + disp('batch mode requested: not launching job interactively'); + disp('launch solution sequence on remote cluster by hand'); + end + + end + %}}} + function Download(cluster,md)% {{{1 + + %some check + if isempty(md.private.runtimename), + error('pfe Download error message: supply runtime name for results to be loaded!'); + end + + %Figure out the directory where all the files are in: + directory=[cluster.executionpath '/' md.private.runtimename '/']; + + %What packages are we picking up from remote cluster + packages={[md.miscellaneous.name '.outlog'],[md.miscellaneous.name '.errlog']}; + %packages={}; + if md.qmu.isdakota, + packages{end+1}=[md.miscellaneous.name '.qmu.err']; + packages{end+1}=[md.miscellaneous.name '.qmu.out']; + if isfield(md.qmu.params,'tabular_graphics_data'), + if md.qmu.params.tabular_graphics_data==true, + packages{end+1}='dakota_tabular.dat'; + end + end + else + packages{end+1}=[md.miscellaneous.name '.outbin']; + end + + %copy files from cluster to present directory + issmscpin(cluster.name, cluster.login, cluster.port, directory, packages); + end %}}} + end +end Index: /issm/trunk/src/m/classes/surfaceforcings.m =================================================================== --- /issm/trunk/src/m/classes/surfaceforcings.m (revision 12292) +++ /issm/trunk/src/m/classes/surfaceforcings.m (revision 12293) @@ -8,4 +8,5 @@ precipitation = NaN; mass_balance = NaN; + ispdd = NaN; end methods @@ -19,4 +20,7 @@ end % }}} function obj = setdefaultparameters(obj) % {{{ + + %pdd methode not use in default mode + obj.ispdd=0; end % }}} @@ -29,4 +33,7 @@ checkfield(md,'surfaceforcings.mass_balance','size',[md.mesh.numberofvertices 1],'NaN',1); end + if (ismember(PrognosticAnalysisEnum,analyses) & md.surfaceforcings.ispdd), + checkfield(md,'surfaceforcings.ispdd','numel',1,'values',[0 1]); + end end % }}} function disp(obj) % {{{ @@ -35,4 +42,5 @@ fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]'); fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]'); + fielddisplay(obj,'ispdd','is pdd activated (0 or 1, default is 0)'); end % }}} @@ -40,4 +48,6 @@ WriteData(fid,'object',obj,'fieldname','precipitation','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','mass_balance','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','ispdd','format','Boolean'); + end % }}} end Index: /issm/trunk/src/m/enum/EnumToModelField.m =================================================================== --- /issm/trunk/src/m/enum/EnumToModelField.m (revision 12293) +++ /issm/trunk/src/m/enum/EnumToModelField.m (revision 12293) @@ -0,0 +1,26 @@ +function string=EnumToModelField(enum) +%ENUMTOMODELFIELD - output string of model field associated to enum +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% string=EnumToModelField(enum) + +switch enum, + + case ThicknessEnum(), string='thickness'; return + case FrictionCoefficientEnum(), string='drag_coefficient'; return + case MaterialsRheologyBEnum(), string='rheology_B'; return + case MaterialsRheologyBbarEnum(), string='rheology_B'; return + case BalancethicknessThickeningRateEnum: string='dhdt'; return + case VxEnum(), string='vx'; return + case InversionVxObsEnum(), string='vx_obs'; return + case VyEnum(), string='vy'; return + case InversionVyObsEnum(), string='vy_obs'; return + case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return + case SurfaceforcingsMassBalanceEnum(), string='surface_mass_balance'; return + otherwise, error(['Enum ' num2str(enum) ' not found associated to any model field']); + +end Index: /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 12292) +++ /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 12293) @@ -58,4 +58,5 @@ %segment on Neumann (Ice Front) pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); + if (md.mesh.dimension==2) pressureload=md.mesh.segments(pos,:); @@ -106,2 +107,3 @@ disp(' no thermal boundary conditions created: no observed temperature found'); end +