[cdftools] 59/228: RD: add the contrib of ssh in cdfhdy and shapiro filter in cdfsmooth

[Alastair McKinstry]

This is an automated email from the git hooks/post-receive script.

mckinstry pushed a commit to branch master
in repository cdftools.

commit 8b1828e83ad3841771f0dd8d4d40198a8577192e
Author: dussin <dussin at 1055176f-818a-41d9-83e1-73fbe5b947c5>
Date:   Wed Jun 23 16:42:27 2010 +0000

    RD: add the contrib of ssh in cdfhdy and shapiro filter in cdfsmooth
    
    
    git-svn-id: http://servforge.legi.grenoble-inp.fr/svn/CDFTOOLS/trunk@335 1055176f-818a-41d9-83e1-73fbe5b947c5
---
 cdfhdy.f90    | 10 ++++---
 cdfsmooth.f90 | 92 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 2 files changed, 96 insertions(+), 6 deletions(-)

diff --git a/cdfhdy.f90 b/cdfhdy.f90
index d52c3c1..15c06dc 100644
--- a/cdfhdy.f90
+++ b/cdfhdy.f90
@@ -100,8 +100,8 @@ PROGRAM cdfhdy
   zsal0(:,:)=35.
 
   zmask(:,:) = getvar(cmask, 'tmask', zlev2, npiglo, npjglo)
+  zssh(:,:)  = getvar(cfilet, 'sossheig', 1, npiglo, npjglo)
   ze3t_1d(:) = getvare3(coordzgr, 'e3t',npk)
-!  PRINT *, 'e3t = ', ze3t_1d
 
   DO jt=1,npt
     PRINT *,' TIME = ', jt, tim(jt)/86400.,' days'
@@ -111,10 +111,14 @@ PROGRAM cdfhdy
 
   DO jk = zlev1, zlev2
 
-!     zdep(:,:)   = getvar(coordzgr, 'e3t_ps', jk,npiglo,npjglo,ldiom=.true.)
+     !zdep(:,:)   = getvar(coordzgr, 'e3t_ps', jk,npiglo,npjglo,ldiom=.true.)
      ! we degrade the computation to smooth the results
      zdep(:,:) = ze3t_1d(jk)
 
+     IF ( zlev1 == 1) THEN
+        zdep(:,:) = zdep(:,:) + zssh(:,:)
+     ENDIF
+
      ! total depth at current level (used for computation of rho in situ)
      zdepth(:,:) = zdepth(:,:) + zdep(:,:)
 
@@ -124,8 +128,6 @@ PROGRAM cdfhdy
      CALL eos_insitu( ztemp0, zsal0, zdepth, npiglo, npjglo, zsig0 ) 
      CALL eos_insitu( ztemp, zsal, zdepth, npiglo, npjglo, zsig ) 
 
-     PRINT *, 'max of ref profile for level', jk ,'is ', MAXVAL(zsig0)
-
      ! we compute the term of the integral : (1/g) *10e4 * sum [ delta * dz ]
      ! with delta = (1/rho - 1/rho0)
      ! 10e4 factor is conversion decibar/pascal
diff --git a/cdfsmooth.f90 b/cdfsmooth.f90
index b30dddc..c37ca96 100644
--- a/cdfsmooth.f90
+++ b/cdfsmooth.f90
@@ -201,7 +201,7 @@ CONTAINS
     CASE ( 2 )
        CALL lishan2d(px,iw,py,ncoup,npiglo,npjglo)
     CASE ( 3 )
-       print *,' not available'
+       CALL lisshapiro1d(px,iw,py,ncoup,npiglo,npjglo)
     CASE ( 4 )
        CALL lisbox(px,kpx,py,npiglo,npjglo,fn,nband,npiglo,npjglo)
     END SELECT
@@ -255,7 +255,9 @@ CONTAINS
   SUBROUTINE initshap(pfn,knj)
     INTEGER, INTENT(in)     :: knj  !: bandwidth
     REAL(KIND=4),INTENT(in) ::  pfn
-    PRINT *,' Init shap not done already' ; STOP
+
+!   nothing to do 
+
   END SUBROUTINE initshap
 
   SUBROUTINE initbox(pfn,knj)
@@ -402,6 +404,92 @@ CONTAINS
 
   END SUBROUTINE lishan2d
 
