subroutine ogw_vert_acceleration(coun,ncom,un1,vn1,tn1,wn1 *,z,UU,Aw,Eps,alpha1) include 'com_topo.fc' c parameters real un1(nb,kgit,igit),tn1(nb,kgit,igit),vn1(nb,kgit,igit) real UU(igit,nb,kgit),Aw(igit,nb,kgit),z(kgit), * alpha1(igit,nb,kgit),Eps(igit,nb,kgit),wn1(nb,kgit+1,igit) integer ncom c other dimension Ros(igit,nb),Fver(igit,nb),Wav(nb,kgit),W_av(nb,45), * lat(nb),lon(igit),alt(kgit),zstep(kgit) dimension bver(igit,nb,kgit),beta(igit,nb,kgit), * Rovert(igit,nb,kgit),aver(igit,nb,kgit),mver(igit,nb,kgit), * Nver(igit,nb,kgit),nuver(igit,nb,kgit) real lon, lat,ZhINT,Ros,Rovert,Nver,Fver,beta,Wav real bver,nuver,aver, zstep,mver Ros = 1.27 c print*,'working on vertical steps...', coun do m=2,kgit-1 if(m.eq.2)then alt(m)=0. zstep(m) = int(z(m-1)*1000.) else alt(m)=int(z(m-2)*1000.) zstep(m)=int(z(m-1)*1000.) - int(z(m-2)*1000.) endif !write(211,*) m,alt(m),zstep(m) enddo c print*, 'MU,QM:',mu(1,1),Qm(1,1) do 44 j=1,nb do 44 i=1,igit do m=2,kgit-1 if(m.eq.2)then Rovert(i,j,m) = 1.27 else Rovert(i,j,m) = Rovert(i,j,m-1) * *exp(-(1/7+abs(tn1(j,m,i)-tn1(j,m-1,i)) * /(zstep(m)*tn1(j,m-1,i)))*zstep(m)) endif enddo 44 continue ! read latitude points in a row do 4 j=1,nb lat(j)=alatb+(j-ilatb)*igs1 ! read longitudes to each latitude do 5 i=1,igit Fver(i,j)=2*7.2E-05*cos(lat(j)*3.1415926535/180.) lon(i)=alonl+(i-ilonl)*igs2 ! print*,'1' if (vn1(j,4,i).ge.0.and.u1_gl(i,j).gt.0) alpha1(i,j,1)= * atan(abs(vn1(j,4,i)/u1_gl(i,j))) if (vn1(j,4,i).gt.0.and.u1_gl(i,j).le.0) alpha1(i,j,1)= * atan(abs(u1_gl(i,j)/vn1(j,4,i)))+3.1415926535/2 if (vn1(j,4,i).le.0.and.u1_gl(i,j).lt.0) alpha1(i,j,1)= * atan(abs(vn1(j,4,i)/u1_gl(i,j)))+3.1415926535 if (vn1(j,4,i).lt.0.and.u1_gl(i,j).ge.0) alpha1(i,j,1)= * atan(abs(u1_gl(i,j)/vn1(j,4,i)))+3*3.1415926535/2 beta(i,j,1)=atan(TT(i,j)/DD(i,j))+3.1415926535 IF(Khsqt(i,j).eq.0.) THEN do m=2,kgit-1 UU(i,j,m-1)=0; aver(i,j,m-1)=0; mver(i,j,m-1)=0 Eps(i,j,m-1)=0; Aw(i,j,m-1)=0 enddo ELSE do 6 m=2,kgit-1 if(m.lt.4)then UU(i,j,m-1)=0 Eps(i,j,m-1)=0 Aw(i,j,m-1)=0 aver(i,j,m-1)=0 Nver(i,j,m-1)=0 nuver(i,j,m-1)=0 elseif(m.eq.4)then Nver(i,j,m-1)=sqrt(9.81*(abs(tn1(j,m,i)-tn1(j,m-1,i)) * /zstep(m)+9.81/1000.)/tn1(j,m-1,i)) aver(i,j,m-1)=(Rovert(i,j,m)*Fver(i,j)**2*Khsqt(i,j))/ *(2*Kh(i,j)**2*Nver(i,j,m-1)* *sqrt(un1(j,m-1,i)**2+vn1(j,m-1,i)**2)*abs(cos(beta(i,j,1)))) UU(i,j,m-1)=0 nuver(i,j,m-1)=0 Eps(i,j,m-1)=0 Aw(i,j,m-1)=0 else nuver(i,j,m-1)=0.01*exp(abs(0.006*log(100.) **alt(m)/1000.)) !====================computing vertical gridsteps usig============================================================================ if (vn1(j,m-1,i).ge.0.and.un1(j,m-1,i).gt.0) alpha1(i,j,m-1)= * atan(abs(vn1(j,m-1,i)/un1(j,m-1,i))) if (vn1(j,m-1,i).gt.0.and.un1(j,m-1,i).le.0) alpha1(i,j,m-1)= * atan(abs(un1(j,m-1,i)/vn1(j,m-1,i)))+3.1415926535/2 if (vn1(j,m-1,i).le.0.and.un1(j,m-1,i).lt.0) alpha1(i,j,m-1)= * atan(abs(vn1(j,m-1,i)/un1(j,m-1,i)))+3.1415926535 if (vn1(j,m-1,i).lt.0.and.un1(j,m-1,i).ge.0) alpha1(i,j,m-1)= * atan(abs(un1(j,m-1,i)/vn1(j,m-1,i)))+3*3.1415926535/2 if (alpha1(i,j,1).lt.teta(i,j)) beta(i,j,m-1)= * abs((3.1415926535-atan(abs(TT(i,j)/DD(i,j)))+ * alpha1(i,j,1))-alpha1(i,j,m-1)) if (alpha1(i,j,1).ge.teta(i,j).and * .alpha1(i,j,1).le.(teta(i,j)+3.1415926535/2)) beta(i,j,m-1)= * abs((3.1415926535+atan(abs(TT(i,j)/DD(i,j))) * +alpha1(i,j,1))-alpha1(i,j,m-1)) if (alpha1(i,j,1).gt.(teta(i,j)+3.1415926535/2).and * .alpha1(i,j,1).lt.(teta(i,j)+3.1415926535)) beta(i,j,m-1)= * abs((3.1415926535-atan(abs(TT(i,j)/DD(i,j)))+alpha1(i,j,1)) * -alpha1(i,j,m-1)) if (alpha1(i,j,1).ge.(teta(i,j)+3.1415926535).and * .alpha1(i,j,1).le.(teta(i,j)+3*3.1415926535/2)) beta(i,j,m-1)= * abs((3.1415926535+atan(abs(TT(i,j)/DD(i,j)))+alpha1(i,j,1)) * -alpha1(i,j,m-1)) if (alpha1(i,j,1).gt.(teta(i,j)+3*3.1415926535/2).and * .alpha1(i,j,1).lt.(teta(i,j)+2*3.1415926535)) beta(i,j,m-1)= * abs((3.1415926535-atan(abs(TT(i,j)/DD(i,j)))+alpha1(i,j,1)) * -alpha1(i,j,m-1)) if (alpha1(i,j,1).ge.(teta(i,j)+2*3.1415926535)) beta(i,j,m-1)= *abs((3.1415926535-atan(abs(TT(i,j)/DD(i,j)))-alpha1(i,j,1)) * -alpha1(i,j,m-1)) if (Kh(i,j)**2*(un1(j,m-2,i)**2+vn1(j,m-2,i)**2)* * (cos(beta(i,j,m-2)))**2.le.Fver(i,j)**2.or * .Kh(i,j)**2*(un1(j,m-1,i)**2+vn1(j,m-1,i)**2)* * (cos(beta(i,j,m-1)))**2.le.Fver(i,j)**2) then do l=m,kgit-1 UU(i,j,l-1)=0; aver(i,j,l-1)=0; mver(i,j,l-1)=0; bver(i,j,l-2)=0 Eps(i,j,l-1)=0; Aw(i,j,l-1)=0 enddo goto 5 else if(abs(un1(j,m-2,i)).lt.4) then do l=m,kgit-1 UU(i,j,l-1)=0; aver(i,j,l-1)=0; mver(i,j,l-1)=0; bver(i,j,l-2)=0 Eps(i,j,l-1)=0; Aw(i,j,l-1)=0 enddo goto 5 endif mver(i,j,m-1)=Kh(i,j)*Nver(i,j,m-2)/sqrt(Kh(i,j)**2* *(un1(j,m-2,i)**2+vn1(j,m-2,i)**2)* *(cos(beta(i,j,m-2)))**2-Fver(i,j)**2) ! write(211,*) mver(i,j,m-1),un1(j,m-2,i),m ! if(mver(i,j,m-1).gt.0.03) then ! mver(i,j,m-1) = 0.03 ! endif ! if(mver(i,j,m-1).lt.0.004) then ! mver(i,j,m-1) = 0.004 ! endif Nver(i,j,m-1)=sqrt(9.81*(abs(tn1(j,m,i)-tn1(j,m-1,i)) * /zstep(m)+9.81/1000.)/tn1(j,m-1,i)) aver(i,j,m-1)=(Rovert(i,j,m)*Fver(i,j)**2*sqrt(Kh(i,j)**2* * (un1(j,m-1,i)**2+vn1(j,m-1,i)**2)* * (cos(beta(i,j,m-1)))**2-Fver(i,j)**2))/ *(2*Kh(i,j)**2*Nver(i,j,m-1)*sqrt(un1(j,m-1,i)**2+vn1(j,m-1,i)**2)) bver(i,j,m-2)=Rovert(i,j,m-1)*mver(i,j,m-1)**2*nuver(i,j,m-1) if(m.eq.5) then UU(i,j,m-1)=2.*Nver(i,j,m-1)*Qm(i,j)/Khsqt(i,j) else UU(i,j,m-1)=aver(i,j,m-2)/aver(i,j,m-1)*UU(i,j,m-2)* * exp(-bver(i,j,m-2)/aver(i,j,m-2)*zstep(m)) endif Eps(i,j,m-1)=mver(i,j,m-1)**2*nuver(i,j,m-1)*UU(i,j,m-1)/1000 if(Eps(i,j,m-1).gt.0.0005) then Eps(i,j,m-1) = 0.0005 endif if(sqrt(un1(j,m-1,i)**2+vn1(j,m-1,i)**2)*cos(beta(i,j,m-1)).gt.0.) * then Aw(i,j,m-1)=(-1000)*Eps(i,j,m-1)*0.5*(1+1/(1.4-1))/ * sqrt(un1(j,m-1,i)**2+vn1(j,m-1,i)**2) else Aw(i,j,m-1)=1000*Eps(i,j,m-1)*0.5*(1+1/(1.4-1))/ * sqrt(un1(j,m-1,i)**2+vn1(j,m-1,i)**2) endif endif !================end computing vertical gridsteps usig MSISe============================================================================ endif !200 format(5f18.10) 6 enddo ENDIF 5 enddo 4 enddo !=============output UU================================================================================================================ do m=2,kgit do j=1,nb www=0;ccc=0 do i=1,igit www = www + wn1(j,m,i) ccc = ccc +1 ! IF(sqrt(UU(i,j,m)).gt.0.and.m.lt.10) then ! write(211,*) i,j,m,sqrt(UU(i,j,m)),ncom,'<---' ! ENDIF enddo Wav(j,m)=www/ccc ! write(211,*) Wav(j,m),j,m,alt(m) enddo enddo do j=1,nb do i=1,igit UU(i,j,4) = 0. Aw(i,j,4) = 0. Eps(i,j,4) = 0. enddo enddo ! do j=1,nb ! numb = 0 ! do 766 m=2,45 ! if(coun.eq.0)then ! write(211,*) i,j,m,so3(j,m),co3(m,j) ! endif ! ZhINT=FLOAT(m-1)*0.25*7*1000. ! write(211,*) alt(m),ZhINT ! do k=2,kgit ! if(ZhINT.ge.alt(k-1).and.ZhINT.lt.alt(k)) then ! numb=numb+1 ! W_av(j,numb) = Wav(j,k-1)+Wav(j,k)-Wav(j,k-1)*(ZhINT-alt(k-1))/ ! *(alt(k)-alt(k-1)) ! write(211,*) j,ZhINT,W_av(j,numb), ! *CO3(numb,j),W_av(j,numb)*CO3(numb,j) ! goto 766 ! endif ! enddo ! 766 continue ! enddo !================================================================================================== !====================================================================================================================================== coun = coun + 1 end subroutine ogw_vert_acceleration