MODULE mrls
USE mkinds
USE mstate
USE mpobs
USE mutil
USE mens
USE mspmd
USE mgrid

IMPLICIT NONE

! Assimilation Parameters

integer(i4),parameter   :: navars=10
REAL(r4)                :: acdepth(navars)
REAL(r4)                :: arlocal(navars)
REAL(r4)                :: acovscl(navars)
character(len=3)        :: atype

REAL(r4)                :: cdpth
REAL(r4)                :: rlocal
REAL(r4)                :: covscl

REAL(r4), ALLOCATABLE   :: hxf(:)
REAL(r4), ALLOCATABLE   :: dmhxf(:)
REAL(r4), ALLOCATABLE   :: oocov(:,:)

CONTAINS

subroutine do_analysis_mgrid(var,ebkg,eobs)
character(:),allocatable :: fname
character(3)             :: var
integer(i4)              :: io,jo,iw,jw,nw,mw,k,m,i,j,maxw,n,ip,jp,ncount,ii,jj,z,n_obs,funit,fid
real(r4),allocatable     :: xa(:,:),w(:),xf(:,:)
real(r4)                 :: xm,bm,wm,var_obs,var_bkg,grdi,grdj
character(:),allocatable :: ncfile
integer(i4)              :: winl,den,l,ivar,mxden
real(r4),allocatable     :: ob_i(:),ob_j(:),dmhxf(:),cov(:,:),mlon(:),mlat(:),varb(:),vard(:)
real(r4)                 :: rlocal,rdist,del,hscl,spval,corr,cdpth,covscl,localr
LOGICAL(bl)              :: fail
real(r4)                 :: ebkg(:,:),eobs(:,:)

CALL get_analysis_params(var,cdpth,rlocal,covscl)



spval=2.00**100  ! fill value
!vard=0.04*0.04   ! variance of the observation errors SLA
!varb=0.125*0.125 ! background/forecast error SLA variance 
!rlocal=24.0     ! localization radius in grid points 40 grid points ~ 160km
winl=rlocal      ! observation search window length 
!n_obs=319       ! total no of observations 

!! allocate necessary arrays

mxden=(2*winl+1)*(2*winl+1)
allocate(dmhxf(mxden),ob_i(mxden),ob_j(mxden),vard(mxden),varb(mxden))


!! now loop over all model grid points calculating an analysis based on
!observations within a box specified by window length


do j=jfp,jlp
 jo=j-winl ; jw =j+winl
 jo=MAX(1,jo)  ; jw=MIN(ny,jw)
  do i=ifp,ilp
  if (mdepth(i,j)<spval .and. mdepth(i,j) >=cdpth) then
   io=i-winl ;  iw=i+winl
   io=MAX(1,io)  ; iw=MIN(nx,iw)

   den=0
   dmhxf=0.0
   ob_i=0.
   ob_j=0.
   del=0.0

!! calculate the no of data points for the model grid point (i,j)   

   do jj=jo,jw
     do ii=io,iw
        if(b(ii,jj,1)<0.5*spval) then
         den=den+1
         ob_i(den)=real(ii,r4)     !data locations on the model grid 
         ob_j(den)=real(jj,r4)     !data locations on the model grid
         rdist=sqrt((i-ob_i(den))**2+(j-ob_j(den))**2) ! distance of the data point to  grid pt (i,j)
         !if(rdist==0.0) rdist=0.1
         dmhxf(den)=b(ii,jj,1)*gslocz(rdist/rlocal) ! innovations (obs-forecast) tapered with localization
         varb(den)=ebkg(ii,jj)      ! background std
         vard(den)=eobs(ii,jj)         ! obs std dev
         endif
      enddo
    enddo

    if(den>0) then
! build the covariance matirx - max size is (mxden,mxden) - local aproximation
allocate(cov(den,den))
DO l=1,den
 DO m=1,den
  IF (l==m) THEN
    cov(l,m) = vard(l)*vard(m)+varb(l)*varb(m)      ! (R+HBHT)
