Index: /trunk/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 5101)
+++ /trunk/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 5102)
@@ -1077,6 +1077,5 @@
 correctly set  ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read 
 using \mdl{fldread}, see \S\ref{SBC_fldread}). 
-The restoring coefficient $\gamma$ is a three-dimensional array initialized by the 
-user in routine \rou{dtacof} also located in module \mdl{tradmp}. 
+The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. The file name is specified by the namelist variable \np{cn\_resto}. The DMP\_TOOLS tool is provided to allow users to generate the netcdf file.
 
 The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} 
@@ -1092,23 +1091,5 @@
 diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity 
 field consistent with the model dynamics whilst having a $T$, $S$ field 
-close to a given climatological field ($T_o$, $S_o$). The time scale 
-associated with $S_o$ is generally not a constant but spatially varying 
-in order to respect other properties. For example, it is usually set to zero 
-in the mixed layer (defined either on a density or $S_o$ criterion) 
-\citep{Madec_al_JPO96} and in the equatorial region 
-\citep{Reverdin1991, Fujio1991, Marti_PhD92} since these two regions 
-have a short time scale of adjustment; while smaller $\gamma$ are used 
-in the deep ocean where the typical time scale is long \citep{Sarmiento1982}. 
-In addition the time scale is reduced (even to zero) along the western 
-boundary to allow the model to reconstruct its own western boundary 
-structure in equilibrium with its physics. 
-The choice of the shape of the Newtonian damping is controlled by two 
-namelist parameters \np{nn\_hdmp} and \np{nn\_zdmp}. The former allows us to specify: the 
-width of the equatorial band in which no damping is applied; a decrease 
-in the vicinity of the coast; and a damping everywhere in the Red and Med Seas.
-The latter sets whether damping should act in the mixed layer or not. 
-The time scale associated with the damping depends on the depth as
-a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as 
-bottom value and a transition depth of \np{rn\_dep}.  
+close to a given climatological field ($T_o$, $S_o$). 
 
 The robust diagnostic method is very efficient in preventing temperature 
@@ -1118,9 +1099,19 @@
 by stabilising the water column too much.
 
-An example of the computation of $\gamma$ for a robust diagnostic experiment 
-with the ORCA2 model is provided in the \mdl{tradmp} module 
-(subroutines \rou{dtacof} and \rou{cofdis} which compute the coefficient 
-and the distance to the bathymetry, respectively). These routines are 
-provided as examples and can be customised by the user. 
+The namelist parameter \np{nn\_zdmp} sets whether the damping should be applied in the whole water column or only below the mixed layer (defined either on a density or $S_o$ criterion). It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here \citep{Madec_al_JPO96}.
+
+\subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS}
+
+DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input. This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
+
+%--------------------------------------------nam_dmp_create-------------------------------------------------
+\namdisplay{nam_dmp_create}
+%-------------------------------------------------------------------------------------------------------
+
+\np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list.
+
+The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations. \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for the ORCA4, ORCA2 and ORCA05 configurations. If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference configurations with previous model versions. \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. This option only has an effect if \np{ln\_full\_field} is true. \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. Finally \np{ln\_custom} specifies that the custom module will be called. This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
+
+The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to the full values of a 10$^{\circ}$ latitud band. This is often used because of the short adjustment time scale in the equatorial region \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.  
 
 % ================================================================
Index: /trunk/DOC/TexFiles/Namelist/nam_dmp_create
===================================================================
--- /trunk/DOC/TexFiles/Namelist/nam_dmp_create	(revision 5102)
+++ /trunk/DOC/TexFiles/Namelist/nam_dmp_create	(revision 5102)
@@ -0,0 +1,23 @@
+&nam_dmp_create
+    cp_cfg = 'orca'            ! Name of model grid (orca and C1D have special options 
+                               ! otherwise ignored)
+    cp_cfz = 'antarctic'       ! Name of zoom configuration (arctic and antarctic have 
+                               ! some special treatment if lzoom=.true.)
+    jp_cfg = 2                 ! Resolution of the model (used for med_red_seas damping)
+    lzoom = .false.            ! Zoom configuration or not
+    ln_full_field = .false.    ! Calculate coefficient over whole of domain
+    ln_med_red_seas = .true.   ! Damping in Med/Red Seas 
+                               ! (or local modifications here if ln_full_field=.true.)
+    ln_old_31_lev_code = .true.   ! Replicate behaviour of old online code for 31 level model 
+                                  ! (Med/Red seas damping based on level number instead of depth)
+    ln_coast = .true.          ! Reduce near to coastlines
+    ln_zero_top_layer = .true. ! No damping in top layer
+    ln_custom = .false.        ! Call "custom" module to apply user modifications to the 
+                               ! damping coefficient field
+    pn_surf = 0.25             ! Surface Relaxation timescale (days)
+    pn_bot = 0.25              ! Bottom relaxation timescale (days)
+    pn_dep = 1000              ! Transition depth from upper to deep ocean
+    nn_hdmp = 10               ! Damp poleward of this latitude (smooth transition up to maximum damping)
+    jperio = 2                 ! Lateral boundary condition (as specified in namelist_cfg for model run).
+/
+
Index: /trunk/DOC/TexFiles/Namelist/namtra_dmp
===================================================================
--- /trunk/DOC/TexFiles/Namelist/namtra_dmp	(revision 5101)
+++ /trunk/DOC/TexFiles/Namelist/namtra_dmp	(revision 5102)
@@ -2,14 +2,9 @@
 &namtra_dmp    !   tracer: T & S newtonian damping
 !-----------------------------------------------------------------------
-   ln_tradmp   =  .true.   !  add a damping termn (T) or not (F)
-   nn_hdmp     =   -1      !  horizontal shape =-1, damping in Med and Red Seas only
-                           !                   =XX, damping poleward of XX degrees (XX>0)
-                           !                      + F(distance-to-coast) + Red and Med Seas
-   nn_zdmp     =    0      !  vertical   shape =0    damping throughout the water column
-                           !                   =1 no damping in the mixing layer (kz  criteria)
-                           !                   =2 no damping in the mixed  layer (rho crieria)
-   rn_surf     =   50.     !  surface time scale of damping   [days]
-   rn_bot      =  360.     !  bottom  time scale of damping   [days]
-   rn_dep      =  800.     !  depth of transition between rn_surf and rn_bot [meters]
-   nn_file     =    0      !  create a damping.coeff NetCDF file (=1) or not (=0)
+   ln_tradmp   =  .true.     !  add a damping termn (T) or not (F)
+   nn_zdmp     =    0        !  vertical   shape =0    damping throughout the water column
+                             !                   =1 no damping in the mixing layer (kz  criteria)
+                             !                   =2 no damping in the mixed  layer (rho crieria)
+   cn_resto    = 'resto.nc'  ! Name of file containing restoration coefficient field (use dmp_tools to create this)
+
 /
Index: /trunk/NEMOGCM/ARCH/OLD/arch-ifort_linux.fcm
===================================================================
--- /trunk/NEMOGCM/ARCH/OLD/arch-ifort_linux.fcm	(revision 5101)
+++ /trunk/NEMOGCM/ARCH/OLD/arch-ifort_linux.fcm	(revision 5102)
@@ -17,5 +17,5 @@
 %NCDF_INC            -I/usr/local/netcdf/include 
 %NCDF_LIB            -L /usr/local/netcdf/lib -lnetcdf   
-%FC                  ifort
+%FC                  ifort 
 %FCFLAGS 	     -r8 -O3  -traceback 
 %FFLAGS 	     -r8 -O3  -traceback 
Index: /trunk/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg
===================================================================
--- /trunk/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 5102)
@@ -241,6 +241,4 @@
 !-----------------------------------------------------------------------
    ln_tradmp   =  .false.  !  add a damping termn (T) or not (F)
-   nn_zdmp     =   -1      !  vertical   shape =0    damping throughout the water column
-   nn_file     =    1      !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg
===================================================================
--- /trunk/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg	(revision 5102)
@@ -245,6 +245,4 @@
 !-----------------------------------------------------------------------
    ln_tradmp   =  .false.  !  add a damping termn (T) or not (F)
-   nn_zdmp     =    1      !  vertical   shape =0    damping throughout the water column
-   nn_file     =    1      !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
===================================================================
--- /trunk/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 5102)
@@ -244,6 +244,4 @@
 !-----------------------------------------------------------------------
    ln_tradmp   =  .false.  !  add a damping termn (T) or not (F)
-   nn_zdmp     =    1      !  vertical   shape =0    damping throughout the water column
-   nn_file     =    1      !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
===================================================================
--- /trunk/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 5102)
@@ -155,6 +155,4 @@
 !-----------------------------------------------------------------------
    ln_tradmp   =  .false.  !  add a damping termn (T) or not (F)
-   nn_zdmp     =    1      !  vertical   shape =0    damping throughout the water column
-   nn_file     =    1      !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/CONFIG/GYRE_XIOS/EXP00/namelist_cfg
===================================================================
--- /trunk/NEMOGCM/CONFIG/GYRE_XIOS/EXP00/namelist_cfg	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/GYRE_XIOS/EXP00/namelist_cfg	(revision 5102)
@@ -216,6 +216,4 @@
 !-----------------------------------------------------------------------
    ln_tradmp   =  .false.  !  add a damping termn (T) or not (F)
-   nn_zdmp     =    1      !  vertical   shape =0    damping throughout the water column
-   nn_file     =    1      !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /trunk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5102)
@@ -791,14 +791,8 @@
 !-----------------------------------------------------------------------
    ln_tradmp   =  .true.   !  add a damping termn (T) or not (F)
