← Back to branch summary

~reducedmodelling/fluidity/ROM_Non-intrusive-ann

« back to all changes in this revision

Viewing changes to climatology/create_climatology_atlas.cpp

  • Committer: tmb1
  • Date: 2010-11-09 14:08:19 UTC
  • Revision ID: svn-v4:5bf5533e-7014-46e3-b1bb-cce4b9d03719:trunk:2423
Changing isntances of Chris' person email address to the group
software email address, and removing some rogue tabs at the same
time to allow me to commit the changes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
climatology/create_climatology_atlas.cpp
8
8
    Department of Earth Science and Engineering
9
9
    Imperial College London
10
10
 
11
 
    C.Pain@Imperial.ac.uk
 
11
    C.Pain@imperial.ac.uk
12
12
    
13
13
    This library is free software; you can redistribute it and/or
14
14
    modify it under the terms of the GNU Lesser General Public
185
185
  for(size_t k=0;k<nlevels;k++){
186
186
    for(size_t j=0; j<(size_t)jdim0; j++){
187
187
      for(size_t i=0; i<(size_t)idim0; i++){
188
 
                if(fabs(variable[0][k][j*idim0+i]-land_val)>1.0){
189
 
                  for(size_t m=0;m<variable.size();m++)
190
 
                        d[m] = variable[m][k][j*idim0+i];
191
 
                  
192
 
                  kilworth_correction(A, d);
193
 
                  
194
 
                  for(size_t m=0;m<variable.size();m++){
195
 
                        if(d[m]>1000.0){
196
 
                          cout<<"A="<<endl;
197
 
                          for(size_t ii=0;ii<(size_t)n;ii++){
198
 
                                for(size_t jj=0;jj<(size_t)n;jj++)
199
 
                                  cout<<A[ii*n+jj]<<" ";
200
 
                                cout<<endl;
201
 
                          }
202
 
                          cout<<"d = ";
203
 
                          for(size_t m=0;m<variable.size();m++)
204
 
                                cout<<variable[m][k][j*idim0+i]<<" ";
205
 
                          cout<<endl;
206
 
                          cout<<"d' = ";
207
 
                          for(size_t m=0;m<variable.size();m++)
208
 
                                cout<<d[m]<<" ";
209
 
                          cout<<endl;
210
 
                        }
211
 
                        max_v = max(max_v, (float)d[m]);
212
 
                        min_v = min(min_v, (float)d[m]);
213
 
                        variable[m][k][j*idim0+i] = d[m];
214
 
                  }
215
 
                }
 
188
    if(fabs(variable[0][k][j*idim0+i]-land_val)>1.0){
 
189
      for(size_t m=0;m<variable.size();m++)
 
190
        d[m] = variable[m][k][j*idim0+i];
 
191
      
 
192
      kilworth_correction(A, d);
 
193
      
 
194
      for(size_t m=0;m<variable.size();m++){
 
195
      if(d[m]>1000.0){
 
196
        cout<<"A="<<endl;
 
197
        for(size_t ii=0;ii<(size_t)n;ii++){
 
198
        for(size_t jj=0;jj<(size_t)n;jj++)
 
199
          cout<<A[ii*n+jj]<<" ";
 
200
        cout<<endl;
 
201
        }
 
202
        cout<<"d = ";
 
203
        for(size_t m=0;m<variable.size();m++)
 
204
        cout<<variable[m][k][j*idim0+i]<<" ";
 
205
        cout<<endl;
 
206
        cout<<"d' = ";
 
207
        for(size_t m=0;m<variable.size();m++)
 
208
        cout<<d[m]<<" ";
 
209
        cout<<endl;
 
210
      }
 
211
      max_v = max(max_v, (float)d[m]);
 
212
      min_v = min(min_v, (float)d[m]);
 
213
      variable[m][k][j*idim0+i] = d[m];
 
214
      }
 
215
    }
216
216
      }
217
217
    }
218
218
  }
