! !Description: ! 化学ポテンシャルをプロットして平衡状態を調べるためのプログラム ! ! Current Code Owner: sugiyama@gfd-dennou.org ! ! Histry: ! ! Version Date Comment ! ------- ---------- -------- ! 1.0 2004-07-26 新規作成 ! ! Code Description: ! Language: Fortran 90 ! Software Standaeds: "気象庁 標準コーディング" ! ! Modules used: ! ! Copyright (C) SUGIYAMA Ko-ichiro, 2000, All rights reserved program ObrChkMu use dcl use interface_netcdf_open use interface_netcdf_close use interface_netcdf_get use interface_display use interface_name_cnd2gas use interface_name2num implicit none real(8), allocatable :: mu_tmp(:,:) real(8), allocatable :: mu(:,:) real, allocatable :: mu_plot(:,:) real, allocatable :: mu_err_plot(:,:) real(8), allocatable :: pres(:) real, allocatable :: pres_plot(:) real(8), allocatable :: vap(:) real(8), allocatable :: svap(:) real(8), allocatable :: molfr(:,:) real(8), allocatable :: molfr_tmp(:,:) real, allocatable :: molfr_plot(:,:) real(8), allocatable :: mol(:,:) real(8), allocatable :: mol_tmp(:,:) real(8), allocatable :: temp(:) real(8), allocatable :: false8(:) real, allocatable :: false(:) real(8), parameter :: zero = 1.0d-40 integer :: x integer :: y integer :: l integer :: ncid integer :: k, k1, k2 integer :: s, s1, s2 integer :: index character(80) :: fname character(1) :: replay character(20) :: name character(20), allocatable :: sname(:) logical :: disp write(*,*) "============================================" write(*,*) "======NetCDF ファイルを指定して下さい=======" write(*,*) "============================================" read(*,*) fname ! fname = "obradbt-1000.nc" write(*,*) "============================================" write(*,*) "======図を 1 枚毎に出力しますか?============" write(*,*) "====== say Y or N ==========================" write(*,*) "============================================" read(*,*) replay if (replay == "Y" .OR. replay == "y") then disp = .false. elseif (replay == "N" .OR. replay == "n") then disp = .true. else write(*,*) "Unknown Answer: ", replay stop end if !!! !!!NetCDF ファイルのオープン !!!次元は pres, spc, elm であることを知っている !!! !ファイルのオープン call netcdf_open(trim(fname), ncid) !変数の取得. name = 'mu' call netcdf_get(ncid, name, mu_tmp) write(*,*) "mu" name = 'pres' call netcdf_get(ncid, name, pres) write(*,*) "pres" name = 'temp' call netcdf_get(ncid, name, temp) write(*,*) "temp" name = 'molfr' call netcdf_get(ncid, name, molfr_tmp) write(*,*) "molfr" name = 'mol' call netcdf_get(ncid, name, mol_tmp) write(*,*) "mol" name = 'sname' call netcdf_get(ncid, name, sname) write(*,*) "sname" !ファイルのクローズ call netcdf_close(ncid) !!! !!!プロットする配列を作る !!! y = size(mu_tmp, 1) x = size(mu_tmp, 2) allocate(mu(x,y), molfr_plot(x,y), & & mu_plot(x,y), mu_err_plot(x,y), pres_plot(x), & & vap(x), svap(x), molfr(x,y), mol(x,y), false(x), false8(x)) mu = transpose(mu_tmp) molfr = transpose(molfr_tmp) mol = transpose(mol_tmp) mu_plot = real(mu, 4) pres_plot = real(pres, 4) molfr_plot = real(molfr, 4) !!! !!!欠損値処理 !!! false = 999.0d0 false8 = 999.0d0 do s = 1, y l = len_trim(sname(s)) name = sname(s) mu_plot(:,s) = merge(mu_plot(:,s), false, mol(:,s) /= zero) end do !!! !!!プロット !!! !グラフのオープン call DclOpenGraphics() if (disp) call DclDivideFrame( 'y', 3, 2 ) if (.NOT. disp) call DclSetParm('WINDOW_HEIGHT', 600) if (.NOT. disp) call DclSetParm('WINDOW_WIDTH', 600) if (.NOT. disp) call DCLSetParm('DRAW_CORNERMARK', .false.) call DclSetFrameMargin( 0.05, 0.05, 0.05, 0.05 ) !!! !!!各化学ポテンシャルをプロット !!! do s = 1, y l = len_trim(sname(s)) name = sname(s) if (name(l:l) /= "g") cycle ! if (minval(mu_plot(:,s), mu_plot(:,s) > 0.0d0) & ! & == maxval(mu_plot(:,s), mu_plot(:,s) > 0.0d0)) cycle if (trim(sname(s)) == "NH4SH-s") cycle name = name(1:l-1) // 'l' call name2num(name, sname, k1) name = name(1:l-1) // 's' call name2num(name, sname, k2) call grfrm call uslset('lyinv', .true.) call gllset( 'LMISS', .TRUE. ) call usspnt(x, mu_plot(:,s), pres_plot) if (k1 /= 0) call usspnt(x, mu_plot(:,k1), pres_plot) if (k2 /= 0) call usspnt(x, mu_plot(:,k2), pres_plot) call grstrn(2) call uspfit call grstrf call ussttl('equilibrium constant', '1', 'Pressure', 'Pa') call usdaxs call sgspli(25) call sgplu(x, mu_plot(:,s), pres_plot) if (k2 /= 0) call sgspli(45) if (k2 /= 0) call sgplu(x, mu_plot(:,k2), pres_plot) if (k1 /= 0) call sgspli(35) if (k1 /= 0) call sgplu(x, mu_plot(:,k1), pres_plot) call uxsttl('T', trim(sname(s)), 0.0) !化学ポテンシャルの相対誤差 if (k1 == 0 .AND. k2 == 0) cycle call grfrm call uslset('lyinv', .true.) call gllset( 'LMISS', .TRUE. ) if (k1 /= 0) then mu(:,k1) = abs((mu(:,k1) - mu(:,s)) / mu(:,s)) mu_err_plot(:,k1) = real(merge(mu(:,k1), false8, mol(:,k1) /= zero),4) where (mu_err_plot(:,k1) <= 0.0d0) mu_err_plot(:,k1) = 1.0d-20 end where end if if (k2 /= 0) then mu(:,k2) = abs((mu(:,k2) - mu(:,s)) / mu(:,s)) mu_err_plot(:,k2) = real(merge(mu(:,k2), false8, mol(:,k2) /= zero),4) where (mu_err_plot(:,k2) <= 0.0d0) mu_err_plot(:,k2) = 1.0d-20 end where end if if (k1 /= 0) call usspnt(x, mu_err_plot(:,k1), pres_plot) if (k2 /= 0) call usspnt(x, mu_err_plot(:,k2), pres_plot) call grstrn(4) call uspfit call grstrf call ussttl('equilibrium constant', '1', 'Pressure', 'Pa') call usdaxs if (k1 /= 0) call sgspli(35) if (k1 /= 0) call sgplu(x, mu_err_plot(:,k1), pres_plot) if (k2 /= 0) call sgspli(45) if (k2 /= 0) call sgplu(x, mu_err_plot(:,k2), pres_plot) call uxsttl('T', trim(sname(s)), 0.0) end do GRAPH2: do s = 1, y if (trim(sname(s)) /= "NH4SH-s") cycle GRAPH2 name = 'NH3-g' call name2num(name, sname, s1) name = 'H2S-g' call name2num(name, sname, s2) vap = log(molfr(:,s1) * molfr(:,s2) * pres **2.0d0) - log(1.0d-2) svap = 61.781d0 - (10834.0d0 / temp) call grfrm call uslset('lyinv', .true.) call usspnt(x, real(vap, 4), pres_plot) call usspnt(x, real(svap, 4), pres_plot) call grstrn(2) call uspfit call grstrf call ussttl('equilibrium constant', '1', 'Pressure', 'Pa') call usdaxs call sgspli(25) call sgplu(x, real(vap, 4), pres_plot) call sgspli(35) call sgplu(x, real(svap, 4), pres_plot) call uxsttl('T', trim(sname(s)), 0.0) end do GRAPH2 call DclCloseGraphics() end program ObrChkMu