subroutine mixsstz(f,stm,h,fsm,zlev,im,jm,nl,tbias)
parameter(kk=9,jj=31,step=5.)
c---------------------------------------------------------------------
c----- this subr. was rewritten 06/17/05 by a.falkovich
c----- to assimilate sst real data in gdem monthly t data for z-levels
c---------------------------------------------------------------------
c
c---- change temperature in the layer 0-125m.
c---- use spline to interpolate from gdem z-levels to even step 10m
c---- here temperature is specified at all levels
c---- (under land and under bottom)
c---- mixed layer depth is the last level n where t(1)-t(k)<0.5deg
c---- dts=stm-t(1); tnew(k)=t(k)+dts for k=1,...,n
c---- if abs(dts)>5deg correct stm
c---- if stm<t(150m) correct stm
c
real stm(im,jm),h(im,jm),fsm(im,jm)
real f(im,jm,nl),zlev(nl)
real dts,t(jj),z(jj),tt(kk),zz(kk),tc(jj),te(kk),dt(kk)
logical plg1,plg2,plg3
c
c--------- stm is used after correction for tbias,
c--------- f temperature before correction for tbias
c
c
c---------------
print *,' mixsstz: im,jm,nl=',im,jm,nl
print *,' tbias=',tbias
print *,'kk,jj=',kk,jj
print *,' zlev, k=1,nl'
c---------------
write(6,202) (zlev(k),k=1,nl)
202 format(10f7.0)
c--------------- for debug
i0=124
j0=92
c--------------- for debug
i1=124
j1=88
j2=98
c---------------
print *,' before sst assim: i1,j1,j2=',i1,j1,j2
print *,' f(i1,j,k), j=j1,j2 nl=',nl
write(6,301) (j,j=j1,j2)
do k=1,kk
write(6,302) zlev(k),(f(i1,j,k),j=j1,j2)
end do
301 format(5x,11i7)
302 format(f8.0,11f7.2)
c---------------
c----- remove correction for tbias
do j=1,jm
do i=1,im
if(fsm(i,j).eq.1.) stm(i,j)=stm(i,j)+tbias
end do
end do
c----------------
write(6,401) i0,j0,stm(i0,j0),h(i0,j0)
401 format('i0,j0,stm(i0,j0),h(i0,j0)=',
* /2i7,f7.2,f7.0)
c---------------
print *,' stm(i1,j), j=j1,j2 '
write(6,303) (stm(i1,j),j=j1,j2)
303 format(8x,11f7.2)
c
c------------------ check difference between stm and f(i,j,1)
c
difmax=0.
difmin=0.
do i=1,im
do j=1,jm
if(fsm(i,j).gt.0.5) then
difp=stm(i,j)-f(i,j,1)
difm=-difp
if(difp.gt.difmax) then
imax=i
jmax=j
difmax=difp
end if
if(difm.gt.difmin) then
imin=i
jmin=j
difmin=difm
end if
end if
end do
end do
c------------------
print *,' mixsstz: max(stm-f(1)),i,j'
print *,' difmax,idifmax,jdifmax=',difmax,imax,jmax
print *,' mixsstz: max(f(1)-stm),i,j'
print *,' difmin,imin,jmin=',difmin,imin,jmin
c------------------ specify z and zz
do j=1,jj
z(j)=step*float(j-1)
end do
do k=1,kk
zz(k)=zlev(k)
end do
c------------------
c----------------------------- loops in i,j
do j=1,jm
do i=1,im
c--------------- for debug
c plg1=i.eq.i0.and.j.eq.j0
plg1=i.eq.i1.and.j.ge.j1.and.j.le.j2
plg2=i.eq.imax.and.j.eq.jmax
plg3=i.eq.imin.and.j.eq.jmin
c---------------
if(fsm(i,j).eq.0.) go to 1000
c------------ send tt
do k=1,kk
tt(k)=f(i,j,k)
end do
c------------ send sst
sst=stm(i,j)
c------------ sp interpolate to even step z
call sp(zz,tt,t,z,kk,jj,kk,jj)
c------------
c if(i.eq.i0.and.j.eq.j0) then
if(plg1.or.plg2.or.plg3) then
c------------
print *,'i,j,sst=',i,j,sst
print *,'initial depth zz'
write(6,101) zz
101 format(10f7.0)
c------------
print *,'initial temperature tt'
write(6,102) tt
102 format(10f7.2)
c------------
print *,' depth for interpolation'
write(6,101) z
c------------
print *,'temperature after spline interpolation'
write(6,102) t
c------------
end if
c
c------------- find n for the depth of mixed layer
n=1
do while ((t(1)-t(n+1)).le.0.5.and.n.lt.(jj))
n=n+1
end do
c------------- if mixed layer more than 80m send 80m (n=17)
if(n.gt.17) n=17
c---------------------- no need to do correction
c---------------------- correct sst if it deviates too much from t(1)
c if(sst-t(1).gt.5.) sst=t(1)+5.
c if(sst-t(1).lt.-5.) sst=t(1)-5.
c---------------------- correct sst if sst<t(jj)
c if(sst.lt.t(jj)) sst=t(jj)
c----------------------
dts=sst-t(1)
c----------------------
do k=1,n
tc(k)=t(k)+dts
end do
c----------------------
m=jj-n
c----------------------
do k=n+1,jj
tc(k)=t(k)+dts*float(jj-k)/float(m)
end do
c----------------------
c if(i.eq.i0.and.j.eq.j0) then
c if(plg2.or.plg3) then
if(plg1.or.plg2.or.plg3) then
print *,'i,j,sst=',i,j,sst
print *,'mix: n,m=',n,m
write(6,201) sst,t(1),dts
201 format('sst,t(1),dts=',3f7.2)
c------------
print *,' depth for interpolation'
write(6,101) z
c------------
print *,'temperature in mix before checking'
write(6,102) tc
end if
c------------- check if(tc(k).lt.t(jj))
do k=2,jj-1
if(tc(k).le.t(jj)) then
n=k-1
m=jj-n
dts=t(jj)-tc(n)
c write(6,104) n,m,dts
104 format('n,m,dts=',2i7,f7.2)
go to 100
end if
end do
go to 200
100 continue
do k=n+1,jj
tc(k)=t(jj)-dts*float(jj-k)/float(m)
end do
c------------
200 continue
c------------- check stability
do k=2,jj
if(tc(k-1).lt.tc(k)) tc(k)=tc(k-1)
end do
c------------ interpolate back to gdem z-levels
call sp(z,tc,te,zz,jj,kk,jj,kk)
c------------
do k=1,kk
dt(k)=te(k)-tt(k)
end do
c------------
c if(i.eq.i0.and.j.eq.j0) then
c if(plg2.or.plg3) then
if(plg1.or.plg2.or.plg3) then
print *,'i,j,sst=',i,j,sst
print *,'temperature in mix after checking'
write(6,102) tc
c------------
print *,'initial depth zz'
write(6,101) zz
c------------
print *,'initial temperature tt'
write(6,102) tt
c------------
print *,'temperature after back interpolation'
write(6,102) te
c------------
print *,'temperature differences after mix'
write(6,102) dt
end if
c---------------------- send back adjusted temperature
do k=1,kk
f(i,j,k)=te(k)
end do
c------------
1000 continue
c--------------------- end loops in i,j
end do
end do
c----------------------
print *,' after sst assim: i1,j1,j2=',i1,j1,j2
print *,' f(i1,j,k), j=j1,j2 nl=',nl
write(6,301) (j,j=j1,j2)
do k=1,kk
write(6,302) zlev(k),(f(i1,j,k),j=j1,j2)
end do
c----------------------
return
end