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