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