Changeset 27646


Ignore:
Timestamp:
03/17/23 16:55:57 (2 years ago)
Author:
inwoo
Message:

CHG: fix SMBsemic - consider QMR(heat from melting or refreezing) in semic..

Location:
issm/trunk-jpl/src/c
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r27510 r27646  
    225225                        iomodel->FindConstant(&ismethod,"md.smb.ismethod");
    226226                        if (ismethod == 1){
    227                                 // initial albedo: TODO: this is not needed.
    228                                 //if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - albedo.\n");
     227                                if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - albedo.\n");
    229228                                iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo",SmbAlbedoInitEnum);
    230229                                iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo_snow",SmbAlbedoSnowInitEnum);
    … …  
    239238                                // assign masking
    240239                                iomodel->FetchDataToInput(inputs,elements,"md.smb.mask",SmbMaskEnum);
     240                                iomodel->FetchDataToInput(inputs,elements,"md.smb.qmr",SmbSemicQmrInitEnum);
    241241                                if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - done.\n");
    242242                        }

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r27586 r27646  
    3131                        IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out);
    3232
    33 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,
    34                         IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *tstic,
     33extern "C" void run_semic_transient_(int *nx, int *ntime, int *nloop,
     34                        IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in,
     35                        IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in,
     36                        IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *qmr_in,
     37                        IssmDouble *tstic,
    3538                        IssmDouble *hcrit, IssmDouble *rcrit,
    3639                        IssmDouble *mask, IssmDouble *hice, IssmDouble *hsnow,
    … …  
    3841                        int *alb_scheme, IssmDouble *alb_smax, IssmDouble *alb_smin, IssmDouble *albi, IssmDouble *albl,
    3942                        IssmDouble *Tamp,
    40                         IssmDouble *tmin, IssmDouble *tmax, IssmDouble *tmid, IssmDouble *mcrit, IssmDouble *wcrit, IssmDouble *tau_a, IssmDouble* tau_f, IssmDouble *afac,
    41                         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);
     43                        IssmDouble *tmin, IssmDouble *tmax, IssmDouble *tmid, IssmDouble *mcrit, IssmDouble *wcrit, IssmDouble *tau_a, IssmDouble* tau_f, IssmDouble *afac, bool *verbose,
     44                        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);
    4245#endif
    4346// _HAVE_SEMIC_
    … …  
    42314234void       Element::SmbSemicTransient(){/*{{{*/
    42324235
     4236        bool isverbose=VerboseSmb();
     4237        if(isverbose && this->Sid()==0){
     4238                _printf0_("smb core: initialize.\n");
     4239        }
    42334240        /*only compute SMB at the surface: */
    42344241        if (!IsOnSurface()) return;
    42354242
    42364243        const int NUM_VERTICES                                  = this->GetNumberOfVertices();
    4237         const int NLOOP = 10; // default internal iteration for SEMIC (Ruckamp et al. 2018).
    42384244
    42394245        // daily forcing inputs
    … …  
    42614267        IssmDouble* hice_in         =xNew<IssmDouble>(NUM_VERTICES);
    42624268        IssmDouble* hsnow_in        =xNew<IssmDouble>(NUM_VERTICES);
     4269        IssmDouble* qmr_in          =xNew<IssmDouble>(NUM_VERTICES);
    42634270
    42644271        // outputs
    … …  
    42714278        IssmDouble* albedo_out =xNew<IssmDouble>(NUM_VERTICES); memset(albedo_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    42724279        IssmDouble* albedo_snow_out =xNew<IssmDouble>(NUM_VERTICES); memset(albedo_snow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4273         IssmDouble* hsnow_out    =xNew<IssmDouble>(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4280        IssmDouble* hsnow_out   =xNew<IssmDouble>(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    42744281        IssmDouble* hice_out    =xNew<IssmDouble>(NUM_VERTICES); memset(hice_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4282        IssmDouble* qmr_out     =xNew<IssmDouble>(NUM_VERTICES); memset(qmr_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    42754283
    42764284        IssmDouble rho_water,rho_ice,desfac,desfacElev,rlaps,rdl;
    … …  
    43254333        GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum);
    43264334
    4327         if(VerboseSmb() && this->Sid()==0){
     4335        if(isverbose && this->Sid()==0){
    43284336                _printf0_("smb core: allocate inputs.\n");
    43294337                _printf0_("smb core: time_yr  : " << time_yr/yts <<"\n");
    … …  
    43454353
    43464354        /*temporal Enum depending on time*/
    4347         int enum_temp =TemperatureSEMICEnum;
    4348         int enum_hice =SmbHIceEnum;
    4349         int enum_hsnow=SmbHSnowEnum;
    4350         int enum_albedo=SmbAlbedoEnum;
     4355        int enum_temp       =TemperatureSEMICEnum;
     4356        int enum_hice       =SmbHIceEnum;
     4357        int enum_hsnow      =SmbHSnowEnum;
     4358        int enum_albedo     =SmbAlbedoEnum;
    43514359        int enum_albedo_snow=SmbAlbedoSnowEnum;
     4360        int enum_qmr        =SmbSemicQmrEnum;
    43524361        if (tstart+dt == time) {
    43534362                /* Load inital value at first time step*/
    … …  
    43574366                enum_albedo=SmbAlbedoInitEnum;
    43584367                enum_albedo_snow=SmbAlbedoSnowInitEnum;
     4368                enum_qmr        =SmbSemicQmrInitEnum;
    43594369        }
     4370        if(isverbose && this->Sid()==0)_printf0_("smb core: assign temp.\n");
    43604371        Input* tsurf_input       = this->GetInput(enum_temp); _assert_(tsurf_in);
     4372        if(isverbose && this->Sid()==0)_printf0_("smb core: assign mask.\n");
    43614373        Input* mask_input        = this->GetInput(SmbMaskEnum); _assert_(mask_input);
     4374        if(isverbose && this->Sid()==0)_printf0_("smb core: assign Tamp.\n");
    43624375        Input* Tamp_input        = this->GetInput(SmbTampEnum); _assert_(Tamp_input);
     4376        if(isverbose && this->Sid()==0)_printf0_("smb core: assign albedo.\n");
    43634377        Input* albedo_input      = this->GetInput(enum_albedo); _assert_(albedo_input);
    43644378        Input* albedo_snow_input = this->GetInput(enum_albedo_snow); _assert_(albedo_snow_input);
    43654379        Input* hice_input        = this->GetInput(enum_hice); _assert_(hice_input);
    43664380        Input* hsnow_input       = this->GetInput(enum_hsnow); _assert_(hsnow_input);
    4367 
    4368         if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: assign inputs.\n");
     4381        if(isverbose && this->Sid()==0)_printf0_("smb core: assign qmr.\n");
     4382        Input* qmr_input         = this->GetInput(enum_qmr); _assert_(qmr_input);
     4383
     4384        if(isverbose && this->Sid()==0)_printf0_("smb core: assign inputs done....\n");
    43694385        for(int iv=0;iv<NUM_VERTICES;iv++){
    43704386                gauss->GaussVertex(iv);
    … …  
    43894405                hsnow_input->GetInputValue(&hsnow_in[iv],gauss);
    43904406                hice_input->GetInputValue(&hice_in[iv],gauss);
     4407                qmr_input->GetInputValue(&qmr_in[iv],gauss);
    43914408
    43924409                /* Surface temperature correction */
    … …  
    44074424                dailydlradiation[iv]=dailydlradiation[iv]+rdl*st[iv];
    44084425        }
    4409         if(VerboseSmb() && this->Sid()==0){
     4426        if(isverbose && this->Sid()==0){
    44104427                _printf0_("smb core: assign tsurf_in        :" << tsurf_in[0] << "\n");
    44114428                _printf0_("smb core: assign dailytemperature:" << dailytemperature[0] << "\n");
    … …  
    44134430                _printf0_("smb core: assign hice            :" << hice_in[0] << "\n");
    44144431                _printf0_("smb core: assign mask            :" << mask_in[0] << "\n");
     4432                _printf0_("smb core: assign Tamp            :" << Tamp_in[0] << "\n");
    44154433                _printf0_("smb core: assign albedo          :" << albedo_in[0] << "\n");
    44164434                _printf0_("smb core: assign albedo_snow     :" << albedo_snow_in[0] << "\n");
    44174435                _printf0_("smb core: assign albedo_scheme   :" << alb_scheme  << "\n");
    4418         }
    4419 
    4420         if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: call semic module.\n");
     4436                _printf0_("smb core: assign qmr             :" << qmr_in[0]  << "\n");
     4437        }
     4438
     4439        if(isverbose && this->Sid()==0)_printf0_("smb core: call run_semic_transient module.\n");
     4440        /* call semic */
     4441        int nx=NUM_VERTICES, ntime=1, nloop=1;
     4442        bool semic_verbose=false; //VerboseSmb();
     4443        run_semic_transient_(&nx, &ntime, &nloop,
     4444                        dailysnowfall,  dailyrainfall, dailydsradiation, dailydlradiation,
     4445                        dailywindspeed, dailypressure, dailyairdensity,  dailyairhumidity, dailytemperature, tsurf_in, qmr_in,
     4446                        &dt,
     4447                        &hcrit, &rcrit,
     4448                        mask_in, hice_in, hsnow_in,
     4449                        albedo_in, albedo_snow_in,
     4450                        &alb_scheme, &alb_smax, &alb_smin, &albi, &albl,
     4451                        Tamp_in,
     4452                        &tmin, &tmax, &tmid, &mcrit, &wcrit, &tau_a, &tau_f, &afac, &semic_verbose,
     4453                        tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out,albedo_out, albedo_snow_out, hsnow_out, hice_out, qmr_out);
     4454
    44214455        for (int iv = 0; iv<NUM_VERTICES; iv++){
    4422                 /* call semic */
    4423                 run_semic_transient_(&dailysnowfall[iv], &dailyrainfall[iv], &dailydsradiation[iv], &dailydlradiation[iv],
    4424                                 &dailywindspeed[iv], &dailypressure[iv], &dailyairdensity[iv], &dailyairhumidity[iv], &dailytemperature[iv], &tsurf_in[iv], &dt,
    4425                                 &hcrit, &rcrit,
    4426                                 &mask_in[iv], &hice_in[iv], &hsnow_in[iv],
    4427                                 &albedo_in[iv], &albedo_snow_in[iv],
    4428                                 &alb_scheme, &alb_smax, &alb_smin, &albi, &albl,
    4429                                 &Tamp_in[iv],
    4430                                 &tmin, &tmax, &tmid, &mcrit, &wcrit, &tau_a, &tau_f, &afac,
    4431                                 &tsurf_out[iv], &smb_out[iv], &smbi_out[iv], &smbs_out[iv], &saccu_out[iv], &smelt_out[iv],&albedo_out[iv], &albedo_snow_out[iv], &hsnow_out[iv], &hice_out[iv]);
    4432 
    44334456                /*
    44344457                 unit conversion: water -> ice
    … …  
    44404463        }
    44414464
    4442         if(VerboseSmb() && this->Sid()==0){
    4443                 _printf0_("smb core: tsurf_out[0] " << tsurf_out[0] << "\n");
    4444                 _printf0_("smb core: hice_out[0]  " << hice_out[0] << "\n");
    4445                 _printf0_("smb core: hsnow_out[0] " << hsnow_out[0] << "\n");
    4446                 _printf0_("smb core: smb_ice[0]   " << smbi_out[0]*yts << "\n");
     4465        if(isverbose && this->Sid()==0){
     4466                _printf0_("smb core: tsurf_out " << tsurf_out[0] << " " << tsurf_out[1] << " " << tsurf_out[2] << "\n");
     4467                _printf0_("smb core: hice_out  " << hice_out[0] << " " <<  hice_out[1] << " " << hice_out[2] << "\n");
     4468                _printf0_("smb core: hsnow_out " << hsnow_out[0] << "\n");
     4469                _printf0_("smb core: smb_ice   " << smbi_out[0]*yts << "\n");
     4470                _printf0_("smb core: smb_ice   " << albedo_out[0] <<" "<<albedo_out[1] << " " << albedo_out[2] << "\n");
    44474471        }
    44484472
    44494473        switch(this->ObjectEnum()){
    44504474                case TriaEnum:
    4451                         this->AddInput(TemperatureSEMICEnum,&tsurf_out[0],P1DGEnum);
     4475                        this->AddInput(TemperatureSEMICEnum,  &tsurf_out[0],P1DGEnum);
    44524476                        // SMBout = SMB_ice + SMB_snow values.
    44534477                        //this->AddInput(SmbMassBalanceTotalEnum,&smb_out[0],P1DGEnum);
    44544478                        // water equivalent SMB ice to ice equivalent.
    4455                         this->AddInput(SmbMassBalanceEnum,&smbi_out[0],P1DGEnum);
     4479                        this->AddInput(SmbMassBalanceEnum,    &smbi_out[0],P1DGEnum);
    44564480                        this->AddInput(SmbMassBalanceSnowEnum,&smbs_out[0],P1DGEnum);
    44574481                        //this->AddInput(SmbMassBalanceSnowEnum,&smbs_out[0],P1DGEnum);
    … …  
    44594483                        this->AddInput(SmbAccumulationEnum,&saccu_out[0],P1DGEnum);
    44604484                        // smelt
    4461                         this->AddInput(SmbMeltEnum,&smelt_out[0],P1DGEnum);
    4462                         this->AddInput(SmbAlbedoEnum,&albedo_out[0],P1DGEnum);
    4463                         this->AddInput(SmbAlbedoSnowEnum,&albedo_snow_out[0],P1DGEnum);
    4464                         this->AddInput(SmbHSnowEnum,&hsnow_out[0],P1DGEnum);
    4465                         this->AddInput(SmbHIceEnum,&hice_out[0],P1DGEnum);
     4485                        this->AddInput(SmbMeltEnum,        &smelt_out[0],P1DGEnum);
     4486                        this->AddInput(SmbAlbedoEnum,      &albedo_out[0],P1DGEnum);
     4487                        this->AddInput(SmbAlbedoSnowEnum,  &albedo_snow_out[0],P1DGEnum);
     4488                        this->AddInput(SmbHSnowEnum,       &hsnow_out[0],P1DGEnum);
     4489                        this->AddInput(SmbHIceEnum,        &hice_out[0],P1DGEnum);
     4490                        this->AddInput(SmbSemicQmrEnum,    &qmr_out[0],P1DGEnum);
    44664491                        break;
    44674492                case PentaEnum:
    … …  
    44984523        xDelete<IssmDouble>(hsnow_out);
    44994524        xDelete<IssmDouble>(hice_out);
     4525        xDelete<IssmDouble>(qmr_out);
    45004526
    45014527        /*for inputs*/
    … …  
    45074533        xDelete<IssmDouble>(albedo_snow_in);
    45084534        xDelete<IssmDouble>(tsurf_in);
     4535        xDelete<IssmDouble>(qmr_in);
    45094536
    45104537        /* for inputs:geometry */

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90

    r27587 r27646  
    1 subroutine run_semic_transient(sf_in, rf_in, swd_in, lwd_in, wind_in, &
    2       sp_in, rhoa_in, qq_in, tt_in, tsurf_in, tstic, &
     1subroutine run_semic_transient(nx, ntime, nloop, sf_in, rf_in, swd_in, lwd_in, wind_in, &
     2      sp_in, rhoa_in, qq_in, tt_in, tsurf_in, qmr_in, &
     3      tstic, &
    34      hcrit, rcrit, &
    45      mask, hice, hsnow, &
    … …  
    67      alb_scheme, alb_smax, alb_smin, albi, albl, &
    78      Tamp, &
    8       tmin, tmax, tmid, mcrit, w_crit, tau_a, tau_f, afac, &
    9       tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, alb_out, alb_snow_out,hsnow_out,hice_out)
     9      tmin, tmax, tmid, mcrit, w_crit, tau_a, tau_f, afac, verbose, &
     10      tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, alb_out, &
     11      alb_snow_out,hsnow_out,hice_out,qmr_out) !{{{
    1012
    1113   use utils
    1214   use surface_physics
    13 
    1415   implicit none
    1516
    … …  
    1718   type(surface_physics_class) :: surface
    1819   ! declare forcing class
    19    type(forc_class) :: forc
     20   !type(forc_class) :: forc
    2021   ! declare validation class
    2122   !type(vali_class) :: vali    ! validation not needed here
    2223
    2324   integer, parameter:: dp=kind(0.d0)  !< define precision (machine specific)
    24    integer :: i, k, nx, nloop, year
    25    integer :: day=1 !< not used value.
     25   integer :: i, k, n
     26   integer :: nnx, nny
     27   integer :: year=0
     28   integer :: day =1 !< not used value.
     29   integer, intent(in) :: nx, ntime, nloop   ! number of grid / number of times
     30   logical, intent(in) :: verbose            ! verbosity
     31   logical :: debug=.false.
    2632
    2733   ! forcing data   
    28    double precision, intent(in) :: sf_in  ! snow fall rate [m/s]
    29    double precision, intent(in) :: rf_in  ! rain fall rate [m/s]
    30    double precision, intent(in) :: swd_in ! downwelling shortwave radiation [W/m2]
    31    double precision, intent(in) :: lwd_in ! downwelling longwave radiation [W/m2]
    32    double precision, intent(in) :: wind_in! surface wind speed [m/s]
    33    double precision, intent(in) :: sp_in  ! surface pressure [Pa]
    34    double precision, intent(in) :: rhoa_in! air density [kg/m3]
    35    double precision, intent(in) :: qq_in  ! air specific humidity [kg/kg]
    36    double precision, intent(in) :: tt_in  ! air temperature [K]
     34   ! input argument format array size (nx, ntime)...
     35   double precision, intent(in), dimension(nx):: sf_in    ! snow fall rate [m/s]
     36   double precision, intent(in), dimension(nx):: rf_in    ! rain fall rate [m/s]
     37   double precision, intent(in), dimension(nx):: swd_in   ! downwelling shortwave radiation [W/m2]
     38   double precision, intent(in), dimension(nx):: lwd_in   ! downwelling longwave radiation [W/m2]
     39   double precision, intent(in), dimension(nx):: wind_in  ! surface wind speed [m/s]
     40   double precision, intent(in), dimension(nx):: sp_in    ! surface pressure [Pa]
     41   double precision, intent(in), dimension(nx):: rhoa_in  ! air density [kg/m3]
     42   double precision, intent(in), dimension(nx):: qq_in    ! air specific humidity [kg/kg]
     43   double precision, intent(in), dimension(nx):: tt_in    ! air temperature [K]
    3744
    3845   ! input data
    3946   double precision, intent(in) :: tstic  ! time step from ISSM [sec].
    40    double precision, intent(in) :: tsurf_in ! input temperature [K]
    4147
    4248   ! output data
    43    double precision, intent(out) :: tsurf_out  ! Ice surface Temperature [K]
    44    double precision, intent(out) :: smb_out    ! surface mass balance=(Accu-Melt) [m/s]
    45    double precision, intent(out) :: smbi_out   ! SMB ice  [water equivalent m/s]
    46    double precision, intent(out) :: smbs_out   ! SMB snow [water equivalent m/s]
    47    double precision, intent(out) :: saccu_out  ! accumulation [m/s]
    48    double precision, intent(out) :: smelt_out  ! ablation [m/s]
    49    double precision, intent(out) :: alb_out    ! grid-averaged albedo [no unit]
    50    double precision, intent(out) :: alb_snow_out 
    51    double precision :: qmr_out   
    52    !double precision, intent(out) :: qmr_out   
    53    double precision, intent(out) :: hice_out   
    54    double precision, intent(out) :: hsnow_out
     49   ! Ice surface Temperature [K]
     50   double precision, intent(out), dimension(nx):: tsurf_out   
     51   ! surface mass balance=(Accu-Melt) [m/s]
     52   double precision, intent(out), dimension(nx):: smb_out     
     53   double precision, intent(out), dimension(nx):: smbi_out     ! SMB ice  [water equivalent m/s]
     54   double precision, intent(out), dimension(nx):: smbs_out     ! SMB snow [water equivalent m/s]
     55   double precision, intent(out), dimension(nx):: saccu_out    ! accumulation [m/s]
     56   double precision, intent(out), dimension(nx):: smelt_out    ! ablation [m/s]
     57   double precision, intent(out), dimension(nx):: alb_out      ! grid-averaged albedo [no unit]
     58   double precision, intent(out), dimension(nx):: alb_snow_out
     59   double precision, intent(out), dimension(nx):: hice_out   
     60   double precision, intent(out), dimension(nx):: hsnow_out   
     61   double precision, intent(out), dimension(nx):: qmr_out     
    5562
    5663   double precision :: total_time, start, finish
    5764
    5865   ! set parameters
    59    character (len=256) :: name         ! not used(?)
    60    character (len=256) :: boundary(30) ! not used(?)
     66   !character (len=256) :: name         ! not used(?)
     67   !character (len=256) :: boundary(30) ! not used(?)
    6168   !character (len=256), intent(in) :: alb_scheme  !< name of albedo scheme
    6269   integer, intent(in)          :: alb_scheme
    63    integer :: n_ksub   
    64    double precision             :: ceff         !< surface heat heat capacity of snow/ice [J/W m2]
    65    double precision, intent(in) :: hsnow
    66    double precision, intent(in) :: hice
     70   !integer :: n_ksub   
     71   !double precision             :: ceff         !< surface heat heat capacity of snow/ice [J/W m2]
     72   double precision, intent(in), dimension(nx):: albedo
     73   double precision, intent(in), dimension(nx):: albedo_snow !< spatial..
     74   double precision, intent(in), dimension(nx):: hsnow
     75   double precision, intent(in), dimension(nx):: hice
     76   double precision, intent(in), dimension(nx):: tsurf_in    !< input temperature [K]
     77   double precision, intent(in), dimension(nx):: qmr_in
     78   double precision, intent(in), dimension(nx):: mask
     79
    6780   double precision, intent(in) :: albi
    6881   double precision, intent(in) :: albl
    69    double precision, intent(in) :: mask
    70    double precision, intent(in) :: albedo 
    71    double precision, intent(in) :: albedo_snow
    7282   double precision, intent(in) :: alb_smax
    7383   double precision, intent(in) :: alb_smin
    … …  
    7585   double precision, intent(in) :: rcrit !< critical snow height fro which refreezing fraction is 50% [m]
    7686   double precision, intent(in) :: Tamp
    77    double precision    :: csh
    78    double precision    :: clh
     87   !double precision    :: csh
     88   !double precision    :: clh
    7989   double precision, intent(in) :: tmin
    8090   double precision, intent(in) :: tmax
    81    double precision    :: tsticsub
     91   !double precision    :: tsticsub
    8292   ! parameters for isba albedo scheme.
    8393   double precision, intent(in) :: tau_a  !< critical liquide water concent for "isba" albedo scheme [kg/m2]
    … …  
    8898   double precision, intent(in) :: tmid !< param for "alex" albedo parameterization [K]
    8999
    90    nx = 1
    91    year = 0
     100   if (debug) then
     101      print*,'   ntime: ', ntime
     102      print*,'   nx   : ', nx
     103   end if
    92104
    93105   ! set vector length
    94106   surface%par%nx = nx
    95107
    96    ! set input (forcing data)
    97    !forc%sf(:,:) = sf_in      ! snowfall flux [m/sec]
    98    !forc%rf(:,:) = rf_in      ! rainfall flux [m/sec]
    99    !forc%swd(:,:) = swd_in    ! short radiation
    100    !forc%lwd(:,:) = lwd_in    ! long radiation downward
    101    !forc%wind(:,:) = wind_in  ! wind speed.
    102    !forc%sp(:,:) = sp_in      ! surface pressure
    103    !forc%rhoa(:,:) = rhoa_in  ! air density [kg/m3]
    104    !forc%qq(:,:) = qq_in      ! air specific humidity [kg/kg]
    105    !forc%tt(:,:) = tt_in      ! 2m air temperature
    106 
    107108   ! FIXME should be user input
    108109   !boundary = "" "" ""
    109    !surface%par%tstic = 86400.0_dp !< time step [s] in one day.
     110   if (debug) then
     111      print*, "run_semic_transient: initialize parameters."
     112   end if
    110113   surface%par%tstic = tstic      !< time step [s]
    111114   surface%par%ceff= 2.0e6_dp     !< surface heat capacity of snow/ice [J/K/m2]
    … …  
    114117   surface%par%alb_smax = alb_smax !0.79_dp !< max snow albedo
    115118   surface%par%alb_smin = alb_smin !0.6_dp  !< min snow albedo
    116    surface%par%albi = albi !0.41_dp     !< albedo for ice
    117    surface%par%albl = albl !0.07_dp     !< albedo for land
     119   surface%par%albi = albi ! 0.41_dp     !< albedo for ice
     120   surface%par%albl = albl ! 0.07_dp     !< albedo for land
    118121   surface%par%tmin = tmin ! -999_dp
    119122   surface%par%tmax = tmax ! 273.15_dp
    … …  
    122125   surface%par%amp   = Tamp !3.0_dp   !< amplitude of diurnal cycle [K]
    123126   if (alb_scheme == 0) then
    124       surface%par%alb_scheme="None"
     127      surface%par%alb_scheme="none"
    125128   else if (alb_scheme == 1) then
    126129      surface%par%alb_scheme = "slater"
    … …  
    129132   else if (alb_scheme == 3) then
    130133      surface%par%alb_scheme = "isba"
     134   else
     135      print*, "ERROR: current albedo scheme is not available."
     136      call exit(1)
    131137   end if
    132138   surface%par%tau_a  = tau_a  !0.008_dp
    … …  
    134140   surface%par%w_crit = w_crit !15.0_dp ! default value
    135141   surface%par%mcrit  = mcrit  !6.0e-8_dp
    136    surface%par%n_ksub = 3.0_dp
     142   surface%par%n_ksub = 3      ! sub ...
    137143   ! snow albedo of alex
    138144   surface%par%afac   = afac
    … …  
    146152
    147153   ! initialise prognostic variables
     154   if (debug) then
     155      print*,"run_semic_transient: initialize variables."
     156   end if
    148157   ! these values will be updated through "surface_energy_and_mass_balance" function.
    149    surface%now%mask(:)     = mask        ! 2.0_dp  !loi_mask(:nx)
    150    surface%now%hsnow(:)    = hsnow       ! initial snow height...
    151    surface%now%hice(:)     = hice        ! initial ice height..
    152    surface%now%tsurf(:)    = tsurf_in
    153    surface%now%alb(:)      = albedo      !< initial albedo for energy balance.
    154    surface%now%alb_snow(:) = albedo_snow !< initial albedo for ISBA albedo method.
    155    surface%now%qmr_res(:)  = 0.0_dp ! residual heat.
    156 
    157    ! initialize with zeros
    158    surface%now%qmr(:) = 0.0_dp
    159 
    160    tsurf_out = 0.0_dp
    161    smb_out   = 0.0_dp
    162    smbi_out  = 0.0_dp
    163    smbs_out  = 0.0_dp
    164    saccu_out = 0.0_dp
    165    smelt_out = 0.0_dp
    166    alb_out      = 0.0_dp
    167    alb_snow_out = 0.0_dp
     158   surface%now%mask    (:) = mask       (:) ! 2.0_dp  !loi_mask(:nx)
     159   if (debug) then
     160      print*,"run_semic_transient: initialize variables: mask"
     161   end if
     162   surface%now%hsnow   (:) = hsnow      (:) ! initial snow height...
     163   surface%now%hice    (:) = hice       (:) ! initial ice height..
     164   surface%now%tsurf   (:) = tsurf_in   (:) !< initial ice surface temperature
     165   surface%now%alb     (:) = albedo     (:) !< initial albedo for energy balance.
     166   surface%now%alb_snow(:) = albedo_snow(:) !< initial albedo for ISBA albedo method.
     167   if (debug) then
     168      print*,"run_semic_transient: initialize variables. DONE."
     169   end if
     170
     171   if (debug) then
     172      !print*, "====== global variable =========="
     173      !print*, "nloop          :", nloop
     174      !print*, "nx             :", surface%par%nx
     175      !print*, "======  parameters ======"
     176      !print*, "csh            :", surface%par%csh
     177      !print*, "clh            :", surface%par%clh
     178      !print*, "albeo scheme   :", surface%par%alb_scheme
     179      !print*, "albeo ice      :", surface%par%albi
     180      !print*, "tstic          :", surface%par%tstic
     181      !print*, "tsticsub       :", surface%par%tsticsub
     182      !print*, "n_ksub         :", surface%par%n_ksub
     183      print*, "====== inputs ========="
     184      print*, "hsnow          :", hsnow
     185      print*, "======  state variables ======"
     186      print*, "hsnow          :", surface%now%hsnow
     187      print*, "hice           :", surface%now%hice
     188      print*, "albeo          :", surface%now%alb
     189      print*, "albeo snow     :", surface%now%alb_snow
     190      print*, "mask           :", surface%now%mask
     191      print*, "tsurf          :", surface%now%tsurf
     192      print*, "sf             :", sf_in
     193   end if
    168194
    169195   ! define boundary conditions (not used, here!)
    … …  
    171197   !call print_boundary_opt(surface%bnd)
    172198
    173    
    174199   ! input with single value
    175    surface%now%sf   = sf_in
    176    surface%now%rf   = rf_in
    177    surface%now%sp   = sp_in
    178    surface%now%lwd  = lwd_in
    179    surface%now%swd  = swd_in
    180    surface%now%wind = wind_in
    181    surface%now%rhoa = rhoa_in
    182    surface%now%t2m  = tt_in
    183    surface%now%qq   = qq_in
    184 
    185    ! calculate prognostic and diagnsotic variables
    186    call cpu_time(start)
    187    call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,day,-1)
    188    call cpu_time(finish)
    189    total_time = total_time + (finish - start)
    190 
    191    tsurf_out=surface%now%tsurf(1)
    192    ! melt - potential surface melt [m/s]
    193    ! smb = SMB_ice + SMB_snow
    194    ! smbi  - SMB_ice  (water equivalent m/sec)
    195    ! smbs  - SMB_snow (water equivalent m/sec)
    196    smb_out   =surface%now%smb(1)      ! smb = smb_snow + smb_ice
    197    smbi_out  =surface%now%smb_ice(1)  ! Csi (snow>ice) - melted_ice + refrezon_rain.
    198    smbs_out  =surface%now%smb_snow(1) ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow
    199    saccu_out =surface%now%acc(1)      ! acc      = snowfall - sublimiation - refreezing
    200    smelt_out =surface%now%melt(1)     ! potential surface melt = melt_ice + melt_snow
    201    alb_out      =surface%now%alb(1)     
    202    alb_snow_out =surface%now%alb_snow(1)
    203    qmr_out    =surface%now%qmr(1)
    204    hsnow_out  =surface%now%hsnow(1)
    205    hice_out   =surface%now%hice(1)
     200   do k =1,nloop
     201      do i =1,ntime
     202         if (debug) then
     203            print*,"run_semic_transient: forcing data: ntime = ", i
     204         end if
     205         surface%now%sf   = sf_in  !(:,i)
     206         surface%now%rf   = rf_in  !(:,i)
     207         surface%now%sp   = sp_in  !(:,i)
     208         surface%now%lwd  = lwd_in !(:,i)
     209         surface%now%swd  = swd_in !(:,i)
     210         surface%now%wind = wind_in!(:,i)
     211         surface%now%rhoa = rhoa_in!(:,i)
     212         surface%now%t2m  = tt_in  !(:,i)
     213         surface%now%qq   = qq_in  !(:,i)
     214         ! qmr_res is used to "energy_balance" in semic.
     215         surface%now%qmr_res = qmr_in
     216
     217         ! calculate prognostic and diagnsotic variables
     218         call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,i,year)
     219         
     220         if (debug) then
     221            print*,"done..."
     222         end if
     223         if (k == nloop) then
     224            tsurf_out         = surface%now%tsurf
     225            ! melt - potential surface melt [m/s]
     226            ! smb = SMB_ice + SMB_snow
     227            ! smbi  - SMB_ice  (water equivalent m/sec)
     228            ! smbs  - SMB_snow (water equivalent m/sec)
     229            smb_out           =surface%now%smb      ! smb = smb_snow + smb_ice
     230            smbi_out          =surface%now%smb_ice  ! Csi (snow>ice) - melted_ice + refrezon_rain.
     231            smbs_out          =surface%now%smb_snow ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow
     232            saccu_out         =surface%now%acc      ! acc      = snowfall - sublimiation - refreezing
     233            smelt_out         =surface%now%melt     ! potential surface melt = melt_ice + melt_snow
     234            alb_out           =surface%now%alb
     235            alb_snow_out      =surface%now%alb_snow
     236            hsnow_out         =surface%now%hsnow
     237            hice_out          =surface%now%hice
     238            qmr_out           =surface%now%qmr_res
     239         end if
     240      end do
     241   end do
    206242
    207243   ! de-allocate surface_physics arrays
    208244   call surface_dealloc(surface%now)
    209245
    210    !write(*,*) 'total time for surface_physics:', nloop, total_time
    211 
    212 end subroutine run_semic_transient
     246end subroutine run_semic_transient ! }}}

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r27510 r27646  
    11181118syn keyword cConstant SmbS0pEnum
    11191119syn keyword cConstant SmbS0tEnum
     1120syn keyword cConstant SmbSemicQmrEnum
     1121syn keyword cConstant SmbSemicQmrInitEnum
    11201122syn keyword cConstant SmbSizeiniEnum
    11211123syn keyword cConstant SmbSmbCorrEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r27510 r27646  
    11151115        SmbS0pEnum,
    11161116        SmbS0tEnum,
     1117        SmbSemicQmrEnum,
     1118        SmbSemicQmrInitEnum,
    11171119        SmbSizeiniEnum,
    11181120        SmbSmbCorrEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r27510 r27646  
    11201120                case SmbS0pEnum : return "SmbS0p";
    11211121                case SmbS0tEnum : return "SmbS0t";
     1122                case SmbSemicQmrEnum : return "SmbSemicQmr";
     1123                case SmbSemicQmrInitEnum : return "SmbSemicQmrInit";
    11221124                case SmbSizeiniEnum : return "SmbSizeini";
    11231125                case SmbSmbCorrEnum : return "SmbSmbCorr";

  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r27510 r27646  
    11111111syn keyword juliaConstC SmbS0pEnum
    11121112syn keyword juliaConstC SmbS0tEnum
     1113syn keyword juliaConstC SmbSemicQmrEnum
     1114syn keyword juliaConstC SmbSemicQmrInitEnum
    11131115syn keyword juliaConstC SmbSizeiniEnum
    11141116syn keyword juliaConstC SmbSmbCorrEnum

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r27510 r27646  
    11471147              else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
    11481148              else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
     1149              else if (strcmp(name,"SmbSemicQmr")==0) return SmbSemicQmrEnum;
     1150              else if (strcmp(name,"SmbSemicQmrInit")==0) return SmbSemicQmrInitEnum;
    11491151              else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    11501152              else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
    … …  
    12421244              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    12431245              else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
    1244               else if (strcmp(name,"VyDebris")==0) return VyDebrisEnum;
    1245               else if (strcmp(name,"Vy")==0) return VyEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
     1249              if (strcmp(name,"VyDebris")==0) return VyDebrisEnum;
     1250              else if (strcmp(name,"Vy")==0) return VyEnum;
     1251              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    12501252              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    12511253              else if (strcmp(name,"VyShear")==0) return VyShearEnum;
    … …  
    13651367              else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
    13661368              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
    1367               else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    1368               else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
     1372              if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
     1373              else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
     1374              else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
    13731375              else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
    13741376              else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum;
    … …  
    14881490              else if (strcmp(name,"Fset")==0) return FsetEnum;
    14891491              else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
    1490               else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum;
    1491               else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
     1495              if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum;
     1496              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
     1497              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    14961498              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    14971499              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    … …  
    16111613              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    16121614              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    1613               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    1614               else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
     1618              if (strcmp(name,"MinVy")==0) return MinVyEnum;
     1619              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
     1620              else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
    16191621              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    16201622              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    … …  
    17341736              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    17351737              else if (strcmp(name,"Tetra")==0) return TetraEnum;
    1736               else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
    1737               else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    17381738         else stage=15;
    17391739   }
    17401740   if(stage==15){
    1741               if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
     1741              if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
     1742              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
     1743              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    17421744              else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    17431745              else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
Note: See TracChangeset for help on using the changeset viewer.