Changeset 19566


Ignore:
Timestamp:
09/19/15 14:54:07 (10 years ago)
Author:
Eric.Larour
Message:

CHG: added verbose of each module + fixed some serious bugs.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r19554 r19566  
    129129
    130130} /*}}}*/
    131 void grainGrowth(IssmDouble* re, IssmDouble* gdn, IssmDouble* gsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx){ /*{{{*/
     131void grainGrowth(IssmDouble* re, IssmDouble* gdn, IssmDouble* gsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx,int sid){ /*{{{*/
    132132
    133133        /*Created by: Alex S. Gardner, University of Alberta
    … …  
    173173        IssmDouble  Q=0;
    174174
     175        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   grain growth module\n");
     176       
    175177        /*only when aIdx = 1 or 2 do we run grainGrowth: */
    176178        if(aIdx!=1 & aIdx!=2){
    … …  
    302304
    303305}  /*}}}*/
    304 void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m) { /*{{{*/
     306void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m,int sid) { /*{{{*/
    305307
    306308        //// Calculates Snow, firn and ice albedo as a function of:
    … …  
    345347        // a = albedo(4, [], [], [], 0.48, 0.85, [0.8 0.5 ... 0.48], ...
    346348        //   [273 272.5 ... 265], [0 0.001 ... 0], 0, 0.01, 15, 15, 7, 3600)
     349
     350        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   albedo module\n");
    347351       
    348352        //some constants:
    … …  
    451455        else if (isnan(a[0])) _error_("albedo == NAN\n");
    452456}  /*}}}*/
    453 void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz) { /*{{{*/
     457void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz,int sid) { /*{{{*/
    454458
    455459        /* ENGLACIAL THERMODYNAMICS*/
    … …  
    527531        IssmDouble EC;
    528532
     533        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
     534
    529535        // INITIALIZE
    530536        CI = 2102;          // heat capacity of snow/ice (J kg-1 k-1)
    … …  
    621627        // must go evenly into one hour or the data frequency if it is smaller
    622628
    623         // all integer factors of the number of second in a day (8600 [s])
     629        // all integer factors of the number of second in a day (86400 [s])
    624630        int f[45] = {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 30, 36, 40, 45, 48, 50, 60,
    625631    72, 75, 80, 90, 100, 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 450, 600, 720, 900, 1200, 1800, 3600};
    626632
    627633        // return the min integer factor that is < dt
    628         max_fdt=0;
     634        max_fdt=f[0];
    629635        for(int i=0;i<45;i++){
    630636                if (f[i]<dt)if(f[i]>=max_fdt)max_fdt=f[i];
    … …  
    773779   
    774780                //SW penetrates surface
    775                 for(int i=0;i<m;i++) T[i] = T[i] + dT_sw[i];
     781                for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
    776782                T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
    777783               
    778784                // temperature diffusion
    779                 for(int i=0;i<m;i++)T0[1+i]=T[i];
    780                 for(int i=0;i<m;i++) Tu[i] = T0[i];
    781                 for(int i=0;i<m;i++) Td[i] = T0[2+i];
    782    
    783                 for(int i=0;i<m;i++) T[i] = (Np[i] * T[i]) + (Nu[i] * Tu[i]) + (Nd[i] * Td[i]);
     785                for(int j=0;j<m;j++)T0[1+j]=T[j];
     786                for(int j=0;j<m;j++) Tu[j] = T0[j];
     787                for(int j=0;j<m;j++) Td[j] = T0[2+j];
     788                for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
    784789               
    785790                // calculate cumulative evaporation (+)/condensation(-)
    … …  
    826831
    827832}  /*}}}*/
    828 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m){ /*{{{*/
     833void shortwave(IssmDouble* swf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
    829834
    830835        // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
    … …  
    848853        //   swf     = absorbed shortwave radiation [W m-2]
    849854
    850         /*outputs:*/
    851         IssmDouble* swf = NULL;
     855        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   shortwave module\n");
    852856
    853857        // SHORTWAVE FUNCTION
    854         // initialize variables
    855         swf = xNewZeroInit<IssmDouble>(m);
    856 
    857858        if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
    858859       
    … …  
    10101011
    10111012                        // add flux absorbed at surface
    1012                         swf[0] = swf[0] + swf_s;
     1013                        swf[0] += swf_s;
    10131014
    10141015                        /*Free ressources:*/
    … …  
    10191020                }
    10201021        }
    1021         /*Assign output pointers: */
    1022         *pswf=swf;
    10231022} /*}}}*/
    1024 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow){ /*{{{*/
     1023void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid){ /*{{{*/
    10251024
    10261025        // Adds precipitation and deposition to the model grid
    … …  
    10621061        IssmDouble* gsp=NULL;
    10631062        int         m;
     1063
     1064        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   accumulation module\n");
    10641065
    10651066        /*Recover pointers: */
    … …  
    11691170        *pm=m;
    11701171} /*}}}*/
    1171 void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin){ /*{{{*/
     1172void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, int sid){ /*{{{*/
    11721173
    11731174        //// MELT ROUTINE
    … …  
    12201221        int         n=*pn;
    12211222       
     1223        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
     1224
    12221225        //// INITIALIZATION
    12231226
    … …  
    15641567
    15651568} /*}}}*/
    1566 void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,int m){ /*{{{*/
     1569void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble dIce, int m, int sid){ /*{{{*/
    15671570
    15681571        //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???]
    … …  
    16061609        //// MAIN FUNCTION
    16071610        // specify constants
    1608         const IssmDouble dIce    = 910;         // density of ice [kg m-3]
    16091611        dt      = dt / 86400;  // convert from [s] to [d]
    16101612        // R     = 8.314        // gas constant [mol-1 K-1]
    … …  
    16151617        /*intermediary: */
    16161618        IssmDouble c0,c1,H;
     1619       
     1620        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
    16171621
    16181622        // initial mass
    … …  
    16791683                dz[i] = mass_init[i] / d[i];
    16801684        }
     1685
    16811686        /*Free ressources:*/
    16821687        xDelete<IssmDouble>(mass_init);
    … …  
    16851690
    16861691} /*}}}*/
    1687 void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz){ /*{{{*/
     1692void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, int sid){ /*{{{*/
    16881693
    16891694        //// TURBULENT HEAT FLUX
    … …  
    17251730        /*output: */
    17261731        IssmDouble shf, lhf, EC;
     1732       
     1733        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
    17271734
    17281735        // calculated air density [kg/m3]
Note: See TracChangeset for help on using the changeset viewer.