Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 27645) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 27646) @@ -225,6 +225,5 @@ iomodel->FindConstant(&ismethod,"md.smb.ismethod"); if (ismethod == 1){ - // initial albedo: TODO: this is not needed. - //if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - albedo.\n"); + if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - albedo.\n"); iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo",SmbAlbedoInitEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo_snow",SmbAlbedoSnowInitEnum); @@ -239,4 +238,5 @@ // assign masking iomodel->FetchDataToInput(inputs,elements,"md.smb.mask",SmbMaskEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.qmr",SmbSemicQmrInitEnum); if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - done.\n"); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27645) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27646) @@ -31,6 +31,9 @@ IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out); -extern "C" void run_semic_transient_(IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in, - IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *tstic, +extern "C" void run_semic_transient_(int *nx, int *ntime, int *nloop, + IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, + IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in, + IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *qmr_in, + IssmDouble *tstic, IssmDouble *hcrit, IssmDouble *rcrit, IssmDouble *mask, IssmDouble *hice, IssmDouble *hsnow, @@ -38,6 +41,6 @@ int *alb_scheme, IssmDouble *alb_smax, IssmDouble *alb_smin, IssmDouble *albi, IssmDouble *albl, IssmDouble *Tamp, - IssmDouble *tmin, IssmDouble *tmax, IssmDouble *tmid, IssmDouble *mcrit, IssmDouble *wcrit, IssmDouble *tau_a, IssmDouble* tau_f, IssmDouble *afac, - IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *smbi_out, IssmDouble *smbs_out, IssmDouble *saccu_out, IssmDouble *smelt_out, IssmDouble *albedo_out, IssmDouble *albedo_snow_out, IssmDouble *hsnow_out, IssmDouble *hice_out); + IssmDouble *tmin, IssmDouble *tmax, IssmDouble *tmid, IssmDouble *mcrit, IssmDouble *wcrit, IssmDouble *tau_a, IssmDouble* tau_f, IssmDouble *afac, bool *verbose, + IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *smbi_out, IssmDouble *smbs_out, IssmDouble *saccu_out, IssmDouble *smelt_out, IssmDouble *albedo_out, IssmDouble *albedo_snow_out, IssmDouble *hsnow_out, IssmDouble *hice_out, IssmDouble *qmr_out); #endif // _HAVE_SEMIC_ @@ -4231,9 +4234,12 @@ void Element::SmbSemicTransient(){/*{{{*/ + bool isverbose=VerboseSmb(); + if(isverbose && this->Sid()==0){ + _printf0_("smb core: initialize.\n"); + } /*only compute SMB at the surface: */ if (!IsOnSurface()) return; const int NUM_VERTICES = this->GetNumberOfVertices(); - const int NLOOP = 10; // default internal iteration for SEMIC (Ruckamp et al. 2018). // daily forcing inputs @@ -4261,4 +4267,5 @@ IssmDouble* hice_in =xNew(NUM_VERTICES); IssmDouble* hsnow_in =xNew(NUM_VERTICES); + IssmDouble* qmr_in =xNew(NUM_VERTICES); // outputs @@ -4271,6 +4278,7 @@ IssmDouble* albedo_out =xNew(NUM_VERTICES); memset(albedo_out, 0., NUM_VERTICES*sizeof(IssmDouble)); IssmDouble* albedo_snow_out =xNew(NUM_VERTICES); memset(albedo_snow_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* hsnow_out =xNew(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble)); + IssmDouble* hsnow_out =xNew(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble)); IssmDouble* hice_out =xNew(NUM_VERTICES); memset(hice_out, 0., NUM_VERTICES*sizeof(IssmDouble)); + IssmDouble* qmr_out =xNew(NUM_VERTICES); memset(qmr_out, 0., NUM_VERTICES*sizeof(IssmDouble)); IssmDouble rho_water,rho_ice,desfac,desfacElev,rlaps,rdl; @@ -4325,5 +4333,5 @@ GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum); - if(VerboseSmb() && this->Sid()==0){ + if(isverbose && this->Sid()==0){ _printf0_("smb core: allocate inputs.\n"); _printf0_("smb core: time_yr : " << time_yr/yts <<"\n"); @@ -4345,9 +4353,10 @@ /*temporal Enum depending on time*/ - int enum_temp =TemperatureSEMICEnum; - int enum_hice =SmbHIceEnum; - int enum_hsnow=SmbHSnowEnum; - int enum_albedo=SmbAlbedoEnum; + int enum_temp =TemperatureSEMICEnum; + int enum_hice =SmbHIceEnum; + int enum_hsnow =SmbHSnowEnum; + int enum_albedo =SmbAlbedoEnum; int enum_albedo_snow=SmbAlbedoSnowEnum; + int enum_qmr =SmbSemicQmrEnum; if (tstart+dt == time) { /* Load inital value at first time step*/ @@ -4357,14 +4366,21 @@ enum_albedo=SmbAlbedoInitEnum; enum_albedo_snow=SmbAlbedoSnowInitEnum; + enum_qmr =SmbSemicQmrInitEnum; } + if(isverbose && this->Sid()==0)_printf0_("smb core: assign temp.\n"); Input* tsurf_input = this->GetInput(enum_temp); _assert_(tsurf_in); + if(isverbose && this->Sid()==0)_printf0_("smb core: assign mask.\n"); Input* mask_input = this->GetInput(SmbMaskEnum); _assert_(mask_input); + if(isverbose && this->Sid()==0)_printf0_("smb core: assign Tamp.\n"); Input* Tamp_input = this->GetInput(SmbTampEnum); _assert_(Tamp_input); + if(isverbose && this->Sid()==0)_printf0_("smb core: assign albedo.\n"); Input* albedo_input = this->GetInput(enum_albedo); _assert_(albedo_input); Input* albedo_snow_input = this->GetInput(enum_albedo_snow); _assert_(albedo_snow_input); Input* hice_input = this->GetInput(enum_hice); _assert_(hice_input); Input* hsnow_input = this->GetInput(enum_hsnow); _assert_(hsnow_input); - - if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: assign inputs.\n"); + if(isverbose && this->Sid()==0)_printf0_("smb core: assign qmr.\n"); + Input* qmr_input = this->GetInput(enum_qmr); _assert_(qmr_input); + + if(isverbose && this->Sid()==0)_printf0_("smb core: assign inputs done....\n"); for(int iv=0;ivGaussVertex(iv); @@ -4389,4 +4405,5 @@ hsnow_input->GetInputValue(&hsnow_in[iv],gauss); hice_input->GetInputValue(&hice_in[iv],gauss); + qmr_input->GetInputValue(&qmr_in[iv],gauss); /* Surface temperature correction */ @@ -4407,5 +4424,5 @@ dailydlradiation[iv]=dailydlradiation[iv]+rdl*st[iv]; } - if(VerboseSmb() && this->Sid()==0){ + if(isverbose && this->Sid()==0){ _printf0_("smb core: assign tsurf_in :" << tsurf_in[0] << "\n"); _printf0_("smb core: assign dailytemperature:" << dailytemperature[0] << "\n"); @@ -4413,22 +4430,28 @@ _printf0_("smb core: assign hice :" << hice_in[0] << "\n"); _printf0_("smb core: assign mask :" << mask_in[0] << "\n"); + _printf0_("smb core: assign Tamp :" << Tamp_in[0] << "\n"); _printf0_("smb core: assign albedo :" << albedo_in[0] << "\n"); _printf0_("smb core: assign albedo_snow :" << albedo_snow_in[0] << "\n"); _printf0_("smb core: assign albedo_scheme :" << alb_scheme << "\n"); - } - - if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: call semic module.\n"); + _printf0_("smb core: assign qmr :" << qmr_in[0] << "\n"); + } + + if(isverbose && this->Sid()==0)_printf0_("smb core: call run_semic_transient module.\n"); + /* call semic */ + int nx=NUM_VERTICES, ntime=1, nloop=1; + bool semic_verbose=false; //VerboseSmb(); + run_semic_transient_(&nx, &ntime, &nloop, + dailysnowfall, dailyrainfall, dailydsradiation, dailydlradiation, + dailywindspeed, dailypressure, dailyairdensity, dailyairhumidity, dailytemperature, tsurf_in, qmr_in, + &dt, + &hcrit, &rcrit, + mask_in, hice_in, hsnow_in, + albedo_in, albedo_snow_in, + &alb_scheme, &alb_smax, &alb_smin, &albi, &albl, + Tamp_in, + &tmin, &tmax, &tmid, &mcrit, &wcrit, &tau_a, &tau_f, &afac, &semic_verbose, + tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out,albedo_out, albedo_snow_out, hsnow_out, hice_out, qmr_out); + for (int iv = 0; iv ice @@ -4440,18 +4463,19 @@ } - if(VerboseSmb() && this->Sid()==0){ - _printf0_("smb core: tsurf_out[0] " << tsurf_out[0] << "\n"); - _printf0_("smb core: hice_out[0] " << hice_out[0] << "\n"); - _printf0_("smb core: hsnow_out[0] " << hsnow_out[0] << "\n"); - _printf0_("smb core: smb_ice[0] " << smbi_out[0]*yts << "\n"); + if(isverbose && this->Sid()==0){ + _printf0_("smb core: tsurf_out " << tsurf_out[0] << " " << tsurf_out[1] << " " << tsurf_out[2] << "\n"); + _printf0_("smb core: hice_out " << hice_out[0] << " " << hice_out[1] << " " << hice_out[2] << "\n"); + _printf0_("smb core: hsnow_out " << hsnow_out[0] << "\n"); + _printf0_("smb core: smb_ice " << smbi_out[0]*yts << "\n"); + _printf0_("smb core: smb_ice " << albedo_out[0] <<" "<ObjectEnum()){ case TriaEnum: - this->AddInput(TemperatureSEMICEnum,&tsurf_out[0],P1DGEnum); + this->AddInput(TemperatureSEMICEnum, &tsurf_out[0],P1DGEnum); // SMBout = SMB_ice + SMB_snow values. //this->AddInput(SmbMassBalanceTotalEnum,&smb_out[0],P1DGEnum); // water equivalent SMB ice to ice equivalent. - this->AddInput(SmbMassBalanceEnum,&smbi_out[0],P1DGEnum); + this->AddInput(SmbMassBalanceEnum, &smbi_out[0],P1DGEnum); this->AddInput(SmbMassBalanceSnowEnum,&smbs_out[0],P1DGEnum); //this->AddInput(SmbMassBalanceSnowEnum,&smbs_out[0],P1DGEnum); @@ -4459,9 +4483,10 @@ this->AddInput(SmbAccumulationEnum,&saccu_out[0],P1DGEnum); // smelt - this->AddInput(SmbMeltEnum,&smelt_out[0],P1DGEnum); - this->AddInput(SmbAlbedoEnum,&albedo_out[0],P1DGEnum); - this->AddInput(SmbAlbedoSnowEnum,&albedo_snow_out[0],P1DGEnum); - this->AddInput(SmbHSnowEnum,&hsnow_out[0],P1DGEnum); - this->AddInput(SmbHIceEnum,&hice_out[0],P1DGEnum); + this->AddInput(SmbMeltEnum, &smelt_out[0],P1DGEnum); + this->AddInput(SmbAlbedoEnum, &albedo_out[0],P1DGEnum); + this->AddInput(SmbAlbedoSnowEnum, &albedo_snow_out[0],P1DGEnum); + this->AddInput(SmbHSnowEnum, &hsnow_out[0],P1DGEnum); + this->AddInput(SmbHIceEnum, &hice_out[0],P1DGEnum); + this->AddInput(SmbSemicQmrEnum, &qmr_out[0],P1DGEnum); break; case PentaEnum: @@ -4498,4 +4523,5 @@ xDelete(hsnow_out); xDelete(hice_out); + xDelete(qmr_out); /*for inputs*/ @@ -4507,4 +4533,5 @@ xDelete(albedo_snow_in); xDelete(tsurf_in); + xDelete(qmr_in); /* for inputs:geometry */ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90 (revision 27645) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90 (revision 27646) @@ -1,4 +1,5 @@ -subroutine run_semic_transient(sf_in, rf_in, swd_in, lwd_in, wind_in, & - sp_in, rhoa_in, qq_in, tt_in, tsurf_in, tstic, & +subroutine run_semic_transient(nx, ntime, nloop, sf_in, rf_in, swd_in, lwd_in, wind_in, & + sp_in, rhoa_in, qq_in, tt_in, tsurf_in, qmr_in, & + tstic, & hcrit, rcrit, & mask, hice, hsnow, & @@ -6,10 +7,10 @@ alb_scheme, alb_smax, alb_smin, albi, albl, & Tamp, & - tmin, tmax, tmid, mcrit, w_crit, tau_a, tau_f, afac, & - tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, alb_out, alb_snow_out,hsnow_out,hice_out) + tmin, tmax, tmid, mcrit, w_crit, tau_a, tau_f, afac, verbose, & + tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, alb_out, & + alb_snow_out,hsnow_out,hice_out,qmr_out) !{{{ use utils use surface_physics - implicit none @@ -17,57 +18,66 @@ type(surface_physics_class) :: surface ! declare forcing class - type(forc_class) :: forc + !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, year - integer :: day=1 !< not used value. + integer :: i, k, n + integer :: nnx, nny + integer :: year=0 + integer :: day =1 !< not used value. + integer, intent(in) :: nx, ntime, nloop ! number of grid / number of times + logical, intent(in) :: verbose ! verbosity + logical :: debug=.false. ! forcing data - double precision, intent(in) :: sf_in ! snow fall rate [m/s] - double precision, intent(in) :: rf_in ! rain fall rate [m/s] - double precision, intent(in) :: swd_in ! downwelling shortwave radiation [W/m2] - double precision, intent(in) :: lwd_in ! downwelling longwave radiation [W/m2] - double precision, intent(in) :: wind_in! surface wind speed [m/s] - double precision, intent(in) :: sp_in ! surface pressure [Pa] - double precision, intent(in) :: rhoa_in! air density [kg/m3] - double precision, intent(in) :: qq_in ! air specific humidity [kg/kg] - double precision, intent(in) :: tt_in ! air temperature [K] + ! input argument format array size (nx, ntime)... + double precision, intent(in), dimension(nx):: sf_in ! snow fall rate [m/s] + double precision, intent(in), dimension(nx):: rf_in ! rain fall rate [m/s] + double precision, intent(in), dimension(nx):: swd_in ! downwelling shortwave radiation [W/m2] + double precision, intent(in), dimension(nx):: lwd_in ! downwelling longwave radiation [W/m2] + double precision, intent(in), dimension(nx):: wind_in ! surface wind speed [m/s] + double precision, intent(in), dimension(nx):: sp_in ! surface pressure [Pa] + double precision, intent(in), dimension(nx):: rhoa_in ! air density [kg/m3] + double precision, intent(in), dimension(nx):: qq_in ! air specific humidity [kg/kg] + double precision, intent(in), dimension(nx):: tt_in ! air temperature [K] ! input data double precision, intent(in) :: tstic ! time step from ISSM [sec]. - double precision, intent(in) :: tsurf_in ! input temperature [K] ! output data - double precision, intent(out) :: tsurf_out ! Ice surface Temperature [K] - double precision, intent(out) :: smb_out ! surface mass balance=(Accu-Melt) [m/s] - double precision, intent(out) :: smbi_out ! SMB ice [water equivalent m/s] - double precision, intent(out) :: smbs_out ! SMB snow [water equivalent m/s] - double precision, intent(out) :: saccu_out ! accumulation [m/s] - double precision, intent(out) :: smelt_out ! ablation [m/s] - double precision, intent(out) :: alb_out ! grid-averaged albedo [no unit] - double precision, intent(out) :: alb_snow_out - double precision :: qmr_out - !double precision, intent(out) :: qmr_out - double precision, intent(out) :: hice_out - double precision, intent(out) :: hsnow_out + ! Ice surface Temperature [K] + double precision, intent(out), dimension(nx):: tsurf_out + ! surface mass balance=(Accu-Melt) [m/s] + double precision, intent(out), dimension(nx):: smb_out + double precision, intent(out), dimension(nx):: smbi_out ! SMB ice [water equivalent m/s] + double precision, intent(out), dimension(nx):: smbs_out ! SMB snow [water equivalent m/s] + double precision, intent(out), dimension(nx):: saccu_out ! accumulation [m/s] + double precision, intent(out), dimension(nx):: smelt_out ! ablation [m/s] + double precision, intent(out), dimension(nx):: alb_out ! grid-averaged albedo [no unit] + double precision, intent(out), dimension(nx):: alb_snow_out + double precision, intent(out), dimension(nx):: hice_out + double precision, intent(out), dimension(nx):: hsnow_out + double precision, intent(out), dimension(nx):: qmr_out double precision :: total_time, start, finish ! set parameters - character (len=256) :: name ! not used(?) - character (len=256) :: boundary(30) ! not used(?) + !character (len=256) :: name ! not used(?) + !character (len=256) :: boundary(30) ! not used(?) !character (len=256), intent(in) :: alb_scheme !< name of albedo scheme integer, intent(in) :: alb_scheme - integer :: n_ksub - double precision :: ceff !< surface heat heat capacity of snow/ice [J/W m2] - double precision, intent(in) :: hsnow - double precision, intent(in) :: hice + !integer :: n_ksub + !double precision :: ceff !< surface heat heat capacity of snow/ice [J/W m2] + double precision, intent(in), dimension(nx):: albedo + double precision, intent(in), dimension(nx):: albedo_snow !< spatial.. + double precision, intent(in), dimension(nx):: hsnow + double precision, intent(in), dimension(nx):: hice + double precision, intent(in), dimension(nx):: tsurf_in !< input temperature [K] + double precision, intent(in), dimension(nx):: qmr_in + double precision, intent(in), dimension(nx):: mask + double precision, intent(in) :: albi double precision, intent(in) :: albl - double precision, intent(in) :: mask - double precision, intent(in) :: albedo - double precision, intent(in) :: albedo_snow double precision, intent(in) :: alb_smax double precision, intent(in) :: alb_smin @@ -75,9 +85,9 @@ double precision, intent(in) :: rcrit !< critical snow height fro which refreezing fraction is 50% [m] double precision, intent(in) :: Tamp - double precision :: csh - double precision :: clh + !double precision :: csh + !double precision :: clh double precision, intent(in) :: tmin double precision, intent(in) :: tmax - double precision :: tsticsub + !double precision :: tsticsub ! parameters for isba albedo scheme. double precision, intent(in) :: tau_a !< critical liquide water concent for "isba" albedo scheme [kg/m2] @@ -88,24 +98,17 @@ double precision, intent(in) :: tmid !< param for "alex" albedo parameterization [K] - nx = 1 - year = 0 + if (debug) then + print*,' ntime: ', ntime + print*,' nx : ', nx + end if ! set vector length surface%par%nx = nx - ! set input (forcing data) - !forc%sf(:,:) = sf_in ! snowfall flux [m/sec] - !forc%rf(:,:) = rf_in ! rainfall flux [m/sec] - !forc%swd(:,:) = swd_in ! short radiation - !forc%lwd(:,:) = lwd_in ! long radiation downward - !forc%wind(:,:) = wind_in ! wind speed. - !forc%sp(:,:) = sp_in ! surface pressure - !forc%rhoa(:,:) = rhoa_in ! air density [kg/m3] - !forc%qq(:,:) = qq_in ! air specific humidity [kg/kg] - !forc%tt(:,:) = tt_in ! 2m air temperature - ! FIXME should be user input !boundary = "" "" "" - !surface%par%tstic = 86400.0_dp !< time step [s] in one day. + if (debug) then + print*, "run_semic_transient: initialize parameters." + end if surface%par%tstic = tstic !< time step [s] surface%par%ceff= 2.0e6_dp !< surface heat capacity of snow/ice [J/K/m2] @@ -114,6 +117,6 @@ surface%par%alb_smax = alb_smax !0.79_dp !< max snow albedo surface%par%alb_smin = alb_smin !0.6_dp !< min snow albedo - surface%par%albi = albi !0.41_dp !< albedo for ice - surface%par%albl = albl !0.07_dp !< albedo for land + surface%par%albi = albi ! 0.41_dp !< albedo for ice + surface%par%albl = albl ! 0.07_dp !< albedo for land surface%par%tmin = tmin ! -999_dp surface%par%tmax = tmax ! 273.15_dp @@ -122,5 +125,5 @@ surface%par%amp = Tamp !3.0_dp !< amplitude of diurnal cycle [K] if (alb_scheme == 0) then - surface%par%alb_scheme="None" + surface%par%alb_scheme="none" else if (alb_scheme == 1) then surface%par%alb_scheme = "slater" @@ -129,4 +132,7 @@ else if (alb_scheme == 3) then surface%par%alb_scheme = "isba" + else + print*, "ERROR: current albedo scheme is not available." + call exit(1) end if surface%par%tau_a = tau_a !0.008_dp @@ -134,5 +140,5 @@ surface%par%w_crit = w_crit !15.0_dp ! default value surface%par%mcrit = mcrit !6.0e-8_dp - surface%par%n_ksub = 3.0_dp + surface%par%n_ksub = 3 ! sub ... ! snow albedo of alex surface%par%afac = afac @@ -146,24 +152,44 @@ ! initialise prognostic variables + if (debug) then + print*,"run_semic_transient: initialize variables." + end if ! these values will be updated through "surface_energy_and_mass_balance" function. - surface%now%mask(:) = mask ! 2.0_dp !loi_mask(:nx) - surface%now%hsnow(:) = hsnow ! initial snow height... - surface%now%hice(:) = hice ! initial ice height.. - surface%now%tsurf(:) = tsurf_in - surface%now%alb(:) = albedo !< initial albedo for energy balance. - surface%now%alb_snow(:) = albedo_snow !< initial albedo for ISBA albedo method. - surface%now%qmr_res(:) = 0.0_dp ! residual heat. - - ! initialize with zeros - surface%now%qmr(:) = 0.0_dp - - tsurf_out = 0.0_dp - smb_out = 0.0_dp - smbi_out = 0.0_dp - smbs_out = 0.0_dp - saccu_out = 0.0_dp - smelt_out = 0.0_dp - alb_out = 0.0_dp - alb_snow_out = 0.0_dp + surface%now%mask (:) = mask (:) ! 2.0_dp !loi_mask(:nx) + if (debug) then + print*,"run_semic_transient: initialize variables: mask" + end if + surface%now%hsnow (:) = hsnow (:) ! initial snow height... + surface%now%hice (:) = hice (:) ! initial ice height.. + surface%now%tsurf (:) = tsurf_in (:) !< initial ice surface temperature + surface%now%alb (:) = albedo (:) !< initial albedo for energy balance. + surface%now%alb_snow(:) = albedo_snow(:) !< initial albedo for ISBA albedo method. + if (debug) then + print*,"run_semic_transient: initialize variables. DONE." + end if + + if (debug) then + !print*, "====== global variable ==========" + !print*, "nloop :", nloop + !print*, "nx :", surface%par%nx + !print*, "====== parameters ======" + !print*, "csh :", surface%par%csh + !print*, "clh :", surface%par%clh + !print*, "albeo scheme :", surface%par%alb_scheme + !print*, "albeo ice :", surface%par%albi + !print*, "tstic :", surface%par%tstic + !print*, "tsticsub :", surface%par%tsticsub + !print*, "n_ksub :", surface%par%n_ksub + print*, "====== inputs =========" + print*, "hsnow :", hsnow + print*, "====== state variables ======" + print*, "hsnow :", surface%now%hsnow + print*, "hice :", surface%now%hice + print*, "albeo :", surface%now%alb + print*, "albeo snow :", surface%now%alb_snow + print*, "mask :", surface%now%mask + print*, "tsurf :", surface%now%tsurf + print*, "sf :", sf_in + end if ! define boundary conditions (not used, here!) @@ -171,42 +197,50 @@ !call print_boundary_opt(surface%bnd) - ! input with single value - surface%now%sf = sf_in - surface%now%rf = rf_in - surface%now%sp = sp_in - surface%now%lwd = lwd_in - surface%now%swd = swd_in - surface%now%wind = wind_in - surface%now%rhoa = rhoa_in - surface%now%t2m = tt_in - surface%now%qq = qq_in - - ! 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) - - tsurf_out=surface%now%tsurf(1) - ! melt - potential surface melt [m/s] - ! smb = SMB_ice + SMB_snow - ! smbi - SMB_ice (water equivalent m/sec) - ! smbs - SMB_snow (water equivalent m/sec) - smb_out =surface%now%smb(1) ! smb = smb_snow + smb_ice - smbi_out =surface%now%smb_ice(1) ! Csi (snow>ice) - melted_ice + refrezon_rain. - smbs_out =surface%now%smb_snow(1) ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow - saccu_out =surface%now%acc(1) ! acc = snowfall - sublimiation - refreezing - smelt_out =surface%now%melt(1) ! potential surface melt = melt_ice + melt_snow - alb_out =surface%now%alb(1) - alb_snow_out =surface%now%alb_snow(1) - qmr_out =surface%now%qmr(1) - hsnow_out =surface%now%hsnow(1) - hice_out =surface%now%hice(1) + do k =1,nloop + do i =1,ntime + if (debug) then + print*,"run_semic_transient: forcing data: ntime = ", i + end if + surface%now%sf = sf_in !(:,i) + surface%now%rf = rf_in !(:,i) + surface%now%sp = sp_in !(:,i) + surface%now%lwd = lwd_in !(:,i) + surface%now%swd = swd_in !(:,i) + surface%now%wind = wind_in!(:,i) + surface%now%rhoa = rhoa_in!(:,i) + surface%now%t2m = tt_in !(:,i) + surface%now%qq = qq_in !(:,i) + ! qmr_res is used to "energy_balance" in semic. + surface%now%qmr_res = qmr_in + + ! calculate prognostic and diagnsotic variables + call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,i,year) + + if (debug) then + print*,"done..." + end if + if (k == nloop) then + tsurf_out = surface%now%tsurf + ! melt - potential surface melt [m/s] + ! smb = SMB_ice + SMB_snow + ! smbi - SMB_ice (water equivalent m/sec) + ! smbs - SMB_snow (water equivalent m/sec) + smb_out =surface%now%smb ! smb = smb_snow + smb_ice + smbi_out =surface%now%smb_ice ! Csi (snow>ice) - melted_ice + refrezon_rain. + smbs_out =surface%now%smb_snow ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow + saccu_out =surface%now%acc ! acc = snowfall - sublimiation - refreezing + smelt_out =surface%now%melt ! potential surface melt = melt_ice + melt_snow + alb_out =surface%now%alb + alb_snow_out =surface%now%alb_snow + hsnow_out =surface%now%hsnow + hice_out =surface%now%hice + qmr_out =surface%now%qmr_res + end if + 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_transient +end subroutine run_semic_transient ! }}} Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27645) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27646) @@ -1118,4 +1118,6 @@ syn keyword cConstant SmbS0pEnum syn keyword cConstant SmbS0tEnum +syn keyword cConstant SmbSemicQmrEnum +syn keyword cConstant SmbSemicQmrInitEnum syn keyword cConstant SmbSizeiniEnum syn keyword cConstant SmbSmbCorrEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27645) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27646) @@ -1115,4 +1115,6 @@ SmbS0pEnum, SmbS0tEnum, + SmbSemicQmrEnum, + SmbSemicQmrInitEnum, SmbSizeiniEnum, SmbSmbCorrEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27645) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27646) @@ -1120,4 +1120,6 @@ case SmbS0pEnum : return "SmbS0p"; case SmbS0tEnum : return "SmbS0t"; + case SmbSemicQmrEnum : return "SmbSemicQmr"; + case SmbSemicQmrInitEnum : return "SmbSemicQmrInit"; case SmbSizeiniEnum : return "SmbSizeini"; case SmbSmbCorrEnum : return "SmbSmbCorr"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27645) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27646) @@ -1111,4 +1111,6 @@ syn keyword juliaConstC SmbS0pEnum syn keyword juliaConstC SmbS0tEnum +syn keyword juliaConstC SmbSemicQmrEnum +syn keyword juliaConstC SmbSemicQmrInitEnum syn keyword juliaConstC SmbSizeiniEnum syn keyword juliaConstC SmbSmbCorrEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27645) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27646) @@ -1147,4 +1147,6 @@ else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum; + else if (strcmp(name,"SmbSemicQmr")==0) return SmbSemicQmrEnum; + else if (strcmp(name,"SmbSemicQmrInit")==0) return SmbSemicQmrInitEnum; else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"VyBase")==0) return VyBaseEnum; - else if (strcmp(name,"VyDebris")==0) return VyDebrisEnum; - else if (strcmp(name,"Vy")==0) return VyEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"VyMesh")==0) return VyMeshEnum; + if (strcmp(name,"VyDebris")==0) return VyDebrisEnum; + else if (strcmp(name,"Vy")==0) return VyEnum; + else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; else if (strcmp(name,"VyObs")==0) return VyObsEnum; else if (strcmp(name,"VyShear")==0) return VyShearEnum; @@ -1365,10 +1367,10 @@ else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; - else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; - else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; else stage=12; } if(stage==12){ - if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; + if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; + else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; + else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum; else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum; @@ -1488,10 +1490,10 @@ else if (strcmp(name,"Fset")==0) return FsetEnum; else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum; - else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum; - else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; + if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum; + else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; + else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; @@ -1611,10 +1613,10 @@ else if (strcmp(name,"MinVel")==0) return MinVelEnum; else if (strcmp(name,"MinVx")==0) return MinVxEnum; - else if (strcmp(name,"MinVy")==0) return MinVyEnum; - else if (strcmp(name,"MinVz")==0) return MinVzEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; + if (strcmp(name,"MinVy")==0) return MinVyEnum; + else if (strcmp(name,"MinVz")==0) return MinVzEnum; + else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; @@ -1734,10 +1736,10 @@ else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else if (strcmp(name,"Tetra")==0) return TetraEnum; - else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; - else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; else stage=15; } if(stage==15){ - if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; + if (strcmp(name,"TetraInput")==0) return TetraInputEnum; + else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; + else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;