+  SUBROUTINE lisshapiro1d(px,kiw,py,order,kpi,kpj)
+    !----------------------------------------------
+    !   px  = input data
+    !   kiw = validity of input data
+    !   py  = output filter
+    !   n=number of input/output data
+    !
+    !   adapted from Mercator code...
+    !--------------------------------------------
+    ! * Arguments
+    INTEGER, INTENT(in) :: kpi, kpj, order
+    INTEGER,DIMENSION(:,:),INTENT(in) :: kiw
+    REAL(KIND=4), DIMENSION(:,:), INTENT(in)  :: px
+    REAL(KIND=4), DIMENSION(:,:), INTENT(out) :: py
+    REAL(KIND=4), PARAMETER :: rp_aniso_diff_XY=2.25 !  anisotrope case
+    REAL(KIND=4)            :: zalphax, zalphay, znum
+
+    ! local
+    REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: ztmp , zpx , zpy, zkiw
+    INTEGER :: imin, imax, halo=0
+    INTEGER :: jj, ji, iorder  !: loop indexes
+
+    LOGICAL :: cycl = .true.
+
+    !PRINT *, 'east-west periodicity is assumed ' , cycl
+
+    IF(cycl) halo=1
+    ! we allocate with an halo
+    ALLOCATE( ztmp(0:kpi+halo,kpj) , zpx(0:kpi+halo,kpj) , zpy(0:kpi+halo,kpj) , zkiw(0:kpi+halo,kpj) )
+
+    IF(cycl) THEN
+       zpx(1:kpi,:) = px(:,:)   ;  zkiw(1:kpi,:) = kiw(:,:)
+       zpx(0,:)     = px(kpi,:) ;  zkiw(0,:)     = kiw(kpi,:)
+       zpx(kpi+1,:) = px(1,:)   ;  zkiw(kpi+1,:) = kiw(1,:)
+    ELSE
+       zpx(:,:) = px(:,:)
+    ENDIF
+
+    zpy(:,:)  = zpx(:,:)  ! init?
+    ztmp(:,:) = zpx(:,:)  ! init
+
+    zalphax=1./2.
+    zalphay=1./2.
+
+    !  Dx/Dy=rp_aniso_diff_XY  , D_ = vitesse de diffusion
+    !  140 passes du fitre, Lx/Ly=1.5, le rp_aniso_diff_XY correspondant est:
+
+    IF ( rp_aniso_diff_XY >=  1. ) zalphay=zalphay/rp_aniso_diff_XY
+    IF ( rp_aniso_diff_XY <   1. ) zalphax=zalphax*rp_aniso_diff_XY
+
+    DO iorder=1,order
+
+    imin=2-halo
+    imax=kpi-1+halo
+
+        DO ji = imin,imax
+           DO jj = 2,kpj-1
+              ! We crop on the coast
+               znum = ztmp(ji,jj)   &
+                      + 0.25*zalphax*(ztmp(ji-1,jj  )-ztmp(ji,jj))*zkiw(ji-1,jj)  &
+                      + 0.25*zalphax*(ztmp(ji+1,jj  )-ztmp(ji,jj))*zkiw(ji+1,jj)  &
+                      + 0.25*zalphay*(ztmp(ji  ,jj-1)-ztmp(ji,jj))*zkiw(ji  ,jj-1)  &
+                      + 0.25*zalphay*(ztmp(ji  ,jj+1)-ztmp(ji,jj))*zkiw(ji  ,jj+1)
+               zpy(ji,jj)=znum*zkiw(ji,jj)+zpx(ji,jj)*(1.-zkiw(ji,jj))
+            ENDDO  ! end loop ji
+        ENDDO  ! end loop jj
+
+    IF(cycl) THEN
+       zpy(0,:)     = zpy(kpi,:) 
+       zpy(kpi+1,:) = zpy(1,:) 
+    ENDIF
+
+    ! update the tmp array
+    ztmp(:,:) = zpy(:,:)
+
+    ENDDO
+
+  ! return this array
+  IF(cycl) THEN
+     py(:,:) = zpy(1:kpi,:)
+  ELSE
+     py(:,:) = zpy(:,:)
+  ENDIF
+
+  END SUBROUTINE lisshapiro1d
+
   SUBROUTINE lisbox(px,kiw,py,knx,kny,pfn,knj,kpi,kpj)
     !----------------------------------------------
     ! perform a box car 2d filtering, of span knj

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/cdftools.git