Commit 7850a717 authored by Leo Demaine's avatar Leo Demaine
Browse files

modification de slope_wind avec les remarques de mardi dernier (non testée)

parent 08c3da38
......@@ -1506,30 +1506,22 @@ c 4.3.5 Computation of high precision wind
if (hireskey.eq.1) then
zref = 10000
slopes_scale = 1
call get_slopes(dset,xlon,xlat,slopes_scale,zradius,
& theta_s,psy_s,ier)
call slope_wind(dset,xlon,xlat,zsurface,zradius,
& temp,temp_gcm,sheight,
& zonwind,merwind,old_zonwind,old_merwind,
& upslope_wind,crossslope_wind)
if (zsurface.le.zref) then
temp_0=207
call slope_wind(dset,xlon,xlat,zsurface,z_0,theta_s,
& psy_s,tl,temp_0,ul,vl,
& old_zonwind,old_merwind,upslope_wind,crossslope_wind)
else
old_zonwind = zonwind
old_merwind = merwind
upslope_wind = 0
crossslope_wind = 0
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
extvar(85) = old_zonwind
endif
endif
c 4.4 Store atmospheric fields (pressure, temperature, density, winds)
......@@ -7696,8 +7688,9 @@ c$$$ write(*,*) 'Ftot:', F
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine slope_wind(dset,xlon,xlat,zsurface,zsurf_0,theta_s,
& psy_s,temp,temp_0,zonwind,merwind,old_zonwind,old_merwind,
subroutine slope_wind(dset,xlon,xlat,zsurface,zradius,
& temp,temp_gcm,sheight,
& zonwind,merwind,old_zonwind,old_merwind,
& upslope_wind,crossslope_wind)
implicit none
......@@ -7710,11 +7703,16 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real,intent(in) :: xlon
real,intent(in) :: xlat
real,intent(in) :: zsurface
real,intent(in) :: zsurf_0
real,intent(in) :: theta_s
real,intent(in) :: psy_s
c$$$ real,intent(in) :: zsurf_0
real,intent(in) :: zradius ! to compute get_slopes
real,intent(in) :: temp
real,intent(in) :: temp_0
real,intent(in) :: temp_gcm(dimlevs)
!inutile pour le moment car temp déjà un imput
c$$$ real,intent(in) :: levweight(2)
c$$$ integer,intent(in) :: levhi
c$$$ integer,intent(in) :: levlow
c$$$ real,intent(in) :: pratio
real,intent(in) :: sheight(dimlevs) ! altitude above surface of sigma levels
real,intent(inout) :: zonwind
real,intent(inout) :: merwind
......@@ -7724,6 +7722,20 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
! 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 :: Float
......@@ -7742,67 +7754,164 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
parameter (degtorad=pi/180.0d0)
write(*,*) ''
write(*,*) '|-----------------------------------------|'
write(*,*) '| Input slope_wind |'
write(*,*) '|-----------------------------------------|'
write(*,*) '|'
write(*,*) '| xlon', xlon
write(*,*) '| xlat', xlat
write(*,*) '| zsurface', zsurface
write(*,*) '| zsurf_0', zsurf_0
write(*,*) '| theta_s', theta_s
write(*,*) '| psy_s', psy_s
write(*,*) '| temp', temp
write(*,*) '| temp_0', temp_0
write(*,*) '| zonwind', zonwind
write(*,*) '| merwind', merwind
write(*,*) '|'
write(*,*) '|-----------------------------------------|'
write(*,*) '|'
slopes_scale = 1
call get_slopes(dset,xlon,xlat,slopes_scale,zradius,
& theta_s,psy_s,ier)
c$$$! TODO: cas où l'altitude est inférieure à la première maille de mcd
c$$$! exemple d'utilisation pour trouver la température
c$$$ elseif ((name.eq.'temp')) then
c$$$! Default linear interpolation between levels:
c$$$ a=profile(levlow)+(profile(levhi)-profile(levlow))*levweight(1)
c$$$! If below first atmospheric level, use linear
c$$$! interpolation with surface temperature
c$$$ if ((levlow.eq.levhi).and.(levlow.eq.1)) then
c$$$ goto 1000
c$$$ endif ! of if ((levlow.eq.levhi).and.(levlow.eq.1))
c$$$ elseif ((name.eq.'u').or.(name.eq.'v')) then
c$$$! Default linear interpolation between levels:
c$$$ a=profile(levlow)+(profile(levhi)-profile(levlow))*levweight(1)
c$$$! If below first atmopheric level, mimic log boundary layer
c$$$! with roughness length z_0; below z_0, set wind velocity to zero
c$$$ if ((levlow.eq.levhi).and.(levlow.eq.1)) then
c$$$ goto 1000
c$$$ endif
open(1, file='Test_slope/TEST_1.txt', status = 'new')
write(1,*) ''
write(1,*) '|-----------------------------------------|'
write(1,*) '| Input slope_wind |'
write(1,*) '|-----------------------------------------|'
write(1,*) '|'
write(1,*) '| xlon', xlon
write(1,*) '| xlat', xlat
write(1,*) '| zsurface', zsurface
write(1,*) '| zsurf_0', zref
write(1,*) '| theta_s', theta_s
write(1,*) '| psy_s', psy_s
write(1,*) '| temp', temp
c$$$ write(1,*) '| temp_0', temp_0
write(1,*) '| zonwind', zonwind
write(1,*) '| merwind', merwind
write(1,*) '|'
write(1,*) '|-----------------------------------------|'
write(1,*) '|'
upslope_wind=0
crossslope_wind=0
! Computation of float terme
Float = g_0*(temp - temp_0)*sin(theta_s*degtorad)/temp_0
write(*,*) '| Float', Float
zref = 10000*sin(theta_s)
if (zsurface.le.zref) then
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! calcul des valeurs pour la ref (sur le modèle de getsi)
!
! pratio_0 non calculé dans les cas extrèmes pour l'instant (et non pris en compte dans les calculs de temp_0
! besion de : R_gcm pour pratio_0
if (zref.lt.sheight(1)) then
! ! below the lowest layer
write(*,*) '/!\ zref trop bas'
levhi_0=1
levlow_0=1
levweight_0=1.
c$$$ levweight_0(2)=0.
c$$$ pratio_0=exp((sheight(1)-zref)*g/(R_gcm(1)*temp_gcm(1)))
pratio_0 = 1.
elseif (zref.ge.sheight(dimlevs)) then
write(*,*) '/!\ zref trop haut'
! ! above the top layer
levhi_0=dimlevs
levlow_0=dimlevs
levweight_0=0.
c$$$ levweight_0(2)=1.
! pratio_0=exp((sheight(dimlevs)-absheight)*g/(R_gcm(dimlevs)
! & *temp_gcm(dimlevs)))
! z = (r1/r0)**2*r0*dz/(r0+dz) where r1 is alt where g was calculated
c$$$ z = (a0+oroheight+absheight)**2/(a0+oroheight+sheight(dimlevs))*
c$$$ & (absheight-sheight(dimlevs))/(a0+oroheight+absheight)
c$$$ pratio_0 = 0.0
c$$$ do igas=1,nbvar3d
c$$$ ! Pressure is computed assuming hydrostatic for each species
c$$$ if(ma3d(igas).lt.0) cycle
c$$$ call profi(vmr_gcm,lon,lat,utime,varname3d(igas),
c$$$ & ierr,itimint,wl,wh,0,dimlevs,dimlevs)
c$$$ pratio_0 = pratio_0+vmr_gcm(dimlevs)*
c$$$ & exp(-z*g/(8314.4598/ma3d(igas)*temp_gcm(dimlevs)))
c$$$ enddo
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))
c$$$ levweight_0(2)=(1.-(sigma(levhi_0)/sigma(levlow_0))**levweight(1))
c$$$ & /(1.-(sigma(levhi_0)/sigma(levlow_0)))
pratio_0=1.
endif
enddo
endif
Cd = (k_vk/log(zsurface/zsurf_0))**2
write(*,*) '| Cd', Cd
temp_0 = temp_gcm(levlow_0)
& + (temp_gcm(levhi_0)-temp_gcm(levlow_0))*levweight_0
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! Computation of float terme
Float = g_0*(temp - temp_0)*sin(theta_s*degtorad)/temp_0
write(1,*) '| Float', Float
do i=0,1000
write(*,*)'|-|------------'
write(*,*)'| | i = ', i
write(*,*)'|-|------------'
uxu = merwind
& + upslope_wind*cos(psy_s*degtorad)
& +crossslope_wind*sin(psy_s*degtorad)
write(*,*)'| | uxu = ', uxu
uyu = zonwind
& + upslope_wind*sin(psy_s*degtorad)
& +crossslope_wind*cos(psy_s*degtorad)
write(*,*)'| | uyu = ', uyu
UU = sqrt( uxu**2 + uyu**2 )
write(*,*)'| | UU = ', uu
Cd = (k_vk/log(zsurface/zref))**2
write(1,*) '| Cd', Cd
upslope_wind = Float*zsurface/(Cd*UU)
write(*,*) '| | Float*zsurface/(Cd*UU)',Float*zsurface/(Cd*UU)
write(*,*)'| | upslope_wind = ', upslope_wind
enddo
upslope_wind = sign(sqrt(Float / Cd), Float/Cd)
c$$$ do i=0,1000
c$$$
c$$$ write(1,*)'|-|------------'
c$$$ write(1,*)'| | i = ', i
c$$$ write(1,*)'|-|------------'
c$$$
c$$$ uxu = merwind
c$$$ & + upslope_wind*cos(psy_s*degtorad)
c$$$ & +crossslope_wind*sin(psy_s*degtorad)
c$$$ write(1,*)'| | uxu = ', uxu
c$$$
c$$$ uyu = zonwind
c$$$ & + upslope_wind*sin(psy_s*degtorad)
c$$$ & +crossslope_wind*cos(psy_s*degtorad)
c$$$ write(1,*)'| | uyu = ', uyu
c$$$
c$$$ UU = sqrt( uxu**2 + uyu**2 )
c$$$ write(1,*)'| | UU = ', uu
c$$$
c$$$ upslope_wind = Float*zsurface/(Cd*UU)
c$$$ write(1,*) '| | Float*zsurface/(Cd*UU)',Float*zsurface/(Cd*UU)
c$$$ write(1,*)'| | upslope_wind = ', upslope_wind
c$$$
c$$$ enddo
endif
old_merwind = merwind
merwind = old_merwind
......@@ -7814,21 +7923,25 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
& + upslope_wind*sin(psy_s*degtorad)
& +crossslope_wind*cos(psy_s*degtorad)
write(*,*) '|'
write(*,*) '|-----------------------------------------|'
write(*,*) '| Output slope_wind |'
write(*,*) '|-----------------------------------------|'
write(*,*) '|'
write(*,*) '| upslope_wind', upslope_wind
write(*,*) '| crossslope_wind', crossslope_wind
write(*,*) '| zonwind', zonwind
write(*,*) '| merwind', merwind
write(*,*) '| old_zonwind', old_zonwind
write(*,*) '| old_merwind', old_merwind
write(*,*) '|'
write(*,*) '|-----------------------------------------|'
write(*,*) ''
write(1,*) '|'
write(1,*) '|-----------------------------------------|'
write(1,*) '| Output slope_wind |'
write(1,*) '|-----------------------------------------|'
write(1,*) '|'
write(1,*) '| upslope_wind', upslope_wind
write(1,*) '| crossslope_wind', crossslope_wind
write(1,*) '| zonwind', zonwind
write(1,*) '| merwind', merwind
write(1,*) '| old_zonwind', old_zonwind
write(1,*) '| old_merwind', old_merwind
write(1,*) '|'
write(1,*) '|-----------------------------------------|'
write(1,*) ''
close(1)
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