source: TOOLS/MOSAIX/nemo.py @ 6174

Last change on this file since 6174 was 6105, checked in by omamce, 2 years ago

MOSAIX, by O.M : adapt to new file format of NEMO4.2

(no halo or peridodicity band in files)
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 21.2 KB
Line 
1# -*- coding: utf-8 -*-
2## ===========================================================================
3##
4##  Warning, to install, configure, run, use any of Olivier Marti's
5##  software or to read the associated documentation you'll need at least
6##  one (1) brain in a reasonably working order. Lack of this implement
7##  will void any warranties (either express or implied).
8##  O. Marti assumes no responsability for errors, omissions,
9##  data loss, or any other consequences caused directly or indirectly by
10##  the usage of his software by incorrectly or partially configured
11##  personal.
12##
13'''
14Utilities to plot NEMO ORCA fields
15Periodicity and other stuff
16
17olivier.marti@lsce.ipsl.fr
18'''
19
20import sys, numpy as np
21try    : import xarray as xr
22except : pass
23
24rpi = np.pi ; rad = rpi / 180.0
25
26nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2]
27
28## SVN information
29__Author__   = "$Author$"
30__Date__     = "$Date$"
31__Revision__ = "$Revision$"
32__Id__       = "$Id$"
33__HeadURL    = "$HeadURL$"
34
35def __guessNperio__ (jpj, jpi, nperio=None) :
36    '''
37    Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details)
38   
39    Inputs
40    jpj    : number of latitudes
41    jpi    : number of longitudes
42    nperio : periodicity parameter
43    '''
44    if nperio == None :
45        ## Values for NEMO version < 4.2
46        if jpi ==  182 : nperio = 4   #   ORCA2. We choose legacy orca2.
47        if jpi ==  362 : nperio = 6   #   eORCA1.
48        if jpi == 1442 : nperio = 6   #   ORCA025.
49
50        if jpj ==  149 : nperio = 4   #   ORCA2. We choose legacy orca2.
51        if jpj ==  332 : nperio = 6   #   eORCA1.
52
53        ## Values for NEMO version >= 4.2. No more halo points
54        if jpi ==  180 : nperio = 4.2 #   ORCA2. We choose legacy orca2.
55        if jpi ==  360 : nperio = 0   #   eORCA1.
56        if jpi == 1440 : nperio = 0   #   ORCA025.
57
58        if jpj ==  148 : nperio = 4.2 #   ORCA2. We choose legacy orca2.
59        if jpj ==  330 : nperio = 0   #   eORCA1.
60
61        #
62        if nperio == None :
63            sys.exit ('in nemo module : nperio not found, and cannot by guessed' )
64        else :
65            sys.ext  ('nperio set as {:} (deduced from jpi={:d}) : nemo.py is not ready for this value'.format (nperio, jpi) )
66           
67            print ('nperio set as {:} (deduced from jpi={:d})'.format (nperio, jpi))
68
69    return nperio
70
71def __guessPoint__ (ptab) :
72    '''
73    Tries to guess the grid point (periodicity parameter. See NEMO documentation for details)
74   
75    For array conforments with xgcm requirements
76
77    Inputs
78    ptab : xarray array
79
80    Credits : who is the original author ?
81    '''
82    gP = None
83    if isinstance (ptab, xr.core.dataarray.DataArray) :
84        if 'x_c' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'T'
85        if 'x_f' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'U'
86        if 'x_c' in ptab.dims and 'y_f' in ptab.dims                        : gP = 'V'
87        if 'x_f' in ptab.dims and 'y_f' in ptab.dims                        : gP = 'F'
88        if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_c' in ptab.dims : gP = 'T'
89        if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'W'
90        if 'x_f' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'U'
91        if 'x_c' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'V'
92        if 'x_f' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'F'
93             
94        if gP == None :
95            sys.exit ('in nemo module : cd_type not found, and cannot by guessed' )
96        else :
97            print ('Grid set as', gP, 'deduced from dims ', ptab.dims)
98            return gP
99    else :
100         sys.exit ('in nemo module : cd_type not found, input is not an xarray data' )
101         #return gP
102
103def fixed_lon (lon) :
104    '''
105    Returns corrected longitudes for nicer plots
106
107    fixed_lon (lon)
108    lon : longitudes of the grid. At least 2D.
109
110    Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7
111    '''
112    fixed_lon = lon.copy ()
113    for i, start in enumerate (np.argmax (np.abs (np.diff (lon, axis=-1)) > 180, axis=-1)) :
114        fixed_lon[..., i, start+1:] += 360
115
116    # Special case for eORCA025
117    if lon.shape[-1] == 1442 : lon [..., -2, :] =  lon [..., -3, :]
118    if lon.shape[-1] == 1440 : lon [..., -1, :] =  lon [..., -2, :]
119               
120    return fixed_lon
121
122def jeq (lat) :
123    '''
124    Returns j index of equator in the grid
125   
126    lat : latitudes of the grid. At least 2D.
127    '''
128    jeq = np.argmin (np.abs (np.float64(lat[:,0])))
129    return jeq
130
131def lon1D (lon, lat=None) :
132    '''
133    Returns 1D longitude for simple plots.
134   
135    lon : longitudes of the grid
136    lat (optionnal) : latitudes of the grid
137    '''
138    if np.max (lat) != None :
139        je     = jeq (lat)
140        lon1D = lon[..., je, :]
141    else :
142        jpj, jpi = lon.shape[-2:]
143        lon1D = lon[..., jpj//3, :]
144
145    start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1)
146    lon1D[..., start+1:] += 360
147       
148    return lon1D
149
150def latreg (lat, diff=0.1) :
151    '''
152    Returns maximum j index where gridlines are along latitude in the northern hemisphere
153   
154    lat : latitudes of the grid (2D)
155    diff [optional] : tolerance
156    '''
157    if diff == None :
158        dy = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2))))
159        diff = dy/100.
160   
161    je     = jeq (lat)
162    jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< diff)[-1][-1] + je
163    latreg = np.float64 (lat[...,jreg,:].mean(axis=-1))
164    JREG   = jreg
165
166    return jreg, latreg
167
168def lat1D (lat) :
169    '''
170    Returns 1D latitudes for zonal means and simple plots.
171
172    lat : latitudes of the grid (2D)
173    '''
174    jpj, jpi = lat.shape[-2:]
175    try :
176        if type (lat) == xr.core.dataarray.DataArray : math = xr
177        else : math = np
178    except : math = np
179
180    dy     = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2))))
181    je     = jeq (lat)
182    lat_eq = np.float64 (lat[...,je,:].mean(axis=-1))
183     
184    jreg, lat_reg = latreg (lat)
185    lat_ave = np.mean (lat, axis=-1)
186
187    if ( np.abs (lat_eq) < dy/100. ) : # T, U or W grid
188        dys    = (90.-lat_reg) / (jpj-jreg-1)*0.5
189        yrange = 90.-dys-lat_reg
190    else                             :  # V or F grid
191        yrange = 90.    -lat_reg
192       
193    lat1D = math.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) )   
194       
195    if math == xr : lat1D.attrs = lat.attrs
196
197    return lat1D
198
199def latlon1D (lat, lon) :
200    '''
201    Returns simple latitude and longitude (1D) for simple plots.
202
203    lat, lon : latitudes and longitudes of the grid (2D)
204    '''
205    return lat1D (lat),  lon1D (lon, lat)
206
207def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) :
208    try :
209        lon = lon.to_masked_array()
210        lat = lat.to_masked_array()
211    except : pass
212           
213    mask = np.logical_and ( np.logical_and(lat>y0, lat<y1), 
214        np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)),
215                                         np.logical_and(lon-360>x0, lon-360<x1)) )
216    try    : tab = xr.where (mask, ptab, np.nan)
217    except : tab = np.where (mask, ptab, np.nan)
218   
219    return tab
220       
221def extend (tab, Lon=False, jplus=25, jpi=None, nperio=4) :
222    '''
223    Returns extended field eastward to have better plots, and box average crossing the boundary
224    Works only for xarray and numpy data (?)
225
226    tab : field to extend.
227    Lon : (optional, default=False) : if True, add 360 in the extended parts of the field
228    jpi : normal longitude dimension of the field. exrtend does nothing it the actual
229        size of the field != jpi (avoid to extend several times)
230    jplus (optional, default=25) : number of points added on the east side of the field
231   
232    '''
233    if tab.shape[-1] == 1 :
234        extend = tab
235
236    else :
237        if jpi == None : jpi = tab.shape[-1]
238   
239        try : 
240            if type (tab) ==  xr.core.dataarray.DataArray : math = xr
241            else : math = np
242        except : math = np
243
244        if Lon : xplus = -360.0
245        else   : xplus =    0.0
246
247        if tab.shape[-1] > jpi :
248            extend = tab
249        else :
250            if nperio == 1 :
251                istart = 0 ; le=jpi+1 ; la=0
252            if nperio > 1 : # OPA case with to halo points for periodicity
253                istart = 1 ; le=jpi-2 ; la=1  # Perfect, except at the pole that should be masked by lbc_plot
254               
255            if math == xr : 
256                extend = xr.concat      ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), dim=tab.dims[-1] )
257            if math == np :
258                extend = np.concatenate ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), axis=-1 )
259           
260    return extend
261
262def orca2reg (ff, lat_name='nav_lat', lon_name='nav_lon', y_name='y', x_name='x') :
263    '''
264    Assign ORCA dataset on a regular grid.
265    For use in the tropical region.
266
267    Inputs :
268      ff : xarray dataset
269      lat_name, lon_name : name of latitude and longitude 2D field in ff
270      y_name, x_name     : namex of dimensions in ff
271
272    Returns : xarray dataset with rectangular grid. Incorrect above 20°N
273    '''
274    # Compute 1D longitude and latitude
275    (lat, lon) = latlon1D (ff[lat_name], ff[lon_name])
276    # Assign lon and lat as dimensions of the dataset
277    if y_name in ff.dims : 
278        lat = xr.DataArray (lat, coords=[lat,], dims=['lat',])     
279        ff  = ff.rename_dims ({y_name: "lat",}).assign_coords (lat=lat)
280    if x_name in ff.dims :
281        lon = xr.DataArray (lon, coords=[lon,], dims=['lon',])
282        ff  = ff.rename_dims ({x_name: "lon",}).assign_coords (lon=lon)
283    # Force dimensions to be in the right order
284    coord_order = ['lat', 'lon']
285    for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z',
286                 'time_counter', 'time', 'tbnds', 
287                 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] :
288        if dim in ff.dims : coord_order.insert (0, dim )
289       
290    ff = ff.transpose (*coord_order)
291    return ff
292
293def lbc_init (ptab, nperio=None, cd_type='T') :
294    """
295    Prepare for all lbc calls
296   
297    Set periodicity on input field
298    nperio    : Type of periodicity
299      0       : No periodicity
300      1, 4, 6 : Cyclic on i dimension (generaly longitudes) with 2 points halo
301      2       : Obsolete (was symmetric condition at southern boundary ?)
302      3, 4    : North fold T-point pivot (legacy ORCA2)
303      5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo)
304    cd_type   : Grid specification : T, U, V or F
305
306    See NEMO documentation for further details
307    """
308    jpj, jpi = ptab.shape[-2:]
309    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio)
310   
311    if nperio not in nperio_valid_range :
312        print ( 'nperio=', nperio, ' is not in the valid range', nperio_valid_range )
313        sys.exit()
314
315    if cd_type == None and isinstance (ptab, xr.core.dataarray.DataArray) :
316        cd_type = __guessPoint__ (ptab)
317
318    if cd_type not in ['T', 'U', 'V', 'F', 'W', 'F' ]:
319        print ( 'cd_type=', cd_type, ' is not in the list T, U, V, F, W, F' )
320        sys.exit ()
321
322    return jpj, jpi, nperio, cd_type
323       
324       
325def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) :
326    """
327    Set periodicity on input field
328    ptab      : Input array (works
329      rank 2 at least : ptab[...., lat, lon]
330    nperio    : Type of periodicity
331    cd_type   : Grid specification : T, U, V or F
332    psgn      : For change of sign for vector components (1 for scalars, -1 for vector components)
333   
334    See NEMO documentation for further details
335    """
336    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type)
337    psgn   = ptab.dtype.type (psgn)
338
339    #
340    #> East-West boundary conditions
341    # ------------------------------
342    if nperio in [1, 4, 6] :
343        # ... cyclic
344        ptab [..., :,  0] = ptab [..., :, -2]
345        ptab [..., :, -1] = ptab [..., :,  1]
346    #
347    #> North-South boundary conditions
348    # --------------------------------
349    if nperio in [3, 4] :  # North fold T-point pivot
350        if cd_type in [ 'T', 'W' ] : # T-, W-point
351            ptab [..., -1, 1:       ] = psgn * ptab [..., -3, -1:0:-1]
352            ptab [..., -1, 0        ] = psgn * ptab [..., -3, 2            ]
353            ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2:0:-1  ]
354               
355        if cd_type == 'U' :
356            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ]       
357            ptab [..., -1,  0       ] = psgn * ptab [..., -3,  1           ]
358            ptab [..., -1, -1       ] = psgn * ptab [..., -3, -2           ]
359               
360            if nemo_4U_bug :
361                ptab [..., -2, jpi//2+1:-1] = psgn * ptab [..., -2, jpi//2-2:0:-1]
362                ptab [..., -2, jpi//2-1   ] = psgn * ptab [..., -2, jpi//2     ]
363            else :
364                ptab [..., -2, jpi//2-1:-1] = psgn * ptab [..., -2, jpi//2:0:-1]
365               
366        if cd_type == 'V' : 
367            ptab [..., -2, 1:       ] = psgn * ptab [..., -3, jpi-1:0:-1   ]
368            ptab [..., -1, 1:       ] = psgn * ptab [..., -4, -1:0:-1      ]   
369            ptab [..., -1, 0        ] = psgn * ptab [..., -4, 2            ]
370               
371        if cd_type == 'F' :
372            ptab [..., -2, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ]
373            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -4, -1:0:-1      ]
374            ptab [..., -1,  0       ] = psgn * ptab [..., -4,  1           ]
375            ptab [..., -1, -1       ] = psgn * ptab [..., -4, -2           ]
376
377    if nperio in [4.2] :  # North fold T-point pivot
378        if cd_type in [ 'T', 'W' ] : # T-, W-point
379            ptab [..., -1, jpi//2:  ] = psgn * ptab [..., -1, jpi//2:0:-1  ]
380               
381        if cd_type == 'U' :
382            ptab [..., -1, jpi//2-1:-1] = psgn * ptab [..., -1, jpi//2:0:-1]
383               
384        if cd_type == 'V' : 
385            ptab [..., -1, 1:       ] = psgn * ptab [..., -2, jpi-1:0:-1   ]
386               
387        if cd_type == 'F' :
388            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -2, -1:0:-1      ]
389         
390       
391    if nperio in [5, 6] :            #  North fold F-point pivot 
392        if cd_type in ['T', 'W']  :
393            ptab [..., -1, 0:       ] = psgn * ptab [..., -2, -1::-1       ]
394               
395        if cd_type == 'U' :
396            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -2, -2::-1       ]       
397            ptab [..., -1, -1       ] = psgn * ptab [..., -2, 0            ] # Bug ?
398               
399        if cd_type == 'V' :
400            ptab [..., -1, 0:       ] = psgn * ptab [..., -3, -1::-1       ]
401            ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2-1::-1 ]
402               
403        if cd_type == 'F' :
404            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -2::-1       ]
405            ptab [..., -1, -1       ] = psgn * ptab [..., -3, 0            ]
406            ptab [..., -2, jpi//2:-1] = psgn * ptab [..., -2, jpi//2-2::-1 ]
407               
408    return ptab
409
410def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) :
411    #
412    """
413    Mask fields on duplicated points
414    ptab      : Input array
415      rank 2 at least : ptab[...., lat, lon]
416    nperio    : Type of periodicity
417    cd_type   : Grid specification : T, U, V or F
418   
419    See NEMO documentation for further details
420    """
421
422    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type)
423         
424    #
425    #> East-West boundary conditions
426    # ------------------------------
427    if nperio in [1, 4, 6] :
428        # ... cyclic
429        ptab [...,:,  0] = sval
430        ptab [...,:, -1] = sval
431
432    #
433    #> South (in which nperio cases ?)
434    # --------------------------------
435    if nperio in [1, 3, 4, 5, 6] :
436        ptab [...,0,:] = sval
437       
438    #
439    #> North-South boundary conditions
440    # --------------------------------
441    if nperio in [3, 4] :  # North fold T-point pivot
442        if cd_type in [ 'T', 'W' ] : # T-, W-point
443            ptab [..., -1, 1:       ] = sval
444            ptab [..., -1, 0        ] = sval
445            ptab [..., -2, jpi//2:  ] = sval
446               
447        if cd_type == 'U' :
448            ptab [..., -1,  :       ] = sval 
449            ptab [..., -2, jpi//2:  ] = sval
450               
451        if cd_type == 'V' :
452            ptab [..., -2, :       ] = sval
453            ptab [..., -1, :       ] = sval   
454               
455        if cd_type == 'F' :
456            ptab [..., -2, :     ] = sval
457            ptab [..., -1, :     ] = sval
458
459    if nperio in [4.2] :  # North fold T-point pivot
460        if cd_type in [ 'T', 'W' ] : # T-, W-point
461            ptab [..., -1, jpi//2:  ] = sval
462               
463        if cd_type == 'U' :
464            ptab [..., -1, jpi//2-1:-1] = sval
465               
466        if cd_type == 'V' : 
467            ptab [..., -1, 1:       ] = sval
468               
469        if cd_type == 'F' :
470            ptab [..., -1, 0:-1     ] = sval
471   
472    if nperio in [5, 6] :            #  North fold F-point pivot
473        if cd_type in ['T', 'W']  :
474            ptab [..., -1, 0:       ] = sval
475               
476        if cd_type == 'U' :
477            ptab [..., -1, 0:-1     ] = sval       
478            ptab [..., -1, -1       ] = sval
479             
480        if cd_type == 'V' :
481            ptab [..., -1, 0:       ] = sval
482            ptab [..., -2, jpi//2:  ] = sval
483                             
484        if cd_type == 'F' :
485            ptab [..., -1, 0:-1     ] = sval
486            ptab [..., -1, -1       ] = sval
487            ptab [..., -2, jpi//2:-1] = sval
488
489    return ptab
490
491
492def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) :
493    """
494    Set periodicity on input field, adapted for plotting for any cartopy projection
495    ptab      : Input array (works
496      rank 2 at least : ptab[...., lat, lon]
497    nperio    : Type of periodicity
498    cd_type   : Grid specification : T, U, V or F
499    psgn      : For change of sign for vector components (1 for scalars, -1 for vector components)
500   
501    See NEMO documentation for further details
502    """
503
504    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type)
505    psgn   = ptab.dtype.type (psgn)
506   
507    #
508    #> East-West boundary conditions
509    # ------------------------------
510    if nperio in [1, 4, 6] :
511        # ... cyclic
512        ptab [..., :,  0] = ptab [..., :, -2]
513        ptab [..., :, -1] = ptab [..., :,  1]
514       
515    #
516    #> North-South boundary conditions
517    # --------------------------------
518    if nperio in [3, 4] :  # North fold T-point pivot
519        if cd_type in [ 'T', 'W' ] : # T-, W-point
520            ptab [..., -1, :] = sval
521            ptab [..., -2, :jpi//2] = sval
522            #ptab [..., -2, :] = sval
523                     
524        if cd_type == 'U' :
525            ptab [..., -1, : ] = sval
526
527        if cd_type == 'V' : 
528            ptab [..., -2, : ] = sval
529            ptab [..., -1, : ] = sval
530           
531        if cd_type == 'F' :
532            ptab [..., -2, : ] = sval
533            ptab [..., -1, : ] = sval
534
535    if nperio in [4.2] :  # North fold T-point pivot
536        if cd_type in [ 'T', 'W' ] : # T-, W-point
537            ptab [..., -1, jpi//2:  ] = sval
538               
539        if cd_type == 'U' :
540            ptab [..., -1, jpi//2-1:-1] = sval
541               
542        if cd_type == 'V' : 
543            ptab [..., -1, 1:       ] = sval
544               
545        if cd_type == 'F' :
546            ptab [..., -1, 0:-1     ] = sval
547     
548    if nperio in [5, 6] :            #  North fold F-point pivot 
549        if cd_type in ['T', 'W']  :
550            ptab [..., -1, : ] = sval
551               
552        if cd_type == 'U' :
553            ptab [..., -1, : ] = sval     
554             
555        if cd_type == 'V' :
556            ptab [..., -1, : ] = sval
557            ptab [..., -2, jpi//2:  ] = sval
558                             
559        if cd_type == 'F' :
560            ptab [..., -1, : ] = sval
561            ptab [..., -2, jpi//2:-1] = sval
562
563    return ptab
564
565def geo2oce (pxx, pyy, pzz, glam, gphi) : 
566    '''
567    Change vector from geocentric to east/north
568
569    Input
570        pxx, pyy, pzz : components on the geoce,tric system
571        glam, gphi : longitue and latitude of the points
572
573    '''
574    gsinlon = np.sin (rad * glam)
575    gcoslon = np.cos (rad * glam)
576    gsinlat = np.sin (rad * gphi)
577    gcoslat = np.cos (rad * gphi)
578         
579    pte = - pxx * gsinlon            + pyy * gcoslon
580    ptn = - pxx * gcoslon * gsinlat  - pyy * gsinlon * gsinlat + pzz * gcoslat
581
582    return pte, ptn
583
584def oce2geo (pte, ptn, glam, gphi) :
585    ##----------------------------------------------------------------------
586    ##                    ***  ROUTINE oce2geo  ***
587    ##     
588    ## ** Purpose :
589    ##
590    ## ** Method  :   Change vector from east/north to geocentric
591    ##
592    ## History :
593    ##        #         (A. Caubel)  oce2geo - Original code
594    ##----------------------------------------------------------------------
595    gsinlon = np.sin (rad * glam)
596    gcoslon = np.cos (rad * glam)
597    gsinlat = np.sin (rad * gphi)
598    gcoslat = np.cos (rad * gphi)
599
600    pxx = - pte * gsinlon - ptn * gcoslon * gsinlat
601    pyy =   pte * gcoslon - ptn * gsinlon * gsinlat
602    pzz =   ptn * gcoslat
603   
604    return pxx, pyy, pzz
605
606## ===========================================================================
607##
608##                               That's all folk's !!!
609##
610## ===========================================================================
Note: See TracBrowser for help on using the repository browser.