Ticket #1761: make_INITICE.py

File make_INITICE.py, 3.0 KB (added by clevy, 4 years ago)
Line 
1#!/opt/local/bin/python
2
3import os,sys
4from netCDF4 import Dataset as netcdf
5import numpy as np
6import matplotlib.pyplot as plt
7from math import exp
8
9resname=''
10
11# input file
12###fcoord='coordinates_'+str(resname)+'.nc'
13#fcoord='mesh_mask_'+str(resname)+'.nc'
14fcoord='mesh_mask.nc'
15
16# output file
17#fflx='initice_'+str(resname)+'.nc'
18fflx='initice.nc'
19
20print '   creating init ice file  ' +fflx
21
22# Reading coordinates file
23nccoord=netcdf(fcoord,'r')
24nav_lon=nccoord.variables['nav_lon']
25nav_lat=nccoord.variables['nav_lat']
26LON1= nav_lon.shape[1]
27LAT1= nav_lon.shape[0]
28print 'nav_lon.shape[1]' ,nav_lon.shape[1]
29print 'LON1 ', LON1
30print 'LAT1 ', LAT1
31
32# Creating INITICE netcdf file
33nc=netcdf(fflx,'w')
34nc.createDimension('y',LAT1)
35nc.createDimension('x',LON1)
36nc.createDimension('time_counter',0)    # Setting dimension size to 0 or None makes it unlimited.
37
38cdflon=nc.createVariable('nav_lon','f',('y','x'))
39cdflat=nc.createVariable('nav_lat','f',('y','x'))
40cdftimecounter=nc.createVariable('time_counter','f',('time_counter'))
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 :
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
49rn_hts_ini=0.3            #  initial real snow thickness (m)
50rn_hti_ini=3.0            #  initial real ice thickness  (m)
51rn_ati_ini=0.9            #  initial ice concentration   (-)
52rn_smi_ini=6.3            #  initial ice salinity     (g/kg)
53rn_tmi_ini=270.           #  initial ice/snw temperature (K)
54rn_tsu_ini=270.           #  initial sea ice temperature (K)
55#
56cdfati=nc.createVariable('ati','f',('time_counter','y','x'))
57cdfati.units='Percentage'
58cdfati.long_name='Fraction of open waters in sea ice'
59cdfhti=nc.createVariable('hti','f',('time_counter','y','x'))
60cdfhti.long_name='Sea ice thickness'
61cdfhti.units='m'
62cdfhts=nc.createVariable('hts','f',('time_counter','y','x'))
63cdfhts.long_name='Snow thickness'
64cdfhts.units='m'
65cdfsmi=nc.createVariable('smi','f',('time_counter','y','x'))
66cdfsmi.long_name=''
67cdfsmi.units=''
68cdftmi=nc.createVariable('tmi','f',('time_counter','y','x'))
69cdftmi.long_name='Sea ice internal temperature'
70cdftmi.units='Kelvin'
71cdftsu=nc.createVariable('tsu','f',('time_counter','y','x'))
72cdftsu.long_name='Sea ice surface temperature'
73cdftsu.units='Kelvin'
74
75cdflon[:,:]=nav_lon[:,:]
76cdflat[:,:]=nav_lat[:,:]
77
78# Fill fields
79#print 'cdfati[:,1]', cdfati[:,1]  -> 32 values
80## Add a gaussian for sea ice thickness here
81cdfhti[:,:,:]=0.
82A=7.     # max of gaussian for sea ice thichness in meters
83a=-0.1
84xshift=5.
85
86x=np.arange(0,LON1,1)
87y = [A*exp(a*t**2)for t in x-xshift]
88
89plt.plot(x,y,'ro')
90
91for t in np.arange(1,LON1,1) :
92    cdfhti[:,:,t] = A*exp(a*(t-xshift)**2)
93   
94#
95# Intialise all other arrays to constant values
96   
97cdfhts[:,:,:]=rn_hts_ini
98cdfati[:,:,:]=rn_ati_ini
99cdfsmi[:,:,:]=rn_smi_ini
100cdftmi[:,:,:]=rn_tmi_ini
101cdftsu[:,:,:]=rn_tsu_ini
102
103
104nc.close()
105nccoord.close()
106
107#sys.exit()