!     cov(l,m) = vard(den)*vard(den) + SUM(ens%var2d(int(ob_i(l)),int(ob_j(l)),:)*ens%var2d(int(ob_i(l)),int(ob_j(l)),:))/REAL(ens%nmembers-1,r4)
  ELSE
    rdist=sqrt((ob_i(l)-ob_i(m))**2+(ob_j(l)-ob_j(m))**2)
    !!if(rdist==0.0) rdist=0.1
    hscl=rdist/rlocal         ! scaled length to the observations - rlocal can be spatially variable  
    corr=(1+hscl)*exp(-hscl)  ! correlation based on second order auto regressive functions - compact support
    cov(l,m) = varb(l)*varb(m)*corr      ! (HBHT)
    !cov(l,m) = SUM(ens%var2d(int(ob_i(l)),int(ob_j(l)),:)*ens%var2d(int(ob_i(m)),int(ob_j(m)),:))/REAL(ens%nmembers-1,r4)
  ENDIF
!  print *,l,m,cov(l,m),scl(l),scl(m)
 ENDDO
ENDDO

!compute the Cholesky factorization of a real symmetric positive definite covariance matrix 

CALL SPOTRF ('U',den,cov,den,fail)
      IF(fail) THEN
      PRINT *, "chol"
      ENDIF
!solve a system of linear equations A*X = B with a symmetric positive definite
!matrix A using the Cholesky fac- torization A = U**T*U

CALL SPOTRS ('U',den,1,cov,den,dmhxf(1:den),den, fail)
      IF(fail) THEN
      PRINT *,"lu"
      ENDIF

! compute corrections to grid point(i,j) based on nearby observations      
   do l=1,den
              rdist=sqrt((i-ob_i(l))**2+(j-ob_j(l))**2)
              !!if(rdist==0.0) rdist=0.1
              hscl=rdist/rlocal
!              print *,scl(l),rossby(i,j),rdist,hscl
              corr=(1+hscl)*exp(-hscl) ! correlations between observation points and grid point
              del=del+dmhxf(l)*varb(l)*ebkg(i,j)*corr
              !del=del+dmhxf(l)*SUM(ens%var2d(int(ob_i(l)),int(ob_j(l)),:)*ens%var2d(i,j,:))/REAL(ens%nmembers-1,r4)
   enddo
work2d(i,j,1)=del
deallocate(cov)



endif  ! if non zero obs



endif  ! model i,j, in water
enddo  !j
enddo  !i
!print *, pid, "....finished analysis...." 
end subroutine 

subroutine do_analysis_mgrid_ens(var,ebkg,eobs)
character(:),allocatable :: fname
character(3)             :: var
integer(i4)              :: io,jo,iw,jw,nw,mw,k,m,i,j,maxw,n,ip,jp,ncount,ii,jj,z,n_obs,funit,fid
real(r4),allocatable     :: xa(:,:),w(:),xf(:,:)
real(r4)                 :: xm,bm,wm,var_obs,var_bkg,grdi,grdj
character(:),allocatable :: ncfile
integer(i4)              :: winl,den,l,ivar,mxden
real(r4),allocatable     :: ob_i(:),ob_j(:),dmhxf(:),cov(:,:),mlon(:),mlat(:),varb(:),vard(:)
real(r4)                 :: rlocal,rdist,del,hscl,spval,corr,cdpth,covscl,localr
LOGICAL(bl)              :: fail
real(r4)                 :: ebkg(:,:),eobs(:,:)

CALL get_analysis_params(var,cdpth,rlocal,covscl)



spval=2.00**100  ! fill value
!vard=0.04*0.04   ! variance of the observation errors SLA
!varb=0.125*0.125 ! background/forecast error SLA variance 
!rlocal=24.0     ! localization radius in grid points 40 grid points ~ 160km
winl=rlocal      ! observation search window length
!n_obs=319       ! total no of observations 

