! ! $Id$ ! SUBROUTINE clc_core_cp(du, dt, dq, t1, q1, & zu, zt, zq, & P,& n_it, mixte, & coeffs, rugosm, rugosh) IMPLICIT NONE ! INTEGER :: i,j,k REAL, INTENT(IN) :: du,dt,dq,t1,q1 REAL, INTENT(IN) :: zu,zt,zq, P INTEGER, INTENT(IN):: n_it LOGICAL, INTENT(IN) :: mixte REAL, DIMENSION(3), INTENT(OUT) :: coeffs real, intent(out) :: rugosm real, intent(out) :: rugosh ! REAL :: du,dt,t1,dq,q1 REAL :: dt_u, dq_u REAL :: cd, ch, cd_rt, cd_n10, ce_n10, ce,& u_N, X,& rho,cpa,le REAL :: usr,tsr,qsr REAL, DIMENSION(3) :: zeta ! dans l'ordre: it courant, it precedente, var REAL, DIMENSION(2,3) :: psi REAL, PARAMETER :: g = 9.81, k=.41, sqrt_k = .640312 !avec g= gravité, k = cste von karman, sqrt_k = racine(k) REAL :: logzu_10,logzt_10,logzq_10,logzt_zu,logzq_zu,& u,cd_n10_rt,ch_n10_rt,ch_n10,chi,z_0m,z_0h,tv INTEGER :: i,j,m1,m2,m3,m4,m5 REAL :: phih,phid,phie REAL, parameter :: alpha=2.7e-3,& !alpha = a1 ds formulation smith beta = 1.42e-4,& !beta = a2 ds la formulation smith gamma = 7.64e-5,& !Large and yaeger 2006 !gamma = a3 ds formulation smith ! gamma = 13.09e-3,& !Large and yaeger 2004 q0 = 1.64474 ! q0 dzfinis mais pas reutilisée ds la suite ?? REAL, dimension(3) :: z integer, dimension(2) :: tmp_shape z = [zu, zt, zq] logzu_10 = LOG(zu/10.) logzt_10 = LOG(zt/10.) logzq_10 = LOG(zq/10.) logzt_zu = logzt_10 - logzu_10 logzq_zu = logzq_10 - logzu_10 cpa = 1004.67 !!! pas reutilise ds nvelle formulation... tv = t1 * (1+.608*q1) !!!!! tv pas reutiliser ds nvelle formulation... rho = p / (287.1 * t1 * (1.+.61*q1)) !!!! rho pas reutiliser ds nvelle formulation.... le = (2.501-.00237*(t1-273.15-dt))*1.e6 !!!!!! le pas reutiliser ds nvelle formulation.... u = max(du,.5) !!!! u = vitesse du vent 1 er niveau u_N = u !!!!!!! u_N utilisée pour le calcul des premier calcul... !!!!! initialisation des parametres et premier calcul de z_0 cd et ch if (mixte) then z_0m = 1.e-4 !!!!! pour l'initialisation on fixe z_0 a 10-4 cd_n10 = k**2 / (log(10/z_0m))**2 !!!!! on calcul ensuite cd a 10 m d'après smith cd_n10_rt = sqrt(cd_n10) else ! Large and yaeger 2006 cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur ! Large and yaeger 2004 ! cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) cd_n10_rt = sqrt(cd_n10) !!!!!! z_0m = 10. * exp(-k/cd_n10_rt) !!!!!!! on calcul le premier z_0 endif if ( dt .gt. 0. ) then ch_n10 = 0.018 * cd_n10_rt !!!!!!! si regime stable calcul du ch ! z_0h = 10. * exp(-k/ch_n10) else ch_n10 = 0.0327 * cd_n10_rt !!!!!!! si regime instable calcul du ch ! z_0h = 10. * exp(-k/ch_n10) endif ce_n10 = 3.46e-2 *cd_n10_rt !!!!!! calcul ensuite de ce dt_u = dt !!! def dq_u = dq !!! def cd = cd_n10 !!! def ch = ch_n10 !!! def ce = ce_n10 !!! def cd_rt = sqrt(cd) !!! def ! open(7,file="err_rel.dat",position="append") !!!!!!!!! début des itérations.... do i = 1,n_it usr = cd_rt * u tsr = ch / cd_rt * dt_u qsr = ce / cd_rt * dq_u zeta = k*g / (usr**2) * tsr / t1& * z zeta = sign( min(abs(zeta),10.0), zeta ) !!!! calcul du zeta pour connaitre le signe et donc en deduire le regime... IF ( zeta(1) .GT. 0 ) THEN !!! si regime stable alors valeurs pr chi et psi chi = 0.018 psi(1,1) = -5.*zeta(1) psi(2,1) = psi(1,1) ELSE !!!!! si regime instable alors valeurs pr psi et chi !!! chi = 0.0327 X = sqrt( sqrt( 1-16.*zeta(1) )) psi(1,1) = 3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X) psi(2,1) = 2.*LOG( (1 + X**2) / 2 ) END IF do j =2,3 IF ( zeta(j) .GT. 0 ) THEN psi(1,j) = -5.*zeta(j) psi(2,j) = psi(1,1) ELSE X = sqrt( sqrt( 1-16.*zeta(j) )) psi(1,j) = 3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X) psi(2,j) = 2.*LOG( (1 + X**2) / 2 ) END IF enddo dt_u = dt - tsr / k * & ( logzt_zu + psi(2,1) - psi(2,2) ) dq_u = dq - qsr / k * & ( logzq_zu + psi(2,1) - psi(2,3) ) u_N = u / ( 1 + cd_n10_rt / k * ( logzu_10 - psi(1,1)) ) if (mixte) then z_0m=0.018*usr**2/9.81 + 0.11*14e-6 / (usr) ! z_0h=(0.11*14e-6 / (usr)) + 1.4e-5 cd_n10 = k**2 / (log(10/z_0m))**2 cd_n10_rt = sqrt(cd_n10) ! ch_n10 = k**2 / ((log(10/z_0m))*(log(10/z_0h))) ! ch_n10_rt = sqrt(ch_n10) ch_n10 = 1.e-3 else ! Large and yaeger 2006 cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur ! Large and yaeger 2004 ! cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) ! cd_n10_dun = -alpha / u_n**2 + gamma cd_n10_rt = sqrt(cd_n10) 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 ! ch_n10 = chi * cd_n10_rt ! z_0h = 10. * exp(-k/ch_n10) ! ch_n10_rt = sqrt(ch_n10) ch_n10 = chi * cd_n10_rt endif ce_n10 = 3.46e-2*cd_n10_rt phid = 1+cd_n10_rt/k * (logzu_10 - psi(1,1)) phih = 1+chi/k * (logzu_10 - psi(2,1)) phie = 1+3.46e-2/k * (logzu_10 - psi(2,1)) cd_rt = cd_n10_rt / abs( phid ) ch = ch_n10 / ( phih * phid ) ce = ce_n10 / ( phie * phid ) cd=cd_rt**2 END DO usr = cd_rt * u tsr = ch / cd_rt * dt_u qsr = ce / cd_rt * dq_u ! coeffs(:,m) = [ & ! rho * usr**2,& ! rho * cpa * usr *tsr,& ! rho * le * usr *qsr] coeffs = [ cd_rt**2, ch , ce ] rugosm = z_0m rugosh = z_0h RETURN END subroutine clc_core_cp