external/io_grib_share/get_region_center.c

   1 #include "get_region_center.h"
   2 #include "gridnav.h"
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6
   7 #include "wrf_projection.h"
   8
   9 int get_gridnav_projection(int wrf_projection);
  10
  11 /****************************************************************************
  12  *
  13  * This function calculates the center lat and lon of a region within a larger
  14  *   domain.  It is useful for calculating the center of the boundary regions
  15  *   in a domain.
  16  *
  17  ****************************************************************************/
  18
  19
  20 int GET_REGION_CENTER(char *MemoryOrderIn, int *projection,
  21                       float *domain_center_lat,
  22                       float *domain_center_lon, int *full_xsize,
  23                       int *full_ysize, float *dx, float *dy,
  24                       float *proj_central_lon,
  25                       int *proj_center_flag, float *truelat1,
  26                       float *truelat2, int *region_xsize, int *region_ysize,
  27                       float *region_center_lat, float *region_center_lon,
  28                       int strlen1)
  29 {
  30
  31   char *MemoryOrder;
  32   int grid_projection;
  33   float full_xcenter, full_ycenter;
  34   float region_xcenter, region_ycenter;
  35   int status;
  36   int orig;
  37   int x_pos, y_pos;
  38   GridNav gridnav;
  39
  40   MemoryOrder = (char *)malloc((strlen1+1)*sizeof(char));
  41   memcpy(MemoryOrder,MemoryOrderIn,strlen1);
  42   MemoryOrder[strlen1] = '\0';
  43
  44   grid_projection = get_gridnav_projection(*projection);
  45
  46   full_xcenter = (*full_xsize - 1) / 2.;
  47   full_ycenter = (*full_ysize - 1) / 2.;
  48   region_xcenter = (*region_xsize - 1) / 2.;
  49   region_ycenter = (*region_ysize - 1) / 2.;
  50
  51   orig = 0;
  52
  53   if (strncmp(MemoryOrder,"XS", 2) == 0)
  54     {
  55       x_pos = region_xcenter;
  56       y_pos = full_ycenter;
  57     }
  58   else if (strncmp(MemoryOrder,"XE", 2) == 0)
  59     {
  60       x_pos = (*full_xsize - 1) - region_xcenter;
  61       y_pos = full_ycenter;
  62     }
  63   else if (strncmp(MemoryOrder,"YS", 2) == 0)
  64     {
  65       x_pos = full_xcenter;
  66       y_pos = region_ycenter;
  67     }
  68   else if (strncmp(MemoryOrder,"YE", 2) == 0)
  69     {
  70       x_pos = full_xcenter;
  71       y_pos = (*full_ysize - 1) - region_ycenter;
  72     }
  73   else
  74     {
  75       orig = 1;
  76     }
  77
  78   if (orig == 1)
  79     {
  80       *region_center_lat = *domain_center_lat;
  81       *region_center_lon = *domain_center_lon;
  82       status = 0;
  83     }
  84   else
  85     {
  86       /* Initialize grid structure */
  87       /*
  88       status = GRID_init(grid_info->center_lat, grid_info->central_lon,
  89                          grid_projection,
  90                          grid_info->latin1, grid_info->latin2,
  91                          grid_info->xpoints, grid_info->ypoints,
  92                          grid_info->Di, grid_info->Dj,
  93                          grid_info->center_lat, grid_info->center_lon,
  94                          x_center, y_center,
  95                          &gridnav);
  96       */
  97       status = GRID_init(*domain_center_lat, *proj_central_lon,
  98                          grid_projection,
  99                          *truelat1, *truelat2,
 100                          *full_xsize, *full_ysize, *dx, *dy,
 101                          *domain_center_lat, *domain_center_lon,
 102                          full_xcenter, full_ycenter,
 103                          &gridnav);
 104       if (!status)
 105         {
 106           fprintf(stderr,"get_region_center: error from GRID_init\n");
 107         }
 108
 109       /* get lat/lon of center of region */
 110       status = GRID_to_latlon(&gridnav, x_pos, y_pos, region_center_lat,
 111                               region_center_lon);
 112       if (!status)
 113         {
 114           fprintf(stderr,
 115                   "get_region_cneter: error from GRID_to_latlon for first lat/lon\n");
 116         }
 117
 118     }
 119
 120   free(MemoryOrder);
 121   return status;
 122
 123 }
 124 /******************************************************************************
 125  * translates the grid projection identifier from the WRF id to the grib id.
 126  *****************************************************************************/
 127
 128 int get_gridnav_projection(int wrf_projection)
 129 {
 130   int gridnav_projection;
 131
 132   /* Set the grid projection in the gridnav units */
 133   switch (wrf_projection)
 134     {
 135     case WRF_LATLON:
 136     case WRF_CASSINI:
 137       gridnav_projection = GRID_LATLON;
 138       break;
 139     case WRF_MERCATOR:
 140       gridnav_projection = GRID_MERCATOR;
 141       break;
 142     case WRF_LAMBERT:
 143       gridnav_projection = GRID_LAMCON;
 144       break;
 145     case WRF_POLAR_STEREO:
 146       gridnav_projection = GRID_POLSTR;
 147       break;
 148     default:
 149       fprintf(stderr,"Error, invalid projection: %d\n",wrf_projection);
 150       gridnav_projection = -1;
 151     }
 152
 153   return gridnav_projection;
 154 }
 155
 156 int GET_LL_LATLON(float *central_lat, float *central_lon, int *projection,
 157                   float *latin1, float *latin2, int *nx, int *ny,
 158                   float *dx, float *dy, float *center_lat, float *center_lon,
 159                   float *LLLa, float *LLLo, float *URLa, float *URLo, int *ierr)
 160 {
 161
 162   int grid_projection;
 163   float x_center;
 164   float y_center;
 165   GridNav gridnav;
 166   int status;
 167
 168   grid_projection = get_gridnav_projection(*projection);
 169
 170   /* Get coords of center of grid */
 171   x_center = (*nx + 1)/2.;
 172   y_center = (*ny + 1)/2.;
 173
 174  /* Initialize grid structure */
 175   status = GRID_init(*central_lat, *central_lon, grid_projection,
 176                      *latin1, *latin2, *nx, *ny, *dx, *dy,
 177                      *center_lat, *center_lon, x_center, y_center,
 178                      &gridnav);
 179   if (!status)
 180     {
 181       fprintf(stderr,"write_grib: error from GRID_init\n");
 182       *ierr = 1;
 183       return(0);
 184     }
 185
 186   /* get lat/lon of lower left corner */
 187   status = GRID_to_latlon(&gridnav, 1, 1, LLLa, LLLo);
 188   if (!status)
 189     {
 190       fprintf(stderr,
 191               "write_grib: error from GRID_to_latlon for first lat/lon\n");
 192       *ierr = 1;
 193       return(0);
 194     }
 195
 196   /* get lat/lon of upper right corner */
 197   status = GRID_to_latlon(&gridnav, *nx, *ny, URLa, URLo);
 198   if (!status)
 199     {
 200       fprintf(stderr,
 201               "write_grib: error from GRID_to_latlon for first lat/lon\n");
 202       *ierr = 1;
 203       return(0);
 204     }
 205
 206   *ierr = 0;
 207   return(0);
 208
 209 }