doc/tec/evaluation: a1_a3.ncl

; Concepts illustrated:
;   - Drawing a scatter plot over a map using the "overlay" procedure
;   - Using gsn_csm_blank_plot to create a scatter plot with filled polygons
;   - Generating dummy data using "random_uniform"
;   - Changing the draw order of filled polygons
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"





begin


;Indizes ÃŒberdenken, besonders bei interpolation (sehe a4_2.f90 an)

nx=335    
ny=83    
nz=64    

dx=2.5
dy=2.5






;Nahbereich a1-1

;Beginn
nxnb=72
nynb=3
nznb=1


;Ende
nxne=172
nyne=ny-4
nzne=16


nyh=floattointeger((ny+1)/2.-1)       ;Mitte y-Richtung

;Erlaubte Abweichungen
W=0.01    ;a1-1, a1-2
D=0.05

Wc1=0.07   ;c1
Dc1=0.25




;Verschiebung des Gitters in m
sdx=230. 
sdy=105.


;Normierungspunkt
xn=-70.
yn=0.
zn=75.


findin=0.
findjn=0.
findkn=0.

upol=new((/nz,ny,nx/),double)
xupol=new((/nz,ny,nx/),double)
yupol=new((/nz,ny,nx/),double)
zupol=new((/nz,ny,nx/),double)

upol4=new((/nz,ny,nx/),double)
xupol4=new((/nz,ny,nx/),double)
yupol4=new((/nz,ny,nx/),double)
zupol4=new((/nz,ny,nx/),double)

upol5=new((/2*nz,2*ny,2*nx/),double)
xupol5=new((/2*nz,2*ny,2*nx/),double)
yupol5=new((/2*nz,2*ny,2*nx/),double)
zupol5=new((/2*nz,2*ny,2*nx/),double)

xspeicher=new((/nx*nz/),double)
ipos=new((/nx*nz/),integer)
kpos=new((/nx*nz/),integer)



;Testfaelle auswaehlen -> 1: ja, 0: nein
test_a1_1 = 1
test_a1_2 = 1    ;wenn test_a1_2=1, dann muss test_a1_1=1 sein
test_a2   = 1    ;wenn test_a2=1, dann muss test_a1_2=1 sein
test_a3_1 = 1    ;wenn test_a3_1=1, dann muss test_a1_2=1 sein
test_a3_2 = 1    ;wenn test_a3_2=1, dann muss test_a3_1=1 sein

test_c1   = 1    ;wenn test_c1=1, dann muss test_a1_2=1 sein

print(" ")
;---------------------------------------------------------------
;----Einlesen der Daten--------------------------------------------
;---------------------------------------------------------------
if(test_a1_1 .ge. 1) then
f10 = addfile("../OUTPUT/a1-1_av_3d.nc","r")
dNames = getfilevardims(f10,"time")
dSizes = getfilevardimsizes(f10,"time")
tmax=dSizes-1
t=tmax
uwind=f10->u(t,:,:,:)       ;u(t,k,j,i)
wwind=f10->w(t,:,:,:)       ;w(t,k,j,i)
DKM=f10->km(t,:,:,:)
hoehe=f10->zu_3d(:)
xu=f10->xu(:)-sdx
y=f10->y(:)-sdy
x=f10->x(:)-sdx
hoehew=f10->zw_3d(:)
end if

if(test_a1_2 .ge. 1) then
f2 = addfile("../OUTPUT/a1-2_av_3d.nc","r")  ;Name der Datei Àndern
dNames = getfilevardims(f2,"time")
dSizes = getfilevardimsizes(f2,"time")
tmax=dSizes-1
t=tmax
uwind1=f2->u(t,:,:,:)       ;u(t,k,j,i)
wwind1=f2->w(t,:,:,:)       ;w(t,k,j,i)
nznea1_2=26  ; bis 75m Hoehe
hoehe1=f2->zu_3d(:)
xu1=f2->xu(:)-sdx
y1=f2->y(:)-sdy
x1=f2->x(:)-sdx
hoehew1=f2->zw_3d(:)

end if


