adh2=ckh(i,j,k)*(v(i-1,j,k)+v(i+1,j,k)+v(i,j-1,k)+v(i,j+1,k)-4.*v(i,j,k))/dd/dd*step
   adh3=ckh(i,j,k)*(q(i-1,j,k)+q(i+1,j,k)+q(i,j-1,k)+q(i,j+1,k)-4.*q(i,j,k))/dd/dd*step  
   adh4=ckh(i,j,k)*(t(i-1,j,k)+t(i+1,j,k)+t(i,j-1,k)+t(i,j+1,k)-4.*t(i,j,k))/dd/dd*step  
        
   adv1=(ckz(i,j,k+1)*(u(i,j,k+1)-u(i,j,k))/(z(k+1)-z(k))-ckz(i,j,k)*(u(i,j,k)-u(i,j,k-1))/(z(k)-z(k-1)))/(z(k+1)-z(k-1))*(s/(s-zg(i,j)))**2*step
   adv2=(ckz(i,j,k+1)*(v(i,j,k+1)-v(i,j,k))/(z(k+1)-z(k))-ckz(i,j,k)*(v(i,j,k)-v(i,j,k-1))/(z(k)-z(k-1)))/(z(k+1)-z(k-1))*(s/(s-zg(i,j)))**2*step
   adv3=(cke(i,j,k+1)*(q(i,j,k+1)-q(i,j,k))/(z(k+1)-z(k))-cke(i,j,k)*(q(i,j,k)-q(i,j,k-1))/(z(k)-z(k-1)))/(z(k+1)-z(k-1))*(s/(s-zg(i,j)))**2*step
   adv4=(ckt(i,j,k+1)*(t(i,j,k+1)-t(i,j,k))/(z(k+1)-z(k))-ckt(i,j,k)*(t(i,j,k)-t(i,j,k-1))/(z(k)-z(k-1)))/(z(k+1)-z(k-1))*(s/(s-zg(i,j)))**2*step
  
   u(i,j,k)=u(i,j,k)+adh1+adv1
   v(i,j,k)=v(i,j,k)+adh2+adv2
   q(i,j,k)=q(i,j,k)+adh3+adv3
   t(i,j,k)=t(i,j,k)+adh4+adv4
enddo
enddo
enddo
!$OMP END DO NOWAIT

!SETP 14 U V smoothing
!call smooth(u,l,m,n,aa22,bb22)
!call smooth(v,l,m,n,aa22,bb22)


!$OMP DO PRIVATE(i,j,k)
do k=2,nn
do j=2,mm
do i=2,ll
   u(i,j,k)=aa22*u(i,j,k)+bb22*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k))
   v(i,j,k)=aa22*v(i,j,k)+bb22*(v(i+1,j,k)+v(i-1,j,k)+v(i,j+1,k)+v(i,j-1,k))
enddo
enddo
enddo
!$OMP END DO NOWAIT

!$OMP END PARALLEL


WRITE(*,"(5F10.4,I6)")U(10,10,2),V(10,10,2),W(10,10,4),T(10,10,3),p(10,10,3),mk
!WRITE(*,"(6F10.4)")QK(10,10),QR(10,10),SHEAT(10,10),LHEAT(10,10),TSOIL(10,10),TT
WRITE(*,"(4F10.4)")CKZ(10,10,2),CKH(10,10,2),Q(10,10,2)
write(*,"(4f10.4)")qk(10,10),qr(10,10),qh(10,10),qe(10,10)



!!!!一次積分完成，得到新的u,v,w,t,p,q，積分步數增加1，返回55，再次積分
mk=mk+1


goto 55

1994	stop

end

!  sub top                              
!  This program is for dealing with absorbing layer in the five highest grid points near the top of the model                            
subroutine top (a,l,m,n)
dimension a(l,m,n)
do k=n-4,n
   tk=k
do j=1,m
do i=1,l
   a(i,j,k)=a(i,j,k)-(a(i,j,k)-a(1,1,k))*(tk-16.)/4.
enddo
enddo
enddo
return
end 

!  sub  pp                          
!  This program is for calculating pressure gradient force term     

subroutine pp(u,v,p,t,zg,z,dd,s,step,l,m,n)
dimension u(l,m,n),v(l,m,n),t(l,m,n),p(l,m,n),z(n),zg(l,m)

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
return
end

