36 (u,v,w,p,t,c_u,c_v,c_w,c_p,c_ff,c_gg,c_hh,c_qu,c_qv,c_qw)
45 real(mytype),
dimension(:,:,:),
intent(in) :: u,v,w,p,t
46 complex(mytype),
dimension(:,:,:),
intent(in) :: c_u,c_v,c_w,c_p
47 complex(mytype),
dimension(:,:,:),
intent(out) :: c_ff,c_gg,c_hh
48 complex(mytype),
dimension(:,:,:),
intent(out) :: c_qu,c_qv,c_qw
51 real(mytype) :: tmp1,tmp2,a1,a2
52 real(mytype) :: coe_sk
70 c_qu(:,:,:)=c_u(:,:,:)
83 tmp1= ( (u(i,j,k+1)-u(i,j,k))/
dz_t(k) - &
84 (u(i,j,k)-u(i,j,k-1))/
dz_b(k) ) /
dz(k)
98 c_qu(:,:,:)=c_qu(:,:,:)+
rk_b*
dt*c_ff(:,:,:)
103 c_ff(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
109 c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*
c_var(:,:,:)
112 c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*
c_var(:,:,:)
118 a1=u(i,j,k+1)*
l_t(k)+u(i,j,k)*(1.d0-
l_t(k))
119 a2=u(i,j,k-1)*
l_b(k)+u(i,j,k)*(1.d0-
l_b(k))
120 tmp1 = ( a1 * w(i,j,k) - &
121 a2 * w(i,j,k-1) ) /
dz(k)
124 tmp2 =
fc*( v(i,j,k) -
vg )
138 c_ff(:,:,:)=c_ff(:,:,:)+coe_sk*
c_var(:,:,:)
142 if (
issk .eq. 1)
then 153 cmplx( u(i,j,k) *
real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
160 c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*
c_var(:,:,:)
171 real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
178 c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*
c_var(:,:,:)
185 a1=u(i,j,k+1)*
l_t(k)+u(i,j,k)*(1.d0-
l_t(k))
186 a2=u(i,j,k-1)*
l_b(k)+u(i,j,k)*(1.d0-
l_b(k))
187 tmp1 = ( a1 - a2 ) /
dz(k)
188 tmp2 = 0.5d0*( w(i,j,k) + w(i,j,k-1) )
195 c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*
c_var(:,:,:)
201 c_qu(:,:,:)=c_qu(:,:,:)+
rk_a*
dt*c_ff(:,:,:)
212 c_qv(:,:,:)=c_v(:,:,:)
225 tmp1 = ( (v(i,j,k+1)-v(i,j,k))/
dz_t(k) - &
226 (v(i,j,k)-v(i,j,k-1))/
dz_b(k) ) /
dz(k)
240 c_qv(:,:,:)=c_qv(:,:,:)+
rk_b*
dt*c_gg(:,:,:)
245 c_gg(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
251 c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*
c_var(:,:,:)
254 c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*
c_var(:,:,:)
260 a1=v(i,j,k+1)*
l_t(k)+v(i,j,k)*(1.d0-
l_t(k))
261 a2=v(i,j,k-1)*
l_b(k)+v(i,j,k)*(1.d0-
l_b(k))
262 tmp1 = ( a1 * w(i,j,k) - &
263 a2 * w(i,j,k-1) ) /
dz(k)
272 cmplx( 1.d0/coe_sk*tmp2-tmp1,0.d0, kind=mytype )
278 c_gg(:,:,:)=c_gg(:,:,:)+coe_sk*
c_var(:,:,:)
283 if (
issk .eq. 1)
then 294 real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
301 c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*
c_var(:,:,:)
312 real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
319 c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*
c_var(:,:,:)
326 a1=v(i,j,k+1)*
l_t(k)+v(i,j,k)*(1.d0-
l_t(k))
327 a2=v(i,j,k-1)*
l_b(k)+v(i,j,k)*(1.d0-
l_b(k))
328 tmp1 = ( a1 - a2 ) /
dz(k)
329 tmp2 = 0.5d0*( w(i,j,k) + w(i,j,k-1) )
336 c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*
c_var(:,:,:)
342 c_qv(:,:,:)=c_qv(:,:,:)+
rk_a*
dt*c_gg(:,:,:)
353 c_qw(:,:,:)=c_w(:,:,:)
366 tmp1 = ( (w(i,j,k+1)-w(i,j,k))/
dz(k+1) - &
367 (w(i,j,k)-w(i,j,k-1))/
dz(k) ) /
dz_t(k)
382 tmp1= ( p(i,j,k+1)-p(i,j,k) ) /
dz_t(k)
393 c_qw(:,:,:)=c_qw(:,:,:)+
rk_b*
dt*c_hh(:,:,:)
398 c_hh(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
404 c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*
c_var(:,:,:)
407 c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*
c_var(:,:,:)
414 a1 = ( 0.5d0*(w(i,j,k)+w(i,j,k+1)) )**2.d0
415 a2 = ( 0.5d0*(w(i,j,k)+w(i,j,k-1)) )**2.d0
416 tmp1= (a1-a2) /
dz_t(k)
424 cmplx( 1.d0/coe_sk*tmp2 - tmp1,0.d0, kind=mytype )
433 c_hh(:,:,:)=c_hh(:,:,:)+coe_sk*
c_var(:,:,:)
437 if (
issk .eq. 1)
then 447 tmp1=
l_t(k)*u(i,j,k+1)+(1.d0-
l_t(k))*u(i,j,k)
449 real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
456 c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*
c_var(:,:,:)
466 tmp1=
l_t(k)*v(i,j,k+1)+(1.d0-
l_t(k))*v(i,j,k)
468 real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
475 c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*
c_var(:,:,:)
482 a1 = 0.5d0*(w(i,j,k)+w(i,j,k+1))
483 a2 = 0.5d0*(w(i,j,k)+w(i,j,k-1))
484 tmp1= (a1-a2) /
dz_t(k)
486 cmplx( w(i,j,k)*tmp1,0.d0, kind=mytype )
492 c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*
c_var(:,:,:)
499 c_qw(:,:,:)=c_qw(:,:,:)+
rk_a*
dt*c_hh(:,:,:)
510 (u,v,w,t,c_t,c_ss,c_qt)
519 real(mytype),
dimension(:,:,:),
intent(in) :: u,v,w,t
520 complex(mytype),
dimension(:,:,:),
intent(in) :: c_t
521 complex(mytype),
dimension(:,:,:),
intent(out) :: c_ss
522 complex(mytype),
dimension(:,:,:),
intent(out) :: c_qt
525 real(mytype) :: tmp1,a1,a2
534 c_qt(:,:,:)=c_t(:,:,:)
547 tmp1= ( (t(i,j,k+1)-t(i,j,k))/
dz_t(k) - &
548 (t(i,j,k)-t(i,j,k-1))/
dz_b(k) ) /
dz(k)
558 c_qt(:,:,:)=c_qt(:,:,:)+
rk_b*
dt*c_ss(:,:,:)
562 c_ss(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
568 c_ss(:,:,:)=c_ss(:,:,:)-
c_var(:,:,:)
571 c_ss(:,:,:)=c_ss(:,:,:)-
c_var(:,:,:)
577 a1=t(i,j,k+1)*
l_t(k)+t(i,j,k)*(1.d0-
l_t(k))
578 a2=t(i,j,k-1)*
l_b(k)+t(i,j,k)*(1.d0-
l_b(k))
579 tmp1 = ( a1 * w(i,j,k) - &
580 a2 * w(i,j,k-1) ) /
dz(k)
587 c_ss(:,:,:)=c_ss(:,:,:)-
c_var(:,:,:)
590 c_qt(:,:,:)=c_qt(:,:,:)+
rk_a*
dt*c_ss(:,:,:)
603 (aa,bb,cc,c_src_in,c_var_out,zstagger,delnull,tbnd,ubnd)
606 use,
intrinsic :: iso_c_binding
615 include
'include/fftw3.f03' 618 real(mytype),
dimension(:,:,:),
intent(inout) :: aa,cc
619 complex(mytype),
dimension(:,:,:),
intent(inout) :: bb
621 complex(mytype),
dimension(:,:,:),
intent(in) :: c_src_in
631 integer,
intent(in) :: zstagger,delnull,tbnd,ubnd
633 complex(mytype),
dimension(:,:,:),
intent(out) :: c_var_out
636 complex(mytype) :: tmp
638 real(mytype) :: aa_t_zstr,cc_t_zend
639 real(mytype) :: aa_u_zstr,cc_u_zend
640 integer :: i,j,k,kmax3
643 c_var_out(:,:,:)=0.d0
644 c_src(:,:,:)=c_src_in(:,:,:)
648 if (zstagger .eq. 1)
then 650 elseif (zstagger .eq. 0)
then 653 print *,
'Invalid zstagger flag!' 659 if (delnull .eq. 1)
then 661 if (
myid .eq. 0)
then 664 c_src(1,1,1)=cmplx(0.d0,0.d0 , kind=mytype )
666 elseif (delnull .eq. 0)
then 669 print *,
'delnull flag error!' 673 if (tbnd .eq. 1)
then 675 if (
myid .eq. 0)
then 683 elseif (tbnd .eq. 0)
then 686 print *,
'tbnd flag error!' 691 if (ubnd .eq. 1)
then 693 if (
myid .eq. 0)
then 703 elseif (ubnd .eq. 0)
then 706 print *,
'ubnd flag error!' 716 tmp=aa(i,j,k)/bb(i,j,k-1)
717 bb(i,j,k)=bb(i,j,k)-cc(i,j,k-1)*tmp
727 c_var_out(i,j,kmax3-
ghst)=
c_src(i,j,kmax3-
ghst)/bb(i,j,kmax3)
732 do k=kmax3-1,
kstr3,-1
737 cc(i,j,k)*c_var_out(i,j,k-
ghst+1) ) / bb(i,j,k)
762 complex(mytype),
dimension(:,:,:),
intent(in) :: c_u,c_v
763 real(mytype),
dimension(:,:,:),
intent(in) :: w
764 complex(mytype),
dimension(:,:,:),
intent(out) :: c_div
772 c_div(:,:,:)=
c_var(:,:,:)
776 c_div(:,:,:)=c_div(:,:,:)+
c_var(:,:,:)
783 dwdz = ( w(i,j,k)-w(i,j,k-1) )/
dz(k)
790 c_div(:,:,:)=c_div(:,:,:)+
c_var(:,:,:)
809 real(mytype),
dimension(:,:,:),
intent(in) :: phi
810 complex(mytype),
dimension(:,:,:),
intent(in) :: c_phi
811 complex(mytype),
dimension(:,:,:),
intent(inout) :: c_u,c_v,c_w
814 real(mytype) :: dphidz
829 dphidz= ( phi(i,j,k+1) - phi(i,j,k) ) /
dz_t(k)
862 real(mytype),
dimension(:,:,:),
intent(in) :: phi
863 complex(mytype),
dimension(:,:,:),
intent(in) :: c_phi
864 complex(mytype),
dimension(:,:,:),
intent(inout) :: c_p
867 real(mytype) :: dphidz2
870 c_p(:,:,:)=c_p(:,:,:) + c_phi(:,:,:)
882 dphidz2= ( (phi(i,j,k+1)-phi(i,j,k))/
dz_t(k) - &
883 (phi(i,j,k)-phi(i,j,k-1))/
dz_b(k) ) /
dz(k)
908 real(mytype),
dimension(:,:,:),
intent(in) :: u,v,w
909 complex(mytype),
dimension(:,:,:),
intent(in) :: c_u,c_v
911 integer,
intent(in) :: isfinal,ishistout
914 real(mytype) :: div_avg, div_max
915 real(mytype) :: div_avg_all,div_max_all
916 real(mytype) :: div_hist_max_all,tmp1
929 tmp1=
real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype)
931 div_avg=div_avg+abs(tmp1)
933 if ( abs(tmp1) .gt. div_max )
then 937 if ( u(i,j,k) .ne. u(i,j,k) .or. &
938 v(i,j,k) .ne. v(i,j,k) .or. &
939 w(i,j,k) .ne. w(i,j,k) )
then 940 print *,
'NaN error! ',i,
' ',j,
' ',k
944 if (isfinal .eq. 1)
then 957 call mpi_reduce(div_avg,div_avg_all,1,real_type,mpi_sum,0, &
960 call mpi_reduce(div_max,div_max_all,1,real_type,mpi_max,0, &
963 call mpi_reduce(
div_hist_max,div_hist_max_all,1,real_type,mpi_max,0, &
967 if (
myid .eq. 0)
then 969 if ( isfinal .eq. 0)
then 971 write (*,104) div_avg_all,div_max_all
972 104
format(
' DIV-ERR INIT : DIV_AVG = ',es11.4,
'; DIV_MAX = ',es11.4)
974 elseif ( isfinal .eq. 1 .and. ishistout .eq. 0)
then 976 write (*,105) div_avg_all,div_max_all
977 105
format(
' DIV-ERR FINAL: DIV_AVG = ',es11.4,
'; DIV_MAX = ',es11.4)
979 elseif ( isfinal .eq. 1 .and. ishistout .eq. 1)
then 981 write (*,115) div_avg_all,div_max_all
982 115
format(
' DIV-ERR FINAL: DIV_AVG = ',es11.4,
'; DIV_MAX = ',es11.4)
983 write (*,106) div_hist_max_all
984 106
format(
' DIV-ERR HIST_MAX = ', es11.4)
989 258
format(
' Ri_b = ', es12.5)
994 259
format(
' Ubar = ', es12.5)
998 print *,
'div out flag error!' 1013 (u,v,w,t,mean_u,mean_v,mean_w,mean_t,mean_uu,mean_vv,mean_ww, &
1014 mean_uw,mean_vw,mean_tt,mean_tw)
1018 real(mytype),
dimension(:,:,:),
intent(in) :: u,v,w,t
1019 real(mytype),
dimension(:,:,:),
intent(out) :: mean_u,mean_v,mean_w,mean_t, &
1020 mean_uu,mean_vv,mean_ww,mean_uw, mean_vw, mean_tt,mean_tw
1029 mean_uu(i,j,k)=mean_uu(i,j,k)+(u(i,j,k)+
u_mrf)**2.d0
1030 mean_vv(i,j,k)=mean_vv(i,j,k)+v(i,j,k)**2.d0
1031 mean_ww(i,j,k)=mean_ww(i,j,k)+w(i,j,k)**2.d0
1032 mean_uw(i,j,k)=mean_uw(i,j,k)+0.5d0*( w(i,j,k)+w(i,j,k-1) ) * &
1034 mean_vw(i,j,k)=mean_vw(i,j,k)+0.5d0*( w(i,j,k)+w(i,j,k-1) ) * v(i,j,k)
1036 mean_tt(i,j,k)=mean_tt(i,j,k)+t(i,j,k)**2.d0
1037 mean_tw(i,j,k)=mean_tw(i,j,k)+t(i,j,k) * 0.5d0 * ( w(i,j,k)+w(i,j,k-1) )
1040 mean_u(i,j,k)=mean_u(i,j,k) + u(i,j,k)+
u_mrf 1041 mean_v(i,j,k)=mean_v(i,j,k) + v(i,j,k)
1042 mean_w(i,j,k)=mean_w(i,j,k) + w(i,j,k)
1043 mean_t(i,j,k)=mean_t(i,j,k) + t(i,j,k)
1067 real(mytype),
dimension(:,:,:),
intent(in) :: r_in1,r_in2
1068 complex(mytype),
dimension(:,:,:),
intent(out) :: c_out3
1069 integer,
intent(in) :: dir,weqn
1071 integer :: i,j,k, l, m
1072 real(mytype) :: tmp1
1079 if (weqn .eq. 1)
then 1084 c_out3(i,j,k)=cmplx( tmp1 * &
1122 if (dir .eq. 1)
then 1124 c_out3(i,j,k) = c_out3(i,j,k)*
iii*2.d0*
pi/
lx*l
1126 elseif (dir .eq. 2)
then 1128 c_out3(i,j,k) = c_out3(i,j,k)*
iii*2.d0*
pi/
ly*m
1132 print *,
'flag error' 1138 if (l .eq. -
nx/2 .or. m .eq. -
ny/2)
then 1139 c_out3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
1163 real(mytype),
dimension(:,:,:),
intent(in) :: r_in
1164 complex(mytype),
dimension(:,:,:),
intent(out) :: c_out
1172 c_out(i,j,k)=cmplx(r_in(i+
ghst,j+
ghst,k+
ghst),0.d0,kind=mytype)
1196 real(mytype),
dimension(:,:,:),
intent(out) :: r_out
1197 complex(mytype),
dimension(:,:,:),
intent(in) :: c_in
1201 c_var(:,:,:)=c_in(:,:,:)
1209 r_out(i+
ghst,j+
ghst,k+
ghst)=
real(c_var(i,j,k),kind=mytype)
1226 complex(mytype),
dimension(:,:,:),
intent(inout) :: c_inout3
1228 integer :: i,j,k,l,m
1248 if (l .ge. int(
nx/3.d0) .or. l .le. -int(
nx/3.d0) )
then 1250 c_inout3(i,j,k) = cmplx(0.d0,0.d0,kind=mytype)
1255 if (m .ge. int(
ny/3.d0) .or. m .le. -int(
ny/3.d0) )
then 1257 c_inout3(i,j,k) = cmplx(0.d0,0.d0,kind=mytype)
1278 complex(mytype),
dimension(:,:,:),
intent(in) :: c_in3
1279 complex(mytype),
dimension(:,:,:),
intent(out) :: c_out3
1280 integer,
intent(in) :: ord,dir
1282 integer :: i,j,k,l,m
1284 c_out3(:,:,:)=c_in3(:,:,:)
1305 if (ord .eq. 1 .and. dir .eq. 1)
then 1307 c_out3(i,j,k) = c_out3(i,j,k)*
iii*2.d0*
pi/
lx*l
1309 elseif (ord .eq. 1 .and. dir .eq. 2)
then 1311 c_out3(i,j,k) = c_out3(i,j,k)*
iii*2.d0*
pi/
ly*m
1313 elseif (ord .eq. 2 .and. dir .eq. 1)
then 1315 c_out3(i,j,k) = -c_out3(i,j,k)*(2.d0*
pi/
lx*l)**2.d0
1317 elseif (ord .eq. 2 .and. dir .eq. 2)
then 1319 c_out3(i,j,k) = -c_out3(i,j,k)*(2.d0*
pi/
ly*m)**2.d0
1323 print *,
'flag error' 1329 if (l .eq. -
nx/2 .or. m .eq. -
ny/2)
then 1330 c_out3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
1351 complex(mytype),
dimension(:,:,:),
intent(inout) :: c_varin
1353 real(mytype) :: zd,sigz
1354 complex(mytype) :: mean_local,mean_all
1368 mean_local=mean_local/zsize(1)/zsize(2)
1372 call mpi_allreduce(mean_local,mean_all,1,real_type,mpi_sum, &
1373 mpi_comm_world,
ierr)
1374 mean_all=mean_all/
nprc 1378 sigz=
cdamp*0.5d0*( 1.d0-cos(
pi*zd) )
integer, save, protected is_ri_var
integer, save, protected myid
real(mytype), save, protected rk_b
real(mytype), dimension(:), allocatable, save, protected l_b
integer, save, protected ochannel
integer, save, protected nprc
real(mytype), save, protected ubar
integer, save, protected ny
real(mytype), save, protected fc
subroutine real_2_cmplx(r_in, c_out)
integer, save, protected issk
complex(mytype), dimension(:,:,:), allocatable, save c_src
real(mytype), parameter pi
real(mytype), save, protected dt
subroutine get_temperature_cnvdiff_spec(u, v, w, t, c_t, c_ss, c_qt)
real(mytype), save, protected rk_a
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ttop
integer, save, protected nz
real(mytype), save, protected ly
real(mytype), save, protected ug
subroutine calculate_mean_fields_spec(u, v, w, t, mean_u, mean_v, mean_w, mean_t, mean_uu, mean_vv, mean_ww, mean_uw, mean_vw, mean_tt, mean_tw)
subroutine spectral_truncation(c_inout3)
real(mytype), dimension(:), allocatable, save, protected dz
integer, save, protected isscalar
real(mytype), save, protected cdamp
real(mytype), save, protected alpha
integer, save, protected istr3
real(mytype), save, protected vg
real(mytype), save, protected t_ref
subroutine get_momentum_cnvdiff_spec(u, v, w, p, t, c_u, c_v, c_w, c_p, c_ff, c_gg, c_hh, c_qu, c_qv, c_qw)
subroutine rayleigh_damping_spec(c_varin)
integer, save, protected jstr3
real(mytype), dimension(:), allocatable, save, protected dz_b
integer, save, protected i_offset
integer, save, protected j_offset
subroutine cmplx_2_real(c_in, r_out)
real(mytype), save, protected u_mrf
complex(mytype), parameter iii
real(mytype), save, protected lx
real(mytype), save, protected dpx_drive
integer, save, protected nzdamp
real(mytype), dimension(:), allocatable, save, protected dz_t
subroutine poisson_solver_fft_spec(aa, bb, cc, c_src_in, c_var_out, zstagger, delnull, tbnd, ubnd)
integer, save, protected nx
subroutine dev_fft_linear_cmplx(c_in3, c_out3, ord, dir)
subroutine update_dpx_drive
subroutine screen_div_error_spec(u, v, w, c_u, c_v, isfinal, ishistout)
real(mytype), save, protected got
subroutine get_div_spec(c_u, c_v, w, c_div)
real(mytype), save div_hist_max
subroutine calculate_new_p_spec(phi, c_phi, c_p)
subroutine calculate_uvw_spec(phi, c_phi, c_u, c_v, c_w)
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), dimension(:), allocatable, save, protected l_t
real(mytype), save, protected tbot
complex(mytype), dimension(:,:,:), allocatable, save c_var
integer, save, protected dp_opt
subroutine dev_fft_nonlinear_cmplx(r_in1, r_in2, c_out3, dir, weqn)
integer, save, protected kstr3
integer, save, protected kend3