!! allocate necessary arrays

mxden=(2*winl+1)*(2*winl+1)
allocate(dmhxf(mxden),ob_i(mxden),ob_j(mxden),vard(mxden),varb(mxden))


!! now loop over all model grid points calculating an analysis based on
!observations within a box specified by window length


do j=jfp,jlp
 jo=j-winl ; jw =j+winl
 jo=MAX(1,jo)  ; jw=MIN(ny,jw)
  do i=ifp,ilp
  if (mdepth(i,j)<spval .and. mdepth(i,j) >=cdpth) then
   io=i-winl ;  iw=i+winl
   io=MAX(1,io)  ; iw=MIN(nx,iw)

   den=0
   dmhxf=0.0
   ob_i=0.
   ob_j=0.
   del=0.0

!! calculate the no of data points for the model grid point (i,j)   

   do jj=jo,jw
     do ii=io,iw
        if(b(ii,jj,1)<0.5*spval) then
         den=den+1
         ob_i(den)=real(ii,r4)     !data locations on the model grid 
         ob_j(den)=real(jj,r4)     !data locations on the model grid
         rdist=sqrt((i-ob_i(den))**2+(j-ob_j(den))**2) ! distance of the data point to  grid pt (i,j)
         !if(rdist==0.0) rdist=0.1
         dmhxf(den)=b(ii,jj,1)*gslocz(rdist/rlocal) ! innovations (obs-forecast) tapered with localization
         varb(den)=ebkg(ii,jj)      ! background std
         vard(den)=eobs(ii,jj)         ! obs std dev
         endif
      enddo
    enddo

    if(den>0) then
! build the covariance matirx - max size is (mxden,mxden) - local aproximation
allocate(cov(den,den))
DO l=1,den
 DO m=1,den
  IF (l==m) THEN
!    cov(l,m) = vard(l)*vard(m)+varb(l)*varb(m)      ! (R+HBHT)
     cov(l,m) = vard(den)*vard(den) + SUM(ens%var2d(int(ob_i(l)),int(ob_j(l)),:)*ens%var2d(int(ob_i(l)),int(ob_j(l)),:))/REAL(ens%nmembers-1,r4)
  ELSE
    !rdist=sqrt((ob_i(l)-ob_i(m))**2+(ob_j(l)-ob_j(m))**2)
    !!if(rdist==0.0) rdist=0.1
    !hscl=rdist/rlocal         ! scaled length to the observations - rlocal can be spatially variable  
    !corr=(1+hscl)*exp(-hscl)  ! correlation based on second order auto regressive functions - compact support
!    cov(l,m) = varb(l)*varb(m)*corr      ! (HBHT)
     cov(l,m) = SUM(ens%var2d(int(ob_i(l)),int(ob_j(l)),:)*ens%var2d(int(ob_i(m)),int(ob_j(m)),:))/REAL(ens%nmembers-1,r4)
  ENDIF
!  print *,l,m,cov(l,m),scl(l),scl(m)
 ENDDO
ENDDO

!compute the Cholesky factorization of a real symmetric positive definite covariance matrix 

CALL SPOTRF ('U',den,cov,den,fail)
      IF(fail) THEN
      PRINT *, "chol"
      ENDIF
!solve a system of linear equations A*X = B with a symmetric positive definite
!matrix A using the Cholesky fac- torization A = U**T*U

CALL SPOTRS ('U',den,1,cov,den,dmhxf(1:den),den, fail)
      IF(fail) THEN
      PRINT *,"lu"
      ENDIF

! compute corrections to grid point(i,j) based on nearby observations      
   do l=1,den
              !rdist=sqrt((i-ob_i(l))**2+(j-ob_j(l))**2)
              !!if(rdist==0.0) rdist=0.1
              !hscl=rdist/rlocal