!  sub  advuu                         
!  this program is for calculating advective terms of momentum  u,v   
subroutine advuu(u,v,w,z,dd,step,l,m,n)
dimension u(l,m,n),v(l,m,n),w(l,m,n),z(n)



do 1 k=2,kk
dz1=z(k+1)-z(k)
dz=z(k)-z(k-1)
do 1 j=2,mm     
do 1 i=2,ll
vv=0.25*(v(i-1,j,k)+v(i-1,j+1,k)+v(i,j,k)+v(i,j+1,k))

if(vv.ge.0.)  then
    u(i,j,k)=u(i,j,k)-vv*(u(i,j,k)-u(i,j-1,k))/dd*step
else
    u(i,j,k)=u(i,j,k)-vv*(u(i,j+1,k)-u(i,j,k))/dd*step
endif

if(u(i,j,k).ge.0.) then
    u(i,j,k)=u(i,j,k)-u(i,j,k)*(u(i,j,k)-u(i-1,j,k))/dd*step
else
    u(i,j,k)=u(i,j,k)-u(i,j,k)*(u(i+1,j,k)-u(i,j,k))/dd*step
endif  

ww=.25*(w(i-1,j,k)+w(i,j,k)+w(i-1,j,k+1)+w(i,j,k+1))

if(ww.ge.0.) then
    u(i,j,k)=u(i,j,k)-ww*(u(i,j,k)-u(i,j,k-1))/dz*step
else
    u(i,j,k)=u(i,j,k)-ww*(u(i,j,k+1)-u(i,j,k))/dz1*step
endif

tt=tt
uu=.25*(u(i,j-1,k)+u(i+1,j-1,k)+u(i,j,k)+u(i+1,j,k))

if(uu.ge.0.) then
    v(i,j,k)=v(i,j,k)-uu*(v(i,j,k)-v(i-1,j,k))/dd*step
else
    v(i,j,k)=v(i,j,k)-uu*(v(i+1,j,k)-v(i,j,k))/dd*step
endif

if(v(i,j,k).ge.0.) then
    v(i,j,k)=v(i,j,k)-v(i,j,k)*(v(i,j,k)-v(i,j-1,k))/dd*step
else
    v(i,j,k)=v(i,j,k)-v(i,j,k)*(v(i,j+1,k)-v(i,j,k))/dd*step
endif

ww=.25*(w(i,j-1,k)+w(i,j,k)+w(i,j-1,k+1)+w(i,j,k+1))

if(ww.ge.0.) then
    v(i,j,k)=v(i,j,k)-ww*(v(i,j,k)-v(i,j,k-1))/dz*step
else
    v(i,j,k)=v(i,j,k)-ww*(v(i,j,k+1)-v(i,j,k))/dz1*step
endif

1     continue
return
end



!  sub adv
!  this program is for caculating advective terms  of scalar variables          
subroutine adv(u,v,w,a,z,dd,step,l,m,n)
dimension u(l,m,n),v(l,m,n),w(l,m,n),a(l,m,n),z(n)

do 1 k=2,nn
dz=z(k)-z(k-1)
dz1=z(k+1)-z(k)
do 1 j=2,mm
do 1 i=2,ll

uu=.5*(u(i,j,k)+u(i+1,j,k))
vv=.5*(v(i,j,k)+v(i,j+1,k))
ww=.5*(w(i,j,k)+w(i,j,k+1))
if (uu.ge.0.) then 
    a(i,j,k)=a(i,j,k)-uu*(a(i,j,k)-a(i-1,j,k))/dd*step
else
    a(i,j,k)=a(i,j,k)-uu*(a(i+1,j,k)-a(i,j,k))/dd*step
endif

if(vv.ge.0.) then 
    a(i,j,k)=a(i,j,k)-vv*(a(i,j,k)-a(i,j-1,k))/dd*step
else
    a(i,j,k)=a(i,j,k)-vv*(a(i,j+1,k)-a(i,j,k))/dd*step
endif

if(ww.ge.0.) then 
    a(i,j,k)=a(i,j,k)-ww*(a(i,j,k)-a(i,j,k-1))/dz*step
else
    a(i,j,k)=a(i,j,k)-ww*(a(i,j,k+1)-a(i,j,k))/dz1*step
endif 

1     continue
return  
end


!  sub smooth                     
!  this progam  is used to smooth calculation noise caused by diference schemes      
subroutine smooth(a,l,m,n,aa11,bb11)
dimension a(l,m,n)