-   nn_hdmp     =   -1      !  horizontal shape =-1, damping in Med and Red Seas only
-                           !                   =XX, damping poleward of XX degrees (XX>0)
-                           !                      + F(distance-to-coast) + Red and Med Seas
    nn_zdmp     =    0      !  vertical   shape =0    damping throughout the water column
                            !                   =1 no damping in the mixing layer (kz  criteria)
                            !                   =2 no damping in the mixed  layer (rho crieria)
-   rn_surf     =   50.     !  surface time scale of damping   [days]
-   rn_bot      =  360.     !  bottom  time scale of damping   [days]
-   rn_dep      =  800.     !  depth of transition between rn_surf and rn_bot [meters]
-   nn_file     =    0      !  create a damping.coeff NetCDF file (=1) or not (=0)
+   cn_resto    = 'resto.nc' ! Name of file containing restoration coefficient field (use dmp_tools to create this)
 /
 
Index: /trunk/NEMOGCM/CONFIG/SHARED/namelist_top_ref
===================================================================
--- /trunk/NEMOGCM/CONFIG/SHARED/namelist_top_ref	(revision 5101)
+++ /trunk/NEMOGCM/CONFIG/SHARED/namelist_top_ref	(revision 5102)
@@ -76,14 +76,8 @@
 &namtrc_dmp    !   passive tracer newtonian damping   
 !-----------------------------------------------------------------------
-   nn_hdmp_tr  =   -1      !  horizontal shape =-1, damping in Med and Red Seas only
-                           !                   =XX, damping poleward of XX degrees (XX>0)
-                           !                      + F(distance-to-coast) + Red and Med Seas
    nn_zdmp_tr  =    1      !  vertical   shape =0    damping throughout the water column
                            !                   =1 no damping in the mixing layer (kz  criteria)
                            !                   =2 no damping in the mixed  layer (rho crieria)
-   rn_surf_tr  =   50.     !  surface time scale of damping   [days]
-   rn_bot_tr   =  360.     !  bottom  time scale of damping   [days]
-   rn_dep_tr   =  800.     !  depth of transition between rn_surf and rn_bot [meters]
-   nn_file_tr  =    0      !  create a damping.coeff NetCDF file (=1) or not (=0)
+   cn_resto_tr  = 'resto_tr.nc'    !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90	(revision 5101)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90	(revision 5102)
@@ -3,6 +3,9 @@
    !!                       ***  MODULE  dyndmp  ***
    !! Ocean dynamics: internal restoring trend on momentum (U and V current)
+   !!                 This should only be used for C1D case in current form
    !!======================================================================
    !! History :  3.5  ! 2013-08  (D. Calvert)  Original code
+   !!            3.6  ! 2014-08  (T. Graham) Modified to use netcdf file of
+   !!                             restoration coefficients supplied to tradmp
    !!----------------------------------------------------------------------
 
@@ -25,4 +28,5 @@
    USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
+   USE iom            ! I/O manager
 
    IMPLICIT NONE
@@ -73,4 +77,5 @@
       NAMELIST/namc1d_dyndmp/ ln_dyndmp
       INTEGER :: ios
+      INTEGER :: imask
       !!----------------------------------------------------------------------
 
@@ -91,10 +96,7 @@
          WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp
          WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters'
-         WRITE(numout,*) '      horizontal damping option       nn_hdmp   = ', nn_hdmp
-         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)'
-         WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf
-         WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot
-         WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep
-         WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file
+         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
+         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp
+         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto
          WRITE(numout,*)
       ENDIF
@@ -104,14 +106,4 @@
          IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
          !
-#if ! defined key_c1d
-         SELECT CASE ( nn_hdmp )             !==   control print of horizontal option   ==!
-         CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   momentum damping in the Med & Red seas only'
-         CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   momentum damping poleward of', nn_hdmp, ' degrees'
-         CASE DEFAULT
-            WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp
-            CALL ctl_stop(ctmp1)
-         END SELECT
-         !
-#endif
          SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==!
          CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column'
@@ -130,13 +122,9 @@
          utrdmp(:,:,:) = 0._wp               ! internal damping trends
          vtrdmp(:,:,:) = 0._wp
-         !                                   !==   Damping coefficients calculation:                           ==!
-         !                                   !==   use tradmp.F90 subroutines dtacof, dtacof_zoom and cofdis   ==!
-         !                                   !!!   NOTE: these need to be altered for use in this module if 
-         !                                   !!!       they are to be used outside the C1D context 
-         !                                   !!!       (use of U,V grid variables)
-         IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto_uv )
-         ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'DYN', resto_uv )
-         ENDIF
-         !
+         !
+         !Read in mask from file
+         CALL iom_open ( cn_resto, imask)
+         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
+         CALL iom_close( imask )
       ENDIF
       !
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 5101)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 5102)
@@ -21,7 +21,4 @@
    !!   tra_dmp       : update the tracer trend with the internal damping
    !!   tra_dmp_init  : initialization, namlist read, parameters control
-   !!   dtacof_zoom   : restoring coefficient for zoom domain
-   !!   dtacof        : restoring coefficient for global domain
-   !!   cofdis        : compute the distance to the coastline
    !!----------------------------------------------------------------------
    USE oce            ! ocean: variables
@@ -39,4 +36,5 @@
    USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -45,17 +43,12 @@
    PUBLIC   tra_dmp      ! routine called by step.F90
    PUBLIC   tra_dmp_init ! routine called by opa.F90
-   PUBLIC   dtacof       ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90
-   PUBLIC   dtacof_zoom  ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90
-
-!!gm  why all namelist variable public????   only ln_tradmp should be sufficient
 
    !                               !!* Namelist namtra_dmp : T & S newtonian damping *
+   ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90
    LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag
-   INTEGER , PUBLIC ::   nn_hdmp     ! = 0/-1/'latitude' for damping over T and S
    INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer
-   REAL(wp), PUBLIC ::   rn_surf     ! surface time scale for internal damping        [days]
-   REAL(wp), PUBLIC ::   rn_bot      ! bottom time scale for internal damping         [days]
-   REAL(wp), PUBLIC ::   rn_dep      ! depth of transition between rn_surf and rn_bot [meters]
-   INTEGER , PUBLIC ::   nn_file     ! = 1 create a damping.coeff NetCDF file 
+   CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field
+   !
+
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   strdmp   !: damping salinity trend (psu/s)
@@ -197,623 +190,60 @@
       !! ** Method  :   read the namtra_dmp namelist and check the parameters
       !!----------------------------------------------------------------------
-      INTEGER  ::   ios   ! Local integer output status for namelist read
-      !!
-      NAMELIST/namtra_dmp/ ln_tradmp, nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file
-      !!----------------------------------------------------------------------
-      !
-      REWIND( numnam_ref )              ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term
+      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto
+      INTEGER ::  ios         ! Local integer for output status of namelist read
+      INTEGER :: imask        ! File handle 
+      !!
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation
       READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp )
       !
-      REWIND( numnam_cfg )              ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term
+      REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation
       READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namtra_dmp )
-      
-      IF( lzoom .AND. .NOT. lk_c1d )   nn_zdmp = 0          ! restoring to climatology at closed north or south boundaries
-
-      IF(lwp) THEN                       ! Namelist print
+
+      IF(lwp) THEN                 !Namelist print
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_dmp_init : T and S newtonian damping'
+         WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation'
          WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '   Namelist namtra_dmp : set damping parameter'
-         WRITE(numout,*) '      add a damping term or not       ln_tradmp = ', ln_tradmp
-         WRITE(numout,*) '      T and S damping option          nn_hdmp   = ', nn_hdmp
-         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)'
-         WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf
-         WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot
-         WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep
-         WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file
+         WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters'
+         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
+         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp
+         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto
          WRITE(numout,*)
       ENDIF
 
-      IF( ln_tradmp ) THEN               ! initialization for T-S damping
-         !
+      IF( ln_tradmp) THEN
+         !
+         !Allocate arrays
          IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' )
-         !
-!!gm  I don't understand the specificities of c1d case......   
-!!gm  to be check with the autor of these lines
-         
-#if ! defined key_c1d
-         SELECT CASE ( nn_hdmp )
-         CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only'
-         CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees'
-         CASE DEFAULT
-            WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp
-            CALL ctl_stop(ctmp1)
+
+         !Check values of nn_zdmp
+         SELECT CASE (nn_zdmp)
+         CASE ( 0 )  ; IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask'
+         CASE ( 1 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline'
+         CASE ( 2 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
          END SELECT
-         !
-#endif
-         SELECT CASE ( nn_zdmp )
-         CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column'
-         CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)'
-         CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
-         CASE DEFAULT
-            WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp
-            CALL ctl_stop(ctmp1)
-         END SELECT
-         !
+
+         !TG: Initialisation of dtatsd - Would it be better to have dmpdta routine
+         !so can damp to something other than intitial conditions files?
          IF( .NOT.ln_tsd_tradmp ) THEN
             CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' )
             CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data
          ENDIF
-         !
-         strdmp(:,:,:) = 0._wp       ! internal damping salinity trend (used in asmtrj)
+
+         !initialise arrays - Are these actually used anywhere else?
+         strdmp(:,:,:) = 0._wp
          ttrdmp(:,:,:) = 0._wp
-         !                          ! Damping coefficients initialization
-         IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto )
-         ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto )
-         ENDIF
-         !
-      ENDIF
-      !
+
+         !Read in mask from file
+         CALL iom_open ( cn_resto, imask)
+         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
+         CALL iom_close( imask )
+       ENDIF
+
    END SUBROUTINE tra_dmp_init
 