!              print *,scl(l),rossby(i,j),rdist,hscl
               !corr=(1+hscl)*exp(-hscl) ! correlations between observation points and grid point
     !         del=del+dmhxf(l)*varb(l)*ebkg(i,j)*corr
              del=del+dmhxf(l)*SUM(ens%var2d(int(ob_i(l)),int(ob_j(l)),:)*ens%var2d(i,j,:))/REAL(ens%nmembers-1,r4)
   enddo
work2d(i,j,1)=del
deallocate(cov)



endif  ! if non zero obs



endif  ! model i,j, in water
enddo  !j
enddo  !i
print *, pid, "....finished analysis...." 
end subroutine 


subroutine do_analysis_mgrid_uvp(var,ebkg,eobs)
character(:),allocatable :: fname
character(3)             :: var
integer(i4)              :: io,jo,iw,jw,nw,mw,k,m,i,j,maxw,n,ip,jp,ncount,ii,jj,z,n_obs,funit,fid,nvalid
real(r4),allocatable     :: xa(:,:),w(:),xf(:,:)
real(r4)                 :: xm,bm,wm,var_obs,var_bkg,grdi,grdj
character(:),allocatable :: ncfile
integer(i4)              :: winl,den,l,ivar,mxden
real(r4),allocatable     :: ob_i(:),ob_j(:),dmhxf(:),cov(:,:),varb(:),vard(:)
real(r4)                 :: rlocal,rdist,del,hscl,spval,corr,cdpth,covscl,localr
LOGICAL(bl)              :: fail
real(r4)                 :: ebkg(:,:,:),eobs(:,:,:)
real(r4)                 :: r2,xx,yy,xy,coruu,corvv,coruv!,varb,vard
real(r4)                 :: pi,f,f0,scl,g,eps

CALL get_analysis_params(var,cdpth,rlocal,covscl)

!
pi     = 2. * asin(1.0)
f0     = 4. * pi / 24. / 3600.
g      = 9.8
scl    = 6000.0 
eps    =0.001
! analysis parameters

spval=2.00**100   ! fill value
!vard=0.05         ! variance of the observation errors SLA
!varb=0.25         ! background/forecast error SLA variance
winl=12!rlocal     ! observation search window length 
rlocal=12         ! localization radius

mxden=2*(2*winl+1)*(2*winl+1)
allocate(dmhxf(mxden),ob_i(mxden),ob_j(mxden),vard(mxden),varb(mxden))


!! now loop over all model grid points calculating an analysis based on
!observations within a box specified by window length
cdpth=0.0

do j=jfp,jlp
do i=ifp,ilp
  if (mdepth(i,j)<spval) then
   jo=j-winl ; jw =j+winl
   jo=MAX(1,jo)  ; jw=MIN(ny,jw)
   io=i-winl ;  iw=i+winl
   io=MAX(1,io)  ; iw=MIN(nx,iw)


   nvalid=0
   den=0
   dmhxf=0.0
   ob_i=0.
   ob_j=0.
   del=0.0

!! calculate the no of data points for the model grid point (i,j)   

   do jj=jo,jw
     do ii=io,iw
        if(b(ii,jj,1)<0.5*spval) then
         den=den+1
         ob_i(den)=real(ii,r4)     !data locations on the model grid 
         ob_j(den)=real(jj,r4)     !data locations on the model grid
         rdist=sqrt((i-ob_i(den))**2+(j-ob_j(den))**2) ! distance of the data point to  grid pt (i,j)
         dmhxf(den)=b(ii,jj,1)*gslocz(rdist/rlocal) ! innovations (obs-forecast) tapered with localization
         varb(den)=ebkg(ii,jj,1)      ! background std
         vard(den)=eobs(ii,jj,1)         ! obs std dev
         endif
     enddo
  enddo
 
