program mount_wave
!**************************
!  0.  COMPILATION 
!
! make
!****************************
!  1.  Declarations
!****************************

!use nr
!Prantl, Richardon number, small parameter, boundary layer depth:
use viscous_m, only: prt=>prk,zj=>jrk, epd=>eps, sld=>bly
use viscous_m, only: viscous,vis1,vis2,vis3,vis4,hankel,hanksol,bab1,bab2,hypsol
use nr, only: hypgeo
implicit none

!  Versions of the programme:
logical :: lout  ; integer :: iout! Control outputs
!  Integer controlling dynamics and boundary condition:
integer :: icond    ! 0 for inviscid, 1 for no slip, 2 for now heat flux
                    ! 3 For free slip (not working well)
character(len=3), dimension(0:3) :: cond=(/'INV','NSL','NFL','FSL'/)
integer iric           !  Control the choice of the Richardon number
real, allocatable :: ric(:)
!real, dimension(9) :: ric=(/0.01,0.5,1.,2.,4.,9.,16.,25.,36./)
!  Looping case:
!real, dimension(40) :: ric=(/0.01,0.02,0.025,0.050,0.075,0.1,0.15,0.175,&
!                           0.2,0.225,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.2, &
!                           1.4,1.6,1.8,2.,2.5,3.,4.,5.,7.5,9., &
!                           11.,13.,16.,19.,22.,25.,28.,31.,34.,36.,40./)
integer islope       ! Control the choice of the Slope (zrid)
real, dimension(9) :: slope=(/0.01,0.02,0.05,0.075,0.1,0.125,0.15,0.175,0.185/)
integer isld         ! Control the shear layer depth sld
real, dimension(7) :: depth=(/0.1,0.25,0.5,1.,2.,4.,10./)
integer :: ibouss                         ! 1 for Boussinesq 0 for hydro 
character(len=4) :: exptype='BXXXX'       ! Nom de l'expérience
character(len=8) :: exppara='JX_SX_DV'    ! Nom de l'expérience
integer, parameter :: iim=512             ! Horizontal grid size
real, parameter ::  xlen=100., zlen=10.    ! Longueur et hauteur du domain.
!integer, parameter :: iim=256         ! Horizontal grid size
!real, parameter ::  xlen=100., zlen=2.    ! Longueur et hauteur du domain.
!debug values:
!integer, parameter :: iim=512            ! Horizontal grid size
!real, parameter ::  xlen=25., zlen=3.0   ! Longueur et hauteur du domain.
integer :: llm                            ! Vertical grid size
real ::             zrid                  ! Hauteur de l'obstacle
real ::  zbl=0.1   !  Boundary layer depth

complex :: cgamma
complex, parameter :: ci=(0.,1.),c1=(1.,0.)
real ::    pi=3.1415927
character(len=36) ::  filename
integer :: ii,jj,ll,kk,jm,jp                       ! looping integers
logical :: logic,lfirst                      ! tests

! Inviscid solution:
complex, dimension(iim) :: aa,bb  !  aa, bb entering in hyp. functions
complex :: cc                     !  cc entering in the hyp. functions
complex :: cwp,cwm,cwdet          !  Temporary arrays                  
real ::    zm
complex :: zmu                !  Vertical wavenumber and sqrt(J-0.25)
complex, dimension(iim) :: oa1,oa2 ! ai's on outer variables
complex, dimension(iim) :: ta1,ta2,ta3,ta4 ! ai's in inner variables
complex, dimension(iim) :: tb3,tb4,tc3,tc4 ! bi's and ci's in inner variables
complex, dimension(iim) :: a1,a2,a3,a4
complex, dimension(iim) :: b1,b2,det,cm ! Matching structure, vertical wavenumber
complex, dimension(3) :: inv

! Grille et montagne
real :: dx,dz,dk,xkmax,dzmin,dzmax       ! Resolutions domaines physique et spectral
integer :: lmax
real, dimension(iim) :: xx,xk ! Grille et nombre d'onde horizontaux
real, dimension(iim) :: del,tiz  ! Haut couche limite, altitude en var. intérieure.
real, dimension(iim) :: zh,zhx       !  Hauteur montagne et sa derive
real, dimension(iim) :: psol,fxsol,fzsol !  Surface pressure and action fluxes
real, dimension(iim) :: usol,wsol,bsol,rsol,uzsol   !  Surface winds and co
real, dimension(iim) :: psth,usth,wsth,bsth,rsth,uzsth,bzsth  ! Theoric surface values
real ::     zth      ! Hauteur montagne en variable interieure
complex, dimension(iim) ::  cfa,cfb,cfc,cfd  ! Fonctions f(k)
real, dimension(12) :: visc
real, allocatable :: zz(:),rr(:)          ! Vertical grid
real, allocatable :: u0(:),u0z(:),u0zz(:),psi0(:) ! Background wind and shears

! Disturbance fields
real, allocatable :: vitw(:,:),vitu(:,:),buoy(:,:),pres(:,:),vort(:,:),shea(:,:),disp(:,:)
real, allocatable :: wwav(:,:),uwav(:,:),bwav(:,:),pwav(:,:),vwav(:,:),swav(:,:),dwav(:,:)
real, allocatable :: acti(:,:),fx(:,:),fz(:,:),psi(:,:),pdtr(:,:)
real, allocatable :: gwstr(:),gwpre(:),gwvis(:)     !  Gravity wave stress:


!  Complex numbers to recover JAS07

complex :: ws1,us1,bs1,ws2,us2,bs2,ws3,us3,bs3,ws4,us4,bs4
complex :: sde,cre,bb1,dd1,bb2,dd2,bb3,dd3,bb4,dd4
complex, dimension(6) :: cv1,cv2,cv3,cv4,cv12
complex, allocatable :: cs1(:,:),cs2(:,:),cs3(:,:),cs4(:,:)
complex, allocatable :: csa(:,:),csb(:,:),csc(:,:),csd(:,:)

!  Complex matrices needed to solve the problem

complex, dimension(iim,iim) :: cw1,cu1,cb1,cw2,cu2,cb2,cw3,cu3,cb3,cw4,cu4,cb4
complex, dimension(iim,iim) :: cc1,cc2,cc3,cc4                                 
complex, dimension(iim,iim) :: cma,cmb,cmc,cmd,cme,cmf,cmh,cm1,cm2,cm3,cm4,cm5
complex, dimension(iim,iim) :: cwa,cua,cba,cwb,cub,cbb,cwc,cuc,cbc,cwd,cud,cbd
complex, dimension(iim,iim) :: cca,ccb,ccc  ! ch=u_z

! Complex vectors to represent the solutions:
complex, dimension(iim) :: cwvis,cwinv,cwvis1,cwvis2,cwvis3,cwvis4,cwmat,cwuni
complex, dimension(iim) :: cuvis,cuinv,cuvis1,cuvis2,cuvis3,cuvis4,cumat,cuuni
complex, dimension(iim) :: cbvis,cbinv,cbvis1,cbvis2,cbvis3,cbvis4,cbmat,cbuni
complex, dimension(iim) :: ccvis,ccinv,ccvis1,ccvis2,ccvis3,ccvis4,ccmat,ccuni
complex, dimension(iim) :: cxvis,cxinv,cxvis1,cxvis2,cxvis3,cxvis4,cxmat,cxuni
complex, dimension(iim) :: cruni,cduni  !  Rotationnal, displacement
complex, dimension(iim) :: cpvis,cpinv,cpvis1,cpvis2,cpvis3,cpvis4,cpmat,cpuni
complex, dimension(iim) :: cwvisa,cwvisb,cwvisc,cwvisd
complex, dimension(iim) :: cuvisa,cuvisb,cuvisc,cuvisd
complex, dimension(iim) :: cbvisa,cbvisb,cbvisc,cbvisd
complex, dimension(iim) :: cpvisa,cpvisb,cpvisc,cpvisd
complex, dimension(iim) :: ccvisa,ccvisb,ccvisc,ccvisd
complex, dimension(iim) :: cxvisa,cxvisb,cxvisc,cxvisd
real, dimension(12) :: visc1,visc2,visc3,visc4
real, allocatable :: vsols(:,:),aruni(:),arvis(:),armat(:),arinv(:)
!
!  Tableaux pour les programmes d'inversion de matrice, FFT, ect
!
complex, dimension(iim,iim) :: cdet,cdet2
complex, dimension(2) :: fdet
complex, dimension(iim)  ::  cvec,cinv 
integer :: ipvt(iim)     !  Integer array needed to inverse CWJK
integer :: info            ! Measure if the inversion is feasible
real :: rcond               ! an estimate of the reciprocal condition of  cdet
real, dimension(2*iim) :: w1,w2   !  Work array for the fft
integer, dimension(15) :: ifac    !  Work array for the fft

! Cosmetics prints, synthetic diagnostics: 
real :: zmin,zmax,vdrag,storm,storp,down,tdrag,tgwdr,actfl,mdrag,gwmin,gwmax
real :: adtr,bdtr

!**************************
! 2. Initialisations
!**************************

! 2.1 Constant
!***************

epd=0.001              ! smallest grid distance
prt=2.0                ! Prandtl number

! 2.2 Interative set-up:
!**********************************
write (*,*)'Enter 1 for extensive output (non-looping option)';lout=.false.
write(*,*)'If first trial type 1'
read(*,*)iout ; if(iout.eq.1)lout=.true.

if(lout)then;allocate(ric(0:9));ric=(/0.01,0.1,0.2,0.3,0.7,1.,1.7,4.,9./)
!if(lout)then;allocate(ric(0:9));ric=(/0.15,0.20,0.24,0.26,0.3,0.35,0.45,0.5,1.,2./)
else;allocate(ric(40))
 ric=(/0.01,0.02,0.025,0.050,0.075,0.1,0.15,0.175,&
       0.2,0.225,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.2, &
       1.4,1.6,1.8,2.,2.5,3.,4.,5.,7.5,9., &
       11.,13.,16.,19.,22.,25.,28.,31.,34.,36.,40./)
endif

if(lout) write(*,*) '0: Hydrostatic constant shear (Part I)'
if(lout) write(*,*) '1: Non-hydro constant shear (Part II)'
if(lout) write(*,*) '2: Non-hydro variable shear (Part III)'
if(lout) write(*,*) '3: Hydrostatic variable shear (Part III)'
if(lout) write (*,*)'Enter 0 for HCSH, 1 for BCSH, 2 for BVSH, 3 for HVSH:'
read(*,*)ibouss

if(lout)write (*,*)' For no slip type 1'
if(lout)write (*,167)'Bound    :',(cond(ii),ii=0,3);167 format(a11,4(1x,a3,1x))
if(lout)write (*,168)'For icond:',(ii,ii=0,3);168 format(a11,9(2x,i2,1x))
if(lout)print *, 'Bound Cond, enter icond:'
read(*,*)icond
if(lout)print *,'Boundary condition:',cond(icond)
if(icond.eq.0) then; print *,'Inviscid case not clean yet!';stop; endif

if(lout)write (*,165)'  J     =',(ric(ii),ii=0,9);165 format(a9,10(1x,f5.2))
if(lout)write (*,166)'For iric=',(ii,ii=0,9);166 format(a9,10(2x,i2,2x))
if(lout)print *, 'Choice of J, enter iric:'
read(*,*)iric;if(iout.eq.0)iric=iric-10
zj=ric(iric); if(lout)print *,'Richardson zj=',zj

if(lout)write (*,169)'   Slope  =',(slope(ii),ii=1,9);169 format(a11,9(1x,f5.3))
if(lout)write (*,170)'For islope=',(ii,ii=1,9);170 format(a11,9(2x,i2,2x))
if(lout)print *, 'Slope enter islope:'
read(*,*)islope
zrid=slope(islope); if(lout)print *,'Slope zrid=',zrid

if(lout)write (*,171)'Shear layer depth D:',(depth(ii),ii=1,7);171 format(a20,7(1x,f6.3))
if(lout)write (*,172)'          For isld=:',(ii,ii=1,7);172 format(a20,7(3x,i1,3X))
if(lout)print *, 'Depth enter isld:' !; isld=4
read(*,*)isld;
sld=depth(isld); if(lout)print *,'Depth=',sld

if(ibouss.ne.2.and.isld.ne.4)then
write(*,*)'Pas de isld variable en hydro ou en cte shear';stop;endif
if(ibouss.ne.0.and.ibouss.ne.1.and.ibouss.ne.2.and.ibouss.ne.3) then
write(*,*)' problem with ibouss';stop;endif

! Change title according to the setup:
 
write(exptype(2:4),186)cond(icond);186 format(a3)
if(ibouss.eq.0)write(exptype(1:1),193)'H';193 format(a1)
if(ibouss.eq.1)write(exptype(1:1),193)'B'
if(ibouss.eq.2)write(exptype(1:1),193)'V'
if(ibouss.eq.3)write(exptype(1:1),193)'W'
if(lout)print *,'exptype=',exptype
if(iric.lt.10) then;
write(exppara(2:2),187)iric;write(exppara(5:5),187)islope
else
write(exppara(2:2),187)0;write(exppara(5:5),187)islope
endif;187 format(i1)
write(exppara(8:8),187)isld
if(lout)print *,'exppara=',exppara

