/** * \file Aeolus-SkillScore_q.cpp * * $Id$ * * Cost function for parameter optimisation for specific humidity of Aeolus when coupled to FMS/POEM. * */ #define DEBUG_SKILLSCORE static HorizontalVariable ** q_1000, /** q_900,*/ ** q_500; ///< specific humidity at 1000, 900 and 500 mb (ca. sealevel, 1km and 5km height) static int *ntsteps; #ifdef DEBUG_SKILLSCORE static NetCDFOutput *scorecdf; #endif static void init_Aeolus_Score(//, const int ntsteps ///< select between monthly, seasonal, annual evaluation ) { const char *monthabbr[12] = { "jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec" }; // allocate arrays for temporal aggregation // 12 monthly means q_1000 = new HorizontalVariable *[12]; //q_900 = new HorizontalVariable *[12]; q_500 = new HorizontalVariable *[12]; ntsteps = new int[12]; for (int i = 0; i < 12; i++) { q_1000[i] = new HorizontalVariable(*grid, SF, ((string)("q_1000_")+monthabbr[i]).c_str(), "kg/kg", "", Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, 0., -1. ); //q_900[i] = new HorizontalVariable(*grid, 1, ((string)("q_900_")+monthabbr[i]).c_str(), "kg/kg", // "", // Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, // 0., -1. ); q_500[i] = new HorizontalVariable(*grid, EBL, ((string)("q_500_")+monthabbr[i]).c_str(), "kg/kg", "", Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, 0., -1. ); q_1000[i]->UniformValue(0.0); //q_900[i]->UniformValue(0.0); q_500[i]->UniformValue(0.0); ntsteps[i] = 0; } #ifdef DEBUG_SKILLSCORE // make time units and calendar compatible with FMS diag manager output scorecdf = new NetCDFOutput(*grid, "debug_skillscore", false, "$Id$", "days since 0001-01-01 00:00:00", "NOLEAP"); for (int i = 0; i < 12; i++) { scorecdf->AddVariable(*(q_1000[i])); scorecdf->AddVariable(*(q_500[i])); } #endif } /// Aggregate to monthly mean values. /// Start at 1-Jan of the 2nd year of simulation. /// With fixed-surface-condition runs, we can assume that already in the 2nd year we have /// reached a steady state. /// \param days_since_0 is days since start of time axis, not day of year! /// \param seconds_of_day is seconds of day static void Aeolus_Score_Add_Tstep(const int year, const int month, const int days_since_0, const int seconds_of_day) { static int prevyear = -1; static int agg_year = 0; cout << __FUNCTION__ << " year " << year << " month " << month << endl; if (-1 != prevyear && year != prevyear) agg_year++; // we are at least in the 2nd year of simulation cout << " prevyear " << prevyear << " aggyear " << agg_year << endl; if (1 == agg_year) { static HorizontalVariable tmpl1(*grid, 0, "tmpl1", "", ""); static HorizontalVariable tmpl3(*grid, 2, "tmpl3", "", ""); cout << __FUNCTION__ << " aggregating " << (q_1000[month-1])->Longname() << endl; //cout << "Types q_1000 " << (q_1000[month-1])->Type() << " q_500 " << (q_500[month-1])->Type() // << " qS " << humidity->qS_().Type() // << " tmpl1 " << tmpl1.Type() << " tmpl3 " << tmpl3.Type() << endl; *(q_1000[month-1]) += humidity->qS_(); //tmpl1.ImportFrom(humidity->q_()); //*(q_900[month-1]) += tmpl1; tmpl3.ImportFrom(humidity->q_()); *(q_500[month-1]) += tmpl3; ntsteps[month-1]++; #ifdef DEBUG_SKILLSCORE scorecdf->OutputToFile(days_since_0+seconds_of_day/(24.0*3600.0)); // approximate date #endif } prevyear = year; } // from Sonja: ///p/projects/climber3/molnos/Taylor_step1_20_05_2016_el_nino_la_nina_new/main_Aeolus_DynamicalCore_TaylorSkillScore.cpp::taylor_skillscore*() static double taylor_skillscore(const HorizontalVariable& xm, ///< model data const HorizontalVariable& xo, ///< observation data const HorizontalVariable& weights ///< weight function ) { assert(xm.Size() == xo.Size()); assert(xm.Size() == weights.Size()); cout << __FUNCTION__ << "Types xm " << xm.Type() << " xo " << xo.Type() << endl << flush; static HorizontalVariable tmp(*grid, 0, "tmp", "", ""); tmp.SetType(weights.Type(), *grid); assert(tmp.Size() == weights.Size()); // weighted Mean values double weights_sum = weights.sum(); tmp = weights; tmp.SetType(xm.Type(), *grid); cout << __FUNCTION__ << "Types xm " << xm.Type() << " tmp " << tmp.Type() << endl << flush; tmp *= xm; double mean_model=tmp.Mean() / weights_sum; tmp.SetType(weights.Type(), *grid); tmp = weights; tmp.SetType(xo.Type(), *grid); cout << __FUNCTION__ << "Types xo " << xo.Type() << " tmp " << tmp.Type() << endl << flush; tmp *= xo; double mean_obs = tmp.Mean() / weights_sum; //Spatial statistical moments double std_model = 0., std_obs = 0., crossCorrelation = 0.; for (unsigned int i = 0; i < xm.Size(); i++) { std_model += weights[i]*(xm[i]-mean_model)*(xm[i]-mean_model); std_obs += weights[i]*(xo[i]-mean_obs)*(xo[i]-mean_obs); crossCorrelation += weights[i]*(xm[i]-mean_model)*(xo[i]-mean_obs); } double sigma_model=sqrt(std_model/weights_sum); double sigma_obs=sqrt(std_obs/weights_sum); //ratio between standard deviation of model and of obs -> 1 ==perfect // <1 == std of model smaller than std of obs // >1 == std of model greater than std of obs double amplitude=sigma_model/sigma_obs; //ratio of the spatial correlation -> mean between obs and model double ratio_spatialCorrelation=crossCorrelation/(weights_sum*sigma_model*sigma_obs); //finally calculate Taylor Skillscore return 0.25*pow(1.0+ratio_spatialCorrelation,4)/((amplitude+1.0/amplitude)*(amplitude+1.0/amplitude)); } static double Aeolus_Score_Calc() { /* first finish aggregation */ for (int i = 0; i < 12; i++) { *(q_1000[i]) /= ntsteps[i]; /*(q_900[i]) /= ntsteps[i];*/ *(q_500[i]) /= ntsteps[i]; } // Init weights // Weight is cell area, between 60S and 70N, else 0 // For specific humidity, the equator region is important. HorizontalVariable weights2d(*grid, SF, "weight2d"," ", "area weights for simple skillscore", 0.0, 1.0); double wsum = 0.0; for (vector::const_iterator it = grid->SurfaceCellsBegin(); it < grid->SurfaceCellsEnd(); it++){ unsigned int id = it->ID(); double lat = it->lat(); double w = 0; if (lat>-60 && lat<70) { w = it->FirstConnectedZInterface()->Area(); } weights2d.Write(id,w);//Level=3 wsum+=w; } weights2d /= wsum; // // Prepare to read observation data (ERA interim reanalysis in this case) // HorizontalVariable q_1000_obs(*grid, SF, "Q","kg kg**-1", "q 1000mbar observation data",Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, 0, 100); //HorizontalVariable q_900_obs(*grid, 1, "Q","kg kg**-1", "q 900mbar observation data",Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, 0, 100); HorizontalVariable q_500_obs(*grid, 3, "Q","kg kg**-1", "q 500mbar observation data",Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, 0, 100); // Sonja does 3 series of experiments: all months, el nino months, la nina months. // Here we do only the all-months experiment char* q_1000_file = "q-ymonmean-1000-90S-90N-1979-1999.nc"; //char* q_900_file = "q-ymonmean-900-90S-90N-1979-1999.nc"; char* q_500_file = "q-ymonmean-500-90S-90N-1979-1999.nc"; //load q at different levels NetCDFInput input_q_1000(*grid, q_1000_file); //NetCDFInput input_q_900(*grid, q_900_file); NetCDFInput input_q_500(*grid, q_500_file); input_q_1000.GetVariable(q_1000_obs); //input_q_900.GetVariable(q_900_obs); input_q_500.GetVariable(q_500_obs); valarray skillscore_vector_q_1000(12), /*skillscore_vector_q_900[12],*/ skillscore_vector_q_500(12); // skill score is computed for each timestep, i.e. monthly mean, separately. for (int month = 0; month < 12; month++) { input_q_1000.InputFromFile(month); //input_q_900.InputFromFile(month); input_q_500.InputFromFile(month); skillscore_vector_q_1000[month] = taylor_skillscore(*(q_1000[month]), q_1000_obs, weights2d); //skillscore_vector_q_900[month] = taylor_skillscore(*(q_900[month]), q_900_obs, weights2d); skillscore_vector_q_500[month] = taylor_skillscore(*(q_500[month]), q_500_obs, weights2d); } //skillscore is mean value of vector skillscore_vector values[t] double skillscore_q_1000 = skillscore_vector_q_1000.sum()/skillscore_vector_q_1000.size(); //double skillscore_q_900 = skillscore_vector_q_900.sum()/skillscore_vector_q_900.size(); double skillscore_q_500 = skillscore_vector_q_500.sum()/skillscore_vector_q_500.size(); cout<<"skillscore_q_1000; "<100000) skillscore = 100000; return skillscore; } /// Clean up. /// \param days is number of days since begin of simulation static void Aeolus_Score_End(double days) { #ifdef DEBUG_SKILLSCORE scorecdf->OutputToFile(days); // approximate date delete scorecdf; #endif // deallocate storage for temporal aggregation for (int i = 0; i < 12; i++) { delete q_1000[i]; //delete q_900[i]; delete q_500[i]; } delete[] q_1000; //delete[] q_900; delete[] q_500; delete[] ntsteps; }