if(test_a2 .ge. 1) then
f3 = addfile("../OUTPUT/a2_av_3d.nc","r")  ;Name der Datei Àndern
dNames = getfilevardims(f3,"time")
dSizes = getfilevardimsizes(f3,"time")
tmax=dSizes-1
t=tmax
uwind2=f3->u(t,:,:,:)       ;u(t,k,j,i)
wwind2=f3->w(t,:,:,:)       ;w(t,k,j,i)
hoehe2=f3->zu_3d(:)
xu2=f3->xu(:)-sdx
y2=f3->y(:)-sdy
x2=f3->x(:)-sdx
hoehew2=f3->zw_3d(:)
nznea1_2=26  ; bis 75m Hoehe

end if


if(test_a3_2 .ge. 1) then
f4 = addfile("../OUTPUT/a3-2_av_3d.nc","r")  ;Name der Datei Àndern
dNames = getfilevardims(f4,"time")
dSizes = getfilevardimsizes(f4,"time")
tmax=dSizes-1
t=tmax
uwind3=f4->u(t,:,:,:)       ;u(t,k,j,i)
wwind3=f4->w(t,:,:,:)       ;w(t,k,j,i)


end if






if(test_c1 .ge. 1) then
ascii_filename1 = "../comparative_data/C1.dat"
seismic1 = asciiread(ascii_filename1,(/651,7/),"float")
 
vdi=651
xvdi=new((/vdi/),double)
yvdi=new((/vdi/),double)
zvdi=new((/vdi/),double)
unormvdi=new((/vdi/),double)
wnormvdi=new((/vdi/),double)

xvdi(:) = seismic1(:,1) 
yvdi(:) = seismic1(:,2) 
zvdi(:) = seismic1(:,3)
unormvdi(:) = seismic1(:,4)
wnormvdi(:)  = seismic1(:,6)
end if
  
;---------------------------------------------------------------
;----Trilineare Interpolation (siehe wikipedia)--------------------------------------------
;---------------------------------------------------------------


;----Am Normierungspunkt--------------------------------------------

;---fuer a1-1-----------
 if(test_a1_1 .eq. 1) then
do i=0,nx

if(findin .le.0.) then
if(xu(i) .ge. xn) then
in=i
findin=findin+1.

end if
end if
end do


do j=0,ny
if(findjn .le. 0.) then
if(y(j) .ge. yn) then
jn=j
findjn=findjn+1.
end if
end if
end do



do k=0,nz
if(findkn .le. 0.) then
if(hoehe(k) .ge. zn) then
kn=k
findkn=findkn+1.
end if
end if
end do





xd=(xn-(xu(in-1)))/(xu(in)-(xu(in-1)))
yd=(yn-(y(jn-1)))/(y(jn)-(y(jn-1)))
zd=(zn-hoehe(kn-1))/(hoehe(kn)-hoehe(kn-1))




 c00=uwind(kn-1,jn-1,in-1)*(1.-xd)+uwind(kn-1,jn-1,in)*xd
 c01=uwind(kn,jn-1,in-1)*(1.-xd)+uwind(kn,jn-1,in)*xd
 c10=uwind(kn-1,jn,in-1)*(1.-xd)+uwind(kn-1,jn,in)*xd
 c11=uwind(kn,jn,in-1)*(1.-xd)+uwind(kn,jn,in)*xd


 c0=c00*(1.-yd)+c10*yd
 c1=c01*(1.-yd)+c11*yd

 c=c0*(1.-zd)+c1*zd
 unorm=c

print("Die Windgeschwindigkeit am Normierungspunkt betraegt bei a1-1: "+c+" m/s.")
print(" ")

end if
 ;---------fuer a1-2--------------------
 if(test_a1_2 .eq. 1) then
 
 findin=0.
 findjn=0.
 findkn=0.
 
 do i=0,nx

if(findin .le.0.) then
if(xu1(i) .ge. xn) then
in=i
findin=findin+1.

end if
end if
end do


do j=0,ny
if(findjn .le. 0.) then
if(y1(j) .ge. yn) then
jn=j
findjn=findjn+1.
end if
end if
end do



do k=0,nz
if(findkn .le. 0.) then
if(hoehe1(k) .ge. zn) then
kn=k
findkn=findkn+1.
end if
end if
end do
 
