Changeset 6042 for TOOLS


Ignore:
Timestamp:
01/13/22 10:45:31 (2 years ago)
Author:
omamce
Message:

O.M. : MOASIX.
RunoffWeights?.py : correction in writing of the output file
CreateWeightsMask?.bash : small changes

Location:
TOOLS/MOSAIX
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/MOSAIX/CreateWeightsMask.bash

    r6040 r6042  
    688688    cp add_varoce.nco add_varoce_$(basename ${rmpFile} .nc)_o_to_a.nco 
    689689    [[ ${o_to_a} = true ]] && ncap2 --history --append --script-file add_varoce_$(basename ${rmpFile} .nc)_o_to_a.nco ${OCE}_coordinates_mask.nc ${rmpFile} 
    690     sed -i -e "s/src_/dst_/g" add_varoce.nco 
     690    sed -i "" -e "s/src_/dst_/g" add_varoce.nco 
    691691    cp add_varoce.nco add_varoce_$(basename ${rmpFile} .nc)_a_to_o.nco 
    692692    [[ ${a_to_o} = true ]] && ncap2 --history --append --script-file add_varoce_$(basename ${rmpFile} .nc)_a_to_o.nco ${OCE}_coordinates_mask.nc ${rmpFile} 
    … …  
    725725        fi 
    726726                     
    727         sed -i -e "s/dst_/src_/g" add_varatm.nco 
     727        sed -i "" -e "s/dst_/src_/g" add_varatm.nco 
    728728        cp add_varatm.nco add_varatm_$(basename ${rmpFile} .nc)_a_to_o.nco 
    729729        if [[ ${a_to_o} = true ]] ; then 
    … …  
    769769        fi 
    770770         
    771         sed -i -e "s/dst/src/g" add_varatm.nco 
     771        sed -i "" -e "s/dst/src/g" add_varatm.nco 
    772772        cp add_varatm.nco add_varatm_$(basename ${rmpFile} .nc)_a_to_o.nco 
    773773        if [[ ${a_to_o} = true ]] ; then 

  • TOOLS/MOSAIX/RunoffWeights.py

    r5915 r6042  
    133133atm_address = np.arange(atm_jpj*atm_jpi) 
    134134 
    135  
    136135(oce_nvertex, oce_jpj, oce_jpi) = gridFile['torc.cla'][:].shape ; jpon=oce_jpj*oce_jpj 
    137136oce_grid_size = oce_jpj*oce_jpi 
    … …  
    149148    if oce_jpi ==  362 : oce_perio = 6 # ORCA 1 
    150149    if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 
    151 print(f"oce_perio = {oce_perio}") 
     150print ("Oce NPERIO parameter : {:d}".format(oce_perio)) 
    152151oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask, (oce_jpj,oce_jpi)), 'T', oce_perio).ravel() 
    153152oce_address = np.arange(oce_jpj*oce_jpi) 
    154153 
    155 ## Fill closed sea with image processing library 
     154print ("Fill closed sea with image processing library") 
    156155oce_grid_imask2D = np.reshape(oce_grid_pmask,(oce_jpj,oce_jpi)) 
    157156oce_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' ) 
    … …  
    173172oceCoast = oceCoast2D.ravel() 
    174173 
    175 print ('Number of points in oceLand  : ' + str(oceLand.sum()) ) 
    176 print ('Number of points in oceOcean : ' + str(oceOcean.sum()) ) 
    177 print ('Number of points in oceCoast : ' + str(oceCoast.sum()) ) 
     174print ('Number of points in oceLand  : {:8d}'.format (oceLand.sum()) ) 
     175print ('Number of points in oceOcean : {:8d}'.format (oceOcean.sum()) ) 
     176print ('Number of points in oceCoast : {:8d}'.format (oceCoast.sum()) ) 
    178177 
    179178# Arrays with coastal points only 
    … …  
    194193## For LMDZ only !! 
    195194if atmDomainType == 'rectilinear' : 
     195    print ("Extend coastal band " ) 
    196196    NNatm = 1+2*atmCoastWidth 
    197197    atmLand2D = np.reshape ( atmLand, ( atm_jpj, atm_jpi) ) 
    … …  
    203203    atmCoast = atmCoast2D.ravel() 
    204204 
    205     print ('Number of points in atmLand  : ' + str(atmLand.sum())  ) 
    206     print ('Number of points in atmOcean : ' + str(atmOcean.sum()) ) 
    207     print ('Number of points in atmCoast : ' + str(atmCoast.sum()) ) 
     205    print ('Number of points in atmLand  : {:8d}'.format (atmLand.sum())  ) 
     206    print ('Number of points in atmOcean : {:8d}'.format (atmOcean.sum()) ) 
     207    print ('Number of points in atmCoast : {:8d}'.format (atmCoast.sum()) ) 
    208208 
    209209else : 
    … …  
    225225 
    226226## Loop on atmosphere coastal points 
    227 for ja in np.arange(len(atmCoast_grid_pmask)) : 
    228     z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) 
    229     z_mask = np.where ( z_dist*ra < searchRadius, True, False) 
    230     num_links = np.int(z_mask.sum()) 
    231     if num_links == 0 : continue 
    232     z_area = oceCoast_grid_area[z_mask].sum() 
    233     poids = np.ones ((num_links),dtype=np.float64) / z_area 
    234     if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] 
    235     if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] 
    236     if  ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print 
    237         print ( 'ja:{:8d}, num_links:{:8d},  z_area:{:8.4e},  atm area:{:8.4e},  weights sum:{:8.4e}  '.format(ja, num_links, z_area, atm_grid_area[ja], poids.sum() ) ) 
    238     # 
    239     matrix_local = poids 
    240     atm_address_local = np.ones(num_links, dtype=np.int32 ) * atmCoast_address[ja] 
    241     # Address on destination grid 
    242     oce_address_local = oceCoast_address[z_mask] 
    243     # Append to global arrays 
    244     remap_matrix = np.append ( remap_matrix, matrix_local      ) 
    245     atm_address  = np.append ( atm_address , atm_address_local ) 
    246     oce_address  = np.append ( oce_address , oce_address_local ) 
    247  
    248 print ('End of loop') 
     227if searchRadius > 0. : 
     228    print ("Loop on atmosphere coastal points") 
     229    for ja in np.arange(len(atmCoast_grid_pmask)) : 
     230        z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) 
     231        z_mask = np.where (z_dist*ra < searchRadius, True, False) 
     232        num_links = int(z_mask.sum()) 
     233        if num_links == 0 : continue 
     234        z_area = oceCoast_grid_area[z_mask].sum() 
     235        poids = np.ones ((num_links),dtype=np.float64) / z_area 
     236        if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] 
     237        if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] 
     238        if  ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print 
     239            print ( 'ja:{:8d}, num_links:{:8d},  z_area:{:8.4e},  atm area:{:8.4e},  weights sum:{:8.4e}  '.format(ja, num_links, z_area, atm_grid_area[ja], poids.sum() ) ) 
     240        # 
     241        matrix_local = poids 
     242        atm_address_local = np.ones(num_links, dtype=np.int32 ) * atmCoast_address[ja] 
     243        # Address on destination grid 
     244        oce_address_local = oceCoast_address[z_mask] 
     245        # Append to global arrays 
     246        remap_matrix = np.append ( remap_matrix, matrix_local      ) 
     247        atm_address  = np.append ( atm_address , atm_address_local ) 
     248        oce_address  = np.append ( oce_address , oce_address_local ) 
     249 
     250    print ('End of loop') 
    249251 
    250252num_links = remap_matrix.shape[0] 
    251253 
    252 ### Output file 
     254print ("Write output file") 
    253255runoff = myargs.output 
    254256f_runoff = netCDF4.Dataset ( runoff, 'w', format=FmtNetcdf ) 
    … …  
    266268f_runoff.title           = runoff 
    267269f_runoff.Program         = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) 
    268 f_runoff.atmCoastWidth   = str(atmCoastWidth) + " grid points" 
    269 f_runoff.oceCoastWidth   = str(oceCoastWidth) + " grid points" 
    270 f_runoff.searchRadius    = str(searchRadius/1000.) + " km" 
     270f_runoff.atmCoastWidth   = "{:d} grid points".format(atmCoastWidth) 
     271f_runoff.oceCoastWidth   = "{:d} grid points".format(oceCoastWidth) 
     272f_runoff.searchRadius    = "{:.0f} km".format(searchRadius/1000.) 
    271273f_runoff.atmQuantity     = myargs.atmQuantity 
    272274f_runoff.oceQuantity     = myargs.oceQuantity 
    … …  
    333335v_atm_grid_center_lon[:] = atm_grid_center_lon[:] 
    334336v_atm_grid_center_lat[:] = atm_grid_center_lat[:] 
    335 v_atm_grid_corner_lon[:] = atm_grid_corner_lon 
    336 v_atm_grid_corner_lat[:] = atm_grid_corner_lat 
     337v_atm_grid_corner_lon[:] = np.transpose(atm_grid_corner_lon) 
     338v_atm_grid_corner_lat[:] = np.transpose(atm_grid_corner_lat) 
    337339v_atm_grid_area      [:] = atm_grid_area[:] 
    338340v_atm_grid_imask     [:] = atm_grid_imask[:] 
    … …  
    360362v_oce_grid_center_lon[:] = oce_grid_center_lon[:] 
    361363v_oce_grid_center_lat[:] = oce_grid_center_lat[:] 
    362 v_oce_grid_corner_lon[:] = oce_grid_corner_lon 
    363 v_oce_grid_corner_lat[:] = oce_grid_corner_lon 
     364v_oce_grid_corner_lon[:] = np.transpose(oce_grid_corner_lon) 
     365v_oce_grid_corner_lat[:] = np.transpose(oce_grid_corner_lat) 
    364366v_oce_grid_area      [:] = oce_grid_area[:] 
    365367v_oce_grid_imask     [:] = oce_grid_imask[:] 
    … …  
    402404    v_atmOceanFrac    = f_runoff.createVariable ( 'atmOceanFrac'   , 'i4', ('src_grid_size',) ) 
    403405     
    404     v_atmLand[:]         = atmLand 
    405     v_atmLandFrac[:]     = atmLandFrac 
    406     v_atmLandFiltered[:] = atmLandFiltered 
    407     v_atmFrac[:]         = atmFrac 
    408     v_atmOcean[:]        = atmOcean 
    409     v_atmOceanFrac[:]    = atmOceanFrac 
    410  
    411 v_atmCoast         = f_runoff.createVariable ( 'atmCoast'       , 'i4', ('src_grid_size',) )  
     406    v_atmLand[:]         = atmLand.ravel() 
     407    v_atmLandFrac[:]     = atmLandFrac.ravel() 
     408    v_atmLandFiltered[:] = atmLandFiltered.ravel() 
     409    v_atmFrac[:]         = atmFrac.ravel() 
     410    v_atmOcean[:]        = atmOcean.ravel() 
     411    v_atmOceanFrac[:]    = atmOceanFrac.ravel() 
     412 
     413v_atmCoast         = f_runoff.createVariable ( 'atmCoast'        , 'i4', ('src_grid_size',) )  
    412414v_atmCoast[:]      = atmCoast 
    413415 
    … …  
    416418v_oceOceanFiltered = f_runoff.createVariable ( 'oceOceanFiltered', 'f4', ('dst_grid_size',) ) 
    417419v_oceCoast         = f_runoff.createVariable ( 'oceCoast'        , 'i4', ('dst_grid_size',) ) 
    418  
    419 v_oceLand[:]      = oceLand 
    420 v_oceOcean[:]     = oceOcean 
    421 v_oceOceanFiltered[:]     = oceOceanFiltered 
    422 v_oceCoast[:]     = oceCoast 
     420  
     421v_oceLand[:]          = oceLand 
     422v_oceOcean[:]         = oceOcean 
     423v_oceOceanFiltered[:] = oceOceanFiltered 
     424v_oceCoast[:]         = oceCoast 
    423425 
    424426f_runoff.sync () 
    425427 
    426428##  
    427 f_runoff.close() 
    428  
    429 print ('The end') 
     429f_runoff.close () 
     430 
     431print ('That''s all folks !') 
    430432 
    431433## ====================================================================================== 
Note: See TracChangeset for help on using the changeset viewer.