mp-sang.f90

一个数值计算程序的OpenMP的并行化代码

!  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