xd=(xn-(xu1(in-1)))/(xu1(in)-(xu1(in-1)))
yd=(yn-(y1(jn-1)))/(y1(jn)-(y1(jn-1)))
zd=(zn-hoehe1(kn-1))/(hoehe1(kn)-hoehe1(kn-1))


 c00=uwind1(kn-1,jn-1,in-1)*(1.-xd)+uwind1(kn-1,jn-1,in)*xd
 c01=uwind1(kn,jn-1,in-1)*(1.-xd)+uwind1(kn,jn-1,in)*xd
 c10=uwind1(kn-1,jn,in-1)*(1.-xd)+uwind1(kn-1,jn,in)*xd
 c11=uwind1(kn,jn,in-1)*(1.-xd)+uwind1(kn,jn,in)*xd


 c0=c00*(1.-yd)+c10*yd
 c1=c01*(1.-yd)+c11*yd

 c=c0*(1.-zd)+c1*zd
 
 unorm1=c
 
 print("Die Windgeschwindigkeit am Normierungspunkt betraegt bei a1-2: "+c+" m/s.")
print(" ")
 
 end if
 
 ;---------fuer a2--------------------
 if(test_a2 .eq. 1) then     ;!!!!!!hoehe,xu,y anpassen 
 
 
xd=(xn-(xu2(in-1)))/(xu2(in)-(xu2(in-1)))
yd=(yn-(y2(jn-1)))/(y2(jn)-(y2(jn-1)))
zd=(zn-hoehe2(kn-1))/(hoehe2(kn)-hoehe2(kn-1))


 c00=uwind2(kn-1,jn-1,in-1)*(1.-xd)+uwind2(kn-1,jn-1,in)*xd
 c01=uwind2(kn,jn-1,in-1)*(1.-xd)+uwind2(kn,jn-1,in)*xd
 c10=uwind2(kn-1,jn,in-1)*(1.-xd)+uwind2(kn-1,jn,in)*xd
 c11=uwind2(kn,jn,in-1)*(1.-xd)+uwind2(kn,jn,in)*xd


 c0=c00*(1.-yd)+c10*yd
 c1=c01*(1.-yd)+c11*yd

 c=c0*(1.-zd)+c1*zd
 
 unorm2=c
 
  
 print("Die Windgeschwindigkeit am Normierungspunkt betraegt bei a2: "+c+" m/s.")
print(" ")
 
 end if






;---------------------------------------------------------------
;----Gitterpunkte fÃŒr VDI Datei C1 bestimmen--------------------------------------------
;---------------------------------------------------------------
if(test_c1 .ge. 1) then
;Initialisierung
U_zaehl=0.
W_zaehl=0.
DKM_zaehl=0.
U_treffer=0.
W_treffer=0.
;----fuer u---------------------


do jj=0,vdi-1 ;Schleife ÃŒber einzelne Vergleichspunkte



findin=0.
findjn=0.
findkn=0.

xn=doubletofloat(xvdi(jj))
yn=doubletofloat(yvdi(jj))
zn=doubletofloat(zvdi(jj))

do i=1,nx
if(findin .le.0.) then
if(xu1(i) .ge. xn) then
in=i
findin=findin+1.

end if
end if
end do




do j=1,ny
if(findjn .le. 0.) then
if(y1(j) .ge. yn) then
jn=j
findjn=findjn+1.
end if
end if
end do




do k=1,nz
if(findkn .le. 0.) then
if(hoehe1(k) .ge. zn) then
kn=k
findkn=findkn+1.
end if
end if
end do





xd=(xn-(xu1(in-1)))/(xu1(in)-(xu1(in-1)))
yd=(yn-(y1(jn-1)))/(y1(jn)-(y1(jn-1)))
zd=(zn-hoehe1(kn-1))/(hoehe1(kn)-hoehe1(kn-1))


 c00=uwind1(kn-1,jn-1,in-1)*(1.-xd)+uwind1(kn-1,jn-1,in)*xd
 c01=uwind1(kn,jn-1,in-1)*(1.-xd)+uwind1(kn,jn-1,in)*xd
 c10=uwind1(kn-1,jn,in-1)*(1.-xd)+uwind1(kn-1,jn,in)*xd
 c11=uwind1(kn,jn,in-1)*(1.-xd)+uwind1(kn,jn,in)*xd


 c0=c00*(1.-yd)+c10*yd
 c1=c01*(1.-yd)+c11*yd

 c=c0*(1.-zd)+c1*zd
 
 

 
 
 

 
