***s/r eole_3d
*
subroutine eole_3d( geobus, ngeop ) 1,14
implicit none
*
**
#include "nesting.cdk"
#include "dynmem.cdk"
#include "tracers.cdk"
#include "nestpnt.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
c#include "yomdyn2.cdk"
#include "consdyn_8.cdk"
#include "levels.cdk"
#include "partopo.cdk"
c***
#include "topo.cdk"
#include "h_geop.cdk"
#include "eole.cdk"
*
integer i,j,k,n,gnix,gnjy,dim,err,kc
real*8 c1,tropo
real tp1,tp2,qp1,s,d,qstar
pointer (patp1, tp1(gni+2*hx,gnj+2*hy,gnk)),
$ (patp2, tp2(gni+2*hx,gnj+2*hy,gnk)),
$ (paqp1, qp1(gni+2*hx,gnj+2*hy,0:gnk)),
$ (pas1 ,s(*)), (pad1 ,d(*))
* to hold the global sby, sbx
real sby1(gni+2*hx,gnj+2*hy),sbx1(gni+2*hx,gnj+2*hy)
real*8 rtstar1
c
integer ngeop
real geobus(ngeop)
c*** use new variables Tprime and Qprime (beginning), Jan. 2003
c
integer km1, kp1
real*8 beta, dh, dqdz, pt5
real hht(minx:maxx,miny:maxy,gnk), hhm(minx:maxx,miny:maxy,0:gnk)
real hhw(minx:maxx,miny:maxy, gnk)
pt5 = 0.5d0
c
c*** use new variables Tprime and Qprime (end), Jan. 2003
c
*
****
* Mise en place des profils
****
print*, 'NOT ADAPTED FOR NEW LATERAL BCS --- ABORT ---'
call mc2stop(-1)
gnix = gni+2*hx
gnjy = gnj+2*hy
if (myproc.eq.0) then
call hpalloc (patp1, gnix*gnjy*gnk, err, 1)
call hpalloc (patp2, gnix*gnjy*gnk, err, 1)
call hpalloc (paqp1, gnix*gnjy*(gnk+1), err, 1)
endif
*
if (myproc.eq.0) then
do k=1,gnk
do j=1,gnjy
do i=1,gnix
tp1(i,j,k) = 0. ! no humidity
end do
end do
end do
endif
call glbdist2 (tp1,hup ,minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk)
*
if (myproc.eq.0) then
do k=1,gnk
if (zm(k).gt.haut4) then
if (k.gt.1) then
tp1(1,1,k)=tp1(1,1,k-1)
tp2(1,1,k)=tp2(1,1,k-1)
else
print*,'zm(1) > haut4'
print*,'Initialisation du profil de vent avec ug et vg'
tp1(1,1,k)=uvg
tp2(1,1,k)=vvg
endif
else
call intcub_2pts_2d2('u',zm(k),tp1(1,1,k))
call intcub_2pts_2d2('v',zm(k),tp2(1,1,k))
endif
print*,'ug (',k,')=',tp1(1,1,k)
print*,'vg (',k,')=',tp2(1,1,k)
end do
*
do k=1,gnk
do j=1,gnjy
do i=1,gnix
tp1(i,j,k) = tp1(1,1,k)
tp2(i,j,k) = tp2(1,1,k)
enddo
enddo
enddo
endif
call glbdist2 (tp1,uup ,minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk)
call glbdist2 (tp2,vvp ,minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk)
*
* get the m2
call glbcolc (sby1,gnix,gnjy,sby,minx,maxx,miny,maxy,1) !verifier dims
*
tropo = 10000.d0
if (myproc.eq.0) then
do k=1,gnk
if (ztr(k).le.tropo) then
if (ztr(k).gt.haut4) then
c MDMD inliner intlin (de grace)
cccccccccccc call intlin(haut3,haut4,tprofil3,tprofil4,
cccccccccccc $ ztr(k),tp1(1,1,k))
else
call intcub_2pts_2d2('t',ztr(k),tp1(1,1,k))
endif
else
if (k.eq.1) then
tp1(1,1,k)=grtstar
else
tp1(1,1,k)=tp1(1,1,k-1)
endif
endif
enddo
endif
*
****
* Calcul des valeurs pleines : P=Pstar+Pprime
****
c*** use new variables Tprime and Qprime (beginning), Jan. 2003
c
my_psol = 100000.0d0
*
if (myproc.eq.0) then
c
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
hhm(i,j,0)=hh0i(i,j,1)
end do
end do
c
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
hhw(i,j,k)=h_geop(zt(k),i,j)
hhm(i,j,k)=h_geop(zm(k),i,j)
hht(i,j,k)=h_geop(ztr(k),i,j)
end do
end do
end do
c
call qntstar(qstr,nssq,gots,orts,hht,hhm,
& (maxx-minx+1)*(maxy-miny+1),0,gnk)
*
c1 = grav_8 / rgasd_8
c
qp1(1,1,0) = dlog(my_psol)
*
qp1(1,1,1) = qp1(1,1,0)-c1/tp1(1,1,1)*dble(zm(1)) !zm(0) n existe pas
*
print*,'k=',0,'t= p=',exp(qp1(1,1,0))
print*,'k=',1,'t=',tp1(1,1,1),'p=',exp(qp1(1,1,1))
*
do k=2,gnk
qp1(1,1,k)= qp1(1,1,k-1)-c1/tp1(1,1,k)*(zm(k)-zm(k-1))
print*,'k=',k,'t=',tp1(1,1,k),'p=',exp(qp1(1,1,k))
end do
*
****
* Calcul des perturbations
****
do k = 0, gnk
qp1(1,1,k) = ( qp1(1,1,k) - qstr(1,1,k) ) / orts(1,1,k)
end do
c
print*,'k=',0,'tprime= pprime=',qp1(1,1,0)
c
do k = 1, gnk
km1=max(k-1,1)
kp1=min(k+1,gnk)
beta = nssq(1,1,k) / grav_8 - gots(1,1,k) / cpd_8
dh = ( hhw(1,1,kp1) - hhw(1,1,km1) ) * pt5
dqdz = ( qp1(1,j,k) - qp1(1,1,k-1) ) / dh
& - ( qp1(1,1,k) + qp1(1,1,k-1) ) * pt5 * beta
tp1(1,1,k) = grav_8 * dqdz / ( grav_8 - dqdz )
c
print*,'k=',k,'tprime=',tp1(1,1,k),'qprime=',qp1(1,1,k),
& 'beta, dh, dqdz=', beta, dh, dqdz
end do
**********************************
*
****
* Champs constants de dq et T en (x,y)
****
do j=1,gnjy
do i=1,gnix
qp1 (i,j,0)= qp1(1,1,0)
end do
end do
*
do k=1,gnk
do j=1,gnjy
do i=1,gnix
tp1 (i,j,k)= tp1(1,1,k)
qp1 (i,j,k)= qp1(1,1,k)
end do
end do
end do
*
endif
*
call glbdist2 (tp1,bbp ,minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk)
call glbdist2 (qp1,ppp ,minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk+1)
*
* do the pressure balance for geostrophic flow (operates on ppp)
print *,'s/r eole_qbal called....'
*
call eole_qbal
*
*
* Vertical motion
*
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
wwp(i,j,k) = 0.
end do
end do
end do
*
c---- A Buble is introduced as a tracer
if (myproc.eq.0) then
kc=0
do k=1,gnk-1
if (zm(k).lt.blb_zp) kc=k
enddo
print *,"kc=",kc," : ",zm(kc)
do k=1,gnk
do j=1,gnjy
do i=1,gnix
tp1(i,j,k)=100./(1.+
$ (real(k-kc)/real(blb_zs))**4.+
$ (real(i-gnix/4)/real(blb_xs))**4.)
c tp1(i,j,k)= min(5,max(5-abs(i-gnix/4),0))
c tp1(i,j,k)= 1.
* tp1(i,j,k)= 0.
end do
end do
end do
endif
* attention ...
* le Bubble tracer doit aller en derniere place des traceurs
call glbdist2 (tp1,trp(minx,miny,1,ntr) ,
$ minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk)
c----
if (ctebcs) then
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppnta(i,j,k) = ppp(i,j,k)
uunta(i,j,k) = uup(i,j,k)
vvnta(i,j,k) = vvp(i,j,k)
wwnta(i,j,k) = wwp(i,j,k)
bbnta(i,j,k) = bbp(i,j,k)
hunta(i,j,k) = hup(i,j,k)
end do
end do
do n=1,ntr
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
trnta(i,j,k,n) = trp(i,j,k,n)
end do
end do
end do
end do
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppnta(i,j,0) = ppp(i,j,0)
end do
end do
endif
*
dim = ndynvar*dim3d+dim2d
pas1 = pappp
pad1 = papp0
do i=1,dim
d(i) = s(i)
end do
*
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppntt(i,j,k) = ppp(i,j,k)
uuntt(i,j,k) = uup(i,j,k)
vvntt(i,j,k) = vvp(i,j,k)
wwntt(i,j,k) = wwp(i,j,k)
bbntt(i,j,k) = bbp(i,j,k)
huntt(i,j,k) = hup(i,j,k)
end do
end do
end do
do n = 1, ntr
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
* trntt(i,j,k,n) = trp(i,j,k,n)
end do
end do
end do
end do
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppntt(i,j,0) = ppp(i,j,0)
end do
end do
*
if (myproc.eq.0) then
write (6,602)
call hpdeallc (patp1 ,err,1)
call hpdeallc (patp2 ,err,1)
call hpdeallc (paqp1 ,err,1)
endif
*
if ( fhalo .and. vraies_mtn.eq.2 ) then
call draglaw ( geobus, ngeop )
endif
c
602 format (/'INITIALIZATION OF EOLE PROBLEM COMPLETED')
*---------------------------------------------------------------------
return
end
*
subroutine eole_qbal 1,1
*
implicit none
*
#include "grd.cdk"
#include "consdyn_8.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
c#include "yomdyn2.cdk"
#include "dynmem.cdk"
#include "levels.cdk"
c***
#include "topo.cdk"
#include "h_geop.cdk"
#include "eole.cdk"
*
integer i, j, k, kk, km, kp, err,nb_instab
real*8 dxi, dyi, one, two, three, p25,
$ c2,c1, vbarxy, ubarxy, tpbarz, dz,
$ vr(minx:maxx,miny:maxy,1:nk),
$ ur(minx:maxx,miny:maxy,1:nk),
$ ppr(minx:maxx,miny:maxy)
*
parameter (p25=0.25, one = 1.0)
real ref_adiab, ref_conv, critere_stab, petit
pointer (pat_ad, ref_adiab(2:gnk)),
$ (pat_cv, ref_conv(2:gnk))
parameter (critere_stab=9.3) !refroidissement adiab : -9.8C/km
parameter (petit=0.0001)
c ***
integer km1, kp1
real*8 beta, dh, dqdz, pt5
real hht(minx:maxx,miny:maxy,gnk), hhm(minx:maxx,miny:maxy,0:gnk)
real hhw(minx:maxx,miny:maxy, gnk)
pt5 = 0.5d0
c ***
* ------------------------------------
*
dxi = one / dble(grdx)
dyi = one / dble(grdy)
*
****
* Calcul de grad (q) : ur=R(Tstar+Tprime)dq/dx ; vr=R(Tstar+Tprime)dq/dy
****
do k=1,gnk-1
do j = 1-hy, ldnj+hy-1
do i = 1-hx+1,ldni+hx
vbarxy = p25 * ( ( vvp(i-1,j+1,k)+vvp(i-1,j,k) ) +
$ ( vvp(i ,j+1,k)+vvp(i ,j,k) ) )
ur(i,j,k) = + vbarxy * ( fcor(i,j+1) + fcor(i,j) ) * pt5
!Coriolis doit rester constant
end do
end do
*
do j = 1-hy+1,ldnj+hy
do i = 1-hx,ldni+hx-1
ubarxy = p25 * ( (uup(i+1,j-1,k)+uup(i,j-1,k) )
$ + (uup(i+1,j ,k)+uup(i, j,k) ) )
vr(i,j,k) = - ubarxy * ( fcor(i+1,j) + fcor(i,j) ) * pt5
end do
end do
*
i=ldni+hx
do j = 1-hy+1,ldnj+hy
vr(i,j,k) = vr(i-1,j,k) !OK car grad q constant
enddo
enddo
*
****
* ppp by integration
* Tout doit etre constant dans l equilibre geostrophique
****
do k=0,gnk
* kk=min(gnk-1,max(1,k))
kk=min(gnk-1,k)
km=max(1,k)
kp=min(gnk,k+1)
tpbarz=(bbp(1,1,km)+bbp(1,1,kp))*pt5
c*** use new variable Qprime=R*Tstar*Qprime(old), (beginning) Jan. 2003
c1 = dble(grdx) / ( 1 + tpbarz / grav_8 )
c*** use new variable Qprime=R*Tstar*Qprime(old), (end) Jan. 2003
c2=grav_8/(rgasd_8*dble(grtstar))
*
* mid point
*
i=ldni/2
j=ldnj/2
ppr(i,j)=ppp(i,j,k) !Calculs en double precision
*
* mid row minus mid point
*
* A k=0, ur pas defini. Si defini, probleme dans runsor pour TT !!?
do i= ldni/2-1,1-hx,-1
if (k.eq.0) then
ppr(i,j)=ppr(i+1,j)-c1*vprofil1*fcor(i,j)
else
ppr(i,j)=ppr(i+1,j)-c1*ur(i+1,j,kk)
endif
enddo
*
do i= ldni/2+1,ldni+hx
if (k.eq.0) then
ppr(i,j)=ppr(i-1,j)+c1*vprofil1*fcor(i,j)
else
ppr(i,j)=ppr(i-1,j)+c1*ur(i,j,kk)
endif
enddo
*
* other rows
*
do j= ldnj/2-1,1-hy,-1
do i= 1-hx,ldni+hx
if (k.eq.0) then
ppr(i,j)=ppr(i,j+1)-c1*(-uprofil1*fcor(i,j))
else
ppr(i,j)=ppr(i,j+1)-c1*vr(i,j+1,kk)
endif
enddo
enddo
*
do j= ldnj/2+1,ldnj+hy
do i= 1-hx,ldni+hx
if (k.eq.0) then
ppr(i,j)=ppr(i,j-1)+c1*(-uprofil1*fcor(i,j))
else
ppr(i,j)=ppr(i,j-1)+c1*vr(i,j,kk)
endif
enddo
enddo
*
do j= 1-hy,ldnj+hy
do i= 1-hx,ldni+hx
ppp(i,j,k)=ppr(i,j)
c ppm(i,j,k)=ppp(i,j,k)
c pp0(i,j,k)=ppp(i,j,k)
enddo
enddo
c ***
print*, 'in eole_qbal: k, ppp = ', k, ppp(3,3,k)
enddo
*
****
* Iteration avec l'equation hydrostatique pour trouver
* un profil de temperature T(x,y,z) plus proche des 2 equilibres
****
c
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
hhm(i,j,0)=hh0i(i,j,1)
end do
end do
c
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
hhw(i,j,k)=h_geop(zt(k),i,j)
hhm(i,j,k)=h_geop(zm(k),i,j)
hht(i,j,k)=h_geop(ztr(k),i,j)
end do
end do
end do
c
call qntstar(qstr,nssq,gots,orts,hht,hhm,
& (maxx-minx+1)*(maxy-miny+1),0,gnk)
c
do k=1,gnk
dz = zm(1)
if(k.ne.1) dz = zm(k) - zm(k-1)
km1=max(k-1,1)
kp1=min(k+1,gnk)
do j= 1-hy,ldnj+hy
do i= 1-hx,ldni+hx
beta = nssq(i,j,k) / grav_8 - gots(i,j,k) / cpd_8
dh = ( hhw(i,j,kp1) - hhw(i,j,km1) ) * pt5
dqdz = ( ppp(i,j,k) - ppp(i,j,k-1) ) / dh
& - ( ppp(i,j,k) + ppp(i,j,k-1) ) * pt5 * beta
bbp(i,j,k)= grav_8 * dqdz / ( grav_8 - dqdz )
c bbm(i,j,k)= bbp(i,j,k)
c bb0(i,j,k)= bbp(i,j,k)
enddo
enddo
print*,'tprime_recalculee(1,1',k,'=)',bbp(1,1,k),
& 'beta, dh, dqdz=', beta, dh, dqdz
enddo
*
if (stabilite_air.eq.1) then
* Attention, les 2 equilibres sont moins bons
nb_instab=0
call hpalloc (pat_ad, gnk-1, err, 1)
call hpalloc (pat_cv, gnk-1, err, 1)
do k=2,gnk
ref_adiab(k)=-critere_stab*(ztr(k)-ztr(k-1))/1000.
enddo
do j= 1-hy,ldnj+hy
do i= 1-hx,ldni+hx
100 continue
do k=2,gnk
ref_conv(k)=(bbp(i,j,k)-bbp(i,j,k-1))
ref_conv(k)=ref_conv(k)-ref_adiab(k)
enddo
*
k=2
110 if (ref_conv(k).lt.(0.-petit)) then
nb_instab=nb_instab+1
bbp(i,j,k)=bbp(i,j,k-1)+ref_adiab(k)
bbp(i,j,k+1)=bbp(i,j,k+1)+ref_conv(k)
print*,' ! Profil T(',i,',',j,',',k,') -
$ T(',i,',',j,',',k-1,') instable '
goto 100
else
k=k+1
if (k.eq.gnk) then !pas d instabilite en gnk
goto 200 !car profil ct de temperature
else
goto 110
endif
endif
200 continue
enddo
enddo
*
call hpdeallc (pat_ad, err, 1)
call hpdeallc (pat_cv, err, 1)
print*,' '
print*,'----------------------------------------------'
print*,'Nombre d instabilites corrigees: ',nb_instab
print*,'----------------------------------------------'
endif
*
*----------------------------------------------------------------------
*
return
end
*
copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
*
***s/r draglaw (geobus,ngeop)
*
subroutine draglaw (geobus,ngeop) 1
implicit none
*
*AUTHOR David Lemarquis Sep 2001
*
*REVISION
*
* Wei Yu (Jan. 2003)
* Use Z0 instead of ZP for surface roughness
*
*IMPLICIT
#include "eole.cdk"
#include "lesbus.cdk"
#include "dynmem.cdk"
#include "consdyn_8.cdk"
#include "yomdyn1.cdk"
#include "levels.cdk"
#include "grd.cdk"
*
* integer i,j,k,ldni,ldnj,ni,nj, nk,k_hcl
integer i, j, k, ni, nj, k_hcl
integer offj,izp,err
integer ngeop
real geobus(ngeop)
* real zzero,ue,tau0,tau,taux,tauy
* pointer (pat_z0, zzero(minx:maxx,miny:maxy)),
* $ (pat_ue, ue(minx:maxx,miny:maxy)),
* $ (pat_tau0, tau0(minx:maxx,miny:maxy)),
* $ (pat_tau, tau(minx:maxx,miny:maxy,gnk)),
* $ (pat_taux, taux(minx:maxx,miny:maxy,gnk)),
* $ (pat_tauy, tauy(minx:maxx,miny:maxy,gnk))
real zzero(minx:maxx,miny:maxy),
$ ue(minx:maxx,miny:maxy),
$ tau0(minx:maxx,miny:maxy)
* $ tau(minx:maxx,miny:maxy,gnk),
* $ taux(minx:maxx,miny:maxy,gnk),
* $ tauy(minx:maxx,miny:maxy,gnk)
real uemin,uemax,eps,ueg,ued,ueopp,deltaue
integer nit,nitmax,signeyg,signeyd,signey
real yg,yd,y,A,B,hcl,cdd,sdd,rho0,tbarz
real *8 ug1,vg1,geostr,costh,sinth,uv1,u1,v1,bx,by,bb,dens,G
logical verb
data verb/.false./
*
* profil de gradient de stress=gprime
*
real gprime,z,h
gprime(z,h)=-2*(1-z/h)/h
*
A=0.9 !These Robert Benoit mars 1073 p21
B=4.5
* uemin=0.0001 !Hauteur couche lim = ue/f
uemin=0.000001 !gaspe testing
* uemax=1. !ue moyen = 0.07
uemax=100. ! uemax=1 maybe too low (gaspesie)
nitmax=100
eps=0.001 !eps=delta ue / ue
G=sqrt(dble(uprofil1**2.)+dble(vprofil1**2.))
hcl=1000. !hauteur de la couche limite (m)
*
print *,'draglaw.gni=',gni,' gnj=',gnj,' ldni,ldnj=',ldni,ldnj,
$ ' ni,nj=',ni,nj
* ldni=gni
* ldnj=gnj
ni=gni+2*hx
nj=gnj+2*hy
* nk=gnk
print *,'draglaw.ni=',ni,' nj=',nj,' gni=',gni,' gnj=',gnj,
$ ' minx=',minx,' maxx=',maxx,' miny=',miny,' maxy=',maxy,
$ 'gnk,nk=',gnk,nk
* stop
*
do k=1,nk
if (ztr(k).gt.hcl) then
goto 10
else
k_hcl=k+1
endif
enddo
*
10 continue
* 10 call hpalloc (pat_z0, ni*nj, err, 1)
* call hpalloc (pat_ue, ni*nj, err, 1)
* call hpalloc (pat_tau0, ni*nj, err, 1)
* call hpalloc (pat_tau, ni*nj*nk, err, 1)
* call hpalloc (pat_taux, ni*nj*nk, err, 1)
* call hpalloc (pat_tauy, ni*nj*nk, err, 1)
*
do i=1,geotop
c Line changed (WY, Jan. 2003)
c if (geonm(i,2).eq.'ZP') izp=i
if (geonm(i,2).eq.'Z0') izp=i
c-----------
enddo
* Lecture de la longueur de rugosite au centre du domaine
do j= 1,ldnj
offj=ldni*(j-1)
do i= 1,ldni
zzero(i,j)=geobus(geopar(izp,1)+offj+i-1)
enddo
enddo
*
* Recopiage de Z0 sur les bandes
do i=minx,0
do j=1,ldnj
zzero(i,j)=zzero(1,j)
enddo
enddo
*
do i=ldni+1,maxx
do j=1,ldnj
zzero(i,j)=zzero(ldni,j)
enddo
enddo
*
do i=1,ldni
do j=miny,0
zzero(i,j)=zzero(i,1)
enddo
enddo
*
do i=1,ldni
do j=ldnj+1,maxy
zzero(i,j)=zzero(i,ldnj)
enddo
enddo
*
* Recopiage de Z0 dans les coins
do i=minx,0
do j=miny,0
zzero(i,j)=zzero(1,1)
enddo
enddo
*
do i=minx,0
do j=ldnj+1,maxy
zzero(i,j)=zzero(1,ldnj)
enddo
enddo
*
do i=ldni+1,maxx
do j=miny,0
zzero(i,j)=zzero(ldni,1)
enddo
enddo
*
do i=ldni+1,maxx
do j=ldnj+1,maxy
zzero(i,j)=zzero(ldni,ldnj)
enddo
enddo
*
* Calcul de Uetoile, le stress de surface par dichotomie
* do i=minx,maxx
do i=minx,maxx-1
do j=miny,maxy-1
* do j=miny,maxy
nit=0
ueg=uemin
ued=uemax
ue(i,j)=(ueg+ued)/2.
* print *,'fcor,i,j=',fcor(i,j),i,j
*
50 nit=nit+1
if (nit.gt.nitmax) then
print*,'Calcul du stress de surface impossible'
print*,'Nombre max d iteration depasse'
stop
endif
yg=ueg/karman_8 * sqrt((alog(ueg/(fcor(i,j)*
$ zzero(i,j))) - A)**2. + B**2) - G
if (yg.lt.0.) then
signeyg=-1
else
signeyg=+1
endif
yd=ued/karman_8 * sqrt((alog(ued/(fcor(i,j)*
$ zzero(i,j))) - A)**2. + B**2) - G
if (yd.lt.0.) then
signeyd=-1
else
signeyd=+1
endif
if (signeyg.eq.signeyd) then
print*,'Calcul du stress de surface impossible'
print *,nit,ueg,ued,i,j,zzero(i,j),g,fcor(i,j),
$ karman_8
stop
endif
*
y=ue(i,j)/karman_8 * sqrt((alog(ue(i,j)/(fcor(i,j)*
$ zzero(i,j))) - A)**2. + B**2) - G
*
if (y.lt.0.) then
signey=-1
else
signey=+1
endif
*
if (signey.ne.signeyg) then
ued=ue(i,j)
ueopp=ueg
else
ueg=ue(i,j)
ueopp=ued
endif
ue(i,j)=(ue(i,j)+ueopp)/2.
*
deltaue=sqrt((ue(i,j)-ueopp)**2.)
if ((deltaue/ue(i,j)).lt.eps) then
goto 100
else
goto 50
endif
*
100 continue
*
* Calcul de Tau0
c rho0=my_psol/(rgasd_8*(bbp(i,j,1)+grtstar))
rho0 = my_psol / (rgasd_8*(1.+bbp(i,j,1)/grav_8)*grtstar)
tau0(i,j)=ue(i,j)**2. * rho0
if (verb)
$ print*,'draglaw...',i,j,bbp(i,j,1)+grtstar,tau0(i,j),zzero(i,j),
$ ue(i,j)
*
* Profil quadratique de tau(z)
* do k=1,k_hcl
* tau(i,j,k)=tau0(i,j) * (ztr(k)**2./hcl**2. -
* $ 2.*ztr(k)/hcl + 1.)
* cdd=uup(i,j,k)/(sqrt(uup(i,j,k)**2.+vvp(i,j,k)**2.))
* sdd=vvp(i,j,k)/(sqrt(uup(i,j,k)**2.+vvp(i,j,k)**2.))
* taux(i,j,k)=tau(i,j,k)*(-cdd)
* tauy(i,j,k)=tau(i,j,k)*(-sdd)
* enddo
*
enddo
enddo
*
* Calcul des nouveaux vents
* a l'entree, uup,vvp contiennent ug,vg (geostrophique)
* premier niveau (special, donne direction profil de stress)
* do j=miny,maxy
* do i=minx,maxx
do j=miny,maxy-1
do i=minx,maxx-1
k=1
ug1=uup(i,j,k)
vg1=vvp(i,j,k)
bx=-fcor(i,j)*vg1
by=+fcor(i,j)*ug1
geostr=sqrt(ug1**2+vg1**2)
bb=sqrt(bx*bx+by*by)
costh=ue(i,j)**2*dble(abs(gprime(ztr(k),hcl)/fcor(i,j)))/geostr
if (costh.gt.0.99) then
print *,'draglaw...costh,i,j,=',costh,i,j
costh=0.99 !avoid 1.0000x value (large z0?)
print *,'draglaw...costh____,i,j,=',costh,i,j
endif
sinth=sqrt(dble(1.0)-costh**2)
uv1=geostr*sinth
u1=uv1*(costh*bx+sinth*by)/bb
v1=uv1*(costh*by-sinth*bx)/bb
uup(i,j,k)=u1
vvp(i,j,k)=v1
uu0(i,j,k) = uup(i,j,k)
vv0(i,j,k) = vvp(i,j,k)
* nest aussi !
uunta(i,j,k) = uup(i,j,k)
vvnta(i,j,k) = vvp(i,j,k)
uuntt(i,j,k) = uup(i,j,k)
vvntt(i,j,k) = vvp(i,j,k)
if (verb)
$ print*,'draglaw...uv1',
$ i,j,ug1,vg1,bx,by,geostr,bb,costh,sinth,uv1,u1,v1,
$ (ug1*u1+vg1*v1)/(geostr*uv1),
$ hcl,ztr(k),gprime(ztr(k),hcl),zzero(i,j),fcor(i,j)
do k=2,k_hcl-1
ug1=uup(i,j,k)
vg1=vvp(i,j,k)
* tbarz=(bbp(i,j,k)+bbp(i,j,k+1))/2.
c dens=1. ! a ajuster ....
dens = my_psol/(rgasd_8*(1.+bbp(i,j,k)/grav_8)*grtstar)
vvp(i,j,k)=vg1
$ -tau0(i,j)*gprime(ztr(k),hcl)*u1/uv1/fcor(i,j)/dens
uup(i,j,k)=ug1
$ +tau0(i,j)*gprime(ztr(k),hcl)*v1/uv1/fcor(i,j)/dens
* nest aussi !
uu0(i,j,k) = uup(i,j,k)
vv0(i,j,k) = vvp(i,j,k)
uunta(i,j,k) = uup(i,j,k)
vvnta(i,j,k) = vvp(i,j,k)
uuntt(i,j,k) = uup(i,j,k)
vvntt(i,j,k) = vvp(i,j,k)
geostr=sqrt(ug1**2+vg1**2)
u1=uup(i,j,k)
v1=vvp(i,j,k)
uv1=sqrt(u1**2+v1**2)
if (verb)
$ print*,'draglaw...uvk',k,ug1,vg1,uup(i,j,k),vvp(i,j,k),
$ geostr,uv1,(ug1*u1+vg1*v1)/(geostr*uv1)
enddo
enddo
enddo
* stop
*
* call hpdeallc (pat_z0, err, 1)
* call hpdeallc (pat_ue, err, 1)
* call hpdeallc (pat_tau0, err, 1)
* call hpdeallc (pat_tau, err, 1)
* call hpdeallc (pat_taux, err, 1)
* call hpdeallc (pat_tauy, err, 1)
*
print *,'draglaw...returning'
return
end
*
***s/r intcub_2pts_2d2
*
subroutine intcub_2pts_2d2(choix,z,f) 3
implicit none
*
#include "eole.cdk"
*
real f1,f2,f3,f4,x1,x2,x3,x4
real f4p,f1pp,f2pp,f3pp
real z,f
real a,b,c,d
character choix
*
*
x1=haut1
x2=haut2
x3=haut3
x4=haut4
*
if (choix.eq.'u') then
f1=uprofil1
f2=uprofil2
f3=uprofil3
f4=uprofil4
endif
if (choix.eq.'v') then
f1=vprofil1
f2=vprofil2
f3=vprofil3
f4=vprofil4
endif
if (choix.eq.'t') then
f1=tprofil1
f2=tprofil2
f3=tprofil3
f4=tprofil4
endif
*
f1pp=0.
f2pp=( ((f3-f2)/(x3-x2)) - ((f2-f1)/(x2-x1)) ) /
$ (0.5*(x3-x1))
f3pp=( ((f4-f3)/(x4-x3)) - ((f3-f2)/(x3-x2)) ) /
$ (0.5*(x4-x2))
if (choix.eq.'t') then
f4p=(f4-f3)/(x4-x3)
else
f4p=0.
endif
*
*
if ((z.gt.x4).or.(z.lt.x1)) then
print*,'intcub_2pts_2d2 : altitude hors du domaine'
stop
endif
*
if (z.le.x2) then
d=(f2pp-f1pp)/(6.*(x2-x1))
c=(f1pp-6.*d*x1)/2.
b=(f2-f1)/(x2-x1) - c*(x2+x1) - d*(x1**2.+x1*x2+x2**2.)
a=f1-b*x1-c*x1**2.-d*x1**3.
endif
*
if ((z.gt.x2).and.(z.le.x3)) then
d=(f3pp-f2pp)/(6.*(x3-x2))
c=(f2pp-6.*d*x2)/2.
b=(f3-f2)/(x3-x2) - c*(x3+x2) - d*(x2**2.+x2*x3+x3**2.)
a=f2-b*x2-c*x2**2.-d*x2**3.
endif
*
if ((z.gt.x3).and.(z.le.x4)) then
c=1./(-x3+x4+(x3**2.+x3*x4-2.*x4**2.)/(3.*x3)) *
$ (f4p-(f4-f3)/(x4-x3) + f3pp/(6.*x3) *
$ (x3**2.+x3*x4-2.*x4**2.))
d=(f3pp-2.*c)/(6.*x3)
b=f4p-2.*c*x4-(f3pp-2.*c)/(2.*x3) * x4**2.
a=f3-b*x3-c*x3**2.-d*x3**3.
endif
*
*
f=a+b*z+c*z**2.+d*z**3.
*
c***debug
if ( choix .eq. 't' ) then
write(78,*) 'x1, x2, x3, x4, z=', x1, x2, x3, x4, z
write(78,*) 'f1, f2, f3, f4, f=', f1, f2, f3, f4, f
endif
c***debug
*
*
return
end