Changeset 151 for altifloat/src

Timestamp:

07/02/15 10:08:35 (9 years ago)

Author:

jbrlod

Message:

load and interp wind data

Location:

altifloat/src

Files:

• floater.h (modified) (6 diffs)
• floater_delta.d (modified) (2 diffs)

altifloat/src/floater.h

@@ -1 +1 @@
 #define UPDATE_BCK
+
+#include <netcdf.h>
 
 //constantes
@@ -21 +23 @@
 YREAL umod[nlon][nlat][jptfl];
 YREAL vmod[nlon][nlat][jptfl];
+YREAL uwind[nlon][nlat][jptfl];
+YREAL vwind[nlon][nlat][jptfl];
 
 //Background
@@ -47 +51 @@
 #endif
 #endif
+
+#define ERR(e) {fprintf(stderr,"Error: %s\n", nc_strerror(e)); exit(EXIT_FAILURE);}
+#define LONDIM "lon_interp"
+#define LATDIM "lat_interp"
+#define UVAR   "u10_f"
+#define VVAR   "v10_f"
 
 //Options du run incr士ental
@@ -61 +71 @@
 void run_inc(int argv, char *argc[]);
 void div_init();
+void init_wind(char fname[],int iter0,int *ptime0);
 
 void appli_start(int argc, char *argv[]){
@@ -195 +206 @@
 
 }
