Ignore:
Timestamp:
11/23/18 13:06:47 (6 years ago)
Author:
dubos
Message:

devel/Python : XIOS output on dual mesh

Location:
codes/icosagcm/devel/Python/test
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py

    r785 r787  
    3131parser.add_argument("--Davies_N1", type=int, default=3) 
    3232parser.add_argument("--Davies_N2", type=int, default=3) 
    33  
     33parser.add_argument("--llm", type=int, default=50) 
     34parser.add_argument("--nx", type=int, default=200) 
     35parser.add_argument("--ny", type=int, default=30) 
     36parser.add_argument("--dt", type=float, default=360., help='Time step in seconds') 
     37parser.add_argument("--hydrostatic", action='store_true') 
     38parser.add_argument("--betaplane", action='store_true') 
    3439 
    3540args = parser.parse_args() 
    … …  
    4449    a     = 6.371229e6 # Radius of the Earth (m) 
    4550    ptop  = 2000. 
    46     y0    = Ly*0.5 
     51    y0    = .5*Ly 
    4752    Cpd   = 1004.5 
    4853    p0    = 1e5 
    … …  
    5156 
    5257    omega  = 7.292e-5  # Angular velocity of the Earth (s^-1)     
    53     phi0   = 90.*np.pi/180.0 # Reference latitude North pi/4 (deg) 
     58    phi0   = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg) 
    5459    f0     = 2*omega*np.sin(phi0)  
    55     beta0  = 2*omega*np.cos(phi0)/a 
    56     fb     = 2*omega*np.sin(phi0) - y0*2*omega*np.cos(phi0)/a 
     60    beta0  = 2*omega*np.cos(phi0)/a if args.betaplane else 0. 
     61    fb     = f0 - beta0*y0 
    5762 
    5863    def Phi_xy(y): 
    … …  
    6469    def ulon(x,y,eta): 
    6570        u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b)))  
    66         u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2)) 
     71#        u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2)) 
    6772        return  u 
    6873 
    … …  
    7176    def p(eta): return p0*eta     # eta  = p/p_s 
    7277 
    73     def eta(alpha) : return (1-(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels 
    74     def coriolis(x,y): return f0+beta0*y 
     78    def eta(alpha) : return (1.-(lap*ztop*alpha/T0))**(g/(Rd*lap)) # roughly equispaced levels 
     79    def coriolis(x,y): return f0+beta0*y # here y is 0 at domain center 
    7580 
    7681    filename = 'cart_%03d_%03d.nc'%(nx,ny) 
    … …  
    8287    mesh  = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 
    8388 
    84  
    8589    alpha_k = (np.arange(llm)  +.5)/llm 
    8690    alpha_l = (np.arange(llm+1)+ 0.)/llm 
    8791    x_ik, alpha_ik = np.meshgrid(mesh.lon_i,       alpha_k, indexing='ij') 
    88     y_ik, alpha_ik = np.meshgrid(mesh.lat_i+.5*Ly, alpha_k, indexing='ij') 
     92    y_ik, alpha_ik = np.meshgrid(mesh.lat_i+y0, alpha_k, indexing='ij') 
    8993    x_il, alpha_il = np.meshgrid(mesh.lon_i,       alpha_l, indexing='ij') 
    90     y_il, alpha_il = np.meshgrid(mesh.lat_i+.5*Ly, alpha_l, indexing='ij') 
     94    y_il, alpha_il = np.meshgrid(mesh.lat_i+y0, alpha_l, indexing='ij') 
    9195    x_ek, alpha_ek = np.meshgrid(mesh.lon_e,       alpha_k, indexing='ij') 
    92     y_ek, alpha_ek = np.meshgrid(mesh.lat_e+.5*Ly, alpha_k, indexing='ij') 
     96    y_ek, alpha_ek = np.meshgrid(mesh.lat_e+y0, alpha_k, indexing='ij') 
    9397 
    9498    print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g)))  
    … …  
    103107    Phi_il = Phi_xyeta(y_il, eta_il) 
    104108    Phi_ik = Phi_xyeta(y_ik, eta_ik) 
    105     p_ik = p(eta_ik) 
     109    p_ik, p_il = p(eta_ik), p(eta_il) 
    106110    T_ik = T(y_ik, eta_ik)  #ik full level(40), il(41) 
    107111 
    … …  
    110114    for l in range(llm):  
    111115        mass_ik[:,l]=(Phi_il[:,l+1]-Phi_il[:,l])/(g*gas.v[:,l]) 
     116#        mass_ik[:,l]=(p_il[:,l]-p_il[:,l+1])/g 
    112117    Sik, Wil  = gas.s*mass_ik, mesh.field_w() 
    113118     
    … …  
    224229 
    225230    g, Lx, Ly = 9.81, 4e7, 6e6 
    226     nx, ny, llm = 200, 30, 30 
     231    nx, ny, llm = args.nx, args.ny, args.llm 
    227232    dx,dy = Lx/nx,Ly/ny 
    228233 
    … …  
    236241 
    237242    T  = args.T 
    238     dt = 100.#180.#360. 
     243    dt = args.dt 
    239244    dz = flow0[3].max()/(params.g*llm) 
    240245    nt = int(math.ceil(T/dt)) 
    … …  
    242247    logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)) 
    243248     
     249    dt=0. # FIXME 
     250 
    244251    caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 
    245252 
    … …  
    252259    else: # time stepping in Fortran 
    253260        scheme = time_step.ARK2(None, dt) 
    254         caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 
     261        if args.hydrostatic: 
     262            caldyn_step = unst.caldyn_step_HPE(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 
     263        else: 
     264            caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 
    255265        davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) 
    256266        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() ) 
    … …  
    292302    z=mesh.field_mass() 
    293303 
    294  
    295304    #T, nslice, dt = 3600., 1, 3600. 
    296305    Nslice = args.N 
    … …  
    305314            # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W 
    306315            gas, w, z = diagnose(Phi,S,m,W) 
     316            ps = caldyn_step.ps 
    307317            curl_vk = curl(mesh, u) 
     318            du_fast, du_slow = caldyn_step.du_fast[0,:,:], caldyn_step.du_slow[0,:,:]  
     319            div_fast, div_slow = div(mesh,du_fast), div(mesh,du_slow) 
    308320            # write to disk 
    309             context.send_field_primal('ps', davies.beta_i) 
     321            context.send_field_primal('ps', ps) 
    310322            context.send_field_primal('temp', gas.T) 
    311323            context.send_field_primal('p', gas.p) 
    … …  
    313325            context.send_field_primal('uz', w) 
    314326            context.send_field_primal('z', z) 
    315             context.send_field_primal('curl', curl_vk) 
     327            context.send_field_primal('div_fast', div_fast) 
     328            context.send_field_primal('div_slow', div_slow) 
     329            print curl_vk.size 
     330            context.send_field_dual('curl', curl_vk) 
    316331 
    317332            print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')) 

  • codes/icosagcm/devel/Python/test/xml/context_icosa_simple.xml

    r622 r787  
    2020 
    2121     <domain_group id="i"> 
    22        <domain id="i" name="mesh"/> 
     22       <domain id="i" name="primal"/> 
     23     </domain_group> 
     24 
     25     <domain_group id="v"> 
     26       <domain id="v" name="dual"/> 
    2327     </domain_group> 
    2428 

  • codes/icosagcm/devel/Python/test/xml/field_def_simple.xml

    r785 r787  
    99    <field id="mid_ap" axis_ref="lev"   long_name="hybrid A coefficient at midpoints" /> 
    1010    <field id="mid_bp" axis_ref="lev"  long_name="hybrid B coefficient at midpoints" /> 
     11 
     12    <field_group domain_ref="v"> 
     13      <field_group axis_ref="lev">  
     14        <field id="curl" /> 
     15      </field_group> 
     16    </field_group> 
    1117     
    1218    <field_group domain_ref="i"> 
    … …  
    2329        <field id="z" /> 
    2430        <field id="theta" /> 
    25         <field id="curl" /> 
     31        <field id="div_fast" /> 
     32        <field id="div_slow" /> 
    2633      </field_group> 
    2734      

  • codes/icosagcm/devel/Python/test/xml/filedef_simple.xml

    r785 r787  
    2121      <field field_ref="theta" /> 
    2222      <field field_ref="curl" /> 
     23      <field field_ref="div_fast" /> 
     24      <field field_ref="div_slow" /> 
    2325      <field field_ref="phis" operation="once" freq_offset="0ts" name="PHIS"   standard_name="surface_geopotential" long_name="Surface geopotential" unit="m2/m2"/>       
    2426    </field_group> 
Note: See TracChangeset for help on using the changeset viewer.