[10515] | 1 | #!/usr/bin/python |
---|
| 2 | |
---|
| 3 | import os,sys |
---|
| 4 | from netCDF4 import Dataset as netcdf |
---|
| 5 | import numpy as np |
---|
| 6 | import matplotlib.pyplot as plt |
---|
| 7 | from math import exp |
---|
| 8 | from math import ceil |
---|
| 9 | |
---|
| 10 | resname='' |
---|
| 11 | |
---|
| 12 | # input file |
---|
| 13 | fcoord='mesh_mask.nc' |
---|
| 14 | |
---|
| 15 | # output file |
---|
| 16 | fflx='initice.nc' |
---|
| 17 | |
---|
| 18 | print ' creating init ice file ' +fflx |
---|
| 19 | |
---|
| 20 | # Reading coordinates file |
---|
| 21 | nccoord=netcdf(fcoord,'r') |
---|
| 22 | nav_lon=nccoord.variables['nav_lon'] |
---|
| 23 | nav_lat=nccoord.variables['nav_lat'] |
---|
| 24 | time_counter=1 |
---|
| 25 | LON1= nav_lon.shape[1] |
---|
| 26 | LAT1= nav_lon.shape[0] |
---|
| 27 | print 'nav_lon.shape[1]' ,nav_lon.shape[1] |
---|
| 28 | print 'LON1 ', LON1 |
---|
| 29 | print 'LAT1 ', LAT1 |
---|
| 30 | |
---|
| 31 | # Creating INITICE netcdf file |
---|
| 32 | nc=netcdf(fflx,'w') |
---|
| 33 | nc.createDimension('y',LAT1) |
---|
| 34 | nc.createDimension('x',LON1) |
---|
| 35 | nc.createDimension('time_counter',None) # Setting dimension size to 0 or None makes it unlimited. |
---|
| 36 | |
---|
| 37 | cdflon=nc.createVariable('nav_lon','f',('y','x')) |
---|
| 38 | cdflat=nc.createVariable('nav_lat','f',('y','x')) |
---|
| 39 | cdftimecounter=nc.createVariable('time_counter','f',('time_counter')) |
---|
| 40 | |
---|
| 41 | # ati : Fraction of open waters in sea ice - units % |
---|
| 42 | # hti : Sea ice thickness - units m |
---|
| 43 | # hts : Snow thickness - units m |
---|
| 44 | # smi : Sea ice salinity: |
---|
| 45 | # tmi : Sea ice internal temperature - units K |
---|
| 46 | # tsu : Sea ice surface temperature - units K |
---|
| 47 | # |
---|
| 48 | # Take constant values from namelist &namiceini of NEMO |
---|
| 49 | rn_hti_ini=2.0 |
---|
| 50 | rn_hts_ini=0.2 # initial real snow thickness (m) |
---|
| 51 | rn_ati_ini=0.9 # initial ice concentration (-) |
---|
| 52 | rn_smi_ini=6.3 # initial ice salinity (g/kg) |
---|
| 53 | rn_tmi_ini=270. # initial ice/snw temperature (K) |
---|
| 54 | rn_tsu_ini=270. # initial sea ice temperature (K) |
---|
| 55 | # |
---|
| 56 | cdfati=nc.createVariable('ati','f',('time_counter','y','x')) |
---|
| 57 | cdfati.units='Percentage' |
---|
| 58 | cdfati.long_name='Sea ice concentration' |
---|
| 59 | cdfhti=nc.createVariable('hti','f',('time_counter','y','x')) |
---|
| 60 | cdfhti.long_name='Sea ice thickness' |
---|
| 61 | cdfhti.units='m' |
---|
| 62 | cdfhts=nc.createVariable('hts','f',('time_counter','y','x')) |
---|
| 63 | cdfhts.long_name='Snow thickness' |
---|
| 64 | cdfhts.units='m' |
---|
| 65 | cdfsmi=nc.createVariable('smi','f',('time_counter','y','x')) |
---|
| 66 | cdfsmi.long_name='Sea ice salinity' |
---|
| 67 | cdfsmi.units='pss' |
---|
| 68 | cdftmi=nc.createVariable('tmi','f',('time_counter','y','x')) |
---|
| 69 | cdftmi.long_name='Sea ice internal temperature' |
---|
| 70 | cdftmi.units='Kelvin' |
---|
| 71 | cdftsu=nc.createVariable('tsu','f',('time_counter','y','x')) |
---|
| 72 | cdftsu.long_name='Sea ice surface temperature' |
---|
| 73 | cdftsu.units='Kelvin' |
---|
| 74 | |
---|
| 75 | cdflon[:,:]=nav_lon[:,:] |
---|
| 76 | cdflat[:,:]=nav_lat[:,:] |
---|
| 77 | cdftimecounter[0]=1 |
---|
| 78 | |
---|
| 79 | # Fill fields |
---|
| 80 | #print 'cdfati[:,1]', cdfati[:,1] -> 32 values |
---|
| 81 | |
---|
| 82 | # Add a gaussian for sea ice thickness here |
---|
| 83 | cdfhti[:,:,:]=0. |
---|
| 84 | cdfhts[:,:,:]=0. |
---|
| 85 | cdfati[:,:,:]=0. |
---|
| 86 | cdfsmi[:,:,:]=0. |
---|
| 87 | cdftmi[:,:,:]=rn_tmi_ini |
---|
| 88 | cdftsu[:,:,:]=rn_tsu_ini |
---|
| 89 | |
---|
| 90 | # -------------------------------------- |
---|
| 91 | # for basin=99x99km with dx=1km ; dy=3km |
---|
| 92 | #sigx=-0.04 |
---|
| 93 | #sigy=-0.04*9. |
---|
| 94 | #xshift=50.-1. |
---|
| 95 | #yshift=17.-1. |
---|
| 96 | #dlat=7 |
---|
| 97 | #dlon=21 |
---|
| 98 | |
---|
| 99 | # -------------------------------------- |
---|
| 100 | # for basin=99x99km with dx=3km ; dy=3km |
---|
| 101 | sigx=-0.04 |
---|
| 102 | sigy=-0.04 |
---|
| 103 | xshift=50.-1. |
---|
| 104 | yshift=50.-1. |
---|
| 105 | dlat=21 |
---|
| 106 | dlon=21 |
---|
| 107 | |
---|
| 108 | # --- gaussian and square experiment --- |
---|
| 109 | for y in np.arange(dlat,LAT1-dlat,1) : |
---|
| 110 | for x in np.arange(dlon,LON1-dlon,1) : |
---|
| 111 | cdfhti[:,y,x] = rn_hti_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2) |
---|
| 112 | cdfhts[:,y,x] = rn_hts_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2) |
---|
| 113 | cdfati[:,y,x] = rn_ati_ini |
---|
| 114 | cdfsmi[:,y,x] = rn_smi_ini |
---|
| 115 | |
---|
| 116 | nc.close() |
---|
| 117 | nccoord.close() |
---|
| 118 | |
---|
| 119 | #sys.exit() |
---|