Index: LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/README
===================================================================
--- LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/README	(revision 4171)
+++ LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/README	(revision 4171)
@@ -0,0 +1,41 @@
+make_sso
+
+routine to compute subfgrid scale orography parameters
+from a spherical harmonics method
+
+Requirements:
+- SPHEREPACK library
+- FFTW ?
+- ifortran
+- netcdf
+- an adequate environment defined in make_sso.env
+
+Options:
+
+ ./make_sso  -i <file>  -v <var>  -res <nlon> <nlat>  -m <mask>  -f <ffac>  -t
+  <tfac>
+ 
+ Where:
+ 
+  <file>:   input file name
+  <var> :   height variable name
+  <nlon>:   number of distinct longitude points of output grid (iim)
+  <nlat>:   number of distinct latitude  points of output grid (jjm+1)
+  <mask>:   can be:  1. a mask file (o2a.nc)
+                     2. "noro": computed with grid_noro
+            default: 3. "spec": computed internally
+  <ffac>:   filtering  width: ffac*Do         (Do:output resol.)
+  <tfac>:   transition width: Di+tfac*(Do-Di) (Di: input resol.)
+ 
+ Note that latitudes grid contains both poles.
+
+
+Example to compile then run make_sso from a high-resolution Relief.nc file to
+a 180x90 horizontal grid. 
+
+./compile
+file=Relief.nc
+./make_sso -i $file -res 180 91 -v z
+
+Output file is make_sso_180x91_${file}
+
Index: LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/compile
===================================================================
--- LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/compile	(revision 4171)
+++ LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/compile	(revision 4171)
@@ -0,0 +1,39 @@
+source make_sso.env
+
+
+if [ ! -f Relief.nc ]; then
+echo Attention missing Relief.nc file
+echo I download it
+wget --no-check-certificate -nv https://web.lmd.jussieu.fr/~lmdz/pub/3DInputData/Limit/Relief.nc
+fi
+
+if [ ! -d SPHEREPACK ]; then
+echo Missing SPHEREPACK library. I download it
+wget --no-check-certificate -nv https://web.lmd.jussieu.fr/~lmdz/pub/import/spherepack3.2.tgz 
+tar -xvf spherepack3.2.tgz
+mv spherepack3.2 SPHEREPACK
+mv SPHEREPACK/make.inc SPHEREPACK/make.inc.old
+
+# copy more adequate make.inc to SPHEREPACK
+mv make.inc SPHEREPACK/.
+# copy corrected Makefile to SPHEREPACK/test
+mv SPHEREPACK/test/Makefile SPHEREPACK/test/Makefile.old
+mv  Makefile.test SPHEREPACK/test/Makefile 
+# clean up
+rm spherepack3.2.tgz
+fi
+
+
+
+if [ ! -f SPHEREPACK/lib/libspherepack.a ]; then
+  cd SPHEREPACK ; gmake ; cd -
+fi
+SRC="make_sso_SpherePack.f90"
+mode=prod
+
+case $mode in
+  prod) OPTS="-O3" ;;
+  debug) OPTS="-check all -g" ;;
+esac
+
+ifort -assume realloc_lhs $OPTS -o make_sso.e make_sso_SpherePack.f90 -I$NETCDFFORTRAN_INCDIR $netcdf -L$NETCDFFORTRAN_LIBDIR -lnetcdff -lnetcdf -L./SPHEREPACK/lib -lspherepack -L$FFTW3_LIBDIR -lfftw3
Index: LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make.inc
===================================================================
--- LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make.inc	(revision 4171)
+++ LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make.inc	(revision 4171)
@@ -0,0 +1,11 @@
+
+LIB=../lib/libspherepack.a
+
+F90 := ifort -I../lib 
+CPP := ifort -E 
+
+
+MAKE := gmake
+AR := /usr/bin/ar
+
+
Index: LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso
===================================================================
--- LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso	(revision 4171)
+++ LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso	(revision 4171)
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+source make_sso.env
+./make_sso.e $*
Index: LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso.env
===================================================================
--- LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso.env	(revision 4171)
+++ LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso.env	(revision 4171)
@@ -0,0 +1,20 @@
+case ${HOSTNAME:0:5} in
+  irene)
+    module unload netcdf-fortran
+    module unload netcdf-c
+    module load netcdf-c
+    module load netcdf-fortran
+    module load fftw3/gnu
+    netcdf="-L$NETCDF_LIBDIR"
+  ;;
+  cicla)
+    module purge
+    module load intel
+    module load netcdf4/4.3.3.1-ifort
+    NETCDF=$(which ncdump); NETCDF=${NETCDF%%/bin*}
+    NETCDFFORTRAN_INCDIR=$NETCDF/include
+    NETCDFFORTRAN_LIBDIR=$NETCDF/lib
+    netcdf=""
+  ;;
+esac
+
Index: LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso_SpherePack.f90
===================================================================
--- LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso_SpherePack.f90	(revision 4171)
+++ LMDZ6/branches/LMDZ-ECRAD/tools/make_sso/make_sso_SpherePack.f90	(revision 4171)
@@ -0,0 +1,1093 @@
+!-------------------------------------------------------------------------------
+!
+PROGRAM make_sso
+!
+!-------------------------------------------------------------------------------
+! Purpose: Project ETOPO file (GMT4 axes conventions) on spherical harmonics.
+!-------------------------------------------------------------------------------
+  USE netcdf
+!  USE sphpack
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Local variables:
+  REAL,    PARAMETER :: xpi=4.*ATAN(1.),         &    !--- pi
+                        hp=xpi/2., d2r=xpi/180., &    !--- pi, pi/2, pi/180  
+                        REarth=2.E7/xpi,         &    !--- Earth radius
+                        eps=1.E-6                     !--- Tolerance for angles
+  !=== Transition length for spectral filters (to avoid side lobes)
+!  REAL,    PARAMETER :: del=2000./REarth              !--- On unit sphere
+  !=== .TRUE./.FALSE.: Fast/slow routine with heavy O(N^3)/light O(N^2) storage
+  LOGICAL, PARAMETER :: lfast=.TRUE.                  !--- Fast/slow routine
+
+!--- From the input arguments
+  CHARACTER(LEN=256), ALLOCATABLE :: args(:)          !--- All the arguments
+  CHARACTER(LEN=2048):: call_seq                      !--- Calling sequence
+  INTEGER            :: nargs                         !--- Input args number
+  CHARACTER(LEN=256) :: f_in                          !--- Input file name
+  CHARACTER(LEN=256) :: vnam                          !--- Height variable name
+  CHARACTER(LEN=256) :: res_in                        !--- Input resolution
+  CHARACTER(LEN=256) :: fmsk                          !--- Mask file or keyword
+  !=== Filter on output field, of resolution=ffactor*output resolution
+  LOGICAL            :: lfilter=.TRUE.                !--- Filter activation
+  REAL               :: ffactor                       !--- Filter resol. factor
+  REAL               :: tfactor                       !--- Transition factor
+  INTEGER            :: nlon_ou, nlat_ou              !--- Output grid resol.
+
+!--- Quantities related to input field grid
+  CHARACTER(LEN=256) :: lonn, latn                    !--- Axes vars names
+  CHARACTER(LEN=256) :: lonu, latu                    !--- Angular units
+  REAL, ALLOCATABLE  :: lon_in(:), lat_in(:)          !--- Cells centers
+  INTEGER            :: nlon_in,   nlat_in            !--- Input  grid dims
+  INTEGER            :: igrid_in(2)                   !--- Grid identification
+  REAL               :: dlon, dlat                    !--- Domain sizes
+  REAL               :: D1, del1                      !--- Resol. at equator
+  REAL               :: al                            !--- Longitude ratio
+  LOGICAL            :: lcyc                          !--- T: redundant longitude
+  LOGICAL            :: lXdec, lYinc                  !--- T: dec/increasing lon/lat
+  REAL, ALLOCATABLE, DIMENSION(:,:) ::v,h , a , b,mski!--- Heights, coeffs, mask
+  REAL, ALLOCATABLE, DIMENSION(:,:,:) ::hs, as, bs, & !--- Same, filtered
+    dxZ, dyZ, dxZ2, dyZ2, dxyZ, Xk, Xl, Xm, Xp, Xq, dH!--- Grads, correlations
+
+!--- Quantities related to output grid
+  CHARACTER(LEN=256) :: f_ou                          !--- Output file name
+  CHARACTER(LEN=256) :: res_ou                        !--- Output resolution
+  INTEGER(KIND=1), ALLOCATABLE :: mask(:,:)           !--- Boolean out land mask
+  REAL, ALLOCATABLE  :: lon_ou(:), lat_ou(:), tmp1(:) !--- Output grid coords
+  INTEGER            :: igrid_ou(2)                   !--- Grid identification
+  REAL               :: D2, del2                      !--- Resol. at equator
+  REAL, ALLOCATABLE, DIMENSION(:,:) ::   msko, &      !--- Fractional land mask
+                                           h0, &      !--- Large scale orography
+                   Zphi, Zsig, Zthe, Zgam, mu, tmp2   !--- Geopo, SSO parameters
+
+!--- Miscellanous Netcdf-related variables
+  INTEGER :: fID, xID, loID, phiID, sigID, theID, picID, dIDs(4), nn
+  INTEGER :: vID, yID, laID, meaID, gamID, mskID, valID, muID, ix180
+  LOGICAL :: ll
+
+!--- Miscellanous SPHEREPACK-specific quantities
+  INTEGER :: l1, idh, mdab, lsh, la1, la2, lwa, lsvmin, lsave, ierr, intl
+  INTEGER :: l2, jdh, ndab, lvh, lb1, lb2, lwb, lwkmin, lwork, ldwork
+  REAL, ALLOCATABLE :: wsha(:), wsave(:), wvh (:)
+  REAL, ALLOCATABLE :: wshs(:), dwork(:), work(:)
+  INTEGER :: k, n, nt, isym
+
+  CHARACTER(LEN=256) :: sub, fnam, msg, arg
+
+!-------------------------------------------------------------------------------
+  sub='make_sso'
+
+!=== DEAL WITH INPUT ARGUMENTS =================================================
+  call_seq='./'//TRIM(sub)//'.e'
+  nargs=COMMAND_ARGUMENT_COUNT(); ALLOCATE(args(nargs))
+!--- MANUAL DISPLAYED IF NO ARGUMENTS ARE PROVIDED
+  IF(nargs==0) THEN
+    WRITE(*,*)''
+    WRITE(*,*)'General calling sequence:'
+    WRITE(*,*)''
+    WRITE(*,*)TRIM(call_seq)//'  -i <file>  -v <var>  -res <nlon> <nlat>'      &
+      //'  -m <mask>  -f <ffac>  -t <tfac>'
+    WRITE(*,*)''
+    WRITE(*,*)'Where:'
+    WRITE(*,*)''
+    WRITE(*,*)' <file>:   input file name'
+    WRITE(*,*)' <var> :   height variable name'
+    WRITE(*,*)' <nlon>:   number of distinct longitude points of output grid'
+    WRITE(*,*)' <nlat>:   number of distinct latitude  points of output grid'
+    WRITE(*,*)' <mask>:   can be:  1. a mask file (o2a.nc)'
+    WRITE(*,*)'                    2. "noro": computed with grid_noro'
+    WRITE(*,*)'           default: 3. "spec": computed internally'
+    WRITE(*,*)' <ffac>:   filtering  width: ffac*Do         (Do:output resol.)'
+    WRITE(*,*)' <tfac>:   transition width: Di+tfac*(Do-Di) (Di: input resol.)'
+    WRITE(*,*)
+    WRITE(*,*)'Note that latitudes grid contains both poles.'
+    WRITE(*,*)
+    STOP
+  END IF
+!--- ARGUMENTS PARSING
+  DO k=1,nargs
+    CALL GET_COMMAND_ARGUMENT(NUMBER=k,VALUE=args(k))
+    call_seq=TRIM(call_seq)//' '//TRIM(args(k))
+  END DO
+  k=1
+  f_in=''; nlon_ou=0; nlat_ou=0; vnam=''; fmsk=''; ffactor=2.0; tfactor=0.
+  DO
+    k=k+1; ll=.FALSE.; arg=args(k-1)
+    SELECT CASE(arg)
+      CASE('-i');   f_in=args(k); msg='Missing file "'//TRIM(f_in)//'".'
+                    ll=NF90_OPEN(f_in,NF90_NOWRITE,fID)/=NF90_NOERR
+                    IF(.NOT.ll) n=NF90_CLOSE(fID)
+      CASE('-res'); nlon_ou=str2int(args(k)); k=k+1
+                    nlat_ou=str2int(args(k))
+      CASE('-v');   vnam=args(k)
+      CASE('-m');   fmsk=args(k)
+      CASE('-f');   ffactor=str2real(args(k))
+      CASE('-t');   tfactor=str2real(args(k))
+      CASE DEFAULT; msg='Unrecognized argument "'//TRIM(arg)//'".'; ll=.TRUE.
+    END SELECT
+    IF(ll) THEN; WRITE(*,*)msg; STOP; END IF
+    k=k+1
+    IF(k>nargs) EXIT
+  END DO
+!--- CHECK OUTPUT GRID IS DEFINED
+  IF(TRIM(fmsk)=='') fmsk='spec'
+  IF(ALL(['noro','spec']/=fmsk)) THEN
+    msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file'
+    CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=NF90_NOERR,msg)
+    CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean")    !--- MASK ID
+    CALL nc(NF90_INQUIRE_VARIABLE(fID,vID,dimids=dIDs))         !--- DIMS IDS
+    CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(1),len=nlon_ou),'x')!--- NB LONG
+    CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(2),len=nlat_ou),'y')!--- NB LAT
+    CALL nc(NF90_CLOSE(fID))
+  END IF
+  IF(nlon_ou<=0.OR.nlat_ou<=0) THEN
+    WRITE(*,*)'Missing or wrong "-res" arguments for output grid'; STOP
+  END IF
+
+!=== READ THE INPUT FIELD =======================================================
+  CALL nc(NF90_OPEN(f_in,NF90_NOWRITE,fID))
+  WRITE(*,*)'>> Reading variable "'//TRIM(vnam)//'" from "'//TRIM(f_in)//'"...'
+
+!--- GET HEIGHT VARIABLE AND ITS DIMENSIONS 
+  CALL nc(NF90_INQ_VARID(fID,vnam,vID))
+  CALL nc(NF90_INQUIRE_VARIABLE (fID,vID,dimids=dIDs))
+  CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(1),len=nlon_in,name=lonn))
+  CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(2),len=nlat_in,name=latn))
+  WRITE(*,*)TRIM(vnam)//' is '//TRIM(lonn)//'('//TRIM(int2str(nlon_in))//')*'   &
+                              //TRIM(latn)//'('//TRIM(int2str(nlat_in))//')'
+  nlon_in=nlon_in-1                                   !--- DISTINCT POINTS NUMBER
+  ALLOCATE(lon_in(nlon_in+1),lat_in(nlat_in),h(nlon_in+1,nlat_in))
+  CALL nc(NF90_INQ_VARID(fID,lonn,loID)        ,lonn)
+  CALL nc(NF90_INQ_VARID(fID,latn,laID)        ,latn)
+  CALL nc(NF90_GET_VAR  (fID,loID,lon_in)      ,lonn)
+  CALL nc(NF90_GET_VAR  (fID,laID,lat_in)      ,latn)
+  CALL nc(NF90_GET_ATT  (fID,loID,'units',lonu),lonn)
+  CALL nc(NF90_GET_ATT  (fID,laID,'units',latu),latn)
+  CALL nc(NF90_GET_VAR  (fID,vID,h(:,:),[1,1],[nlon_in+1,nlat_in]),vnam)
+  CALL nc(NF90_CLOSE(fID))
+
+!--- CHECK WETHER GRID IS CORRECT (GLOBAL DOMAIN, IDENTIFIED UNITS...)
+  dlon=ABS(lon_in(nlon_in+1)-lon_in(1)); dlat=ABS(lat_in(nlat_in)-lat_in(1))
+  CALL err(lonu(1:6)/='degree','Invalid longitudes unit: "'//TRIM(lonu)//'"')
+  CALL err(latu(1:6)/='degree','Invalid latitudes unit: "' //TRIM(latu)//'"')
+  CALL err(ABS(dlon-360.)>eps,'Longitudes domain is not global.')
+  CALL err(ABS(dlat-180.)>eps, 'Latitudes domain is not global.')
+
+!--- ENSURE LONGITUDES ARE INCREASING AND LATITUDES ARE DECREASING
+  lXdec=lon_in(1)>lon_in(nlon_in)
+  lYinc=lat_in(1)<lat_in(nlat_in)
+  IF(lXdec) THEN; lon_in=lon_in(nlon_in+1:1:-1); h(:,:)=h(nlon_in+1:1:-1,:); END IF
+  IF(lYinc) THEN; lat_in=lat_in(nlat_in  :1:-1); h(:,:)=h(:,nlat_in  :1:-1); END IF
+
+!--- ENSURE LONGITUDES ARE BETWEEN -180 and 180deg WITHOUT REDUNDANT END POINT
+  DO ix180=1,nlon_in+1; IF(lon_in(ix180)>=180.) EXIT; END DO
+  IF(ix180>=nlon_in+1) THEN
+    lon_in(:)=[lon_in(ix180:nlon_in)-360.,lon_in(1:ix180)]
+    ALLOCATE(tmp2(nlon_in+1,nlat_in)); nn=nlon_in-ix180+1
+    tmp2(1:nn,:)=h(ix180:nlon_in,:); tmp2(nn+1:nlon_in+1,:)=h(1:ix180,:)
+    DEALLOCATE(tmp2)
+  END IF
+
+!--- REMOVE OCEAN ; FLIP FIELD (SPHEREPACK: LAT FIRST) ; INPUT BOOLEAN MASK
+  WHERE(h(:,:)<0.) h(:,:)=0
+  idh=nlat_in; jdh=nlon_in; CALL flip(h)
+  ALLOCATE(mski(idh,jdh)); mski(:,:)=0.; WHERE(h(:,:)>0.) mski(:,:)=1.
+
+!=== SPHERICAL HARMONICS ANALYSIS ===============================================
+!--- Allocate work vector (min. length depends on used routines)
+  nt=2                                                        ! MAX nt
+  l1=MIN(nlat_in,(nlon_in+2)/2)
+  IF(MODULO(nlon_in,2)==1) l1=MIN(nlat_in,(nlon_in+1)/2)
+  l2=(nlat_in+1)/2
+  IF(lfast) THEN
+    lwork=2*(nlat_in+1)                                       ! shaeci/shseci
+    lwork=MAX(lwork,nlat_in*(nt*nlon_in+MAX(3*l2,nlon_in)))   ! shaec /shaes
+    lwork=MAX(lwork,4*(nlat_in+1))                            ! vhseci
+    lwork=MAX(lwork,nlat_in*(2*nt*nlon_in+MAX(6*l2,nlon_in))+nlat_in*(2*l1*nt+1)) ! gradec
+    lsh=2*nlat_in*l2+3*(  (l1-2  )*(nlat_in+nlat_in-l1-1))/2+nlon_in+15 ! shaec /shsec
+    lvh=4*nlat_in*l2+3*MAX(l1-2,0)*(nlat_in+nlat_in-l1-1)   +nlon_in+15 ! vhseci/gradec
+  ELSE
+    lwork=5*nlat_in*l2+3*((l1-2)*(nlat_in+nlat_in-l1-1))/2    ! shaesi/shsesi
+    lwork=MAX(lwork,(nt+1)*nlat_in*nlon_in)                   ! shaes /shses
+    lwork=MAX(lwork,3*(MAX(l1-2,0)*(nlat_in+nlat_in-l1-1))/2+5*l2*nlat_in)     ! vhsesi
+    lwork=MAX(lwork,nlat_in*((2*nt+1)*nlon_in+2*l1*nt+1))     ! grades
+    lsh=(l1*l2*(nlat_in+nlat_in-l1+1))/2+nlon_in+15           ! shaes /shses
+    lvh= l1*l2*(nlat_in+nlat_in-l1+1)   +nlon_in+15           ! vhsesi/grades
+  END IF
+  ALLOCATE(work(lwork))
+  ALLOCATE(wsha(lsh),wshs(lsh),wvh(lvh))
+
+!--- Allocate spectral coefficients vectors
+  mdab=l1; ndab=nlat_in; IF(.NOT.lfast) mdab=MIN(nlat_in,(nlon_in+2)/2)
+  ALLOCATE(a(mdab,ndab),b(mdab,ndab))
+  isym=0; nt=1                                                !--- Symetries, fields nb
+
+!--- FIELD PROJECTION ON THE SPHERICAL HARMONICS BASE.
+  WRITE(*,*)'>> Projecting fields on the spherical harmonics...'
+  IF(lfast) THEN
+    CALL shaeci(nlat_in,nlon_in,wsha,lsh,work,lwork,ierr)
+    CALL error('shaeci',ierr)
+    CALL shaec(nlat_in,nlon_in,isym,nt,h,idh,jdh,a,b,mdab,ndab,wsha,lsh,work,lwork,ierr)
+    CALL error('shaec',ierr)
+  ELSE
+    CALL shaesi(nlat_in,nlon_in,wsha,lsh,work,lwork,ierr)
+    CALL error('shaesi',ierr)
+    CALL shaes(nlat_in,nlon_in,isym,nt,h,idh,jdh,a,b,mdab,ndab,wsha,lsh,work,lwork,ierr)
+    CALL error('shaes',ierr)
+  END IF
+
+!--- ALIASING DAMPLING FILTER
+!  WRITE(*,*)'>> Bi-Laplacian filtering to avoid Gibbs phenomenon...'
+!  CALL diffusion_filter(a,b)
+
+!=== INPUT/OUTPUT RESOLUTIONS AT EQUATOR (normalized with Earth radius)
+!--- The area of a cell at (lam,phi) of size(dl,dp) is: dl.SIN(dp).COS(phi)
+!--- At equator, Res=SQRT(dl.SIN(dp)) ; with dl=2.pi/nlon and dp=pi/nlat:
+!---             Res=SQRT(2.pi.SIN(pi/nlat)/nlon)~pi*SQRT(2/(nlon*nlat))
+!--- Note: nlat is the slices number = latitude points number -1
+  D1=SQRT(2.*xpi*SIN(xpi/(nlat_in-1))/nlon_in)   !--- Input  grid
+  D2=SQRT(2.*xpi*SIN(xpi/(nlat_ou-1))/nlon_ou)   !--- Output grid
+
+!=== FILTER FACTOR EQUAL TO ZERO => NO FILTERING
+  lfilter=ffactor/=0.
+
+!=== DEFAULT TRANSITION SCALE: INPUT GRID RESOLUTION
+  WRITE(*,*)'tfactor=',tfactor
+  del1=D1; del2=D1+(D2-D1)*tfactor
+  WRITE(*,*)'2*D1,del1,D2,del2=',2*D1,del1,D2,del2
+
+!=== SCALES SEPARATION AT THE OUTPUT GRID RESOLUTION ============================
+  nt=2
+  WRITE(*,*)'>> Scales separation at '//TRIM(real2str(D2/d2r))//' degrees (deta'&
+      //'ils below '//TRIM(real2str(60*D1/d2r))//' minutes removed)...'
+  ALLOCATE(as(mdab,ndab,nt),bs(mdab,ndab,nt))
+  DO n=1,nt; as(:,:,n)=a; bs(:,:,n)=b; END DO
+  CALL  lowpass_filter(as(:,:,1),bs(:,:,1),          D2,del2)   !--- LARGE SCALES
+  CALL bandpass_filter(as(:,:,2),bs(:,:,2),2*D1,del1,D2,del2)   !--- SMALL SCALES
+
+!--- BACKWARD TRANSFORM FOR LARGE AND SUB-CELL SCALES OROGRAPHYS
+  WRITE(*,*)'>> Backward transform of large and sub-cell scales orographies...'
+  ALLOCATE(hs(idh,jdh,nt))
+  IF(lfast) THEN
+    CALL shseci(nlat_in,nlon_in,wshs,lsh,work,lwork,ierr)
+    CALL error('shseci',ierr)
+    CALL shsec(nlat_in,nlon_in,isym,nt,hs,idh,jdh,as,bs,mdab,ndab,wshs,lsh,work,lwork,ierr)
+    CALL error('shsec',ierr)
+  ELSE
+    CALL shsesi(nlat_in,nlon_in,wshs,lsh,work,lwork,ierr)
+    CALL error('shsesi',ierr)
+    CALL shses(nlat_in,nlon_in,isym,nt,hs,idh,jdh,as,bs,mdab,ndab,wshs,lsh,work,lwork,ierr)
+    CALL error('shses',ierr)
+  END IF
+  WHERE(mski==0.); hs(:,:,1)=0.; hs(:,:,2)=0.; END WHERE
+  hs(:,:,2)=hs(:,:,2)**2
+
+!=== COMPUTE GRADIENT ; n=1: large scale ; n=2: sub-cell grid ===================
+  WRITE(*,*)'>> Computing large and sub-cell scales gradients...'
+  ALLOCATE(dxZ(idh,jdh,nt),dyZ(idh,jdh,nt))
+  IF(lfast) THEN
+    CALL vhseci(nlat_in,nlon_in,wvh,lvh,work,lwork,ierr)
+    CALL error('vhseci',ierr)
+    DO n=1,2
+      CALL gradec(nlat_in,nlon_in,isym,1,dyZ(:,:,n),dxZ(:,:,n),idh,jdh,         &
+                           as(:,:,n),bs(:,:,n),mdab,ndab,wvh,lvh,work,lwork,ierr)
+    END DO
+    CALL error('gradec',ierr)
+  ELSE
+    CALL vhsesi(nlat_in,nlon_in,wvh,lvh,work,lwork,ierr)
+    CALL error('vhsesi',ierr)
+    DO n=1,2
+      CALL gradec(nlat_in,nlon_in,isym,1,dyZ(:,:,n),dxZ(:,:,n),idh,jdh,as(:,:,n),&
+                                   bs(:,:,n),mdab,ndab,wvh,lvh,work,lwork,ierr)
+    END DO
+    CALL error('grades',ierr)
+  END IF
+  DEALLOCATE(as,bs)
+  dxZ=dxZ/REarth; dyZ=dyZ/REarth
+  IF(lat_in(nlat_in)>lat_in(1)) dyZ=-dyZ
+  WHERE(mski==0.)
+    dxZ(:,:,1)=0.; dyZ(:,:,1)=0.
+    dxZ(:,:,2)=0.; dyZ(:,:,2)=0.
+  END WHERE
+  ALLOCATE(dH(idh,jdh,6))
+  dH(:,:,1)=dxZ(:,:,1)*dxZ(:,:,1); dH(:,:,4)=dxZ(:,:,2)*dxZ(:,:,2)
+  dH(:,:,2)=dyZ(:,:,1)*dyZ(:,:,1); dH(:,:,5)=dyZ(:,:,2)*dyZ(:,:,2)
+  dH(:,:,3)=dxZ(:,:,1)*dyZ(:,:,1); dH(:,:,6)=dxZ(:,:,2)*dyZ(:,:,2)
+  DEALLOCATE(dxZ,dyZ)
+
+!===============================================================================
+  IF(lfilter) THEN    !=== FILTER TO REMOVE ALIASING ===========================
+!===============================================================================
+    WRITE(*,*)'>> Smoothing large and sub-cell scales fields at '//            &
+             &TRIM(real2str(ffactor*D2/d2r))//' degrees...'
+    CALL spatial_filter_m(dH       ,ffactor*D2,del2)
+    CALL spatial_filter_1(hs(:,:,2),ffactor*D2,del2)
+  END IF
+  DEALLOCATE(work)
+!===============================================================================
+
+!=== DEAL WITH OUTPUT GRID =====================================================
+  intl=0
+!--- OUTPUT GRID (CONTAINS A LONGITUDE REDUNDANT POINT FOR NOW)
+  lon_ou = [(360.*DBLE(n)/DBLE(nlon_ou  ), n = 0, nlon_ou)] + lon_in(1)
+  lat_ou = [(180.*DBLE(n)/DBLE(nlat_ou-1), n = nlat_ou-1, 0, -1)] - 90.
+  la1 = MIN(nlat_in,(nlon_in+1)/2); la2 = (nlat_in+1)/2
+  lb1 = MIN(nlat_ou,(nlon_ou+1)/2); lb2 = (nlat_ou+1)/2
+  lwa = 4*nlat_in*la2+3*MAX(la1-2,0)*(2*nlat_in-la1-1)+la2+nlon_in+15
+  lwb = 4*nlat_ou*lb2+3*MAX(lb1-2,0)*(2*nlat_ou-lb1-1)+    nlon_ou+15
+  lsave=lwa+lwb
+  l1 = MAX(nlat_in,nlat_ou)
+  l2 = MAX(nlon_in,nlon_ou)
+  lwork  = 2*l1*(8*MIN(l1,(l2+2)/2)+4*l2+3)
+  ldwork = 2*l1*(l1+1)+1
+  ALLOCATE(wsave(lsave),work(lwork),dwork(ldwork))
+
+!--- Latitude grid: equally spaced partition of [-pi/2,pi/2]
+  igrid_in(:)=[-1,1]                                  !--- NP->SP, nlat*nlon
+  igrid_ou(:)=[-1,0]                                  !--- NP->SP, nlon*nlat
+
+!--- BUILD FRACTIONAL LAND MASK msko AND GEOPOTENTIAL ('noro' CASE)
+  SELECT CASE(fmsk)
+    CASE('noro')
+      WRITE(*,*)'>> Computing mask & geop (grid_noro compatible)...'
+      ALLOCATE(msko(nlon_ou+1,nlat_ou), Zphi(nlon_ou+1,nlat_ou))
+      CALL flip(h)
+      CALL grid_noro0(lon_in*d2r, lat_in*d2r, h,lon_ou*d2r, lat_ou*d2r, Zphi, msko)
+      msko = msko(1:nlon_ou,:)
+      Zphi = Zphi(1:nlon_ou,:)
+    CASE('spec')
+      WRITE(*,*)'>> Computing mask in spectral space...'
+      ALLOCATE(msko(nlon_ou,nlat_ou))
+      CALL trssph(intl, igrid_in, nlon_in, nlat_in, mski,                      &
+                        igrid_ou, nlon_ou, nlat_ou, msko, wsave, lsave, lsvmin,&
+                                     work, lwork, lwkmin, dwork, ldwork, ierr)
+      intl=1
+    CASE DEFAULT
+      WRITE(*,*)'>> Trying to read mask from "'//TRIM(fmsk)//'" file'
+      ALLOCATE(msko(nlon_ou,nlat_ou))
+      msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file'
+      CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=NF90_NOERR,msg)
+      CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean")    !--- MASK ID
+      CALL nc(NF90_GET_VAR(fID,vID,msko(:,:)))                    !--- MASK
+      CALL nc(NF90_CLOSE(fID))
+      msko(:,:)=1.0-msko(:,:)
+  END SELECT
+  DEALLOCATE(mski,h)
+
+!--- GET RID OF REDUNDANT LONGITUDE (USELESS FOR SPHEREPACK)
+  lon_ou=lon_ou(1:nlon_ou)
+
+!=== REGRIDDING ON THE TARGET OUTPUT GRID ======================================
+  WRITE(*,*)'>> Regridding on the target grid...'
+
+!--- COMPUTE BOOLEAN OUTPUT MASK
+  ALLOCATE(mask(nlon_ou,nlat_ou)); mask=0_1; WHERE(msko>=0.5) mask=1_1
+
+!--- MEAN OROGRAPHY
+  ALLOCATE(h0(nlon_ou,nlat_ou))
+  CALL trssph(intl, igrid_in, nlon_in, nlat_in, hs(:,:,1),                     &
+                    igrid_ou, nlon_ou, nlat_ou, h0,                            &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  WHERE(h0<0..OR.mask==0_1) h0(:,:)=0.
+
+!--- STANDARD DEVIATION OF SUB-CELL SCALES OROGRAPHY
+  ALLOCATE(mu(nlon_ou,nlat_ou))
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, hs(:,:,2),                        &
+                 igrid_ou, nlon_ou, nlat_ou, mu,                               &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  DEALLOCATE(hs)
+  WHERE(mu(:,:)<0..OR.mask==0_1) mu(:,:)=0.; mu=SQRT(mu)
+
+!--- SQUARE OF LARGE AND SUB-CELL SCALES ZONAL SLOPE
+  ALLOCATE(dxZ2(nlon_ou,nlat_ou,2))
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, dH  (:,:,1),                      &
+                 igrid_ou, nlon_ou, nlat_ou, dxZ2(:,:,1),                      &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, dH  (:,:,4),                      &
+                 igrid_ou, nlon_ou, nlat_ou, dxZ2(:,:,2),                      &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  WHERE(mask==0_1); dxZ2(:,:,1)=0.; dxZ2(:,:,2)=0.; END WHERE
+
+!--- SQUARE OF LARGE AND SUB-CELL SCALES MERIDIONAL SLOPE
+  ALLOCATE(dyZ2(nlon_ou,nlat_ou,2))
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, dH  (:,:,2),                      &
+                 igrid_ou, nlon_ou, nlat_ou, dyZ2(:,:,1),                      &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, dH  (:,:,5),                      &
+                 igrid_ou, nlon_ou, nlat_ou, dyZ2(:,:,2),                      &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  WHERE(mask==0_1); dyZ2(:,:,1)=0.; dyZ2(:,:,2)=0.; END WHERE
+
+!--- PRODUCT OF LARGE AND SUB-CELL SCALES MERIDIONAL AND ZONAL SLOPES
+  ALLOCATE(dxyZ(nlon_ou,nlat_ou,2))
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, dH  (:,:,3),                      &
+                 igrid_ou, nlon_ou, nlat_ou, dxyZ(:,:,1),                      &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  CALL trssph(1, igrid_in, nlon_in, nlat_in, dH  (:,:,6),                      &
+                 igrid_ou, nlon_ou, nlat_ou, dxyZ(:,:,2),                      &
+    wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr)
+  WHERE(mask==0_1); dxyZ(:,:,1)=0.; dxyZ(:,:,2)=0.; END WHERE
+  DEALLOCATE(dH)
+
+!--- COMPUTE LARGE AND SUB-CELL SCALES DERIVED QUANTITIES
+  ALLOCATE(Xk(nlon_ou,nlat_ou,2)); Xk=(dxZ2+dyZ2)/2.
+  ALLOCATE(Xl(nlon_ou,nlat_ou,2)); Xl=(dxZ2-dyZ2)/2.
+  DEALLOCATE(dxZ2,dyZ2)
+  ALLOCATE(Xm(nlon_ou,nlat_ou,2)); Xm=dxyZ
+  DEALLOCATE(dxyZ)
+  ALLOCATE(Xp(nlon_ou,nlat_ou,2),Xq(nlon_ou,nlat_ou,2))
+  Xp=Xk; Xq=Xk; Xk=SQRT(Xl**2+Xm**2); Xp=Xp-Xk; Xq=Xq+Xk
+  DEALLOCATE(Xk)
+  WHERE(ABS(Xm)<=1e-16) Xm=1e-16*SIGN(1.,Xm)
+  ALLOCATE(Zthe(nlon_ou,nlat_ou)); Zthe=0.5*ATAN2(Xm(:,:,1),Xl(:,:,1))/d2r
+  DEALLOCATE(Xm,Xl)
+  ALLOCATE(Zgam(nlon_ou,nlat_ou)); Zgam=0.
+  WHERE(Xq(:,:,1)>0..AND.Xp(:,:,1)>1.E-8) Zgam(:,:)=Xp(:,:,1)/Xq(:,:,1)
+  DEALLOCATE(Xp)
+  ALLOCATE(Zsig(nlon_ou,nlat_ou)); Zsig=SQRT(Xq(:,:,2))
+  WHERE(Xq(:,:,2)<=1.E-16) Zsig(:,:)=1.E-8
+  DEALLOCATE(Xq)
+  WHERE(mask==0_1); Zthe=0.; Zgam=0.; Zsig=0.; END WHERE
+
+!--- Save the data
+  WRITE(*,*)'>> Saving fields in a netcdf file...'
+  res_in=TRIM(int2str(nlon_in))//'x'//TRIM(int2str(nlat_in))
+  res_ou=TRIM(int2str(nlon_ou))//'x'//TRIM(int2str(nlat_ou))
+  f_ou='make_sso_'//TRIM(res_ou)//'_'//TRIM(f_in)
+  fnam=f_ou
+  CALL nc(NF90_CREATE(f_ou,NF90_CLOBBER,fID))
+
+  CALL nc(NF90_DEF_DIM(fID,'x',nlon_ou,xID))
+  CALL nc(NF90_DEF_VAR(fID,'x',NF90_REAL,xID,loID)      ,'x')
+  CALL nc(NF90_PUT_ATT(fID,loID,'long_name','Longitude'),'x')
+  CALL nc(NF90_PUT_ATT(fID,loID,'units','degrees_east') ,'x')
+
+  CALL nc(NF90_DEF_DIM(fID,'y',nlat_ou,yID))
+  CALL nc(NF90_DEF_VAR(fID,'y',NF90_REAL,yID,laID)      ,'y')
+  CALL nc(NF90_PUT_ATT(fID,laID,'long_name','Latitude') ,'y')
+  CALL nc(NF90_PUT_ATT(fID,laID,'units','degrees_north'),'y')
+
+  CALL nc(NF90_DEF_VAR(fID,'mask',NF90_REAL,[xID,yID],mskID),'mask')
+  IF(fmsk=='noro') &
+  CALL nc(NF90_DEF_VAR(fID,'Zphi',NF90_REAL,[xID,yID],phiID),'Zphi')
+  CALL nc(NF90_DEF_VAR(fID,'Zmea',NF90_REAL,[xID,yID],meaID),'Zmea')
+  CALL nc(NF90_DEF_VAR(fID,'mu'  ,NF90_REAL,[xID,yID], muID),'mu'  )
+  CALL nc(NF90_DEF_VAR(fID,'Zsig',NF90_REAL,[xID,yID],sigID),'Zsig')
+  CALL nc(NF90_DEF_VAR(fID,'Zgam',NF90_REAL,[xID,yID],gamID),'Zgam')
+  CALL nc(NF90_DEF_VAR(fID,'Zthe',NF90_REAL,[xID,yID],theID),'Zthe')
+  CALL nc(NF90_DEF_VAR(fID,'Zpic',NF90_REAL,[xID,yID],picID),'Zpic')
+  CALL nc(NF90_DEF_VAR(fID,'Zval',NF90_REAL,[xID,yID],valID),'Zval')
+
+  CALL nc(NF90_PUT_ATT(fID,mskID,'long_name','Fractional land mask'                      ),'mask')
+  IF(fmsk=='noro') &
+  CALL nc(NF90_PUT_ATT(fID,phiID,'long_name','Geopotential'                              ),'Zphi')
+  CALL nc(NF90_PUT_ATT(fID,meaID,'long_name','Mean orography'                            ),'Zmea')
+  CALL nc(NF90_PUT_ATT(fID, muID,'long_name','Std deviation of sub-cell scales orography'),'mu'  )
+  CALL nc(NF90_PUT_ATT(fID,sigID,'long_name','Slope along principal axis'                ),'Zsig')
+  CALL nc(NF90_PUT_ATT(fID,gamID,'long_name','Anisotropy (aspect ratio)'                 ),'Zgam')
+  CALL nc(NF90_PUT_ATT(fID,theID,'long_name','Orientation (principal axis)'              ),'Zthe')
+  CALL nc(NF90_PUT_ATT(fID,picID,'long_name','Maximum height'                            ),'Zpic')
+  CALL nc(NF90_PUT_ATT(fID,valID,'long_name','Minimum height'                            ),'Zval')
+
+  CALL nc(NF90_PUT_ATT(fID,mskID,'units','none'   ),'mask')
+  IF(fmsk=='noro') &
+  CALL nc(NF90_PUT_ATT(fID,phiID,'units','m'      ),'Zphi')
+  CALL nc(NF90_PUT_ATT(fID,meaID,'units','m'      ),'Zmea')
+  CALL nc(NF90_PUT_ATT(fID, muID,'units','m'      ),'mu'  )
+  CALL nc(NF90_PUT_ATT(fID,sigID,'units','m/m'    ),'Zsig')
+  CALL nc(NF90_PUT_ATT(fID,gamID,'units','none'   ),'Zgam')
+  CALL nc(NF90_PUT_ATT(fID,theID,'units','degrees'),'Zthe')
+  CALL nc(NF90_PUT_ATT(fID,picID,'units','m'      ),'Zpic')
+  CALL nc(NF90_PUT_ATT(fID,valID,'units','m'      ),'Zval')
+
+  CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Conventions','COARDS/CF-1.0'))
+  CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Initial_Grid',TRIM(res_in)))
+  CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'history',TRIM(call_seq)))
+  CALL nc(NF90_ENDDEF(fID))
+
+  CALL nc(NF90_PUT_VAR(fID, loID,lon_ou),'x' )
+  CALL nc(NF90_PUT_VAR(fID, laID,lat_ou),'y' )
+  CALL nc(NF90_PUT_VAR(fID,mskID,msko),'mask')
+  IF(fmsk=='noro') &
+  CALL nc(NF90_PUT_VAR(fID,phiID,Zphi),'Zphi')
+  CALL nc(NF90_PUT_VAR(fID,meaID,h0  ),'Zmea')
+  CALL nc(NF90_PUT_VAR(fID, muID,mu  ),'mu'  )
+  CALL nc(NF90_PUT_VAR(fID,sigID,Zsig),'Zsig')
+  CALL nc(NF90_PUT_VAR(fID,gamID,Zgam),'Zgam')
+  CALL nc(NF90_PUT_VAR(fID,theID,Zthe),'Zthe')
+  CALL nc(NF90_PUT_VAR(fID,picID,       h0+2*mu ),'Zpic')
+  CALL nc(NF90_PUT_VAR(fID,valID,MAX(0.,h0-2*mu)),'Zval')
+  CALL nc(NF90_CLOSE(fID))
+  WRITE(*,*)'Finished.'
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE diffusion_filter(a,b)
+!
+!-------------------------------------------------------------------------------
+! Purpose: ECMWF bi-laplacian diffusion filter to limit aliasing.
+! Note that one third of the frequencies is skipped to allow some accuracy for
+! derived quantities computation.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: a(:,:), b(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: m
+  REAL, ALLOCATABLE :: fn(:)
+!-------------------------------------------------------------------------------
+  ALLOCATE(fn(ndab)); fn=[(REAL(m*(m-1)),m=1,ndab)]
+  a(1,:)=a(1,:)/(1+4.*fn(:)/fn(ndab))**2
+  b(1,:)=b(1,:)/(1+4.*fn(:)/fn(ndab))**2
+  DO m=1,ndab
+    IF(m<=2*ndab/3) THEN
+      a(m,:)=a(m,:)/(1+5.*fn(:)/fn(ndab))**2
+      b(m,:)=b(m,:)/(1+5.*fn(:)/fn(ndab))**2
+    ELSE
+      a(m,:)=0.
+      b(m,:)=0.
+    END IF
+  END DO
+
+END SUBROUTINE diffusion_filter
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE lowpass_filter(a,b,d,dd)
+!
+!-------------------------------------------------------------------------------
+! Purpose: ECMWF lowpass filter.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: a(:,:), b(:,:)
+  REAL, INTENT(IN)    :: d, dd
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: m, na
+  REAL, ALLOCATABLE  :: h(:)
+!-------------------------------------------------------------------------------
+  na=SIZE(a,1)
+  IF(ALLOCATED(h)) THEN; IF(SIZE(h)/=na) DEALLOCATE(h); END IF
+  IF(.NOT.ALLOCATED(h)) ALLOCATE(h(na))
+  CALL lowpass(d,dd,h)
+  DO m=1,SIZE(h); a(m,:)=a(m,:)*h; b(m,:)=b(m,:)*h; END DO
+
+END SUBROUTINE lowpass_filter
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE bandpass_filter(a,b,d1,dd1,d2,dd2)
+!
+!-------------------------------------------------------------------------------
+! Purpose: ECMWF bandpass filter. Calls lowpass filter with 2 different scales.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: a(:,:), b(:,:)
+  REAL, INTENT(IN)    :: d1, dd1, d2, dd2
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: m, na
+  REAL, ALLOCATABLE  :: h(:), hh(:)
+!-------------------------------------------------------------------------------
+  na=SIZE(a,1)
+  IF(ALLOCATED(hh)) THEN; IF(SIZE(hh)/=na) DEALLOCATE(hh); END IF
+  IF(ALLOCATED(h )) THEN; IF(SIZE(h )/=na) DEALLOCATE(h ); END IF
+  IF(.NOT.ALLOCATED(hh)) ALLOCATE(hh(na))
+  IF(.NOT.ALLOCATED(h )) ALLOCATE(h (na))
+  CALL lowpass(d1,dd1,hh)
+  CALL lowpass(d2,dd2,h ); hh=hh-h
+  DO m=1,SIZE(h); a(m,:)=a(m,:)*hh; b(m,:)=b(m,:)*hh; END DO
+
+END SUBROUTINE bandpass_filter
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE lowpass(d,dd,h)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)    :: d, dd   ! d is DELTA/REarth ; dd is delta/REarth
+  REAL, INTENT(INOUT) :: h(:)
+!-------------------------------------------------------------------------------
+  INTEGER :: m, n
+  REAL    :: hd
+!-------------------------------------------------------------------------------
+  n=SIZE(h); hd=0.5*d
+  h(:)=[((SIN(m*(hd-dd))/m+(COS(hp+m*hd)*SIN(hp+m*dd))/(hp/dd+m)               &
+        + SIN(m*(hd+dd))/m+(COS(hp-m*hd)*SIN(hp-m*dd))/(hp/dd-m))**2/d/d,m=1,n)]
+
+END SUBROUTINE lowpass
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE spatial_filter_1(u,D,dd)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Spatial filter ; using spectral lowpass filter.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: u(:,:)
+  REAL, INTENT(IN)    :: D, dd
+!-------------------------------------------------------------------------------
+! Local variables:
+  REAL, ALLOCATABLE  :: uu(:,:,:)
+!-------------------------------------------------------------------------------
+  ALLOCATE(uu(SIZE(u,1),SIZE(u,2),1)); uu(:,:,1)=u
+  CALL spatial_filter_m(uu,D,dd)
+  u=uu(:,:,1)
+  DEALLOCATE(uu)
+
+END SUBROUTINE spatial_filter_1
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE spatial_filter_m(u,D,dd)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Spatial filter ; using spectral lowpass filter.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: u(:,:,:)
+  REAL, INTENT(IN)    :: D, dd
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: ku
+  REAL, ALLOCATABLE  :: aa(:,:,:),bb(:,:,:)
+!-------------------------------------------------------------------------------
+  ku=SIZE(u,3)                                   !--- NUMBER OF FIELDS TO SMOOTH
+  ALLOCATE(aa(mdab,ndab,ku),bb(mdab,ndab,ku))
+
+!--- PROJECTION OF THE SUB-CELL SCALES FIELDS ON THE SPHERICAL HARMONICS BASE.
+  IF(lfast) THEN
+    CALL shaec(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wsha,lsh,work,lwork,ierr)
+    CALL error('shaec',ierr)
+  ELSE
+    CALL shaes(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wsha,lsh,work,lwork,ierr)
+    CALL error('shaes',ierr)
+  END IF
+
+!--- LOW PASS FILTER AT THE OUTPUT GRID RESOLUTION
+  DO n=1,ku; CALL lowpass_filter(aa(:,:,n),bb(:,:,n),D,dd); END DO
+
+!--- BACKWARD TRANSFORM FOR MEAN SUB-CELL SCALES PARAMETERS
+  IF(lfast) THEN
+    CALL shsec(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wshs,lsh,work,lwork,ierr)
+    CALL error('shsec',ierr)
+  ELSE
+    CALL shses(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wshs,lsh,work,lwork,ierr)
+    CALL error('shses',ierr)
+  END IF
+  DEALLOCATE(aa,bb)
+
+END SUBROUTINE spatial_filter_m
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE grid_noro0(xin,yin,zin,xou,you,zphi,mask)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Sub-cell scales orographic parameters coomputation. Angles in radians
+! Notes: Input coordinates contain a longitude periodic redundant point.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: xin(:), yin(:)   !--- INPUT  COORDINATES (nxi)   (nyi)
+  REAL, INTENT(IN)  :: zin(:,:)         !--- INPUT  FIELD       (nxi,    nyi)
+  REAL, INTENT(IN)  :: xou(:), you(:)   !--- OUTPUT COORDINATES (nxo+1) (nyo)
+  REAL, INTENT(OUT) :: zphi(:,:)        !--- GEOPOTENTIAL       (nxo+1,  nyo)
+  REAL, INTENT(OUT) :: mask(:,:)        !--- FRACTIONAL MASK    (nxo+1,  nyo)
+!-------------------------------------------------------------------------------
+! Local variables:
+  REAL, ALLOCATABLE :: xi(:), xo_w(:), yo_s(:), zi(:,:), num_tot(:,:)
+  REAL, ALLOCATABLE :: yi(:), xo_e(:), yo_n(:),  w(:,:), num_lan(:,:)
+  REAL    :: dx, dy, wx, wy, zx, zy, ex, ey
+  INTEGER :: nxi, nxo, i, ii, nn
+  INTEGER :: nyi, nyo, j, jj, di
+!-------------------------------------------------------------------------------
+  nxi=SIZE(xin); nxo=SIZE(xou)-1; ex=SIGN(1.,xin(nxi)-xin(1))
+  nyi=SIZE(yin); nyo=SIZE(you);   ey=SIGN(1.,yin(nyi)-yin(1))
+
+!--- MARGIN TO ENSURE OUTPUT CELLS HAVE MATCHING INPUT CELLS
+  di=CEILING(1.+MAXVAL([xou(2:nxo+1)-xou(1:nxo)])/MINVAL([xin(2:nxi)-xin(1:nxi-1)]))
+
+!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES
+  ALLOCATE(xi(nxi+2*di)); xi(:)=[xin(nxi-di:nxi-1)-2.*xpi,xin,xin(2:di+1)+2.*xpi]
+  ALLOCATE(yi(nyi+2));    yi(:)=[2.*yin(1)-yin(2),yin,2.*yin(nyi)-yin(nyi-1)]
+
+  ALLOCATE(zi(nxi+2*di,nyi+2))
+  zi(1      +di:nxi+  di,2:nyi+1)=zin(:,:)
+  zi(1         :      di,2:nyi+1)=zin(nxi-di:nxi-1,:)
+  zi(1+nxi  +di:nxi+2*di,2:nyi+1)=zin(2     :di+1 ,:)
+  zi(1         :nxi/2+di,      1)=zi (1+nxi/2:nxi  +di,    2)
+  zi(1+nxi/2+di:nxi+2*di,      1)=zi (1      :nxi/2+di,    2)
+  zi(1         :nxi/2+di,  nyi+2)=zi (1+nxi/2:nxi  +di,nyi+1)
+  zi(1+nxi/2+di:nxi+2*di,  nyi+2)=zi (1      :nxi/2+di,nyi+1)
+
+!  write(*,'(a,2f10.5,a,2f10.5)')'xin=',xin(1:2),'...',xin(nxi-1:nxi)
+!  write(*,'(a,2f10.5,a,2f10.5)')'yin=',yin(1:2),'...',yin(nyi-1:nyi)
+!  write(*,'(a,2f10.5,a,2f10.5)')'xou=',xou(1:2),'...',xou(nxo:nxo+1)
+!  write(*,'(a,2f10.5,a,2f10.5)')'you=',you(1:2),'...',you(nyo-1:nyo)
+!  write(*,*)'xi=',xi(1:di+1),'...',xi(nxi+di:nxi+2*di)
+!  write(*,'(a,3f10.5,a,3f10.5)')'yi=',yi(1:3),'...',yi(nyi:nyi+2)
+
+!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) bdac
+  ALLOCATE(xo_w(nxo+1)); xo_w=(xou+[2.*xou(    1)-xou(  2),xou(1:nxo  )])/2.0
+  ALLOCATE(yo_s(nyo  )); yo_s=(you+[2.*you(    1)-you(  2),you(1:nyo-1)])/2.0
+  ALLOCATE(xo_e(nxo+1)); xo_e=(xou+[xou(2:nxo+1),2.*xou(nxo+1)-xou(nxo)])/2.0
+  ALLOCATE(yo_n(nyo  )); yo_n=(you+[you(2:nyo  ),2.*you(nyo)-you(nyo-1)])/2.0
+
+!--- INITIALIZATIONS:
+  ALLOCATE(w(nxo+1,nyo)); w(:,:)=0.0; zphi(:,:)=0.0
+
+!--- SUMMATION OVER GRIDPOINT AREA
+  zx=2.*xpi/REAL(nxi)                            !--- LON CELL LENGTH (1-SPHERE)
+  dx=zx/2.                                       !--- HALF "
+  ALLOCATE(num_tot(nxo+1,nyo)); num_tot(:,:)=0.
+  ALLOCATE(num_lan(nxo+1,nyo)); num_lan(:,:)=0.
+  DO ii=1,nxo+1
+    DO jj=1,nyo
+      DO j=2,nyi+1
+        zy=xpi/REAL(nyi)*COS(yi(j))              !--- LAT CELL LENGTH (1-SPHERE)
+        dy=zy/2.
+        wy=MAX(0.,MIN(zy,MINVAL(dy+ey*[yo_n(jj)-yi(j),yi(j)-yo_s(jj)])))
+        IF(wy==0.) CYCLE
+        DO i=2,nxi+2*di-1
+          wx=MAX(0.,MIN(zx,MINVAL(dx+ex*[xo_e(ii)-xi(i),xi(i)-xo_w(ii)])*COS(yi(j))))
+          IF(wx==0.) CYCLE
+          num_tot(ii,jj)=num_tot(ii,jj)+1.0
+          IF(zi(i,j)>=1.) num_lan(ii,jj)=num_lan(ii,jj)+1.0
+          w   (ii,jj)=   w(ii,jj)+wx*wy
+          zphi(ii,jj)=zphi(ii,jj)+wx*wy*zi(i,j)
+        END DO
+      END DO
+    END DO
+  END DO
+
+!--- COMPUTE FRACTIONAL MASK AND GEOPOTENTIAL
+  WHERE(w(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
+  nn=COUNT(w(:,:)==0.0)
+  IF(nn/=0) WRITE(*,*)'Problem with weight ; vanishing occurrences: ',nn
+  WHERE(w/=0.0) zphi=zphi/w
+  WHERE(w==0.0.OR.mask<0.1) zphi=0.0
+  zphi(nxo+1,:)=zphi(1,:)
+  zphi(:,  1)=SUM(w(1:nxo,  1)*zphi(1:nxo,  1))/SUM(w(1:nxo,  1))
+  zphi(:,nyo)=SUM(w(1:nxo,nyo)*zphi(1:nxo,nyo))/SUM(w(1:nxo,nyo))
+
+END SUBROUTINE grid_noro0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE nc(ncres,var)
+!
+!-------------------------------------------------------------------------------
+! Purpose: NetCDF errors handling.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,          INTENT(IN)           :: ncres
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: var
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: msg
+!-------------------------------------------------------------------------------
+  IF(ncres/=NF90_NoErr) THEN
+    msg='Error in routine '//TRIM(sub)
+    IF(fnam/='')     msg=TRIM(msg)//' for file "'//TRIM(fnam)//'"'
+    IF(PRESENT(var)) msg=TRIM(msg)//'" and variable "'//TRIM(var)//'"'
+    WRITE(*,*)TRIM(msg)//': '//NF90_STRERROR(ncres); STOP
+  END IF
+
+END SUBROUTINE nc
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+ELEMENTAL FUNCTION int2str(val)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: convert an integer into a string.
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=128)  :: int2str
+  INTEGER, INTENT(IN) :: val
+!-------------------------------------------------------------------------------
+! Local variable:
+  CHARACTER(LEN=128) :: s1
+!-------------------------------------------------------------------------------
+  WRITE(s1,*)val; int2str=TRIM(ADJUSTL(s1))
+
+END FUNCTION int2str
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE err(ierr,str)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Argument:
+  LOGICAL,          INTENT(IN) :: ierr
+  CHARACTER(LEN=*), INTENT(IN) :: str
+!-------------------------------------------------------------------------------
+  IF(.NOT.ierr) RETURN
+  WRITE(*,*)TRIM(str)
+  STOP
+
+END SUBROUTINE err
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE error(sub,ierr)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: sub
+  INTEGER,          INTENT(IN) :: ierr
+!-------------------------------------------------------------------------------
+  IF(ierr==0) RETURN
+  WRITE(*,*)'In routine "'//TRIM(sub)//'": error in the specification of argume&
+      &nt nr. '//TRIM(int2str(ierr))
+  STOP
+END SUBROUTINE error
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+ELEMENTAL FUNCTION real2str(val)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: convert an integer into a string.
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=128) :: real2str
+  REAL,   INTENT(IN) :: val
+  INTEGER :: r, s
+!-------------------------------------------------------------------------------
+                s=NINT(val); r=NINT(1000.*(val-REAL(s)))
+  IF(r<0) THEN; s= INT(val); r= INT(1000.*(val-REAL(s))); END IF
+  real2str=TRIM(ADJUSTL(int2str(s)))//'.'//TRIM(ADJUSTL(int2str(r)))
+
+END FUNCTION real2str
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE flip(a)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, ALLOCATABLE, INTENT(INOUT) :: a(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: k, n
+  REAL, ALLOCATABLE :: tmp(:,:)
+!-------------------------------------------------------------------------------
+  n=SIZE(a,2)
+  ALLOCATE(tmp(n,SIZE(a,1))); DO k=1,n; tmp(k,:)=a(:,k); END DO
+  CALL MOVE_ALLOC(FROM=tmp,TO=a)
+
+END SUBROUTINE flip
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+ELEMENTAL LOGICAL FUNCTION is_int(str) RESULT(out)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: tests if str is a valid integer number
+!-------------------------------------------------------------------------------
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: k, k0
+!-------------------------------------------------------------------------------
+  out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN
+  k0=1; IF(INDEX('+-',str(1:1))/=0) k0=2
+  DO k=k0,LEN_TRIM(str); IF(INDEX('0123456789',str(k:k))==0) out=.FALSE.; END DO
+END FUNCTION is_int
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+ELEMENTAL LOGICAL FUNCTION is_dec(str) RESULT(out)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: tests if str is a valid decimal number or an integer
+!-------------------------------------------------------------------------------
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: ix
+!-------------------------------------------------------------------------------
+  out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN
+
+  !--- TEST IF THE NUMBER IS A PURE INTEGER
+  IF(is_int(str)) RETURN
+
+  !--- CHECK FOR POINT SEPARATOR (SHOULD HAVE A STRING INT BEFORE AND/OR AFTER
+  ix=INDEX(str,'.'); out=ix>0.AND.LEN_TRIM(str)/=1; IF(.NOT.out) RETURN
+
+  !--- CHECK SURROUNDING STRINGS ARE INTEGERS
+  out=is_int(str(1:ix-1)); IF(ix==LEN_TRIM(str)) RETURN
+  out=out.AND.is_int(str(ix+1:LEN_TRIM(str)))
+END FUNCTION is_dec
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+ELEMENTAL LOGICAL FUNCTION is_num(str) RESULT(out)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: test if str is a valid integer/decimal/scientific notation number
+!-------------------------------------------------------------------------------
+  CHARACTER(LEN=*),  INTENT(IN) :: str
+  INTEGER :: ix, ls
+!-------------------------------------------------------------------------------
+  out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN; ls=LEN_TRIM(str)
+
+  !--- TEST IF THE NUMBER IS A PURE INTEGER OR A PURE DECIMAL NUMBER
+  IF(is_int(str)) RETURN
+  IF(is_dec(str)) RETURN
+
+  !--- CHECK FOR "D" OR "E" SEPARATOR
+  ix=MAX(INDEX(strLow(str),'d'),INDEX(strLow(str),'e')); out=ix>0
+  IF(.NOT.out) RETURN
+
+  !--- CHECK SURROUNDING STRINGS ARE NUMBERS
+  out=is_dec(str(1:ix-1)).AND.is_int(str(ix+1:ls))
+END FUNCTION is_num
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+INTEGER FUNCTION str2int(str) RESULT(out)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: convert a string into an integer
+!-------------------------------------------------------------------------------
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: k
+!-------------------------------------------------------------------------------
+  CALL err(.NOT.is_int(str),'Invalid integer "'//TRIM(str)//'"')
+  out=0
+  DO k=1,LEN_TRIM(str); out=10*out+INT(INDEX('0123456789',str(k:k))-1); END DO
+END FUNCTION str2int
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+REAL FUNCTION str2real(str) RESULT(out)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Convert a string into an real
+!-------------------------------------------------------------------------------
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: ix, l
+!-------------------------------------------------------------------------------
+  CALL err(.NOT.is_num(str),'Invalid number "'//TRIM(str)//'"')
+  ix=MAX(INDEX(strLow(str),'d'),INDEX(strLow(str),'e')); l=LEN_TRIM(str)
+  IF(ix==0) out=dec2real(str)
+  IF(ix/=0) out=dec2real(str(1:ix-1))*10.**str2int(str(ix+1:l))
+END FUNCTION str2real
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+ELEMENTAL REAL FUNCTION dec2real(str) RESULT(out)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: convert a valid (already checked) string into an real
+!-------------------------------------------------------------------------------
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: k, k0, ix
+!-------------------------------------------------------------------------------
+  out=0.
+  k0=1; IF(INDEX('+-',str(1:1))/=0) k0=2; ix=INDEX(str,'.')     !--- CHECK SIGN
+  DO k=k0,LEN_TRIM(str); IF(k==ix) CYCLE                        !--- SKIP .+-
+    out=out*10.+REAL(INDEX('0123456789',str(k:k))-1)
+  END DO
+  IF(ix/=0) out=out*10.**(ix-LEN_TRIM(str))                     !--- COMA
+  IF(k0==2) out=-out                                            !--- SIGN
+END FUNCTION dec2real
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+!
+ELEMENTAL CHARACTER(LEN=256) FUNCTION strLow(str) RESULT(out)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Lower case conversion.
+!-------------------------------------------------------------------------------
+! Parameters:
+  CHARACTER(LEN=*), INTENT(IN) :: str
+!-------------------------------------------------------------------------------
+! Local variable:
+  INTEGER :: k, ix
+!-------------------------------------------------------------------------------
+  out=str
+  DO k=1,LEN(str)
+    ix=IACHAR(str(k:k)); IF(64<ix.AND.ix<91) out(k:k)=ACHAR(ix+32)
+  END DO
+
+END FUNCTION strLow
+!
+!-------------------------------------------------------------------------------
+
+END PROGRAM make_sso
+!
+!-------------------------------------------------------------------------------