223
223
int diffuse_boundaries(deque< deque< vector<float> > > &dat, float minv, int iterations){
224
224
  // Iterate through time
225
225
  for(deque< deque< vector<float> > >::iterator tdat=dat.begin(); tdat!=dat.end(); tdat++){
226
 
        // Iterate through levels
 
226
  // Iterate through levels
227
227
    for(size_t k=0;k<tdat->size();k++){
228
 
          
229
 
          // Find which points are null
230
 
          vector<bool> null_points(idim0*jdim0, false);
231
 
          for(size_t j=0; j<(size_t)jdim0; j++){
232
 
                for(size_t i=0; i<(size_t)idim0; i++){
 
228
    
 
229
    // Find which points are null
 
230
    vector<bool> null_points(idim0*jdim0, false);
 
231
    for(size_t j=0; j<(size_t)jdim0; j++){
 
232
    for(size_t i=0; i<(size_t)idim0; i++){
233
233
          if(fabs((*tdat)[k][j*idim0+i]-land_val)<1.0){
234
 
                        null_points[j*idim0+i] = true;
235
 
                        (*tdat)[k][j*idim0+i] = minv;
 
234
      null_points[j*idim0+i] = true;
 
235
      (*tdat)[k][j*idim0+i] = minv;
236
236
          }
237
 
                }
238
 
      }
239
 
          
240
 
          // Start diffusing
241
 
          for(size_t its=0;its<iterations;its++){
242
 
                for(size_t j=0; j<(size_t)jdim0; j++){ 
243
 
                  for(size_t i=0; i<(size_t)idim0; i++){
244
 
                        if(null_points[j*idim0+i]){
245
 
                          float sum=(*tdat)[k][j*idim0+i];
246
 
                          int npts=1;
247
 
                          
248
 
                          if(j>0){
249
 
                                sum+=(*tdat)[k][(j-1)*idim0+i];
250
 
                                npts++;
251
 
                          }
252
 
 
253
 
                          if(j+1<jdim0){
254
 
                                sum+=(*tdat)[k][(j+1)*idim0+i];
255
 
                                npts++;
256
 
                          }
257
 
 
258
 
                          sum+=(*tdat)[k][j*idim0+(i+idim0-1)%idim0];
259
 
                          npts++;
260
 
 
261
 
                          sum+=(*tdat)[k][j*idim0+(i+1)%idim0];
262
 
                          npts++;
263
 
                          
264
 
                          if(k>0){
265
 
                                sum+=(*tdat)[k-1][j*idim0+i];
266
 
                                npts++;
267
 
                          }
268
 
 
269
 
                          (*tdat)[k][j*idim0+i] = sum/npts;
270
 
                        }
271
 
                  }
272
 
                }
273
 
          }
274
 
        }
 
237
    }
 
238
      }
 
239
    
 
240
    // Start diffusing
 
241
    for(size_t its=0;its<iterations;its++){
 
242
    for(size_t j=0; j<(size_t)jdim0; j++){ 
 
243
       for(size_t i=0; i<(size_t)idim0; i++){
 
244
      if(null_points[j*idim0+i]){
 
245
        float sum=(*tdat)[k][j*idim0+i];
 
246
        int npts=1;
 
247
        
 
248
        if(j>0){
 
249
        sum+=(*tdat)[k][(j-1)*idim0+i];
 
250
        npts++;
 
251
        }
 
252
 
 
253
         if(j+1<jdim0){
 
254
        sum+=(*tdat)[k][(j+1)*idim0+i];
 
255
        npts++;
 
256
        }
 
257
 
 
258
        sum+=(*tdat)[k][j*idim0+(i+idim0-1)%idim0];
 
259
        npts++;
 
260
 
 
261
        sum+=(*tdat)[k][j*idim0+(i+1)%idim0];
 
262
        npts++;
 
263
        
 
264
        if(k>0){
 
265
        sum+=(*tdat)[k-1][j*idim0+i];
 
266
        npts++;
 
267
        }
 
268
 
 
269
        (*tdat)[k][j*idim0+i] = sum/npts;
 
270
      }
 
271
      }
 
272
    }
 
273
    }
 
274
  }
275
275
  }
276
276
 
277
277
  return 0;
278
278
}
279
279
 
