source: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/make_INITICE.py @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 3 years ago

Merge of dev_CNRS_2017 into branch

  • Property svn:executable set to *
File size: 4.6 KB
Line 
1#!/usr/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
8from math import ceil
9
10resname=''
11
12# input file
13fcoord='mesh_mask.nc'
14
15# output file
16fflx='initice.nc'
17
18print '   creating init ice file  ' +fflx
19
20# Reading coordinates file
21nccoord=netcdf(fcoord,'r')
22nav_lon=nccoord.variables['nav_lon']
23nav_lat=nccoord.variables['nav_lat']
24time_counter=1
25LON1= nav_lon.shape[1]
26LAT1= nav_lon.shape[0]
27print 'nav_lon.shape[1]' ,nav_lon.shape[1]
28print 'LON1 ', LON1
29print 'LAT1 ', LAT1
30
31# Creating INITICE netcdf file
32nc=netcdf(fflx,'w')
33nc.createDimension('y',LAT1)
34nc.createDimension('x',LON1)
35nc.createDimension('time_counter',None)    # Setting dimension size to 0 or None makes it unlimited.
36
37cdflon=nc.createVariable('nav_lon','f',('y','x'))
38cdflat=nc.createVariable('nav_lat','f',('y','x'))
39cdftimecounter=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
49rn_hti_ini=2.0
50rn_hts_ini=0.2            #  initial real snow 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='Sea ice concentration'
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='Sea ice salinity'
67cdfsmi.units='pss'
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[:,:]
77cdftimecounter[0]=1
78
79# Fill fields
80#print 'cdfati[:,1]', cdfati[:,1]  -> 32 values
81
82# Add a gaussian for sea ice thickness here
83cdfhti[:,:,:]=0.
84cdfhts[:,:,:]=0.
85cdfati[:,:,:]=0.
86cdfsmi[:,:,:]=0.
87cdftmi[:,:,:]=rn_tmi_ini
88cdftsu[:,:,:]=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=1km ; dy=1km
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# --- Lipscomb 2004 experiment ---
117#cdfhti[:,:,:]=1.
118#cdfati[:,:,:]=0.1
119#for y in np.arange(0,LAT1,1) :
120#    for x in np.arange(0,LON1,1) :
121#        if (x > ceil(0.5*xshift) and x < xshift): # and (y > ceil(0.5*yshift) and y < yshift):
122#            cdfati[:,y,x] = rn_ati_ini / (ceil(0.5*xshift)*ceil(0.5*yshift)) * (x - ceil(0.5*xshift)) * (y - ceil(0.5*yshift))
123#            cdfsmi[:,y,x] = rn_smi_ini / (ceil(0.5*xshift)*ceil(0.5*yshift)) * (x - ceil(0.5*xshift)) * (y - ceil(0.5*yshift))
124#            cdfhti[:,y,x] = rn_hti_ini
125#            cdfhts[:,y,x] = rn_hts_ini
126#        if (x >= 25. and x < 50.):
127#            cdfati[:,y,x] = 0.9 * (x - 25.) / 25.
128#            cdfhti[:,y,x] = 1.
129#        elif (x >= 50. and x <= 75.):
130#            cdfati[:,y,x] = 0.9
131#            cdfhti[:,y,x] = 1.
132#        if (x > 30. and x < 70.):
133#            cdfhti[:,y,x] = 0.2
134           
135           
136# ----------------------------------------------
137# for basin=99x99km with dx=1km ; dy=1km + AGRIF
138sigx=-0.012
139sigy=-0.012
140xshift=20.-1.
141yshift=50.-1.
142dlat=18
143dlon=18
144for y in np.arange(32,66,1) :
145    for x in np.arange(2,36,1) :
146        cdfhti[:,y,x] = rn_hti_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2)
147        cdfhts[:,y,x] = rn_hts_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2)
148        cdfati[:,y,x] = rn_ati_ini   
149        cdfsmi[:,y,x] = rn_smi_ini
150# ------------------------------------------------
151
152nc.close()
153nccoord.close()
154
155#sys.exit()
Note: See TracBrowser for help on using the repository browser.