| 1 | subroutine my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh) |
|---|
| 2 | |
|---|
| 3 | *********************************************************************** |
|---|
| 4 | ******* FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) ******* |
|---|
| 5 | ** q2 au interfaces entre mailles. |
|---|
| 6 | *********************************************************************** |
|---|
| 7 | |
|---|
| 8 | |
|---|
| 9 | ************** DECLARATIONS ******************************************* |
|---|
| 10 | |
|---|
| 11 | |
|---|
| 12 | integer imax,kmax |
|---|
| 13 | |
|---|
| 14 | parameter (impmax=5) |
|---|
| 15 | |
|---|
| 16 | real kappa,khmin,kmmin,kqmin,longc |
|---|
| 17 | |
|---|
| 18 | parameter (kappa=0.4) |
|---|
| 19 | parameter (a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08, |
|---|
| 20 | a e1=1.8, e2=1.33) |
|---|
| 21 | parameter (khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5, |
|---|
| 22 | a q2min=0.001, q2lmin=0.001) |
|---|
| 23 | parameter (ghmax=0.023, ghmin=-0.28) |
|---|
| 24 | |
|---|
| 25 | real cd(imax),teta(imax,kmax),u(imax,kmax),v(imax,kmax), |
|---|
| 26 | a z(imax,kmax),zi(imax,kmax+1) |
|---|
| 27 | real kh(imax,kmax+1),km(imax,kmax+1),q2(imax,kmax+1), |
|---|
| 28 | a long(imax,kmax+1) |
|---|
| 29 | real unsdz(imax,kmax),unsdzi(imax,kmax+1) |
|---|
| 30 | real kq(imax,kmax), |
|---|
| 31 | a m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1) |
|---|
| 32 | real a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax), |
|---|
| 33 | a alph(imax,kmax) |
|---|
| 34 | real ksdz2inf,ksdz2sup |
|---|
| 35 | |
|---|
| 36 | c save q2,q2l |
|---|
| 37 | |
|---|
| 38 | save imp |
|---|
| 39 | |
|---|
| 40 | data imp/0/ |
|---|
| 41 | |
|---|
| 42 | *********************************************************************** |
|---|
| 43 | |
|---|
| 44 | imp=imp+1 |
|---|
| 45 | |
|---|
| 46 | ************** INCREMENTS VERTICAUX *********************************** |
|---|
| 47 | |
|---|
| 48 | do 9 i=1,imax |
|---|
| 49 | zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax)) |
|---|
| 50 | 9 continue |
|---|
| 51 | do 10 k=1,kmax |
|---|
| 52 | do 10 i=1,imax |
|---|
| 53 | unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k)) |
|---|
| 54 | 10 continue |
|---|
| 55 | do 11 k=2,kmax |
|---|
| 56 | do 11 i=1,imax |
|---|
| 57 | unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1)) |
|---|
| 58 | 11 continue |
|---|
| 59 | do 12 i=1,imax |
|---|
| 60 | unsdzi(i,1)=0.5/(z(i,1)-zi(i,1)) |
|---|
| 61 | unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax)) |
|---|
| 62 | 12 continue |
|---|
| 63 | |
|---|
| 64 | *********************************************************************** |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | ************** DIFFUSIVITES KH, KM et KQ ****************************** |
|---|
| 68 | * Ci-dessous, une premiere estimation des diffusivites turbulentes km * |
|---|
| 69 | * et kh est effectuee pour utilisation dans les taux de production * |
|---|
| 70 | * et de destruction de q2 et q2l. On calcule aussi kq. * |
|---|
| 71 | |
|---|
| 72 | do 100 k=2,kmax |
|---|
| 73 | do 100 i=1,imax |
|---|
| 74 | beta=2.0/(teta(i,k)+teta(i,k-1)) |
|---|
| 75 | n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1)) |
|---|
| 76 | n2(i,k)=amax1(0.0,n2(i,k)) |
|---|
| 77 | du=unsdzi(i,k)*(u(i,k)-u(i,k-1)) |
|---|
| 78 | dv=unsdzi(i,k)*(v(i,k)-v(i,k-1)) |
|---|
| 79 | m2(i,k)=du*du+dv*dv |
|---|
| 80 | ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10) |
|---|
| 81 | ri(i,k)=amax1(-0.1,min(4.0,ri(i,k))) |
|---|
| 82 | 100 continue |
|---|
| 83 | |
|---|
| 84 | do 110 k=2,kmax |
|---|
| 85 | do 110 i=1,imax |
|---|
| 86 | vt=kappa*(zi(i,k)-zi(i,1)) |
|---|
| 87 | long(i,k)=vt/(1.0+vt/160.0) |
|---|
| 88 | if(n2(i,k).gt.0.0) then |
|---|
| 89 | long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) |
|---|
| 90 | endif |
|---|
| 91 | gh=amax1(ghmin, |
|---|
| 92 | a min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k))) |
|---|
| 93 | sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh* |
|---|
| 94 | a ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/ |
|---|
| 95 | a ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh)) |
|---|
| 96 | km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) |
|---|
| 97 | sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) |
|---|
| 98 | kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) |
|---|
| 99 | 110 continue |
|---|
| 100 | do 111 i=1,imax |
|---|
| 101 | us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1))) |
|---|
| 102 | vt1=(b1**0.666667)*us*us |
|---|
| 103 | vt2=(b1**0.6666667)*kappa*kappa* |
|---|
| 104 | a m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1)) |
|---|
| 105 | c q2(i,1)=amax1(q2min,vt1,vt2) |
|---|
| 106 | q2(i,1)=amax1(q2min,vt1) |
|---|
| 107 | long(i,1)=0.0 |
|---|
| 108 | long(i,kmax+1)=long(i,kmax) |
|---|
| 109 | sq=0.2 |
|---|
| 110 | kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq) |
|---|
| 111 | 111 continue |
|---|
| 112 | |
|---|
| 113 | do 120 k=2,kmax |
|---|
| 114 | do 120 i=1,imax |
|---|
| 115 | longc=0.5*(long(i,k)+long(i,k+1)) |
|---|
| 116 | q2c=0.5*(q2(i,k)+q2(i,k+1)) |
|---|
| 117 | sq=0.2 |
|---|
| 118 | kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq) |
|---|
| 119 | 120 continue |
|---|
| 120 | |
|---|
| 121 | *********************************************************************** |
|---|
| 122 | |
|---|
| 123 | |
|---|
| 124 | ************** CALCUL DE Q2 ******************************************* |
|---|
| 125 | |
|---|
| 126 | do 200 k=2,kmax |
|---|
| 127 | do 200 i=1,imax |
|---|
| 128 | prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k))) |
|---|
| 129 | dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k)) |
|---|
| 130 | a +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k))) |
|---|
| 131 | if(k.lt.kmax) then |
|---|
| 132 | ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k) |
|---|
| 133 | else |
|---|
| 134 | ksdz2sup=0.0 |
|---|
| 135 | endif |
|---|
| 136 | ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1) |
|---|
| 137 | b(i,k)=-ksdz2inf*dt |
|---|
| 138 | a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup) |
|---|
| 139 | c(i,k)=-ksdz2sup*dt |
|---|
| 140 | f(i,k)=q2(i,k)+dt*prod |
|---|
| 141 | 200 continue |
|---|
| 142 | do 201 i=1,imax |
|---|
| 143 | f(i,2)=f(i,2) |
|---|
| 144 | a +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1) |
|---|
| 145 | 201 continue |
|---|
| 146 | |
|---|
| 147 | do 210 i=1,imax |
|---|
| 148 | alph(i,2)=a(i,2) |
|---|
| 149 | 210 continue |
|---|
| 150 | do 211 k=3,kmax |
|---|
| 151 | do 211 i=1,imax |
|---|
| 152 | bet=b(i,k)/alph(i,k-1) |
|---|
| 153 | alph(i,k)=a(i,k)-bet*c(i,k-1) |
|---|
| 154 | f(i,k)=f(i,k)-bet*f(i,k-1) |
|---|
| 155 | 211 continue |
|---|
| 156 | |
|---|
| 157 | do 220 i=1,imax |
|---|
| 158 | q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax)) |
|---|
| 159 | q2(i,kmax+1)=q2(i,kmax) |
|---|
| 160 | 220 continue |
|---|
| 161 | do 221 k=kmax-1,2,-1 |
|---|
| 162 | do 221 i=1,imax |
|---|
| 163 | q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k)) |
|---|
| 164 | 221 continue |
|---|
| 165 | do 222 i=1,imax |
|---|
| 166 | q2(i,2)=amax1(q2(i,2),q2(i,1)) |
|---|
| 167 | 222 continue |
|---|
| 168 | |
|---|
| 169 | *********************************************************************** |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | ************** EVALUATION FINALE DE KH ET KM ************************** |
|---|
| 173 | |
|---|
| 174 | do 400 k=2,kmax |
|---|
| 175 | do 400 i=1,imax |
|---|
| 176 | if(n2(i,k).gt.0.0) then |
|---|
| 177 | long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) |
|---|
| 178 | endif |
|---|
| 179 | gh=amax1(ghmin, |
|---|
| 180 | a min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k))) |
|---|
| 181 | sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh* |
|---|
| 182 | a ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/ |
|---|
| 183 | a ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh)) |
|---|
| 184 | km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) |
|---|
| 185 | sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) |
|---|
| 186 | kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) |
|---|
| 187 | 400 continue |
|---|
| 188 | do 401 i=1,imax |
|---|
| 189 | km(i,1)=kmmin |
|---|
| 190 | km(i,kmax+1)=km(i,kmax) |
|---|
| 191 | kh(i,1)=khmin |
|---|
| 192 | kh(i,kmax+1)=kh(i,kmax) |
|---|
| 193 | 401 continue |
|---|
| 194 | |
|---|
| 195 | *********************************************************************** |
|---|
| 196 | |
|---|
| 197 | c if(imp.eq.impmax) then |
|---|
| 198 | am=1.0/float(imax) |
|---|
| 199 | imp=0 |
|---|
| 200 | do 1000 k=kmax,1,-1 |
|---|
| 201 | au=0.0 |
|---|
| 202 | ateta=0.0 |
|---|
| 203 | az=0.0 |
|---|
| 204 | adz=0.0 |
|---|
| 205 | akq=0.0 |
|---|
| 206 | acd=0.0 |
|---|
| 207 | do 1001 i=1,imax |
|---|
| 208 | au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k)) |
|---|
| 209 | ateta=ateta+am*teta(i,k) |
|---|
| 210 | az=az+am*z(i,k) |
|---|
| 211 | adz=adz+am*(1.0/unsdz(i,k)) |
|---|
| 212 | akq=akq+am*kq(i,k) |
|---|
| 213 | acd=acd+am*cd(i) |
|---|
| 214 | 1001 continue |
|---|
| 215 | c write(*,2000) k,az,adz,au,ateta,akq,acd*1000.0 |
|---|
| 216 | 2000 format(2x,i3,2x,6(f9.2,2x)) |
|---|
| 217 | 1000 continue |
|---|
| 218 | |
|---|
| 219 | write(*,*) |
|---|
| 220 | write(*,*) |
|---|
| 221 | |
|---|
| 222 | do 1002 k=kmax+1,1,-1 |
|---|
| 223 | azi=0.0 |
|---|
| 224 | adzi=0.0 |
|---|
| 225 | aq2=0.0 |
|---|
| 226 | al=0.0 |
|---|
| 227 | akm=0.0 |
|---|
| 228 | akh=0.0 |
|---|
| 229 | am2=0.0 |
|---|
| 230 | al0=0.0 |
|---|
| 231 | do 1003 i=1,imax |
|---|
| 232 | azi=azi+am*zi(i,k) |
|---|
| 233 | adzi=adzi+am*(1.0/unsdzi(i,k)) |
|---|
| 234 | aq2=aq2+am*q2(i,k) |
|---|
| 235 | al=al+am*long(i,k) |
|---|
| 236 | akm=akm+am*km(i,k) |
|---|
| 237 | akh=akh+am*kh(i,k) |
|---|
| 238 | am2=am2+am*m2(i,k) |
|---|
| 239 | c al0=al0+am*long0d(i) |
|---|
| 240 | 1003 continue |
|---|
| 241 | c write(*,2001) k,azi,aq2,al,akm,akh,am2*1.0e5 |
|---|
| 242 | 2001 format(2x,i3,6(2x,f9.3)) |
|---|
| 243 | 1002 continue |
|---|
| 244 | c endif |
|---|
| 245 | |
|---|
| 246 | return |
|---|
| 247 | end |
|---|
