Changeset 6064


Ignore:
Timestamp:
02/24/22 09:58:49 (2 years ago)
Author:
omamce
Message:

O.M. : changes in MOSAIX

  • Correct ocean mask to remove periodic duplicated points, to have full conservation of run-off
  • Suppress creation of corc mask, which is not used
  • Put SVN keywords in make_mosaic
  • Update nemo.py with somùe additionnal utilities
  • Adapt RunOffWeights?.py to new nemo module
  • Adapt CalvingWeights?.py to new ORCA configurations
Location:
TOOLS/MOSAIX
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/MOSAIX/CalvingWeights.py

    r5915 r6064  
    122122    if dst_Name == 'eORCA025.1' : nperio_dst = 6 
    123123print ('nperio_dst: ' + str(nperio_dst) ) 
    124 dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' ) 
     124dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T', sval=0 ) 
    125125print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) 
    126126 
    … …  
    131131# Preparation by closing some straits 
    132132# ----------------------------------- 
    133 if dst_Name in [ 'eORCA025', 'eORCA025.1' ] : 
     133if dst_Name in [ 'eORCA025', 'eORCA025.1', 'eORCA025.4' ] : 
    134134    print ('Filling some seas for ' + dst_Name ) 
    135135    # Set Gibraltar strait to 0 to fill Mediterranean sea 
    … …  
    144144    dst_mskutil[997:1012,907] = 0 
    145145     
    146 if dst_Name == 'eORCA1.2' : 
    147     print ('Filling some seas for eORCA1.2') 
     146if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4' ] : 
     147    print ('Filling some seas for ' + dst_Name) 
    148148    # Set Gibraltar strait to 0 to fill Mediterrean sea 
    149149    dst_mskutil[240, 283] = 0 
    … …  
    157157    dst_mskutil[284,222:225] = 0 
    158158     
    159 if dst_Name == 'ORCA2.3' : 
    160     print ('Filling some seas for ORCA2.3') 
     159if dst_Name in [ 'ORCA2.3', 'ORCA2.4' ] : 
     160    print ('Filling some seas for ' + dst_Name) 
    161161    # Set Gibraltar strait to 0 to fill Mediterrean sea 
    162162    dst_mskutil[101,139] = 0 
    … …  
    175175# Fill closed seas with image processing library 
    176176# ---------------------------------------------- 
    177 dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), nperio=nperio_dst, cd_type='T' ) 
     177dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), nperio=nperio_dst, cd_type='T', sval=0) 
    178178 
    179179# Surfaces 

  • TOOLS/MOSAIX/CreateOasisGrids.bash

    r5159 r6064  
    5151if [[ $(hostname) = irene* ]]    ; then arch=irene ; center=tgcc ; fi 
    5252if [[ $(hostname) = lsce3005* ]] ; then arch=spip  ; center=spip ; fi 
     53if [[ $(hostname) = lsce5138* ]] ; then arch=spip  ; center=spip ; fi 
    5354 
    5455PROGRAM=$(basename ${0}) 
    … …  
    152153    ncks -C --history --append --variable area_grid_${OCEGRID} ${OCE}_coordinates_mask.nc areas_${CplModel}.nc 
    153154    # Inverts mask values and switch to integer 
    154     ncks --history -C --variable mask_${OCEGRID} ${OCE}_coordinates_mask.nc mask_${OCEGRID}_tmp.nc 
    155     ncatted --history \ 
    156             --attribute coordinates,mask_${OCEGRID},d,,         \ 
    157             --attribute online_operation,mask_${OCEGRID},d,,    \ 
    158             --attribute cell_measures,mask_${OCEGRID},d,,       \ 
     155    ncks --history -C --variable maskutil_${OCEGRID} ${OCE}_coordinates_mask.nc mask_${OCEGRID}_tmp.nc 
     156    ncatted --history \ 
     157            --attribute coordinates,maskutil_${OCEGRID},d,,         \ 
     158            --attribute online_operation,maskutil_${OCEGRID},d,,    \ 
     159            --attribute cell_measures,maskutil_${OCEGRID},d,,       \ 
    159160    mask_${OCEGRID}_tmp.nc 
    160161     
    161     ncap2 --history --append --script "mask_${OCEGRID}=int(1-mask_${OCEGRID});"    mask_${OCEGRID}_tmp.nc masks_${CplModel}.nc 
     162    ncap2 --history --append --script "maskutil_${OCEGRID}=int(1-maskutil_${OCEGRID});"    mask_${OCEGRID}_tmp.nc masks_${CplModel}.nc 
    162163    rm mask_${OCEGRID}_tmp.nc 
    163164    ncatted --history \ 
    164             --attribute long_name,mask_${OCEGRID},o,c,"Land-sea mask" \ 
    165             --attribute units,mask_${OCEGRID},o,c,"Land:1, Ocean:0"    masks_${CplModel}.nc 
     165            --attribute long_name,maskutil_${OCEGRID},o,c,"Land-sea mask" \ 
     166            --attribute units,maskutil_${OCEGRID},o,c,"Land:1, Ocean:0"    masks_${CplModel}.nc 
    166167    # Change order of dimensions 
    167168    mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc 
    … …  
    327328 
    328329 
    329 echo ${Titre}"Generates grid c${atm}, 'o' meaning 'one'"${Norm} 
    330 # same as t${atm} grid, with surfaces set to grid area 
    331 # and mask to 0 (ocean everywhere, to compute integral over the whole grid)) 
    332 # ---------------------------------------------------------------------------- 
    333 mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc 
    334 ncap2 --history --script "c${atm}_lon=alon ; c${atm}_lat=alat ; bounds_c${atm}_lon=bounds_lon ; bounds_c${atm}_lat=bounds_lat ; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc 
    335  
    336 mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc 
    337 ncap2 --history --script "c${atm}_aire=aire"   areas_${CplModel}_tmp.nc areas_${CplModel}.nc 
    338  
    339 mv masks_${CplModel}.nc masks_${CplModel}_tmp.nc 
    340 ncap2 --history --script "c${atm}_mask=int(Oce2AtmMask)*0+0 ; " masks_${CplModel}_tmp.nc masks_${CplModel}.nc 
    341  
    342 rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc masks_${CplModel}_tmp.nc 
     330# echo ${Titre}"Generates grid c${atm}, 'o' meaning 'one'"${Norm} 
     331# # same as t${atm} grid, with surfaces set to grid area 
     332# # and mask to 0 (ocean everywhere, to compute integral over the whole grid)) 
     333# # ---------------------------------------------------------------------------- 
     334# mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc 
     335# ncap2 --history --script "c${atm}_lon=alon ; c${atm}_lat=alat ; bounds_c${atm}_lon=bounds_lon ; bounds_c${atm}_lat=bounds_lat ; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc 
     336 
     337# mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc 
     338# ncap2 --history --script "c${atm}_aire=aire"   areas_${CplModel}_tmp.nc areas_${CplModel}.nc 
     339 
     340# mv masks_${CplModel}.nc masks_${CplModel}_tmp.nc 
     341# ncap2 --history --script "c${atm}_mask=int(Oce2AtmMask)*0+0 ; " masks_${CplModel}_tmp.nc masks_${CplModel}.nc 
     342 
     343# rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc masks_${CplModel}_tmp.nc 
    343344 
    344345echo ${Titre}"Generates grid o${oce}, 'o' meaning 'one'"${Norm} 
    345346# same as t${oce} grid, with surfaces set to 1 
    346347# and mask to 0 (ocean everywhere, to compute integral over the whole grid)) 
    347 # This grid is used when field are quantities instead of fluxes (i.e river flow) 
     348# This grid is used when field are quantities instead of fluxes 
    348349# -------------------------------------------------------------------------------------------------------- 
    349350mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc 
    … …  
    354355 
    355356ncks -C --history --overwrite -v maskutil_T ${OCE}_coordinates_mask.nc maskutil_T.nc 
    356 ncap2 --history --append --script "mask_O=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc 
     357ncap2 --history --append --script "maskutil_O=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc 
    357358 
    358359rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc  
    359360 
    360361 
    361 echo ${Titre}"Creates OCEAN C grid : redundant points removed to compute proper integrals"${Norm} 
    362 # -------------------------------------------------------------------------------------------------------- 
    363  
    364 mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc 
    365 ncap2 --history --script "nav_lon_grid_C=nav_lon_grid_T; nav_lat_grid_C=nav_lat_grid_T; bounds_lon_grid_C=bounds_lon_grid_T; bounds_lat_grid_C=bounds_o${oce}_lat=bounds_lat_grid_T; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc 
    366  
    367 mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc 
    368 ncap2 --history --script "area_grid_C=area_grid_T ; "   areas_${CplModel}_tmp.nc areas_${CplModel}.nc 
    369  
    370 ncap2 --history --append --script "mask_C=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc 
    371  
    372 rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc maskutil_T.nc 
     362# echo ${Titre}"Creates OCEAN C grid : redundant points removed to compute proper integrals"${Norm} 
     363# # -------------------------------------------------------------------------------------------------------- 
     364 
     365# mv grids_${CplModel}.nc grids_${CplModel}_tmp.nc 
     366# ncap2 --history --script "nav_lon_grid_C=nav_lon_grid_T; nav_lat_grid_C=nav_lat_grid_T; bounds_lon_grid_C=bounds_lon_grid_T; bounds_lat_grid_C=bounds_o${oce}_lat=bounds_lat_grid_T; " grids_${CplModel}_tmp.nc grids_${CplModel}.nc 
     367 
     368# mv areas_${CplModel}.nc areas_${CplModel}_tmp.nc 
     369# ncap2 --history --script "area_grid_C=area_grid_T ; "   areas_${CplModel}_tmp.nc areas_${CplModel}.nc 
     370 
     371# ncap2 --history --append --script "mask_C=maskutil_T; " maskutil_T.nc masks_${CplModel}.nc 
     372 
     373# rm grids_${CplModel}_tmp.nc areas_${CplModel}_tmp.nc maskutil_T.nc 
    373374 
    374375# 
    375376echo ${Titre}"Final renaming"${Norm} 
    376377# ---------------------------------------------------------------------------- 
    377 for OCEGRID in T U V O C 
     378for OCEGRID in T U V O 
    378379do 
    379380    ocegrid=${OCEGRID~} # To lowercase 
    … …  
    382383    ncrename --history --variable bounds_lon_grid_${OCEGRID},${ocegrid}${oce}.clo  grids_${CplModel}.nc 
    383384    ncrename --history --variable bounds_lat_grid_${OCEGRID},${ocegrid}${oce}.cla  grids_${CplModel}.nc 
    384     ncrename --history --variable mask_${OCEGRID},${ocegrid}${oce}.msk             masks_${CplModel}.nc 
     385    ncrename --history --variable maskutil_${OCEGRID},${ocegrid}${oce}.msk         masks_${CplModel}.nc 
    385386    ncrename --history --variable area_grid_${OCEGRID},${ocegrid}${oce}.srf        areas_${CplModel}.nc 
    386387     
    … …  
    399400ncrename --history --variable bounds_lat,t${atm}.cla          grids_${CplModel}.nc 
    400401 
    401 for ATMGRID in O C ; do 
     402for ATMGRID in O ; do 
    402403    atmgrid=${ATMGRID~} # To lowercase 
    403404    ncrename --history --variable ${atmgrid}${atm}_lon,${atmgrid}${atm}.lon         grids_${CplModel}.nc 
    … …  
    409410done 
    410411     
    411 for ATMGRID in T O C ; do 
     412for ATMGRID in T O ; do 
    412413    atmgrid=${ATMGRID~} # To lowercase 
    413414    ncatted --history \ 

  • TOOLS/MOSAIX/CreateWeightsMask.bash

    r6045 r6064  
    6868echo ${Titre}"Defines model"${Norm} 
    6969# ================================= 
    70 CplModel=ORCA2.3xLMD9695 
    71 #CplModel=ORCA2.3xICO30 
     70#CplModel=ORCA2.3xLMD9695 
     71CplModel=ORCA2.3xICO30 
    7272#CplModel=ORCA2.3xICO40 
    7373#CplModel=eORCA1.2xLMD144142 
    … …  
    7979#Version="v0" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" 
    8080#Version="v1" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" 
    81 Version="NoSearchRadius" ; Comment="For testing new routing" 
     81#Version="NoSearchRadius" ; Comment="For testing new routing" 
     82Version="v2" ;  Comment="Correction of ORCA masks to have a perfect conservation of run-off" 
    8283 
    8384# If available, get model name from job name 
    … …  
    104105 
    105106# More specific cases 
    106 [[ ${CplModel} = eORCA1.2xLMD144142 ]] && { atmCoastWidth=2 ; oceCoastWidth=2 ; searchRadius=550.0 ; Version="v1" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" ; } 
     107[[ ${CplModel} = eORCA1.2xLMD144142 ]] && { atmCoastWidth=2 ; oceCoastWidth=2 ; searchRadius=550.0 ; } 
    107108 
    108109if [[ ${CplModel} = eORCA1.2xLMD144142 && Version = "NoSearchRadius" ]] ; then 
    … …  
    178179## 
    179180## =========================================================================== 
    180 Tag="MOSAIX" 
     181#Tag="MOSAIX" 
    181182SUBMIT_DIR=$(pwd) 
    182183FMT_XIOS=netcdf4 
    … …  
    589590            --attribute NCO,global,o,c,"NCO netCDF Operator ${NCO} http://nco.sourceforge.net" \ 
    590591            --attribute Python,global,o,c,"Python3 version ${PYTHON_VER}"          \ 
    591             --attribute OS,global,o,c,"$(uname -o)"                               \ 
    592592            --attribute release,global,o,c,"$(uname -r)"                          \ 
    593593            --attribute directory,global,o,c,"$(pwd)"                             \ 
    … …  
    800800              --atmQuantity=${runOff_atmQuantity} --oceQuantity=${runOff_oceQuantity} --ocePerio=${OcePerio} 
    801801fi 
    802  
     802exit 
    803803## 
    804804echo ${Titre}"Calving weights"${Norm} 
    … …  
    909909NCO             : NCO netCDF Operator ${NCO} http://nco.sourceforge.net 
    910910Python version  : ${PYTHON_VER} 
    911 OS              : $(uname -o) 
    912911release         : $(uname -r) 
    913 hardware        : $(uname -i) 
    914912EOF 
    915913 

  • TOOLS/MOSAIX/RunoffWeights.py

    r6042 r6064  
    149149    if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 
    150150print ("Oce NPERIO parameter : {:d}".format(oce_perio)) 
    151 oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask, (oce_jpj,oce_jpi)), 'T', oce_perio).ravel() 
     151oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() 
    152152oce_address = np.arange(oce_jpj*oce_jpi) 
    153153 
    154154print ("Fill closed sea with image processing library") 
    155155oce_grid_imask2D = np.reshape(oce_grid_pmask,(oce_jpj,oce_jpi)) 
    156 oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T' ) 
     156oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T', sval=0 ) 
    157157oce_grid_imask = oce_grid_imask2D.ravel() 
    158158##  
    … …  
    165165oceOceanFiltered2D = ndimage.uniform_filter(oceOcean2D.astype(float), size=NNocean) 
    166166oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D 
    167 oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D,(oce_jpj,oce_jpi)), oce_perio, 'T').ravel() 
     167oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D,(oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T').ravel() 
    168168 
    169169oceOceanFiltered = oceOceanFiltered2D.ravel() 

  • TOOLS/MOSAIX/nemo.py

    r5948 r6064  
    1818''' 
    1919 
    20 import sys 
    21  
    22 try    : import numpy as np 
    23 except : pass 
    24  
     20import sys, numpy as np 
    2521try    : import xarray as xr 
    2622except : pass 
     23 
     24rpi = np.pi ; rad = rpi / 180.0 
     25 
     26nperio_valid_range = [0,1,4,5,6] 
    2727 
    2828## SVN information 
    … …  
    3333__HeadURL    = "$HeadURL$" 
    3434 
    35 def __guessNperio__ (jpi, nperio=None) : 
     35def __guessNperio__ (jpj, jpi, nperio=None) : 
    3636    ''' 
    3737    Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details) 
    … …  
    4141    nperio : periodicity parameter 
    4242    ''' 
    43     
    4443    if nperio == None : 
    4544        if jpi ==  182 : nperio = 4 #   ORCA2. We choose legacy orca2. 
    46         if jpi ==  362 : nperio = 6 #   ORCA1. 
    47         if jpi == 1442 : nperio = 6 # ORCA025. 
     45        if jpi ==  362 : nperio = 6 #   eORCA1. 
     46        if jpi == 1442 : nperio = 6 #   ORCA025. 
     47 
     48        if jpj ==  149 : nperio = 4 #   ORCA2. We choose legacy orca2. 
     49        if jpj ==  332 : nperio = 6 #   eORCA1. 
    4850        # 
    4951        if nperio == None : 
    50             sys.exit ('in nemo.lbc : nperio not found, and cannot by guessed' ) 
     52            sys.exit ('in nemo module : nperio not found, and cannot by guessed' ) 
    5153        else : 
    52             print ('nperio set as {:d} (deduced from jpi={:d})'.format(nperio, jpi)) 
     54            print ('nperio set as {:d} (deduced from jpi={:d})'.format (nperio, jpi)) 
    5355 
    5456    return nperio 
     57 
     58def __guessPoint__ (ptab) : 
     59    ''' 
     60    Tries to guess the grid point (periodicity parameter. See NEMO documentation for details) 
     61     
     62    For array conforments with xgcm requirements 
     63 
     64    Inputs 
     65    ptab : xarray array 
     66 
     67    Credits : who is the original author ? 
     68    ''' 
     69    gP = None 
     70    if isinstance (ptab, xr.core.dataarray.DataArray) : 
     71        if 'x_c' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'T' 
     72        if 'x_f' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'U' 
     73        if 'x_c' in ptab.dims and 'y_f' in ptab.dims                        : gP = 'V' 
     74        if 'x_f' in ptab.dims and 'y_f' in ptab.dims                        : gP = 'F' 
     75        if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_c' in ptab.dims : gP = 'T' 
     76        if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'W' 
     77        if 'x_f' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'U' 
     78        if 'x_c' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'V' 
     79        if 'x_f' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'F' 
     80              
     81        if gP == None : 
     82            sys.exit ('in nemo module : cd_type not found, and cannot by guessed' ) 
     83        else : 
     84            print ('Grid set as', gP, 'deduced from dims ', ptab.dims) 
     85            return gP 
     86    else : 
     87         sys.exit ('in nemo module : cd_type not found, input is not an xarray data' ) 
     88         #return gP 
    5589 
    5690def fixed_lon (lon) : 
    … …  
    6094    fixed_lon (lon) 
    6195    lon : longitudes of the grid. At least 2D. 
     96 
     97    Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7 
    6298    ''' 
    6399    fixed_lon = lon.copy () 
    64100    for i, start in enumerate (np.argmax (np.abs (np.diff (lon, axis=-1)) > 180, axis=-1)) : 
    65         fixed_lon[...,i, start+1:] += 360 
     101        fixed_lon[..., i, start+1:] += 360 
     102 
     103    # Special case for eORCA025 
     104    if lon.shape[-1] == 1442 : lon [..., -2, :] =  lon [..., -3, :] 
     105         
    66106    return fixed_lon 
    67107 
    … …  
    82122    lat (optionnal) : latitudes of the grid 
    83123    ''' 
    84      
    85124    if np.max (lat) != None : 
    86125        je     = jeq (lat) 
    … …  
    95134    return lon1D 
    96135 
    97 def latreg (lat) : 
    98     ''' 
    99     Returns maximum j index wehre gridlines are along latitude in the norethern hemisphere 
    100      
    101     lon : longitudes of the grid 
    102     lat (optionnal) : latitudes of the grid 
    103     ''' 
    104     dy     = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2)))) 
     136def latreg (lat, diff=0.1) : 
     137    ''' 
     138    Returns maximum j index where gridlines are along latitude in the northern hemisphere 
     139     
     140    lat : latitudes of the grid (2D) 
     141    diff [optional] : tolerance 
     142    ''' 
     143    if diff == None : 
     144        dy = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2)))) 
     145        diff = dy/100. 
     146     
    105147    je     = jeq (lat) 
    106     jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< dy/100.)[-1][-1] + je 
     148    jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< diff)[-1][-1] + je 
    107149    latreg = np.float64 (lat[...,jreg,:].mean(axis=-1)) 
    108     JREG = jreg 
     150    JREG   = jreg 
    109151 
    110152    return jreg, latreg 
    … …  
    114156    Returns 1D latitudes for zonal means and simple plots. 
    115157 
    116     lat : latitudes of the grid 
     158    lat : latitudes of the grid (2D) 
    117159    ''' 
    118160    jpj, jpi = lat.shape[-2:] 
    119161    try : 
    120         if type (lat) ==  xr.core.dataarray.DataArray : math = xr 
     162        if type (lat) == xr.core.dataarray.DataArray : math = xr 
    121163        else : math = np 
    122164    except : math = np 
    123165 
    124     # jreg : index of last uniform latitude, north of equator 
    125  
    126     dy = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2)))) 
     166    dy     = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2)))) 
    127167    je     = jeq (lat) 
    128     lateq = np.float64 (lat[...,je,:].mean(axis=-1)) 
     168    lat_eq = np.float64 (lat[...,je,:].mean(axis=-1)) 
    129169       
    130     jreg   = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< dy/100.)[-1][-1] + je 
    131     latreg = np.float64 (lat[...,jreg,:].mean(axis=-1)) 
    132     latave = np.mean (lat, axis=-1) 
    133  
    134     if ( np.abs (lateq) <  dy/100. ) : # T, U or W grid 
    135         dys = (90.-latreg)//(jpj-jreg-1)*0.5 
    136         yrange = 90.-dys-latreg 
    137     else                          :  # V or F grid 
    138         yrange = 90.    -latreg 
    139          
    140     lat1D = math.where (latave<latreg, latave, latreg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) )    
    141          
    142     if math == xr : 
    143         lat1D.attrs = lat.attrs 
     170    jreg, lat_reg = latreg (lat) 
     171    lat_ave = np.mean (lat, axis=-1) 
     172 
     173    if ( np.abs (lat_eq) < dy/100. ) : # T, U or W grid 
     174        dys    = (90.-lat_reg) / (jpj-jreg-1)*0.5 
     175        yrange = 90.-dys-lat_reg 
     176    else                             :  # V or F grid 
     177        yrange = 90.    -lat_reg 
     178         
     179    lat1D = math.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) )    
     180         
     181    if math == xr : lat1D.attrs = lat.attrs 
    144182 
    145183    return lat1D 
    … …  
    148186    ''' 
    149187    Returns simple latitude and longitude (1D) for simple plots. 
    150     ''' 
    151     laty = lat1D (lat) 
    152     lonx = lon1D (lon, lat) 
    153  
    154     return laty, lonx 
    155          
    156 def extend (tab, Lon=False, jplus=25) : 
    157     ''' 
    158     Returns extended field eastward to have better plots 
     188 
     189    lat, lon : latitudes and longitudes of the grid (2D) 
     190    ''' 
     191    return lat1D (lat),  lon1D (lon, lat) 
     192 
     193def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) : 
     194    try : 
     195        lon = lon.to_masked_array() 
     196        lat = lat.to_masked_array() 
     197    except : pass 
     198             
     199    mask = np.logical_and ( np.logical_and(lat>y0, lat<y1),  
     200        np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)), 
     201                                         np.logical_and(lon-360>x0, lon-360<x1)) ) 
     202    try    : tab = xr.where (mask, ptab, np.nan) 
     203    except : tab = np.where (mask, ptab, np.nan) 
     204     
     205    return tab 
     206         
     207def extend (tab, Lon=False, jplus=25, jpi=None, nperio=4) : 
     208    ''' 
     209    Returns extended field eastward to have better plots, and box average crossing the boundary 
    159210    Works only for xarray and numpy data (?) 
    160211 
    161212    tab : field to extend. 
    162213    Lon : (optional, default=False) : if True, add 360 in the extended parts of the field 
     214    jpi : normal longitude dimension of the field. exrtend does nothing it the actual 
     215        size of the field != jpi (avoid to extend several times) 
    163216    jplus (optional, default=25) : number of points added on the east side of the field 
    164217     
    165218    ''' 
    166     jpi   = tab.shape[-1] 
    167     JPLUS = jplus 
    168      
    169     try :  
    170         if type (tab) ==  xr.core.dataarray.DataArray : math = xr 
    171         else : math = np 
    172     except : math = np 
    173  
    174     if Lon : xplus = -360.0 
    175     else   : xplus =    0.0 
    176  
    177     if math == xr :  
    178         extend = xr.concat      ( (tab[..., 2:jpi] + xplus, tab[..., 2:jplus]), dim=tab.dims[-1] ) 
    179     if math == np : 
    180         extend = np.concatenate ( (tab[..., 2:jpi] + xplus, tab[..., 2:jplus]), axis=-1 ) 
     219    if tab.shape[-1] == 1 : 
     220        extend = tab 
     221 
     222    else : 
     223        if jpi == None : jpi = tab.shape[-1] 
     224     
     225        try :  
     226            if type (tab) ==  xr.core.dataarray.DataArray : math = xr 
     227            else : math = np 
     228        except : math = np 
     229 
     230        if Lon : xplus = -360.0 
     231        else   : xplus =    0.0 
     232 
     233        if tab.shape[-1] > jpi : 
     234            extend = tab 
     235        else : 
     236            if nperio == 1 : 
     237                istart = 0 ; le=jpi+1 ; la=0 
     238            if nperio > 1 : # OPA case with to halo points for periodicity 
     239                istart = 1 ; le=jpi-2 ; la=1  # Perfect, except at the pole that should be masked by lbc_plot 
     240                 
     241            if math == xr :  
     242                extend = xr.concat      ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), dim=tab.dims[-1] ) 
     243            if math == np : 
     244                extend = np.concatenate ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), axis=-1 ) 
    181245             
    182246    return extend 
    183247 
     248def orca2reg (ff, lat_name='nav_lat', lon_name='nav_lon', y_name='y', x_name='x') : 
     249    ''' 
     250    Assign ORCA dataset on a regular grid. 
     251    For use in the tropical region. 
     252 
     253    Inputs :  
     254      ff : xarray dataset 
     255      lat_name, lon_name : name of latitude and longitude 2D field in ff 
     256      y_name, x_name     : namex of dimensions in ff 
     257 
     258    Returns : xarray dataset with rectangular grid. Incorrect above 20°N 
     259    ''' 
     260    # Compute 1D longitude and latitude 
     261    (lat, lon) = latlon1D (ff[lat_name], ff[lon_name]) 
     262    # Assign lon and lat as dimensions of the dataset 
     263    if y_name in ff.dims :  
     264        lat = xr.DataArray (lat, coords=[lat,], dims=['lat',])      
     265        ff  = ff.rename_dims ({y_name: "lat",}).assign_coords (lat=lat) 
     266    if x_name in ff.dims : 
     267        lon = xr.DataArray (lon, coords=[lon,], dims=['lon',]) 
     268        ff  = ff.rename_dims ({x_name: "lon",}).assign_coords (lon=lon) 
     269    # Force dimensions to be in the right order 
     270    coord_order = ['lat', 'lon'] 
     271    for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z', 
     272                 'time_counter', 'time', 'tbnds',  
     273                 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] : 
     274        if dim in ff.dims : coord_order.insert (0, dim ) 
     275         
     276    ff = ff.transpose (*coord_order) 
     277    return ff 
     278 
     279def lbc_init (ptab, nperio=None, cd_type='T') : 
     280    """ 
     281    Prepare for all lbc calls 
     282     
     283    Set periodicity on input field 
     284    nperio    : Type of periodicity 
     285      0       : No periodicity 
     286      1, 4, 6 : Cyclic on i dimension (generaly longitudes) with 2 points halo 
     287      2       : Obsolete (was symmetric condition at southern boundary ?) 
     288      3, 4    : North fold T-point pivot (legacy ORCA2) 
     289      5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo) 
     290    cd_type   : Grid specification : T, U, V or F 
     291 
     292    See NEMO documentation for further details 
     293    """ 
     294    jpj, jpi = ptab.shape[-2:] 
     295    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) 
     296     
     297    if nperio not in nperio_valid_range : 
     298        print ( 'nperio=', nperio, ' is not in the valid range', nperio_valid_range ) 
     299        sys.exit() 
     300 
     301    if cd_type == None and isinstance (ptab, xr.core.dataarray.DataArray) : 
     302        cd_type = __guessPoint__ (ptab) 
     303 
     304    if cd_type not in ['T', 'U', 'V', 'F', 'W', 'F' ]: 
     305        print ( 'cd_type=', cd_type, ' is not in the list T, U, V, F, W, F' ) 
     306        sys.exit () 
     307 
     308    return jpj, jpi, nperio, cd_type 
     309         
    184310         
    185311def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) : 
    … …  
    189315      rank 2 at least : ptab[...., lat, lon] 
    190316    nperio    : Type of periodicity 
    191       1, 4, 6 : Cyclic on i dimension (generaly longitudes) 
    192       2       : Obsolete (was symmetric condition at southern boundary ?) 
    193       3, 4    : North fold T-point pivot (legacy ORCA2) 
    194       5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo) 
    195317    cd_type   : Grid specification : T, U, V or F 
    196318    psgn      : For change of sign for vector components (1 for scalars, -1 for vector components) 
    … …  
    198320    See NEMO documentation for further details 
    199321    """ 
    200  
    201     jpi    = ptab.shape[-1] 
    202     nperio = __guessNperio__ ( jpi, nperio) 
    203     psgn   = ptab.dtype.type(psgn) 
    204      
     322    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) 
     323    psgn   = ptab.dtype.type (psgn) 
     324 
    205325    # 
    206326    #> East-West boundary conditions 
    … …  
    210330        ptab [..., :,  0] = ptab [..., :, -2] 
    211331        ptab [..., :, -1] = ptab [..., :,  1] 
    212     
    213332    # 
    214333    #> North-South boundary conditions 
    … …  
    216335    if nperio in [3, 4] :  # North fold T-point pivot 
    217336        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    218             ptab [..., -1, 1:       ] = psgn * ptab [..., -3, -1:0:-1      ]       
     337            ptab [..., -1, 1:       ] = psgn * ptab [..., -3, -1:0:-1] 
    219338            ptab [..., -1, 0        ] = psgn * ptab [..., -3, 2            ] 
    220339            ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2:0:-1  ] 
    221                        
     340                 
    222341        if cd_type == 'U' : 
    223342            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ]        
    224343            ptab [..., -1,  0       ] = psgn * ptab [..., -3,  1           ] 
    225344            ptab [..., -1, -1       ] = psgn * ptab [..., -3, -2           ] 
    226  
     345                 
    227346            if nemo_4U_bug : 
    228347                ptab [..., -2, jpi//2+1:-1] = psgn * ptab [..., -2, jpi//2-2:0:-1] 
    … …  
    230349            else : 
    231350                ptab [..., -2, jpi//2-1:-1] = psgn * ptab [..., -2, jpi//2:0:-1] 
    232              
     351                 
    233352        if cd_type == 'V' :  
    234353            ptab [..., -2, 1:       ] = psgn * ptab [..., -3, jpi-1:0:-1   ] 
    235354            ptab [..., -1, 1:       ] = psgn * ptab [..., -4, -1:0:-1      ]    
    236355            ptab [..., -1, 0        ] = psgn * ptab [..., -4, 2            ] 
    237              
     356                 
    238357        if cd_type == 'F' : 
    239358            ptab [..., -2, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ] 
    … …  
    241360            ptab [..., -1,  0       ] = psgn * ptab [..., -4,  1           ] 
    242361            ptab [..., -1, -1       ] = psgn * ptab [..., -4, -2           ]  
    243        
     362         
    244363    if nperio in [5, 6] :            #  North fold F-point pivot   
    245364        if cd_type in ['T', 'W']  : 
    … …  
    248367        if cd_type == 'U' : 
    249368            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -2, -2::-1       ]        
    250             ptab [..., -1, -1       ] = psgn * ptab [..., -2, 0            ] 
    251               
     369            ptab [..., -1, -1       ] = psgn * ptab [..., -2, 0            ] # Bug ? 
     370                 
    252371        if cd_type == 'V' : 
    253372            ptab [..., -1, 0:       ] = psgn * ptab [..., -3, -1::-1       ] 
    254373            ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2-1::-1 ] 
    255                               
     374                 
    256375        if cd_type == 'F' : 
    257376            ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -2::-1       ] 
    258377            ptab [..., -1, -1       ] = psgn * ptab [..., -3, 0            ] 
    259378            ptab [..., -2, jpi//2:-1] = psgn * ptab [..., -2, jpi//2-2::-1 ] 
    260  
     379                 
    261380    return ptab 
    262381 
    263 def lbc_mask (ptab, nperio=None, cd_type='T', sval=None) : 
     382def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) : 
    264383    # 
    265384    """ 
    … …  
    268387      rank 2 at least : ptab[...., lat, lon] 
    269388    nperio    : Type of periodicity 
    270       1, 4, 6 : Cyclic on i dimension (generaly longitudes) 
    271       2       : Obsolete (was symmetric condition at southern boundary ?) 
    272       3, 4    : North fold T-point pivot (legacy ORCA2) 
    273       5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo) 
    274389    cd_type   : Grid specification : T, U, V or F 
    275390     
    … …  
    277392    """ 
    278393 
    279     jpi    = ptab.shape[-1] 
    280     nperio = __guessNperio__ (jpi, nperio) 
    281     if sval == None : sval = ptab.dtype.type (0) 
    282  
     394    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) 
     395          
    283396    # 
    284397    #> East-West boundary conditions 
    … …  
    303416            ptab [..., -1, 0        ] = sval  
    304417            ptab [..., -2, jpi//2:  ] = sval 
    305                        
     418                 
    306419        if cd_type == 'U' : 
    307             ptab [..., -1, 0:-1     ] = sval   
    308             ptab [..., -1,  0       ] = sval 
    309             ptab [..., -1, -1       ] = sval  
    310             ptab [..., -2, jpi//2  :] = sval 
    311              
     420            ptab [..., -1,  :       ] = sval   
     421            ptab [..., -2, jpi//2:  ] = sval 
     422                 
    312423        if cd_type == 'V' : 
    313             ptab [..., -2, 1:       ] = sval 
    314             ptab [..., -1, 1:       ] = sval    
    315             ptab [..., -1, 0        ] = sval 
    316              
     424            ptab [..., -2, :       ] = sval 
     425            ptab [..., -1, :       ] = sval    
     426                 
    317427        if cd_type == 'F' : 
    318             ptab [..., -2, 0:-1     ] = sval 
    319             ptab [..., -1, 0:-1     ] = sval 
    320             ptab [..., -1,  0       ] = sval 
    321             ptab [..., -1, -1       ] = sval 
    322        
     428            ptab [..., -2, :     ] = sval 
     429            ptab [..., -1, :     ] = sval 
     430     
    323431    if nperio in [5, 6] :            #  North fold F-point pivot 
    324432        if cd_type in ['T', 'W']  : 
    … …  
    340448    return ptab 
    341449 
     450 
     451def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) : 
     452    """ 
     453    Set periodicity on input field, adapted for plotting for any cartopy projection 
     454    ptab      : Input array (works  
     455      rank 2 at least : ptab[...., lat, lon] 
     456    nperio    : Type of periodicity 
     457    cd_type   : Grid specification : T, U, V or F 
     458    psgn      : For change of sign for vector components (1 for scalars, -1 for vector components) 
     459     
     460    See NEMO documentation for further details 
     461    """ 
     462 
     463    jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) 
     464    psgn   = ptab.dtype.type (psgn) 
     465     
     466    # 
     467    #> East-West boundary conditions 
     468    # ------------------------------ 
     469    if nperio in [1, 4, 6] : 
     470        # ... cyclic 
     471        ptab [..., :,  0] = ptab [..., :, -2] 
     472        ptab [..., :, -1] = ptab [..., :,  1] 
     473         
     474    # 
     475    #> North-South boundary conditions 
     476    # -------------------------------- 
     477    if nperio in [3, 4] :  # North fold T-point pivot 
     478        if cd_type in [ 'T', 'W' ] : # T-, W-point 
     479            ptab [..., -1, :] = sval 
     480            ptab [..., -2, :jpi//2] = sval 
     481            #ptab [..., -2, :] = sval 
     482                       
     483        if cd_type == 'U' : 
     484            ptab [..., -1, : ] = sval 
     485 
     486        if cd_type == 'V' :  
     487            ptab [..., -2, : ] = sval 
     488            ptab [..., -1, : ] = sval 
     489             
     490        if cd_type == 'F' : 
     491            ptab [..., -2, : ] = sval 
     492            ptab [..., -1, : ] = sval 
     493       
     494    if nperio in [5, 6] :            #  North fold F-point pivot   
     495        if cd_type in ['T', 'W']  : 
     496            ptab [..., -1, : ] = sval 
     497                 
     498        if cd_type == 'U' : 
     499            ptab [..., -1, : ] = sval       
     500              
     501        if cd_type == 'V' : 
     502            ptab [..., -1, : ] = sval 
     503            ptab [..., -2, jpi//2:  ] = sval 
     504                              
     505        if cd_type == 'F' : 
     506            ptab [..., -1, : ] = sval 
     507            ptab [..., -2, jpi//2:-1] = sval 
     508 
     509    return ptab 
     510 
     511def geo2oce (pxx, pyy, pzz, glam, gphi) :  
     512    ''' 
     513    Change vector from geocentric to east/north 
     514 
     515    Input 
     516        pxx, pyy, pzz : components on the geoce,tric system 
     517        glam, gphi : longitue and latitude of the points 
     518 
     519    ''' 
     520    gsinlon = np.sin (rad * glam) 
     521    gcoslon = np.cos (rad * glam) 
     522    gsinlat = np.sin (rad * gphi) 
     523    gcoslat = np.cos (rad * gphi) 
     524           
     525    pte = - pxx * gsinlon            + pyy * gcoslon 
     526    ptn = - pxx * gcoslon * gsinlat  - pyy * gsinlon * gsinlat + pzz * gcoslat 
     527 
     528    return pte, ptn 
     529 
     530def oce2geo (pte, ptn, glam, gphi) : 
     531    ##---------------------------------------------------------------------- 
     532    ##                    ***  ROUTINE oce2geo  *** 
     533    ##       
     534    ## ** Purpose : 
     535    ## 
     536    ## ** Method  :   Change vector from east/north to geocentric 
     537    ## 
     538    ## History : 
     539    ##        #         (A. Caubel)  oce2geo - Original code 
     540    ##---------------------------------------------------------------------- 
     541    gsinlon = np.sin (rad * glam) 
     542    gcoslon = np.cos (rad * glam) 
     543    gsinlat = np.sin (rad * gphi) 
     544    gcoslat = np.cos (rad * gphi) 
     545 
     546    pxx = - pte * gsinlon - ptn * gcoslon * gsinlat 
     547    pyy =   pte * gcoslon - ptn * gsinlon * gsinlat 
     548    pzz =   ptn * gcoslat 
     549     
     550    return pxx, pyy, pzz 
     551 
    342552## =========================================================================== 
    343553## 
Note: See TracChangeset for help on using the changeset viewer.