if(icond.eq.0)zbl=zbl  ! A value for invisci stuff to behave as viscous (on work)
if(zj.gt.100.)then; print*,'Hankel not good for zj>100'; stop;endif

!  2.2)  Grid construction orography
!************************************

if(zj.ge.0.25) then;zmu=sqrt(abs(zj-0.25));else;zmu=ci*sqrt(abs(zj-0.25));endif
dx=xlen/float(iim);dk=2.*pi/xlen
do ii=1,iim 
xx(ii)=(ii-float(iim+1)/2.)*dx
zh(ii)=zrid/(1.+xx(ii)**2/2.)
zhx(ii)=-xx(ii)*zrid/(1.+xx(ii)**2/2.)**2
!Gaussian
zh(ii)=zrid*exp(-xx(ii)**2/2.)
zhx(ii)=-xx(ii)*zrid*exp(-xx(ii)**2/2.)
!Twin peacks
!zh(ii)=zrid*exp(-xx(ii)**2/2.+1)*xx(ii)**2./2.
!zhx(ii)=zrid*exp(-xx(ii)**2/2.+1)*(xx(ii)-xx(ii)**3./2.)
enddo
do ii=1,iim/2+1;xk(ii)=(ii-1)*dk;enddo;xk(1)=epd*dk
do ii=iim/2+2,iim;xk(ii)=-xk(iim+2-ii);enddo



if(icond.gt.0)then
  do ii=2,iim; del(ii)=zbl/(abs(xk(ii)))**(1./3.);enddo;del(1)=del(2)
else
  do ii=2,iim; del(ii)=(zbl/xk(ii)+zbl*xk(ii))/2.; enddo;del(1)=del(2)
endif
  zm=100.; do jj=1,iim; zm=amin1(zm,del(jj)); enddo
if(lout) print *,'min del',zm,'h over del:',zrid/zm
  zm=0.; do jj=1,iim; zm=amax1(zm,del(jj)); enddo
if(lout) then;print *,'max del',zm;print *,xk(10),del(10);endif

! Vertical grid with zoom at low level, determination of llm and allocations

!dzmax=zlen/100.; dzmin=min(zrid/10.,dzmax); zmin =0.; lmax=0
!dzmax=zlen/500.; dzmin=dzmax  ; zmin =0.; lmax=0
dzmax=zlen/100.; dzmin=zrid/10.  ; zmin =0.; lmax=0
dzmax=zlen/500.; dzmin=min(dzmax,zrid/10.)     ; zmin =0.; lmax=0
dzmax=zlen/500.; dzmin=dzmax    ; zmin =0.; lmax=0
do ll=2,1000;dz=(1.-tanh((zmin/zrid-20.)/10.))/2.*dzmin &
               +(1.+tanh((zmin/zrid-20.)/10.))/2.*dzmax ; zmax =zmin+dz 
if(zmax.ge.zlen.and.zmin.lt.zlen) then
if(lout)print *,'llm should equal',ll,zmin;lmax=ll;endif;zmin=zmax;enddo

if(lmax.eq.0) then ; print *,'llm should  be greater than 1000!'; stop;endif
llm=lmax

! Array allocation once vertical  dimension is evaluated
!*******************************************************

allocate(zz(llm),rr(llm),u0(llm),u0z(llm),u0zz(llm),psi0(llm))                      
zz(1)=0.;
do ll=2,llm;dz=(1.-tanh((zz(ll-1)/zrid-20.)/10.))/2.*dzmin&
              +(1.+tanh((zz(ll-1)/zrid-20.)/10.))/2.*dzmax            
!zz(ll)=zz(ll-1)+dz;write(134,*),dz,(zz(ll)+zz(ll-1))/2.;enddo
zz(ll)=zz(ll-1)+dz;enddo; rr(:)=tanh(zz(:)/sld)**2
if(ibouss.ge.2) then
u0(:)=sld*tanh(zz(:)/sld);u0z(:)=1.-rr(:);u0zz(:)=-2./sld**2*u0(:)*u0z(:)
psi0(:)=sld**2*log(cosh(zz(:)/sld))
else;u0(:)=zz(:);u0z(:)=1.;u0zz(:)=0.;psi0(:)=zz(:)**2/2.;endif


allocate(vitw(iim,llm),vitu(iim,llm),buoy(iim,llm),wwav(iim,llm),dwav(iim,llm))
allocate(uwav(iim,llm),bwav(iim,llm),pres(iim,llm),psi(iim,llm),vort(iim,llm))
allocate(pwav(iim,llm),swav(iim,llm),vwav(iim,llm))
allocate(acti(iim,llm),fx(iim,llm),fz(iim,llm),shea(iim,llm),pdtr(iim,llm))
allocate(cs1(6,llm),cs2(6,llm),cs3(6,llm),cs4(6,llm))
allocate(csa(6,llm),csb(6,llm),csc(6,llm),csd(6,llm))
allocate(vsols(12,llm),gwstr(llm),gwpre(llm),gwvis(llm))
allocate(aruni(llm),arvis(llm),armat(llm),arinv(llm))
!***************************************
! 3. Inviscid and viscous solutions used
!***************************************
! 3.1) a1 and a2 associated with Hankel or Hypergeometric functions:
!*******************************************************************

if(ibouss.eq.0) then ! Hydrostatic:
  if(zj.le.0.25) then;write(*,*)'No Hydrostatic when J<0.25'; stop; endif
  do ii=2,iim/2+1;oa2(ii)=1.;oa1(ii)=0.;ta1(ii)=0.
  ta2(ii)=oa2(ii)/xk(ii)*del(ii)**(-0.5+ci*zmu);enddo
endif
if(ibouss.eq.1)then  ! Boussinesq
  do 3 ii=2,iim/2+1
  oa1(ii)=-ci*sqrt(pi)/sinh(zmu*pi)/cgamma(1.-ci*zmu)*(xk(ii)/2.)**(0.5-ci*zmu)
  oa2(ii)=conjg(oa1(ii))
  ta1(ii)=oa1(ii)/xk(ii)*del(ii)**(-0.5-ci*zmu)
  ta2(ii)=+conjg(ta1(ii))
  if(zj.lt.0.25)then
  oa2(ii)=+ci*sqrt(pi)/sinh(zmu*pi)/cgamma(1.+ci*zmu)*(xk(ii)/2.)**(0.5+ci*zmu)
  ta2(ii)=oa2(ii)/xk(ii)*del(ii)**(-0.5+ci*zmu)
  endif
3 continue; endif 
if(ibouss.eq.2.or.ibouss.eq.3) then ! Boussinesq variable shear
do ii=2,iim/2+1
if(ibouss.eq.2) zm=sld**2*xk(ii)**2-zj;if(ibouss.eq.3) zm=-zj
if(zm.ge.0)then; cm(ii)=sqrt(zm); else; cm(ii)=-sign(1.,xk(ii))*ci*sqrt(-zm);endif
aa(ii)=+5./4.+ci*zmu/2.-cm(ii)/2.
bb(ii)=-1./4.+ci*zmu/2.-cm(ii)/2.
cc=1.+ci*zmu
call test_neg_int(bb(ii)-cc+1.,logic);if(logic.and.ii.gt.1)goto 111
call test_neg_int(aa(ii)+bb(ii)-cc,logic);if(logic.and.ii.gt.1)goto 111
call test_neg_int(bb(ii),logic);if(logic.and.ii.gt.1)goto 111
call test_neg_int(aa(ii),logic);if(logic.and.ii.gt.1)goto 111
call test_neg_int(aa(ii)-cc+1.,logic);if(logic.and.ii.gt.1)goto 111
a1(ii)=cgamma(cc)*cgamma(cc-aa(ii)-bb(ii))/cgamma(cc-aa(ii))/cgamma(cc-bb(ii))
a2(ii)=cgamma(2.-cc)*cgamma(cc-aa(ii)-bb(ii))/cgamma(1.-aa(ii))/cgamma(1.-bb(ii))
a3(ii)=cgamma(cc)*cgamma(aa(ii)+bb(ii)-cc)/cgamma(aa(ii))/cgamma(bb(ii))
a4(ii)=cgamma(2.-cc)*cgamma(aa(ii)+bb(ii)-cc)/cgamma(aa(ii)-cc+1.)/cgamma(bb(ii)-cc+1.)
det(ii)=a1(ii)*a4(ii)-a2(ii)*a3(ii)
b1(ii)=+a1(ii)/det(ii)*2**(-cm(ii));b2(ii)=-a2(ii)/det(ii)*2**(-cm(ii))
goto 112 ; 111 continue
det(ii)=det(ii-1); b1(ii)=(999.9999,999.999); b2(ii)=(999.9999,999.999)
112 continue; enddo
do ii=3,iim/2+1
if(b1(ii).eq.(999.9999,999.999))then;b1(ii)=b1(ii-1); b2(ii)=b2(ii-1); endif
enddo
do ii=2,iim/2+1
oa1(ii)=b1(ii)*sld**(-0.5+ci*zmu) ; oa2(ii)=b2(ii)*sld**(-0.5-ci*zmu)
ta1(ii)=oa1(ii)/xk(ii)*del(ii)**(-0.5-ci*zmu)
ta2(ii)=oa2(ii)/xk(ii)*del(ii)**(-0.5+ci*zmu)
enddo
endif
if(zj.lt.0.25) then; ta1(:)=ta1(:)/(1.,1.)*sqrt(2.)
                     ta2(:)=ta2(:)/(1.,-1.)*sqrt(2.); endif

  if(lout)then
  filename='../xmvg/xxxx_tais_kall_xxxxxxxx.dat'
  write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
  330 format(a4);331 format(a8)
  open(32,file=filename);do ii=1,iim/2+1; write(32,301)xk(ii),ta1(ii),&
  sqrt(real(ta1(ii)*conjg(ta1(ii)))),ta2(ii),sqrt(real(ta2(ii)*conjg(ta2(ii))))
  enddo;close(32);endif;301  format(1x,f7.3,4(2(3X,E10.4,2X)))

! 3.2) a3, a4, b3, b4, c3, c4 associated with the viscous solutions 
!*****************************************************************

!chi:1,2 ; u:3,4 ; xsi:5,6 ; b:7,8 ; w:9,10 ; p:11,12

if(icond.gt.0) then

  ! zth used here has an inner altitude varying between 0 and 10

  do ll=llm,1,-1;lfirst=.true.;zth=10.*float(ll-1)/float(llm)
  visc1=viscous(zth,vis1,lfirst)
  do ii=1,6;cs1(ii,ll)=visc1(2*ii-1)+ci*visc1(2*ii);enddo;enddo
  do ll=llm,1,-1;lfirst=.true.;zth=10.*float(ll-1)/float(llm)
  visc2=viscous(zth,vis2,lfirst)
  do ii=1,6;cs2(ii,ll)=visc2(2*ii-1)+ci*visc2(2*ii);enddo;enddo
  do ll=llm,1,-1;lfirst=.true.;zth=10.*float(ll-1)/float(llm)
  visc3=viscous(zth,vis3,lfirst)
  do ii=1,6;cs3(ii,ll)=visc3(2*ii-1)+ci*visc3(2*ii);enddo;enddo
  do ll=llm,1,-1;lfirst=.true.;zth=10.*float(ll-1)/float(llm)
  visc4=viscous(zth,vis4,lfirst)
  do ii=1,6;cs4(ii,ll)=visc4(2*ii-1)+ci*visc4(2*ii);enddo;enddo

  cv1(:)=cs1(:,1);cv2(:)=cs2(:,1);cv3(:)=cs3(:,1);cv4(:)=cs4(:,1)

! The values in ii=1 should not be used except for check:
  if(ibouss.eq.0)then;ta1(1)=0.0;ta2(1)=(1.,1.)/sqrt(2.);endif
  if(ibouss.gt.1)then;ta1(1)=(1.,1.)/sqrt(2.);ta2(1)=conjg(ta1(1));endif
! print *,'cv12:',ta1(1)*cv1+ta2(1)*cv2
  do ii=1,iim/2+1
  cv12=ta1(ii)*cv1+ta2(ii)*cv2
  if(icond.eq.1)then
! Viscous sol with chi=0; xsi=0
    ta3(ii)= -(cv12(1)*cv4(3)-cv12(3)*cv4(1))/(cv3(1)*cv4(3)-cv3(3)*cv4(1))
    ta4(ii)= -(cv3(1)*cv12(3)-cv3(3)*cv12(1))/(cv3(1)*cv4(3)-cv3(3)*cv4(1))
! Viscous sol with chi=0, p=0
    tb3(ii)= -(cv12(1)*cv4(6)-cv12(6)*cv4(1))/(cv3(1)*cv4(6)-cv3(6)*cv4(1))
    tb4(ii)= -(cv3(1)*cv12(6)-cv3(6)*cv12(1))/(cv3(1)*cv4(6)-cv3(6)*cv4(1))
! Viscous sol with xsi=0, p=0
    tc3(ii)= -(cv12(3)*cv4(6)-cv12(6)*cv4(3))/(cv3(3)*cv4(6)-cv3(6)*cv4(3))
    tc4(ii)= -(cv3(3)*cv12(6)-cv3(6)*cv12(3))/(cv3(3)*cv4(6)-cv3(6)*cv4(3))
  endif
  if(icond.eq.2)then
! Viscous sol with chi=0; b=0
    ta3(ii)= -(cv12(1)*cv4(4)-cv12(4)*cv4(1))/(cv3(1)*cv4(4)-cv3(4)*cv4(1))
    ta4(ii)= -(cv3(1)*cv12(4)-cv3(4)*cv12(1))/(cv3(1)*cv4(4)-cv3(4)*cv4(1))
! Viscous sol with chi=0, p=0
    tb3(ii)= -(cv12(1)*cv4(6)-cv12(6)*cv4(1))/(cv3(1)*cv4(6)-cv3(6)*cv4(1))
    tb4(ii)= -(cv3(1)*cv12(6)-cv3(6)*cv12(1))/(cv3(1)*cv4(6)-cv3(6)*cv4(1))
! Viscous sol with b=0, p=0
    tc3(ii)= -(cv12(4)*cv4(6)-cv12(6)*cv4(4))/(cv3(4)*cv4(6)-cv3(6)*cv4(4))
    tc4(ii)= -(cv3(4)*cv12(6)-cv3(6)*cv12(4))/(cv3(4)*cv4(6)-cv3(6)*cv4(4))
  endif
  if(icond.eq.3)then
! Viscous sol with u=0; b=0
    ta3(ii)= -(cv12(2)*cv4(4)-cv12(4)*cv4(2))/(cv3(2)*cv4(4)-cv3(4)*cv4(2))
    ta4(ii)= -(cv3(2)*cv12(4)-cv3(4)*cv12(2))/(cv3(2)*cv4(4)-cv3(4)*cv4(2))
! Viscous sol with u=0, p=0
    tb3(ii)= -(cv12(2)*cv4(6)-cv12(6)*cv4(2))/(cv3(2)*cv4(6)-cv3(6)*cv4(2))
    tb4(ii)= -(cv3(2)*cv12(6)-cv3(6)*cv12(2))/(cv3(2)*cv4(6)-cv3(6)*cv4(2))
! Viscous sol with b=0, p=0
    tc3(ii)= -(cv12(4)*cv4(6)-cv12(6)*cv4(4))/(cv3(4)*cv4(6)-cv3(6)*cv4(4))
    tc4(ii)= -(cv3(4)*cv12(6)-cv3(6)*cv12(4))/(cv3(4)*cv4(6)-cv3(6)*cv4(4))
  endif
enddo

kk=iim/10
do ii=2,iim/2; if(xk(ii+1).gt.1.and.xk(ii).le.1.)kk=ii;enddo
csa(:,:)=ta1(kk)*cs1(:,:)+ta2(kk)*cs2(:,:)+ta3(kk)*cs3(:,:)+ta4(kk)*cs4(:,:)
csb(:,:)=ta1(kk)*cs1(:,:)+ta2(kk)*cs2(:,:)+tb3(kk)*cs3(:,:)+tb4(kk)*cs4(:,:)
csc(:,:)=ta1(kk)*cs1(:,:)+ta2(kk)*cs2(:,:)+tc3(kk)*cs3(:,:)+tc4(kk)*cs4(:,:)

  if(lout) then
   filename='../xmvg/xxxx_visc_sols_xxxxxxxx.dat'
   write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
   open(12,file=filename)
   do ii=1,6; do ll=1,llm;zth=10.*float(ll-1)/float(llm)
   write(12,*)sqrt(real(csa(ii,ll)*conjg(csa(ii,ll)))),zth    
   enddo;write(12,*)' ';enddo
   do ii=1,6; do ll=1,llm;zth=10.*float(ll-1)/float(llm)
   write(12,*)sqrt(real(csb(ii,ll)*conjg(csb(ii,ll)))),zth   
   enddo;write(12,*)' ';enddo
   do ii=1,6; do ll=1,llm;zth=10.*float(ll-1)/float(llm)
   write(12,*)sqrt(real(csc(ii,ll)*conjg(csc(ii,ll)))),zth    
   enddo;write(12,*)' ';enddo;close(12)
  endif
endif

!**********************************************************************
! 4. Surface values evaluated with RK protocole, harmonics by harmonics
!**********************************************************************

if(icond.gt.0) then

  if(icond.ne.1.and.icond.ne.2.and.icond.ne.3) then;print*,'fix icond';stop;endif
!  print *,'viscous 1:'
  cw1=0.;cu1=0.;cb1=0.;lfirst=.true.;do jj=2,iim/2+1;do ii=1,iim
  zth=zh(ii)/del(jj);visc=viscous(zth,vis1,lfirst)
  cw1(ii,jj)=visc(9)+ci*visc(10);cu1(ii,jj)=visc(3)+ci*visc(4)
  if(icond.eq.1)cb1(ii,jj)=visc(7)+ci*visc(8)
  if(icond.eq.2.or.icond.eq.3)cb1(ii,jj)=visc(5)+ci*visc(6)             
  if(icond.eq.3)cc1(ii,jj)=visc(1)+ci*visc(2);enddo;enddo

!  print *,'viscous 2:'
  cw2=0.;cu2=0.;cb2=0.;lfirst=.true.;do jj=2,iim/2+1;do ii=1,iim
  zth=zh(ii)/del(jj);visc=viscous(zth,vis2,lfirst)
  cw2(ii,jj)=visc(9)+ci*visc(10);cu2(ii,jj)=visc(3)+ci*visc(4)
  if(icond.eq.1)cb2(ii,jj)=visc(7)+ci*visc(8)
  if(icond.eq.2.or.icond.eq.3)cb2(ii,jj)=visc(5)+ci*visc(6)
  if(icond.eq.3)cc2(ii,jj)=visc(1)+ci*visc(2);enddo;enddo

!  print *,'viscous 3:'
  cw3=0.;cu3=0.;cb3=0.;lfirst=.true.;do jj=2,iim/2+1;do ii=1,iim
  zth=zh(ii)/del(jj);visc=viscous(zth,vis3,lfirst);cw3(ii,jj)=visc(9)+ci*visc(10)
  cu3(ii,jj)=visc(3)+ci*visc(4)
  if(icond.eq.1)cb3(ii,jj)=visc(7)+ci*visc(8)
  if(icond.eq.2.or.icond.eq.3)cb3(ii,jj)=visc(5)+ci*visc(6)
  if(icond.eq.3)cc3(ii,jj)=visc(1)+ci*visc(2);enddo;enddo

!  print *,'viscous 4:'
  cw4=0.;cu4=0.;cb4=0.;lfirst=.true.;do jj=2,iim/2+1;do ii=1,iim
  zth=zh(ii)/del(jj);visc=viscous(zth,vis4,lfirst);cw4(ii,jj)=visc(9)+ci*visc(10)
  cu4(ii,jj)=visc(3)+ci*visc(4)
  if(icond.eq.1) cb4(ii,jj)=visc(7)+ci*visc(8)
  if(icond.eq.2.or.icond.eq.3) cb4(ii,jj)=visc(5)+ci*visc(6)
  if(icond.eq.3)cc4(ii,jj)=visc(1)+ci*visc(2);enddo;enddo

! Build solutions with reasonable amplitudes near the surface.
do jj=2,iim/2+1
cua(:,jj)=ta1(jj)*cu1(:,jj)+ta2(jj)*cu2(:,jj)+ta3(jj)*cu3(:,jj)+ta4(jj)*cu4(:,jj)
cub(:,jj)=ta1(jj)*cu1(:,jj)+ta2(jj)*cu2(:,jj)+tb3(jj)*cu3(:,jj)+tb4(jj)*cu4(:,jj)
cuc(:,jj)=ta1(jj)*cu1(:,jj)+ta2(jj)*cu2(:,jj)+tc3(jj)*cu3(:,jj)+tc4(jj)*cu4(:,jj)
cwa(:,jj)=ta1(jj)*cw1(:,jj)+ta2(jj)*cw2(:,jj)+ta3(jj)*cw3(:,jj)+ta4(jj)*cw4(:,jj)
cwb(:,jj)=ta1(jj)*cw1(:,jj)+ta2(jj)*cw2(:,jj)+tb3(jj)*cw3(:,jj)+tb4(jj)*cw4(:,jj)
cwc(:,jj)=ta1(jj)*cw1(:,jj)+ta2(jj)*cw2(:,jj)+tc3(jj)*cw3(:,jj)+tc4(jj)*cw4(:,jj)
cba(:,jj)=ta1(jj)*cb1(:,jj)+ta2(jj)*cb2(:,jj)+ta3(jj)*cb3(:,jj)+ta4(jj)*cb4(:,jj)
cbb(:,jj)=ta1(jj)*cb1(:,jj)+ta2(jj)*cb2(:,jj)+tb3(jj)*cb3(:,jj)+tb4(jj)*cb4(:,jj)
cbc(:,jj)=ta1(jj)*cb1(:,jj)+ta2(jj)*cb2(:,jj)+tc3(jj)*cb3(:,jj)+tc4(jj)*cb4(:,jj)
enddo
if(icond.eq.3)then;do jj=2,iim/2+1
cca(:,jj)=ta1(jj)*cc1(:,jj)+ta2(jj)*cc2(:,jj)+ta3(jj)*cc3(:,jj)+ta4(jj)*cc4(:,jj)
ccb(:,jj)=ta1(jj)*cc1(:,jj)+ta2(jj)*cc2(:,jj)+tb3(jj)*cc3(:,jj)+tb4(jj)*cc4(:,jj)
ccc(:,jj)=ta1(jj)*cc1(:,jj)+ta2(jj)*cc2(:,jj)+tc3(jj)*cc3(:,jj)+tc4(jj)*cc4(:,jj)
enddo;endif

! Expressed in term of outer variables:
do jj=2,iim/2+1; do ii=1,iim; cwa(ii,jj)=xk(jj)*del(jj)*cwa(ii,jj)
cwb(ii,jj)=xk(jj)*del(jj)*cwb(ii,jj); cwc(ii,jj)=xk(jj)*del(jj)*cwc(ii,jj)
enddo; enddo
if (icond.eq.2.or.icond.eq.3)then;do jj=2,iim/2+1; do ii=1,iim
cba(ii,jj)=cba(ii,jj)/del(jj)
cbb(ii,jj)=cbb(ii,jj)/del(jj); cbc(ii,jj)=cbc(ii,jj)/del(jj)
enddo; enddo; endif
if (icond.eq.3)then;do jj=2,iim/2+1; do ii=1,iim
cca(ii,jj)=cca(ii,jj)/del(jj)
ccb(ii,jj)=ccb(ii,jj)/del(jj); ccc(ii,jj)=ccc(ii,jj)/del(jj)
enddo; enddo; endif

! Negative non-infinite wavelength
do jj=iim/2+2,iim
cwb(:,jj)=conjg(cwb(:,iim+2-jj))
cwc(:,jj)=conjg(cwc(:,iim+2-jj));cwa(:,jj)=conjg(cwa(:,iim+2-jj))
cub(:,jj)=conjg(cub(:,iim+2-jj));cuc(:,jj)=conjg(cuc(:,iim+2-jj))
cua(:,jj)=conjg(cua(:,iim+2-jj))
cbb(:,jj)=conjg(cbb(:,iim+2-jj));cbc(:,jj)=conjg(cbc(:,iim+2-jj))
cba(:,jj)=conjg(cba(:,iim+2-jj))
if (icond.eq.3) then;ccb(:,jj)=conjg(ccb(:,iim+2-jj))
ccc(:,jj)=conjg(ccc(:,iim+2-jj));cca(:,jj)=conjg(cca(:,iim+2-jj)); endif
enddo

!Reasonable values on jj=1:
cua(:,1)=cua(1,2);cub(:,1)=cub(1,2);cuc(:,1)=cuc(1,2)
cwa(:,1)=cwa(1,2);cwb(:,1)=cwb(1,2);cwc(:,1)=cwc(1,2)
cba(:,1)=cba(1,2);cbb(:,1)=cbb(1,2);cbc(:,1)=cbc(1,2)
!cua(:,1)=(1.,1.);cub(:,1)=(1.,1.);cuc(:,1)=(1.,1.)
!wa(:,1)=(1.,1.);cwb(:,1)=(1.,1.);cwc(:,1)=(1.,1.)
!ba(:,1)=(1.,1.);cbb(:,1)=(1.,1.);cbc(:,1)=(1.,1.)
if (icond.eq.3) then
cca(:,1)=cca(1,2);ccb(:,1)=ccb(1,2);ccc(:,1)=ccc(1,2) 
endif

if(lout)print *,' viscous 1 to 4  done'

else  ! End of viscous case

if(lout)print *,' Here you are in the inviscid case....'

do jj=2,iim/2+1;do ii=1,iim
! Booker and Bretherton type:
! Matching function
visc1=bab1(xk(jj),c1*zh(ii))
visc2=bab2(xk(jj),c1*zh(ii))  
cw1(ii,jj)=visc1(9)+ci*visc1(10);cu1(ii,jj)=visc1(3)+ci*visc1(4)
cw2(ii,jj)=visc2(9)+ci*visc2(10);cu2(ii,jj)=visc2(3)+ci*visc2(4)
cwc(ii,jj)=oa1(jj)*cw1(ii,jj)+oa2(jj)*cw2(ii,jj)
cuc(ii,jj)=oa1(jj)*cu1(ii,jj)+oa2(jj)*cu2(ii,jj)

! Damped solution
visc1=bab1(xk(jj),c1*zh(ii)-ci*del(jj))
visc2=bab2(xk(jj),c1*zh(ii)-ci*del(jj))
cw1(ii,jj)=visc1(9)+ci*visc1(10);cu1(ii,jj)=visc1(3)+ci*visc1(4)
cw2(ii,jj)=visc2(9)+ci*visc2(10);cu2(ii,jj)=visc2(3)+ci*visc2(4)
cwd(ii,jj)=oa1(jj)*cw1(ii,jj)+oa2(jj)*cw2(ii,jj)
cud(ii,jj)=oa1(jj)*cu1(ii,jj)+oa2(jj)*cu2(ii,jj)

! Exact inviscid
if(zh(ii).gt.epd)then
if(ibouss.eq.0) then;visc2=bab2(xk(jj),c1*zh(ii))
cw2(ii,jj)=visc2(9)+ci*visc2(10);cu2(ii,jj)=visc2(3)+ci*visc2(4)
cwb(ii,jj)=oa2(jj)*cw2(ii,jj);cub(ii,jj)=oa2(jj)*cu2(ii,jj); endif
if(ibouss.eq.1) then;visc1=hanksol(xk(jj),zh(ii))
cwb(ii,jj)=visc1(9)+ci*visc1(10);cub(ii,jj)=visc1(3)+ci*visc1(4); endif
if(ibouss.ge.2) then; visc1=hypsol(aa(jj),bb(jj),cc,b1(jj),b2(jj),xk(jj),zh(ii))
cwb(ii,jj)=visc1(9)+ci*visc1(10);cub(ii,jj)=visc1(3)+ci*visc1(4); endif
else
cwb(ii,jj)=cwc(ii,jj);cub(ii,jj)=cuc(ii,jj)
endif
!visc1=bab1(xk(jj),c1*zh(ii)-ci*del(jj))
!Uniform approximation:
visc=visc1+visc4-visc3
cwa(ii,jj)=cwd(ii,jj)+cwb(ii,jj)-cwc(ii,jj)
cua(ii,jj)=cud(ii,jj)+cub(ii,jj)-cuc(ii,jj)
enddo;enddo
! Negative non-infinite wavelength
do jj=iim/2+2,iim
cwa(:,jj)=conjg(cwa(:,iim+2-jj));cua(:,jj)=conjg(cua(:,iim+2-jj))
enddo
!Reasonable values on jj=1:
cua(:,1)=cua(1,2);cwa(:,1)=cwa(1,2);cba(:,1)=cba(1,2)

endif ! End of inviscid case

!****************************************
! 5. Inversion of the boundary condition
!****************************************

cfa=(0.,0.);cfb=(0.,0.);cfc=(0.,0.)

if(icond.gt.0)then

  if(icond.eq.3) then
  cua(:,1)=-cwa(1,2); cub(:,1)=-cwb(1,2);cuc(:,1)=-cwc(1,2)
  do jj=2,iim;
  cua(:,jj)=cua(:,jj)*zhx(:)-cwa(:,jj);cub(:,jj)=cub(:,jj)*zhx(:)-cwb(:,jj)
  cuc(:,jj)=cuc(:,jj)*zhx(:)-cwc(:,jj);enddo
  cwa=cca;cwb=ccb;cwc=ccc;endif

  do ii=1,iim;do jj=1,iim
  cua(ii,jj)=cua(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cub(ii,jj)=cub(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cuc(ii,jj)=cuc(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cwa(ii,jj)=cwa(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cwb(ii,jj)=cwb(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cwc(ii,jj)=cwc(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cba(ii,jj)=cba(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cbb(ii,jj)=cbb(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  cbc(ii,jj)=cbc(ii,jj)*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  enddo;enddo

! print *,'build the complete matrix....cua:',cua
  cu2=cua; call cgeco(cu2,iim,iim, ipvt,rcond,cvec); call cgedi(cu2,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cu3=cub; call cgeco(cu3,iim,iim, ipvt,rcond,cvec); call cgedi(cu3,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cu4=cuc; call cgeco(cu4,iim,iim, ipvt,rcond,cvec); call cgedi(cu4,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  
  cw2=cwa; call cgeco(cw2,iim,iim, ipvt,rcond,cvec); call cgedi(cw2,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cw3=cwb; call cgeco(cw3,iim,iim, ipvt,rcond,cvec); call cgedi(cw3,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cw4=cwc; call cgeco(cw4,iim,iim, ipvt,rcond,cvec); call cgedi(cw4,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cb2=cba; call cgeco(cb2,iim,iim, ipvt,rcond,cvec); call cgedi(cb2,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cb3=cbb; call cgeco(cb3,iim,iim, ipvt,rcond,cvec); call cgedi(cb3,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cb4=cbc; call cgeco(cb4,iim,iim, ipvt,rcond,cvec); call cgedi(cb4,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))

  if(lout) print *,'Calculation of FA:'
  cma=matmul(cw3,cwa)-matmul(cu3,cua);cmb=matmul(cw3,cwc)-matmul(cu3,cuc)
  cmc=matmul(cw3,cwa)-matmul(cb3,cba);cmd=matmul(cw3,cwc)-matmul(cb3,cbc)
  cm2=cmb;call cgeco(cm2,iim,iim,ipvt,rcond,cvec); call cgedi(cm2,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cm4=cmd;call cgeco(cm4,iim,iim,ipvt,rcond,cvec); call cgedi(cm4,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cme=matmul(cm2,cma)-matmul(cm4,cmc);
  cm5=cme;call cgeco(cm5,iim,iim,ipvt,rcond,cvec); call cgedi(cm5,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
if(icond.eq.1)  cmf=matmul(cm2,cu3)-matmul(cm4,cb3)*zj
if(icond.ge.2)  cmf=matmul(cm2,cu3)
if(icond.lt.3) then; cvec=zh; else; cvec=zh*zhx; endif
cmh=matmul(cm5,cmf);cfa=matmul(cmh,cvec)

  if(lout) print *,'Calculation of FB:'
  cma=matmul(cw2,cwb)-matmul(cu2,cub);cmb=matmul(cw2,cwc)-matmul(cu2,cuc)
  cmc=matmul(cw2,cwb)-matmul(cb2,cbb);cmd=matmul(cw2,cwc)-matmul(cb2,cbc)
  cm2=cmb;call cgeco(cm2,iim,iim,ipvt,rcond,cvec); call cgedi(cm2,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cm4=cmd;call cgeco(cm4,iim,iim,ipvt,rcond,cvec); call cgedi(cm4,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cme=matmul(cm2,cma)-matmul(cm4,cmc);
  cm5=cme;call cgeco(cm5,iim,iim,ipvt,rcond,cvec); call cgedi(cm5,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
if(icond.eq.1)  cmf=matmul(cm2,cu2)-matmul(cm4,cb2)*zj
if(icond.ge.2)  cmf=matmul(cm2,cu2)                      
if(icond.lt.3) then; cvec=zh; else; cvec=zh*zhx; endif
cmh=matmul(cm5,cmf);cfb=matmul(cmh,cvec)

  if(lout) print *,'Calculation of FC:'
  cma=matmul(cw2,cwb)-matmul(cu2,cub);cmb=matmul(cw2,cwc)-matmul(cu2,cuc)
  cmc=matmul(cw2,cwb)-matmul(cb2,cbb);cmd=matmul(cw2,cwc)-matmul(cb2,cbc)
  cm1=cma;call cgeco(cm1,iim,iim,ipvt,rcond,cvec); call cgedi(cm1,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cm3=cmc;call cgeco(cm3,iim,iim,ipvt,rcond,cvec); call cgedi(cm3,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cme=matmul(cm1,cmb)-matmul(cm3,cmd);
  cm5=cme;call cgeco(cm5,iim,iim,ipvt,rcond,cvec); call cgedi(cm5,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  if(icond.eq.1) cmf=matmul(cm1,cu2)-matmul(cm3,cb2)*zj
  if(icond.ge.2) cmf=matmul(cm1,cu2)                        
if(icond.lt.3) then; cvec=zh; else; cvec=zh*zhx; endif
cmh=matmul(cm5,cmf);cfc=matmul(cmh,cvec)

  if(lout) print *,'All viscous about inversion?'

 else
  if(lout) print *,'build the matrixes'
  do ii=1,iim;do jj=1,iim
  cm1(ii,jj)=(cwa(ii,jj)-cua(ii,jj)*zhx(ii))*exp(ci*2.*pi*float((ii-1)*(jj-1))/float(iim))
  enddo; enddo
  if(lout) print *,'end of the matrixes'
  
!print *,'cua:',cua
  call cgeco(cm1,iim,iim,ipvt,rcond,cvec); call cgedi(cm1,iim,iim, ipvt,fdet,cvec,11)
  if(lout) print *,'fdet:',real(sqrt(fdet(1)*conjg(fdet(1)))),real(fdet(2))
  cvec=zh*zhx;cfa=matmul(cm1,cvec)
  if(lout) print *,'All Inviscid inversion?'
 endif
 
 if(lout)then;filename='../xmvg/xxxx_fofk_func_xxxxxxxx.dat'
 write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
 open(42,file=filename);401 format(1x,f10.5,2(3X,e10.3,1X,E10.3)); do ii=1,iim 
 write(42,401)xk(ii),sqrt(real(cfa(ii)*conjg(cfa(ii)))),&
 sqrt(real(cfb(ii)*conjg(cfb(ii)))),sqrt(real(cfc(ii)*conjg(cfc(ii))))
 enddo;close(42); endif

!True verification:

if(lout)print *,'Enter true verif.....'
if(lout.and.icond.gt.0)then; filename='../xmvg/xxxx_uwbpr_sol_xxxxxxxx.dat'
write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
open(41,file=filename)

do ii=1,iim;tiz(:)=zh(ii)/del(:)

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis1,lfirst)
  cwvis1(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis1(jj)=visc(3)+ci*visc(4);cbvis1(jj)=(visc(7)+ci*visc(8))
  ccvis1(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis1(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis1(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis2,lfirst)
  cwvis2(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis2(jj)=visc(3)+ci*visc(4);cbvis2(jj)=(visc(7)+ci*visc(8))
  ccvis2(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis2(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis2(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis3,lfirst)
  cwvis3(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis3(jj)=visc(3)+ci*visc(4);cbvis3(jj)=(visc(7)+ci*visc(8))
  ccvis3(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis3(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis3(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis4,lfirst)
  cwvis4(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis4(jj)=visc(3)+ci*visc(4);cbvis4(jj)=(visc(7)+ci*visc(8))
  ccvis4(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis4(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis4(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  do jj=2,iim/2+1
  cwvisd(jj)=ta1(jj)*cwvis1(jj)+ta2(jj)*cwvis2(jj)
  cwvisa(jj)=cwvisd(jj)+ta3(jj)*cwvis3(jj)+ta4(jj)*cwvis4(jj)
  cwvisb(jj)=cwvisd(jj)+tb3(jj)*cwvis3(jj)+tb4(jj)*cwvis4(jj)
  cwvisc(jj)=cwvisd(jj)+tc3(jj)*cwvis3(jj)+tc4(jj)*cwvis4(jj)
  cuvisd(jj)=ta1(jj)*cuvis1(jj)+ta2(jj)*cuvis2(jj)
  cuvisa(jj)=cuvisd(jj)+ta3(jj)*cuvis3(jj)+ta4(jj)*cuvis4(jj)
  cuvisb(jj)=cuvisd(jj)+tb3(jj)*cuvis3(jj)+tb4(jj)*cuvis4(jj)
  cuvisc(jj)=cuvisd(jj)+tc3(jj)*cuvis3(jj)+tc4(jj)*cuvis4(jj)
  cbvisd(jj)=ta1(jj)*cbvis1(jj)+ta2(jj)*cbvis2(jj)
  cbvisa(jj)=cbvisd(jj)+ta3(jj)*cbvis3(jj)+ta4(jj)*cbvis4(jj)
  cbvisb(jj)=cbvisd(jj)+tb3(jj)*cbvis3(jj)+tb4(jj)*cbvis4(jj)
  cbvisc(jj)=cbvisd(jj)+tc3(jj)*cbvis3(jj)+tc4(jj)*cbvis4(jj)
  cpvisd(jj)=ta1(jj)*cpvis1(jj)+ta2(jj)*cpvis2(jj)
  cpvisa(jj)=cpvisd(jj)+ta3(jj)*cpvis3(jj)+ta4(jj)*cpvis4(jj)
  cpvisb(jj)=cpvisd(jj)+tb3(jj)*cpvis3(jj)+tb4(jj)*cpvis4(jj)
  cpvisc(jj)=cpvisd(jj)+tc3(jj)*cpvis3(jj)+tc4(jj)*cpvis4(jj)
  ccvisd(jj)=ta1(jj)*ccvis1(jj)+ta2(jj)*ccvis2(jj)
  ccvisa(jj)=ccvisd(jj)+ta3(jj)*ccvis3(jj)+ta4(jj)*ccvis4(jj)
  ccvisb(jj)=ccvisd(jj)+tb3(jj)*ccvis3(jj)+tb4(jj)*ccvis4(jj)
  ccvisc(jj)=ccvisd(jj)+tc3(jj)*ccvis3(jj)+tc4(jj)*ccvis4(jj)
  cxvisd(jj)=ta1(jj)*cxvis1(jj)+ta2(jj)*cxvis2(jj)
  cxvisa(jj)=cxvisd(jj)+ta3(jj)*cxvis3(jj)+ta4(jj)*cxvis4(jj)
  cxvisb(jj)=cxvisd(jj)+tb3(jj)*cxvis3(jj)+tb4(jj)*cxvis4(jj)
  cxvisc(jj)=cxvisd(jj)+tc3(jj)*cxvis3(jj)+tc4(jj)*cxvis4(jj)
  enddo

!Restore values on jj=1

if(icond.eq.1)then
  cwvisa(1)=cwa(1,1);cwvisb(1)=cwb(1,1);cwvisc(1)=cwc(1,1)
  cuvisa(1)=cua(1,1);cuvisb(1)=cub(1,1);cuvisc(1)=cuc(1,1)
  cbvisa(1)=cba(1,1);cbvisb(1)=cbb(1,1);cbvisc(1)=cbc(1,1)
endif
if(icond.eq.2)then
  cwvisa(1)=cwa(1,1);cwvisb(1)=cwb(1,1);cwvisc(1)=cwc(1,1)
  cuvisa(1)=cua(1,1);cuvisb(1)=cub(1,1);cuvisc(1)=cuc(1,1)
  cxvisa(1)=cba(1,1);cxvisb(1)=cbb(1,1);cxvisc(1)=cbc(1,1)
endif
if(icond.eq.3)then
  cuvisa(1)=cua(1,1);cuvisb(1)=cub(1,1);cuvisc(1)=cuc(1,1)
  cwvisa(1)=(0.,0.) ;cwvisb(1)=(0.,0.) ;cwvisc(1)=(0.,0.) 
  ccvisa(1)=cca(1,1);ccvisb(1)=ccb(1,1);ccvisc(1)=ccc(1,1)
  cxvisa(1)=cba(1,1);cxvisb(1)=cbb(1,1);cxvisc(1)=cbc(1,1)
endif

  do jj=1,iim/2+1
  cwvis(jj)=cfb(jj)*cwvisb(jj)+cfc(jj)*cwvisc(jj)+cfa(jj)*cwvisa(jj)
  cuvis(jj)=cfb(jj)*cuvisb(jj)+cfc(jj)*cuvisc(jj)+cfa(jj)*cuvisa(jj)
  cbvis(jj)=cfb(jj)*cbvisb(jj)+cfc(jj)*cbvisc(jj)+cfa(jj)*cbvisa(jj)
  cpvis(jj)=cfb(jj)*cpvisb(jj)+cfc(jj)*cpvisc(jj)+cfa(jj)*cpvisa(jj)
  ccvis(jj)=cfb(jj)*ccvisb(jj)+cfc(jj)*ccvisc(jj)+cfa(jj)*ccvisa(jj)
  cxvis(jj)=cfb(jj)*cxvisb(jj)+cfc(jj)*cxvisc(jj)+cfa(jj)*cxvisa(jj)
  enddo
  do jj=iim/2+2,iim
  cwvis(jj)=conjg(cwvis(iim-jj+2));cbvis(jj)=conjg(cbvis(iim-jj+2))
  cuvis(jj)=conjg(cuvis(iim-jj+2));cpvis(jj)=conjg(cpvis(iim-jj+2))
  ccvis(jj)=conjg(ccvis(iim-jj+2));cxvis(jj)=conjg(cxvis(iim-jj+2))
  enddo

call cffti1(iim,w2,ifac)   !  Needed to initialize the fft (first call)
call cfftb1(iim,cwvis,w1,w2,ifac);wsth(ii)=real(cwvis(ii))
call cfftb1(iim,cuvis,w1,w2,ifac);usth(ii)=real(cuvis(ii))
call cfftb1(iim,cbvis,w1,w2,ifac);bsth(ii)=real(cbvis(ii))
call cfftb1(iim,cpvis,w1,w2,ifac);psth(ii)=real(cpvis(ii))
call cfftb1(iim,ccvis,w1,w2,ifac);uzsth(ii)=real(ccvis(ii))
call cfftb1(iim,cxvis,w1,w2,ifac);bzsth(ii)=real(cxvis(ii))

enddo

do ii=1,iim;write(41,*)xx(ii),usth(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),wsth(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),bsth(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),psth(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),uzsth(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),bzsth(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),-zh(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),-zj *zh(ii);enddo;write(41,*)
do ii=1,iim;write(41,*)xx(ii),(zh(ii)+usth(ii))*zhx(ii)-wsth(ii);enddo

print *,'Exit true verif.....'
close(41);endif

!cfb(iim/2+1)=0.;cfa(iim/2+1)=0.;cfc(iim/2+1)=0.

!***********************
!  6:  Reconstruction  *
!***********************

if(lout) print *,'Debut de la reconstruction:'

call cffti1(iim,w2,ifac)   !  Needed to initialize the fft (first call)

do ll=1,llm     !   Boucle sur l'altitude

if(icond.gt.0)then ! Reconstruct the viscous solutions:'

  tiz(:)=zz(ll)/del(:)     !   Boucle sur l'altitude

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis1,lfirst)
  cwvis1(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis1(jj)=visc(3)+ci*visc(4);cbvis1(jj)=(visc(7)+ci*visc(8))
  ccvis1(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis1(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis1(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis2,lfirst)
  cwvis2(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis2(jj)=visc(3)+ci*visc(4);cbvis2(jj)=(visc(7)+ci*visc(8))  
  ccvis2(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis2(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis2(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis3,lfirst)
  cwvis3(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis3(jj)=visc(3)+ci*visc(4);cbvis3(jj)=(visc(7)+ci*visc(8))  
  ccvis3(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis3(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis3(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  lfirst=.true.;do jj=2,iim/2+1; zth=tiz(jj);visc=viscous(zth,vis4,lfirst)
  cwvis4(jj)=abs(xk(jj))*del(jj)*(visc(9)+ci*visc(10))
  cuvis4(jj)=visc(3)+ci*visc(4);cbvis4(jj)=(visc(7)+ci*visc(8))  
  ccvis4(jj)=(visc(1)+ci*visc(2))/del(jj)
  cxvis4(jj)=(visc(5)+ci*visc(6))/del(jj)
  cpvis4(jj)=del(jj)*(visc(11)+ci*visc(12));enddo

  do jj=2,iim/2+1
  cwvisd(jj)=ta1(jj)*cwvis1(jj)+ta2(jj)*cwvis2(jj)
  cwvisa(jj)=cwvisd(jj)+ta3(jj)*cwvis3(jj)+ta4(jj)*cwvis4(jj)
  cwvisb(jj)=cwvisd(jj)+tb3(jj)*cwvis3(jj)+tb4(jj)*cwvis4(jj)
  cwvisc(jj)=cwvisd(jj)+tc3(jj)*cwvis3(jj)+tc4(jj)*cwvis4(jj)
  cuvisd(jj)=ta1(jj)*cuvis1(jj)+ta2(jj)*cuvis2(jj)
  cuvisa(jj)=cuvisd(jj)+ta3(jj)*cuvis3(jj)+ta4(jj)*cuvis4(jj)
  cuvisb(jj)=cuvisd(jj)+tb3(jj)*cuvis3(jj)+tb4(jj)*cuvis4(jj)
  cuvisc(jj)=cuvisd(jj)+tc3(jj)*cuvis3(jj)+tc4(jj)*cuvis4(jj)
  cbvisd(jj)=ta1(jj)*cbvis1(jj)+ta2(jj)*cbvis2(jj)
  cbvisa(jj)=cbvisd(jj)+ta3(jj)*cbvis3(jj)+ta4(jj)*cbvis4(jj)
  cbvisb(jj)=cbvisd(jj)+tb3(jj)*cbvis3(jj)+tb4(jj)*cbvis4(jj)
  cbvisc(jj)=cbvisd(jj)+tc3(jj)*cbvis3(jj)+tc4(jj)*cbvis4(jj)
  cpvisd(jj)=ta1(jj)*cpvis1(jj)+ta2(jj)*cpvis2(jj)
  cpvisa(jj)=cpvisd(jj)+ta3(jj)*cpvis3(jj)+ta4(jj)*cpvis4(jj)
  cpvisb(jj)=cpvisd(jj)+tb3(jj)*cpvis3(jj)+tb4(jj)*cpvis4(jj)
  cpvisc(jj)=cpvisd(jj)+tc3(jj)*cpvis3(jj)+tc4(jj)*cpvis4(jj)
  ccvisd(jj)=ta1(jj)*ccvis1(jj)+ta2(jj)*ccvis2(jj)
  ccvisa(jj)=ccvisd(jj)+ta3(jj)*ccvis3(jj)+ta4(jj)*ccvis4(jj)
  ccvisb(jj)=ccvisd(jj)+tb3(jj)*ccvis3(jj)+tb4(jj)*ccvis4(jj)
  ccvisc(jj)=ccvisd(jj)+tc3(jj)*ccvis3(jj)+tc4(jj)*ccvis4(jj)
  cxvisd(jj)=ta1(jj)*cxvis1(jj)+ta2(jj)*cxvis2(jj)
  cxvisa(jj)=cxvisd(jj)+ta3(jj)*cxvis3(jj)+ta4(jj)*cxvis4(jj)
  cxvisb(jj)=cxvisd(jj)+tb3(jj)*cxvis3(jj)+tb4(jj)*cxvis4(jj)
  cxvisc(jj)=cxvisd(jj)+tc3(jj)*cxvis3(jj)+tc4(jj)*cxvis4(jj)
  enddo

  if(icond.eq.1)then
    cwvisa(1)=cwa(1,1);cwvisb(1)=cwb(1,1);cwvisc(1)=cwc(1,1)
    cuvisa(1)=cua(1,1);cuvisb(1)=cub(1,1);cuvisc(1)=cuc(1,1)
    cbvisa(1)=cba(1,1);cbvisb(1)=cbb(1,1);cbvisc(1)=cbc(1,1)
  endif
  if(icond.eq.2)then
    cwvisa(1)=cwa(1,1);cwvisb(1)=cwb(1,1);cwvisc(1)=cwc(1,1)
    cuvisa(1)=cua(1,1);cuvisb(1)=cub(1,1);cuvisc(1)=cuc(1,1)
    cxvisa(1)=cba(1,1);cxvisb(1)=cbb(1,1);cxvisc(1)=cbc(1,1)
  endif
  if(icond.eq.3)then
    cuvisa(1)=cua(1,1);cuvisb(1)=cub(1,1);cuvisc(1)=cuc(1,1)
    cwvisa(1)=(0.,0.) ;cwvisb(1)=(0.,0.) ;cwvisc(1)=(0.,0.) 
    ccvisa(1)=cca(1,1);ccvisb(1)=ccb(1,1);ccvisc(1)=ccc(1,1)
    cxvisa(1)=cba(1,1);cxvisb(1)=cbb(1,1);cxvisc(1)=cbc(1,1)
  endif

  do jj=1,iim/2+1
  cwvis(jj)=cfb(jj)*cwvisb(jj)+cfc(jj)*cwvisc(jj)+cfa(jj)*cwvisa(jj)
  cuvis(jj)=cfb(jj)*cuvisb(jj)+cfc(jj)*cuvisc(jj)+cfa(jj)*cuvisa(jj)
  cbvis(jj)=cfb(jj)*cbvisb(jj)+cfc(jj)*cbvisc(jj)+cfa(jj)*cbvisa(jj)
  cpvis(jj)=cfb(jj)*cpvisb(jj)+cfc(jj)*cpvisc(jj)+cfa(jj)*cpvisa(jj)
  ccvis(jj)=cfb(jj)*ccvisb(jj)+cfc(jj)*ccvisc(jj)+cfa(jj)*ccvisa(jj)
  cxvis(jj)=cfb(jj)*cxvisb(jj)+cfc(jj)*cxvisc(jj)+cfa(jj)*cxvisa(jj)
  enddo
  do ii=iim/2+2,iim
  cwvis(ii)=conjg(cwvis(iim-ii+2));cbvis(ii)=conjg(cbvis(iim-ii+2))
  cuvis(ii)=conjg(cuvis(iim-ii+2));cpvis(ii)=conjg(cpvis(iim-ii+2))
  ccvis(ii)=conjg(ccvis(iim-ii+2));cxvis(ii)=conjg(cxvis(iim-ii+2))
  enddo

  if(zz(ll).gt.epd)then
! Inviscid solution:
    do jj=2,iim/2+1
    if(ibouss.eq.0)then; visc=bab2(xk(jj),zz(ll)*c1);endif
    if(ibouss.eq.1)then; visc=hanksol(xk(jj),zz(ll));endif
    if(ibouss.ge.2)then
    visc=hypsol(aa(jj),bb(jj),cc,b1(jj),b2(jj),xk(jj),zz(ll));endif
    cwinv(jj)=(cfa(jj)+cfb(jj)+cfc(jj))*(visc(9)+ci*visc(10))
    cuinv(jj)=(cfa(jj)+cfb(jj)+cfc(jj))*(visc(3)+ci*visc(4))
    cbinv(jj)=(cfa(jj)+cfb(jj)+cfc(jj))*(visc(7)+ci*visc(8))
    cpinv(jj)=(cfa(jj)+cfb(jj)+cfc(jj))*(visc(11)+ci*visc(12))
    ccinv(jj)=(cfa(jj)+cfb(jj)+cfc(jj))*(visc(1)+ci*visc(2))
    cxinv(jj)=(cfa(jj)+cfb(jj)+cfc(jj))*(visc(5)+ci*visc(6))
    enddo
    do ii=iim/2+2,iim
    cwinv(ii)=conjg(cwinv(iim-ii+2));cbinv(ii)=conjg(cbinv(iim-ii+2))
    cuinv(ii)=conjg(cuinv(iim-ii+2));cpinv(ii)=conjg(cpinv(iim-ii+2))
    ccinv(ii)=conjg(ccinv(iim-ii+2));cxinv(ii)=conjg(cxinv(iim-ii+2))
    enddo
! Matching solution:
    do jj=2,iim/2+1
    visc=bab1(xk(jj),c1*zz(ll))
    cwmat(jj)=oa1(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(9)+ci*visc(10))
    cumat(jj)=oa1(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(3)+ci*visc(4))
    cbmat(jj)=oa1(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(7)+ci*visc(8))
    cpmat(jj)=oa1(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(11)+ci*visc(12))
    ccmat(jj)=oa1(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(1)+ci*visc(2))
    cxmat(jj)=oa1(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(5)+ci*visc(6))
    visc=bab2(xk(jj),c1*zz(ll))
    cwmat(jj)=oa2(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(9)+ci*visc(10))+cwmat(jj)
    cumat(jj)=oa2(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(3)+ci*visc(4))+cumat(jj)
    cbmat(jj)=oa2(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(7)+ci*visc(8))+cbmat(jj)
    cpmat(jj)=oa2(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(11)+ci*visc(12))+cpmat(jj)
    ccmat(jj)=oa2(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(1)+ci*visc(2))+ccmat(jj)
    cxmat(jj)=oa2(jj)*(cfa(jj)+cfb(jj)+cfc(jj))*(visc(5)+ci*visc(6))+cxmat(jj)
    enddo
    do ii=iim/2+2,iim
    cwmat(ii)=conjg(cwmat(iim-ii+2));cbmat(ii)=conjg(cbmat(iim-ii+2))
    cumat(ii)=conjg(cumat(iim-ii+2));cpmat(ii)=conjg(cpmat(iim-ii+2))
    ccmat(ii)=conjg(ccmat(iim-ii+2));cxmat(ii)=conjg(cxmat(iim-ii+2))
    enddo
  else
    cwinv=0.;cuinv=0.;cbinv=0.;cpinv=0.;ccinv=0.;cxinv=0.
    cwmat=0.;cumat=0.;cbmat=0.;cpmat=0.;ccmat=0.;cxmat=0.
  endif
  cuinv(1)=cuvis(1);cwinv(1)=cwvis(1);cbinv(1)=cbvis(1)
  cpinv(1)=cpvis(1);ccinv(1)=ccvis(1);cxinv(1)=cxvis(1)
  cumat(1)=cuvis(1);cwmat(1)=cwvis(1);cbmat(1)=cbvis(1)
  cpmat(1)=cpvis(1);ccmat(1)=ccvis(1);cxmat(1)=cxvis(1)
  do ii=2,iim; if(abs(xk(ii)*zz(ll)).le.0.025)then
     cwinv(ii)=cwmat(ii);cuinv(ii)=cumat(ii);cbinv(ii)=cbmat(ii);
     cpinv(ii)=cpmat(ii);ccinv(ii)=ccmat(ii);cxinv(ii)=cxmat(ii);endif
  enddo

else  !  End of viscous; beginning of inviscid (icond=0)

  do jj=2,iim/2+1
! Booker and Bretherton type:
! Matching function
visc1=bab1(xk(jj),c1*zz(ll));visc2=bab2(xk(jj),c1*zz(ll))  
cwvis1(jj)=visc1(9)+ci*visc1(10);cuvis1(jj)=visc1(3)+ci*visc1(4)
cbvis1(jj)=visc1(7)+ci*visc1(8);cpvis1(jj)=visc1(11)+ci*visc1(12)
cwvis2(jj)=visc2(9)+ci*visc2(10);cuvis2(jj)=visc2(3)+ci*visc2(4)
cbvis2(jj)=visc2(7)+ci*visc2(8);cpvis2(jj)=visc2(11)+ci*visc2(12)
ccvis2(jj)=visc2(1)+ci*visc2(2);ccvis1(jj)=visc1(1)+ci*visc1(2)
cwmat(jj)=cfa(jj)*(oa1(jj)*cwvis1(jj)+oa2(jj)*cwvis2(jj))
cumat(jj)=cfa(jj)*(oa1(jj)*cuvis1(jj)+oa2(jj)*cuvis2(jj))
cbmat(jj)=cfa(jj)*(oa1(jj)*cbvis1(jj)+oa2(jj)*cbvis2(jj))
cpmat(jj)=cfa(jj)*(oa1(jj)*cpvis1(jj)+oa2(jj)*cpvis2(jj))
ccmat(jj)=cfa(jj)*(oa1(jj)*ccvis1(jj)+oa2(jj)*ccvis2(jj))

! Damped solution
visc1=bab1(xk(jj),c1*zz(ll)-ci*del(jj));visc2=bab2(xk(jj),c1*zz(ll)-ci*del(jj))  
cwvis1(jj)=visc1(9)+ci*visc1(10);cuvis1(jj)=visc1(3)+ci*visc1(4)
cbvis1(jj)=visc1(7)+ci*visc1(8);cpvis1(jj)=visc1(11)+ci*visc1(12)
ccvis1(jj)=visc1(1)+ci*visc1(2)
cwvis2(jj)=visc2(9)+ci*visc2(10);cuvis2(jj)=visc2(3)+ci*visc2(4)
cbvis2(jj)=visc2(7)+ci*visc2(8);cpvis2(jj)=visc2(11)+ci*visc2(12)
ccvis2(jj)=visc2(1)+ci*visc2(2)
cwvis(jj)=cfa(jj)*(oa1(jj)*cwvis1(jj)+oa2(jj)*cwvis2(jj))
cuvis(jj)=cfa(jj)*(oa1(jj)*cuvis1(jj)+oa2(jj)*cuvis2(jj))
cbvis(jj)=cfa(jj)*(oa1(jj)*cbvis1(jj)+oa2(jj)*cbvis2(jj))
cpvis(jj)=cfa(jj)*(oa1(jj)*cpvis1(jj)+oa2(jj)*cpvis2(jj))
ccvis(jj)=cfa(jj)*(oa1(jj)*ccvis1(jj)+oa2(jj)*ccvis2(jj))


! Exact inviscid (Hankel Functions)
! if(abs(zz(ll)*xk(jj)).gt.*epd.and.ibouss.eq.1)then
  if(ibouss.eq.0)then; cwinv(jj)=cwmat(jj);cuinv(jj)=cumat(jj)
    cbinv(jj)=cbmat(jj);cpinv(jj)=cpmat(jj);ccinv(jj)=ccmat(jj);endif
  if(ibouss.eq.1)then; visc1=hanksol(xk(jj),zz(ll))
    cwinv(jj)=cfa(jj)*(visc1(9)+ci*visc1(10))
    cuinv(jj)=cfa(jj)*(visc1(3)+ci*visc1(4))
    cbinv(jj)=cfa(jj)*(visc1(7)+ci*visc1(8))
    cpinv(jj)=cfa(jj)*(visc1(11)+ci*visc1(12))
    ccinv(jj)=cfa(jj)*(visc1(1)+ci*visc1(2)); endif
  if(ibouss.ge.2)then
    visc1=hypsol(aa(jj),bb(jj),cc,b1(jj),b2(jj),xk(jj),zz(ll))
    cwinv(jj)=cfa(jj)*(visc1(9)+ci*visc1(10))
    cuinv(jj)=cfa(jj)*(visc1(3)+ci*visc1(4))
    cbinv(jj)=cfa(jj)*(visc1(7)+ci*visc1(8))
    cpinv(jj)=cfa(jj)*(visc1(11)+ci*visc1(12))
    ccinv(jj)=cfa(jj)*(visc1(1)+ci*visc1(2)); endif

  enddo

endif  ! Fin de l'inviscid

! Uniform approximation:
cwuni=cwvis+cwinv-cwmat
cuuni=cuvis+cuinv-cumat
cbuni=cbvis+cbinv-cbmat
cpuni=cpvis+cpinv-cpmat
ccuni=ccvis+ccinv-ccmat
! Rotational
if(ibouss.eq.0.or.ibouss.eq.3)then ; cruni(:)=ccuni(:) ; else
cruni(:)=ccuni(:)-ci*xk(:)*cwuni(:); endif
! Displacement due to the waves:
cduni=0.;if(zz(ll).ge.zrid.and.ll.gt.2)cduni(:)=-ci*cwuni(:)/u0(ll)/xk(:)

!cwuni=cwvis;cuuni=cuvis;cbuni=cbvis

cwvis=cwuni ; call cfftb1(iim,cwvis,w1,w2,ifac);vitw(:,ll)=real(cwvis(:))
cwuni(1)=0. ; call cfftb1(iim,cwuni,w1,w2,ifac);wwav(:,ll)=real(cwuni(:))
cpvis=cpuni ; call cfftb1(iim,cpvis,w1,w2,ifac);pres(:,ll)=real(cpvis(:))
cpuni(1)=0. ; call cfftb1(iim,cpuni,w1,w2,ifac);pwav(:,ll)=real(cpuni(:))
cuvis=cuuni ; call cfftb1(iim,cuvis,w1,w2,ifac);vitu(:,ll)=real(cuvis(:))+u0(ll)
cuuni(1)=0. ; call cfftb1(iim,cuuni,w1,w2,ifac);uwav(:,ll)=real(cuuni(:))
cbvis=cbuni ; call cfftb1(iim,cbvis,w1,w2,ifac);buoy(:,ll)=real(cbvis(:))+zz(ll)*zj
cbuni(1)=0. ; call cfftb1(iim,cbuni,w1,w2,ifac);bwav(:,ll)=real(cbuni(:))
cxvis=cruni ; call cfftb1(iim,cxvis,w1,w2,ifac);vort(:,ll)=real(cxvis(:))
cruni(1)=0. ; call cfftb1(iim,cruni,w1,w2,ifac);vwav(:,ll)=real(cruni(:))
ccvis=ccuni ; call cfftb1(iim,ccvis,w1,w2,ifac);shea(:,ll)=real(ccvis(:))
ccuni(1)=0. ; call cfftb1(iim,ccuni,w1,w2,ifac);swav(:,ll)=real(ccuni(:))
cduni(1)=0. ; call cfftb1(iim,cduni,w1,w2,ifac);dwav(:,ll)=real(cduni(:))
acti(:,ll)=bwav(:,ll)*vwav(:,ll)/zj-u0zz(ll)*bwav(:,ll)**2/2./zj**2
fx(:,ll)=u0(ll)*acti(:,ll)+bwav(:,ll)**2/2./zj+uwav(:,ll)**2/2.
if(ibouss.eq.1.or.ibouss.eq.2)fx(:,ll)=fx(:,ll)-wwav(:,ll)**2/2.
fz(:,ll)=wwav(:,ll)*uwav(:,ll)

do ii=1,iim
if(zz(ll).gt.zh(ii))then
psi(ii,ll)=psi(ii,ll-1)+vitu(ii,ll)*(zz(ll)-zz(ll-1))
else;psi(ii,ll)=0.;endif
enddo

enddo    !  Fin boucle sur l'altitude....

! Surface values:
zmax=-100000.; zmin=+100000.;
do ii=1,iim
uzsol(ii)=shea(ii,1)+(shea(ii,2)-shea(ii,1))/(zz(2)-zz(1))*(zh(ii)-zz(1))
rsol(ii)=vort(ii,1)+(vort(ii,2)-vort(ii,1))/(zz(2)-zz(1))*(zh(ii)-zz(1))
psol(ii)=pres(ii,1)+(buoy(ii,2)+buoy(ii,1))/2.*(zh(ii)-zz(1))
fxsol(ii)=fx(ii,1)+(fx(ii,2)-fx(ii,1))/(zz(2)-zz(1))*(zh(ii)-zz(1))
fzsol(ii)=fz(ii,1)+(fz(ii,2)-fz(ii,1))/(zz(2)-zz(1))*(zh(ii)-zz(1))
do ll=1,llm-1
if(zz(ll).lt.zh(ii).and.zz(ll+1).ge.zh(ii)) then
uzsol(ii)=shea(ii,ll)+ &
 (zh(ii)-zz(ll))*(shea(ii,ll+1)-shea(ii,ll))/(zz(ll+1)-zz(ll))
rsol(ii)=vort(ii,ll)+ &
(zh(ii)-zz(ll))/(zz(ll+1)-zz(ll))*(vort(ii,ll+1)-vort(ii,ll))
psol(ii)=pres(ii,ll)+ &
(zh(ii)-zz(ll))*(buoy(ii,ll+1)+buoy(ii,ll))/2.
fxsol(ii)=fx(ii,ll)+ &
(zh(ii)-zz(ll))/(zz(ll+1)-zz(ll))*(fx(ii,ll+1)-fx(ii,ll))
fzsol(ii)=fz(ii,ll)+ &
(zh(ii)-zz(ll))/(zz(ll+1)-zz(ll))*(fz(ii,ll+1)-fz(ii,ll))
endif
enddo
enddo

!Theoretical drag:
tdrag=sqrt(zj)*zrid**2*zbl/2.    ! Theoretical value for wind intensity at zb
if(ibouss.ge.1) then
tdrag=zrid**2*zbl/2.*max(sqrt(zj),1.)
if(icond.eq.3) tdrag=zrid**3/2.*sqrt(zj)
endif
tgwdr=tdrag
if(icond.eq.3) tgwdr=zrid**4/2./zbl*sqrt(zj)

!Mountain drag
!mdrag=0;do ii=2,iim; mdrag=mdrag+(zh(ii)-zh(ii-1))* &
!(psol(ii)+psol(ii-1))/2. ; enddo ; mdrag=mdrag/tdrag
mdrag=0;do ii=2,iim; mdrag=mdrag-(zh(ii)+zh(ii-1))* &
(psol(ii)-psol(ii-1))/2. ; enddo ; mdrag=mdrag/tdrag

!Viscous drag
!actfl=0.;do ii=1,iim;actfl=actfl+uzsol(ii)/float(iim);enddo;uzsol(:)=uzsol(:)-actfl
vdrag=0.;do ii=2,iim;vdrag=vdrag+(uzsol(ii)+uzsol(ii-1))/2.*(xx(ii)-xx(ii-1));enddo
actfl=0.;do ii=2,iim; actfl=actfl-(zh(ii)-zh(ii-1))* &
(fxsol(ii)+fxsol(ii-1))/2.+(fzsol(ii)+fzsol(ii-1))/2.*(xx(ii)-xx(ii-1))
enddo ; vdrag=vdrag/tdrag*zbl**3 ; actfl=actfl/tdrag

!GWstress:
gwstr(:)=0.
do ll=2,llm; gwstr(ll)=0.;  do ii=2,iim
gwstr(ll)=gwstr(ll)+(xx(ii)-xx(ii-1))* &
(wwav(ii,ll)*uwav(ii,ll)+wwav(ii-1,ll)*uwav(ii-1,ll))/2.
enddo; enddo ; gwstr(:)=gwstr(:)/tdrag
do ll=llm-1,1,-1; if(zz(ll).le.zrid) gwstr(ll)=gwstr(ll+1) ; enddo
!Min  wave drag:
gwmin=0.;gwmax=0.;
do ll=1,llm; if(zz(ll).gt.zrid)then
gwmin=amin1(gwmin,gwstr(ll)); gwmax=amax1(gwmax,gwstr(ll)); endif; enddo 

! Storm index
storm=0. ; storp=0 ;  do ii=2,iim ; do ll=1,llm
if(xx(ii).gt.0.and.xx(ii).lt.2.and.zz(ll).lt.zrid*0.75.and.zz(ll).gt.zh(ii))&
storp=amax1(storp,sqrt(vitu(ii,ll)**2+vitw(ii,ll)**2))
!if(xx(ii).lt.0.and.xx(ii).gt.-2.and.zz(ll).lt.zrid)&
if(xx(ii).lt.0.and.xx(ii).gt.-2.and.zz(ll).lt.zrid*0.75.and.zz(ll).gt.zh(ii))&
storm=amax1(storm,sqrt(vitu(ii,ll)**2+vitw(ii,ll)**2))
enddo;enddo; down=storp/storm

if(.not.lout)then
write(44,*)zj,-mdrag     
write(45,*)zj,gwstr(llm)*tdrag/tgwdr
write(46,*)zj,-vdrag
!if(gwmax.gt.-gwmin)then;write(47,*)zj,gwmax*tdrag/tgwdr
!else;write(47,*)zj,gwmin*tdrag/tgwdr;endif    
write(47,*)zj,gwmin*tdrag/tgwdr
write(57,*)zj,gwmax*tdrag/tgwdr
write(48,*)zj,down  
write(49,*)zj,actfl*tdrag/tgwdr  
endif

if(lout)then
  print *,'Pressure drag:',mdrag    
  print *,'Viscous drag:',vdrag
  print *,'Wave drag:',gwstr(llm)
  print *,'Min gwdrag',gwmin   
  print *,'Max gwdrag',gwmax   
  print *,'Storm index=',down
  print *,'Surf Act flux=',actfl
  filename='../xmvg/yyyy_gwstress_xxxxxxxx.dat'
  write(filename(23:30),331)exppara;write(filename(9:12),330)exptype
  open(76,file=filename)
        do ll=1,llm;write(76,*)gwstr(ll),zz(ll);enddo;write(76,*)' '
! Other stresses:
! Pressure force acting through material surface:
gwpre=0.;do ll=3,llm-2; do ii=2,iim; gwpre(ll)=gwpre(ll)- &
(dwav(ii,ll)-dwav(ii-1,ll))* &
(pwav(ii,ll)+pwav(ii-1,ll))/2. ; enddo ; gwpre(ll)=gwpre(ll)/tdrag
enddo
do ll=3,llm-2;if(zz(ll).gt.zrid)write(76,*)gwpre(ll),zz(ll);enddo;write(76,*)' '
! Viscous force through material surface:
gwvis=0.;do ll=3,llm-2; do ii=2,iim; gwvis(ll)=gwvis(ll)-( &
 dwav(ii,ll)*((uwav(ii,ll+1)-uwav(ii,ll))/(zz(ll+1)-zz(ll))&
             -(uwav(ii,ll)-uwav(ii,ll-1))/(zz(ll)-zz(ll-1)))/(zz(ll+1)-zz(ll-1)) &
+dwav(ii+1,ll)*((uwav(ii+1,ll+1)-uwav(ii+1,ll))/(zz(ll+1)-zz(ll))&
       -(uwav(ii+1,ll)-uwav(ii+1,ll-1))/(zz(ll)-zz(ll-1)))/(zz(ll+1)-zz(ll-1)))&
*(xx(ii+1)-xx(ii))*zbl**3 ; enddo ; gwvis(ll)=gwvis(ll)/tdrag
enddo
do ll=3,llm-2;if(zz(ll-1).gt.zrid)write(76,*)gwvis(ll),zz(ll);enddo;write(76,*)' '
do ll=3,llm-2;if(zz(ll-1).gt.zrid)write(76,*)gwvis(ll)+gwpre(ll),zz(ll);enddo;write(76,*)' '

close(76);endif

! Theoretical values:
!Mountain drag
mdrag=0;do ii=2,iim; mdrag=mdrag+(zh(ii)-zh(ii-1))* &
(psth(ii)+psth(ii-1))/2. ; enddo ; mdrag=mdrag/tdrag

!Viscous drag
vdrag=0.;do ii=2,iim;vdrag=vdrag+(uzsth(ii)+uzsth(ii-1))/2.*(xx(ii)-xx(ii-1))
enddo ; vdrag=vdrag/tdrag*zbl**3 

if(lout)print *,'Theory for mountain: ',mdrag,'  viscous: ',vdrag
! Mountain drag dtrend!
pdtr=pres; call dtrend(pdtr,iim,llm)
if(lout) then
do ii=1,iim
psol(ii)=pdtr(ii,1)+(pdtr(ii,2)-pdtr(ii,1))/(zz(2)-zz(1))*zh(ii)
do ll=1,llm-1
if(zz(ll).le.zh(ii).and.zz(ll+1).gt.zh(ii)) then
psol(ii)=pdtr(ii,ll)+ &
(zh(ii)-zz(ll))/(zz(ll+1)-zz(ll))*(pdtr(ii,ll+1)-pdtr(ii,ll))
endif
enddo; enddo
! Mountain drag dtrend!
pdtr=shea; call dtrend(pdtr,iim,llm)
do ii=1,iim
uzsol(ii)=pdtr(ii,1)+(pdtr(ii,2)-pdtr(ii,1))/(zz(2)-zz(1))*zh(ii)
do ll=1,llm-1
if(zz(ll).le.zh(ii).and.zz(ll+1).gt.zh(ii)) then
uzsol(ii)=pdtr(ii,ll)+ &
(zh(ii)-zz(ll))/(zz(ll+1)-zz(ll))*(pdtr(ii,ll+1)-pdtr(ii,ll))
endif
enddo; enddo
!Mountain drag
mdrag=0;do ii=2,iim; mdrag=mdrag+(zh(ii)-zh(ii-1))* &
(psol(ii)+psol(ii-1))/2. ; enddo ; mdrag=mdrag/tdrag
!Viscous drag
vdrag=0.;do ii=2,iim;vdrag=vdrag+(uzsol(ii)+uzsol(ii-1))/2.*(xx(ii)-xx(ii-1))
enddo ; vdrag=vdrag/tdrag*zbl**3 
print *,'M and V drag detrended:',mdrag,vdrag
endif
if(lout)then
do ii=1,iim; write(91,*)xx(ii),uzsol(ii); enddo;  write(91,*)' '
do ii=1,iim; write(91,*)xx(ii),uzsth(ii); enddo;  write(91,*)' '
do ii=1,iim; write(91,*)xx(ii),psol(ii); enddo;  write(91,*)' '
do ii=1,iim; write(91,*)xx(ii),psth(ii); enddo;  write(91,*)' '
do ii=1,iim; write(91,*)xx(ii),psth(ii); enddo;  write(91,*)' '
do ii=1,iim; write(91,*)xx(ii),pdtr(ii,2); enddo;  write(91,*)' '
do ii=1,iim; write(91,*)xx(ii),shea(ii,2); enddo;  write(91,*)' '
endif

!***********************************
!  6.9: Lagrangian calculations
!***********************************

if(lout)then

filename='../xmvg/YYYY_Lagr_Stre_xxxxxxxx.dat'
write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
open(85,file=filename)

allocate(disp(iim,llm))

disp(1,:)=zz(:)
do ll=2,llm-1; do ii=1,iim
do kk=2,llm;if(psi(ii,kk).ge.psi0(ll).and.psi(ii,kk-1).le.psi0(ll)) then
disp(ii,ll)=zz(kk-1)+(psi0(ll)-psi(ii,kk-1))*(zz(kk)-zz(kk-1))/(psi(ii,kk)-psi(ii,kk-1));endif
enddo
enddo;enddo
do ii=1,iim; write(89,*)xx(ii),disp(ii,5);enddo

!Print *, 'Lagrangian continues'

do jj=1,3; do ll=2,llm-1

do ii=1,iim 
!psol(ii)=pres(ii,1)+(pres(ii,2)-pres(ii,1))/(zz(2)-zz(1))*(disp(ii,ll)-zz(1))
!uzsol(ii)=shea(ii,1)+(shea(ii,2)-shea(ii,1))/(zz(2)-zz(1))*(disp(ii,ll)-zz(1))
do kk=2,llm ; if(disp(ii,ll).le.zz(kk).and.disp(ii,ll).ge.zz(kk-1)) then
psol(ii)=pres(ii,kk-1)  + &
(disp(ii,ll)-zz(kk-1))/(zz(kk)-zz(kk-1))*(pres(ii,kk)-pres(ii,kk-1)) 
uzsol(ii)=shea(ii,kk-1)  + &
(disp(ii,ll)-zz(kk-1))/(zz(kk)-zz(kk-1))*(shea(ii,kk)-shea(ii,kk-1)) 
endif ;
enddo;enddo

!Mountain drag
mdrag=0;do ii=2,iim; mdrag=mdrag+(disp(ii,ll)-disp(ii-1,ll))* &
(psol(ii)+psol(ii-1))/2. ; enddo ; mdrag=mdrag/tdrag
vdrag=0.;do ii=2,iim;vdrag=vdrag+(uzsol(ii)+uzsol(ii-1))/2.*(xx(ii)-xx(ii-1));enddo
vdrag=vdrag/tdrag*zbl**3
!print *,' A mountain drag at:',disp(1,ll),mdrag
if(jj.eq.1)write(85,*)-mdrag,disp(1,ll)
if(jj.eq.2)write(85,*)-vdrag,disp(1,ll)
if(jj.eq.3)write(85,*)-mdrag-vdrag,disp(1,ll)

enddo;write(85,*)' ';enddo

close(85)

endif

! Vertical and horizontal PM Fluxes:
if(lout)then;   print *,'Action flux budget'
filename='../xmvg/YYYY_Fz_and_Fx_xxxxxxxx.dat'
write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
open(76,file=filename)
filename='../xmvg/YYYY_Fx_downSt_xxxxxxxx.dat'
write(filename(24:31),331)exppara;write(filename(9:12),330)exptype
open(77,file=filename)

kk=3
do jj=iim/2+1,3*iim/4; gwstr(:)=0.;gwvis(1)=0.
  do ll=2,llm; do ii=2,jj               
  gwstr(ll)=gwstr(ll)+(xx(ii)-xx(ii-1))*(fz(ii,ll)+fz(ii-1,ll))/2.
  enddo; enddo ; gwstr(:)=gwstr(:)/tdrag
  do ll=llm-1,1,-1; if(zz(ll).le.zrid) gwstr(ll)=gwstr(ll+1) ; enddo
  do ll=2,llm;gwvis(ll)=gwvis(ll-1)+(zz(ll)-zz(ll-1))*(fx(jj,ll-1)+fx(jj,ll))/2.
  enddo; gwvis(:)=gwvis(:)/tdrag
  if(xx(jj+1).ge.float(kk).and.xx(jj).lt.float(kk))then
  print *,'profiles in xx=',xx(jj);kk=kk+1
  do ll=1,llm;write(76,*)gwstr(ll),zz(ll);enddo;write(76,*)' '
  do ll=1,llm;write(76,*)gwvis(ll),zz(ll);enddo;write(76,*)' '
  do ll=1,llm;write(76,*)gwstr(ll)+gwvis(ll),zz(ll);enddo;write(76,*)' '
  endif
  write(77,*)gwstr(llm-1)+gwvis(llm-1)
enddo 

close(76);close(77); print *,'End of action  flux budget'
endif
! Same diags but over a 0.5x4 box
do ii=1,llm-1;if(zz(ii+1).gt..5.and.zz(ii).le..5) kk=ii+1;enddo
do ii=1,iim-1;if(xx(ii+1).gt.3..and.xx(ii).le.3.)   jp=ii+1;enddo
do ii=2,iim-1;if(xx(ii+1).gt.-3..and.xx(ii).le.-3.) jm=ii;enddo
gwstr(1)=0.;gwvis(1)=0.
do ii=jm,jp; gwstr(1)=gwstr(1)+(xx(ii)-xx(ii-1))*(fz(ii,kk)+fz(ii-1,kk))/2.; enddo
gwstr(1)=gwstr(1)/tdrag
do ll=2,kk;gwvis(1)=gwvis(1)+(zz(ll)-zz(ll-1))*(fx(jp,ll-1)+fx(jp,ll))/2.
           gwvis(1)=gwvis(1)-(zz(ll)-zz(ll-1))*(fx(jm,ll-1)+fx(jm,ll))/2.
enddo;gwvis(1)=gwvis(1)/tdrag
if (lout) then; print *,' PMF out of the box:',zj,gwstr(1)+gwvis(1) 
print *,'FZ and FX:',gwstr(1),gwvis(1); print *,zz(kk),xx(jm),xx(jp)
else;  write(50,*)zj,gwstr(1)+gwvis(1); endif 
!GWstress:

!***********************************
!  7: Storage
!***********************************

do ll=1,llm
wwav(:,ll)=wwav(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
vitw(:,ll)=vitw(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
uwav(:,ll)=uwav(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
bwav(:,ll)=bwav(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
vitu(:,ll)=vitu(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
buoy(:,ll)=buoy(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
pres(:,ll)=pres(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
acti(:,ll)=acti(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
vort(:,ll)=vort(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
fx(:,ll)=fx(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
fz(:,ll)=fz(:,ll)*(1.+sign(1.,zz(ll)-zh(:)))/2.
enddo


if(lout)then; filename='../idl_vg/xxxx_xxxxxxxx.dat'
write(filename(11:14),330)exptype;write(filename(16:23),331)exppara
if(lout) print *,'filename:',filename
!write(filename(15:22),331)exppara;331 format(a8)
open(60,file=filename)
write(60,*)zj; write(60,*)zrid;write(60,*)zbl;write(60,*)sld
write(60,*)iim;write(60,*)llm
do ii=1,iim;write(60,*)xx(ii);enddo;do ii=1,iim;write(60,*)zh(ii);enddo
do ll=1,llm;write(60,*)REAL(zz(ll));enddo
do ll=1,llm;do ii=1,iim;write(60,*)vitw(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)vitu(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)buoy(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)wwav(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)uwav(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)bwav(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)psi(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)pres(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)acti(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)fx(ii,ll);enddo;enddo
do ll=1,llm;do ii=1,iim;write(60,*)fz(ii,ll);enddo;enddo
close(60)

filename='../xmvg/wav_flds_xxxx_xxxxxxxx.dat'
write(filename(18:21),330)exptype;write(filename(23:30),331)exppara
if(lout) print *,'filename:',filename
!write(filename(15:22),331)exppara;331 format(a8)
open(60,file=filename)
ii=iim/2  
!do ll=1,llm;write(60,*)vitw(ii,ll),zz(ll);enddo;write(60,*)' '
!do ll=1,llm;write(60,*)vitu(ii,ll),zz(ll);enddo;write(60,*)' '
!do ll=1,llm;write(60,*)buoy(ii,ll),zz(ll);enddo;write(60,*)' '
do ll=1,llm;write(60,*)wwav(ii,ll),zz(ll);enddo;write(60,*)' '
do ll=1,llm;write(60,*)uwav(ii,ll),zz(ll);enddo;write(60,*)' '
do ll=1,llm;write(60,*)bwav(ii,ll),zz(ll);enddo;write(60,*)' '
!do ll=1,llm;write(60,*)psi(ii,ll),zz(ll);enddo;write(60,*)' '
do ll=1,llm;write(60,*)pres(ii,ll),zz(ll);enddo;write(60,*)' '
do ll=1,llm;write(60,*)shea(ii,ll),zz(ll);enddo;write(60,*)' '
do ll=1,llm;write(60,*)vort(ii,ll),zz(ll);enddo;write(60,*)' '
close(60)


!*******************************************************
!  This is the end, 
write(*,*)'Name: ',exptype,'_',exppara
write (*,*) 'Slope:',islope,slope(islope),zrid
write(*,*)'Zbl=',zbl
write(*,*)'ZJ=',zj

endif  ! End of storage

if(.not.lout)print *,'slope:',zrid,' Ric:',zj,' Shear layer depth:',sld

contains

!***********************************
! A: Subroutines
!***********************************
   subroutine test_neg_int(comp,logic)
!  Test that a complex is not to near a negative integer
   implicit none
   complex, intent(in)::  comp
   logical, intent(out):: logic
   real, parameter :: eps=0.002
!  real, parameter :: eps=0.01 
   integer, parameter :: itest=100
   integer :: it
   logic=.false.
   if(aimag(comp).ge.-eps/100..and.aimag(comp).le.eps/100.)then

     do it=1,itest
        if(real(comp).gt.-float(it)-eps.and.real(comp).lt.-float(it)+eps)logic=.true.
     enddo

   endif

   end subroutine test_neg_int

!  Detrend the pressure man!
   subroutine dtrend(field,iim,llm)
   implicit none
   integer, intent(in):: iim,llm
   real, intent(inout) :: field(iim,llm)
   integer:: ii,ll
   real :: xisum,pxsum,psum,det,pa,pb
!   iim=size(field,dim=1);llm=size(field,dim=2); print *,'In dtrend:',iim,llm
   do ll=1,llm;xisum=0.;pxsum=0.;psum=0.;det=0.
   do ii=1,iim
   xisum=xisum+float(ii)/float(iim);psum = psum+field(ii,ll)/float(iim)
   pxsum=pxsum+field(ii,ll)*float(ii)/float(iim);det=xisum*xisum-1.
   pa=(xisum*pxsum-psum)/det;pb=(-pxsum+psum*xisum)/det
   enddo
   do ii=1,iim; field(ii,ll)=field(ii,ll)-pa-pb*ii;enddo
   enddo
   end subroutine dtrend

   function hypsolold(aa,bb,cc,kk,cm,b1,b2,zj,zmu,zz)
   use nr
   implicit none
   complex :: hypsolold(3)
   complex, intent(in) :: aa,bb,cc,cm,b1,b2,zmu
   real, intent(in) :: kk,zz,zj  ! Hor. wavenumber, altitude, Rich number
   real ss,zp
   complex rr
   real, parameter :: eps=0.01
   complex :: visu,visd,visp,vism
   complex, parameter :: ci=(0.,1.)
  
   zp=zz+eps;ss=tanh(zp/sqrt(zj));rr=ss**2;visu=(0.,0.);visd=(0.,0.)
   if(ss.ge.0.05) visu=2**(-cm)*ss**(0.5+ci*zmu)*cosh(zp/sqrt(zj))**(-cm)*&
             hypgeo(cc-bb,cc-aa,cc-aa-bb+1.,1.-rr)
   if(ss.le.0.95) visd=cosh(zp/sqrt(zj))**cm*(b1*ss**(0.5-ci*zmu)& 
 * hypgeo(aa-cc+1,bb-cc+1,2.-cc,rr) + b2*ss**(0.5+ci*zmu)*hypgeo(aa,bb,cc,rr))   

   hypsolold(1)=visu*(1.+tanh((zp/sqrt(zj)-0.5)/0.1))/2.&
            +visd*(1.-tanh((zp/sqrt(zj)-0.5)/0.1))/2.
   hypsolold(2)=visu*(1.+tanh((zp/sqrt(zj)-0.5)/0.1))/2.&
            +visd*(1.-tanh((zp/sqrt(zj)-0.5)/0.1))/2.

   zp=zz-eps;ss=tanh(zp/sqrt(zj));rr=ss**2
   if(ss.ge.0.05) visu=2**(-cm)*ss**(0.5+ci*zmu)*cosh(zp/sqrt(zj))**(-cm)*&
             hypgeo(cc-bb,cc-aa,cc-aa-bb+1.,1.-rr)
   if(ss.le.0.95) visd=cosh(zp/sqrt(zj))**cm*(b1*ss**(0.5-ci*zmu)& 
 * hypgeo(aa-cc+1,bb-cc+1,2.-cc,rr) + b2*ss**(0.5+ci*zmu)*hypgeo(aa,bb,cc,rr))   

   hypsolold(1)=(hypsolold(1)+visu*(1.+tanh((zp/sqrt(zj)-0.5)/0.1))/2.&
             +visd*(1.-tanh((zp/sqrt(zj)-0.5)/0.1))/2.)/2.
   hypsolold(2)=ci*(hypsolold(2)-visu*(1.+tanh((zp/sqrt(zj)-0.5)/0.1))/2.&
             -visd*(1.-tanh((zp/sqrt(zj)-0.5)/0.1))/2.)/2./eps/kk
   hypsolold(3)=ci*hypsolold(1)/kk/tanh(max(zz,eps)/sqrt(zj))

   if(zz.le.eps)print *,'Its not very safe to use invis when z<',eps
   
   end function hypsolold
!inv=hypsol(aa(jj),bb(jj),cc,xk(jj),cm,b1(jj),b2(jj),zmu,zz(ll))

end program mount_wave

