!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine ogw_stress(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 Const(igit,nb),Ros(igit,nb), * psi(igit,nb), Gsl(igit,nb),alph(igit,nb), * beta(igit,nb),F(igit,nb),N(igit,nb) dimension temp(5) real lon, lat,Const,Ros,psi,Gsl,alph,N,F,beta, * temp,temps c print*,'working on wave stress...' Gsl=1.23 Ros=1.14 do 557 i=1,nb do 558 j=1,igit temp = 0 count1 = 3 count2 = 2 if(i.eq.1) then temp(1)=un1(i,4,j) temp(2)=un1(igit,4,j) temp(3)=un1(i+1,4,j) elseif(i.eq.igit) then temp(1)=un1(i,4,j) temp(2)=un1(i-1,4,j) temp(3)=un1(1,4,j) else temp(1)=un1(i,4,j) temp(2)=un1(i-1,4,j) temp(3)=un1(i+1,4,j) endif if(j.ne.1.and.j.ne.nb) then if(i.eq.1) then temp(4)=un1(i,4,j-1) temp(5)=un1(i,4,j+1) elseif(i.eq.igit) then temp(4)=un1(i,4,j-1) temp(5)=un1(i,4,j+1) else temp(4)=un1(i,4,j-1) temp(5)=un1(i,4,j+1) endif count1 = 5 count2 = 3 endif if(j.eq.1) then if(i.eq.1) then temp(4)=un1(i,4,2) elseif(i.eq.igit) then temp(4)=un1(i,4,2) else temp(4)=un1(i,4,2) endif count1 = 4 count2 = 2 endif if(j.eq.nb) then if(i.eq.1) then temp(4)=un1(i,4,j-1) elseif(i.eq.igit) then temp(4)=un1(i,4,j-1) else temp(4)=un1(i,4,j-1) endif count1 = 4 count2 = 2 endif pass = 1 1001 continue sorted = 1 do 1005 l = 1,count1-pass if(temp(l) .gt. temp(l+1)) then temps = temp(l) temp(l) = temp(l+1) temp(l+1) = temps sorted = 0 endif 1005 continue pass = pass +1 if(sorted .eq. 0) goto 1001 u1_gl(j,i) = temp(count2) 558 continue 557 continue !------------------------read latitude points in a row----------------------------------------- do 4 j=1,nb lat=alatb+(j-ilatb)*igs1 !---------------------------read longitudes to each latitude-------------------------------------- do 5 i=1,igit lon=alonl+(i-ilonl)*igs2 f(i,j)=2*7.2E-05*cos(lat*3.1415926535/180) N(i,j)=sqrt(9.81*(-0.007+9.81/1000)/tn1(j,3,i)) if (vn1(j,4,i).ge.0.and.u1_gl(i,j).gt.0) alph(i,j)= *atan(abs(vn1(j,4,i)/u1_gl(i,j))) if (vn1(j,4,i).gt.0.and.u1_gl(i,j).le.0) alph(i,j)= *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) alph(i,j)= *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) alph(i,j)= *atan(abs(u1_gl(i,j)/vn1(j,4,i)))+3*3.1415926535/2 if (alph(i,j).lt.teta(i,j)) psi(i,j)=teta(i,j)-alph(i,j) if (alph(i,j).ge.teta(i,j).and *.alph(i,j).le.(teta(i,j)+3.1415926535/2)) * psi(i,j)=alph(i,j)-teta(i,j) if (alph(i,j).gt.(teta(i,j)+3.1415926535/2).and *.alph(i,j).lt.(teta(i,j)+3.1415926535)) psi(i,j)= * abs(alph(i,j)-3.1415926535)+teta(i,j) if (alph(i,j).ge.(teta(i,j)+3.1415926535).and *.alph(i,j).le.(teta(i,j)+3*3.1415926535/2)) * psi(i,j)=abs(alph(i,j)-3.1415926535)-teta(i,j) if (alph(i,j).gt.(teta(i,j)+3*3.1415926535/2).and *.alph(i,j).lt.(teta(i,j)+2*3.1415926535)) * psi(i,j)=abs(alph(i,j)-2*3.1415926535)+teta(i,j) if (alph(i,j).ge.(teta(i,j)+2*3.1415926535)) * psi(i,j)=-teta(i,j)+(alph(i,j)-2*3.1415926535) c IF(B(i,j).gt.0) then c write(211,*) i,j,'B=',B(i,j),'C=',C(i,j),'str' c ENDIF DD(i,j)=Gsl(i,j)*(B(i,j)*cos(psi(i,j))**2+C(i,j)*sin(psi(i,j))**2) TT(i,j)=Gsl(i,j)*(B(i,j)-C(i,j))*sin(psi(i,j))*cos(psi(i,j)) beta(i,j)=atan(abs(TT(i,j)/DD(i,j))) Const(i,j)=Ros(i,j)*sqrt(u1_gl(i,j)**2+vn1(j,4,i)**2) **cos(beta(i,j))*N(i,j)*mu(i,j)*sigma(i,j) !--------------computing gravity wave stress/Qm/& its log/ZZZ/----------------------------- Qm(i,j)=Const(i,j)*sqrt(DD(i,j)**2+TT(i,j)**2) !--------------computing horizontal wave number (Kh)-------------------------- if (cos(teta(i,j))*mu(i,j).eq.0.or.abs(u1_gl(i,j)).lt.4.or *.abs(mu(i,j)).lt.1) then Khsqt(i,j)=0 Kh(i,j)=0 else Khsqt(i,j)=(2*Gsl(i,j)*sigma(i,j)*sqrt(DD(i,j)**2+TT(i,j)**2)* *sqrt(u1_gl(i,j)**2+vn1(j,4,i)**2))*cos(beta(i,j))/ *(0.5*cos(teta(i,j))*mu(i,j)) Kh(i,j)=sqrt((Khsqt(i,j)**2+F(i,j)**2)/ *((u1_gl(i,j)**2+vn1(j,4,i)**2)*cos(beta(i,j))*cos(beta(i,j)))) end if !------------------------------------------------------------------------------------------------ ! if(Kh(i,j).gt.1.) then ! Write(211,*) '1=',u1_gl(i,j),'2=',vn1(j,4,i),'3=',beta(i,j),'4=', ! *Khsqt(i,j),'5=',mu(i,j), ! *'6=',sqrt(DD(i,j)**2+TT(i,j)**2),'7=',Kh(i,j) ! endif 5 continue 4 continue ! print*,'stop working on wave stress' call ogw_vert_acceleration(coun,ncom,un1,vn1,tn1,wn1 *,z,UU,Aw,Eps,alpha1) end subroutine ogw_stress