?? mp-sang.f90
字號:
! main program mesoscale model
! this program is used to simulate three-dimensional flow fields over urban boundary layer
! z -- height in terrain following coordinate system
! u,v,w -- wind speeds
! p -- pressure ( exner function )
! t -- potential temperature
! q -- turbulence energy
! ckh -- horizontal diffusion coefficient
! ckz -- vertical diffusion coefficient for momentum
! ckt -- vertical diffusion coefficient for heat
! cke -- vertical diffusion coefficient for turbulent energy
! zg -- topography
! h -- height of mixed layer
! uu,vv,ww -- wind velocities for adpic
! qkh -- horizontal diffusion coefficient in adpic
! zgg -- topography in adpic
! qkv -- vertical diffusion coefficient in adpic
! a,b -- auxiliary variables
program pgl
include 'omp_lib.h'
!integer,EXTERNAL :: OMP_SET_NUM_THREADS
integer :: Itera_steps, Threadsnum1, Threadsnum2
parameter (l=223,m=271,n=20) !設置變量維數
parameter(Itera_steps=100,Threadsnum1=4,Threadsnum2=4)
parameter(ll=l-1,mm=m-1,kk=n-1,nn=n-1)
parameter(dd=500. ,cp=1005.)
parameter(dep=300. ,dtl=0.01 ,dtu=0.005)
parameter(step=2 ) !積分步長
parameter(aa22=0.5 ) !平滑系數
parameter(bb22=0.125 ) !平滑系數
!dimension u(l,m,n),v(l,m,n),w(l,m,n),p(l,m,n),t(l,m,n) !定義風速、氣壓和位溫
!dimension z(n),zg(l,m),q(l,m,n) !z表示垂直高度，zg表示地形高度，q為湍流動能
!dimension a(l,m,n),b(l,m,n) !臨時變量
!dimension ckz(l,m,n),ckt(l,m,n),cke(l,m,n),ckh(l,m,n) !湍流交換系數
!dimension tsoil(l,m) !土壤溫度
!dimension trees(l,m),grass(l,m),soil(l,m),water(l,m),bldmean(l,m) !樹木、草地、土壤、水占網格的百分比，建筑高度
!dimension ct(l,m),alf(l,m),z0(l,m) !地表參數：熱力系數、反照率、粗糙度
!dimension qk(l,m),qr(l,m),qh(l,m),qe(l,m) !短波輻射、長波輻射、感熱、潛熱
dimension z(n);
real,dimension(:,:,:),allocatable::u,v,w,p,t,q
real,dimension(:,:),allocatable::zg
real,dimension(:,:,:),allocatable::a,b
real,dimension(:,:,:),allocatable::ckz,ckt,cke,ckh
real,dimension(:,:),allocatable::tsoil
real,dimension(:,:),allocatable::trees,grass,soil,water,bldmean
real,dimension(:,:),allocatable::ct,alf,z0
real,dimension(:,:),allocatable::qk,qr,qh,qe
allocate(u(l,m,n));allocate(v(l,m,n));allocate(w(l,m,n));allocate(p(l,m,n));allocate(t(l,m,n));
allocate(zg(l,m));allocate(q(l,m,n));
allocate(a(l,m,n));allocate(b(l,m,n));
allocate(ckz(l,m,n));allocate(ckt(l,m,n));allocate(cke(l,m,n));allocate(ckh(l,m,n));
allocate(tsoil(l,m));
allocate(trees(l,m));allocate(grass(l,m));allocate(soil(l,m));allocate(water(l,m));allocate(bldmean(l,m));
allocate(ct(l,m));allocate(alf(l,m));allocate(z0(l,m));
allocate(qk(l,m));allocate(qr(l,m));allocate(qh(l,m));allocate(qe(l,m));
! l*m*n -- grid point
! z() -- vertical coordinate in terrain-following system
!!!定義垂直層高度
z(1)=0.
z(2)=10.
z(3)=20.
z(4)=50.
z(5)=80.
z(6)=100.
z(7)=150.
z(8)=200.
z(9)=250.
z(10)=300.
z(11)=400.
z(12)=500.
z(13)=750.
z(14)=1000.
z(15)=1250.
z(16)=1500.
z(17)=2000.
z(18)=2500.
z(19)=3000.
z(20)=4000.
!!!!調用子程序計算地表參數和地表分類：ct,alf,z0,zg,bldmean
call landclas(ct,alf,z0,zg,bldmean,soil,water,grass,trees,l,m)
! s -- top height of the model
! dd -- horizontal grid interval
s=z(20)
! initial conditions of u,v,t
! dep --- depth of the lower layer
! dtl --- vertical gradient of potential temperature in lower layer
! dtu --- vertical gradient of potential temperature in upper layer
!!!!!!初始化溫度、風速和湍流動能，t,u,v,q
call OMP_SET_NUM_THREADS(Threadsnum1)
!$OMP PARALLEL PRIVATE(TID,i,j,k)
TID = OMP_GET_THREAD_NUM()
PRINT *, 'Hello World from thread = ', TID
!$OMP DO
do k=1,n
do j=1,m
do i=1,l
tsoil(i,j)=290.
u(i,j,k)=-1.*(z(k)/10.)**0.12
v(i,j,k)=-1.5*(z(k)/10.)**0.14
zz=zg(i,j)+(s-zg(i,j))/s*z(k)
if (zz.le.dep) then
t(i,j,k)=290.+dtl*zz
else
t(i,j,k)=293.+dtu*(zz-dep)
end if
! qs(i,j)=1.*bldmean(i,j)
if(k.le.14) then
q(i,j,k)=0.1
else
q(i,j,k)=0.05
endif
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO
! p(x,y,s) is the pressure at the top of model
!!!!!!定義頂層氣壓p
do j=1,m
do i=1,l
p(i,j,20)=875.
enddo
enddo
!$OMP END DO NOWAIT
!$OMP END PARALLEL
! step -- time step every circulating
! mk -- step number of model integration
! aa22 and bb22 are coef of smooth
mk=0 !積分步數
tt=6. !初始積分時間
mk=mk+1
! integrating hydrostatic equation to obtain pressure fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!積分過程!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
55 if(mk.gt.Itera_steps) goto 1994
!ierr=OMP_SET_NUM_THREADS(Threadsnum2)
!!$OMP PARALLEL PRIVATE(TID)
!TID = OMP_GET_THREAD_NUM()
!PRINT *, 'PHASE II from thread = ', TID
!!$OMP END PARALLEL
!!!!step 1:計算各層的氣壓p
!$OMP PARALLEL
!$OMP DO PRIVATE(i,j,k)
do j=1,m
do k=kk,1,-1
do i=1,l
p(i,j,k)=p(i,j,k+1)+9.8*(s-zg(i,j))/s/(t(i,j,k)+t(i,j,k+1))*(z(k+1)-z(k))*2.
enddo
enddo
enddo
!$OMP END DO NOWAIT
! absorbing layer for prevention of wave reflection from model top
!!!!step2:頂層濾波u,v,t
!call top(u,l,m,n)
!call top(v,l,m,n)
!call top(t,l,m,n)
!$OMP DO PRIVATE(i,j,k,tk)
do k=n-4,n
tk=(k-16.)/4.
do j=1,m
do i=1,l
u(i,j,k)=u(i,j,k)*(1-tk)+u(1,1,k)*tk
v(i,j,k)=v(i,j,k)*(1-tk)+v(1,1,k)*tk
t(i,j,k)=t(i,j,k)*(1-tk)+t(1,1,k)*tk
enddo
enddo
enddo
!$OMP END DO NOWAIT
! time filtering for noise damping at boundary
!!!!step3:處理t的側邊界，同時將u,v放入臨時變量a,b
!$OMP DO PRIVATE(i,j,k)
do k=1,n
do j=1,m
t(1,j,k)=t(2,j,k)
t(l,j,k)=t(l-1,j,k)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO PRIVATE(i,j,k)
do k=1,n
do i=1,l
t(i,1,k)=t(i,2,k)
t(i,m,k)=t(i,m-1,k)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO PRIVATE(i,j,k)
do k=1,n
do j=1,m
do i=1,l
a(i,j,k)=u(i,j,k)
b(i,j,k)=v(i,j,k)
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP SINGLE
!!!!step4:計算u,v的平流過程，更新u,v
call advuu(u,v,w,z,dd,step,l,m,n)
!$OMP END SINGLE
!!!!step5:計算u,v的氣壓梯度過程，更新u,v
!call pp(u,v,p,t,zg,z,dd,s,step,l,m,n)
!$OMP DO PRIVATE(i,j,k)
do k=2,kk
do j=2,mm
do i=2,ll
u(i,j,k)=u(i,j,k)-(p(i,j,k)-p(i-1,j,k))/dd*step*0.5*(t(i,j,k)+t(i-1,j,k))+(zg(i,j)-zg(i-1,j))/dd*9.8*(z(k)-s)/s*step
v(i,j,k)=v(i,j,k)-(p(i,j,k)-p(i,j-1,k))/dd*step*0.5*(t(i,j,k)+t(i,j-1,k))+(zg(i,j)-zg(i,j-1))/dd*9.8*(z(k)-s)/s*step
enddo
enddo
enddo
!$OMP END DO NOWAIT
!!!!step5:用新的u,v和原來的u,v結合，平滑擾動，同時將t和q放入臨時變量t,q
!$OMP DO PRIVATE(i,j,k)
do k=1,n
do j=1,m
do i=1,l
u(i,j,k)=0.6*u(i,j,k)+0.4*a(i,j,k)
v(i,j,k)=0.6*v(i,j,k)+0.4*b(i,j,k)
a(i,j,k)=t(i,j,k)
b(i,j,k)=q(i,j,k)
enddo
enddo
enddo
!$OMP SECTIONS
!!!!step6:計算t,q的平流項，更新t,q
!$OMP SECTION
call adv(u,v,w,t,z,dd,step,l,m,n)
!$OMP SECTION
call adv(u,v,w,q,z,dd,step,l,m,n)
!$OMP END SECTIONS
!!!!step7:計算垂直速度w
!call verti(u,v,w,dd,step,s,z,zg,l,m,n,ll,mm,kk)
!$OMP DO PRIVATE(i,j,k)
do j=2,mm
do k=2,kk
do i=2,ll
w(i,j,k)=w(i,j,k-1)-(u(i+1,j,k-1)-u(i,j,k-1)+v(i,j+1,k-1)-v(i,j,k-1))/dd*(z(k)-z(k-1))+.5*((u(i+1,j,k-1)+u(i,j,k-1))*(zg(i+1,j)-zg(i-1,j))+(v(i,j+1,k-1)+v(i,j,k-1))*(zg(i,j+1)-zg(i,j-1)))/dd/2.*(z(k)-z(k-1))/(s-zg(i,j))
enddo
enddo
enddo
!$OMP END DO NOWAIT
!!!!step8:處理u,v,t的邊界
!$OMP SECTIONS
!$OMP SECTION
call bound(u,u,v,l,m,n)
!$OMP SECTION
call bound(v,u,v,l,m,n)
!$OMP SECTION
call bound(t,u,v,l,m,n)
!$OMP SECTION
tt=tt+step/3600.
!$OMP END SECTIONS
!!!!step9:調用interface子程序，計算地表溫度t(i,j,1)
!$OMP DO PRIVATE(i,j,k) PRIVATE(z2,t1,t2,r2,tem,estar,rstar,hu)
do j=1,m
do i=1,l
z2=z(2)
t1=t(i,j,1)
t2=t(i,j,2)
u2=(u(i,j,2)*u(i,j,2)+v(i,j,2)*v(i,j,2))**0.5
r2=0.010
tem=t(i,j,1)*p(i,j,1)/cp
estar=6.112*exp(17.67*(tem-273.15)/(tem-29.65))
rstar=0.622*estar/1000.
hu=0.5
call interface(qk(i,j),qr(i,j),qh(i,j),qe(i,j),u2,t1,t2,r2,ct(i,j),z0(i,j),alf(i,j),z2,step,qa,tt,tem,rstar,hu,soil(i,j),water(i,j),grass(i,j),trees(i,j),tsoil(i,j),i,j)
t(i,j,1)=t1
enddo
enddo
!$OMP END DO NOWAIT
!!!!step10:調用turbu子程序計算湍流擴散系數ckz,ckt,cke
!call turbu(q,ckz,ckt,cke,u,v,t,z,step,l,m,n,tt,s,zg)
al8=35.
!$OMP DO PRIVATE(i,j,k,dz,ri,rif,sm,alf,all,ti,tk)
do 1 k=2,n-1
dz=z(k+1)-z(k-1)
do 1 j=1,m
do 1 i=1,l
if(u(i,j,k+1).eq.u(i,j,k-1)) goto 10
ri=9.8/t(i,j,k)*(t(i,j,k+1)-t(i,j,k-1))*dz/((u(i,j,k+1)-u(i,j,k-1))**2+(v(i,j,k+1)-v(i,j,k-1))**2)
if(ri.ge.0.195) goto 10
rif=.6588*(ri+.1776-(ri*ri-.3221*ri+.03156)**.5)
goto 11
10 rif=.191
11 if(rif.ge.0.16) then
sm=.085
temp_alf=1.25
else
sm=1.96*(.1912-rif)*(.2341-rif)/(1.-rif)/(.2331-rif)
temp_alf=1.318*(.2331-rif)/(.2341-rif)
end if
all=.4*z(k)/(1.+.4*z(k)/al8)
!ti=i
!tk=k
!if(q(i,j,k).le.0.) then
if(q(i,j,k).le.0.00001) then
q(i,j,k)=.01
else
q(i,j,k)=q(i,j,k)+ckz(i,j,k)*(((u(i,j,k+1)-u(i,j,k-1))/dz)**2+((v(i,j,k+1)-v(i,j,k-1))/dz)**2)*step*(s/(s-zg(i,j)))**2-(s/(s-zg(i,j)))*ckt(i,j,k)*9.8*temp_alf/t(i,j,k)*(t(i,j,k+1)-t(i,j,k-1))/dz*step*0.5-(2*q(i,j,k))**1.5/10.6/all*step
end if
! if(q(i,j,k).le.0.00001) q(i,j,k)=0.01
ckz(i,j,k)=sm*all*(2*q(i,j,k))**.5
cke(i,j,k)=.2*all*(2*q(i,j,k))**.5
ckt(i,j,k)=ckz(i,j,k)*temp_alf
!if(i.ne.12.or.k.ne.3) goto 1
1 continue
!$OMP END DO NOWAIT
! up and down boundery conditions
!$OMP DO PRIVATE(i,j)
do j=1,m
do i=1,l
ckt(i,j,n)=ckt(i,j,n-1)
enddo
enddo
!$OMP END DO NOWAIT
!!!!step11:平滑t,q,濾掉擾動
!$OMP DO PRIVATE(i,j,k)
do k=1,n
do j=1,m
do i=1,l
t(i,j,k)=0.6*t(i,j,k)+0.4*a(i,j,k)
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO PRIVATE(i,j,k)
do k=1,n
do j=1,m
do i=1,l
q(i,j,k)=0.6*q(i,j,k)+0.4*b(i,j,k)
enddo
enddo
enddo
!$OMP END DO NOWAIT
!!!!step12:計算水平擴散系數ckh
!$OMP DO PRIVATE(i,j,k)
do k=2,n
do j=2,mm
do i=2,ll
def=(v(i+1,j,k)-v(i-1,j,k)+u(i,j+1,k)-u(i,j-1,k))**2+(u(i+1,j,k)-u(i-1,j,k)-v(i,j+1,k)+v(i,j-1,k))**2
def=def**0.5
ckh(i,j,k)=0.049*dd*def
enddo
enddo
enddo
!$OMP END DO NOWAIT
!!!!step13:利用湍流擴散系數計算u,v,q,t的湍流擴散項，更新u,v,q,t
!call diff(ckz,ckh,z,zg,dd,step,s,l,m,n,u)
!call diff(ckz,ckh,z,zg,dd,step,s,l,m,n,v)
!call diff(cke,ckh,z,zg,dd,step,s,l,m,n,q)
!call diff(ckt,ckh,z,zg,dd,step,s,l,m,n,t)
!$OMP DO PRIVATE(i,j,k)
do k=2,nn
do j=2,mm
do i=2,ll
adh1=ckh(i,j,k)*(u(i-1,j,k)+u(i+1,j,k)+u(i,j-1,k)+u(i,j+1,k)-4.*u(i,j,k))/dd/dd*step
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -