Commit 11771146 authored by ldemaine's avatar ldemaine
Browse files

série de write(*,*) pour trouver débuger

parent 981c58b6
......@@ -517,7 +517,8 @@ c Hires wind
real crossslope_wind
real zref
real upslope_wind_cd
real crossslope_wind_cd
......@@ -1502,38 +1503,37 @@ c******************************
endif
endif ! of if (perturkey.eq.5)
c (re-)compute atmospheric pressure (deplacer depuis 4.4)
pl=ps*(sigma(levlow)+(sigma(levhi)-sigma(levlow))*levweight(2))
pl=pratio*pl
c 4.3.5 Computation of high precision wind
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_old(dset,xlon,xlat,zsurface,zradius,z_0,
& tl,temp_gcm,sheight,
& ul,vl,old_zonwind,old_merwind,
& upslope_wind_cd,crossslope_wind_cd)
write(*,*) 'aa' , upslope_wind_cd
call slope_wind(dset,xlon,xlat,zsurface,zradius,z_0,
& tl,temp_gcm,tsurf,p_pgcm,ps_hr,
& tl,temp_gcm,tsurf,pl,p_pgcm,ps_hr,
& levhi,levlow,sheight,
& utime,itimint,wl,wh,
& ul,vl,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
c 4.4 Store atmospheric fields (pressure, temperature, density, winds)
c (re-)compute atmospheric pressure
pl=ps*(sigma(levlow)+(sigma(levhi)-sigma(levlow))*levweight(2))
pl=pratio*pl
t=tl ! temperature
p=pl ! pressure
......@@ -2416,6 +2416,9 @@ c copy to output
extvar(83) = crossslope_wind
extvar(84) = old_merwind
extvar(85) = old_zonwind
extvar(86) = upslope_wind_cd
extvar(87) = crossslope_wind_cd
return
c Error handling : all the outputs are set to a missing data value
......@@ -7982,7 +7985,7 @@ c$$$ close(1)
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine slope_wind(dset,xlon,xlat,zsurface,zradius,zsurf_0,
& temp,temp_gcm,tsurf,p_pgcm,ps_hr,
& temp,temp_gcm,tsurf,Pl,p_pgcm,ps_hr,
& levhi,levlow,sheight,
& utime,itimint,wl,wh,
& zonwind,merwind,old_zonwind,old_merwind,
......@@ -8008,6 +8011,7 @@ c$$$ close(1)
real,intent(in) :: zradius ! to compute get_slopes
real,intent(in) :: zsurf_0
real,intent(in) :: temp
real,intent(in) :: Pl
real,intent(in) :: temp_gcm(dimlevs)
real,intent(in) :: tsurf
real,intent(in) :: p_pgcm(dimlevs),ps_hr
......@@ -8072,6 +8076,23 @@ c$$$
parameter (degtorad=pi/180.0d0)
write(*,*) ' '
write(*,*) '|--------------------------------------------------|'
write(*,*) '| Input slope_wind |'
write(*,*) '|--------------------------------------------------|'
write(*,*) '| xlon, xlat:', xlon, xlat
write(*,*) '| zsurface, zradius, zsurf_O:',zsurface,zradius,
& zsurf_0
write(*,*) '| levhi, levlow:', levhi, levlow
write(*,*) '| utime:', utime
write(*,*) '| itimint, wl, wh:', itimint, wl, wh
write(*,*) '| zonwind, merwind:', zonwind, merwind
write(*,*) '| temp, tsurf:', temp, tsurf
write(*,*) '| temp_gcm:', temp_gcm
write(*,*) '| sheight:', sheight
write(*,*) '|--------------------------------------------------|'
write(*,*) '| '
c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
upslope_wind=0
......@@ -8080,13 +8101,27 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
old_merwind=merwind
old_zonwind = zonwind
write(*,*) '| upslope_wind, crossslope_wind:'
& ,upslope_wind, crossslope_wind
write(*,*) '| old_merwind, merwind, old_zonwind, zonwind:'
& ,old_merwind, merwind, old_zonwind, zonwind
write(*,*) '| '
slopes_scale = 1
call get_slopes(dset,xlon,xlat,slopes_scale,zradius,
& theta_s,psy_s,ier)
zref = 10000*sin(theta_s*degtorad)
write(*,*) '| slopes_scale, theta_s, psy_s:',
& slopes_scale, theta_s, psy_s
write(*,*) '| zref:', zref
write(*,*) '| '
if (zsurface.le.zref.and.zsurface.gt.0) then
write(*,*) '| | if (zsurface.le.zref.and.zsurface.gt.0)'
write(*,*) '|-| - - - - - - - - - - - - - - - - - - - - - - -'
write(*,*) '| |'
......@@ -8097,21 +8132,27 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
! (et non pris en compte dans les calculs de temp_0
if (zref.lt.sheight(1)) then
write(*,*) '| | | if (zref.lt.sheight(1))'
write(*,*) '|-|-|- - - - - - - - - - - - - - - - - - - - -'
! ! 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
write(*,*) '| | | elseif (zref.ge.sheight(dimlevs))'
write(*,*) '|-|-|- - - - - - - - - - - - - - - - - - - - -'
!c$$$ write(*,*) 'ATTENTION! zref trop haut'
! ! above the top layer
levhi_0=dimlevs
levlow_0=dimlevs
levweight_0=0.
pratio_0 = 1.
else
else
write(*,*) '| | | else'
write(*,*) '|-|-|- - - - - - - - - - - - - - - - - - - - -'
! ! general case: sheight(1)<=absheight<=sheight(dimlevs)
do l=1,dimlevs-1
if ((zref.ge.sheight(l))
......@@ -8123,24 +8164,35 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
pratio_0=1.
endif
enddo
endif
endif
temp_0 = temp_gcm(levlow_0)
& + (temp_gcm(levhi_0)-temp_gcm(levlow_0))*levweight_0
write(*,*) '| | | levhi_0, levlow_0, levweight_0, pratio_0:',
& levhi_0,levlow_0,levweight_0,pratio_0
write(*,*) '| | temp_0:', temp_0
write(*,*) '| |'
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
if(levhi.eq.1) then
if (levhi.eq.1) then
write(*,*) '| | | if (levhi.eq.1)'
write(*,*) '| | |- - - - - - - - - - - - - - - - - - - - -'
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)
write(*,*) '| | | Flotation, Cd, upslope_wind:',
& Flotation, Cd, upslope_wind
write(*,*) '| | |'
else !(levhi.eq.1)
write(*,*) '| | | if not (levhi.eq.1)'
write(*,*) '| | |- - - - - - - - - - - - - - - - - - - - -'
! 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
......@@ -8150,6 +8202,10 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
! Calcul du premier niveau de vent
upwind_gcm(1)=sign(sqrt(abs(Flotation / Cd)), Flotation)
write(*,*) '| | | Flotation, Cd, upwind_gcm(1):',
& Flotation, Cd, upwind_gcm(1)
write(*,*) '| | |'
! get wind profile at same location
tmpname="u"
call profi(u_gcm,xlon,xlat,utime,tmpname,ierr,itimint,
......@@ -8157,18 +8213,29 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
tmpname="v"
call profi(v_gcm,xlon,xlat,utime,tmpname,ierr,itimint,
& wl,wh,0,0,0)
write(*,*) '| | | u_gcm:', u_gcm
write(*,*) '| | |'
write(*,*) '| | | v_gcm:', v_gcm
write(*,*) '| | |'
! 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))
write(*,*) '| | | pKv(1):', pKv(1)
write(*,*) '| | |'
if (levlow.eq.1) then
write(*,*) '| | | | if levlow.eq.1'
write(*,*) '|-|-|-| - - - - - - - - - - - - - - - - - -'
! 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)
call vdif_k(sheight(2), sheight(1),
& U_gcm(2), U_gcm(1),
& V_gcm(2), 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
......@@ -8182,14 +8249,24 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
upslope_wind = upwind_gcm(levlow) + Flotabis + Sciz
write(*,*) '| | | | pKv_loc', pKv_loc
write(*,*) '| | | | Flotation, Flotabis, Sciz:',
& Flotation, Flotabis, Sciz
write(*,*) '| | | | upslope_wind:', upslope_wind
write(*,*) '| | | | '
else !(levlow.eq.1)
write(*,*) '| | | | if not levlow.eq.1'
write(*,*) '|-|-|-| - - - - - - - - - - - - - - - - - -'
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,
write(*,*) '| | | | | do it =2,levlow'
write(*,*) '|-|-|-|-|- - - - - - - - - - - - - - - -'
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)
......@@ -8217,15 +8294,22 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
upwind_gcm(it) = upwind_gcm(it-1) + Flotabis + Sciz
write(*,*) '| | | | | pKv(it):',pKv(it)
write(*,*) '| | | | | Flotation, Flotabis, Sciz:',
& Flotation, Flotabis, Sciz
write(*,*) '| | | | | upwind_gcm(it):',upwind_gcm(it)
write(*,*) '| | | | |'
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)
call vdif_k(zsurface, sheight(levlow),
& merwind, U_gcm(levlow),
& zonwind, V_gcm(levlow),
& temp, Temp_gcm(levlow),
& Pl, P_pgcm(levlow), ps_hr,
& zsurf_0,
& pKv_loc)
Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)
& /temp_0
......@@ -8241,6 +8325,12 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
upslope_wind = upwind_gcm(levlow) + Flotabis + Sciz
write(*,*) '| | | | pKv_loc', pKv_loc
write(*,*) '| | | | Flotation, Flotabis, Sciz:',
& Flotation, Flotabis, Sciz
write(*,*) '| | | | upslope_wind:', upslope_wind
write(*,*) '| | | | '
endif !(levlow.eq.1)
endif !(levhi.eq.1)
......@@ -8254,14 +8344,19 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
& + upslope_wind*sin(psy_s*degtorad)
& + crossslope_wind*cos(psy_s*degtorad)
write(*,*) '| merwind, zonwind:', merwind, zonwind
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,
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
......@@ -8292,6 +8387,15 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
parameter (karman=0.4)
parameter (g_0=3.72)
write(*,*) ' > ********************************************'
write(*,*) ' > Input vdif_k'
write(*,*) ' > ********************************************'
write(*,*) ' > z_0, z_hi, z_low, U_low, U_hi, V_low, V_hi:',
& z_0, z_hi, z_low, U_low, U_hi, V_low, V_hi
write(*,*) ' > temp_low, Temp_hi, P_low, P_hi, Psurf:',
& temp_low, Temp_hi, P_low, P_hi, Psurf
write(*,*) ' >'
lmixmin = 30.
emin_turb = 1.e-8
......@@ -8303,6 +8407,9 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
zdv = V_hi - V_low
zdz = z_hi - z_low
zdvodz2 = (zdu**2 + zdv**2)/(zdz**2)
write(*,*) ' > z1, lmix, zdu, zdv, zdz, zdvodz2:',
& z1, lmix, zdu, zdv, zdz, zdvodz2
write(*,*) ' >'
! 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
......@@ -8310,13 +8417,31 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
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
write(*,*) ' > T_pot_hi, T_pot_low, T_pot_mid, dT_pot_dz',
& T_pot_hi, T_pot_low, T_pot_mid, dT_pot_dz
write(*,*) ' >'
write(*,*) ' >'
write(*,*) ' > zdvodz2', zdvodz2
write(*,*) ' > lmix', lmix
write(*,*) ' >'
write(*,*) ' > dT_pot_dz', dT_pot_dz
write(*,*) ' > g_0', g_0
write(*,*) ' >'
if(zdvodz2.lt.1.e-5) then
pKv = lmix*sqrt(emin_turb)
else
Ri = g_0 * dT_pot_dz /(T_pot_mid*zdvodz2)
write(*,*) ' > * Ri', Ri
write(*,*) ' > 1 - 2.5 Ri', 1.-2.5*Ri
write(*,*) ' >'
pKv = lmix*sqrt(max(lmix**2*zdvodz2*(1-2.5*Ri),emin_turb))
endif
write(*,*) ' > pKv', pKv
write(*,*) ' > ********************************************'
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