-
-   SUBROUTINE dtacof_zoom( presto )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dtacof_zoom  ***
-      !!
-      !! ** Purpose :   Compute the damping coefficient for zoom domain
-      !!
-      !! ** Method  : - set along closed boundary due to zoom a damping over
-      !!                6 points with a max time scale of 5 days.
-      !!              - ORCA arctic/antarctic zoom: set the damping along
-      !!                south/north boundary over a latitude strip.
-      !!
-      !! ** Action  : - resto, the damping coeff. for T and S
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::   presto   ! restoring coeff. (s-1)
-      !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar
-      REAL(wp), DIMENSION(6)  ::   zfact               ! 1Dworkspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'dtacof_zoom')
-      !
-
-      zfact(1) =  1._wp
-      zfact(2) =  1._wp
-      zfact(3) = 11._wp / 12._wp
-      zfact(4) =  8._wp / 12._wp
-      zfact(5) =  4._wp / 12._wp
-      zfact(6) =  1._wp / 12._wp
-      zfact(:) = zfact(:) / ( 5._wp * rday )    ! 5 days max restoring time scale
-
-      presto(:,:,:) = 0._wp
-
-      ! damping along the forced closed boundary over 6 grid-points
-      DO jn = 1, 6
-         IF( lzoom_w )   presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : )                    = zfact(jn)   ! west  closed
-         IF( lzoom_s )   presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : )                    = zfact(jn)   ! south closed 
-         IF( lzoom_e )   presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn)   ! east  closed 
-         IF( lzoom_n )   presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn)   ! north closed
-      END DO
-
-      !                                           ! ====================================================
-      IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN   !  ORCA configuration : arctic or antarctic zoom
-         !                                        ! ====================================================
-         IF(lwp) WRITE(numout,*)
-         IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom'
-         IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom'
-         IF(lwp) WRITE(numout,*)
-         !
-         !                          ! Initialization : 
-         presto(:,:,:) = 0._wp
-         zlat0 = 10._wp                     ! zlat0 : latitude strip where resto decreases
-         zlat1 = 30._wp                     ! zlat1 : resto = 1 before zlat1
-         zlat2 = zlat1 + zlat0              ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
-         z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days
-
-         DO jk = 2, jpkm1           ! Compute arrays resto ; value for internal damping : 5 days
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zlat = ABS( gphit(ji,jj) )
-                  IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
-                     presto(ji,jj,jk) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  ) 
-                  ELSEIF( zlat < zlat1 ) THEN
-                     presto(ji,jj,jk) = z1_5d
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-      !                             ! Mask resto array
-      presto(:,:,:) = presto(:,:,:) * tmask(:,:,:)
-      !
-      IF( nn_timing == 1 )  CALL timing_stop( 'dtacof_zoom')
-      !
-   END SUBROUTINE dtacof_zoom
-
-
-   SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep,  &
-      &               kn_file, cdtype , presto           )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dtacof  ***
-      !!
-      !! ** Purpose :   Compute the damping coefficient
-      !!
-      !! ** Method  :   Arrays defining the damping are computed for each grid
-      !!                point for temperature and salinity (resto)
-      !!                Damping depends on distance to coast, depth and latitude
-      !!
-      !! ** Action  : - resto, the damping coeff. for T and S
-      !!----------------------------------------------------------------------
-      USE iom
-      USE ioipsl
-      !!
-      INTEGER                         , INTENT(in   )  ::  kn_hdmp    ! damping option
-      REAL(wp)                        , INTENT(in   )  ::  pn_surf    ! surface time scale (days)
-      REAL(wp)                        , INTENT(in   )  ::  pn_bot     ! bottom time scale (days)
-      REAL(wp)                        , INTENT(in   )  ::  pn_dep     ! depth of transition (meters)
-      INTEGER                         , INTENT(in   )  ::  kn_file    ! save the damping coef on a file or not
-      CHARACTER(len=3)                , INTENT(in   )  ::  cdtype     ! =TRA, TRC or DYN (tracer/dynamics indicator)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::  presto     ! restoring coeff. (s-1)
-      !
-      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1          ! local integers
-      INTEGER  ::   inum0, icot                 !   -       -
-      REAL(wp) ::   zinfl, zlon                 ! local scalars
-      REAL(wp) ::   zlat, zlat0, zlat1, zlat2   !   -      -
-      REAL(wp) ::   zsdmp, zbdmp                !   -      -
-      CHARACTER(len=20)                   :: cfile
-      REAL(wp), POINTER, DIMENSION(:    ) :: zhfac 
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zmrs 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdct 
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dtacof')
-      !
-      CALL wrk_alloc( jpk, zhfac          )
-      CALL wrk_alloc( jpi, jpj, zmrs      )
-      CALL wrk_alloc( jpi, jpj, jpk, zdct )
-#if defined key_c1d
-      !                                   ! ====================
-      !                                   !  C1D configuration : local domain
-      !                                   ! ====================
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '              dtacof : C1D 3x3 local domain'
-      IF(lwp) WRITE(numout,*) '              -----------------------------'
-      !
-      presto(:,:,:) = 0._wp
-      !
-      zsdmp = 1._wp / ( pn_surf * rday )
-      zbdmp = 1._wp / ( pn_bot  * rday )
-      DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !   ONLY vertical variation from zsdmp (sea surface) to zbdmp (bottom)
-               presto(ji,jj,jk) = zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)
-            END DO
-         END DO
-      END DO
-      !
-      presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)
-#else
-      !                                   ! ====================
-      !                                   !  ORCA configuration : global domain
-      !                                   ! ====================
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '              dtacof : Global domain of ORCA'
-      IF(lwp) WRITE(numout,*) '              ------------------------------'
-      !
-      presto(:,:,:) = 0._wp
-      !
-      IF( kn_hdmp > 0 ) THEN      !  Damping poleward of 'nn_hdmp' degrees  !
-         !                        !-----------------------------------------!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '              Damping poleward of ', kn_hdmp, ' deg.'
-         !
-         CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )
-         !
-         IF( icot > 0 ) THEN          ! distance-to-coast read in file
-            CALL iom_get  ( icot, jpdom_data, 'Tcoast', zdct )
-            CALL iom_close( icot )
-         ELSE                         ! distance-to-coast computed and saved in file (output in zdct)
-            CALL cofdis( zdct )
-         ENDIF
-
-         !                            ! Compute arrays resto 
-         zinfl = 1000.e3_wp                ! distance of influence for damping term
-         zlat0 = 10._wp                    ! latitude strip where resto decreases
-         zlat1 = REAL( kn_hdmp )           ! resto = 0 between -zlat1 and zlat1
-         zlat2 = zlat1 + zlat0             ! resto increases from 0 to 1 between |zlat1| and |zlat2|
-
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zlat = ABS( gphit(ji,jj) )
-               IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
-                  presto(ji,jj,1) = 0.5_wp * (  1._wp - COS( rpi*(zlat-zlat1)/zlat0 )  )
-               ELSEIF ( zlat > zlat2 ) THEN
-                  presto(ji,jj,1) = 1._wp
-               ENDIF
-            END DO
-         END DO
-
-         IF ( kn_hdmp == 20 ) THEN       ! North Indian ocean (20N/30N x 45E/100E) : resto=0
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zlat = gphit(ji,jj)
-                  zlon = MOD( glamt(ji,jj), 360._wp )
-                  IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45._wp < zlon .AND. zlon < 100._wp ) THEN
-                     presto(ji,jj,1) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         ENDIF
-
-         zsdmp = 1._wp / ( pn_surf * rday )
-         zbdmp = 1._wp / ( pn_bot  * rday )
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )
-                  !   ... Decrease the value in the vicinity of the coast
-                  presto(ji,jj,jk) = presto(ji,jj,1 ) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj,jk)/zinfl)  )
-                  !   ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)
-                  presto(ji,jj,jk) = presto(ji,jj,jk) * (  zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)  )
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-
-      !                                  ! =========================
-      !                                  !  Med and Red Sea damping    (ORCA configuration only)
-      !                                  ! =========================
-      IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN
-         IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas'
-         !
-         zmrs(:,:) = 0._wp
-         !
-         SELECT CASE ( jp_cfg )
-         !                                           ! =======================
-         CASE ( 4 )                                  !  ORCA_R4 configuration 
-            !                                        ! =======================
-            ij0 =  50   ;   ij1 =  56                    ! Mediterranean Sea
-
-            ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            ij0 =  50   ;   ij1 =  55
-            ii0 =  75   ;   ii1 =  80   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            ij0 =  52   ;   ij1 =  53
-            ii0 =  70   ;   ii1 =  74   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
-            DO jk = 1, 17
-               zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday
-            END DO
-            DO jk = 18, jpkm1
-               zhfac (jk) = 1._wp / rday
-            END DO
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration 
-            !                                        ! =======================
-            ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea
-            ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            ij0 = 100   ;   ij1 = 110
-            ii0 = 144   ;   ii1 = 156   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            ij0 = 100   ;   ij1 = 103
-            ii0 = 139   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            !
-            ij0 = 101   ;   ij1 = 102                    ! Decrease before Gibraltar Strait
-            ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp
-            ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp
-            ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp
-            ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp
-            !
-            ij0 =  87   ;   ij1 =  96                    ! Red Sea
-            ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            !
-            ij0 =  91   ;   ij1 =  91                    ! Decrease before Bab el Mandeb Strait
-            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80_wp
-            ij0 =  90   ;   ij1 =  90
-            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp
-            ij0 =  89   ;   ij1 =  89
-            ii0 = 158   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp
-            ij0 =  88   ;   ij1 =  88
-            ii0 = 160   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp
-            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
-            DO jk = 1, 17
-               zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday
-            END DO
-            DO jk = 18, jpkm1
-               zhfac (jk) = 1._wp / rday
-            END DO
-            !                                        ! =======================
-         CASE ( 05 )                                 !  ORCA_R05 configuration
-            !                                        ! =======================
-            ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea
-            ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            ii0 = 575   ;   ii1 = 658
-            ij0 = 314   ;   ij1 = 366   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            !
-            ii0 = 641   ;   ii1 = 651                    ! Black Sea (remaining part
-            ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            !
-            ij0 = 324   ;   ij1 = 333                    ! Decrease before Gibraltar Strait
-            ii0 = 565   ;   ii1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp
-            ii0 = 566   ;   ii1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp
-            ii0 = 567   ;   ii1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp
-            !
-            ii0 = 641   ;   ii1 = 665                    ! Red Sea
-            ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
-            !
-            ii0 = 666   ;   ii1 = 675                    ! Decrease before Bab el Mandeb Strait
-            ij0 = 270   ;   ij1 = 290   
-            DO ji = mi0(ii0), mi1(ii1)
-               zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) )
-            END DO 
-            zsdmp = 1._wp / ( pn_surf * rday )
-            zbdmp = 1._wp / ( pn_bot  * rday )
-            DO jk = 1, jpk
-               zhfac(jk) = (  zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep )  )
-            END DO
-            !                                       ! ========================
-         CASE ( 025 )                               !  ORCA_R025 configuration 
-            !                                       ! ========================
-            CALL ctl_stop( ' Not yet implemented in ORCA_R025' )
-            !
-         END SELECT
-
-         DO jk = 1, jpkm1
-            presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1._wp - zmrs(:,:) ) * presto(:,:,jk)
-         END DO
-
-         ! Mask resto array and set to 0 first and last levels
-         presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)
-         presto(:,:, 1 ) = 0._wp
-         presto(:,:,jpk) = 0._wp
-         !                         !--------------------!
-      ELSE                         !     No damping     !
-         !                         !--------------------!
-         CALL ctl_stop( 'Choose a correct value of nn_hdmp or put ln_tradmp to FALSE' )
-      ENDIF
-#endif
-
-      !                            !--------------------------------!
-      IF( kn_file == 1 ) THEN      !  save damping coef. in a file  !
-         !                         !--------------------------------!
-         IF(lwp) WRITE(numout,*) '              create damping.coeff.nc file'
-         IF( cdtype == 'TRA' ) cfile = 'damping.coeff'
-         IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc'
-         IF( cdtype == 'DYN' ) cfile = 'damping.coeff.dyn'
-         cfile = TRIM( cfile )
-         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
-         CALL iom_rstput( 0, 0, inum0, 'Resto', presto )
-         CALL iom_close ( inum0 )
-      ENDIF
-      !
-      CALL wrk_dealloc( jpk, zhfac)
-      CALL wrk_dealloc( jpi, jpj, zmrs )
-      CALL wrk_dealloc( jpi, jpj, jpk, zdct )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dtacof')
-      !
-   END SUBROUTINE dtacof
-
-
-   SUBROUTINE cofdis( pdct )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE cofdis  ***
-      !!
-      !! ** Purpose :   Compute the distance between ocean T-points and the
-      !!      ocean model coastlines. Save the distance in a NetCDF file.
-      !!
-      !! ** Method  :   For each model level, the distance-to-coast is 
-      !!      computed as follows : 
-      !!       - The coastline is defined as the serie of U-,V-,F-points
-      !!      that are at the ocean-land bound.
-      !!       - For each ocean T-point, the distance-to-coast is then 
-      !!      computed as the smallest distance (on the sphere) between the 
-      !!      T-point and all the coastline points.
-      !!       - For land T-points, the distance-to-coast is set to zero.
-      !!      C A U T I O N : Computation not yet implemented in mpp case.
-      !!
-      !! ** Action  : - pdct, distance to the coastline (argument)
-      !!              - NetCDF file 'dist.coast.nc' 
-      !!----------------------------------------------------------------------
-      USE ioipsl      ! IOipsl librairy
-      !!
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   pdct   ! distance to the coastline
-      !!
-      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
-      INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers
-      CHARACTER (len=32) ::   clname                     ! local name
-      REAL(wp) ::   zdate0                               ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask
-      REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace
-      LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('cofdis')
-      !
-      CALL wrk_alloc( jpi, jpj , zxt, zyt, zzt, zmask    )
-      CALL wrk_alloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     )
-      ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj)  )
-      !
-      IF( lk_mpp    )   CALL mpp_sum( ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'cofdis: requested local arrays unavailable')
-
-      ! 0. Initialization
-      ! -----------------
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'
-      IF(lwp) WRITE(numout,*) '~~~~~~'
-      IF(lwp) WRITE(numout,*)
-      IF( lk_mpp ) &
-           & CALL ctl_stop('         Computation not yet implemented with key_mpp_...', &
-           &               '         Rerun the code on another computer or ', &
-           &               '         create the "dist.coast.nc" file using IDL' )
-
-      pdct(:,:,:) = 0._wp
-      zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )
-      zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )
-      zzt(:,:) = SIN( rad * gphit(:,:) )
-
-
-      ! 1. Loop on vertical levels
-      ! --------------------------
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         ! Define the coastline points (U, V and F)
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zmask(ji,jj) =  ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
-                   &           + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) )
-               llcotu(ji,jj) = ( tmask(ji,jj,  jk) + tmask(ji+1,jj  ,jk) == 1._wp ) 
-               llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1._wp ) 
-               llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp )
-            END DO
-         END DO
-
-         ! Lateral boundaries conditions
-         llcotu(:, 1 ) = umask(:,  2  ,jk) == 1
-         llcotu(:,jpj) = umask(:,jpjm1,jk) == 1
-         llcotv(:, 1 ) = vmask(:,  2  ,jk) == 1
-         llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1
-         llcotf(:, 1 ) = fmask(:,  2  ,jk) == 1
-         llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1
-
-         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
-            llcotu( 1 ,:) = llcotu(jpim1,:)
-            llcotu(jpi,:) = llcotu(  2  ,:)
-            llcotv( 1 ,:) = llcotv(jpim1,:)
-            llcotv(jpi,:) = llcotv(  2  ,:)
-            llcotf( 1 ,:) = llcotf(jpim1,:)
-            llcotf(jpi,:) = llcotf(  2  ,:)
-         ELSE
-            llcotu( 1 ,:) = umask(  2  ,:,jk) == 1
-            llcotu(jpi,:) = umask(jpim1,:,jk) == 1
-            llcotv( 1 ,:) = vmask(  2  ,:,jk) == 1
-            llcotv(jpi,:) = vmask(jpim1,:,jk) == 1
-            llcotf( 1 ,:) = fmask(  2  ,:,jk) == 1
-            llcotf(jpi,:) = fmask(jpim1,:,jk) == 1
-         ENDIF
-         IF( nperio == 3 .OR. nperio == 4 ) THEN
-            DO ji = 1, jpim1
-               iju = jpi - ji + 1
-               llcotu(ji,jpj  ) = llcotu(iju,jpj-2)
-               llcotf(ji,jpjm1) = llcotf(iju,jpj-2)
-               llcotf(ji,jpj  ) = llcotf(iju,jpj-3)
-            END DO
-            DO ji = jpi/2, jpim1
-               iju = jpi - ji + 1
-               llcotu(ji,jpjm1) = llcotu(iju,jpjm1)
-            END DO
-            DO ji = 2, jpi
-               ijt = jpi - ji + 2
-               llcotv(ji,jpjm1) = llcotv(ijt,jpj-2)
-               llcotv(ji,jpj  ) = llcotv(ijt,jpj-3)
-            END DO
-         ENDIF
-         IF( nperio == 5 .OR. nperio == 6 ) THEN
-            DO ji = 1, jpim1
-               iju = jpi - ji
-               llcotu(ji,jpj  ) = llcotu(iju,jpjm1)
-               llcotf(ji,jpj  ) = llcotf(iju,jpj-2)
-            END DO
-            DO ji = jpi/2, jpim1
-               iju = jpi - ji
-               llcotf(ji,jpjm1) = llcotf(iju,jpjm1)
-            END DO
-            DO ji = 1, jpi
-               ijt = jpi - ji + 1
-               llcotv(ji,jpj  ) = llcotv(ijt,jpjm1)
-            END DO
-            DO ji = jpi/2+1, jpi
-               ijt = jpi - ji + 1
-               llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)
-            END DO
-         ENDIF
-
-         ! Compute cartesian coordinates of coastline points
-         ! and the number of coastline points
-         icoast = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( llcotf(ji,jj) ) THEN
-                  icoast = icoast + 1
-                  zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )
-                  zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )
-                  zzc(icoast) = SIN( rad*gphif(ji,jj) )
-               ENDIF
-               IF( llcotu(ji,jj) ) THEN
-                  icoast = icoast+1
-                  zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )
-                  zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )
-                  zzc(icoast) = SIN( rad*gphiu(ji,jj) )
-               ENDIF
-               IF( llcotv(ji,jj) ) THEN
-                  icoast = icoast+1
-                  zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )
-                  zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )
-                  zzc(icoast) = SIN( rad*gphiv(ji,jj) )
-               ENDIF
-            END DO
-         END DO
-
-         ! Distance for the T-points
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( tmask(ji,jj,jk) == 0._wp ) THEN
-                  pdct(ji,jj,jk) = 0._wp
-               ELSE
-                  DO jl = 1, icoast
-                     zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   &
-                        &     + ( zyt(ji,jj) - zyc(jl) )**2   &
-                        &     + ( zzt(ji,jj) - zzc(jl) )**2
-                  END DO
-                  pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
-               ENDIF
-            END DO
-         END DO
-         !                                                ! ===============
-      END DO                                              !   End of slab
-      !                                                   ! ===============
-
-
-      ! 2. Create the  distance to the coast file in NetCDF format
-      ! ----------------------------------------------------------    
-      clname = 'dist.coast'
-      itime  = 0
-      CALL ymds2ju( 0     , 1       , 1     , 0._wp , zdate0 )
-      CALL restini( 'NONE', jpi     , jpj   , glamt, gphit ,   &
-         &          jpk   , gdept_1d, clname, itime, zdate0,   &
-         &          rdt   , icot                         )
-      CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
-      CALL restclo( icot )
-      !
-      CALL wrk_dealloc( jpi, jpj , zxt, zyt, zzt, zmask    )
-      CALL wrk_dealloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     )
-      DEALLOCATE( llcotu, llcotv, llcotf  )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('cofdis')
-      !
-   END SUBROUTINE cofdis
-   !!======================================================================
 END MODULE tradmp
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 5101)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 5102)
@@ -24,4 +24,5 @@
    USE trdtra
    USE trd_oce
