!スペクトルデータpsi, theta_s, chiから様々なグリッドデータに変換するプログラム ! !使い方 ! *ソース内の設定を適宜変更する。(パラメータ、切り出し開始点・終了点など) ! *コンパイルして実行ファイル(a.out)をつくる ! *第1引数にスペクトルデータが入れられたnetcdfファイル(data3D.nc)、第2・3引数に出力ファイル名を指定して実行 ! (例) a.out data3D.nc g1.nc g2.nc ! この場合g1.ncに半整数グリッド上のデータ、g2.ncに整数グリッド上のデータが書き出される ! !注意 ! *GEPフラックス = 一般化したEPフラックス ! !このコードには多くのムダやムラがあります。お気づきの点があれば適宜修正して下さい。 !Mac OSX 10.5.6 Macbook Core Duo上でg95を用いると正常に動作しない不具合を確認しています。 ! !2009.3.27 山本 博基 (hiroki@kugi.kyoto-u.ac.jp) ! !----------------------------------------------------------------- include "module_BPM3D_omp.f90" include "wrap_sj.f90" program uvwtheta implicit none call main !メインルーチン end !---------------------------------------------------------- subroutine main use gt4_history use constants use wrap_sj implicit none !****************************適宜変更する integer, parameter :: initdays=0 !切り出し始め位置 integer, parameter :: enddays=100 !切り出し終わり位置 integer, parameter :: ddays= 1 !データ抽出間隔 !**************************** real(8) :: psi((NN+1)*(NN+1),KM+2), theta_s((NN+1)*(NN+1),KM+2) real(8) :: chival((NN+1)*(NN+1),KM+2), chi((NN+1)*(NN+1),KM+1) real(8) :: z_half(KM),z_int(KM+1) real(8) :: theta_g(IM,JM,KM+2),theta_gval(IM,JM,KM) real(8) :: U(IM,JM,KM+2), V(IM,JM,KM+2) real(8) :: Uval(IM,JM,KM), Vval(IM,JM,KM), M(IM,JM,KM) real(8) :: w_g(IM,JM,KM+1), mpsi(IM,JM,KM+1), w_gval(IM,JM,KM) real(8) :: mpsival(IM,JM,KM) real(8) :: lon(IM), lat(JM) real(8) :: w_s((NN+1)*(NN+1),KM+2) real(8) :: chi1z((NN+1)*(NN+1),KM+2) real(8) :: theta_lambda(IM,JM,KM+2), theta_phi(IM,JM,KM+2) real(8) :: PV(IM, JM, KM) real(8) :: zeta(IM,JM,KM+2) !-------GEP計算用 real(8) :: NF_phi(JM,KM+1), NF_z(JM,KM+1), divNF(JM,KM), NF_phi_phi(JM,KM+1),divNF_mean(JM,KM) real(8) :: NFP(JM,KM+1), NF_phi_phi_half(JM,KM), theta_z(IM,JM,KM+1) real(8) :: U_mean(JM,KM+1), theta_mean(JM,KM+1), w_half(IM,JM,KM+1), Uval_mean(JM,KM+1) real(8) :: V_mean(JM,KM+1), w_mean(JM,KM+1), Uval_mean_phi(JM,KM+1), U_half_mean(JM,KM+2) real(8) :: theta_phi_mean(JM,KM+1), theta_z_mean(JM,KM+1), theta_phi_int(IM,JM,KM+1) real(8) :: vtheta_dash(IM,JM,KM+1), wtheta_dash(IM,JM,KM+1), uv_dash(IM,JM,KM+1), uw_dash(IM,JM,KM+1) real(8) :: vtheta_dash_mean(JM,KM+1), wtheta_dash_mean(JM,KM+1), uv_dash_mean(JM,KM+1), uw_dash_mean(JM,KM+1) real(8) :: mpsi_mean(JM,KM+1), resmpsi(JM,KM+1) real(8) :: U_int(IM,JM,KM+1),V_int(IM,JM,KM+1),theta_int(IM,JM,KM+1) real(8) :: M_mean(JM,KM+1), Mv_bar(JM,KM+1), Mw_bar(JM,KM+1), Mv_dag(JM,KM+1), Mw_dag(JM,KM+1) real(8) :: w_dag(JM,KM+1), v_dag(JM,KM+1) real(8) :: Mv_bar_phi(JM,KM+1), Mv_bar_phi_half(JM,KM) real(8) :: Mv_dag_phi(JM,KM+1), Mv_dag_phi_half(JM,KM) real(8) :: div_M_bar(JM,KM), div_M_dag(JM,KM) type(GT_HISTORY) :: hst_data, hst_EP character*80 filename,outname1,outname2, ctt character*5 :: time,cttt character*80 :: cuttime integer :: i,k, ct, j, days call SJ_init call philambda call getarg(1,filename) call getarg(2,outname1) call getarg(3,outname2) w_gval=0.0 time='time=' do ct=initdays, enddays, ddays write(ctt,*) ct cttt = adjustl(ctt) cuttime = time//cttt print *, cuttime call HistoryGet(file=filename, varname='psi,', array=psi, range=cuttime) call HistoryGet(file=filename, varname='theta_s', array=theta_s,range=cuttime) call HistoryGet(file=filename, varname='chi', array=chival, range=cuttime) chi = chival(1:(NN+1)*(NN+1),1:(KM+1)) do k=1, KM+1 call S_S2S_lap(chi(:,k),w_s(:,k)) ! chi から w_s を求める。 w_s(:,k) = -w_s(:,k) call S_S2G(w_s(:,k),w_g(:,:,k)) ! w_s から w_g を求める(内点のみ) chi1z(:,k)=(chi(:,k)-chi(:,k-1))*indz ! chi から chi1z を求める。 call S_S2G(theta_s(:,k), theta_g(:,:,k)) enddo do i=1 , IM do j=1 , JM theta_g(i,j,KM+2) = theta_g(i,j,KM+1) !境界条件 theta_g(i,j,1) = theta_g(i,j,2) w_g(i,j,KM+1) = 0.0 w_g(i,j,1) = 0.0 enddo enddo theta_gval(:,:,:)=theta_g(:,:,2:KM+1) call S2GUV(psi,chi1z,U,V) ! psi chi1z から U V を求める。 do i=1 , IM do j=1 , JM do k=2, KM+1 M(i,j,k-1) = U(i,j,k)*R + R*R*Omega*cphi(j)*cphi(j) Uval(i,j,k-1)=U(i,j,k)/cphi(j) Vval(i,j,k-1)=V(i,j,k)/cphi(j) enddo enddo enddo call makempsi(chi,mpsi) mpsival(:,:,:) = mpsi(:,:,1:KM) do k=1,KM+1 do j=1,JM mpsi_mean(j,k) = sum(mpsi(:,j,k))/real(IM) enddo enddo ! PVを計算---------- call Grad(theta_s,theta_lambda,theta_phi) call ZETAplus(psi,zeta) PV = 0.0 do k=2,KM+1 do j=1, JM do i=1, IM PV(i,j,k-1) = -theta_lambda(i,j,k)*(V(i,j,k+1)-V(i,j,k-1))/(cphi(j)*2.0*dz*R) & & +theta_phi(i,j,k)*(U(i,j,k+1)-U(i,j,k-1))/(cphi(j)*2.0*dz*R) & & +(zeta(i,j,k)+2.0*Omega*sphi(j))*(theta_g(i,j,k+1)-theta_g(i,j,k-1))/(2.0*dz) enddo enddo enddo ! GEPフラックスの計算--------------------------------------------------------------------- !Uint, Vint, thetaint生成 do k=1, KM+1 do j=1, JM do i=1, IM U_int(i,j,k) = (U(i,j,k+1)+U(i,j,k))/2.0 V_int(i,j,k) = (V(i,j,k+1)+V(i,j,k))/2.0 theta_int(i,j,k) = (theta_g(i,j,k+1)+theta_g(i,j,k))/2.0 theta_phi_int(i,j,k) = (theta_phi(i,j,k+1)+theta_phi(i,j,k))/2.0 enddo enddo enddo !meanの生成 do k=1, KM+1 do j=1, JM V_mean(j,k) = sum(V_int(:,j,k))/real(IM) U_mean(j,k) = sum(U_int(:,j,k))/real(IM) Uval_mean(j,k) = U_mean(j,k)/cphi(j) theta_mean(j,k) = sum(theta_int(:,j,k))/real(IM) w_mean(j,k) = sum(w_g(:,j,k))/real(IM) enddo enddo do k=1, KM+2 do j=1, JM U_half_mean(j,k) = sum(U(:,j,k))/real(IM) enddo enddo ! dashの生成 do k=1,KM+1 do j=1,JM do i=1, IM vtheta_dash(i,j,k) = (theta_int(i,j,k)-theta_mean(j,k))*(V_int(i,j,k)-V_mean(j,k))/cphi(j) wtheta_dash(i,j,k) = (theta_int(i,j,k)-theta_mean(j,k))*(w_g(i,j,k)-w_mean(j,k)) uv_dash(i,j,k) = (U_int(i,j,k)-U_mean(j,k))*(V_int(i,j,k)-V_mean(j,k))/(cphi(j)*cphi(j)) uw_dash(i,j,k) = (U_int(i,j,k)-U_mean(j,k))*(w_g(i,j,k)-w_mean(j,k))/cphi(j) enddo enddo enddo ! dashdasubarの生成 do k=1, KM+1 do j=1, JM vtheta_dash_mean(j,k) = sum(vtheta_dash(:,j,k))/real(IM) wtheta_dash_mean(j,k) = sum(wtheta_dash(:,j,k))/real(IM) uv_dash_mean(j,k) = sum(uv_dash(:,j,k))/real(IM) uw_dash_mean(j,k) = sum(uw_dash(:,j,k))/real(IM) enddo enddo !theta_phi_mean theta_z_meanの生成 do k=1, KM+1 do j=1, JM theta_phi_mean(j,k) = (sum(theta_phi_int(:,j,k)))/real(IM) enddo enddo do k=1, KM+1 do j=1, JM do i=1, IM theta_z(i,j,k) = (theta_g(i,j,k+1)-theta_g(i,j,k))/dz enddo theta_z_mean(j,k) = sum(theta_z(:,j,k))/real(IM) enddo enddo !NFP生成 do k=1, KM+1 do j=1, JM !GEPフラックスのとき用 NFP(j,k) = (vtheta_dash_mean(j,k)*theta_z_mean(j,k)-wtheta_dash_mean(j,k)*theta_phi_mean(j,k)/R)& &/(theta_phi_mean(j,k)*theta_phi_mean(j,k)/(R*R)+theta_z_mean(j,k)*theta_z_mean(j,k)) !伝統的なEPフラックス用 !NFP(j,k) = vtheta_dash_mean(j,k)/theta_z_mean(j,k) enddo enddo !print *, maxval(NFP), minval(NFP) !NF_phi, NF_z生成 call divphi(Uval_mean, Uval_mean_phi) do j=1, JM do k=1,KM+1 NF_phi(j,k) = R*cphi(j)*(((U_half_mean(j,k+1)-U_half_mean(j,k))/& &(dz*cphi(j)))*NFP(j,k)-uv_dash_mean(j,k)) NF_z(j,k) = R*cphi(j)*((2.0*Omega*sphi(j)-Uval_mean_phi(j,k)/R)& &*NFP(j,k)-uw_dash_mean(j,k)) enddo !NF_phi(j,1) = -R*cphi(j)*uv_dash_mean(j,1) !NF_z(j,1) = 0.0 !NF_phi(j,KM+1) = -R*cphi(j)*uv_dash_mean(j,KM+1) !NF_z(j,KM+1) = 0.0 enddo ! div F の生成 call divphi(NF_phi,NF_phi_phi) do k=1,KM do j=1,JM NF_phi_phi_half(j,k) = (NF_phi_phi(j,k+1)+NF_phi_phi(j,k))/2.0 enddo enddo do k=1,KM do j=1, JM divNF(j,k) = NF_phi_phi_half(j,k)/R+(NF_z(j,k+1)-NF_z(j,k))/(dz) enddo enddo !残差循環 resmpsi = mpsi_mean - NFP(:,1:KM+1) !---------------------------------------------- !M_mean作成 do k=1, KM+1 do j=1, JM M_mean(j,k) = U_mean(j,k)*R + Omega*R*R*cphi(j)*cphi(j) enddo enddo !w_dag作成 call divphi(NFP,w_dag) w_dag = w_dag/R !v_dag作成 !内点 do k=2, KM do j=1, JM v_dag(j,k) = -(NFP(j,k+1)-NFP(j,k-1))/(2.0*dz) enddo enddo !上端 do j=1, JM v_dag(j,KM+1) = -(NFP(j,KM+1)-NFP(j,KM))/(dz) enddo !下端 do j=1,JM v_dag(j,1) = -0.5*(1+alpha)*(NFP(j,2)-NFP(j,1))/dz enddo !Mv_bar, Mw_bar, Mv_dag, Mw_dag作成 do k=1, KM+1 do j=1, JM Mv_bar(j,k) = M_mean(j,k)*V_mean(j,k)/cphi(j) enddo enddo Mw_bar = M_mean*w_mean Mv_dag = M_mean*v_dag Mw_dag = M_mean*w_dag !div_M_bar作成 call divphi(Mv_bar,Mv_bar_phi) do k=1,KM do j=1,JM Mv_bar_phi_half(j,k) = (Mv_bar_phi(j,k+1)+Mv_bar_phi(j,k))/2.0 enddo enddo do k=1,KM do j=1, JM div_M_bar(j,k) = Mv_bar_phi_half(j,k)/R+(Mw_bar(j,k+1)-Mw_bar(j,k))/(dz) enddo enddo !div_M_dag作成 call divphi(Mv_dag,Mv_dag_phi) do k=1,KM do j=1,JM Mv_dag_phi_half(j,k) = (Mv_dag_phi(j,k+1)+Mv_dag_phi(j,k))/2.0 enddo enddo do k=1,KM do j=1, JM div_M_dag(j,k) = Mv_dag_phi_half(j,k)/R+(Mw_dag(j,k+1)-Mw_dag(j,k))/(dz) enddo enddo if (ct==initdays) then call HistoryCreate(file=outname1,title='Data',source='DataConversion.f90',& &institution='user',& & dims=(/'lon ','lat ','z ','time'/),dimsizes=(/IM,JM,KM,0/),& &longnames=(/'lon ','lat ','z ','time'/), & & units=(/'deg ','deg ','m ','days'/), origin=real(initdays),interval=1.0,history=hst_data) do i=1,IM lon(i) = lambda(i)*180.0/pi enddo do j=1, JM lat(j) = phi(j)*180.0/pi enddo do k=1,KM z_half(k) = real(k-1)*dz+dz/2.0 enddo call HistoryPut('lon',lon,hst_data) call HistoryPut('lat',lat,hst_data) call HistoryPut('z',z_half,hst_data) call HistoryAddVariable(varname='u',dims=(/'lon ','lat ','z ','time'/),& &longname='u',units='m/s',xtype='float',history=hst_data) call HistoryAddVariable(varname='v',dims=(/'lon ','lat ','z ','time'/),& &longname='v',units='m/s',xtype='float',history=hst_data) call HistoryAddVariable(varname='theta',dims=(/'lon ','lat ','z ','time'/),& &longname='theta',units='K',xtype='float',history=hst_data) call HistoryAddVariable(varname='M',dims=(/'lon ','lat ','z ','time'/),& &longname='angular_momentum',units='m/s',xtype='float',history=hst_data) call HistoryAddVariable(varname='PV',dims=(/'lon ','lat ','z ','time'/),& &longname='potential vorticity',units='K/ms',xtype='float',history=hst_data) call HistoryAddVariable(varname='epflx_div',dims=(/'lat ','z ','time'/),& &longname='GEPflux_div',units='m s-2',xtype='float',history=hst_data) call HistoryAddVariable(varname='M_bar_div',dims=(/'lat ','z ','time'/),& &longname='M_bar_div',units='m2 s-2',xtype='float',history=hst_data) call HistoryAddVariable(varname='M_dag_div',dims=(/'lat ','z ','time'/),& &longname='M_dag_div',units='m2 s-2',xtype='float',history=hst_data) call HistoryCreate(file=outname2,title='Data',source='DataConversion.f90',& &institution='user',& & dims=(/'lon ','lat ','z ','time'/),dimsizes=(/IM,JM,KM+1,0/),& &longnames=(/'lon ','lat ','z ','time'/), & & units=(/'deg ','deg ','m ','days'/), origin=real(initdays),interval=1.0,history=hst_EP) do k=1,KM+1 z_int(k) = real(k-1)*dz enddo call HistoryPut('lon',lon,hst_EP) call HistoryPut('lat',lat,hst_EP) call HistoryPut('z',z_int,hst_EP) call HistoryAddVariable(varname='epflx_y',dims=(/'lat ','z ','time'/),& &longname='GEPflux_y',units='m2 s-2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='epflx_z',dims=(/'lat ','z ','time'/),& &longname='GEPflux_z',units='m2 s-2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='mpsi',dims=(/'lon ','lat ','z ','time'/),& &longname='meridional_psi',units='m2/s',xtype='float',history=hst_EP) call HistoryAddVariable(varname='rmpsi',dims=(/'lat ','z ','time'/),& &longname='r_meridional_psi',units='m2/s',xtype='float',history=hst_EP) call HistoryAddVariable(varname='w',dims=(/'lon ','lat ','z ','time'/),& &longname='w',units='m/s',xtype='float',history=hst_EP) call HistoryAddVariable(varname='uv_dash_mean',dims=(/'lat ','z ','time'/),& &longname='uv_dash_mean',units='m2/s2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='uw_dash_mean',dims=(/'lat ','z ','time'/),& &longname='uw_dash_mean',units='m2/s2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='vtheta_dash_mean',dims=(/'lat ','z ','time'/),& &longname='vtheta_dash_mean',units='K m/s',xtype='float',history=hst_EP) call HistoryAddVariable(varname='wtheta_dash_mean',dims=(/'lat ','z ','time'/),& &longname='wtheta_dash_mean',units='K m/s',xtype='float',history=hst_EP) call HistoryAddVariable(varname='Mv_bar',dims=(/'lat ','z ','time'/),& &longname='Mv_bar',units='m3/s2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='Mw_bar',dims=(/'lat ','z ','time'/),& &longname='Mw_bar',units='m3/s2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='Mv_dag',dims=(/'lat ','z ','time'/),& &longname='Mv_dag',units='m3/s2',xtype='float',history=hst_EP) call HistoryAddVariable(varname='Mw_dag',dims=(/'lat ','z ','time'/),& &longname='Mw_dag',units='m3/s2',xtype='float',history=hst_EP) endif call HistoryPut('u', Uval,hst_data) call HistoryPut('v', Vval,hst_data) call HistoryPut('w', w_g,hst_EP) call HistoryPut('theta', theta_gval,hst_data) call HistoryPut('M',M,hst_data) call HistoryPut('mpsi',mpsi,hst_EP) call HistoryPut('rmpsi',resmpsi,hst_EP) call HistoryPut('PV',PV,hst_data) call HistoryPut('epflx_y',NF_phi,hst_EP) call HistoryPut('epflx_z',NF_z,hst_EP) call HistoryPut('epflx_div',divNF,hst_data) call HistoryPut('uv_dash_mean',uv_dash_mean,hst_EP) call HistoryPut('uw_dash_mean',uw_dash_mean,hst_EP) call HistoryPut('vtheta_dash_mean',vtheta_dash_mean,hst_EP) call HistoryPut('wtheta_dash_mean',wtheta_dash_mean,hst_EP) call HistoryPut('M_bar_div',div_M_bar,hst_data) call HistoryPut('M_dag_div',div_M_dag,hst_data) call HistoryPut('Mv_bar',Mv_bar,hst_EP) call HistoryPut('Mw_bar',Mw_bar,hst_EP) call HistoryPut('Mv_dag',Mv_dag,hst_EP) call HistoryPut('Mw_dag',Mw_dag,hst_EP) enddo call HistoryClose(hst_data) call HistoryClose(hst_EP) endsubroutine !--------------------------------------------------------------------------- subroutine defineTE(theta_es) !ニュートン加熱冷却の基準温位場生成 use constants use wrap_sj implicit none real(8) :: theta_es((NN+1)*(NN+1),KM+2) real(8) :: TE(IM,JM,KM+2) integer :: i,j,k do k=1, KM+2 do i=1,IM do j=1,JM TE(i,j,k)=theta_0*(-0.33333*(sphi(j)*sphi(j)-1.0/3.0)+0.125*((real(k-1,8)-0.5)*dz/H-0.5)) !TE(i,j,k)=theta_0*0.0 enddo enddo call S_G2S(TE(:,:,k), theta_es(:,k)) enddo endsubroutine !------------------------------------------------------------------ subroutine zbibun1(x,x1) ! スペクトルデータの整数グリッドの値 x (k=0 ~ KM)から半整数グリッドのdx/dz (k=1/2 ~ KM-1/2)をもとめる。 use constants use wrap_sj implicit none real(8) :: x((NN+1)*(NN+1),KM+1) real(8) :: x1((NN+1)*(NN+1),KM+2) integer :: n,k do n=1, (NN+1)*(NN+1) do k=1, KM x1(n,k+1)=(x(n,k+1)-x(n,k))/dz enddo enddo endsubroutine !------------------------------------------------------------------- subroutine S2GUV(x,y,a,b) !スペクトルデータからU,Vのグリッドデータを計算 use constants use wrap_sj implicit none real(8) :: x((NN+1)*(NN+1),KM+2) real(8) :: y((NN+1)*(NN+1),KM+2) real(8) :: a(IM,JM,KM+2) real(8) :: b(IM,JM,KM+2) real(8) :: sx1b((NN+1)*(NN+1), KM+2) real(8) :: sx2b((NN+1)*(NN+1), KM+2) real(8) :: g1a(IM,JM,KM+2) , g2a(IM,JM,KM+2) real(8) :: g1b(IM,JM,KM+2) , g2b(IM,JM,KM+2) integer :: i , k ,j ,n !real(8) :: Q(JM/2*11), WS1(2*(NN+2)),WS2(2*(NN+2)),W1((JM+1)*IM),W2((JM+1)*IM),WG((IM+2)*JM) !$omp parallel do do k = 2, KM+1 call S_S2G_dlat_dual(x(:,k),y(:,k),g1a(:,:,k),g2a(:,:,k)) call S_S2S_dlon(x(:,k),sx1b(:,k)) call S_S2S_dlon(y(:,k),sx2b(:,k)) call S_S2G_dual(sx1b(:,k),sx2b(:,k),g1b(:,:,k),g2b(:,:,k)) !call SJCS2Y(NN,s1a,sy1a,SJC) !call SJCS2Y(NN,s2a,sy2a,SJC) !call SJMS2G(NN,NM,NN+1,IM,JM,sy1a,sy2a,g1a,g2a,SJIT,SJT,SJP,Q,SJR,WS1,WS2,WG,W1,W2,1) !call SJCS2X(NN,s1b,sx1b) !call SJCS2X(NN,s2b,sx2b) !call SJMS2G(NN,NM,NN,IM,JM,sx1b,sx2b,g1b,g2b,SJIT,SJT,SJP,Q,SJR,WS1,WS2,WG,W1,W2,0) do i=1 , IM do j=1 , JM a(i,j,k) = -g1a(i,j,k)*R*cphi(j) +g2b(i,j,k)*R ! Uの生成 b(i,j,k) = g1b(i,j,k)*R +g2a(i,j,k)*R*cphi(j) ! Vの生成 enddo enddo enddo !$omp end parallel do !$omp parallel do do i=1,IM do j=1,JM a(i,j,1)=alpha*a(i,j,2) !境界条件 a(i,j,KM+2)=a(i,j,KM+1) b(i,j,1) = alpha*b(i,j,2) !境界条件 b(i,j,KM+2) = b(i,j,KM+1) enddo enddo !$omp end parallel do endsubroutine !------------------------------------------------------------------------------------- subroutine ZETAplus(psi,zeta) ! lap psi (- 2/a^2 psi) を求める use constants use wrap_sj implicit none real(8) :: psi((NN+1)*(NN+1),KM+2) real(8) :: lap_psi((NN+1)*(NN+1),KM+2), zeta(IM,JM,KM+2) integer :: i,k do k=1,KM+2 do i=1,(NN+1)*(NN+1) lap_psi(i,k)=-psi(i,k)*(SJn(i))*(SJn(i)+1.0) enddo call S_S2G(lap_psi(:,k),zeta(:,:,k)) enddo endsubroutine !---------------------------------------------------------------------------------- subroutine lap2(x,y) ! y = -nabla^2 x つまり y = n(n+1)xをもとめる。(wを求める) use constants use wrap_sj implicit none real(8) :: x((NN+1)*(NN+1),KM+1) real(8) :: y((NN+1)*(NN+1),KM+1) integer :: i,k do i=1,(NN+1)*(NN+1) do k=1,KM+1 y(i,k) = x(i,k)*(SJn(i))*(SJn(i)+1.0) !print *, n , m enddo enddo endsubroutine !----------------------------------------------------------------------------------- subroutine Grad(s,g1,g2) ! スペクトル → グリッド 勾配 use constants use wrap_sj implicit none real(8) :: s((NN+1)*(NN+1),KM+2), slon((NN+1)*(NN+1),KM+2) real(8) :: g1(IM,JM,KM+2), g2(IM,JM,KM+2) integer :: i,j,k !IPOW = 1 !IFLAG = -1 !call SJTS2G(NN,IM,IM,JM,JM,KM+2,s,g1,snIT,snT,snY,snIPK,snPK,snRK,snIA,snA,snQ,snWS,snWW,IPOW,IFLAG) !IFLAG = 1 !call SJTS2G(NN,IM,IM,JM,JM,KM+2,s,g2,snIT,snT,snY,snIPK,snPK,snRK,snIA,snA,snQ,snWS,snWW,IPOW,IFLAG) do k=1,KM+2 call S_S2G_dlat(s(:,k), g2(:,:,k)) call S_S2S_dlon(s(:,k), slon(:,k)) call S_S2G(slon(:,k),g1(:,:,k)) do j=1, JM do i = 1, IM g1(i,j,k) = g1(i,j,k)*incphi(j) g2(i,j,k) = g2(i,j,k) enddo enddo enddo endsubroutine !----------------------------------------------------------------------------- subroutine divphi(g1,g2) ! 2次元グリッド → 2次元グリッド phi発散 use constants use wrap_sj implicit none real(8) :: s((NN+1)*(NN+1),KM+1) real(8) :: g1(JM,KM+1), g2(JM,KM+1) real(8) :: g31(IM,JM,KM+1), g32(IM,JM,KM+1) integer :: i, k do i=1, IM g31(i,:,:) = g1(:,:) enddo !print *, maxval(g31), minval(g31) do k =1 ,KM+1 !IPOW = 1 !IFLAG = 1 !call SJTG2S(NN,IM,IM,JM,JM,KM+1,g31,s,snIT,snT,snY,isnIPK,isnPK,isnRK,snIA,snA,isnQ,isnWS,isnWW,IPOW,IFLAG) call S_G2S_div(g31(:,:,k), s(:,k)) call S_S2G(s(:,k),g32(:,:,k)) enddo g2(:,:) = g32(1,:,:) !print *, maxval(g32), minval(g32) endsubroutine !----------------------------------------------------------------------------- subroutine zbibun2(x,x1) ! グリッドデータの半整数グリッドの値 x から整数グリッドのdx/dz をもとめる。 use constants use wrap_sj implicit none real(8) :: x(IM,JM,KM+2) real(8) :: x1(IM,JM,KM+1) integer :: i,j,k do i=1, IM do j=1, JM do k=2, KM x1(i,j,k)=(x(i,j,k+1)-x(i,j,k))/dz enddo enddo enddo endsubroutine !------------------------------------------------------------------- !-------------------------------------------------------------------- subroutine makempsi(chi,mpsi) use constants use wrap_sj implicit none real(8) :: chi((NN+1)*(NN+1),KM+1) real(8) :: mpsi(IM,JM,KM+1) integer :: i,j,k real(8) :: s((NN+1)*(NN+1),KM+1),g(IM,JM,KM+1) !IFLAG = 1 !IPOW = 1 !call SJTS2G(NN,IM,IM,JM,JM,KM+1,chi,g,snIT,snT,snY,isnIPK,isnPK,isnRK,snIA,snA,isnQ,isnWS,isnWW,IPOW,IFLAG) do k =1, KM+1 call S_S2G_dlat(chi(:,k), g(:,:,k)) do i=1,IM do j=1,JM mpsi(i,j,k) = g(i,j,k)*R enddo enddo enddo endsubroutine !---------------------------------------------------------------------