# -*- coding: euc-jp -*- # Title: Ruby script drawing contour map for deepconv/arare5 output data # # History: 2011/09/27 (Masatsugu Odaka) # require "numru/ggraph" module NumRu class GPhys def deriv(dim) dim = dim_index(dim) # input dim can be a String or an Integer x = axis(dim).to_gphys a = [true]*dim + [1..-1,false] b = [true]*dim + [0..-2,false] dydx = ( self[*a] - self[*b] ) / ( x[1..-1] - x[0..-2] ) xi = ( x[1..-1] + x[0..-2] ) / 2 dydx.axis(dim).set_pos(xi.coord(0)) dydx.long_name = "d #{name}/d #{x.name}" dydx end end end include NumRu filedir = "time_000180000-000280000" filehead = "BS1998_" cuttime = 180000..266400 puttime = 50000 zcut = 10000 kappa = 155.0 comax = 4.0e-3 comin = -4.0e-3 ins = 2 cpdry = 891.0 grav = 8.87 #comax =0.000010 #comin = -0.0000030 #file0 = 'BS1998_restart_rank000000.nc' dens_file = 'restart_long-DensBZ.nc' ptemp_file = 'restart_long-PTempBZ.nc' ptemp0_file = "../#{filedir}/#{filehead}PTemp_rank000000.nc" ptemp1_file = "../#{filedir}/#{filehead}PTemp_rank000001.nc" ptemp2_file = "../#{filedir}/#{filehead}PTemp_rank000002.nc" ptemp3_file = "../#{filedir}/#{filehead}PTemp_rank000003.nc" ptemp4_file = "../#{filedir}/#{filehead}PTemp_rank000004.nc" ptemp5_file = "../#{filedir}/#{filehead}PTemp_rank000005.nc" ptemp0_file = "../#{filedir}/#{filehead}PTemp_rank000000.nc" ptemp1_file = "../#{filedir}/#{filehead}PTemp_rank000001.nc" ptemp2_file = "../#{filedir}/#{filehead}PTemp_rank000002.nc" ptemp3_file = "../#{filedir}/#{filehead}PTemp_rank000003.nc" ptemp4_file = "../#{filedir}/#{filehead}PTemp_rank000004.nc" ptemp5_file = "../#{filedir}/#{filehead}PTemp_rank000005.nc" adv0_file = "../#{filedir}/#{filehead}PTempAdv_rank000000.nc" adv1_file = "../#{filedir}/#{filehead}PTempAdv_rank000001.nc" adv2_file = "../#{filedir}/#{filehead}PTempAdv_rank000002.nc" adv3_file = "../#{filedir}/#{filehead}PTempAdv_rank000003.nc" adv4_file = "../#{filedir}/#{filehead}PTempAdv_rank000004.nc" adv5_file = "../#{filedir}/#{filehead}PTempAdv_rank000005.nc" turb0_file = "../#{filedir}/#{filehead}PTempTurb_rank000000.nc" turb1_file = "../#{filedir}/#{filehead}PTempTurb_rank000001.nc" turb2_file = "../#{filedir}/#{filehead}PTempTurb_rank000002.nc" turb3_file = "../#{filedir}/#{filehead}PTempTurb_rank000003.nc" turb4_file = "../#{filedir}/#{filehead}PTempTurb_rank000004.nc" turb5_file = "../#{filedir}/#{filehead}PTempTurb_rank000005.nc" velxturb0_file = "../#{filedir}/#{filehead}VelXTurb_rank000000.nc" velxturb1_file = "../#{filedir}/#{filehead}VelXTurb_rank000001.nc" velxturb2_file = "../#{filedir}/#{filehead}VelXTurb_rank000002.nc" velxturb3_file = "../#{filedir}/#{filehead}VelXTurb_rank000003.nc" velxturb4_file = "../#{filedir}/#{filehead}VelXTurb_rank000004.nc" velxturb5_file = "../#{filedir}/#{filehead}VelXTurb_rank000005.nc" velzturb0_file = "../#{filedir}/#{filehead}VelZTurb_rank000000.nc" velzturb1_file = "../#{filedir}/#{filehead}VelZTurb_rank000001.nc" velzturb2_file = "../#{filedir}/#{filehead}VelZTurb_rank000002.nc" velzturb3_file = "../#{filedir}/#{filehead}VelZTurb_rank000003.nc" velzturb4_file = "../#{filedir}/#{filehead}VelZTurb_rank000004.nc" velzturb5_file = "../#{filedir}/#{filehead}VelZTurb_rank000005.nc" sfc0_file = "../#{filedir}/#{filehead}PTempSfc_rank000000.nc" sfc1_file = "../#{filedir}/#{filehead}PTempSfc_rank000001.nc" sfc2_file = "../#{filedir}/#{filehead}PTempSfc_rank000002.nc" sfc3_file = "../#{filedir}/#{filehead}PTempSfc_rank000003.nc" sfc4_file = "../#{filedir}/#{filehead}PTempSfc_rank000004.nc" sfc5_file = "../#{filedir}/#{filehead}PTempSfc_rank000005.nc" rad0_file = "../#{filedir}/#{filehead}PTempRad_rank000000.nc" rad1_file = "../#{filedir}/#{filehead}PTempRad_rank000001.nc" rad2_file = "../#{filedir}/#{filehead}PTempRad_rank000002.nc" rad3_file = "../#{filedir}/#{filehead}PTempRad_rank000003.nc" rad4_file = "../#{filedir}/#{filehead}PTempRad_rank000004.nc" rad5_file = "../#{filedir}/#{filehead}PTempRad_rank000005.nc" exner0_file = "../#{filedir}/#{filehead}Exner_rank000000.nc" exner1_file = "../#{filedir}/#{filehead}Exner_rank000001.nc" exner2_file = "../#{filedir}/#{filehead}Exner_rank000002.nc" exner3_file = "../#{filedir}/#{filehead}Exner_rank000003.nc" exner4_file = "../#{filedir}/#{filehead}Exner_rank000004.nc" exner5_file = "../#{filedir}/#{filehead}Exner_rank000005.nc" velx0_file = "../#{filedir}/#{filehead}VelX_rank000000.nc" velx1_file = "../#{filedir}/#{filehead}VelX_rank000001.nc" velx2_file = "../#{filedir}/#{filehead}VelX_rank000002.nc" velx3_file = "../#{filedir}/#{filehead}VelX_rank000003.nc" velx4_file = "../#{filedir}/#{filehead}VelX_rank000004.nc" velx5_file = "../#{filedir}/#{filehead}VelX_rank000005.nc" velz0_file = "../#{filedir}/#{filehead}VelZ_rank000000.nc" velz1_file = "../#{filedir}/#{filehead}VelZ_rank000001.nc" velz2_file = "../#{filedir}/#{filehead}VelZ_rank000002.nc" velz3_file = "../#{filedir}/#{filehead}VelZ_rank000003.nc" velz4_file = "../#{filedir}/#{filehead}VelZ_rank000004.nc" velz5_file = "../#{filedir}/#{filehead}VelZ_rank000005.nc" velzadv0_file = "../#{filedir}/#{filehead}VelZAdv_rank000000.nc" velzadv1_file = "../#{filedir}/#{filehead}VelZAdv_rank000001.nc" velzadv2_file = "../#{filedir}/#{filehead}VelZAdv_rank000002.nc" velzadv3_file = "../#{filedir}/#{filehead}VelZAdv_rank000003.nc" velzadv4_file = "../#{filedir}/#{filehead}VelZAdv_rank000004.nc" velzadv5_file = "../#{filedir}/#{filehead}VelZAdv_rank000005.nc" velxadv0_file = "../#{filedir}/#{filehead}VelXAdv_rank000000.nc" velxadv1_file = "../#{filedir}/#{filehead}VelXAdv_rank000001.nc" velxadv2_file = "../#{filedir}/#{filehead}VelXAdv_rank000002.nc" velxadv3_file = "../#{filedir}/#{filehead}VelXAdv_rank000003.nc" velxadv4_file = "../#{filedir}/#{filehead}VelXAdv_rank000004.nc" velxadv5_file = "../#{filedir}/#{filehead}VelXAdv_rank000005.nc" ptemp = GPhys::IO.open([ptemp0_file, \ ptemp1_file, \ ptemp2_file, \ ptemp3_file, \ ptemp4_file, \ ptemp5_file ], \ 'PTemp') adv = GPhys::IO.open([adv0_file, \ adv1_file, \ adv2_file, \ adv3_file, \ adv4_file, \ adv5_file ], \ 'PTempAdv') turb = GPhys::IO.open([turb0_file, \ turb1_file, \ turb2_file, \ turb3_file, \ turb4_file, \ turb5_file ], \ 'PTempTurb') velxturb = GPhys::IO.open([velxturb0_file, \ velxturb1_file, \ velxturb2_file, \ velxturb3_file, \ velxturb4_file, \ velxturb5_file ], \ 'VelXTurb') velzturb = GPhys::IO.open([velzturb0_file, \ velzturb1_file, \ velzturb2_file, \ velzturb3_file, \ velzturb4_file, \ velzturb5_file ], \ 'VelZTurb') velxadv = GPhys::IO.open([velxadv0_file, \ velxadv1_file, \ velxadv2_file, \ velxadv3_file, \ velxadv4_file, \ velxadv5_file ], \ 'VelXAdv') velzadv = GPhys::IO.open([velzadv0_file, \ velzadv1_file, \ velzadv2_file, \ velzadv3_file, \ velzadv4_file, \ velzadv5_file ], \ 'VelZAdv') rad = GPhys::IO.open([rad0_file, \ rad1_file, \ rad2_file, \ rad3_file, \ rad4_file, \ rad5_file ], \ 'PTempRad') exner = GPhys::IO.open([exner0_file, \ exner1_file, \ exner2_file, \ exner3_file, \ exner4_file, \ exner5_file ], \ 'Exner') sfc = GPhys::IO.open([sfc0_file, \ sfc1_file, \ sfc2_file, \ sfc3_file, \ sfc4_file, \ sfc5_file ], \ 'PTempSfc') velx = GPhys::IO.open([velx0_file, \ velx1_file, \ velx2_file, \ velx3_file, \ velx4_file, \ velx5_file ], \ 'VelX') velz = GPhys::IO.open([velz0_file, \ velz1_file, \ velz2_file, \ velz3_file, \ velz4_file, \ velz5_file ], \ 'VelZ') densBZ = GPhys::IO.open(dens_file, 'DensBZ') ptempBZ = GPhys::IO.open(ptemp_file, 'PTempBZ') turb = turb + sfc # dt_rho_K kin = ( velx**2.0 + velz**2.0 ) / 2.0 rho_K = densBZ * kin dt_rho_K = rho_K.deriv('t') dt_rho_K_XYZmean = dt_rho_K.mean('x').mean('y').mean('z') # rho_adv dx_K = kin.deriv('x') dz_K = kin.deriv('z') rho_adv = densBZ[1..1199,0..-1,1..99] * ( velx[1..1199,0..-1,1..99,0..-1] * dx_K[0..-1,0..-1,1..99,0..-1] \ + velz[1..1199,0..-1,1..99,0..-1] * dz_K[1..1199,0..-1,0..-1,0..-1] \ ) rho_adv_XYZmean = rho_adv.mean('x').mean('y').mean('z') # rho_adv = densBZ * (velxadv + velzadv) # rho_adv_XYZmean = rho_adv.mean('x').mean('y').mean('z') # press dx_exner = exner.deriv('x') dz_exner = exner.deriv('z') grad_exner = ( velx[1..1199,0..-1,1..99,0..-1] * dx_exner[0..-1,0..-1,1..99,0..-1] \ + velz[1..1199,0..-1,1..99,0..-1] * dz_exner[1..1199,0..-1,0..-1,0..-1] \ ) press_term = - cpdry * ptempBZ[1..1199,0..-1,1..99] * densBZ[1..1199,0..-1,1..99] * grad_exner press_XYZmean = press_term.mean('x').mean('y').mean('z') # Dterm =begin dxdx_u = velx.deriv('x').deriv('x') dzdz_u = velx.deriv('z').deriv('z') dxdz_w = velz.deriv('x').deriv('z') dif_x = kappa * ( \ 2.0 * dxdx_u[0..-1,0..-1,2..99,0..-1] \ + dzdz_u[2..1199,0..-1,0..-1,0..-1] \ + dxdz_w[1..1198,0..-1,1..98,0..-1] \ ) dxdx_w = velz.deriv('x').deriv('x') dxdz_u = velx.deriv('x').deriv('z') dz_w = velz.deriv('z') rho_dz_w = densBZ[0..-1,0..-1,1..99] * dz_w dz_rho_dz_w = rho_dz_w.deriv('z') dif_w = kappa * ( \ dxdx_w[0..-1,0..-1,2..99,0..-1] \ + dxdz_u[1..1198,0..-1,1..98,0..-1] \ ) + 2.0 * kappa * dz_rho_dz_w[2..1199,0..-1,0..-1,0..-1] / densBZ[2..1199,0..-1,2..99] dif = densBZ[2..1199,0..-1,2..99] * ( \ velx[2..1199,0..-1,2..99,0..-1] * dif_x \ + velx[2..1199,0..-1,2..99,0..-1] * dif_w \ ) dif_XYZmean = dif.mean('x').mean('y').mean('z') =end turb_x = densBZ * velx * velxturb turb_z = densBZ * velz * velzturb turb_term = turb_x + turb_z turb_XYZmean = turb_term.mean('x').mean('y').mean('z') # bouy bouy = densBZ * ptemp * grav * velz / ptempBZ bouy_XYZmean = bouy.mean('x').mean('y').mean('z') dif = bouy_XYZmean - rho_adv_XYZmean + press_XYZmean + turb_XYZmean # dif_XYZmean = dif.mean('x').mean('y').mean('z') DCL.gropn(ins) GGraph.set_axes('xtitle'=>'time','ytitle'=>'d Dens K / d t') GGraph.line( dt_rho_K_XYZmean.cut('t'=>cuttime), true, 'exchange'=>false , 'index'=>12, 'type'=>1, #'label'=>'rhoK', 'title'=>"Kinetic energy balance", 'max'=>comax, 'min'=>comin) GGraph.line( - rho_adv_XYZmean.cut('t'=>cuttime), false, 'exchange'=>false , # true, 'exchange'=>false , # 'title'=>"z=#{zcut}", 'index'=>22, 'type'=>1, #'label'=>'adv', 'max'=>comax, 'min'=>comin) GGraph.line( press_XYZmean.cut('t'=>cuttime), false, 'exchange'=>false , 'index'=>42, 'type'=>1, #'label'=>'exner', 'max'=>comax, 'min'=>comin) #GGraph.line( dif.cut('t'=>cuttime), # false, 'exchange'=>false , # 'index'=>62, 'type'=>1, #'label'=>'dif', # 'max'=>comax, 'min'=>comin) GGraph.line( turb_XYZmean.cut('t'=>cuttime), false, 'exchange'=>false , 'index'=>33, 'type'=>1, #'label'=>'D', 'max'=>comax, 'min'=>comin) GGraph.line( bouy_XYZmean.cut('t'=>cuttime), false, 'exchange'=>false , 'index'=>64, 'type'=>1, #'label'=>'bouy', 'max'=>comax, 'min'=>comin) #GGraph.line( all_cut.cut('t'=>cuttime), # false, 'exchange'=>false , # 'index'=>62, 'type'=>1, 'label'=>'all', # 'max'=>comax, 'min'=>comin) DCL.grcls