!----------------------------------------------------------------- include "module_BPM2D.f90" include "wrap_lj.f90" program busi12 use constants implicit none call main !メインルーチン呼び出し contains !---------------------------------------------------------- subroutine main !メインルーチン use gt4_history use constants use wrap_lj implicit none ! スペクトルデータ real(8) :: psi((NN+1),KM+2) , oldpsi((NN+1),KM+2) real(8) :: kpsi((NN+1),KM+2,4) !ルンゲクッタ用 real(8) :: theta_s((NN+1),KM+2), oldtheta_s((NN+1),KM+2) real(8) :: ktheta_s((NN+1),KM+2,4) !ルンゲクッタ用 real(8) :: chi2z((NN+1),KM+1) , oldchi2z((NN+1),KM+1) real(8) :: kchi2z((NN+1),KM+1,4) !ルンゲクッタ用 real(8) :: chi((NN+1),KM+1) real(8) :: zeta_s((NN+1),KM+2) real(8) :: E_s((NN+1),KM+2) real(8) :: wtheta_s((NN+1),KM+2),Vtheta_s((NN+1),KM+2),Utheta_s((NN+1),KM+2) real(8) :: chi1z((NN+1),KM+2) real(8) :: w_s((NN+1),KM+1) real(8) :: theta_es((NN+1),KM+2) ! グリッドデータ real(8) :: psigrid(JM,KM+2) real(8) :: chigrid(JM,KM+1) real(8) :: theta_g(JM,KM+2) real(8) :: U(JM,KM+2), V(JM,KM+2) real(8) :: zeta_g(JM,KM+2) real(8) :: w_g(JM,KM+1) real(8) :: A_g(JM,KM+2) real(8) :: A1z(JM,KM+1) real(8) :: B_g(JM,KM+2) real(8) :: B1z(JM,KM+1) real(8) :: E_g(JM,KM+2) real(8) :: Utheta_g(JM,KM+2) real(8) :: Utheta1z(JM,KM+1) real(8) :: Vtheta_g(JM,KM+2) real(8) :: Vtheta1z(JM,KM+1) real(8) :: wtheta_g(JM,KM+2) real(8) :: rdm real(8) :: totalMold, totalMnew real(8) :: spct(NN+1) real(8) :: chi2zval((NN+1),KM+2) integer :: i,j,k,n,tstep, time ! tstep = ルンゲクッタのカウンター type(GT_HISTORY) :: hst_data1, hst_data2, hst_restart character(len=80) :: filename, filename2 real(8) :: apsi(MM+1),achi(MM+1),atheta(MM+1) real(8) :: ex_apsi(MM+1,KM+2),ex_achi(MM+1,KM+2),ex_atheta(MM+1,KM+2) totalMold = 0.0 call LJ_init call philambda call defineTE(theta_es) do i=1,(NN+1) theta_es(i,KM+2) = theta_es(i,KM+1) theta_es(i,1) = theta_es(i,2) enddo print *, "Theta_e was defined" call RANDOM_SEED if(initial_condition == 1)then do k=2, KM+1 do j=1, JM psigrid(j,k) = - Omega*((1+((real(k-1,8)-0.5)*dz)*(2.0*ga*delta_H)/& &(R*R*Omega*Omega))**0.5-1)*sphi(j) enddo enddo psigrid(:,KM+2) = psigrid(:,KM+1) psigrid(:,1) = psigrid(:,2) do k =1, KM+2 call L_G2S(psigrid(:,k), psi(:,k)) enddo do k=1, KM+2 call L_S2G(theta_es(:,k), theta_g(:,k)) !do j =1, JM !call RANDOM_NUMBER(rdm) !theta_g(j,k) = theta_g(j,k)*(1+0.01*(rdm-0.5)) !enddo call L_G2S(theta_g(:,k), theta_s(:,k)) enddo print *, "Initial condition was set with thermal wind in theta_e " chi2z = 0.0 else !omp parallel do private(rdm) do k=1,KM+2 do j=1,JM call RANDOM_NUMBER(rdm) !theta_g(j,k) = rdm*0.1 theta_g(j,k) = 0.0 enddo call L_G2S(theta_g(:,k),theta_s(:,k)) enddo !omp end prallel do print *, "Initial condition was set with constant PT in rest." endif !print *, "Initial condition was set." !-------------------リスタートファイルの読み込み:この部分を有効にするとリスタートファイルを読む if(restart == 1) then filename = 'restart2D.nc' call HistoryGet(file=filename, varname='psi', array=psi) call HistoryGet(file=filename, varname='theta_s', array=theta_s) call HistoryGet(file=filename, varname='chi2z', array=chi2zval) do k=1,KM+1 do n=1,MM+1 chi2z(n,k)=chi2zval(n,k) enddo enddo endif !------------------------------------------ !!!!!!!!!!!!!!!!!!!!!!!!!水平の力学チェック P^1_1(mu)=root3 cos(phi) Y^1_1=root3 cos(phi)cos(lambda) !do i=1,IM !do j=1,JM !do k=1, KM+2 !psigrid(j,k) = sqrt(3.0)*cos(phi(j))*cos(lambda(i))/R !enddo !enddo !enddo !do k=1,KM+2 !call L_G2S(psigrid(:,k),psi(:,k)) !enddo !*********************** do n=1,MM+1 !硬い方程式の 解核行列用 apsi(n) = -nu_H*(LJn(n)*(LJn(n)+1.0)-2.0)*inR*inR achi(n) = -nu_H*(2.0*LJn(n)*(LJn(n)+1.0)-2.0)*inR*inR atheta(n) = -kappa_H*(LJn(n)*(LJn(n)+1.0))*inR*inR do k = 1, KM+2 ex_apsi(n,k) = dexp(apsi(n)*dt*0.5) ex_atheta(n,k) = dexp(atheta(n)*dt*0.5) ex_achi(n,k) = dexp(achi(n)*dt*0.5) enddo enddo !********************** do time=1,maxdaystep !------------time step !-------- oldpsi = psi oldtheta_s = theta_s oldchi2z = chi2z do tstep = 0, 3 if (tstep == 1) then do i=1,(MM+1) do k=2,KM+1 psi(i,k) = ex_apsi(i,k)*(oldpsi(i,k) + kpsi(i,k,1)*0.5) theta_s(i,k) = ex_atheta(i,k)*(oldtheta_s(i,k) + ktheta_s(i,k,1)*0.5) chi2z(i,k) = ex_achi(i,k)*(oldchi2z(i,k) + kchi2z(i,k,1)*0.5) enddo enddo endif if (tstep == 2) then do i=1,(MM+1) do k=2,KM+1 psi(i,k) = ex_apsi(i,k)*oldpsi(i,k) + kpsi(i,k,2)*0.5 theta_s(i,k) = ex_atheta(i,k)*oldtheta_s(i,k) + ktheta_s(i,k,2)*0.5 chi2z(i,k) = ex_achi(i,k)*oldchi2z(i,k) + kchi2z(i,k,2)*0.5 enddo enddo endif if (tstep == 3) then do i=1,(MM+1) do k=2,KM+1 psi(i,k) = ex_apsi(i,k)*(ex_apsi(i,k)*oldpsi(i,k) + kpsi(i,k,3)) theta_s(i,k) = ex_atheta(i,k)*& &(ex_atheta(i,k)*oldtheta_s(i,k) + ktheta_s(i,k,3)) chi2z(i,k) = ex_achi(i,k)*(ex_achi(i,k)*oldchi2z(i,k) + kchi2z(i,k,3)) enddo enddo endif call useLAPACK(chi2z,chi) ! chi2z から chi を求める。 !$omp parallel do do k=2, KM+1 call L_S2S_lap(psi(:,k),zeta_s(:,k)) ! psi から zeta_s を求める。 call L_S2G_dual(theta_s(:,k),zeta_s(:,k),theta_g(:,k),zeta_g(:,k)) ! theta_s から theta_g , zeta_s から zeta_g を求める(内点のみ) call L_S2S_lap(chi(:,k),w_s(:,k)) ! chi から w_s を求める。 w_s(:,k) = -w_s(:,k) call L_S2G(w_s(:,k),w_g(:,k)) ! w_s から w_g を求める(内点のみ) chi1z(:,k)=(chi(:,k)-chi(:,k-1))*indz ! chi から chi1z を求める。 enddo !$omp end parallel do call S2GUV(psi,chi1z,U,V) ! psi chi1z から U V を求める。 !$omp parallel do do i = 1, (NN+1) !境界条件 zeta_s(i,KM+2) = zeta_s(i,KM+1) psi(i,KM+2) = psi(i,KM+1) zeta_s(i,1) = alpha*zeta_s(i,2) psi(i,1) = alpha*psi(i,2) theta_s(i,KM+2) = theta_s(i,KM+1) theta_s(i,1) = theta_s(i,2) enddo !$omp end parallel do !$omp parallel do do j=1 , JM theta_g(j,KM+2) = theta_g(j,KM+1) !境界条件 theta_g(j,1) = theta_g(j,2) zeta_g(j,KM+2) = zeta_g(j,KM+1) zeta_g(j,1) = alpha*zeta_g(j,2) w_g(j,KM+1) = 0.0 w_g(j,1) = 0.0 enddo !$omp end parallel do ! 非線形項の計算。(グリッド)-------------------------------------- !$omp parallel do do k=2, KM+1 do j=1,JM A_g(j,k) = (zeta_g(j,k)+2.0*Omega*sphi(j))*U(j,k) + (w_g(j,k)+w_g(j,k-1))*0.5*(V(j,k+1)-V(j,k-1))*0.5*indz B_g(j,k) = (zeta_g(j,k)+2.0*Omega*sphi(j))*V(j,k) - (w_g(j,k)+w_g(j,k-1))*0.5*(U(j,k+1)-U(j,k-1))*0.5*indz enddo enddo !$omp end parallel do !$omp parallel do do k=1, KM+2 do j=1,JM E_g(j,k) = (U(j,k)*U(j,k)+V(j,k)*V(j,k))*0.5*incphi(j)*incphi(j) Utheta_g(j,k) = U(j,k)*theta_g(j,k) Vtheta_g(j,k) = V(j,k)*theta_g(j,k) enddo enddo !$omp end parallel do !$omp parallel do do k=1, KM+1 do j=1,JM wtheta_g(j,k) = w_g(j,k)*(theta_g(j,k+1)+theta_g(j,k))*0.5 enddo enddo !$omp end parallel do ! 非線形項のスペクトル変換 !$omp parallel do do k =1, KM+1 call L_G2S_dual(E_g(:,k),wtheta_g(:,k),E_s(:,k),wtheta_s(:,k)) enddo !$omp end parallel do ! グリッドデータの半整数グリッドの値 x から整数グリッドのdx/dz をもとめる。 !$omp parallel do do k=2, KM do j=1, JM A1z(j,k)=(A_g(j,k+1)-A_g(j,k))*indz B1z(j,k)=(B_g(j,k+1)-B_g(j,k))*indz Utheta1z(j,k)=(Utheta_g(j,k+1)-Utheta_g(j,k))*indz Vtheta1z(j,k)=(Vtheta_g(j,k+1)-Vtheta_g(j,k))*indz enddo enddo !$omp end parallel do !$omp parallel do do i=1, (NN+1) !境界条件 chi2z(i,KM+1) = 0.0 chi2z(i,1) = chi(i,2)*(1.0-alpha)*indz*indz enddo !$omp end parallel do if (tstep < 4)then call nextpsi(psi,B_g,kpsi,tstep) call nextchi2z(chi2z,A1z,theta_s,E_s,kchi2z,tstep) call nexttheta(theta_s,Vtheta_g,wtheta_s,theta_es,ktheta_s,tstep) endif !---------------ouput-------------- if(mod(time,daystep*INT(output_interval))==0)then if(tstep==0)then print *, time/daystep+INT(init_day), 'days' !print *, "u,v", maxval(U), maxval(V) call output(time,U,V,theta_g,w_g,chi,theta_s,psi,chi2z,hst_data1,hst_data2,hst_restart,TotalMnew) !open(8,FILE='psiwave.txt') !数値粘性が不十分でないかの確認用に出力する !k=KM+1 !do i=1,(NN+1) !spct(i) = spct(i) + psi(i,k)*psi(i,k) !enddo !do i=1, (NN) !write(8,*) i, spct(i)!*(i*(i+1)) !enddo !close(8) if(abs(totalMold - totalMnew) < 0.0001)then !定常状態の判定。 call HistoryClose(hst_data1) call HistoryClose(hst_data2) call HistoryClose(hst_restart) stop elseif(totalMnew == totalMold)then call HistoryClose(hst_data1) call HistoryClose(hst_data2) call HistoryClose(hst_restart) stop else totalMold = totalMnew endif endif endif enddo !$omp parallel do do i=1,(NN+1) do k=2,KM+1 psi(i,k) = ex_apsi(i,k)*(ex_apsi(i,k)*(oldpsi(i,k) + kpsi(i,k,1)/6.0) & &+kpsi(i,k,2)/3.0 + kpsi(i,k,3)/3.0) + kpsi(i,k,4)/6.0 theta_s(i,k) = ex_atheta(i,k)*(ex_atheta(i,k)*(oldtheta_s(i,k) & &+ ktheta_s(i,k,1)/6.0) + ktheta_s(i,k,2)/3.0 & &+ ktheta_s(i,k,3)/3.0) + ktheta_s(i,k,4)/6.0 chi2z(i,k) = ex_achi(i,k)*(ex_achi(i,k)*(oldchi2z(i,k) & &+ kchi2z(i,k,1)/6.0) + kchi2z(i,k,2)/3.0 & &+ kchi2z(i,k,3)/3.0) + kchi2z(i,k,4)/6.0 enddo enddo !$omp end parallel do enddo call HistoryClose(hst_data1) call HistoryClose(hst_data2) call HistoryClose(hst_restart) endsubroutine !--------------------------------------------------------------------------------- subroutine defineTE(theta_es) !ニュートン加熱冷却の基準温位場生成 use constants use wrap_lj implicit none real(8) :: theta_es((NN+1),KM+2) real(8) :: TE(JM,KM+2) integer :: j,k do k=1, KM+2 do j=1,JM TE(j,k)=theta_0*(-delta_H*(sphi(j)*sphi(j)-1.0/3.0)+0.125*((real(k-1,8)-0.5)*dz/H-0.5)) !TE(j,k)=theta_0*0.0 enddo call L_G2S(TE(:,k), theta_es(:,k)) enddo endsubroutine !--------------------------------------------------------------------------------- subroutine useLAPACK(x2,x) !三重対角行列で表される連立1次方程式を解く use constants use wrap_lj implicit none real(8) :: x2((NN+1),KM+1) real(8) :: x((NN+1),KM+1) real(8) :: D(KM-1) real(8) :: DU(KM-2),DL(KM-2),DU2(KM-3) real(8) :: B(KM-1,1) integer :: info , ii , k, IPIV(KM-1),num character :: TRANS !初期化 D = -2.0 DU = 1.0 DL = 1.0 TRANS = 'N' ! ここでLAPACKの三重対格行列を解くサブルーチンをコールする。 call DGTTRF(KM-1,DL,D,DU,DU2,IPIV,info) !$omp parallel do firstprivate(D,DU,DL,DU2,B,IPIV,num,info,TRANS) do ii=1, (NN+1) do k=2,KM B(k-1,1) = x2(ii,k)*dz*dz enddo num =1 call DGTTRS(TRANS,KM-1,num,DL,D,DU,DU2,IPIV,B,KM-1,info) if(info /= 0)then print *,'eroor at lapack' endif do k=2,KM x(ii,k) = B(k-1,1) enddo x(ii,1)=0.0 !境界条件 chi(1)=chi(KM+1)=0 x(ii,KM+1)=0.0 enddo !$omp end parallel do endsubroutine !------------------------------------------------------------------- subroutine S2GUV(x,y,a,b) !スペクトルデータからU,Vのグリッドデータを計算 use constants use wrap_lj implicit none real(8) :: x((NN+1),KM+2) real(8) :: y((NN+1),KM+2) real(8) :: a(JM,KM+2) real(8) :: b(JM,KM+2) real(8) :: g1a(JM,KM+2) , g2a(JM,KM+2) integer :: i , k ,j ,n !$omp parallel do do k = 2, KM+1 call L_S2G_dlat_dual(x(:,k),y(:,k),g1a(:,k),g2a(:,k)) !call LJCS2Y(NN,s1a,sy1a,LJC) !call LJCS2Y(NN,s2a,sy2a,LJC) !call LJMS2G(NN,NM,NN+1,JM,sy1a,sy2a,g1a,g2a,LJIT,LJT,LJP,Q,LJR,WS1,WS2,WG,W1,W2,1) !call LJCS2X(NN,s1b,sx1b) !call LJCS2X(NN,s2b,sx2b) !call LJMS2G(NN,NM,NN,JM,sx1b,sx2b,g1b,g2b,LJIT,LJT,LJP,Q,LJR,WS1,WS2,WG,W1,W2,0) do j=1 , JM a(j,k) = -g1a(j,k)*R*cphi(j) ! Uの生成 b(j,k) = g2a(j,k)*R*cphi(j) ! Vの生成 enddo enddo !$omp end parallel do !$omp parallel do do j=1,JM a(j,1)=alpha*a(j,2) !境界条件 a(j,KM+2)=a(j,KM+1) b(j,1) = alpha*b(j,2) !境界条件 b(j,KM+2) = b(j,KM+1) enddo !$omp end parallel do endsubroutine !------------------------------------------------------------------------------------- subroutine nextpsi(psi,B,kpsi,tstep) use constants use wrap_lj implicit none real(8) :: psi((NN+1),KM+2),kpsi((NN+1),KM+2,4) real(8) :: B(JM,KM+2) real(8) :: lB(JM,KM+2) integer :: tstep real(8) :: s2((NN+1),KM+2) integer :: i , k ,j !$omp parallel do do k=2 , KM+1 do j =1, JM lB(j,k) = B(j,k)*incphi(j) enddo call L_G2S_div(lB(:,k),s2(:,k)) do i=2 , (NN+1) !通常の粘性 kpsi(i,k,tstep+1) = ((s2(i,k))/(R*LJn(i)*(LJn(i)+1.0)) & !& -nu_H*(LJn(i)*(LJn(i)+1.0)-2.0)*psi(i,k)*inR*inR & & +nu_V*(psi(i,k+1)-2.*psi(i,k)+psi(i,k-1))*indz*indz)*dt !4階の超粘性 !kpsi(i,k,tstep+1) = ((s2(i,k))/(R*LJn(i)*(LJn(i)+1.)) - & !&nu_H*(LJn(i)*(LJn(i)+1.0)*LJn(i)*(LJn(i)+1.0)-2.0*LJn(i)*(LJn(i)+1.0))*psi(i,k)*inR*inR +& !&nu_V*(psi(i,k+1)-2.*psi(i,k)+psi(i,k-1))*indz*indz)*dt enddo enddo !$omp end parallel do endsubroutine !--------------------------------------------------------------------------------- subroutine nextchi2z(chi2z,A1z,theta_s,E_s,kchi2z,tstep) use constants use wrap_lj implicit none real(8) :: theta_s((NN+1),KM+2),chi2z((NN+1),KM+1),E_s((NN+1),KM+2) real(8) :: A1z(JM,KM+1),kchi2z((NN+1),KM+1,4) integer :: tstep real(8) :: lA(JM,KM+1) real(8) :: s1((NN+1),KM+1) integer :: i , k ,j !$omp parallel do do k=2 , KM do j =1, JM lA(j,k) = A1z(j,k)*incphi(j) enddo call L_G2S_div(lA(:,k), s1(:,k)) do i=2 , (NN+1) !通常の粘性 kchi2z(i,k,tstep+1) = ((s1(i,k))/(R*LJn(i)*(LJn(i)+1.)) & & -((ga/theta_0)*(theta_s(i,k+1)+theta_s(i,k))*0.5+(E_s(i,k+1)-E_s(i,k))*indz)*inR*inR & !& - nu_H*(2.0*LJn(i)*(LJn(i)+1.)-2.0)*chi2z(i,k)*inR*inR & & +nu_V*(chi2z(i,k+1)-2.*chi2z(i,k)+chi2z(i,k-1))*indz*indz)*dt !4階の超粘性 !kchi2z(i,k,tstep+1) = ((s1(i,k)-s2(i,k))/(R*LJn(i)*(LJn(i)+1.)) & !& -((ga/theta_0)*(theta_s(i,k+1)+theta_s(i,k))*0.5+(E_s(i,k+1)-E_s(i,k))*indz)*inR*inR & !& - nu_H*(2.0*LJn(i)*(LJn(i)+1.)*LJn(i)*(LJn(i)+1.)-2.0*LJn(i)*(LJn(i)+1.))& !&*chi2z(i,k)*inR*inR +& !&nu_V*(chi2z(i,k+1)-2.*chi2z(i,k)+chi2z(i,k-1))*indz*indz)*dt enddo enddo !$omp end parallel do endsubroutine !--------------------------------------------------------------------------------- subroutine nexttheta(theta_s,Vtheta_g,wtheta_s,theta_es,ktheta_s,tstep) use constants use wrap_lj implicit none real(8) :: wtheta_s((NN+1),KM+1) real(8) :: theta_s((NN+1),KM+2) real(8) :: ktheta_s((NN+1),KM+2,4) real(8) :: Vtheta_g(JM,KM+2) real(8) :: lV(JM,KM+2) real(8) :: theta_es((NN+1),KM+2) integer :: tstep real(8) :: s2((NN+1),KM+2) integer :: i , k ,j !$omp parallel do do k=2 , KM+1 do j =1, JM lV(j,k) = Vtheta_g(j,k)*incphi(j) enddo call L_G2S_div(lV(:,k), s2(:,k)) do i=1 , (NN+1) ktheta_s(i,k,tstep+1) = (-(s2(i,k))*inR & & -((wtheta_s(i,k)-wtheta_s(i,k-1))*indz)-(theta_s(i,k)-theta_es(i,k))*intau & !& - kappa_H*LJn(i)*(LJn(i)+1.)*theta_s(i,k)*inR*inR& & +kappa_V*(theta_s(i,k+1)-2.*theta_s(i,k)+theta_s(i,k-1))*indz*indz)*dt enddo enddo !$omp end parallel do endsubroutine !--------------------------------------------------------------------------------- subroutine output(time,U,V,theta_g,w_g,chi,theta_s,psi,chi2z,hst_data1, hst_data2,hst_restart,TotalM) !出力ルーチン use constants use gt4_history use wrap_lj implicit none real(8) :: U(JM,KM+2), lat(JM), hight(KM+2), Uval(JM,KM+2) real(8) :: V(JM,KM+2), Vval(JM,KM+2), M(JM,KM+2) real(8) :: theta_g(JM,KM+2), w_g(JM, KM+1), PV(JM,KM+2) real(8) :: theta_s((NN+1),KM+2), chi(NN+1, KM+1) real(8) :: totalPT,totalM,totalDM,mpsi(JM,KM+1) real(8) :: psi((NN+1),KM+2), chi2z((NN+1),KM+1), chi2zval((NN+1),KM+2) integer :: time, j,k real(8) :: z_half(KM+2), z_int(KM+1) !character :: name,fname type(GT_HISTORY) :: hst_data1, hst_restart, hst_data2 if (time==daystep*INT(output_interval)) then !netcdfに書き込まれるメタ情報。適当に書き換える。 !g95を用いるときの注意:文字列配列の要素文字列の長さをそろえる必要がある。 call HistoryCreate(file='data2D-g1.nc',title='BPM2D',source='BPM2D.f90',institution='user',& & dims=(/'lat ','z ','time'/),dimsizes=(/JM,KM+2,0/),longnames=(/'lat ','z ','time'/), & & units=(/'deg ','m ','days'/), origin=init_day,interval=output_interval, history=hst_data1) do j=1, JM lat(j) = phi(j)*180.0/pi enddo do k=1 ,KM+2 z_half(k) = dz*real(k-1,8) - 0.5*dz enddo call HistoryPut('lat ',lat,hst_data1) call HistoryPut('z ',z_half,hst_data1) call HistoryAddVariable(varname='u',dims=(/'lat ','z ','time'/),& &longname='u',units='m/s',xtype='float',history=hst_data1) call HistoryAddVariable(varname='v',dims=(/'lat ','z ','time'/),& &longname='v',units='m/s',xtype='float',history=hst_data1) call HistoryAddVariable(varname='theta',dims=(/'lat ','z ','time'/),& &longname='theta',units='K',xtype='float',history=hst_data1) call HistoryAddVariable(varname='M',dims=(/'lat ','z ','time'/),& &longname='absolute_angular_momentum',units='m/s',xtype='float',history=hst_data1) call HistoryAddVariable(varname='PV',dims=(/'lat ','z ','time'/),& &longname='potential vorticity',units='K/ms',xtype='float',history=hst_data1) call HistoryAddVariable(varname='totalPT',dims=(/'time'/),& &longname='totalPT',units='K',xtype='float',history=hst_data1) call HistoryAddVariable(varname='totalM',dims=(/'time'/),& &longname='totalM',units='m/s',xtype='float',history=hst_data1) call HistoryAddVariable(varname='totalDM',dims=(/'time'/),& &longname='totalDM',units='MKS',xtype='float',history=hst_data1) call HistoryCreate(file='data2D-g2.nc',title='BPM2D',source='BPM2D.f90',institution='user',& & dims=(/'lat ','z ','time'/),dimsizes=(/JM,KM+1,0/),longnames=(/'lat ','z ','time'/), & & units=(/'deg ','m ','days'/), origin=init_day,interval=output_interval, history=hst_data2) do k=1,KM+1 z_int(k) = real(k-1)*dz enddo call HistoryPut('lat',lat,hst_data2) call HistoryPut('z',z_int,hst_data2) call HistoryAddVariable(varname='w',dims=(/'lat ','z ','time'/),& &longname='w',units='m/s',xtype='float',history=hst_data2) call HistoryAddVariable(varname='mpsi',dims=(/'lat ','z ','time'/),& &longname='meridional_psi',units='m2/s',xtype='float',history=hst_data2) endif call conserves(theta_s,U,totalPT,totalM,totalDM) call makempsi(chi, mpsi) call HistoryPut('totalPT',totalPT,hst_data1) call HistoryPut('totalM',totalM,hst_data1) call HistoryPut('totalDM',totalDM,hst_data1) do j = 1, JM Uval(j,:) = U(j,:)*incphi(j) M(j,:) = U(j,:)*R + R*Omega*Omega*cphi(j)+cphi(j) Vval(j,:) = V(j,:)*incphi(j) enddo call makePV(Uval, theta_g, theta_s, psi, PV) call HistoryPut('u', Uval,hst_data1) call HistoryPut('v', Vval,hst_data1) call HistoryPut('w', w_g,hst_data2) call HistoryPut('theta', theta_g,hst_data1) call HistoryPut('M',M,hst_data1) call HistoryPut('mpsi',mpsi,hst_data2) call HistoryPut('PV',PV,hst_data1) !----------------------------------------- !リスタートファイル if (time==daystep*INT(output_interval)) then call HistoryCreate(file='restart2D.nc',title='BPM2D',source='BPM2D-s.f90',institution='user',& & dims=(/'wn','z '/),dimsizes=(/(NN+1),KM+2/),longnames=(/'wn','z '/), & & units=(/'1/m','m '/),history=hst_restart) call HistoryPut('wn',LJn,hst_restart) call HistoryPut('z ',z_half,hst_restart) call HistoryAddVariable(varname='psi',dims=(/'wn','z '/),& &longname='psi',units='MKS',xtype='double',history=hst_restart) call HistoryAddVariable(varname='chi2z',dims=(/'wn','z '/),& &longname='theta_s',units='MKS',xtype='double',history=hst_restart) call HistoryAddVariable(varname='theta_s',dims=(/'wn','z '/),& &longname='chi2z',units='MKS',xtype='double',history=hst_restart) endif do j=1,(NN+1) do k=1,KM+1 chi2zval(j,k)= chi2z(j,k) enddo chi2zval(j,KM+2)=chi2z(j,KM+1) enddo call HistoryPut('psi',psi,hst_restart) call HistoryPut('chi2z',chi2zval,hst_restart) call HistoryPut('theta_s',theta_s,hst_restart) endsubroutine !-------------------------------------------------------------------- subroutine conserves(theta_s,U,totalPT,totalM,totalDM) use constants use wrap_lj implicit none real(8) :: theta_s((NN+1),KM+2) real(8) :: U(JM,KM+2),U_s((NN+1),KM+2) real(8) :: totalPT,totalM,totalDM integer :: k totalPT = 0.0 totalM = 0.0 !$omp parallel do do k=2,KM+1 call L_G2S(U(:,k), U_s(:,k)) enddo !omp end parallel do do k=2, KM+1 totalPT = totalPT + 2.0*theta_s(1,k)*dz totalM = totalM + 2.0*U_s(1,k)*dz enddo call L_G2S(U(:,1),U_s(:,1)) totalDM = -2.0*C*(U_s(1,2)-U_s(1,1))*indz print *, totalPT,totalM,totalDM endsubroutine !------------------------------------------------------------------------------------- subroutine makempsi(chi,mpsi) !子午面流線関数の計算 use constants use wrap_lj implicit none real(8) :: chi((MM+1),KM+1) real(8) :: mpsi(JM ,KM+1) integer :: i,j,k do k=1,KM+1 call L_S2G_dlat(chi(:,k), mpsi(:,k)) do j=1,JM mpsi(j,k) = mpsi(j,k)*R enddo enddo endsubroutine !------------------------------------------------------------------------------- subroutine makePV(Uval, theta_g,theta_s,psi,PV) ! PVの計算 use constants use wrap_lj implicit none real(8) :: theta_g(JM, KM+2), theta_s(NN+1, KM+2), psi(NN+1, KM+2), PV(JM, KM+2) real(8) :: theta_phi(JM,KM+2), zeta_g(JM,KM+2), zeta_s(NN+1, KM+2), Uval(JM,KM+1) integer :: j, k PV = 0.0 do k =2 , KM+1 call L_S2G_dlat(theta_s(:,k), theta_phi(:,k)) call L_S2S_lap(psi(:,k), zeta_s(:,k)) call L_S2G(zeta_s(:,k), zeta_g(:,k)) do j=1, JM PV(j,k-1) = theta_phi(j,k)*(Uval(j,k+1)-Uval(j,k-1))*0.5*indz*inR & & +(zeta_g(j,k)+2.0*Omega*sphi(j))*(theta_g(j,k+1)-theta_g(j,k-1))*0.5*indz enddo enddo do j =1, JM PV(j,1) = PV(j,2) PV(j,KM+2) = PV(j,KM+1) enddo endsubroutine endprogram