Ignore:
Timestamp:
01/17/24 17:08:12 (14 months ago)
Author:
sehrenfe
Message:

CHG: GlaDS- new input parameter rheology_B_base and new optional output HydrologySheetDischarge added to hydrology model

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r27937 r28050  
    151151        iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
    152152        iomodel->FetchDataToInput(inputs,elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);
     153        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.rheology_B_base",HydrologyRheologyBBaseEnum);
     154
    153155        if(iomodel->domaintype==Domain2DhorizontalEnum){
    154156                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
    … …  
    249251        IssmDouble g         = element->FindParam(ConstantsGEnum);
    250252        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
    251         Input* hr_input      = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
     253        Input* hr_input  = element->GetInput(HydrologyBumpHeightEnum);       _assert_(hr_input);
    252254        Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
    253         Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
    254         Input* h_input   = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input);
    255         Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
    256         Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
    257         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    258         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     255        Input* phi_input = element->GetInput(HydraulicPotentialEnum);        _assert_(phi_input);
     256        Input* h_input   = element->GetInput(HydrologySheetThicknessEnum);   _assert_(h_input);
     257        Input* H_input   = element->GetInput(ThicknessEnum);                _assert_(H_input);
     258        Input* b_input   = element->GetInput(BedEnum);                      _assert_(b_input);
     259        Input* B_input   = element->GetInput(HydrologyRheologyBBaseEnum);    _assert_(B_input);
     260        Input* n_input   = element->GetInput(MaterialsRheologyNEnum);        _assert_(n_input);
    259261
    260262        /* Start  looping on the number of gaussian points: */
    … …  
    275277                b_input->GetInputValue(&b,gauss);
    276278                H_input->GetInputValue(&H,gauss);
    277 
    278                 /*Hard code B*/
    279                 B = Cuffey(273.15-2);
    280279               
    281280                /*Get norm of gradient of hydraulic potential and make sure it is >0*/
    … …  
    364363        IssmDouble g         = element->FindParam(ConstantsGEnum);
    365364        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
    366         Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    367         Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
    368         Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
    369         Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
    370         Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
     365        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);                _assert_(hr_input);
     366        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);            _assert_(h_input);
     367        Input* H_input      = element->GetInput(ThicknessEnum);                          _assert_(H_input);
     368        Input* b_input      = element->GetInput(BedEnum);                                _assert_(b_input);
     369        Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);        _assert_(G_input);
    371370        Input* melt_input   = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);_assert_(melt_input);
    372371        Input* RO_input     = NULL;
    373         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    374         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    375         Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
    376         Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
     372        Input* B_input      = element->GetInput(HydrologyRheologyBBaseEnum);             _assert_(B_input);
     373        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);                 _assert_(n_input);
     374        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);              _assert_(phiold_input);
     375        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);                 _assert_(phi_input);
    377376
    378377        /*Build friction element, needed later: */
    … …  
    397396                melt_input->GetInputValue(&melt,gauss);
    398397
    399                 /*Hard code B*/
    400                 B = Cuffey(273.15-2);
    401 
    402398                /*Get basal velocity*/
    403399                friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss);
    … …  
    457453        element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
    458454
    459         /*Compute Hydrology Vx and Vy for time stepping purposes (These inputs do not affect GlaDS)*/
     455        /*Compute Hydrology Vx and Vy for time stepping purposes, and Sheet Discharge as an optional output (These inputs do not affect GlaDS)*/
    460456
    461457        /*Intermediaries*/
    462458   IssmDouble  dphi[3],h,k,phi;
    463         IssmDouble  h_r;
     459        IssmDouble  h_r;
    464460        IssmDouble  oceanLS,iceLS;
    465461        IssmDouble* xyz_list = NULL;
    … …  
    468464        int numvertices = element->GetNumberOfVertices();
    469465
    470         /*Initialize new sheet thickness*/
     466        /*Initialize water sheet velocity and discharge*/
    471467        IssmDouble* vx = xNew<IssmDouble>(numvertices);
    472468        IssmDouble* vy = xNew<IssmDouble>(numvertices);
     469        IssmDouble* d = xNew<IssmDouble>(numvertices);
    473470
    474471        /*Set to 0 if inactive element*/
    … …  
    476473                for(int iv=0;iv<numvertices;iv++) vx[iv] = 0.;
    477474                for(int iv=0;iv<numvertices;iv++) vy[iv] = 0.;
     475                for(int iv=0;iv<numvertices;iv++) d[iv] = 0.;
    478476                element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
    479477                element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
     478                element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum);
    480479                xDelete<IssmDouble>(vx);
    481480                xDelete<IssmDouble>(vy);
     481                xDelete<IssmDouble>(d);
    482482                return;
    483483        }
    … …  
    499499        Input *iceLS_input   = element->GetInput(MaskIceLevelsetEnum);            _assert_(iceLS_input);
    500500
    501         /* Start  looping on the number of gaussian points: */
     501        /* Start looping on the number of gaussian points: */
    502502        Gauss* gauss=element->NewGauss();
    503503        for(int iv=0;iv<numvertices;iv++){
    … …  
    508508      phi_input->GetInputValue(&phi,gauss);
    509509      h_input->GetInputValue(&h,gauss);
    510         hr_input->GetInputValue(&h_r,gauss);
     510      hr_input->GetInputValue(&h_r,gauss);
    511511      k_input->GetInputValue(&k,gauss);
    512512                oceanLS_input->GetInputValue(&oceanLS,gauss);
    513513                iceLS_input->GetInputValue(&iceLS,gauss);
    514514
    515                 /*Set sheet thickness to zero if floating or no ice*/
     515                /*Set to zero if floating or no ice*/
    516516                if(oceanLS<0. || iceLS>0.){
    517517                        vx[iv] = 0.;
    518518         vy[iv] = 0.;
     519                        d[iv] = 0.;
    519520                }
    520521                else{
    … …  
    526527         /*If omega is zero, use standard model, otherwise transition model*/
    527528         IssmDouble nu = mu_water/rho_water;
    528         IssmDouble coeff;
    529         if(istransition==1 && omega>=AEPS){
    530                 IssmDouble hratio = fabs(h/h_r);
    531                 IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
    532                 coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi/max(AEPS,h);  // divide by h to get speed instead of discharge
    533         }
    534         else {
    535                 coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h);  // divide by h to get speed instead of discharge
    536         }
    537 
    538 
    539                         vx[iv] = -coeff*dphi[0];
    540                         vy[iv] = -coeff*dphi[1];
     529                        IssmDouble coeff;
     530                        if(istransition==1 && omega>=AEPS){
     531                                IssmDouble hratio = fabs(h/h_r);
     532                                IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
     533                                coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi;  // coeff gives discharge; divide by h to get speed instead of discharge
     534                        }
     535                        else {
     536                        coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);  // coeff gives discharge; divide by h to get speed instead of discharge
     537                        }
     538
     539                        vx[iv] = -coeff/max(AEPS,h)*dphi[0];
     540                        vy[iv] = -coeff/max(AEPS,h)*dphi[1];
     541                               
     542                        d[iv] = coeff*normgradphi;
    541543                }
    542544        }
    … …  
    544546        element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
    545547        element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
     548        element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum);
    546549
    547550        /*Clean up and return*/
    … …  
    549552        xDelete<IssmDouble>(vx);
    550553        xDelete<IssmDouble>(vy);
     554        xDelete<IssmDouble>(d);
    551555        delete gauss;
    552556}/*}}}*/
    … …  
    646650        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
    647651        IssmDouble g         = element->FindParam(ConstantsGEnum);
    648         Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    649         Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input);
    650         Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input);
    651         Input* H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
    652         Input* b_input  = element->GetInput(BedEnum); _assert_(b_input);
    653         Input* hold_input  = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
    654         Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    655         Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     652        Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);         _assert_(hr_input);
     653        Input* vx_input = element->GetInput(VxBaseEnum);                      _assert_(vx_input);
     654        Input* vy_input = element->GetInput(VyBaseEnum);                      _assert_(vy_input);
     655        Input* H_input = element->GetInput(ThicknessEnum);                    _assert_(H_input);
     656        Input* b_input = element->GetInput(BedEnum);                          _assert_(b_input);
     657        Input* hold_input = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
     658        Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum);       _assert_(B_input);
     659        Input* n_input = element->GetInput(MaterialsRheologyNEnum);           _assert_(n_input);
    656660        Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
    657         Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
    658         Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
    659 
     661        Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum);      _assert_(oceanLS_input);
     662        Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum);          _assert_(iceLS_input);
    660663
    661664        /* Start  looping on the number of gaussian points: */
    … …  
    677680                iceLS_input->GetInputValue(&iceLS,gauss);
    678681
    679                 /*Hard code B*/
    680                 B = Cuffey(273.15-2);
    681 
    682682                /*Set sheet thickness to zero if floating or no ice*/
    683683                if(oceanLS<0. || iceLS>0.){
Note: See TracChangeset for help on using the changeset viewer.