+   USE iom
 
    IMPLICIT NONE
@@ -302,13 +303,8 @@
       !!----------------------------------------------------------------------
       !
+      INTEGER :: imask  !local file handle
+
       IF( nn_timing == 1 )  CALL timing_start('trc_dmp_init')
       !
-      SELECT CASE ( nn_hdmp_tr )
-      CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only'
-      CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp_tr, ' degrees'
-      CASE DEFAULT
-         WRITE(ctmp1,*) '          bad flag value for nn_hdmp_tr = ', nn_hdmp_tr
-         CALL ctl_stop(ctmp1)
-      END SELECT
 
       IF( lzoom )   nn_zdmp_tr = 0           ! restoring to climatology at closed north or south boundaries
@@ -325,9 +321,9 @@
          &   CALL ctl_stop( 'passive trace damping need key_tradmp to compute damping coef.' )
       !
-      !                          ! Damping coefficients initialization
-      IF( lzoom ) THEN   ;   CALL dtacof_zoom( restotr )
-      ELSE               ;   CALL dtacof( nn_hdmp_tr, rn_surf_tr, rn_bot_tr, rn_dep_tr,  &
-                             &            nn_file_tr, 'TRC'     , restotr                )
-      ENDIF
+      !                          ! Read damping coefficients from file
+      !Read in mask from file
+      CALL iom_open ( cn_resto_tr, imask)
+      CALL iom_get  ( imask, jpdom_autoglo, 'resto', restotr)
+      CALL iom_close( imask )
       !
       IF( nn_timing == 1 )  CALL timing_stop('trc_dmp_init')
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90	(revision 5101)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90	(revision 5102)
@@ -51,10 +51,6 @@
    !                                                 !!: ** newtonian damping namelist (nam_trcdmp) **
    !                          !!* Namelist namtrc_dmp : passive tracer newtonian damping *
