Changeset 23846


Ignore:
Timestamp:
04/15/19 19:44:09 (6 years ago)
Author:
schlegel
Message:

CHG: add GEMB output for surface fluxes

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

Legend:

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

    r23844 r23846  
    32413241        int        swIdx=0;
    32423242        IssmDouble cldFrac,t0wet, t0dry, K;
    3243         IssmDouble ulw=0.0;
    3244         IssmDouble netSW=0.0;
    3245         IssmDouble netLW=0.0;
    32463243        IssmDouble lhf=0.0;
    32473244        IssmDouble shf=0.0;
     
    32743271        IssmDouble* gsp = NULL;
    32753272        IssmDouble  EC = 0.0;
     3273        IssmDouble  ulw = 0.0;
     3274        IssmDouble  netSW=0.0;
     3275        IssmDouble  netLW=0.0;
     3276        IssmDouble  netULW=0.0;
     3277        IssmDouble  netLHF=0.0;
     3278        IssmDouble  netSHF=0.0;
    32763279        IssmDouble* W = NULL;
    32773280        IssmDouble* a = NULL;
     
    35233526
    35243527                /*Calculate net shortwave [W m-2]*/
    3525                 netSW = cellsum(swf,m);
     3528                netSW = netSW + cellsum(swf,m);
    35263529
    35273530                /*Thermal profile computation:*/
    3528                 if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
     3531                if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
    35293532
    35303533                /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
     
    35443547                /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
    35453548                 * sub-time step in thermo equations*/
    3546                 ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
     3549                //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
    35473550
    35483551                /*Calculate net longwave [W m-2]*/
    3549                 netLW = dlw - ulw;
     3552                netULW = netULW + ulw;
     3553                netLW = netLW + dlw - ulw;
    35503554
    35513555                /*Calculate turbulent heat fluxes [W m-2]*/
     
    35643568                                                << "gsp[" << cellsum(gsp,m)  << "] "
    35653569                                                << "swf[" << netSW << "] "
     3570                                                << "lwf[" << netLW << "] "
    35663571                                                << "a[" << a << "] "
    35673572                                                << "te[" << teValue << "] "
    35683573                                                << "\n");
    35693574                } /*}}}*/
     3575
     3576                netLHF = netLHF + lhf;
     3577                netSHF = netSHF + shf;
    35703578
    35713579                /*Sum component mass changes [kg m-2]*/
     
    36203628        this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/dt/rho_ice));
    36213629        this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/dt/rho_ice));
     3630        this->AddInput(new DoubleInput(SmbNetULWEnum,netULW));
     3631        this->AddInput(new DoubleInput(SmbNetLWEnum,netLW));
     3632        this->AddInput(new DoubleInput(SmbNetSWEnum,netSW));
     3633        this->AddInput(new DoubleInput(SmbNetLHFEnum,netLHF));
     3634        this->AddInput(new DoubleInput(SmbNetSHFEnum,netSHF));
    36223635        this->AddInput(new DoubleInput(SmbDzAddEnum,sumdz_add));
    36233636        this->AddInput(new DoubleInput(SmbMAddEnum,sumMassAdd/dt));
  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r23468 r23846  
    509509
    510510}  /*}}}*/
    511 void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/
     511void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/
    512512
    513513        /* ENGLACIAL THERMODYNAMICS*/
     
    536536        // T: grid cell temperature [k]
    537537        // EC: evaporation/condensation [kg]
     538        // ulwrf: upward longwave radiation flux [W m-2]
    538539
    539540        /*intermediary: */
     
    583584        IssmDouble EC=0.0;
    584585        IssmDouble* T=*pT;
     586        IssmDouble ulwrf=0.0;
    585587
    586588        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
     
    610612
    611613        // Bulk-transfer coefficient for turbulent fluxes
    612         An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
     614        //An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
     615        An =  pow(0.4,2) / (log(Tz/z0)*log(Vz/z0));     // Bulk-transfer coefficient
    613616        C = An * dAir * V;                        // shf & lhf common coefficient
    614617
     
    820823
    821824                // upward longwave contribution
    822                 ulw = - (SB * pow(Ts,4.0)* teValue) * dt ; //+20
     825                ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here
     826                ulwrf = ulwrf - ulw;
     827
    823828                dT_ulw = ulw / TCs;
    824829
     
    875880        *pEC=EC;
    876881        *pT=T;
     882        *pulwrf=ulwrf;
    877883
    878884}  /*}}}*/
     
    19391945                        case 6: // Ligtenberg and others (2011) [semi-emperical], Antarctica
    19401946                                // common variable
    1941                                 H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
     1947                                // From literature: H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
     1948                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
    19421949                                c0arth = 0.07 * H;
    19431950                                c1arth = 0.03 * H;
     
    19521959                        case 7: // Kuipers Munneke and others (2015) [semi-emperical], Greenland
    19531960                                // common variable
    1954                                 H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
     1961                                // From literature: H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
     1962                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
    19551963                                c0arth = 0.07 * H;
    19561964                                c1arth = 0.03 * H;
     
    20682076
    20692077        //// Bulk-transfer coefficient
    2070         An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
     2078        //An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
     2079        An =  pow(0.4,2) / (log(Tz/z0)*log(Vz/z0));     // Bulk-transfer coefficient
    20712080        C = An * d_air * V;             // shf & lhf common coefficient
    20722081
  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r23814 r23846  
    3131void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
    3232void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
    33 void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid);
     33void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid);
    3434void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
    3535void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid);
  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r23843 r23846  
    11351135        SmbDzAddEnum,
    11361136        SmbFACEnum,
     1137        SmbNetULWEnum,
     1138        SmbNetLWEnum,
     1139        SmbNetSWEnum,
     1140        SmbNetLHFEnum,
     1141        SmbNetSHFEnum,
    11371142        SMBforcingEnum,
    11381143        SMBgcmEnum,
  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r23843 r23846  
    11391139                case SmbDzAddEnum : return "SmbDzAdd";
    11401140                case SmbFACEnum : return "SmbFAC";
     1141                case SmbNetULWEnum : return "SmbNetULW";
     1142                case SmbNetLWEnum : return "SmbNetLW";
     1143                case SmbNetSWEnum : return "SmbNetSW";
     1144                case SmbNetLHFEnum : return "SmbNetLHF";
     1145                case SmbNetSHFEnum : return "SmbNetSHF";
    11411146                case SMBforcingEnum : return "SMBforcing";
    11421147                case SMBgcmEnum : return "SMBgcm";
  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r23843 r23846  
    11661166              else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
    11671167              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
     1168              else if (strcmp(name,"SmbNetULW")==0) return SmbNetULWEnum;
     1169              else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
     1170              else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum;
     1171              else if (strcmp(name,"SmbNetLHF")==0) return SmbNetLHFEnum;
     1172              else if (strcmp(name,"SmbNetSHF")==0) return SmbNetSHFEnum;
    11681173              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    11691174              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
     
    12391244              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
    12401245              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    1241               else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    12421250              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    12431251              else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
    12441252              else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
    12451253              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
     1254              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    12501255              else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    12511256              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
Note: See TracChangeset for help on using the changeset viewer.