Changeset 4203 for dynamico_lmdz/simple_physics/phyparam
- Timestamp:
- Dec 20, 2019, 3:25:27 PM (5 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/physics/mellor_yamada.F90
r4202 r4203 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 1 MODULE mellor_yamada 2 3 IMPLICIT NONE 4 5 CONTAINS 6 7 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 15 INTEGER, INTENT(IN) :: imax,kmax 16 REAL, INTENT(IN) :: dt, gp 17 REAL, INTENT(IN) :: z(imax,kmax), & 18 u(imax,kmax), v(imax,kmax), teta(imax,kmax), cd(imax) 19 REAL, INTENT(INOUT) :: zi(imax,kmax+1), q2(imax,kmax+1) 20 REAL, INTENT(OUT) :: long(imax,kmax+1), km(imax,kmax+1), kh(imax,kmax+1) 21 22 !************* DECLARATIONS ******************************************* 23 24 INTEGER, PARAMETER :: impmax=5 25 REAL, PARAMETER :: kappa=0.4, & 26 a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08, & 27 e1=1.8, e2=1.33, & 28 khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5, & 29 q2min=0.001, q2lmin=0.001, & 30 ghmax=0.023, ghmin=-0.28 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
Note: See TracChangeset
for help on using the changeset viewer.