-   INTEGER , PUBLIC ::   nn_hdmp_tr    ! = 0/-1/'latitude' for damping over passive tracer
    INTEGER , PUBLIC ::   nn_zdmp_tr    ! = 0/1/2 flag for damping in the mixed layer
-   REAL(wp), PUBLIC ::   rn_surf_tr    ! surface time scale for internal damping        [days]
-   REAL(wp), PUBLIC ::   rn_bot_tr     ! bottom time scale for internal damping         [days]
-   REAL(wp), PUBLIC ::   rn_dep_tr     ! depth of transition between rn_surf and rn_bot [meters]
-   INTEGER , PUBLIC ::   nn_file_tr    ! = 1 create a damping.coeff NetCDF file
+   CHARACTER(LEN=200) , PUBLIC :: cn_resto_tr    !File containing restoration coefficient
 
    !!----------------------------------------------------------------------
@@ -82,6 +78,5 @@
       NAMELIST/namtrc_zdf/ ln_trczdf_exp  , nn_trczdf_exp
       NAMELIST/namtrc_rad/ ln_trcrad
-      NAMELIST/namtrc_dmp/ nn_hdmp_tr, nn_zdmp_tr, rn_surf_tr, &
-        &                  rn_bot_tr , rn_dep_tr , nn_file_tr
+      NAMELIST/namtrc_dmp/ nn_zdmp_tr , cn_resto_tr
       !!----------------------------------------------------------------------
 
@@ -184,10 +179,6 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist namtrc_dmp : set damping parameter'
-         WRITE(numout,*) '      tracer damping option          nn_hdmp_tr = ', nn_hdmp_tr
          WRITE(numout,*) '      mixed layer damping option     nn_zdmp_tr = ', nn_zdmp_tr, '(zoom: forced to 0)'