;----Abweichung von vdi-Wert-----------------------------------------------------
if(.not. ismissing(c/unorm1) .and. .not. ismissing(unormvdi(jj)) .and. unormvdi(jj) .ne. 0.) then

Dc=abs((c/unorm1-unormvdi(jj))/unormvdi(jj))
Wc=abs(c/unorm1-unormvdi(jj))

U_zaehl=U_zaehl+1

if(Dc .le. Dc1 .or. Wc .le. Wc1) then
       U_treffer=U_treffer+1.
end if
end if


end do



print("Trefferquote fuer U bei c1: " +U_treffer/U_zaehl*100.+" %")


;----fuer w  (anderes Gitter)---------------------


do jj=0,vdi-1 ;Schleife ÃŒber einzelne Vergleichspunkte


findin=0.
findjn=0.
findkn=0.

xn=doubletofloat(xvdi(jj))
yn=doubletofloat(yvdi(jj))
zn=doubletofloat(zvdi(jj))

do i=1,nx
if(findin .le.0.) then
if(x1(i) .ge. xn) then
in=i
findin=findin+1.

end if
end if
end do




do j=1,ny
if(findjn .le. 0.) then
if(y1(j) .ge. yn) then
jn=j
findjn=findjn+1.
end if
end if
end do




do k=1,nz
if(findkn .le. 0.) then
if(hoehew1(k) .ge. zn) then
kn=k
findkn=findkn+1.
end if
end if
end do


xd=(xn-(x1(in-1)))/(x1(in)-(x1(in-1)))
yd=(yn-(y1(jn-1)))/(y1(jn)-(y1(jn-1)))
zd=(zn-hoehew1(kn-1))/(hoehew1(kn)-hoehew1(kn-1))



 c00=wwind1(kn-1,jn-1,in-1)*(1.-xd)+wwind1(kn-1,jn-1,in)*xd
 c01=wwind1(kn,jn-1,in-1)*(1.-xd)+wwind1(kn,jn-1,in)*xd
 c10=wwind1(kn-1,jn,in-1)*(1.-xd)+wwind1(kn-1,jn,in)*xd
 c11=wwind1(kn,jn,in-1)*(1.-xd)+wwind1(kn,jn,in)*xd


 c0=c00*(1.-yd)+c10*yd
 c1=c01*(1.-yd)+c11*yd

 c=c0*(1.-zd)+c1*zd
 


 
;----Abweichung von vdi-Wert-----------------------------------------------------

if(.not. ismissing(c/unorm1) .and. .not. ismissing(wnormvdi(jj)) .and. wnormvdi(jj) .ne. 0.) then

Dc=abs((c/unorm1-wnormvdi(jj))/wnormvdi(jj))
Wc=abs(c/unorm1-wnormvdi(jj))

W_zaehl=W_zaehl+1

if(Dc .le. Dc1 .or. Wc .le. Wc1) then
       W_treffer=W_treffer+1.
end if

end if

end do

print("Trefferquote fuer W bei c1: " +W_treffer/W_zaehl*100.+" %")

;-----Testfall bestanden?----------------------------

if(U_treffer/U_zaehl*100..gt. 66.  .and. W_treffer/W_zaehl*100..gt. 66. ) then
print("Testfall c1 bestanden")
end if

print(" ")
end if
;------------------------------------------------------------------









;---------------------------------------------------------------
;----Vergleich Homogen--------------------------------------------
;---------------------------------------------------------------

if(test_a1_1 .eq. 1) then    ;!!!!!!!!!!!!!!!!!!TEST
;Initialisierung
U_zaehl=0.
W_zaehl=0.
DKM_zaehl=0.
U_treffer=0.
W_treffer=0.
DKM_treffer=0.



do k=1,nz
do j=1,ny
do i=1,nx


if(j .ne. nyh) then
;----u-Komponente-----------
if(.not. ismissing(uwind(k,j,i)) .and. .not. ismissing(uwind(k,nyh,i))) then

Dc=abs((uwind(k,j,i)-uwind(k,nyh,i))/uwind(k,nyh,i))
Wc=abs(uwind(k,j,i)-uwind(k,nyh,i))