do jj=jo,jw
       do ii=io,iw
         if(b(ii,jj,2)<0.5*spval) then
         den=den+1
         ob_i(den)=real(ii,r4)     !data locations on the model grid 
         ob_j(den)=real(jj,r4)     !data locations on the model grid
         rdist=sqrt((i-ob_i(den))**2+(j-ob_j(den))**2) ! distance of the data point to  grid pt (i,j)
         dmhxf(den)=b(ii,jj,2)*gslocz(rdist/rlocal) ! innovations (obs-forecast) tapered with localization
         varb(den)=ebkg(ii,jj,2)      ! background std
         vard(den)=eobs(ii,jj,2)         ! obs std dev
         endif
      enddo
    enddo




if(den>0) then
! build the covariance matirx - max size is (mxden,mxden) - local aproximation
allocate(cov(den,den))
DO l=1,den
 DO m=1,den

   xx=(ob_i(l)-ob_i(m))**2 + eps
   yy=(ob_j(l)-ob_j(m))**2 + eps
   xy=(ob_i(l)-ob_i(m))*(ob_j(l)-ob_j(m)) + eps
   rdist=sqrt(xx+yy) + eps

   coruu=FG(rdist,0.0) + (FF(rdist,0.0)-FG(rdist,0.0))*xx/(rdist*rdist)
   corvv=FG(rdist,0.0) + (FF(rdist,0.0)-FG(rdist,0.0))*yy/(rdist*rdist)
   coruv=(FF(rdist,0.0)-FG(rdist,0.0))*xy/(rdist*rdist)

   if(l<=den/2 .and. m<=den/2 ) then
           if(l==m) then
           cov(l,m)=vard(l)*vard(m)+varb(l)*varb(m)*coruu
            !cov(l,m)=vard*vard+varb*varb*coruu
           else
           cov(l,m)=varb(l)*varb(m)*coruu
           ! cov(l,m)=varb*varb*coruu
           endif

 else if(l>den/2 .and. m>den/2) then
           if(l==m) then
           cov(l,m)=vard(l)*vard(m)+varb(l)*varb(m)*corvv
           !cov(l,m)=vard*vard+varb*varb*corvv
           else
           !cov(l,m)=varb*varb*corvv
           cov(l,m)=varb(l)*varb(m)*corvv
           endif
 else
           cov(l,m)=varb(l)*varb(m)*coruv
 endif

  ENDDO
ENDDO

!compute the Cholesky factorization of a real symmetric positive definite covariance matrix 

CALL SPOTRF ('U',den,cov,den,fail)
      IF(fail) THEN
      PRINT *, "chol"
      ENDIF
!solve a system of linear equations A*X = B with a symmetric positive definite
!matrix A using the Cholesky fac- torization A = U**T*U

CALL SPOTRS ('U',den,1,cov,den,dmhxf(1:den),den, fail)
      IF(fail) THEN
      PRINT *,"lu"
      ENDIF

! compute corrections to grid point(i,j) based on nearby observations      
   do l=1,den
              f = f0 * sin (mlat(i,j) * pi / 180. )
              xx=(i-ob_i(l))**2 + eps
              yy=(j-ob_j(l))**2 + eps
              xy=(i-ob_i(l))*(j-ob_j(l)) + eps 
              rdist=sqrt(xx+yy) + eps
              r2=rdist*rdist
              coruu=FG(rdist,0.0) + (FF(rdist,0.0)-FG(rdist,0.0))*xx/r2
              corvv=FG(rdist,0.0) + (FF(rdist,0.0)-FG(rdist,0.0))*yy/r2
              coruv=(FF(rdist,0.0)-FG(rdist,0.0))*xy/r2
              if(l<=den/2) then
               work2d(i,j,1)=work2d(i,j,1)+dmhxf(l)*varb(l)*eobs(i,j,1)*coruu
               work2d(i,j,2)=work2d(i,j,2)+dmhxf(l)*varb(l)*eobs(i,j,2)*coruv
               work2d(i,j,3)=work2d(i,j,3) + dmhxf(l)*varb(l)*eobs(i,j,3)*FF(rdist,0.0)*sqrt(yy)*scl*f/g
              else
               work2d(i,j,1)=work2d(i,j,1)+dmhxf(l)*varb(l)*eobs(i,j,1)*coruv
               work2d(i,j,2)=work2d(i,j,2)+dmhxf(l)*varb(l)*eobs(i,j,2)*corvv
               work2d(i,j,3)=work2d(i,j,3) - dmhxf(l)*varb(l)*eobs(i,j,3)*FF(rdist,0.0)*sqrt(xx)*scl*f/g
              endif

     enddo
deallocate(cov)

endif  ! if non zero obs
endif  ! model i,j, in water
!print *, i,j,del
enddo  !j
enddo  !i


end subroutine do_analysis_mgrid_uvp

subroutine do_gridded_mgrid(var,b,ebkg,eobs)
character(len=*)      :: var
real(r4)              :: b(:,:),ebkg(:,:),eobs(:,:)
INTEGER(i4)           :: i,J
REAL(r4)              :: varb,vard
REAL(r4)              :: wgt

CALL get_analysis_params(var,cdpth,rlocal,covscl)

DO j=1,ny  !jfp,jlp
DO  i=1,nx !ifp,ilp
  IF (fcst(i,j,1) <2.00**99 .and. mdepth(i,j) .ge. cdpth) then

           varb=ebkg(i,j)*ebkg(i,j)
           vard=eobs(i,j)*eobs(i,j)

          if(varb >0. .and. vard >0.) then
                 wgt=varb/(varb+vard)
                 work2d(i,j,1)=wgt*(b(i,j)-fcst(i,j,1))
          endif
  endif
 enddo
enddo

end subroutine        

REAL(KIND=r4) FUNCTION gslocz(r)
!Gaspari-Cohn taper function. r should be positive, and normalized so taper = 0 at r = 1
!very close to exp(-(r/c)**2), where c = 0.388. See -- Gaspari and Cohn, 1999, QJRMS, p. 723-757, eqn 4.10.
!2011-04-30  Initial version.

IMPLICIT NONE
REAL(r4), INTENT(IN) :: r
REAL(r4)             :: rr

 rr = 2.*r

 IF(r <= 0.5) THEN

    gslocz = ( ( ( -0.25*rr +0.5 )*rr +0.625 )*rr -5.0/3.0 )*rr**2 + 1.0

 ELSE IF (r < 1.0) THEN

    gslocz = ( ( ( ( rr/12.0 -0.5 )*rr +0.625 )*rr +5.0/3.0 )*rr -5.0 )*rr &
            + 4.0 - 2.0 / (3.0 * rr)
 ELSE

    gslocz = 0.

 ENDIF

END FUNCTION gslocz

REAL(KIND=r4) FUNCTION FF(r,t)
    use mkinds
    IMPLICIT NONE
    REAL(r4), INTENT(IN) :: r,t
    REAL(r4)             :: rr
    REAL(r4),parameter   :: a=1./8.0
    REAL(r4),parameter   :: TT=10.0
    REAL(r4)             :: hscl,rlocal

    hscl=a*r
    FF=(1. + hscl -0.25*(hscl*hscl))*exp(-hscl)!*exp((-t/TT)**2)

END FUNCTION FF

REAL(KIND=r4) FUNCTION FG(r,t)
    use mkinds
    IMPLICIT NONE
    REAL(r4), INTENT(IN) :: r,t
    REAL(r4)             :: rr
    REAL(r4),parameter   :: a=1./8.0
    REAL(r4),parameter   :: TT=10.0
    REAL(r4)             :: hscl,rlocal
    hscl=a*r
    FG=(1. + hscl -7.0*0.25*(hscl*hscl) +0.25*(hscl)**3)*exp(-hscl)!*exp((-t/TT)**2)

END FUNCTION FG

REAL(KIND=r4) FUNCTION FE(r,t)
    use mkinds
    IMPLICIT NONE
    REAL(r4), INTENT(IN) :: r,t
    REAL(r4)             :: rr
    REAL(r4),parameter   :: c1=1.0/.0,c2=0.9,a=1./8.
    REAL(r4),parameter   :: TT=10.0
    REAL(r4)             :: hscl,rlocal

    rlocal=8.
    hscl=r/rlocal

    FE=(1+hscl)*exp(-hscl)

    FE=(1.+a*r+((a*r)**2)/6.0-((a*r)**3)/6.0)*exp(-a*r)
    !FE= 1-c2+c2*(1+c1*r)*exp(-c1*r)

END FUNCTION FE

SUBROUTINE adjustStateLilo()
USE mhycom
USE mgrid

IMPLICIT NONE
!print *, minval(work2d(:,:,3)), maxval(work2d(:,:,3)),sum(ists)
CALL lilo(start,count,imask,ists,mdepth,cdpth,thref,thbase,dp0k,athk,aden,fmld,work2d(:,:,1),fathk)
!fthk=fathk
!!Alex CALL getssh(count(1),count(2),count(3),kapref,pref,thbase,thref,imask,asal,atem,fathk,aden-thbase,fpba,fpsikk,fthkk,fmgp,assh)
CALL getssh(count(1),count(2),count(3),kapref,pref,thbase,thref,imask,asal,atem,fathk,aden-thbase,fpba,fpsikk,fthkk,amgp,assh)
!print *, minval(assh),maxval(assh,mask=assh<2.00**100)
!!Alex test uncomment next line
!amgp=amgp-fmgp

END SUBROUTINE


SUBROUTINE read_analysis_parameters()
IMPLICIT NONE

NAMELIST /anlp/acdepth,arlocal,acovscl
! read analysis parameters
OPEN(21,FILE='tsis.nlist')
READ(21,NML=anlp)
CLOSE(21)
END SUBROUTINE

SUBROUTINE get_analysis_params(var,cdpth,rlocal,covscl)
IMPLICIT NONE
CHARACTER(len=*)         :: var
REAL(i4)                 :: cdpth
REAL(i4)                 :: rlocal
REAL(i4)                 :: covscl


if(var .eq. 'ssh') then
   cdpth=acdepth(1)
   rlocal=arlocal(1)
   covscl=acovscl(1)
!   work2d=>fssh
else if(var .eq. 'sst') then
   cdpth=acdepth(2)
   rlocal=arlocal(2)
   covscl=acovscl(2)
!   work2d=>fsst
else if(var .eq. 'sss') then
   cdpth=acdepth(3)
   rlocal=arlocal(3)
   covscl=acovscl(3)
else if(var .eq. 'tem') then
   cdpth=acdepth(4)
   rlocal=arlocal(4)
   covscl=acovscl(4)
else if(var .eq. 'sal') then
   cdpth=acdepth(5)
   rlocal=arlocal(5)
   covscl=acovscl(5)
else if(var .eq. 'thk') then
   cdpth=acdepth(6)
   rlocal=arlocal(6)
   covscl=acovscl(6)
else if(var .eq. 'pin') then
   cdpth=acdepth(7)
   rlocal=arlocal(7)
   covscl=acovscl(7)
else if(var .eq. 'den') then
   cdpth=acdepth(8)
   rlocal=arlocal(8)
   covscl=acovscl(8)
else if(var .eq. 'uvl') then
   cdpth=acdepth(9)
   rlocal=arlocal(9)
   covscl=acovscl(9)
else if(var .eq. 'vvl') then
   cdpth=acdepth(10)
   rlocal=arlocal(10)
   covscl=acovscl(10)
endif

END SUBROUTINE


END MODULE mrls