source: lmdz_wrf/trunk/WRFV3/external/io_grib_share/gridnav.c @ 616

Last change on this file since 616 was 1, checked in by lfita, 11 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 16.2 KB
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4#include "gridnav.h"
5
6#define MISSING -999
7#define PI 3.1415927
8#define RAD_TO_DEG 57.29577951
9/* #define EARTH_RAD 6371.200 */
10#define EARTH_RAD 6370.000
11
12int fill_proj_parms(GridNav *gridnav);
13
14/*****************************************************************************
15 * Fills up the gridnav structure using arguments input to the function
16 *
17 * input:
18 *   central_lat - the central latitude of the projection, in degrees
19 *   central_lon - the central longitude of the projection, in degrees
20 *   projection  - the projection number (defined by #define's in gridnav.h)
21 *   truelat1    - the first "true" latitude.  Only has an effect with
22 *                 polar stereographic and lambert projections
23 *   truelat2    - the second true latitude.  Only has an effect with lambert
24 *                 projection
25 *   num_columns - number of columns in grid
26 *   num_rows    - number of rows in grid
27 *   dx,dy       - east-west (dx) and north-south (dy) distance between grid
28 *                 points, in degrees for latlon (i.e., cylindrical
29 *                 equidistant) projection, in km for all other projections
30 *                 (including mercator).
31 *   lat_origin  - latitude of grid point defined in origin_column, origin_row
32 *   lon_origin  - longitude of grid point defined in origin_column, origin_row
33 *   origin_column - column for lat_origin, long_origin pair
34 *   origin_row  - row for lat_origin, long_origin pair
35 * 
36 * output:
37 *   gridnav - filled gridnav structure
38 *
39 *****************************************************************************/
40
41int GRID_init(float central_lat, float central_lon, int projection, 
42              float truelat1, float truelat2, int num_columns, 
43              int num_rows, float dx, float dy, float lat_origin, 
44              float lon_origin, float origin_column, float origin_row, 
45              GridNav *gridnav)
46{
47
48  int status = 1;
49
50  gridnav->proj.central_lon = central_lon;
51  gridnav->proj.central_lat = central_lat;
52  gridnav->proj.map_proj = projection;
53  gridnav->proj.truelat1 = truelat1;
54  gridnav->proj.truelat2 = truelat2;
55
56  gridnav->grid.num_columns = num_columns;
57  gridnav->grid.num_rows = num_rows;
58  gridnav->grid.dx = dx;
59  gridnav->grid.dy = dy;
60  gridnav->grid.lat_origin = lat_origin;
61  gridnav->grid.lon_origin = lon_origin;
62  gridnav->grid.origin_column = origin_column;
63  gridnav->grid.origin_row = origin_row;
64
65  fill_proj_parms(gridnav);
66
67  return status;
68}
69
70/*****************************************************************************
71 * Outputs the latitude and longitude of a grid point.
72 *
73 * input:
74 *   *gridnav - a pointer to a filled gridnav structure.  The gridnav
75 *              structure should have been filled using one of the init
76 *              functions.
77 *   column   - the column in the grid to retrieve.
78 *   row      - the row in the grid to retrieve.
79 *
80 * output:
81 *   *lat - pointer to output latitude
82 *   *lon - pointer to output longitude
83 *
84 *****************************************************************************/
85
86int GRID_to_latlon(GridNav *gridnav, float column, float row, float *lat, 
87                   float *lon)
88{
89    /*
90     * Calculate the latitude and longitude for an input grid point location
91     */
92    double X, Y, R;
93    int status = 1;
94    double eps = 0.00001;
95
96    switch (gridnav->proj.map_proj)
97        {
98        case GRID_MERCATOR:
99            *lat = 2 * RAD_TO_DEG *
100                atan(exp((gridnav->proj_transform.parm2 + gridnav->grid.dx *
101                          (row - gridnav->grid.origin_row)) /
102                         gridnav->proj_transform.parm1 )) - 90;
103            *lon = gridnav->grid.lon_origin + RAD_TO_DEG * gridnav->grid.dx *
104                (column - gridnav->grid.origin_column) /
105                gridnav->proj_transform.parm1;
106            break;
107
108        case GRID_LAMCON:
109            X = (column - gridnav->grid.origin_column) * gridnav->grid.dx +
110                gridnav->proj_transform.parm6;
111            Y = (row - gridnav->grid.origin_row) * gridnav->grid.dy +
112                gridnav->proj_transform.parm3 + gridnav->proj_transform.parm7;
113            R = sqrt(X*X + Y*Y);
114            *lat = gridnav->proj_transform.parm5 * 90 - 
115                2 * RAD_TO_DEG * atan(gridnav->proj_transform.parm4 * 
116                                      pow(R, 1 / 
117                                          gridnav->proj_transform.parm2));
118            *lon = gridnav->proj.central_lon + 
119                RAD_TO_DEG / gridnav->proj_transform.parm2 *
120                atan(X / (gridnav->proj_transform.parm5 * -Y));
121            while (*lon > 180) *lon -= 360;
122            while (*lon <= -180) *lon += 360;
123            break;
124
125        case GRID_POLSTR:
126            X = (column - gridnav->grid.origin_column) * 
127                gridnav->proj_transform.parm3 + 
128                gridnav->proj_transform.parm1;
129            Y = (row - gridnav->grid.origin_row) * 
130                gridnav->proj_transform.parm3 +
131                gridnav->proj_transform.parm2 + eps;
132            *lon = gridnav->proj_transform.parm5 * -1 * 
133                atan(X / Y) * RAD_TO_DEG + gridnav->proj.central_lon;
134            *lat = 90 - 2 * RAD_TO_DEG *
135                atan(X / (2 * EARTH_RAD * 
136                          sin((*lon - gridnav->proj.central_lon ) / 
137                              RAD_TO_DEG) + eps) );
138            while (*lon > 180) *lon -= 360;
139            while (*lon <= -180) *lon += 360;
140            break;
141     
142        case GRID_LATLON:
143            *lat = (row - gridnav->grid.origin_row)*gridnav->grid.dy +
144                gridnav->grid.lat_origin;
145            *lon = (column - gridnav->grid.origin_column)*gridnav->grid.dx +
146                gridnav->grid.lon_origin;
147            break;
148
149        default:
150            /*
151             *  Unsupported map projection: set lat-lon to no-data values.
152             */ 
153            fprintf(stderr,"GRID_to_latlon: Unsupport map projection type %d\n", 
154                    gridnav->proj.map_proj); 
155            *lon = -9999; 
156            *lat = -9999; 
157            status = 0; 
158            break; 
159
160        } /* end of switch on map projection type */
161
162
163    return status;
164}
165
166
167/*****************************************************************************
168 * Outputs grid point indices, given lat/lon location.
169 *
170 * input:
171 *   *gridnav - a pointer to a filled gridnav structure.  The gridnav
172 *              structure should have been filled using one of the init
173 *              functions.
174 *   lat      - latitude for location
175 *   lon      - longitude for location
176 * output:
177 *   *column  - pointer to value for column
178 *   *row     - pointer to value for row
179 *
180 *****************************************************************************/
181
182int GRID_from_latlon(GridNav *gridnav, float lat, float lon, float *column, 
183                     float *row)
184{
185    double X, Y, Rs;
186    double lat_rad;
187    int status = 1;
188
189    switch (gridnav->proj.map_proj) 
190        {
191        case GRID_MERCATOR:
192            X = gridnav->proj_transform.parm1 *
193                ((lon - gridnav->grid.lon_origin) / RAD_TO_DEG);
194            *column = gridnav->grid.origin_column + X / gridnav->grid.dx;
195            lat_rad = lat / RAD_TO_DEG;
196            Y = gridnav->proj_transform.parm1 * log((1 + sin(lat_rad)) / 
197                                                    cos(lat_rad));
198            *row = gridnav->grid.origin_row + 
199                ((Y-gridnav->proj_transform.parm2) / gridnav->grid.dy);;
200            break;
201
202        case GRID_LAMCON:
203            Rs = (EARTH_RAD / gridnav->proj_transform.parm2) *
204                sin(gridnav->proj_transform.parm1) *
205                pow(tan((gridnav->proj_transform.parm5 * 90 - lat) / 
206                        (2 * RAD_TO_DEG))/
207                    tan(gridnav->proj_transform.parm1 / 2),
208                    gridnav->proj_transform.parm2);
209            *row = gridnav->grid.origin_row - 
210                (1 / gridnav->grid.dy) * 
211                (gridnav->proj_transform.parm3 + 
212                 Rs * cos(gridnav->proj_transform.parm2*
213                          (lon - gridnav->proj.central_lon) / RAD_TO_DEG)) -
214                gridnav->proj_transform.parm7 / gridnav->grid.dy;
215            *column = gridnav->grid.origin_column + 
216                ( gridnav->proj_transform.parm5*
217                ((Rs / gridnav->grid.dx)*
218                 sin(gridnav->proj_transform.parm2 * 
219                     (lon - gridnav->proj.central_lon) / RAD_TO_DEG) - 
220                  gridnav->proj_transform.parm6 / gridnav->grid.dx));
221            break;
222
223        case GRID_POLSTR:           
224            Rs = 2 * EARTH_RAD * tan ( (45 - fabs(lat)/2) / RAD_TO_DEG );
225            Y = gridnav->proj_transform.parm5 * -1 * Rs * 
226                cos ((lon - gridnav->proj.central_lon) / RAD_TO_DEG);
227            X = Rs * sin ((lon - gridnav->proj.central_lon) / RAD_TO_DEG);
228
229            *row = (Y - gridnav->proj_transform.parm2) / 
230                gridnav->proj_transform.parm3
231                + gridnav->grid.origin_row;
232            *column = (X - gridnav->proj_transform.parm1) / 
233                gridnav->proj_transform.parm3
234                + gridnav->grid.origin_column;
235            break;
236   
237        case GRID_LATLON:
238            *row = ((lat - gridnav->grid.lat_origin) / gridnav->grid.dy) +
239                gridnav->grid.origin_row;
240            /* If lon is negative, make it positive */
241            while (lon < 0) lon += 360.;
242            *column = ((lon - gridnav->grid.lon_origin) / gridnav->grid.dx) +
243                gridnav->grid.origin_column;
244            break;
245
246        default:
247            fprintf(stderr,"GRID_from_latlon: Unsupported map projection type %d\n",
248                    gridnav->proj.map_proj);
249            *column = -9999;
250            *row = -9999;
251            status = 0;
252            break;
253
254        } /* End of switch on map projection type */
255
256    return status;
257}
258
259int fill_proj_parms(GridNav *gridnav)
260{
261    double orig_lat_rad;
262    double R_orig;
263    double r_not;
264    double r_truelat1;
265    int hemifactor;
266 
267    switch (gridnav->proj.map_proj) 
268        {
269        case GRID_MERCATOR:
270            gridnav->proj_transform.parm1 = 
271                EARTH_RAD * cos(gridnav->proj.truelat1 / RAD_TO_DEG);
272            orig_lat_rad = (gridnav->grid.lat_origin) / RAD_TO_DEG;
273            gridnav->proj_transform.parm2 = 
274                gridnav->proj_transform.parm1 *
275                log((1. + sin(orig_lat_rad)) / cos(orig_lat_rad));
276            gridnav->proj_transform.parm3 = MISSING;
277            gridnav->proj_transform.parm4 = MISSING;
278            gridnav->proj_transform.parm5 = MISSING;
279            break;
280        case GRID_LAMCON:
281            if (gridnav->proj.truelat1 >= 0) 
282                {
283                    hemifactor = 1;
284                } 
285            else 
286                {
287                    hemifactor = -1;
288                }
289            /* This is Psi1 in MM5 speak */
290            gridnav->proj_transform.parm1 = 
291                hemifactor*(PI/2 - fabs(gridnav->proj.truelat1) / RAD_TO_DEG);
292            /* This is Kappa in MM5 speak */
293            if (fabs(gridnav->proj.truelat1 - gridnav->proj.truelat2) < .01) 
294                {
295                    gridnav->proj_transform.parm2 = 
296                        fabs(sin(gridnav->proj.truelat1 / RAD_TO_DEG));
297                }
298            else 
299                {
300                    gridnav->proj_transform.parm2 = 
301                        (log10(cos(gridnav->proj.truelat1 / RAD_TO_DEG)) - 
302                         log10(cos(gridnav->proj.truelat2 / RAD_TO_DEG))) /
303                        (log10(tan((45 - fabs(gridnav->proj.truelat1) / 2) / RAD_TO_DEG)) - 
304                         log10(tan((45 - fabs(gridnav->proj.truelat2) / 2) / RAD_TO_DEG)));
305                }
306            /* This is Yc in MM5 speak */
307            gridnav->proj_transform.parm3 = 
308                -EARTH_RAD/gridnav->proj_transform.parm2 * 
309                sin(gridnav->proj_transform.parm1) *
310                pow(
311                    (tan((hemifactor * 90 - gridnav->grid.lat_origin) / 
312                         (RAD_TO_DEG * 2)) /
313                     tan(gridnav->proj_transform.parm1 / 2)),
314                    gridnav->proj_transform.parm2);
315            gridnav->proj_transform.parm4 = 
316                tan(gridnav->proj_transform.parm1 / 2)*
317                pow(hemifactor * (gridnav->proj_transform.parm2)/
318                    (EARTH_RAD * sin(gridnav->proj_transform.parm1)),
319                    1 / gridnav->proj_transform.parm2);
320            gridnav->proj_transform.parm5 = hemifactor;
321            R_orig = (EARTH_RAD / gridnav->proj_transform.parm2) *
322                sin(gridnav->proj_transform.parm1) *
323                pow(tan((gridnav->proj_transform.parm5 * 90 - 
324                         gridnav->grid.lat_origin) / 
325                        (2 * RAD_TO_DEG))/
326                    tan(gridnav->proj_transform.parm1 / 2),
327                    gridnav->proj_transform.parm2);
328            /* X origin */
329            gridnav->proj_transform.parm6 = 
330                R_orig * sin(gridnav->proj_transform.parm2 * 
331                             (gridnav->grid.lon_origin - 
332                              gridnav->proj.central_lon) / RAD_TO_DEG);
333
334            /* Y origin */
335            gridnav->proj_transform.parm7 = -(gridnav->proj_transform.parm3 + 
336                R_orig * cos(gridnav->proj_transform.parm2 * 
337                             (gridnav->grid.lon_origin - 
338                              gridnav->proj.central_lon) / RAD_TO_DEG));
339            break;
340        case GRID_POLSTR:
341            if (gridnav->proj.central_lat > 0) 
342                {
343                    hemifactor = 1;
344                } 
345            else 
346                {
347                    hemifactor = -1;
348                }
349
350            /* Calculate X for origin */
351            r_not = 2 * EARTH_RAD * 
352                tan((45 - fabs(gridnav->grid.lat_origin) / 2) / RAD_TO_DEG);
353            gridnav->proj_transform.parm1 = 
354                r_not * 
355                sin ( (gridnav->grid.lon_origin - gridnav->proj.central_lon) /
356                      RAD_TO_DEG);
357            /* Calculate Y for origin */
358            gridnav->proj_transform.parm2 = 
359                hemifactor * -1 * r_not * 
360                cos ( (gridnav->grid.lon_origin - gridnav->proj.central_lon) / 
361                      RAD_TO_DEG);
362            /* Calculate grid spacing at pole */
363            r_truelat1 = 2 * EARTH_RAD * 
364                tan((45 - fabs(gridnav->proj.truelat1) / 2) / RAD_TO_DEG);
365            gridnav->proj_transform.parm3 = 
366                gridnav->grid.dx * r_truelat1 / 
367                ( EARTH_RAD * cos (gridnav->proj.truelat1 / RAD_TO_DEG));
368            gridnav->proj_transform.parm4 = MISSING;
369            gridnav->proj_transform.parm5 = hemifactor;
370            break;
371        case GRID_LATLON:
372          gridnav->proj_transform.parm1 = MISSING;
373          gridnav->proj_transform.parm2 = MISSING;
374          gridnav->proj_transform.parm3 = MISSING;
375          gridnav->proj_transform.parm4 = MISSING;
376          gridnav->proj_transform.parm5 = MISSING;
377          break;
378
379        default:
380            fprintf(stderr,"GRID_init_mm5data: Invalid Projection\n");
381            return 0;
382        }
383    return 1;
384}
385
386/*****************************************************************************
387 * Rotates u and v components of wind, from being relative to earth, to
388 *  relative to grid.
389 *
390 * input:
391 *   *gridnav - a pointer to a filled gridnav structure.  The gridnav
392 *              structure should have been filled using one of the init
393 *              functions.
394 *   lon      - longitude for location
395 *   u_earth  - the u component of the wind in earth (north-relative)
396 *              coordinates.
397 *   v_earth  - the v component of the wind in earth (north-relative)
398 *              coordinates.
399 * output:
400 *   *u_grid  - pointer to value for u in grid coordinates
401 *   *v_grid  - pointer to value for v in grid coordinates
402 *
403 *****************************************************************************/
404
405int GRID_rotate_from_earth_coords(GridNav *gridnav, float lon, float u_earth, 
406                                  float v_earth, float *u_grid, float *v_grid)
407{
408    float speed, dir;
409    float dir_grid;
410
411    /* Calculate Speed and Direction from u,v */
412    switch (gridnav->proj.map_proj) 
413        {
414        case GRID_MERCATOR:
415            *u_grid = u_earth;
416            *v_grid = v_earth;
417            break;
418        case GRID_POLSTR: case GRID_LAMCON:
419            speed = sqrt(u_earth * u_earth + v_earth * v_earth);
420            dir = RAD_TO_DEG * atan2(-u_earth,-v_earth);
421            while (dir >= 360) dir -= 360;
422            while (dir < 0) dir += 360;
423           
424            dir_grid = dir - (lon - gridnav->proj.central_lon) * 
425                gridnav->proj_transform.parm2;
426            while (dir_grid >= 360) dir_grid -= 360;
427            while (dir_grid < 0) dir_grid += 360;
428           
429            *u_grid = -1. * speed * sin(dir_grid / RAD_TO_DEG);
430            *v_grid = -1. * speed * cos(dir_grid / RAD_TO_DEG);
431            break;
432        default:
433            fprintf(stderr,
434                    "GRID_rotate_from_earth_coords: Invalid Projection\n");
435            return 0;
436        }  /* End of switch projection */
437
438    return 1;
439}
440
441/*****************************************************************************
442 * Rotates u and v components of wind, from being relative to earth, to
443 *  relative to grid.
444 *
445 * input:
446 *   *gridnav - a pointer to a filled gridnav structure.  The gridnav
447 *              structure should have been filled using one of the init
448 *              functions.
449 *   lon      - longitude for location
450 *   u_grid   - the u component of the wind in grid coordinates
451 *   v_grid   - the v component of the wind in grid coordinates
452 *
453 * output:
454 *   *u_earth - pointer to value for u in earth-relative coordinates
455 *   *v_grid  - pointer to value for v in earth-relative coordinates
456 *
457 *****************************************************************************/
458
459int GRID_rotate_to_earth_coords(GridNav *gridnav, float lon, float u_grid, 
460                                float v_grid, float *u_earth, float *v_earth)
461{
462    float speed, dir_grid;
463    float dir_earth;
464   
465    /* Calculate Speed and Direction from u,v */
466    switch (gridnav->proj.map_proj) 
467        {
468        case GRID_MERCATOR:
469            *u_earth = u_grid;
470            *v_earth = v_grid;
471            break;
472        case GRID_POLSTR: case GRID_LAMCON:
473            speed = sqrt(u_grid * u_grid + v_grid * v_grid);
474            dir_grid = RAD_TO_DEG * atan2(-u_grid, -v_grid);
475            while (dir_grid >= 360) dir_grid -= 360;
476            while (dir_grid < 0) dir_grid += 360;
477           
478            dir_earth = dir_grid + (lon - gridnav->proj.central_lon) *
479                gridnav->proj_transform.parm2;
480           
481            while (dir_earth >= 360) dir_earth -= 360;
482            while (dir_earth < 0) dir_earth += 360;
483           
484            *u_earth = -1. * speed * sin(dir_earth / RAD_TO_DEG);
485            *v_earth = -1. * speed * cos(dir_earth / RAD_TO_DEG);
486            break;
487        default:
488            fprintf(stderr,
489                    "GRID_rotate_to_earth_coords: Invalid Projection\n");
490            return 0;
491        }  /* End of switch projection */
492    return 1;
493}
Note: See TracBrowser for help on using the repository browser.