U_zaehl=U_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       U_treffer=U_treffer+1.
end if

end if


;;----w-Komponente---------
if(.not. ismissing(wwind(k,j,i)) .and. .not. ismissing(wwind(k,nyh,i))) then

Dc=abs((wwind(k,j,i)-wwind(k,nyh,i))/wwind(k,nyh,i))
Wc=abs(wwind(k,j,i)-wwind(k,nyh,i))

W_zaehl=W_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       W_treffer=W_treffer+1.
end if

end if

;----DKM--------------
if(.not. ismissing(DKM(k,j,i)) .and. .not. ismissing(DKM(k,nyh,i)) .and. DKM(k,nyh,i) .ne. 0.) then

Dc=abs((DKM(k,j,i)-DKM(k,nyh,i))/DKM(k,nyh,i))
Wc=abs(DKM(k,j,i)-DKM(k,nyh,i))

DKM_zaehl=DKM_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       DKM_treffer=DKM_treffer+1.
end if

end if


end if

end do
end do
end do



print("Trefferquote fuer U bei a1-1: " +U_treffer/U_zaehl*100.+" %")


print("Trefferquote fuer W bei a1-1: " +W_treffer/W_zaehl*100.+" %")


print("Trefferquote fuer DKM bei a1-1: " +DKM_treffer/DKM_zaehl*100.+" %")


if(U_treffer/U_zaehl*100..gt. 95. .and. W_treffer/W_zaehl*100..gt. 95. .and. DKM_treffer/DKM_zaehl*100..gt. 95.) then
print("Testfall a1-1 bestanden")
end if

print(" ")
end if


;---------------------------------------------------------------
;----Vergleich Skalierbarkeit a1-2--------------------------------------------
;---------------------------------------------------------------

if(test_a1_2 .eq. 1) then
;Initialisierung
U_zaehl=0.
W_zaehl=0.
DKM_zaehl=0.
U_treffer=0.
W_treffer=0.

do k=nznb,nznea1_2
do i=nxnb,nxne




;u-Komponente
if(.not. ismissing(uwind(k,nyh,i)/unorm) .and. .not. ismissing(uwind1(k,nyh,i)/unorm1)) then


unorma1_1=doubletofloat(uwind(k,nyh,i)/unorm)
unorma1_2=doubletofloat(uwind1(k,nyh,i)/unorm1)



Dc=abs((unorma1_2-unorma1_1)/unorma1_1)
Wc=abs(unorma1_2-unorma1_1)

U_zaehl=U_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       U_treffer=U_treffer+1.
end if

end if


;;w-Komponente
if(.not. ismissing(wwind(k,nyh,i)/unorm) .and. .not. ismissing(wwind1(k,nyh,i)/unorm1)) then


wnorma1_1=doubletofloat(wwind(k,nyh,i)/unorm)
wnorma1_2=doubletofloat(wwind1(k,nyh,i)/unorm1)

Dc=abs((wnorma1_2-wnorma1_1)/wnorma1_1)
Wc=abs(wnorma1_2-wnorma1_1)

W_zaehl=W_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       W_treffer=W_treffer+1.
end if

end if

end do
end do




print("Trefferquote fuer U bei a1-2: " +U_treffer/U_zaehl*100.+" %")



print("Trefferquote fuer W bei a1-2: " +W_treffer/W_zaehl*100.+" %")


if(U_treffer/U_zaehl*100..gt. 95. .and. W_treffer/W_zaehl*100..gt. 95.) then
print("Testfall a1-2 bestanden")
end if
print(" ")
end if


;---------------------------------------------------------------
;----Vergleich Stationaritaet a2--------------------------------------------
;---------------------------------------------------------------

if(test_a2 .eq. 1) then
;Initialisierung
U_zaehl=0.
W_zaehl=0.
DKM_zaehl=0.
U_treffer=0.
W_treffer=0.


unorm1=1.
unorm2=1.

do k=nznb,nznea1_2
do i=nxnb,nxne




;u-Komponente
if(.not. ismissing(uwind1(k,nyh,i)/unorm1) .and. .not. ismissing(uwind2(k,nyh,i)/unorm2)) then


unorma2=doubletofloat(uwind2(k,nyh,i)/unorm2)
unorma1_2=doubletofloat(uwind1(k,nyh,i)/unorm1)



