File indexing completed on 2018-03-02 18:45:48 UTC
view on githubraw file Latest commit 0ba3967d on 2014-06-24 08:00:16 UTC
0ba3967dec Mart*0001
0002
0003
0004 import numpy as np
0005
0006 kwr, kprt =1, 0
0007 nx, ny, nr, nt = 80, 42, 3, 1
0008
0009 xc = np.arange(nx)
0010 xc = xc - xc.mean()
0011 yc = np.arange(ny)+0.5
0012
0013
0014
0015 def writefield(fname,data):
0016 import sys
0017 print 'write to file: '+fname
0018 if sys.byteorder == 'little': data.byteswap(True)
0019 fid = open(fname,"wb")
0020 data.tofile(fid)
0021 fid.close()
0022
0023
0024 windx = 10.
0025 H0 = -100.
0026
0027
0028 namf='channel.bin'
0029 depth=H0*np.ones((ny,nx),dtype='float64'); depth[0,:] = 0.
0030
0031
0032
0033 namf='bathy_3c.bin';
0034 depth=H0*np.ones((ny,nx),dtype='float64'); depth[0,:]=0.;
0035 msk=np.tile(np.abs(xc),(ny,1)) + np.tile(yc,(nx,1)).T
0036 depth[msk < 24]=0.;
0037 y2d=np.tile(yc,(nx,1)).T
0038 depth[y2d > ny/2.]=H0;
0039 if kwr > 0: writefield(namf,depth)
0040
0041
0042 namf='windx.bin';
0043 wnd=windx*np.ones((nt,ny,nx),dtype='float64')
0044 if kwr > 0: writefield(namf,wnd)
0045
0046
0047 namf='const_00.bin';
0048 fld=np.zeros((nt,ny,nx),dtype='float64')
0049 if kwr > 0: writefield(namf,fld)
0050
0051 namf='const100.bin'; w0=1.;
0052 var=w0*np.ones((nt,ny,nx),dtype='float64')
0053 if kwr > 0: writefield(namf,var)
0054
0055 namf='const+20.bin'; w0=0.2;
0056 var=w0*np.ones((nt,ny,nx),dtype='float64')
0057
0058
0059 namf='heff_quartic.bin'
0060 hf_y = 8.*(yc/ny)**4
0061 hf_y.resize((1,ny,1))
0062 hf=np.tile(hf_y,(nt,1,nx))
0063 if kwr > 0: writefield(namf,hf)
0064
0065
0066
0067 namf='ice0_area.bin'; iceC0=1.;
0068 iceConc=iceC0*np.ones((ny,nx),dtype='float64'); iceConc[0,:]=0;
0069 iceConc[1,:] =0.00*iceC0;
0070 iceConc[2,:] =0.10*iceC0;
0071 iceConc[-1,:]=0.00*iceC0;
0072 iceConc[-2,:]=0.01*iceC0;
0073
0074
0075 namf='ice0_heff.bin'; iceH0=0.2;
0076 iceVol=iceConc*iceH0;
0077
0078
0079
0080
0081 dsw0=100;
0082 namf='dsw_'+str(dsw0)+'.bin'
0083 fld=dsw0*np.ones((nt,ny,nx),dtype='float64')
0084 if kwr > 0: writefield(namf,fld)
0085
0086 dlw0=250;
0087 namf='dlw_'+str(dlw0)+'.bin'
0088 fld=dlw0*np.ones((nt,ny,nx),dtype='float64')
0089 if kwr > 0: writefield(namf,fld)
0090
0091 cel2K=273.15; dtx=4;
0092 ta_x=cel2K + dtx*np.sin(np.pi*(1+2*xc/nx));
0093 ta=np.tile(ta_x,(nt,ny,1))
0094 namf='tair_'+str(dtx)+'x.bin'
0095 if kwr > 0: writefield(namf,ta)
0096
0097 cvapor_fac = 640380.000
0098 cvapor_exp = 5107.400
0099 atmrho = 1.200
0100 rh=70;
0101 tmpbulk = cvapor_fac*np.exp(-cvapor_exp/ta_x);
0102 qa_x = (rh/100.)*tmpbulk/atmrho
0103 qa=np.tile(qa_x,(nt,ny,1))
0104 namf='qa'+str(rh)+'_'+str(dtx)+'x.bin'
0105 if kwr > 0: writefield(namf,qa)
0106
0107
0108 sCst=30
0109 so=sCst*np.ones((nt,ny,nx),dtype='float64')
0110 namf='socn.bin'
0111
0112
0113 muTf = 5.4e-2;
0114 tfreeze=-muTf*sCst;
0115 print 'T-freeze = {0:10.6f}'.format(tfreeze)
0116
0117 to_y=(yc-3.5)/(ny-4)
0118 to_y=tfreeze+0.5-to_y*to_y
0119 mnV=to_y.min(); MxV=to_y.max(); Avr=to_y[1:].mean()
0120 print ' SST* av,mn,Mx: {0:9.6f} , {1:9.6f} , {2:9.6f} , {3:9.6f}'.format(Avr,mnV,MxV,MxV-mnV)
0121 to_y.resize((1,ny,1))
0122 to=np.tile(to_y,(nt,1,nx))
0123 namf='tocn.bin';
0124 if kwr > 0: writefield(namf,to)
0125
0126
0127
0128 import matplotlib.pyplot as plt
0129 hScal=np.asarray([-1.1, 0.1])*np.abs(H0);
0130 xg = xc-.5
0131 yg = yc-.5
0132 plt.figure(1); plt.clf()
0133 plt.subplot(211)
0134 var=depth
0135 plt.pcolormesh(xg,yg,var)
0136 plt.colorbar()
0137 plt.title('Depth [m]')
0138
0139 plt.subplot(413)
0140 var=depth
0141 j1=1
0142 j2=ny/2-1
0143 j3=j2+1
0144 plt.plot(xc,var[j1,:],'k-',label=str(j1))
0145 plt.plot(xc,var[j2,:],'ro-',label=str(j2))
0146 plt.plot(xc,var[j3,:],'b-',label=str(j3))
0147 plt.axis([-nx/2, nx/2, hScal[0], hScal[1]])
0148 plt.grid()
0149 plt.legend()
0150 plt.title('Depth @ j= cst');
0151
0152 plt.subplot(414);
0153 i=nx/2-1
0154 plt.plot(yc,var[:,i],'k-')
0155 plt.axis([0,ny,H0*1.1,-H0*.1]);
0156 plt.grid()
0157 plt.title(['Depth @ i=',str(i)]);
0158
0159
0160 dewPt=(qa_x*atmrho)/cvapor_fac;
0161 dewPt=-cvapor_exp/np.log(dewPt);
0162
0163 plt.figure(2);plt.clf();
0164 plt.subplot(311)
0165 plt.plot(xc,ta_x-cel2K,'r-',label='ta',linewidth = 1)
0166 plt.plot(xc,dewPt-cel2K,'b-',label='dew',linewidth = 1)
0167 plt.plot(xc,tfreeze*np.ones((nx,1)),'k-',linewidth = 1)
0168 AA=plt.axis()
0169 plt.axis([-nx/2, nx/2, AA[2], AA[3]]);
0170 plt.legend()
0171 plt.grid()
0172 plt.xlabel('X')
0173 plt.title('Air Temp ($^\circ$C): del-Temp-X = '+str(dtx)+' , RH= '+str(rh));
0174 plt.subplot(312)
0175 plt.plot(yc,to_y.flatten(),'b-',linewidth = 1)
0176 plt.plot(yc,tfreeze*np.ones((ny,1)),'k-',linewidth = 1)
0177 AA=plt.axis()
0178 plt.plot([1,1],AA[2:],'k',linewidth=2.)
0179 plt.axis([0, ny, AA[2], AA[3]])
0180 plt.grid()
0181 plt.xlabel('Y')
0182 plt.title('Ocean Temp $^\circ$C');
0183
0184 plt.subplot(313)
0185 var=iceConc[:,0]
0186 plt.semilogy(yc,var,'b-x',linewidth = 1,label='iceC')
0187 var=iceVol[:,0]
0188 plt.semilogy(yc,var,'r-x',linewidth = 1,label='hEff')
0189 AA=plt.axis()
0190 plt.plot([1,1],AA[2:],'k',linewidth=2.)
0191 plt.axis([0, ny, 0, 2*iceC0])
0192 plt.grid()
0193 plt.xlabel('Y')
0194 plt.legend(loc='lower center')
0195 plt.title('Initial ice in Channel : y-section');
0196
0197 if kprt == 1:
0198 f=2
0199 namfig='forcing_{0:02g}'.format(f)
0200 print ' print fig= {0:2g} to file: '.format(f) + namfig
0201 print(f,'-depsc2',namfig); fprintf('\n');
0202
0203
0204 plt.figure(3);plt.clf();
0205 plt.subplot(311)
0206 plt.pcolormesh(xg,yg,iceConc,vmin=-0.1, vmax=1.2)
0207 plt.colorbar()
0208 plt.title('Ice Concentration in Channel');
0209
0210 plt.subplot(312)
0211 plt.pcolormesh(xg,yg,iceVol,vmin=-1./50., vmax=12./50.)
0212 plt.colorbar()
0213 plt.title('Effective ice thickness in Channel');
0214
0215 plt.subplot(313)
0216 plt.semilogy(yc,iceConc[:,0],'b-x',label='iceC')
0217 plt.semilogy(yc,iceVol[:,0],'r-x',label='hEff')
0218 AA=plt.axis(); plt.axis([0,ny,0,2*iceC0]);
0219 plt.grid()
0220 plt.legend(loc='lower center')
0221 plt.title('Initial ice in Channel : y-section');
0222
0223 plt.show()