-         WRITE(numout,*) '      surface time scale (days)      rn_surf_tr = ', rn_surf_tr
-         WRITE(numout,*) '      bottom time scale (days)       rn_bot_tr  = ', rn_bot_tr
-         WRITE(numout,*) '      depth of transition (meters)   rn_dep_tr  = ', rn_dep_tr
-         WRITE(numout,*) '      create a damping.coeff file    nn_file_tr = ', nn_file_tr
+         WRITE(numout,*) '      Restoration coeff file    cn_resto_tr = ', cn_resto_tr
       ENDIF
       !
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/README
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/README	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/README	(revision 5102)
@@ -0,0 +1,5 @@
+DMP_TOOLS should be used to create a netcdf file called resto.nc containing restoration coefficients for use with the tra_dmp module in NEMO. Further instructions for it's use are available in the NEMO users guide.
+
+The tool can be compiled using the maketools script in the NEMOGCM/TOOLS directory as follows:
+./maketools -m $ARCH -n DMP_TOOLS
+where $ARCH indicates the arch file to be used from the directory NEMOGCM/ARCH. For example to use NEMOGCM/ARCH/arch-PW7_MONSOON.fcm, $ARCH would be PW7_MONSOON.
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/namelist
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/namelist	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/namelist	(revision 5102)
@@ -0,0 +1,24 @@
+&nam_dmp_create
+    cp_cfg = 'orca'            ! Name of model grid (orca and C1D have special options - otherwise ignored)
+    cp_cfz = 'antarctic'       ! Name of zoom configuration (arctic and antarctic have some special treatment if lzoom=.true.)
+    jp_cfg = 2                 ! Resolution of the model (used for med_red_seas damping)
+    lzoom = .false.            ! Zoom configuration or not
+    ln_full_field = .false.    ! Calculate coefficient over whole of domain 
+    ln_med_red_seas = .true.   ! Damping in Med/Red Seas (or local modifications here if ln_full_field=.true.)
+    ln_old_31_lev_code = .true.   ! Replicate behaviour of old online code for 31 level model (Med/Red seas damping based on level number instead of depth)
+    ln_coast = .true.          ! Reduce near to coastlines 
+    ln_zero_top_layer = .true. ! No damping in top layer
+    ln_custom = .false.        ! Call "custom" module to apply user modifications to the damping coefficient field
+    nn_hdmp = 10               ! Damp poleward of this latitude (smooth transition up to maximum damping)
+    pn_surf = 0.25             ! Surface Relaxation timescale (days) 
+    pn_bot = 0.25              ! Bottom relaxation timescale (days)
+    pn_dep = 1000              ! Transition depth from upper to deep ocean
+    jperio = 2                 ! Lateral boundary condition (as specified in namelist_cfg for model run).
+/
+
+&nam_zoom_dmp
+    lzoom_n = .false.     ! Open boundary had northern edge?
+    lzoom_e = .false.     ! Open boundary at eastern edge?
+    lzoom_w = .false.     ! Open boundary at western edge?
+    lzoom_s = .false.     ! Open boundary at southern edge?
+/
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90	(revision 5102)
@@ -0,0 +1,220 @@
+MODULE coastdist
+
+   USE utils
+   USE netcdf
+
+   IMPLICIT NONE
+   PUBLIC 
+
+   CONTAINS
+
+   SUBROUTINE coast_dist_weight( presto ) 
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE coast_dist_weight ***
+      !!
+      !! ** Purpose: Weight restoration coefficient by distance to coast
+      !!
+      !! ** Method: 1) Calculate distance to coast
+      !!            2) Reduce resto with 1000km of coast
+      !!
+      IMPLICIT NONE
+      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: presto
+      REAL(wp), DIMENSION(jpi,jpj) :: zdct
+      REAL(wp) :: zinfl = 1000.e3_wp  ! Distance of influence of coast line (could be
+                                  ! a namelist setting)
+      INTEGER :: jj, ji           ! dummy loop indices
+      
+
+      CALL cofdis( zdct )
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zdct(ji,jj) = MIN( zinfl, zdct(ji,jj) )
+            presto(ji,jj) = presto(ji, jj) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj)/zinfl) )
+         END DO
+      END DO
+
+   END SUBROUTINE coast_dist_weight
+
+
+   SUBROUTINE cofdis( pdct )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE cofdis  ***
+      !!
+      !! ** Purpose :   Compute the distance between ocean T-points and the
+      !!      ocean model coastlines.
+      !!
+      !! ** Method  :   For each model level, the distance-to-coast is 
+      !!      computed as follows : 
+      !!       - The coastline is defined as the serie of U-,V-,F-points
+      !!      that are at the ocean-land bound.
+      !!       - For each ocean T-point, the distance-to-coast is then 
+      !!      computed as the smallest distance (on the sphere) between the 
+      !!      T-point and all the coastline points.
+      !!       - For land T-points, the distance-to-coast is set to zero.
+      !!
+      !! ** Action  : - pdct, distance to the coastline (argument)
+      !!              - NetCDF file 'dist.coast.nc' 
+      !!----------------------------------------------------------------------
+      !!
+      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) ::   pdct   ! distance to the coastline
+      !!
+      INTEGER ::   ji, jj, jl   ! dummy loop indices
+      INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers
+      CHARACTER (len=32) ::   clname                     ! local name
+      REAL(wp) ::   zdate0                               ! local scalar
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask
+      REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace
+      LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace
+
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( zxt(jpi,jpj) , zyt(jpi,jpj) , zzt(jpi,jpj) , zmask(jpi,jpj)    )
+      ALLOCATE(zxc(3*jpi*jpj), zyc(3*jpi*jpj), zzc(3*jpi*jpj), zdis(3*jpi*jpj)     )
+      ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj)  )
+      ALLOCATE( gphiu(jpi,jpj), gphiv(jpi,jpj), gphif(jpi,jpj)  )
+      ALLOCATE( glamu(jpi,jpj), glamv(jpi,jpj), glamf(jpi,jpj), glamt(jpi,jpj)  )
+      ALLOCATE( umask(jpi,jpj), vmask(jpi,jpj), fmask(jpi,jpj)  )
+      !
+
+      CALL check_nf90( nf90_get_var( ncin, gphit_id, gphit, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, gphiu_id, gphiu, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, gphiv_id, gphiv, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, gphif_id, gphif, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, glamt_id, glamt, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, glamu_id, glamu, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, glamv_id, glamv, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, glamf_id, glamf, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, tmask_id, tmask, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, umask_id, umask, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, vmask_id, vmask, (/ 1,1 /), (/ jpi, jpj /) ) )
+      CALL check_nf90( nf90_get_var( ncin, fmask_id, fmask, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+      pdct(:,:) = 0._wp
+      zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )
+      zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )
+      zzt(:,:) = SIN( rad * gphit(:,:) )
+
+
+         ! Define the coastline points (U, V and F)
+         DO jj = 2, jpj-1
+            DO ji = 2, jpi-1
+               zmask(ji,jj) =  ( tmask(ji,jj+1) + tmask(ji+1,jj+1) &
+                   &           + tmask(ji,jj  ) + tmask(ji+1,jj  ) )
+               llcotu(ji,jj) = ( tmask(ji,jj ) + tmask(ji+1,jj  ) == 1._wp ) 
+               llcotv(ji,jj) = ( tmask(ji,jj  ) + tmask(ji  ,jj+1) == 1._wp ) 
+               llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp )
+            END DO
+         END DO
+
+         ! Lateral boundaries conditions
+         llcotu(:, 1 ) = umask(:,  2  ) == 1
+         llcotu(:,jpj) = umask(:,jpj-1) == 1
+         llcotv(:, 1 ) = vmask(:,  2  ) == 1
+         llcotv(:,jpj) = vmask(:,jpj-1) == 1
+         llcotf(:, 1 ) = fmask(:,  2  ) == 1
+         llcotf(:,jpj) = fmask(:,jpj-1) == 1
+
+         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
+            llcotu( 1 ,:) = llcotu(jpi-1,:)
+            llcotu(jpi,:) = llcotu(  2  ,:)
+            llcotv( 1 ,:) = llcotv(jpi-1,:)
+            llcotv(jpi,:) = llcotv(  2  ,:)
+            llcotf( 1 ,:) = llcotf(jpi-1,:)
+            llcotf(jpi,:) = llcotf(  2  ,:)
+         ELSE
+            llcotu( 1 ,:) = umask(  2  ,:) == 1
+            llcotu(jpi,:) = umask(jpi-1,:) == 1
+            llcotv( 1 ,:) = vmask(  2  ,:) == 1
+            llcotv(jpi,:) = vmask(jpi-1,:) == 1
+            llcotf( 1 ,:) = fmask(  2  ,:) == 1
+            llcotf(jpi,:) = fmask(jpi-1,:) == 1
+         ENDIF
+         IF( jperio == 3 .OR. jperio == 4 ) THEN
+            DO ji = 1, jpi-1
+               iju = jpi - ji + 1
+               llcotu(ji,jpj  ) = llcotu(iju,jpj-2)
+               llcotf(ji,jpj-1) = llcotf(iju,jpj-2)
+               llcotf(ji,jpj  ) = llcotf(iju,jpj-3)
+            END DO
+            DO ji = jpi/2, jpi-1
+               iju = jpi - ji + 1
+               llcotu(ji,jpj-1) = llcotu(iju,jpj-1)
+            END DO
+            DO ji = 2, jpi
+               ijt = jpi - ji + 2
+               llcotv(ji,jpj-1) = llcotv(ijt,jpj-2)
+               llcotv(ji,jpj  ) = llcotv(ijt,jpj-3)
+            END DO
+         ENDIF
+         IF( jperio == 5 .OR. jperio == 6 ) THEN
+            DO ji = 1, jpi-1
+               iju = jpi - ji
+               llcotu(ji,jpj  ) = llcotu(iju,jpj-1)
+               llcotf(ji,jpj  ) = llcotf(iju,jpj-2)
+            END DO
+            DO ji = jpi/2, jpi-1
+               iju = jpi - ji
+               llcotf(ji,jpj-1) = llcotf(iju,jpj-1)
+            END DO
+            DO ji = 1, jpi
+               ijt = jpi - ji + 1
+               llcotv(ji,jpj  ) = llcotv(ijt,jpj-1)
+            END DO
+            DO ji = jpi/2+1, jpi
+               ijt = jpi - ji + 1
+               llcotv(ji,jpj-1) = llcotv(ijt,jpj-1)
+            END DO
+         ENDIF
+
+         ! Compute cartesian coordinates of coastline points
+         ! and the number of coastline points
+         icoast = 0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( llcotf(ji,jj) ) THEN
+                  icoast = icoast + 1
+                  zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )
+                  zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )
+                  zzc(icoast) = SIN( rad*gphif(ji,jj) )
+               ENDIF
+               IF( llcotu(ji,jj) ) THEN
+                  icoast = icoast+1
+                  zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )
+                  zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )
+                  zzc(icoast) = SIN( rad*gphiu(ji,jj) )
+               ENDIF
+               IF( llcotv(ji,jj) ) THEN
+                  icoast = icoast+1
+                  zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )
+                  zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )
+                  zzc(icoast) = SIN( rad*gphiv(ji,jj) )
+               ENDIF
+            END DO
+         END DO
+
+         ! Distance for the T-points
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj) == 0._wp ) THEN
+                  pdct(ji,jj) = 0._wp
+               ELSE
+                  DO jl = 1, icoast
+                     zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   &
+                        &     + ( zyt(ji,jj) - zyc(jl) )**2   &
+                        &     + ( zzt(ji,jj) - zzc(jl) )**2
+                  END DO
+                  pdct(ji,jj) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
+               ENDIF
+            END DO
+         END DO
+
+      DEALLOCATE( zxt , zyt , zzt , zmask    )
+      DEALLOCATE(zxc, zyc, zzc, zdis    )
+      DEALLOCATE( llcotu, llcotv, llcotf  )
+      DEALLOCATE( gphiu, gphiv, gphif  )
+      DEALLOCATE( glamu, glamv, glamf, glamt  )
+      DEALLOCATE( umask, vmask, fmask  )
+
+   END SUBROUTINE cofdis
+
+END MODULE coastdist
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/custom.F90
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/custom.F90	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/custom.F90	(revision 5102)
@@ -0,0 +1,22 @@
+MODULE custom
+
+   USE utils
+
+   IMPLICIT NONE
+   PUBLIC
+
+   CONTAINS
+
+   SUBROUTINE custom_resto( presto )
+      !!---------------------------------
+      !!     **ROUTINE: custom_resto
+      !!
+      !!  ** Purpose: Module to be edited by users to create custom restoration
+      !!              coefficient files (e.g. regional damping).
+      !!
+      !!-------------------------------------
+      REAL(wp), DIMENSION(jpi,jpk), INTENT(inout) :: presto
+
+   END SUBROUTINE custom_resto
+
+END MODULE custom
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90	(revision 5102)
@@ -0,0 +1,123 @@
+PROGRAM make_dmp_file
+  !================================================================================
+  !               *** PROGRAM make_dmp_file ****
+  !================================================================================
+  !
+  !  Purpose: Create a file containing a spacially varying
+  !           restoration coefficient to be used by TRADMP
+  !
+  !  Method:  1) Read in tmask from mesh_mask file to use as a template
+  !           2) Calculate restoration coefficients according to options
+  !           specified in the namelist. The user may modify custom.F90 to
+  !           specify specific damping options e.g. to mask certain regions only).
+  !           3) Write the array to output file
+  !
+  !  History: Original code: Tim Graham (Jul 2014) - some code moved from
+  !                            old tradmp.F90 module to this tool (as part of NEMO
+  !                            simplification process).
+  !-------------------------------------------------------------------------------
+
+  ! Declare variables
+  USE netcdf 
+  USE utils
+  USE coastdist
+  USE med_red_seas
+  USE zoom
+  USE custom
+
+  IMPLICIT NONE
+  INTEGER  :: ji, jj, jk                         ! dummpy loop variables
+  REAL(wp) :: zsdmp, zbdmp                     ! Surface and bottom damping coeff
+  CHARACTER(LEN=200) :: meshfile = 'mesh_mask.nc'   ! mesh file
+  CHARACTER(LEN=200) :: outfile = 'resto.nc'     ! output file
+  REAL(wp) :: zlat, zlat2, zlat0
+
+  ! Read namelist
+  OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
+  READ( numnam, nam_dmp_create )
+  CLOSE( numnam )
+  
+  IF ( ln_full_field .AND. lzoom ) THEN
+    WRITE(numerr,*) 'Only one of ln_full_field and lzoom can be .true.'
+    STOP
+  ENDIF
+
+  CALL grid_info(meshfile)
+  WRITE(numout, *) 'jpi = ',jpi
+  WRITE(numout, *) 'jpj = ',jpj
+  WRITE(numout, *) 'jpk = ',jpk
+
+  ALLOCATE( resto(jpi, jpj) )
+
+  !Create output file
+  CALL make_outfile( outfile )
+  
+  CALL check_nf90( nf90_get_var( ncin, gphit_id, gphit, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+  !Calculate surface and bottom damping coefficients
+  zsdmp = 1._wp / ( pn_surf * rday )
+  zbdmp = 1._wp / ( pn_bot  * rday )
+
+  !Loop through levels and read in tmask for each level as starting point for
+  !coefficient array
+  DO jk = 1, jpk-1
+     resto(:,:) = 0._wp
+     
+     IF (.NOT. (jk == 1 .AND. ln_zero_top_layer) ) THEN 
+        !Read in tmask depth for this level
+        CALL check_nf90( nf90_get_var( ncin, tmask_id, tmask, (/ 1,1,jk /), (/ jpi, jpj,1 /) ) )
+        CALL check_nf90( nf90_get_var( ncin, gdept_id, gdept, (/ 1,1,jk /), (/ jpi, jpj,1 /) ) )
+     
+
+        IF ( ln_full_field ) THEN
+           !Set basic value of resto
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 resto(ji,jj) = tmask(ji, jj) * (zbdmp + (zsdmp-zbdmp) * EXP(-gdept(ji,jj)/pn_dep))
+              END DO
+           END DO
+           IF ((nn_hdmp > 0)) THEN
+              zlat0 = 10. !width of latitude strip where resto decreases
+              zlat2 = nn_hdmp + zlat0
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    zlat = ABS(gphit(ji,jj))
+                    IF ( nn_hdmp <= zlat .AND. zlat <= zlat2 ) THEN
+                       resto(ji,jj) = resto(ji,jj) * 0.5_wp * (  1._wp - COS( rpi*(zlat-nn_hdmp)/zlat0 ) )
+                    ELSE IF ( zlat < nn_hdmp ) THEN
+                       resto(ji,jj) = 0._wp
+                    ENDIF
+                 END DO
+              END DO
+           ENDIF
+    
+           IF (ln_coast) THEN
+              ! Reduce damping in vicinity of coastlines
+              CALL coast_dist_weight(resto)
+           ENDIF
+        ENDIF
+
+        ! Damping in Med/Red Seas (or local modifications if full field is set)
+        IF (ln_med_red_seas .AND. (cp_cfg == 'orca') .AND. (.NOT. lzoom)) THEN
+           CALL med_red_dmp(resto, jk, ln_old_31_lev_code)
+        ENDIF
+
+        IF ( lzoom ) THEN
+              CALL dtacof_zoom(resto, tmask)
+        ENDIF
+        
+        !Any user modifications can be added in the custom module
+        IF ( ln_custom ) THEN
+              CALL custom_resto( resto )
+        ENDIF
+     ENDIF
+
+     ! Write out resto for this level
+     CALL check_nf90( nf90_put_var( ncout, resto_id, resto, (/ 1,1,jk /), (/ jpi, jpj,1 /) ) )
+
+  END DO
+  
+  ! Close the output file
+  CALL check_nf90( nf90_close(ncout) )
+
+END PROGRAM make_dmp_file
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90	(revision 5102)
@@ -0,0 +1,135 @@
+MODULE med_red_seas
+
+   USE utils
+   
+   IMPLICIT NONE
+   PUBLIC
+
+   CONTAINS  
+
+   SUBROUTINE med_red_dmp(presto, jk, ln_31_lev)
+      !!------------------------------------
+      !!    **ROUTINE: med_red_dmp
+      !! 
+      !! **Purpose: Apply specific modifications to damping coefficients on ORCA
+      !!            grids in Med and Red Seas
+      !!
+      !!-----------------------------------
+      INTEGER :: ij0,ij1,ii0,ii1,ji,jj     
+      INTEGER, INTENT(in) :: jk
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmrs
+      REAL(wp) :: zhfac, zsdmp, zbdmp
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: presto
+      LOGICAL, INTENT(in), OPTIONAL :: ln_31_lev
+      LOGICAL :: l_31_lev
+
+      WRITE(numout,*) 'ORCA Med and Red Seas Damping'
+      
+      IF ( PRESENT(ln_31_lev)) THEN
+         l_31_lev = ln_31_lev
+      ELSE
+         l_31_lev = .false.
+      ENDIF
+      
+      ALLOCATE( zmrs(jpi, jpj) )
+         !
+         zmrs(:,:) = 0._wp
+         !
+         SELECT CASE ( jp_cfg )
+         !                                           ! =======================
+         CASE ( 4 )                                  !  ORCA_R4 configuration 
+            !                                        ! =======================
+            ij0 =  50   ;   ij1 =  56                    ! Mediterranean Sea
+
+            ii0 =  81   ;   ii1 =  91   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1.
+            ij0 =  50   ;   ij1 =  55
+            ii0 =  75   ;   ii1 =  80   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1.
+            ij0 =  52   ;   ij1 =  53
+            ii0 =  70   ;   ii1 =  74   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1.
+           !
+           !                                        ! =======================
+         CASE ( 2 )                                  !  ORCA_R2 configuration 
+            !                                        ! =======================
+            ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea
+            ii0 = 157   ;   ii1 = 181   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            ij0 = 100   ;   ij1 = 110
+            ii0 = 144   ;   ii1 = 156   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            ij0 = 100   ;   ij1 = 103
+            ii0 = 139   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            !
+            ij0 = 101   ;   ij1 = 102                    ! Decrease before Gibraltar Strait
+            ii0 = 139   ;   ii1 = 141   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0._wp
+            ii0 = 142   ;   ii1 = 142   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp / 90._wp
+            ii0 = 143   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40_wp
+            ii0 = 144   ;   ii1 = 144   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75_wp
+            !
+            ij0 =  87   ;   ij1 =  96                    ! Red Sea
+            ii0 = 147   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            !
+            ij0 =  91   ;   ij1 =  91                    ! Decrease before Bab el Mandeb Strait
+            ii0 = 153   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.80_wp
+            ij0 =  90   ;   ij1 =  90
+            ii0 = 153   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40_wp
+            ij0 =  89   ;   ij1 =  89
+            ii0 = 158   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp / 90._wp
+            ij0 =  88   ;   ij1 =  88
+            ii0 = 160   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0._wp
+            !
+            !                                        ! =======================
+         CASE ( 05 )                                 !  ORCA_R05 configuration
+            !                                        ! =======================
+            ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea
+            ij0 = 324   ;   ij1 = 333   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            ii0 = 575   ;   ii1 = 658
+            ij0 = 314   ;   ij1 = 366   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            !
+            ii0 = 641   ;   ii1 = 651                    ! Black Sea (remaining part
+            ij0 = 367   ;   ij1 = 372   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            !
+            ij0 = 324   ;   ij1 = 333                    ! Decrease before Gibraltar Strait
+            ii0 = 565   ;   ii1 = 565   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp / 90._wp
+            ii0 = 566   ;   ii1 = 566   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40_wp
+            ii0 = 567   ;   ii1 = 567   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75_wp
+            !
+            ii0 = 641   ;   ii1 = 665                    ! Red Sea
+            ij0 = 270   ;   ij1 = 310   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp
+            !
+            ii0 = 666   ;   ii1 = 675                    ! Decrease before Bab el Mandeb Strait
+            ij0 = 270   ;   ij1 = 290   
+            DO ji = ii0, ii1
+               zmrs( ji , ij0:ij1 ) = 0.1_wp * ABS( FLOAT(ji - ii1) )
+            END DO 
+            !                                       ! ========================
+         CASE ( 025 )                               !  ORCA_R025 configuration 
+            !                                       ! ========================
+            WRITE(numerr,*) ' Mediterranean and Red Sea damping option not implemented for ORCA_R025'
+            WRITE(numerr,*) ' Set ln_med_red = .false.'
+            STOP
+            !
+         END SELECT
+
+         zsdmp = 1._wp / ( pn_surf * rday )
+         zbdmp = 1._wp / ( pn_bot  * rday )
+
+         ! The l_31_lev option is used to reproduce the old behaviour of
+         ! defining the restoration coefficient based on the level number.
+         ! This is included to allow damping coefficients for reference
+         ! configurations to be kept the same.
+         IF (l_31_lev) THEN
+            IF (jk <= 17) THEN
+               zhfac = 0.5_wp * (  1. - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday
+            ELSE
+               zhfac = 1._wp / rday
+            ENDIF
+         ELSE
+            zhfac = (  zbdmp + (zsdmp-zbdmp) * EXP( -gdept(1,1)/pn_dep )  )
+         ENDIF
+
+         presto(:,:) = zmrs(:,:) * zhfac + ( 1._wp - zmrs(:,:) ) * presto(:,:)
+
+         DEALLOCATE( zmrs )         
+
+   END SUBROUTINE med_red_dmp
+
+
+END MODULE med_red_seas
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90	(revision 5102)
@@ -0,0 +1,130 @@
+MODULE utils
+
+  USE netcdf
+
+  IMPLICIT NONE
+  PUBLIC
+  
+  INTEGER, PUBLIC, PARAMETER   :: dp=8 , sp=4, wp=dp
+  INTEGER :: tmask_id, umask_id, vmask_id, fmask_id
+  INTEGER :: gdept_id
+  INTEGER :: gphit_id, gphiv_id, gphiu_id, gphif_id        ! Variable ids
+  INTEGER :: glamt_id, glamv_id, glamu_id, glamf_id        ! Variable ids
+  INTEGER :: resto_id                           ! Variable ID for output
+  INTEGER :: jpi, jpj, jpk                      ! Size of domain
+  INTEGER  :: ncin, ncout                              ! File handles for netCDF files
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphit, glamt
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphiu, glamu
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphiv, glamv
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphif, glamf
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: tmask, umask, vmask, fmask
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gdept
+  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: resto
+
+  INTEGER,PARAMETER :: numout = 6
+  INTEGER,PARAMETER :: numerr = 0
+  INTEGER,PARAMETER :: numnam = 11
+  REAL(wp),PARAMETER :: rday = 86400           ! seconds in a day
+  REAL(wp),PARAMETER ::   rpi = 3.141592653589793
+  REAL(wp),PARAMETER ::   rad = 3.141592653589793/180.
+  REAL(wp),PARAMETER ::   ra =  6371229.
+
+  ! Namelist variables
+  CHARACTER(LEN=30) :: cp_cfg = 'ORCA'
+  CHARACTER(LEN=30) :: cp_cfz = 'No zoom'
+  INTEGER :: jp_cfg = 2
+  REAL(KIND=8) :: pn_surf = 1
+  REAL(KIND=8) :: pn_bot = 1
+  REAL(KIND=8) :: pn_dep = 1000         
+  INTEGER  :: nn_hdmp = 0                          ! damping option
+  INTEGER  :: jperio = 0                          ! damping option
+  LOGICAL :: lzoom = .false.
+  LOGICAL :: ln_coast = .false.
+  LOGICAL :: ln_full_field = .true.
+  LOGICAL :: ln_med_red_seas = .false.
+  LOGICAL :: ln_old_31_lev_code = .false.
+  LOGICAL :: ln_zero_top_layer = .false.
+  LOGICAL :: ln_custom = .false.
+
+  NAMELIST/nam_dmp_create/cp_cfg, cp_cfz, jp_cfg, lzoom, ln_full_field, &
+                          ln_med_red_seas, ln_old_31_lev_code, ln_coast, &
+                          ln_zero_top_layer, ln_custom, &
+                          pn_surf, pn_bot, pn_dep, nn_hdmp, jperio
+
+  CONTAINS 
+
+  SUBROUTINE grid_info(mesh)
+     CHARACTER(LEN=*),INTENT(in) :: mesh
+
+     ! Open meshfile
+     CALL check_nf90( nf90_open(mesh, NF90_NOWRITE, ncin), 'Error opening mesh_mask file' )
+  
+     ! Get size of grid from meshfile
+     CALL dimlen( ncin, 'x', jpi )
+     CALL dimlen( ncin, 'y', jpj )
+     CALL dimlen( ncin, 'z', jpk )
+
+     ALLOCATE( tmask(jpi, jpj), gdept(jpi, jpj), gphit(jpi,jpj) )
+
+     !Get ID of tmask in meshfile
+     CALL check_nf90( nf90_inq_varid( ncin, 'tmask', tmask_id ), 'Cannot get variable ID for tmask')
+     CALL check_nf90( nf90_inq_varid( ncin, 'umask', umask_id ), 'Cannot get variable ID for umask')
+     CALL check_nf90( nf90_inq_varid( ncin, 'vmask', vmask_id ), 'Cannot get variable ID for vmask')
+     CALL check_nf90( nf90_inq_varid( ncin, 'fmask', fmask_id ), 'Cannot get variable ID for fmask')
+     CALL check_nf90( nf90_inq_varid( ncin, 'gdept_0', gdept_id ), 'Cannot get variable ID for gdept_0')
+     CALL check_nf90( nf90_inq_varid( ncin, 'gphit', gphit_id ), 'Cannot get variable ID for gphit')
+     CALL check_nf90( nf90_inq_varid( ncin, 'gphiu', gphiu_id ), 'Cannot get variable ID for gphiu')
+     CALL check_nf90( nf90_inq_varid( ncin, 'gphiv', gphiv_id ), 'Cannot get variable ID for gphiv')
+     CALL check_nf90( nf90_inq_varid( ncin, 'gphif', gphif_id ), 'Cannot get variable ID for gphif')
+     CALL check_nf90( nf90_inq_varid( ncin, 'glamt', glamt_id ), 'Cannot get variable ID for glamt')
+     CALL check_nf90( nf90_inq_varid( ncin, 'glamu', glamu_id ), 'Cannot get variable ID for glamu')
+     CALL check_nf90( nf90_inq_varid( ncin, 'glamv', glamv_id ), 'Cannot get variable ID for glamv')
+     CALL check_nf90( nf90_inq_varid( ncin, 'glamf', glamf_id ), 'Cannot get variable ID for glamf')
+  
+  END SUBROUTINE grid_info
+
+  SUBROUTINE dimlen( ncid, dimname, len )
+     ! Determine the length of dimension dimname
+     INTEGER, INTENT(in)          :: ncid
+     CHARACTER(LEN=*), INTENT(in) :: dimname
+     INTEGER, INTENT(out)         :: len
+     ! Local variables
+     INTEGER :: id_var, istatus
+
+     id_var = 1
+     CALL check_nf90( nf90_inq_dimid(ncid, dimname, id_var), 'Dimension not found in file')
+     CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len))
+
+  END SUBROUTINE dimlen
+
+  SUBROUTINE make_outfile( filename )
+     ! Create the output file
+     ! Define dimensions and resto variable
+     CHARACTER(LEN=*), INTENT(in) :: filename
+     INTEGER :: id_x, id_y, id_z
+
+     CALL check_nf90( nf90_create(filename, NF90_CLOBBER, ncout), 'Could not create output file')
+     CALL check_nf90( nf90_def_dim(ncout, 'x', jpi, id_x) )
+     CALL check_nf90( nf90_def_dim(ncout, 'y', jpj, id_y) )
+     CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) )
+
+     CALL check_nf90( nf90_def_var(ncout, 'resto', nf90_double, (/id_x,id_y,id_z/), resto_id ) )
+     CALL check_nf90( nf90_enddef(ncout) )
+
+  END SUBROUTINE make_outfile
+
+
+  SUBROUTINE check_nf90( istat, message )
+     !Check for netcdf errors
+     INTEGER, INTENT(in) :: istat
+     CHARACTER(LEN=*), INTENT(in), OPTIONAL :: message
+     
+     IF (istat /= nf90_noerr) THEN
+        WRITE(numerr,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
+        IF ( PRESENT(message) ) THEN ; WRITE(numerr,*) message ; ENDIF
+        STOP
+     ENDIF
+
+  END SUBROUTINE check_nf90
+
+END MODULE utils
Index: /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90
===================================================================
--- /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90	(revision 5102)
+++ /trunk/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90	(revision 5102)
@@ -0,0 +1,90 @@
+MODULE zoom
+   
+   USE utils
+
+   CONTAINS
+
+   SUBROUTINE dtacof_zoom( presto, mask)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dtacof_zoom  ***
+      !!
+      !! ** Purpose :   Compute the damping coefficient for zoom domain
+      !!
+      !! ** Method  : - set along closed boundary due to zoom a damping over
+      !!                6 points with a max time scale of 5 days.
+      !!              - ORCA arctic/antarctic zoom: set the damping along
+      !!                south/north boundary over a latitude strip.
+      !!
+      !! ** Action  : - resto, the damping coeff. for T and S
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  ::   presto   ! restoring coeff. (s-1)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::   mask   ! restoring coeff. (s-1)
+      !
+      INTEGER  ::   ji, jj, jn   ! dummy loop indices
+      REAL(wp) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar
+      REAL(wp), DIMENSION(6)  ::   zfact               ! 1Dworkspace
+
+      !Namelist variables
+      LOGICAL :: lzoom_w, lzoom_e, lzoom_n, lzoom_s 
+      NAMELIST/nam_zoom_dmp/lzoom_n,lzoom_e,lzoom_w,lzoom_s
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start( 'dtacof_zoom')
+      !
+      
+      ! Read namelist
+      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
+      READ( numnam, nam_dmp_create )
+      CLOSE( numnam )
+
+      zfact(1) =  1._wp
+      zfact(2) =  1._wp
+      zfact(3) = 11._wp / 12._wp
+      zfact(4) =  8._wp / 12._wp
+      zfact(5) =  4._wp / 12._wp
+      zfact(6) =  1._wp / 12._wp
+      zfact(:) = zfact(:) / ( 5._wp * rday )    ! 5 days max restoring time scale
+
+      presto(:,:) = 0._wp
+
+      ! damping along the forced closed boundary over 6 grid-points
+      DO jn = 1, 6
+         IF( lzoom_w )   presto( jn, : )                    = zfact(jn)   ! west  closed
+         IF( lzoom_s )   presto( : , jn )                    = zfact(jn)   ! south closed 
+         IF( lzoom_e )   presto( jpi+1-jn , : ) = zfact(jn)   ! east  closed 
+         IF( lzoom_n )   presto( : , jpj+1-jn ) = zfact(jn)   ! north closed
+      END DO
+
+      !                                           ! ====================================================
+      IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN   !  ORCA configuration : arctic or antarctic zoom
+         !                                        ! ====================================================
+         WRITE(numout,*)
+         IF(cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom'
+         IF(cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom'
+         WRITE(numout,*)
+         !
+         !                          ! Initialization : 
+         presto(:,:) = 0._wp
+         zlat0 = 10._wp                     ! zlat0 : latitude strip where resto decreases
+         zlat1 = 30._wp                     ! zlat1 : resto = 1 before zlat1
+         zlat2 = zlat1 + zlat0              ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
+         z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zlat = ABS( gphit(ji,jj) )
+               IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
+                  presto(ji,jj) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  ) 
+               ELSEIF( zlat < zlat1 ) THEN
+                  presto(ji,jj) = z1_5d
+               ENDIF
+            END DO
+         END DO
+         !
+      ENDIF
+      !                             ! Mask resto array
+      presto(:,:) = presto(:,:) * mask(:,:)
+
+   END SUBROUTINE dtacof_zoom
+
+END MODULE zoom
