Commit 981c58b6 authored by ldemaine's avatar ldemaine
Browse files

Vent de pente par le calcul de K (non testé)

parent e1d0d719
......@@ -1507,16 +1507,22 @@ c$$$ write(*,*) 'temp', temp
if (hireskey.eq.1) then
name='tsurf'
call var2d(tsurf,lon,lat,utime,name,ierr,itimint,wl,wh,
& ps_psgcm)
c$$$ call slope_wind_old(dset,xlon,xlat,zsurface,zradius,z_0,
c$$$ & tl,temp_gcm,sheight,
c$$$ & ul,vl,old_zonwind,old_merwind,
c$$$ & upslope_wind,crossslope_wind)
call slope_wind(dset,xlon,xlat,zsurface,zradius,z_0,
& tl,temp_gcm,sheight,
& tl,temp_gcm,tsurf,p_pgcm,ps_hr,
& levhi,levlow,sheight,
& utime,itimint,wl,wh,
& ul,vl,old_zonwind,old_merwind,
& upslope_wind,crossslope_wind)
c$$$ call slope_wind(dset,xlon,xlat,zsurface,z_0,theta_s,
c$$$ & psy_s,tl,temp_0,ul,vl,
c$$$ & old_zonwind,old_merwind,upslope_wind,crossslope_wind)
extvar(82) = upslope_wind
extvar(83) = crossslope_wind
extvar(84) = old_merwind
......@@ -7711,7 +7717,7 @@ c$$$ write(*,*) 'Ftot:', F
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine slope_wind(dset,xlon,xlat,zsurface,zradius,zsurf_0,
subroutine slope_wind_old(dset,xlon,xlat,zsurface,zradius,zsurf_0,
& temp,temp_gcm,sheight,
& zonwind,merwind,old_zonwind,old_merwind,
& upslope_wind,crossslope_wind)
......@@ -7769,8 +7775,8 @@ c$$$ real,intent(in) :: pratio
real :: uxu
real :: uyu
real :: k_vk,g_0
parameter (k_vk=0.4)
real :: karman,g_0
parameter (karman=0.4)
parameter (g_0=3.72)
real pi,degtorad
......@@ -7906,7 +7912,7 @@ c$$$ write(1,*) '| Float', Float
write(*,*) 'zsurface.le.zsurf_0'
endif
Cd = (k_vk/log(zsurface/zsurf_0))**2
Cd = (karman/log(zsurface/zsurf_0))**2
c$$$ write(1,*) '| Cd', Cd
write(*,*) '| Cd', Cd
......@@ -7973,3 +7979,345 @@ c$$$ close(1)
return
end
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine slope_wind(dset,xlon,xlat,zsurface,zradius,zsurf_0,
& temp,temp_gcm,tsurf,p_pgcm,ps_hr,
& levhi,levlow,sheight,
& utime,itimint,wl,wh,
& zonwind,merwind,old_zonwind,old_merwind,
& upslope_wind,crossslope_wind)
implicit none
include "constants_mcd.inc"
! TODO: profil de la température
! TODO: calcul de K pour chaque couche
! Inputs
character*(*),intent(in) :: dset ! Path to MCD datafiles
real,intent(in) :: xlon
real,intent(in) :: xlat
real,intent(in) :: zsurface
real,intent(in) :: zradius ! to compute get_slopes
real,intent(in) :: zsurf_0
real,intent(in) :: temp
real,intent(in) :: temp_gcm(dimlevs)
real,intent(in) :: tsurf
real,intent(in) :: p_pgcm(dimlevs),ps_hr
c$$$ !inutile pour le moment car temp déjà un imput
c$$$ real,intent(in) :: levweight(2)
integer,intent(in) :: levhi
integer,intent(in) :: levlow
c$$$ real,intent(in) :: pratio
real,intent(in) :: sheight(dimlevs) ! altitude above surface of sigma levels
real,intent(in) :: utime !Universal time (0. to 24. hrs)=local time at lon=0
integer,intent(in) :: itimint !seasonal interpolation flags
real,intent(in) :: wl,wh !seasonal interpolation weights
real,intent(inout) :: zonwind
real,intent(inout) :: merwind
real,intent(out) :: old_zonwind
real,intent(out) :: old_merwind
real,intent(out) :: upslope_wind
real,intent(out) :: crossslope_wind
integer :: ier ! need intent(out)
real :: theta_s
real :: psy_s
real :: slopes_scale = 1
c$$$
c$$$ ! varibles pour la référence
real :: zref
integer :: levhi_0
integer :: levlow_0
real :: levweight_0
real :: pratio_0
integer :: l
real temp_0
real u_gcm(dimlevs),v_gcm(dimlevs)
character*16 tmpname
integer :: ierr
real :: Cd
real :: Flotation
real :: pKv(dimlevs),pKv_loc
real Flotabis,Sciz
real upwind_gcm(dimlevs),crosswind_gcm(dimlevs)
integer it
c$$$ real :: UU
c$$$ integer :: i
c$$$ real :: uxu
c$$$ real :: uyu
c$$$
real :: karman,g_0
parameter (karman=0.4)
parameter (g_0=3.72)
real pi,degtorad
parameter (pi=3.14159265358979d0)
parameter (degtorad=pi/180.0d0)
c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
upslope_wind=0
crossslope_wind=0
old_merwind=merwind
old_zonwind = zonwind
slopes_scale = 1
call get_slopes(dset,xlon,xlat,slopes_scale,zradius,
& theta_s,psy_s,ier)
zref = 10000*sin(theta_s*degtorad)
if (zsurface.le.zref.and.zsurface.gt.0) then
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! calcul des valeurs pour la ref (sur le modèle de getsi)
!
! ATENTION pratio_0 non calculé dans les cas extrèmes pour l'instant
! (et non pris en compte dans les calculs de temp_0
if (zref.lt.sheight(1)) then
! ! below the lowest layer
!c$$$ write(*,*) 'ATTENTION! zref trop bas'
levhi_0=1
levlow_0=1
levweight_0=1.
pratio_0 = 1.
elseif (zref.ge.sheight(dimlevs)) then
!c$$$ write(*,*) 'ATTENTION! zref trop haut'
! ! above the top layer
levhi_0=dimlevs
levlow_0=dimlevs
levweight_0=0.
pratio_0 = 1.
else
! ! general case: sheight(1)<=absheight<=sheight(dimlevs)
do l=1,dimlevs-1
if ((zref.ge.sheight(l))
& .and.(zref.lt.sheight(l+1))) then
levhi_0=l+1
levlow_0=l
levweight_0=(zref-sheight(levlow_0))
& /(sheight(levhi_0)-sheight(levlow_0))
pratio_0=1.
endif
enddo
endif
temp_0 = temp_gcm(levlow_0)
& + (temp_gcm(levhi_0)-temp_gcm(levlow_0))*levweight_0
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
if(levhi.eq.1) then
Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)/temp_0
Cd = (karman/log(zsurface/zsurf_0))**2
upslope_wind = sign(sqrt(abs(Flotation / Cd)), Flotation)
else !(levhi.eq.1)
! Calcul des valeurs au premier niveau gcm avec le coef Cd
Flotation=g_0*(temp_gcm(1) - temp_0)*sin(theta_s*degtorad)
& /temp_0
Cd = (karman/log(sheight(1)/zsurf_0))**2
! Calcul du premier niveau de vent
upwind_gcm(1)=sign(sqrt(abs(Flotation / Cd)), Flotation)
! get wind profile at same location
tmpname="u"
call profi(u_gcm,xlon,xlat,utime,tmpname,ierr,itimint,
& wl,wh,0,0,0)
tmpname="v"
call profi(v_gcm,xlon,xlat,utime,tmpname,ierr,itimint,
& wl,wh,0,0,0)
! Calcul de la première valeur de K (K(1))
call vdif_k(sheight(1),0.,U_gcm(1),0.,V_gcm(1),0.,
& Temp_gcm(1),tsurf,p_pgcm(1),ps_hr,ps_hr,
& zsurf_0,pKv(1))
if (levlow.eq.1) then
! K entre la dernière couche et la couche actuelle (pKv_loc)
call vdif_k(sheight(2),sheight(1),
& U_gcm(2),U_gcm(1),V_gcm(1),V_gcm(1),
& Temp_gcm(2),Temp_gcm(1),P_pgcm(2),P_pgcm(1),ps_hr,
& zsurf_0,pKv_loc)
Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)
& /temp_0
Flotabis = (zsurface - sheight(1))*(zsurface)*Flotation
& /(2*pKv_loc)
Sciz=(pkv(1)/pkv_loc)
& *((zsurface-sheight(1))/(sheight(1)))
& *(upwind_gcm(1))
upslope_wind = upwind_gcm(levlow) + Flotabis + Sciz
else !(levlow.eq.1)
do it=2,levlow
! K(it)
call vdif_k(sheight(it),sheight(it-1),
& U_gcm(it),U_gcm(it-1),V_gcm(it),V_gcm(it-1),
& Temp_gcm(it),Temp_gcm(it-1)
& ,P_pgcm(it),P_pgcm(it-1),ps_hr,
& zsurf_0,pKv(it))
Flotation = g_0*(temp_gcm(it)-temp_0)
& *sin(theta_s*degtorad)/temp_0
if (it.eq.2) then
Sciz = (pkv(1)/pkv(2))
& *((sheight(2)-sheight(1))/sheight(1))
& *upwind_gcm(1)
Flotabis = (sheight(2)-sheight(1))*sheight(2)
& *Flotation/(2*pkv(2))
else !(it.eq.2)
Sciz = (pkv(it-1)/pkv(it))
& *((sheight(it)-sheight(it-1))
& /(sheight(it-1)-sheight(it-2)))
& *(upwind_gcm(it-1)-upwind_gcm(it-2))
Flotabis = (sheight(it)-sheight(it-1))
& *(sheight(it)-sheight(it-2))
& *Flotation/(2*pkv(it))
endif !(it.eq.2)
upwind_gcm(it) = upwind_gcm(it-1) + Flotabis + Sciz
enddo !it=2,levlow
! K entre la dernière couche et la couche actuelle (pKv_loc)
call vdif_k(zsurface,sheight(levlow),
& U_gcm(levlow),U_gcm(levlow-1),
& V_gcm(levlow),V_gcm(levlow-1),
& Temp_gcm(levlow),Temp_gcm(levlow-1),
& P_pgcm(levlow),P_pgcm(levlow-1),ps_hr,
& zsurf_0,pKv_loc)
Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)
& /temp_0
Flotabis = (zsurface - sheight(levlow))
& *(zsurface - sheight(levlow-1))
& *Flotation/(2*pKv_loc)
Sciz = (pkv(levlow)/pkv_loc)
& *((zsurface-sheight(levlow))
& /(sheight(levlow)-sheight(levlow-1)))
& *(upwind_gcm(levlow)-upwind_gcm(levlow-1))
upslope_wind = upwind_gcm(levlow) + Flotabis + Sciz
endif !(levlow.eq.1)
endif !(levhi.eq.1)
endif !(zsurface.le.zref.and.zsurface.gt.0)
merwind = old_merwind
& + upslope_wind*cos(psy_s*degtorad)
& + crossslope_wind*sin(psy_s*degtorad)
zonwind = old_zonwind
& + upslope_wind*sin(psy_s*degtorad)
& + crossslope_wind*cos(psy_s*degtorad)
return
end
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine vdif_k(z_hi,z_low,U_hi,u_low,V_hi,V_low,
& temp_hi,temp_low,P_hi,P_low,Psurf,z_0,
& pKv)
real,intent(in) :: z_hi,z_low
real,intent(in) :: U_hi,U_low
real,intent(in) :: V_hi,V_low
real,intent(in) :: Temp_low,Temp_hi
real,intent(in) :: P_hi,P_low
real,intent(in) :: Psurf
real,intent(in) :: z_0
real,intent(out) :: pKv
real :: zdu,zdv,zdz
real :: zdvodz2
real :: rcp
real :: T_pot_low,T_pot_hi
real :: T_pot_mid,dT_pot_dz
real :: z1
real :: Ri
real :: lmix,lmixmin
real :: emin_turb
real :: karman, g_0
parameter (karman=0.4)
parameter (g_0=3.72)
lmixmin = 30.
emin_turb = 1.e-8
z1 = (z_hi + z_low)/2 + zsurf_0
lmix = karman*z1/(1.+karman*z1/lmixmin)
zdu = U_hi - U_low
zdv = V_hi - V_low
zdz = z_hi - z_low
zdvodz2 = (zdu**2 + zdv**2)/(zdz**2)
! need declaration
rcp=0.259 ! on considère que gamma = 1.35 avec rcp = (1-gamma)/gamma, ce qui est valable à 3% près en dessous de 5000m
T_pot_hi=Temp_hi*(Psurf/P_hi)**rcp
T_pot_low=Temp_low*(Psurf/P_low)**rcp
T_pot_mid = (T_pot_hi + T_pot_low)/2.
dT_pot_dz = (T_pot_hi - T_pot_low)/zdz
if(zdvodz2.lt.1.e-5) then
pKv = lmix*sqrt(emin_turb)
else
Ri = g_0 * dT_pot_dz /(T_pot_mid*zdvodz2)
pKv = lmix*sqrt(max(lmix**2*zdvodz2*(1-2.5*Ri),emin_turb))
endif
return
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment