Back to home page

MITgcm

 
 

    


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 # simple python script to generate input data
                0002 # this script is very similar to $ROOTDIR/verification/offline_exf_seaice/input/gendata.m
                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 # pure channel
                0028 namf='channel.bin'
                0029 depth=H0*np.ones((ny,nx),dtype='float64'); depth[0,:] = 0.
                0030 #if kwr > 0: writefield(namf,depth)
                0031 
                0032 # triangular obstacle
                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 # wind field
                0042 namf='windx.bin';
                0043 wnd=windx*np.ones((nt,ny,nx),dtype='float64')
                0044 if kwr > 0: writefield(namf,wnd)
                0045 
                0046 #- file name convention: "const_{xx}.bin" <-> uniform value = xx (in percent)
                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 #if kwr > 0: writefield(namf,var)
                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 #if kwr > 0: writefield(namf,iceConc)
                0074 
                0075 namf='ice0_heff.bin'; iceH0=0.2;
                0076 iceVol=iceConc*iceH0;
                0077 #if kwr > 0: writefield(namf,iceVol)
                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; #- dtx = amplitude of air temp variations in X-dir
                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; #- specific humid <--> 70.% relative humid
                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 #- salinity
                0108 sCst=30
                0109 so=sCst*np.ones((nt,ny,nx),dtype='float64')
                0110 namf='socn.bin'
                0111 #if kwr > 0: writefield(namf,so)
                0112 
                0113 muTf = 5.4e-2;
                0114 tfreeze=-muTf*sCst;
                0115 print 'T-freeze = {0:10.6f}'.format(tfreeze)
                0116 #- parabolic profile in Y, max @ j=4, min @ j=ny, amplitude=1.K
                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 #-- make some plots to check: ----------------
                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()