+
+void load_wind(int argv, char *argc[]) {
+  /* load_wind windname iter0 time0
+     load a wind netcdfile ( from ERAI)
+     time0 is the time of the first iteration considered in the file
+     (by default time0 is time(0) field in the netcdf file)
+     iter0 is the timestep corresponding to time0.
+     (by default iter0 is 0)
+   */
+  if ((argv<2) | (argv >4)) {
+    fprintf(stderr,"(load_wind): wrong number of arguments\n");
+    return;
+  }
+
+  //time0 is a pointer to pass the NULL address in case time0 is
+  //not specified in the call of the function
+  int *time0=NULL;
+
+  int iter0=0;
+  char fname[50];
+  sprintf(fname,argc[1]);
+  if (argv > 2) {
+    iter0 = atoi(argc[2]);
+    if (argv > 3) 
+      *time0 = atoi(argc[3]);
+  }
+  init_wind(fname,iter0,time0);
+}
+
+void init_wind(char fname[],int iter0,int *ptime0) {
+  int ncid,retval;
+  int dlatid,dlonid,dtimeid;
+  int latid,lonid,timeid,u10id,v10id;
+  size_t Nlon,Nlat,Ntime;
+  if ((retval = nc_open(fname,NC_NOWRITE, &ncid)))
+    ERR(retval);
+  if ((retval = nc_inq_dimid(ncid,LONDIM,&dlonid)))
+    ERR(retval);
+  if ((retval = nc_inq_dimid(ncid,LATDIM,&dlatid)))
+    ERR(retval);
+  if ((retval = nc_inq_dimid(ncid,"time",&dtimeid)))
+    ERR(retval);
+  if ((retval = nc_inq_dimlen(ncid,dlonid,&Nlon)))
+    ERR(retval);
+  if ((retval = nc_inq_dimlen(ncid,dlatid,&Nlat)))
+    ERR(retval);
+  if ((retval = nc_inq_dimlen(ncid,dtimeid,&Ntime)))
+    ERR(retval);
+  
+  if ((Nlon != nlon) | (Nlat != nlat)) {
+   fprintf(stderr,
+           "(init_wind) Wrong dimension (nlat=%d,nlon=%d) for the wind_field %s\n",
+           (int)Nlat, (int)Nlon, fname);
+  }
+  if ((retval = nc_inq_varid(ncid,LATDIM,&latid)))
+    ERR(retval);
+  if ((retval = nc_inq_varid(ncid,LONDIM,&lonid)))
+    ERR(retval);
+  if ((retval = nc_inq_varid(ncid,"time",&timeid)))
+    ERR(retval);
+  if ((retval = nc_inq_varid(ncid,UVAR,&u10id)))
+    ERR(retval);
+  if ((retval = nc_inq_varid(ncid,VVAR,&v10id)))
+    ERR(retval);
+  
+  int time[Ntime];
+  float u10[Ntime][nlat][nlon];
+  float v10[Ntime][nlat][nlon];
+  
+  if ((retval = nc_get_var_int(ncid,timeid,&time[0])))
+    ERR(retval);
+  
+  if ((retval = nc_get_var_float(ncid,u10id,&u10[0][0][0])))
+    ERR(retval);
+  
+  if ((retval = nc_get_var_float(ncid,v10id,&v10[0][0][0])))
+    ERR(retval);
+
+  //Load time dimensions
+  if (ptime0 == NULL)
+    {
+      ptime0 = (int *)malloc(sizeof(int));
+      *ptime0 = time[0];
+    }
+
+  //Copy arrays
+  int nt=0;
+  int realtime=*ptime0; //time value from time (in hour)
+  int timeindex=iter0; //time index for uwind and vwind
+  int it;
+  YREAL alpha; //interp coefficient
+  for (it=0;it<(int)Ntime-1;it++) {
+    
+// ionterp whil relatime is between 2 step time
+    while (time[it]<=realtime && time[it+1]>realtime) {
+      
+
+      if (timeindex >= jptfl) {
+        fprintf(stderr,"(init_wind) Warning: time values in %s don't fill the jptfl\n",fname);
+        return;
+      }
+      alpha = ((double)realtime - time[it])/(time[it+1]-time[it]);
+      for (int ilon=0;ilon<nlon;ilon++)
+        for (int ilat=0;ilat<nlat;ilat++) {
+          uwind[ilon][ilat][timeindex] = (1-alpha)*u10[it][ilat][ilon] + alpha*u10[it+1][ilat][ilon];
+          vwind[ilon][ilat][timeindex] = (1-alpha)*v10[it][ilat][ilon] + alpha*v10[it+1][ilat][ilon];
+         
+        }
+      //To avoid cumulative rounding mistake, realtime is calculate with :
+      nt++;
+      realtime = *ptime0 + nt*(double)rdtflo/3600;
+      timeindex ++;
+
+    } //while time[it]<=realtime
+  } //for it
+  if (time[it]==realtime && timeindex < jptfl  ) //last step time
+    {
+      for (int ilon=0;ilon<nlon;ilon++)
+        for (int ilat=0;ilat<nlat;ilat++) {
+          uwind[ilon][ilat][timeindex] = u10[it][ilat][ilon] ;
+          vwind[ilon][ilat][timeindex] = v10[it][ilat][ilon] ;
+          
+        }
+    }
+  
+}
+
+
 
 void read_mask(int argv, char *argc[]) {
@@ -446 +585 @@
       fprintf(fid,"%d %d %f %f\n",it,j,piret[it][j],pjret[it][j]);
         
+                                             
+  fclose(fid);
+
+}
+
+
+void save_output_wind (int argc, char *argv[]) {
+  //Function used only for test
+FILE *fid;
+  fid=fopen(argv[1],"w");
+  if (fid==NULL) {
+    printf("\nfailed to open %s",argv[1]);
+    exit(3);
+  }
+  int j,k,it;
+    for (it=0;it<jptfl;it++)
+      for (j=0;j<nlon;j++)
+        for (k=0;k<nlat;k++) {
+          fprintf(fid,"%d %d %d %f %f\n",it,j,k,uwind[j][k][it],vwind[j][k][it]);
+        }
                                              
   fclose(fid);

altifloat/src/floater_delta.d

@@ -82 +82 @@
 modul u       space S_euler       noward    output 1 tempo //noward veut dire pas de calcul
 modul v       space S_euler       noward    output 1 tempo
+modul uwind   space S_euler       noward    output 1 tempo
+modul vwind   space S_euler       noward    output 1 tempo
 modul locate  space S_eulerlocate input  2  output 1 tempo
 modul mask    space S_euler       noward    output 1
@@ -280 +282 @@
 insert_fct multirun
 insert_fct arg save_output_uv
+insert_fct arg save_output_wind
 insert_fct arg save_output_rfloat
 insert_fct arg read_bck