subroutine my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh) *********************************************************************** ******* FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) ******* ** q2 au interfaces entre mailles. *********************************************************************** ************** DECLARATIONS ******************************************* integer imax,kmax parameter (impmax=5) real kappa,khmin,kmmin,kqmin,longc parameter (kappa=0.4) parameter (a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08, a e1=1.8, e2=1.33) parameter (khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5, a q2min=0.001, q2lmin=0.001) parameter (ghmax=0.023, ghmin=-0.28) real cd(imax),teta(imax,kmax),u(imax,kmax),v(imax,kmax), a z(imax,kmax),zi(imax,kmax+1) real kh(imax,kmax+1),km(imax,kmax+1),q2(imax,kmax+1), a long(imax,kmax+1) real unsdz(imax,kmax),unsdzi(imax,kmax+1) real kq(imax,kmax), a m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1) real a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax), a alph(imax,kmax) real ksdz2inf,ksdz2sup c save q2,q2l save imp data imp/0/ *********************************************************************** imp=imp+1 ************** INCREMENTS VERTICAUX *********************************** do 9 i=1,imax zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax)) 9 continue do 10 k=1,kmax do 10 i=1,imax unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k)) 10 continue do 11 k=2,kmax do 11 i=1,imax unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1)) 11 continue do 12 i=1,imax unsdzi(i,1)=0.5/(z(i,1)-zi(i,1)) unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax)) 12 continue *********************************************************************** ************** DIFFUSIVITES KH, KM et KQ ****************************** * Ci-dessous, une premiere estimation des diffusivites turbulentes km * * et kh est effectuee pour utilisation dans les taux de production * * et de destruction de q2 et q2l. On calcule aussi kq. * do 100 k=2,kmax do 100 i=1,imax beta=2.0/(teta(i,k)+teta(i,k-1)) n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1)) n2(i,k)=amax1(0.0,n2(i,k)) du=unsdzi(i,k)*(u(i,k)-u(i,k-1)) dv=unsdzi(i,k)*(v(i,k)-v(i,k-1)) m2(i,k)=du*du+dv*dv ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10) ri(i,k)=amax1(-0.1,min(4.0,ri(i,k))) 100 continue do 110 k=2,kmax do 110 i=1,imax vt=kappa*(zi(i,k)-zi(i,1)) long(i,k)=vt/(1.0+vt/160.0) if(n2(i,k).gt.0.0) then long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) endif gh=amax1(ghmin, a min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k))) sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh* a ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/ a ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh)) km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) 110 continue do 111 i=1,imax us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1))) vt1=(b1**0.666667)*us*us vt2=(b1**0.6666667)*kappa*kappa* a m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1)) c q2(i,1)=amax1(q2min,vt1,vt2) q2(i,1)=amax1(q2min,vt1) long(i,1)=0.0 long(i,kmax+1)=long(i,kmax) sq=0.2 kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq) 111 continue do 120 k=2,kmax do 120 i=1,imax longc=0.5*(long(i,k)+long(i,k+1)) q2c=0.5*(q2(i,k)+q2(i,k+1)) sq=0.2 kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq) 120 continue *********************************************************************** ************** CALCUL DE Q2 ******************************************* do 200 k=2,kmax do 200 i=1,imax prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k))) dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k)) a +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k))) if(k.lt.kmax) then ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k) else ksdz2sup=0.0 endif ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1) b(i,k)=-ksdz2inf*dt a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup) c(i,k)=-ksdz2sup*dt f(i,k)=q2(i,k)+dt*prod 200 continue do 201 i=1,imax f(i,2)=f(i,2) a +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1) 201 continue do 210 i=1,imax alph(i,2)=a(i,2) 210 continue do 211 k=3,kmax do 211 i=1,imax bet=b(i,k)/alph(i,k-1) alph(i,k)=a(i,k)-bet*c(i,k-1) f(i,k)=f(i,k)-bet*f(i,k-1) 211 continue do 220 i=1,imax q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax)) q2(i,kmax+1)=q2(i,kmax) 220 continue do 221 k=kmax-1,2,-1 do 221 i=1,imax q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k)) 221 continue do 222 i=1,imax q2(i,2)=amax1(q2(i,2),q2(i,1)) 222 continue *********************************************************************** ************** EVALUATION FINALE DE KH ET KM ************************** do 400 k=2,kmax do 400 i=1,imax if(n2(i,k).gt.0.0) then long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) endif gh=amax1(ghmin, a min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k))) sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh* a ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/ a ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh)) km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) 400 continue do 401 i=1,imax km(i,1)=kmmin km(i,kmax+1)=km(i,kmax) kh(i,1)=khmin kh(i,kmax+1)=kh(i,kmax) 401 continue *********************************************************************** c if(imp.eq.impmax) then am=1.0/float(imax) imp=0 do 1000 k=kmax,1,-1 au=0.0 ateta=0.0 az=0.0 adz=0.0 akq=0.0 acd=0.0 do 1001 i=1,imax au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k)) ateta=ateta+am*teta(i,k) az=az+am*z(i,k) adz=adz+am*(1.0/unsdz(i,k)) akq=akq+am*kq(i,k) acd=acd+am*cd(i) 1001 continue c write(*,2000) k,az,adz,au,ateta,akq,acd*1000.0 2000 format(2x,i3,2x,6(f9.2,2x)) 1000 continue write(*,*) write(*,*) do 1002 k=kmax+1,1,-1 azi=0.0 adzi=0.0 aq2=0.0 al=0.0 akm=0.0 akh=0.0 am2=0.0 al0=0.0 do 1003 i=1,imax azi=azi+am*zi(i,k) adzi=adzi+am*(1.0/unsdzi(i,k)) aq2=aq2+am*q2(i,k) al=al+am*long(i,k) akm=akm+am*km(i,k) akh=akh+am*kh(i,k) am2=am2+am*m2(i,k) c al0=al0+am*long0d(i) 1003 continue c write(*,2001) k,azi,aq2,al,akm,akh,am2*1.0e5 2001 format(2x,i3,6(2x,f9.3)) 1002 continue c endif return end