do k=2,nn
do j=2,mm
do i=2,ll
   a(i,j,k)=aa11*a(i,j,k)+bb11*(a(i+1,j,k)+a(i-1,j,k)+a(i,j+1,k)+a(i,j-1,k))
enddo
enddo
enddo
return
end


!  sub verti                        
!  this program is for calculating continuity equation      
subroutine verti(u,v,w,dd,step,s,z,zg,l,m,n,ll,mm,kk)
dimension u(l,m,n),v(l,m,n),w(l,m,n),z(n),zg(l,m)
do k=2,kk
do j=2,mm
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
return
end

!  sub diff                         
!  this program is for calculating diffusion terms             
subroutine diff(ck,ckh,z,zg,dd,step,s,l,m,n,a)
dimension a(l,m,n),z(n),zg(l,m),ck(l,m,n),ckh(l,m,n)


!$OMP PARALLEL 
!$OMP DO PRIVATE(i,j,k)
do k=2,nn
do j=2,mm
do i=2,ll
   adh=ckh(i,j,k)*(a(i-1,j,k)+a(i+1,j,k)+a(i,j-1,k)+a(i,j+1,k)-4.*a(i,j,k))/dd/dd*step
   adv=(ck(i,j,k+1)*(a(i,j,k+1)-a(i,j,k))/(z(k+1)-z(k))-ck(i,j,k)*(a(i,j,k)-a(i,j,k-1))/(z(k)-z(k-1)))/(z(k+1)-z(k-1))*(s/(s-zg(i,j)))**2*step
   a(i,j,k)=a(i,j,k)+adh+adv
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP END PARALLEL

return 
end

!  sub bound                             
!  this program is for dealing with the lateral boundaries         
subroutine bound(a,u,v,l,m,n)
dimension a(l,m,n),u(l,m,n),v(l,m,n)
do 1 k=2,n
do 1 j=1,m
if(u(1,j,k).le.0.) then
    a(1,j,k)=a(2,j,k)
else
    a(l,j,k)=a(l-1,j,k)
endif
1   continue


do 3 k=2,n
do 3 i=1,l
if(v(i,2,k).lt.0.) then
    a(i,1,k)=a(i,2,k)
elseif(v(i,m-1,k).gt.0.) then
    a(i,m,k)=a(i,m-1,k)
endif
3   continue
return
end



! calculating turbulence energy      
subroutine turbu(q,ckz,ckt,cke,u,v,t,z,step,l,m,n,tm,s,zg)
dimension  q(l,m,n),ckz(l,m,n),u(l,m,n),v(l,m,n),t(l,m,n),zg(l,m)
dimension  ckt(l,m,n),cke(l,m,n)
dimension  z(n)
al8=35.


!$OMP PARALLEL
!$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
alf=1.25
else
sm=1.96*(.1912-rif)*(.2341-rif)/(1.-rif)/(.2331-rif)
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)))*alf*ckt(i,j,k)*9.8/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)*alf
!if(i.ne.12.or.k.ne.3) goto 1
1   continue
!$OMP END DO NOWAIT
!$OMP END PARALLEL

!  up and down boundery conditions
do j=1,m
do i=1,l
ckt(i,j,n)=ckt(i,j,n-1)
enddo
enddo

return
end

subroutine interface(qk,qr,qh,qe,u2,t1,t2,r2,ct,z0,alf,z2,step,qa,tm,tem,rstar,hu,soil,water,grass,trees,tsoil,i,j)

parameter(as=0.2,taus=0.9	,b = 9.4,pai=3.1416,rs=40.,laitre=3.,laigra=2.,pa=11.,fai=32.5/180.*3.1416)

s0 = 1353!

delta=-23.5*3.1416/180*cos(2*3.1416*(173+16)/365.)


!	richardson number
ri = 9.8 / t2 * z2 * (t2 - t1) / u2 / u2
a2=0.16/(log(z2/z0))**2
c=7.4*a2*9.4*(z2/z0)**.5
if(ri.le.0.) then
fm=1.-9.4*ri/(1.+c*(abs(ri))**.5)
else
fm=1./(1.+4.7*ri)**2
endif
us=(a2)**0.5*u2*fm**.5
c=5.3*a2*b*(z2/z0)**0.5
if(ri.le.0.) then
fh=1.-b*ri/(1.+c*(abs(ri))**0.5)
else
fh=fm
endif

qh = a2/.74*u2*(t2-t1)*fh
conduc=a2/0.74*u2*fh

