source: NEMO/trunk/tests/ICE_AGRIF/EXPREF/make_INITICE.py @ 10516

Last change on this file since 10516 was 10516, checked in by clem, 19 months ago

change demo config name from ICEDYN to ICE_AGRIF

  • Property svn:executable set to *
File size: 3.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=1km
92#sigx=-0.04
93#sigy=-0.04
94#xshift=50.-1.
95#yshift=50.-1.
96#dlat=21
97#dlon=21
98
99# --- gaussian and square experiment ---
100#for y in np.arange(dlat,LAT1-dlat,1) :
101#    for x in np.arange(dlon,LON1-dlon,1) :
102#        cdfhti[:,y,x] = rn_hti_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2)
103#        cdfhts[:,y,x] = rn_hts_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2)
104#        cdfati[:,y,x] = rn_ati_ini
105#        cdfsmi[:,y,x] = rn_smi_ini
106#
107           
108# ----------------------------------------------
109# for basin=300x300km with dx=3km ; dy=3km + AGRIF
110sigx=-0.012
111sigy=-0.012
112xshift=20.-1.
113yshift=50.-1.
114dlat=18
115dlon=18
116for y in np.arange(32,66,1) :
117    for x in np.arange(2,36,1) :
118        cdfhti[:,y,x] = rn_hti_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2)
119        cdfhts[:,y,x] = rn_hts_ini*exp(sigx*(x-xshift)**2)*exp(sigy*(y-yshift)**2)
120        cdfati[:,y,x] = rn_ati_ini   
121        cdfsmi[:,y,x] = rn_smi_ini
122# ------------------------------------------------
123
124nc.close()
125nccoord.close()
126
127#sys.exit()
Note: See TracBrowser for help on using the repository browser.