; 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