280
280
int write_netcdf_variable(int ncid, const int* dimids, const char* name, const char* long_name, const char *units,
281
 
                          deque< deque< vector<float> > >& variable){
 
281
        deque< deque< vector<float> > >& variable){
282
282
#ifdef HAVE_NETCDF
283
283
 
284
284
  float max_v, min_v;
320
320
    for(size_t k=0;k<variable[m].size();k++){
321
321
      start[1] = k;
322
322
      for(size_t j=0; j<(size_t)jdim0; j++)
323
 
                for(size_t i=0; i<(size_t)idim0; i++){
324
 
                  // After applying a diffusion filter we need to ensure the
325
 
                  // values are still within bounds.
326
 
                  float v = min(max_v, variable[m][k][j*idim0+i]);
327
 
                  v = max(v, min_v);
328
 
                  
329
 
                  frame[j*idim0+i] = lround((v - add_offset)/scale_factor);
330
 
                }
 
323
    for(size_t i=0; i<(size_t)idim0; i++){
 
324
      // After applying a diffusion filter we need to ensure the
 
325
      // values are still within bounds.
 
326
      float v = min(max_v, variable[m][k][j*idim0+i]);
 
327
      v = max(v, min_v);
 
328
      
 
329
      frame[j*idim0+i] = lround((v - add_offset)/scale_factor);
 
330
    }
331
331
      ncvarput(ncid, varid, start, count, &(frame[0]));
332
332
    }
333
333
  }
460
460
    for(size_t k=0; k<24; k++){
461
461
      // cout<<"Reading level "<<k<<endl;
462
462
      for(size_t j=0; j<(size_t)jdim0; j++)
463
 
        for(size_t i=0; i<(size_t)idim0; i++){
464
 
          fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
465
 
        }
 
463
  for(size_t i=0; i<(size_t)idim0; i++){
 
464
    fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
 
465
  }
466
466
    }
467
467
    fclose(fp);
468
468
  }
489
489
    for(size_t k=0; k<24; k++){
490
490
      // cout<<"Reading level "<<k<<endl;
491
491
      for(size_t j=0; j<(size_t)jdim0; j++)
492
 
        for(size_t i=0; i<(size_t)idim0; i++){
493
 
          fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
494
 
        }
 
492
  for(size_t i=0; i<(size_t)idim0; i++){
 
493
    fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
 
494
  }
495
495
    }
496
496
    fclose(fp);
497
497
  }
526
526
    for(size_t k=0; k<33; k++){
527
527
      // cout<<"Reading level "<<k<<endl;
528
528
      for(size_t j=0; j<(size_t)jdim0; j++)
529
 
        for(size_t i=0; i<(size_t)idim0; i++){
530
 
          float v;
531
 
          fscanf(fp, "%f", &v);
532
 
          if(k>=24){
533
 
            dat[m][k-24][j*idim0+i] = v;
534
 
          }
535
 
        }
 
529
  for(size_t i=0; i<(size_t)idim0; i++){
 
530
    float v;
 
531
    fscanf(fp, "%f", &v);
 
532
    if(k>=24){
 
533
      dat[m][k-24][j*idim0+i] = v;
 
534
    }
 
535
  }
536
536
    }
537
537
    fclose(fp);
538
538
  }
558
558
    for(size_t k=0; k<33; k++){
559
559
      // cout<<"Reading level "<<k<<endl;
560
560
      for(size_t j=0; j<(size_t)jdim0; j++)
561
 
        for(size_t i=0; i<(size_t)idim0; i++){
562
 
          float v;
563
 
          fscanf(fp, "%f", &v);
564
 
          if(k>=24){
565
 
            dat[m][k-24][j*idim0+i] = v;
566
 
          }
567
 
        }
 
561
  for(size_t i=0; i<(size_t)idim0; i++){
 
562
    float v;
 
563
    fscanf(fp, "%f", &v);
 
564
    if(k>=24){
 
565
      dat[m][k-24][j*idim0+i] = v;
 
566
    }
 
567
  }
568
568
    }
569
569
    fclose(fp);
570
570
  }