ra=1./conduc

	 
qe =conduc*water*(r2-rstar)+grass*(r2-rstar)/(ra+rs/laigra)+trees*(r2-rstar)/(ra+rs/laitre)
if(r2.le.hu*rstar) then
qe=qe+conduc*soil*(r2-rstar)
else
endif
ts=qh/us
if(ts.eq.0.) then
al=-100!
else
al=t2*us*us/9.8/ts/.4
end if
qh=1.29*1005.*qh
qe=1.29*2500000.*qe


ea=1.24*(pa/t2)**.14
qr=(ea*t2**4-0.95*t1**4)*5.67/10.**8
tk=tm
h=3.1416*(12.-tk)/12.  !Sang orig
!h=(tk-12)*15.*3.1416/180.


!delta=0.409*cos(2.*3.1416*(180.-173.)/365.) !Sang orig

csz=sin(fai)*sin(delta)+cos(fai)*cos(delta)*cos(h)
!qk=(1.-alf)*s0*csz*taus**(1./csz) !Sang orig
if(csz.ge.0.000001) then
!qk=s0*csz*(1-5*(0.28/(1-6.43*csz)))*(1.-alf) !xujm
qk=(1.-alf)*s0*csz*taus**(1./csz)
else
qk=0.
endif
t1=t1+ct*(qk+qr+qh+qe+qa)*step-2.*pai/24./3600.*(t1-tsoil)*step
if(i.eq.10 .and. j.eq.10) then
!write(*,*)"qk=",qk,"t1=",t1,"tk=",tk
endif
tsoil=tsoil+(t1-tsoil)/24./3600.*step
return
end

subroutine landclas(ct,alf,z0,zg,bldmean,soil,water,grass,trees,l,m)
dimension ct(l,m),alf(l,m),z0(l,m),zg(l,m),bldmean(l,m)
dimension soil(l,m),water(l,m),grass(l,m),trees(l,m)
!dimension landuse_orig(400,223,271)

integer l,m
do j=1,m
do i=1,l
   num_water=10
   num_busin=100
   num_inhab=100
   num_trans=90
   num_green=40
   num_garden=40
   num_farm=20
   num_orch=0
   num_tree=0
   num_beach=0
   num_lack=0
   
   alf(i,j)=(0.09*num_water+0.10*num_busin+0.12*num_inhab+0.2*num_green+0.19*num_garden+0.19*num_farm+0.18*num_orch+0.18*num_tree+0.16*num_beach)/(400-num_lack)
   z0(i,j)=(0.001*num_water+2.5*num_busin+2*num_inhab+0.5*num_green+0.5*num_garden+0.1*num_farm+0.1*num_orch+3.*num_tree+0.005*num_beach)/(400-num_lack)
   zg(i,j)=(0.2*num_water+3.5*num_busin+3.5*num_inhab+3.5*num_green+3.5*num_garden+4.*num_farm+3.5*num_orch+3.5*num_tree+1.*num_beach)/(400-num_lack)
   ct(i,j)=(0.08*num_water+0.35*num_busin+0.3*num_inhab+0.2*num_green+0.16*num_garden+0.12*num_farm+0.13*num_orch+0.12*num_tree+0.14*num_beach)/100000/(400-num_lack)
   bldmean(i,j)=(0.0*num_water+18.*num_busin+15.*num_inhab+0.0*num_green+1.*num_garden+1.*num_farm+0.*num_orch+0.*num_tree+0.*num_beach)/(400-num_lack)
   water(i,j)=(1.0*num_water+0.*num_busin+0.*num_inhab+0.*num_green+0.*num_garden+0.*num_farm+0.*num_orch+0.*num_tree+0.4*num_beach)/(400-num_lack)
   trees(i,j)=(0.*num_water+0.1*num_busin+0.2*num_inhab+0.3*num_green+0.5*num_garden+0.*num_farm+1.*num_orch+1.*num_tree+0.*num_beach)/(400-num_lack)
   grass(i,j)=(0.*num_water+0.1*num_busin+0.2*num_inhab+1.*num_green+0.8*num_garden+0.3*num_farm+0.3*num_orch+0.*num_tree+0.*num_beach)/(400-num_lack)
   soil(i,j)=(0.*num_water+0.9*num_busin+0.6*num_inhab+0.*num_green+0.2*num_garden+0.8*num_farm+0.*num_orch+0.*num_tree+0.5*num_beach)/(400-num_lack)

enddo
enddo

end