Dc=abs((unorma2-unorma1_2)/unorma1_2)
Wc=abs(unorma2-unorma1_2)

U_zaehl=U_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       U_treffer=U_treffer+1.
end if

end if


;;w-Komponente
if(.not. ismissing(wwind1(k,nyh,i)/unorm1) .and. .not. ismissing(wwind2(k,nyh,i)/unorm2)) then


wnorma2=doubletofloat(wwind2(k,nyh,i)/unorm2)
wnorma1_2=doubletofloat(wwind1(k,nyh,i)/unorm1)

Dc=abs((wnorma2-wnorma1_2)/wnorma1_2)
Wc=abs(wnorma2-wnorma1_2)

W_zaehl=W_zaehl+1

if(Dc .le. D .or. Wc .le. W) then
       W_treffer=W_treffer+1.
end if

end if

end do
end do




print("Trefferquote fuer U bei a2: " +U_treffer/U_zaehl*100.+" %")



print("Trefferquote fuer W bei a2: " +W_treffer/W_zaehl*100.+" %")


if(U_treffer/U_zaehl*100..gt. 95. .and. W_treffer/W_zaehl*100..gt. 95.) then
print("Testfall a2 bestanden")
end if
print(" ")
end if



;---------------------------------------------------------------
;----Advektion, Turbulent, a3-1--------------------------------------------
;---------------------------------------------------------------

if(test_a3_1 .eq. 1) then

gebaeudeh = 25.
wirbell = 0.
l=0

do k=nznb,nznea1_2
do i=nxnb,nxne

;----Bestimmen der x-Koordinaten, bei denen u negativ ist
if(.not. ismissing(uwind1(k,nyh,i))) then 
if(uwind1(k,nyh,i) .lt. 0..and. xu1(i) .ge. 12.5 ) then
l=l+1
xspeicher(l)=xu(i)
ipos(l)=i 
kpos(l)=k
end if
end if

end do
end do



;----das groesste x raussuchen
maxu=0.
do m=1,l-1
if(xspeicher(m) .gt. maxu) then
maxu=doubletofloat(xspeicher(m))
ind_i=ipos(m)
ind_k=kpos(m)
end if
end do



diff=uwind1(ind_k,nyh,ind_i+1)-uwind1(ind_k,nyh,ind_i)




wirbell=maxu-12.5 + dx/(diff/abs(uwind1(ind_k,nyh,ind_i))) ;halbe Gebaeudebreite abziehen


print("Die LÀnge des Nachlaufwirbels betrÀgt bei a3-1: " +wirbell+" m")


if(wirbell .ge. 4.* gebaeudeh .and. wirbell .le. 5. * gebaeudeh) then
print("Testfall a3-1 bestanden")
end if

print(" ")
end if



;---------------------------------------------------------------
;----Advektion, Turbulent, a3-2--------------------------------------------
;---------------------------------------------------------------

if(test_a3_2 .eq. 1) then


wirbell2 = 0.
l=0
ipos=0
kpos=0

do k=nznb,nznea1_2
do i=nxnb,nxne

;----Bestimmen der x-Koordinaten, bei denen u negativ ist
if(.not. ismissing(uwind3(k,nyh,i))) then 
if(uwind3(k,nyh,i) .lt. 0..and. xu(i) .ge. 12.5 ) then
l=l+1
xspeicher(l)=xu(i)
ipos(l)=i 
kpos(l)=k
end if
end if

end do
end do



;----das groesste x raussuchen
maxu=0.
do m=1,l-1
if(xspeicher(m) .gt. maxu) then
maxu=doubletofloat(xspeicher(m))
ind_i=ipos(m)
ind_k=kpos(m)
end if
end do



diff=uwind3(ind_k,nyh,ind_i+1)-uwind3(ind_k,nyh,ind_i)



wirbell2=maxu-12.5+dx/(diff/abs(uwind3(ind_k,nyh,ind_i)))   ;halbe Gebaeudebreite abziehen


print("Die LÀnge des Nachlaufwirbels betrÀgt bei a3-2: " +wirbell2+" m")


if(wirbell2 .gt. wirbell) then
print("Testfall a3-2 bestanden")
end if

print(" ")

end if








end