| 1 | ! |
|---|
| 2 | ! $Id$ |
|---|
| 3 | ! |
|---|
| 4 | SUBROUTINE clc_core_cp(du, dt, dq, t1, q1, & |
|---|
| 5 | zu, zt, zq, & |
|---|
| 6 | P,& |
|---|
| 7 | n_it, mixte, & |
|---|
| 8 | coeffs, rugosm, rugosh) |
|---|
| 9 | |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | |
|---|
| 12 | |
|---|
| 13 | ! INTEGER :: i,j,k |
|---|
| 14 | REAL, INTENT(IN) :: du,dt,dq,t1,q1 |
|---|
| 15 | REAL, INTENT(IN) :: zu,zt,zq, P |
|---|
| 16 | INTEGER, INTENT(IN):: n_it |
|---|
| 17 | LOGICAL, INTENT(IN) :: mixte |
|---|
| 18 | REAL, DIMENSION(3), INTENT(OUT) :: coeffs |
|---|
| 19 | real, intent(out) :: rugosm |
|---|
| 20 | real, intent(out) :: rugosh |
|---|
| 21 | |
|---|
| 22 | |
|---|
| 23 | ! REAL :: du,dt,t1,dq,q1 |
|---|
| 24 | REAL :: dt_u, dq_u |
|---|
| 25 | |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | REAL :: cd, ch, cd_rt, cd_n10, ce_n10, ce,& |
|---|
| 29 | u_N, X,& |
|---|
| 30 | rho,cpa,le |
|---|
| 31 | |
|---|
| 32 | REAL :: usr,tsr,qsr |
|---|
| 33 | REAL, DIMENSION(3) :: zeta ! dans l'ordre: it courant, it precedente, var |
|---|
| 34 | REAL, DIMENSION(2,3) :: psi |
|---|
| 35 | REAL, PARAMETER :: g = 9.81, k=.41, sqrt_k = .640312 !avec g= gravité, k = cste von karman, sqrt_k = racine(k) |
|---|
| 36 | REAL :: logzu_10,logzt_10,logzq_10,logzt_zu,logzq_zu,& |
|---|
| 37 | u,cd_n10_rt,ch_n10_rt,ch_n10,chi,z_0m,z_0h,tv |
|---|
| 38 | INTEGER :: i,j,m1,m2,m3,m4,m5 |
|---|
| 39 | |
|---|
| 40 | |
|---|
| 41 | REAL :: phih,phid,phie |
|---|
| 42 | |
|---|
| 43 | |
|---|
| 44 | REAL, parameter :: alpha=2.7e-3,& !alpha = a1 ds formulation smith |
|---|
| 45 | beta = 1.42e-4,& !beta = a2 ds la formulation smith |
|---|
| 46 | gamma = 7.64e-5,& !Large and yaeger 2006 !gamma = a3 ds formulation smith |
|---|
| 47 | ! gamma = 13.09e-3,& !Large and yaeger 2004 |
|---|
| 48 | q0 = 1.64474 ! q0 dzfinis mais pas reutilisée ds la suite ?? |
|---|
| 49 | |
|---|
| 50 | REAL, dimension(3) :: z |
|---|
| 51 | integer, dimension(2) :: tmp_shape |
|---|
| 52 | z = [zu, zt, zq] |
|---|
| 53 | |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | logzu_10 = LOG(zu/10.) |
|---|
| 57 | logzt_10 = LOG(zt/10.) |
|---|
| 58 | logzq_10 = LOG(zq/10.) |
|---|
| 59 | logzt_zu = logzt_10 - logzu_10 |
|---|
| 60 | logzq_zu = logzq_10 - logzu_10 |
|---|
| 61 | |
|---|
| 62 | cpa = 1004.67 !!! pas reutilise ds nvelle formulation... |
|---|
| 63 | |
|---|
| 64 | tv = t1 * (1+.608*q1) !!!!! tv pas reutiliser ds nvelle formulation... |
|---|
| 65 | rho = p / (287.1 * t1 * (1.+.61*q1)) !!!! rho pas reutiliser ds nvelle formulation.... |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | le = (2.501-.00237*(t1-273.15-dt))*1.e6 !!!!!! le pas reutiliser ds nvelle formulation.... |
|---|
| 69 | |
|---|
| 70 | u = max(du,.5) !!!! u = vitesse du vent 1 er niveau |
|---|
| 71 | u_N = u !!!!!!! u_N utilisée pour le calcul des premier calcul... |
|---|
| 72 | |
|---|
| 73 | |
|---|
| 74 | !!!!! initialisation des parametres et premier calcul de z_0 cd et ch |
|---|
| 75 | |
|---|
| 76 | |
|---|
| 77 | if (mixte) then |
|---|
| 78 | z_0m = 1.e-4 !!!!! pour l'initialisation on fixe z_0 a 10-4 |
|---|
| 79 | cd_n10 = k**2 / (log(10/z_0m))**2 !!!!! on calcul ensuite cd a 10 m d'après smith |
|---|
| 80 | cd_n10_rt = sqrt(cd_n10) |
|---|
| 81 | else |
|---|
| 82 | ! Large and yaeger 2006 |
|---|
| 83 | cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur |
|---|
| 84 | ! Large and yaeger 2004 |
|---|
| 85 | ! cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) |
|---|
| 86 | cd_n10_rt = sqrt(cd_n10) !!!!!! |
|---|
| 87 | z_0m = 10. * exp(-k/cd_n10_rt) !!!!!!! on calcul le premier z_0 |
|---|
| 88 | endif |
|---|
| 89 | |
|---|
| 90 | |
|---|
| 91 | if ( dt .gt. 0. ) then |
|---|
| 92 | ch_n10 = 0.018 * cd_n10_rt !!!!!!! si regime stable calcul du ch |
|---|
| 93 | ! z_0h = 10. * exp(-k/ch_n10) |
|---|
| 94 | else |
|---|
| 95 | ch_n10 = 0.0327 * cd_n10_rt !!!!!!! si regime instable calcul du ch |
|---|
| 96 | ! z_0h = 10. * exp(-k/ch_n10) |
|---|
| 97 | endif |
|---|
| 98 | |
|---|
| 99 | ce_n10 = 3.46e-2 *cd_n10_rt !!!!!! calcul ensuite de ce |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | dt_u = dt !!! def |
|---|
| 103 | dq_u = dq !!! def |
|---|
| 104 | |
|---|
| 105 | cd = cd_n10 !!! def |
|---|
| 106 | ch = ch_n10 !!! def |
|---|
| 107 | ce = ce_n10 !!! def |
|---|
| 108 | |
|---|
| 109 | cd_rt = sqrt(cd) !!! def |
|---|
| 110 | |
|---|
| 111 | ! open(7,file="err_rel.dat",position="append") |
|---|
| 112 | |
|---|
| 113 | |
|---|
| 114 | !!!!!!!!! début des itérations.... |
|---|
| 115 | |
|---|
| 116 | do i = 1,n_it |
|---|
| 117 | |
|---|
| 118 | usr = cd_rt * u |
|---|
| 119 | tsr = ch / cd_rt * dt_u |
|---|
| 120 | qsr = ce / cd_rt * dq_u |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | zeta = k*g / (usr**2) * tsr / t1& |
|---|
| 124 | * z |
|---|
| 125 | |
|---|
| 126 | zeta = sign( min(abs(zeta),10.0), zeta ) |
|---|
| 127 | |
|---|
| 128 | |
|---|
| 129 | !!!! calcul du zeta pour connaitre le signe et donc en deduire le regime... |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | IF ( zeta(1) .GT. 0 ) THEN !!! si regime stable alors valeurs pr chi et psi |
|---|
| 135 | |
|---|
| 136 | chi = 0.018 |
|---|
| 137 | |
|---|
| 138 | |
|---|
| 139 | psi(1,1) = -5.*zeta(1) |
|---|
| 140 | psi(2,1) = psi(1,1) |
|---|
| 141 | |
|---|
| 142 | ELSE !!!!! si regime instable alors valeurs pr psi et chi !!! |
|---|
| 143 | |
|---|
| 144 | chi = 0.0327 |
|---|
| 145 | |
|---|
| 146 | X = sqrt( sqrt( 1-16.*zeta(1) )) |
|---|
| 147 | psi(1,1) = 3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X) |
|---|
| 148 | psi(2,1) = 2.*LOG( (1 + X**2) / 2 ) |
|---|
| 149 | |
|---|
| 150 | END IF |
|---|
| 151 | |
|---|
| 152 | do j =2,3 |
|---|
| 153 | IF ( zeta(j) .GT. 0 ) THEN |
|---|
| 154 | psi(1,j) = -5.*zeta(j) |
|---|
| 155 | psi(2,j) = psi(1,1) |
|---|
| 156 | |
|---|
| 157 | ELSE |
|---|
| 158 | X = sqrt( sqrt( 1-16.*zeta(j) )) |
|---|
| 159 | psi(1,j) = 3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X) |
|---|
| 160 | psi(2,j) = 2.*LOG( (1 + X**2) / 2 ) |
|---|
| 161 | |
|---|
| 162 | END IF |
|---|
| 163 | enddo |
|---|
| 164 | |
|---|
| 165 | |
|---|
| 166 | dt_u = dt - tsr / k * & |
|---|
| 167 | ( logzt_zu + psi(2,1) - psi(2,2) ) |
|---|
| 168 | |
|---|
| 169 | dq_u = dq - qsr / k * & |
|---|
| 170 | ( logzq_zu + psi(2,1) - psi(2,3) ) |
|---|
| 171 | |
|---|
| 172 | |
|---|
| 173 | |
|---|
| 174 | u_N = u / ( 1 + cd_n10_rt / k * ( logzu_10 - psi(1,1)) ) |
|---|
| 175 | |
|---|
| 176 | if (mixte) then |
|---|
| 177 | z_0m=0.018*usr**2/9.81 + 0.11*14e-6 / (usr) |
|---|
| 178 | ! z_0h=(0.11*14e-6 / (usr)) + 1.4e-5 |
|---|
| 179 | cd_n10 = k**2 / (log(10/z_0m))**2 |
|---|
| 180 | cd_n10_rt = sqrt(cd_n10) |
|---|
| 181 | ! ch_n10 = k**2 / ((log(10/z_0m))*(log(10/z_0h))) |
|---|
| 182 | ! ch_n10_rt = sqrt(ch_n10) |
|---|
| 183 | ch_n10 = 1.e-3 |
|---|
| 184 | else |
|---|
| 185 | ! Large and yaeger 2006 |
|---|
| 186 | cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur |
|---|
| 187 | ! Large and yaeger 2004 |
|---|
| 188 | ! cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) |
|---|
| 189 | ! cd_n10_dun = -alpha / u_n**2 + gamma |
|---|
| 190 | cd_n10_rt = sqrt(cd_n10) |
|---|
| 191 | z_0m = 10. * exp(-k/cd_n10_rt) !!!!! c'est sur cette ligne la qu'il y a un bug... ou sur la ligne du u_N |
|---|
| 192 | ! ch_n10 = chi * cd_n10_rt |
|---|
| 193 | ! z_0h = 10. * exp(-k/ch_n10) |
|---|
| 194 | ! ch_n10_rt = sqrt(ch_n10) |
|---|
| 195 | ch_n10 = chi * cd_n10_rt |
|---|
| 196 | endif |
|---|
| 197 | |
|---|
| 198 | |
|---|
| 199 | ce_n10 = 3.46e-2*cd_n10_rt |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | |
|---|
| 203 | phid = 1+cd_n10_rt/k * (logzu_10 - psi(1,1)) |
|---|
| 204 | phih = 1+chi/k * (logzu_10 - psi(2,1)) |
|---|
| 205 | phie = 1+3.46e-2/k * (logzu_10 - psi(2,1)) |
|---|
| 206 | |
|---|
| 207 | |
|---|
| 208 | cd_rt = cd_n10_rt / abs( phid ) |
|---|
| 209 | ch = ch_n10 / ( phih * phid ) |
|---|
| 210 | ce = ce_n10 / ( phie * phid ) |
|---|
| 211 | |
|---|
| 212 | cd=cd_rt**2 |
|---|
| 213 | |
|---|
| 214 | END DO |
|---|
| 215 | |
|---|
| 216 | usr = cd_rt * u |
|---|
| 217 | tsr = ch / cd_rt * dt_u |
|---|
| 218 | qsr = ce / cd_rt * dq_u |
|---|
| 219 | |
|---|
| 220 | |
|---|
| 221 | ! coeffs(:,m) = [ & |
|---|
| 222 | ! rho * usr**2,& |
|---|
| 223 | ! rho * cpa * usr *tsr,& |
|---|
| 224 | ! rho * le * usr *qsr] |
|---|
| 225 | |
|---|
| 226 | coeffs = [ cd_rt**2, ch , ce ] |
|---|
| 227 | |
|---|
| 228 | rugosm = z_0m |
|---|
| 229 | rugosh = z_0h |
|---|
| 230 | |
|---|
| 231 | RETURN |
|---|
| 232 | END subroutine clc_core_cp |
|---|
