Changeset 771 for codes/icosagcm
- Timestamp:
- 11/08/18 17:29:34 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r769 r771 18 18 parser = argparse.ArgumentParser() 19 19 20 parser.add_argument("--mpi_ni", type=int, 21 default=64, choices=None, 20 parser.add_argument("--mpi_ni", type=int, default=1, 22 21 help="number of x processors") 23 parser.add_argument("--mpi_nj", type=int, 24 default=64, choices=None, 22 parser.add_argument("--mpi_nj", type=int, default=1, 25 23 help="number of y processors") 26 24 parser.add_argument("--T", type=float, default=5., 27 25 help="Length of time slice in seconds") 26 parser.add_argument("--Davies_N1", type=int, default=5) 27 parser.add_argument("--Davies_N2", type=int, default=5) 28 28 29 29 30 args = parser.parse_args() … … 43 44 44 45 omega = 7.292e-5 # Angular velocity of the Earth (s^-1) 45 phi0 = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg)46 phi0 = 90.*np.pi/180.0 # Reference latitude North pi/4 (deg) 46 47 f0 = 2*omega*np.sin(phi0) 47 48 beta0 = 2*omega*np.cos(phi0)/a … … 63 64 filename = 'cart_%03d_%03d.nc'%(nx,ny) 64 65 print 'Reading Cartesian mesh ...' 65 def coriolis( lon,lat): return f0+0.*lon66 def coriolis(x,y): return f0+beta0*(y+.5*Ly) 66 67 meshfile = meshes.DYNAMICO_Format(filename) 67 68 pmesh = meshes.Unstructured_PMesh(comm,meshfile) … … 131 132 return gas, w, z 132 133 134 class Davies: 135 def __init__(self,N1,N2,x_i,y_i,x_e,y_e): 136 self.N1, self.N2 = N1, N2 137 self.beta_i = self.mask(x_i,y_i) 138 self.beta_e = self.mask(x_e,y_e) 139 def mask0(self,c,c0,delta): # 1D building block for Davies relaxation 140 N1, N2 = self.N1, self.N2 141 N3=N1+N2 142 m = np.zeros(c.size) 143 144 for i in range(c.size): 145 ci=c[i] 146 m[i] = (1.+np.cos((ci-c0+N3*delta)*np.pi/(N2*delta)))/2.0 147 if ci < c0-N3*delta : m[i]=1. 148 if ci > c0-N1*delta : m[i]=0. 149 return m 150 def relax(self, llm, step, flow): 151 beta_i, beta_e = self.beta_i, self.beta_e 152 m,S,u,Phi,W=flow 153 for l in range(llm): 154 step.mass[:,l] = beta_i*step.mass[:,l] + (1.-beta_i)*m[:,l] 155 step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l] 156 step.u[:,l] = beta_e*step.u[:,l] + (1.-beta_e)*u[:,l] 157 for l in range(llm+1): 158 step.geopot[:,l] = beta_i*step.geopot[:,l] + (1.-beta_i)*Phi[:,l] 159 step.W[:,l] = beta_i*step.W[:,l] + (1.-beta_i)*W[:,l] 160 161 162 class myDavies(Davies): 163 def mask(self,x,y): 164 # return self.mask0(y,Ly,dy) 165 return self.mask0(y,-.5*Ly,dy)*self.mask0(-y,-.5*Ly,dy) 166 133 167 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator 134 168 comm = client.comm … … 165 199 else: # time stepping in Fortran 166 200 scheme = time_step.ARK2(None, dt) 167 caldyn_step = unst.caldyn_step_NH(mesh,scheme, nt, caldyn_thermo,caldyn_eta,168 thermo,params,params.g)201 caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 202 davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) 169 203 def next_flow(m,S,u,Phi,W): 170 204 # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 171 205 caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u 172 206 caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W 173 caldyn_step.next() 207 for i in range(nt): 208 caldyn_step.next() 209 davies.relax(llm, caldyn_step, flow0) 174 210 return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), 175 caldyn_step.geopot.copy(), caldyn_step.W.copy())211 caldyn_step.geopot.copy(), caldyn_step.W.copy()) 176 212 177 213 m,S,u,Phi,W=flow0 … … 186 222 w=mesh.field_mass() 187 223 z=mesh.field_mass() 188 xx,yy = mesh.lat_i, mesh.lon_i189 224 190 225
Note: See TracChangeset
for help on using the changeset viewer.