Index: /issm/trunk-jpl/externalpackages/semic/install.sh =================================================================== --- /issm/trunk-jpl/externalpackages/semic/install.sh (revision 23539) +++ /issm/trunk-jpl/externalpackages/semic/install.sh (revision 23540) @@ -3,9 +3,9 @@ #Some cleanup -#rm -rf install src -#mkdir install +rm -rf install src +mkdir install #Download latest version -#git clone https://github.com/mkrapp/semic.git src +git clone https://github.com/mkrapp/semic.git src if which ifort >/dev/null; then Index: /issm/trunk-jpl/m4/issm_options.m4 =================================================================== --- /issm/trunk-jpl/m4/issm_options.m4 (revision 23539) +++ /issm/trunk-jpl/m4/issm_options.m4 (revision 23540) @@ -1697,4 +1697,5 @@ AC_SUBST([SEMICLIB]) fi + AM_CONDITIONAL([SEMIC],[test x$HAVE_SEMIC = xyes]) dnl }}} dnl spai{{{ Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 23539) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 23540) @@ -539,5 +539,11 @@ endif #}}} - +#SEMIC sources {{{ +if SEMIC +if FORTRAN +issm_sources += ./modules/SurfaceMassBalancex/run_semic.f90 +endif +endif +#}}} #Wrapper sources #Kml sources {{{ @@ -645,5 +651,5 @@ if !WINDOWS if !STANDALONE_LIBRARIES -libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(M1QN3LIB) $(SEMICLIB)$(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB) $(ADOLCLIB) $(AMPILIB) $(ADJOINTMPILIB) $(METEOIOLIB) $(SNOWPACKLIB) +libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(M1QN3LIB) $(SEMICLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB) $(ADOLCLIB) $(AMPILIB) $(ADJOINTMPILIB) $(METEOIOLIB) $(SNOWPACKLIB) if FORTRAN libISSMCore_la_LIBADD += $(FLIBS) $(FORTRANLIB) Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23539) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23540) @@ -158,4 +158,17 @@ iomodel->FetchDataToInput(elements,"md.smb.runoffgrad",SmbRunoffgradEnum); break; + case SMBsemicEnum: + iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum); + iomodel->FetchDataToInput(elements,"md.smb.s0gcm",SmbS0gcmEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailypressure",SmbDailypressureEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum); + iomodel->FetchDataToInput(elements,"md.smb.dailytemperature",SmbDailytemperatureEnum); + break; default: _error_("Surface mass balance model "<DeleteData(temp,"md.smb.runoffref"); break; + case SMBsemicEnum: + /*Nothing to add to parameters*/ + break; default: _error_("Surface mass balance model "<GetNumberOfVertices(); + + IssmDouble* s=xNew(numvertices); + IssmDouble* s0gcm=xNew(numvertices); + IssmDouble* st=xNew(numvertices); + + // daily forcing inputs + IssmDouble* dailyrainfall=xNew(365*numvertices); + IssmDouble* dailysnowfall=xNew(365*numvertices); + IssmDouble* dailydlradiation=xNew(365*numvertices); + IssmDouble* dailydsradiation=xNew(365*numvertices); + IssmDouble* dailywindspeed=xNew(365*numvertices); + IssmDouble* dailypressure=xNew(365*numvertices); + IssmDouble* dailyairdensity=xNew(365*numvertices); + IssmDouble* dailyairhumidity=xNew(365*numvertices); + IssmDouble* dailytemperature=xNew(365*numvertices); + // daily outputs + IssmDouble* tsurf_out=xNew(numvertices); memset(tsurf_out, 0., numvertices*sizeof(IssmDouble)); + IssmDouble* smb_out=xNew(numvertices); memset(smb_out, 0., numvertices*sizeof(IssmDouble)); + IssmDouble* saccu_out=xNew(numvertices); memset(saccu_out, 0., numvertices*sizeof(IssmDouble)); + IssmDouble* smelt_out=xNew(numvertices); memset(smelt_out, 0., numvertices*sizeof(IssmDouble)); + + IssmDouble rho_water,rho_ice,desfac,rlaps,rdl; + IssmDouble time,yts,time_yr; + + /* Get time: */ + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /*Get material parameters :*/ + rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + desfac=this->matpar->GetMaterialParameter(SmbDesfacEnum); + rlaps=this->matpar->GetMaterialParameter(SmbRlapsEnum); + rdl=this->matpar->GetMaterialParameter(SmbRdlEnum); + + /* Retrieve inputs: */ + Input* dailysnowfall_input=this->GetInput(SmbDailysnowfallEnum); _assert_(dailysnowfall_input); + Input* dailyrainfall_input=this->GetInput(SmbDailyrainfallEnum); _assert_(dailyrainfall_input); + Input* dailydlradiation_input=this->GetInput(SmbDailydlradiationEnum); _assert_(dailydlradiation_input); + Input* dailydsradiation_input=this->GetInput(SmbDailydsradiationEnum); _assert_(dailydsradiation_input); + Input* dailywindspeed_input=this->GetInput(SmbDailywindspeedEnum); _assert_(dailywindspeed_input); + Input* dailypressure_input=this->GetInput(SmbDailypressureEnum); _assert_(dailypressure_input); + Input* dailyairdensity_input=this->GetInput(SmbDailyairdensityEnum); _assert_(dailyairdensity_input); + Input* dailyairhumidity_input=this->GetInput(SmbDailyairhumidityEnum); _assert_(dailyairhumidity_input); + Input* dailytemperature_input=this->GetInput(SmbDailytemperatureEnum); _assert_(dailytemperature_input); + + /* Recover info at the vertices: */ + GetInputListOnVertices(&s[0],SurfaceEnum); + GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum); + + /* loop over vertices and days */ //FIXME account for leap years (365 -> 366) + Gauss* gauss=this->NewGauss(); + for (int iday = 0; iday < 365; iday++){ + for(int iv=0;ivGaussVertex(iv); + /* get forcing */ + dailyrainfall_input->GetInputValue(&dailyrainfall[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailysnowfall_input->GetInputValue(&dailysnowfall[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailydlradiation_input->GetInputValue(&dailydlradiation[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailydsradiation_input->GetInputValue(&dailydsradiation[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailywindspeed_input->GetInputValue(&dailywindspeed[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailypressure_input->GetInputValue(&dailypressure[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailyairdensity_input->GetInputValue(&dailyairdensity[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailyairhumidity_input->GetInputValue(&dailyairhumidity[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + dailytemperature_input->GetInputValue(&dailytemperature[iv*365+iday],gauss,time_yr+(iday+1)/365.*yts); + + /* Surface temperature correction */ + st[iv]=(s[iv]-s0gcm[iv])/1000.; + dailytemperature[iv*365+iday]=dailytemperature[iv*365+iday]-rlaps *st[iv]; + + /* Precipitation correction (Vizcaino et al. 2010) */ + if (s0gcm[iv] < 2000.0) { + dailysnowfall[iv*365+iday] = dailysnowfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-2000.0)); + dailyrainfall[iv*365+iday] = dailyrainfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-2000.0)); + }else{ + dailysnowfall[iv*365+iday] = dailysnowfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-s0gcm[iv])); + dailyrainfall[iv*365+iday] = dailyrainfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-s0gcm[iv])); + } + + /* downward longwave radiation correction (Marty et al. 2002) */ + st[iv]=(s[iv]-s0gcm[iv])/1000.; + dailydlradiation[iv*365+iday]=dailydlradiation[iv*365+iday]+rdl*st[iv]; + } + } + + for (int iv = 0; ivObjectEnum()){ + case TriaEnum: + this->inputs->AddInput(new TriaInput(TemperatureSEMICEnum,&tsurf_out[0],P1Enum)); // TODO add TemperatureSEMICEnum to EnumDefinitions + this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb_out[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbAccumulationEnum,&saccu_out[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbMeltEnum,&smelt_out[0],P1Enum)); + break; + case PentaEnum: + // TODO + break; + case TetraEnum: + // TODO + break; + default: _error_("Not implemented yet"); + } + + /*clean-up*/ + delete gauss; + xDelete(dailysnowfall); + xDelete(dailyrainfall); + xDelete(dailydlradiation); + xDelete(dailydsradiation); + xDelete(dailywindspeed); + xDelete(dailypressure); + xDelete(dailyairdensity); + xDelete(dailyairhumidity); + xDelete(dailypressure); + xDelete(dailytemperature); + xDelete(smb_out); + xDelete(smelt_out); + xDelete(saccu_out); + xDelete(tsurf_out); + xDelete(s); + xDelete(st); + xDelete(s0gcm); +} +/*}}}*/ +#endif // _HAVE_SEMIC_ int Element::Sid(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23539) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23540) @@ -154,4 +154,5 @@ void ResultToVector(Vector* vector,int output_enum); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); + void SmbSemic(); int Sid(); void SmbGemb(); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23539) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23540) @@ -37,4 +37,5 @@ rlaps = 0; rlapslgm = 0; + rdl = 0; dpermil = 0; rheology_law = 0; @@ -133,4 +134,9 @@ case SMBgradientscomponentsEnum: /*Nothing to add*/ + break; + case SMBsemicEnum: + iomodel->FindConstant(&this->desfac,"md.smb.desfac"); + iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); + iomodel->FindConstant(&this->rdl,"md.smb.rdl"); break; default: @@ -520,4 +526,7 @@ this->rlapslgm=constant; break; + case SmbRdlEnum: + this->rdl=constant; + break; case SmbDpermilEnum: this->dpermil=constant; @@ -639,4 +648,5 @@ case SmbRlapsEnum: return this->rlaps; case SmbRlapslgmEnum: return this->rlapslgm; + case SmbRdlEnum: return this->rdl; case SmbDpermilEnum: return this->dpermil; case MaterialsLithosphereShearModulusEnum: return this->lithosphere_shear_modulus; Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23539) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23540) @@ -438,2 +438,14 @@ }/*}}}*/ +#ifdef _HAVE_SEMIC_ +void SmbSemicx(FemModel* femmodel){/*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SmbSemic(); + } + +}/*}}}*/ +#else +void SmbSemicx(FemModel* femmodel){_error_("SEMIC not installed");} +#endif //_HAVE_SEMIC_ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23539) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23540) @@ -21,4 +21,6 @@ void SmbMeltComponentsx(FemModel* femmodel); void SmbGradientsComponentsx(FemModel* femmodel); +/* SEMIC: */ +void SmbSemicx(FemModel* femmodel); /*GEMB: */ void Gembx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic.f90 (revision 23540) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic.f90 (revision 23540) @@ -0,0 +1,185 @@ +subroutine run_semic(sf_in, rf_in, swd_in, lwd_in, wind_in, & + sp_in, rhoa_in, qq_in, tt_in, tsurf_out, & + smb_out, saccu_out, smelt_out) + + use utils + use surface_physics + + implicit none + + ! declare surface physics class + type(surface_physics_class) :: surface + ! declare forcing class + type(forc_class) :: forc + ! declare validation class + !type(vali_class) :: vali ! validation not needed here + + integer, parameter:: dp=kind(0.d0) !< define precision (machine specific) + integer :: i, k, nx, nloop, ntime, year, day + + ! forcing data + double precision, intent(in),dimension(1,365) :: sf_in ! snow fall rate [m/s] + double precision, intent(in),dimension(1,365) :: rf_in ! rain fall rate [m/s] + double precision, intent(in),dimension(1,365) :: swd_in ! downwelling shortwave radiation [W/m2] + double precision, intent(in),dimension(1,365) :: lwd_in ! downwelling longwave radiation [W/m2] + double precision, intent(in),dimension(1,365) :: wind_in! surface wind speed [m/s] + double precision, intent(in),dimension(1,365) :: sp_in ! surface pressure [Pa] + double precision, intent(in),dimension(1,365) :: rhoa_in! air density [kg/m3] + double precision, intent(in),dimension(1,365) :: qq_in ! air specific humidity [kg/kg] + double precision, intent(in),dimension(1,365) :: tt_in ! air temperature [K] + + ! output data + double precision :: tsurf_out ! Ice surface Temperature [K] + double precision :: smb_out ! surface mass balance=(Accu-Melt) [m/s] + double precision :: saccu_out ! accumulation [m/s] + double precision :: smelt_out ! ablation [m/s] + + double precision :: total_time, start, finish + + ! set parameters + character (len=256) :: name ! not used(?) + character (len=256) :: boundary(30) ! not used(?) + character (len=256) :: alb_scheme + integer :: n_ksub + double precision :: ceff + double precision :: albi + double precision :: albl + double precision :: alb_smax + double precision :: alb_smin + double precision :: hcrit + double precision :: rcrit + double precision :: amp + double precision :: csh + double precision :: clh + double precision :: tmin + double precision :: tmax + double precision :: tstic + double precision :: tsticsub + double precision :: tau_a + double precision :: tau_f + double precision :: w_crit + double precision :: mcrit + double precision :: afac + double precision :: tmid + + nloop = 10 + nx = 1 + ntime = 365 + year = 0 + + ! set vector length + surface%par%nx = nx + + ! set input (forcing data) + allocate(forc%sf(nx,ntime)) + allocate(forc%rf(nx,ntime)) + allocate(forc%swd(nx,ntime)) + allocate(forc%lwd(nx,ntime)) + allocate(forc%wind(nx,ntime)) + allocate(forc%sp(nx,ntime)) + allocate(forc%rhoa(nx,ntime)) + allocate(forc%tt(nx,ntime)) + allocate(forc%qq(nx,ntime)) + +! write(*,*) sf_in + forc%sf(:,:) = sf_in + forc%rf(:,:) = rf_in + forc%swd(:,:) = swd_in + forc%lwd(:,:) = lwd_in + forc%wind(:,:) = wind_in + forc%sp(:,:) = sp_in + forc%rhoa(:,:) = rhoa_in + forc%qq(:,:) = qq_in + forc%tt(:,:) = tt_in + + ! FIXME should be user input + !boundary = "" "" "" + surface%par%tstic = 86400.0_dp + surface%par%ceff= 2.0e6_dp + surface%par%csh = 2.0e-3_dp + surface%par%clh = 5.0e-4_dp + surface%par%alb_smax = 0.79_dp + surface%par%alb_smin = 0.6_dp + surface%par%albi = 0.41_dp + surface%par%albl = 0.07_dp + surface%par%tmin = -999_dp + surface%par%tmax = 273.15_dp + surface%par%hcrit = 0.028_dp + surface%par%rcrit = 0.85_dp + surface%par%amp = 3.0_dp + surface%par%alb_scheme = "None" + surface%par%tau_a = 0.008_dp + surface%par%tau_f = 0.24_dp + surface%par%w_crit = 15.0_dp + surface%par%mcrit = 6.0e-8_dp + surface%par%n_ksub = 3.0_dp + + ! initialize sub-daily time step tsticsub + surface%par%tsticsub = surface%par%tstic / dble(surface%par%n_ksub) + + ! allocate necessary arrays for surface_physics module + call surface_alloc(surface%now,surface%par%nx) + + ! initialise prognostic variables + surface%now%mask(:) = 2.0_dp !loi_mask(:nx) + surface%now%hsnow(:) = 1.0_dp + surface%now%hice(:) = 0.0_dp + surface%now%alb(:) = 0.8_dp + surface%now%tsurf(:) = 260.0_dp + surface%now%alb_snow(:) = 0.8_dp + !surface%now%acc(:) = 0.0_dp + !surface%now%smb(:) = 0.0_dp + !surface%now%melt(:) = 0.0_dp + surface%now%qmr_res(:) = 0.0_dp + + tsurf_out = 0.0_dp + smb_out = 0.0_dp + saccu_out = 0.0_dp + smelt_out = 0.0_dp + + ! define boundary conditions (not used, here!) + call surface_boundary_define(surface%bnd,surface%par%boundary) + !call print_boundary_opt(surface%bnd) + + do k=1,nloop ! re-iterate 'nloop' times + + day = 1 + + do i=1,ntime ! loop over one year + + + ! read input for i-th day of year + surface%now%sf = forc%sf(1,day) + surface%now%rf = forc%rf(1,day) + surface%now%sp = forc%sp(1,day) + surface%now%lwd = forc%lwd(1,day) + surface%now%swd = forc%swd(1,day) + surface%now%wind = forc%wind(1,day) + surface%now%rhoa = forc%rhoa(1,day) + surface%now%t2m = forc%tt(1,day) + surface%now%qq = forc%qq(1,day) + + ! calculate prognostic and diagnsotic variables + call cpu_time(start) + call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,day,-1) + call cpu_time(finish) + total_time = total_time + (finish - start) + + if (k==nloop) then + tsurf_out=tsurf_out+surface%now%tsurf(1)*1.0_dp/365.0_dp + smb_out=smb_out+surface%now%smb(1)*1.0_dp/365.0_dp + saccu_out=saccu_out+surface%now%alb(1)*1.0_dp/365.0_dp + smelt_out=smelt_out+surface%now%melt(1)*1.0_dp/365.0_dp + endif + day = day + 1 + + end do + + end do + + ! de-allocate surface_physics arrays + call surface_dealloc(surface%now) + + !write(*,*) 'total time for surface_physics:', nloop, total_time + +end subroutine run_semic Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23539) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23540) @@ -560,4 +560,13 @@ SmbBPosEnum, SmbCEnum, + SmbDailysnowfallEnum, + SmbDailyrainfallEnum, + SmbDailydsradiationEnum, + SmbDailydlradiationEnum, + SmbDailywindspeedEnum, + SmbDailypressureEnum, + SmbDailyairdensityEnum, + SmbDailyairhumidityEnum, + SmbDailytemperatureEnum, SmbDEnum, SmbDiniEnum, @@ -597,4 +606,5 @@ SmbS0pEnum, SmbS0tEnum, + SmbS0gcmEnum, SmbSizeiniEnum, SmbSmbrefEnum, @@ -1093,4 +1103,5 @@ SigmaVMEnum, SmbAnalysisEnum, + SMBsemicEnum, SMBcomponentsEnum, SMBd18opddEnum, @@ -1110,4 +1121,5 @@ SMBpddSicopolisEnum, SMBgradientscomponentsEnum, + SmbRdlEnum, SmbRlapsEnum, SmbRlapslgmEnum, @@ -1141,4 +1153,5 @@ TaylorHoodEnum, TemperaturePDDEnum, + TemperatureSEMICEnum, TetraEnum, TetraInputEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23539) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23540) @@ -566,4 +566,13 @@ case SmbBPosEnum : return "SmbBPos"; case SmbCEnum : return "SmbC"; + case SmbDailysnowfallEnum : return "SmbDailysnowfall"; + case SmbDailyrainfallEnum : return "SmbDailyrainfall"; + case SmbDailydsradiationEnum : return "SmbDailydsradiation"; + case SmbDailydlradiationEnum : return "SmbDailydlradiation"; + case SmbDailywindspeedEnum : return "SmbDailywindspeed"; + case SmbDailypressureEnum : return "SmbDailypressure"; + case SmbDailyairdensityEnum : return "SmbDailyairdensity"; + case SmbDailyairhumidityEnum : return "SmbDailyairhumidity"; + case SmbDailytemperatureEnum : return "SmbDailytemperature"; case SmbDEnum : return "SmbD"; case SmbDiniEnum : return "SmbDini"; @@ -603,4 +612,5 @@ case SmbS0pEnum : return "SmbS0p"; case SmbS0tEnum : return "SmbS0t"; + case SmbS0gcmEnum : return "SmbS0gcm"; case SmbSizeiniEnum : return "SmbSizeini"; case SmbSmbrefEnum : return "SmbSmbref"; @@ -1097,4 +1107,5 @@ case SigmaVMEnum : return "SigmaVM"; case SmbAnalysisEnum : return "SmbAnalysis"; + case SMBsemicEnum : return "SMBsemic"; case SMBcomponentsEnum : return "SMBcomponents"; case SMBd18opddEnum : return "SMBd18opdd"; @@ -1114,4 +1125,5 @@ case SMBpddSicopolisEnum : return "SMBpddSicopolis"; case SMBgradientscomponentsEnum : return "SMBgradientscomponents"; + case SmbRdlEnum : return "SmbRdl"; case SmbRlapsEnum : return "SmbRlaps"; case SmbRlapslgmEnum : return "SmbRlapslgm"; @@ -1145,4 +1157,5 @@ case TaylorHoodEnum : return "TaylorHood"; case TemperaturePDDEnum : return "TemperaturePDD"; + case TemperatureSEMICEnum : return "TemperatureSEMIC"; case TetraEnum : return "Tetra"; case TetraInputEnum : return "TetraInput"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23539) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23540) @@ -578,4 +578,13 @@ else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum; else if (strcmp(name,"SmbC")==0) return SmbCEnum; + else if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum; + else if (strcmp(name,"SmbDailyrainfall")==0) return SmbDailyrainfallEnum; + else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum; + else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; + else if (strcmp(name,"SmbDailywindspeed")==0) return SmbDailywindspeedEnum; + else if (strcmp(name,"SmbDailypressure")==0) return SmbDailypressureEnum; + else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; + else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum; + else if (strcmp(name,"SmbDailytemperature")==0) return SmbDailytemperatureEnum; else if (strcmp(name,"SmbD")==0) return SmbDEnum; else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum; @@ -615,9 +624,13 @@ else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum; + else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum; else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; - else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; @@ -629,8 +642,5 @@ else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"SmbV")==0) return SmbVEnum; + else if (strcmp(name,"SmbV")==0) return SmbVEnum; else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; else if (strcmp(name,"SmbW")==0) return SmbWEnum; @@ -742,5 +752,8 @@ else if (strcmp(name,"Contour")==0) return ContourEnum; else if (strcmp(name,"Contours")==0) return ContoursEnum; - else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"ControlInputValues")==0) return ControlInputValuesEnum; else if (strcmp(name,"ControlInputMins")==0) return ControlInputMinsEnum; @@ -752,8 +765,5 @@ else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"DataSet")==0) return DataSetEnum; + else if (strcmp(name,"DataSet")==0) return DataSetEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; @@ -865,5 +875,8 @@ else if (strcmp(name,"Intersect")==0) return IntersectEnum; else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; - else if (strcmp(name,"IntInput")==0) return IntInputEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; @@ -875,8 +888,5 @@ else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; + else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; @@ -988,5 +998,8 @@ else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; - else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; @@ -998,8 +1011,5 @@ else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; + else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; @@ -1111,5 +1121,8 @@ else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; else if (strcmp(name,"SedimentHeadStacked")==0) return SedimentHeadStackedEnum; - else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"Seg")==0) return SegEnum; else if (strcmp(name,"SegInput")==0) return SegInputEnum; @@ -1121,8 +1134,6 @@ else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; + else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; + else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; @@ -1141,4 +1152,5 @@ else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; @@ -1172,4 +1184,5 @@ else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; + else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum; else if (strcmp(name,"Tetra")==0) return TetraEnum; else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23539) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23540) @@ -185,4 +185,5 @@ case 10: return SMBpddSicopolisEnum; case 11: return SMBgradientscomponentsEnum; + case 12: return SMBsemicEnum; default: _error_("Marshalled SMB code \""<=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1); + md = checkfield(md,'fieldname','smb.rdl','>=',0,'numel',1); + md = checkfield(md,'fieldname','smb.dailysnowfall','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailyrainfall','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailydsradiation','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailydlradiation','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailywindspeed','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailypressure','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailyairdensity','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailyairhumidity','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.dailytemperature','timeseries',1,'NaN',1,'Inf',1); + end + md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); + + end % }}} + function disp(self) % {{{ + disp(sprintf(' surface forcings parameters:')); + + disp(sprintf(' Interface for coupling GCM data to the energy balance model SEMIC (Krapp et al (2017) https://doi.org/10.5194/tc-11-1519-2017).')); + disp(sprintf(' The implemented coupling uses daily mean GCM input to calculate yearly mean smb, accumulation, ablation, and surface temperature.')); + disp(sprintf(' smb and temperatures are updated every year')); + disp(sprintf('\n SEMIC parameters:')); + fielddisplay(self,'dailysnowfall','daily surface dailysnowfall [m/s]'); + fielddisplay(self,'dailyrainfall','daily surface dailyrainfall [m/s]'); + fielddisplay(self,'dailydsradiation','daily downwelling shortwave radiation [W/m2]'); + fielddisplay(self,'dailydlradiation','daily downwelling longwave radiation [W/m2]'); + fielddisplay(self,'dailywindspeed','daily surface wind speed [m/s]'); + fielddisplay(self,'dailypressure','daily surface pressure [Pa]'); + fielddisplay(self,'dailyairdensity','daily air density [kg/m3]'); + fielddisplay(self,'dailyairhumidity','daily air specific humidity [kg/kg]'); + fielddisplay(self,'dailytemperature','daily surface air temperature [K]'); + fielddisplay(self,'rlaps','present day lapse rate (default is 7.4 [degree/km]; )Erokhina et al. 2017)'); + fielddisplay(self,'desfac','desertification elevation factor (default is -log(2.0)/1000 [1/km]; Vizcaino et al. 2010)'); + fielddisplay(self,'rdl','longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)'); + fielddisplay(self,'s0gcm','GCM reference elevation; (default is 0) [m]'); + fielddisplay(self,'requested_outputs','additional outputs requested'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + WriteData(fid,prefix,'name','md.smb.model','data',12,'format','Integer'); + + WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0gcm','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','rdl','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailysnowfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyrainfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydsradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydlradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailywindspeed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailypressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairdensity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairhumidity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + + %process requested outputs + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray'); + + end % }}} + end +end