Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/CONFIG/SHARED/namelist_top_ref
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/CONFIG/SHARED/namelist_top_ref	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/CONFIG/SHARED/namelist_top_ref	(revision 4739)
@@ -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  = 'resto_tr.nc'    !  create a damping.coeff NetCDF file (=1) or not (=0)
 /
 !-----------------------------------------------------------------------
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90	(revision 4739)
@@ -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
    !!----------------------------------------------------------------------
 
@@ -27,4 +30,5 @@
    USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
+   USE iom            ! I/O manager
 
    IMPLICIT NONE
@@ -75,4 +79,5 @@
       NAMELIST/namc1d_dyndmp/ ln_dyndmp
       INTEGER :: ios
+      INTEGER :: imask
       !!----------------------------------------------------------------------
 
@@ -93,10 +98,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
@@ -106,14 +108,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'
@@ -132,13 +124,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: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 4739)
@@ -36,4 +36,5 @@
    USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -46,7 +47,8 @@
    LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag
    INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer
-   CHARACTER(LEN=200) :: cn_resto      ! name of netcdf file containing restoration coefficient field
+   CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field
    LOGICAL , PUBLIC ::   ln_miss      ! check for missing data in T/S data file (slow?)
    REAL(wp), PUBLIC ::   rn_miss      ! Value of missing data
+   !
 
 
@@ -189,7 +191,7 @@
       !!----------------------------------------------------------------------
 
-      NAMELIST/namtra_dmp/ ln_tradmp, rn_dmp, nn_zdmp, cn_resto
+      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto
       INTEGER ::  ios         ! Local integer for output status of namelist read
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dmp_mask   ! 3D mask
+      INTEGER :: imask        ! File handle 
       !!----------------------------------------------------------------------
 
@@ -224,4 +226,5 @@
          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
 
          !Initialisation of dtatsd - Would it be better to have dmpdta routine
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 4739)
@@ -302,13 +302,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 +320,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: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90	(revision 4739)
@@ -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 , nn_file_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: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/namelist
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/namelist	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/namelist	(revision 4739)
@@ -3,7 +3,8 @@
     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 = .true.             ! Zoom configuration or not
+    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
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90	(revision 4739)
@@ -19,7 +19,7 @@
       !!
       IMPLICIT NONE
-      REAL(8), DIMENSION(jpi,jpj), INTENT( inout ) :: presto
-      REAL(8), DIMENSION(jpi,jpj) :: zdct
-      REAL(8) :: zinfl = 1000.e3  ! Distance of influence of coast line (could be
+      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
@@ -30,5 +30,5 @@
          DO ji = 1, jpi
             zdct(ji,jj) = MIN( zinfl, zdct(ji,jj) )
-            presto(ji,jj) = presto(ji, jj) * 0.5 * (  1. - COS( rpi*zdct(ji,jj)/zinfl) )
+            presto(ji,jj) = presto(ji, jj) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj)/zinfl) )
          END DO
       END DO
@@ -57,12 +57,12 @@
       !!----------------------------------------------------------------------
       !!
-      REAL(8), DIMENSION(jpi,jpj), INTENT( out ) ::   pdct   ! distance to the coastline
+      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(8) ::   zdate0                               ! local scalar
-      REAL(8), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask
-      REAL(8), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace
+      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
 
@@ -90,5 +90,5 @@
       CALL check_nf90( nf90_get_var( ncin, fmask_id, fmask, (/ 1,1 /), (/ jpi, jpj /) ) )
 
-      pdct(:,:) = 0.
+      pdct(:,:) = 0._wp
       zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )
       zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )
@@ -101,7 +101,7 @@
                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. ) 
-               llcotv(ji,jj) = ( tmask(ji,jj  ) + tmask(ji  ,jj+1) == 1. ) 
-               llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )
+               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
@@ -196,6 +196,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( tmask(ji,jj) == 0. ) THEN
-                  pdct(ji,jj) = 0.
+               IF( tmask(ji,jj) == 0._wp ) THEN
+                  pdct(ji,jj) = 0._wp
                ELSE
                   DO jl = 1, icoast
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/custom.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/custom.F90	(revision 4739)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/custom.F90	(revision 4739)
@@ -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: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90	(revision 4739)
@@ -24,11 +24,12 @@
   USE med_red_seas
   USE zoom
+  USE custom
 
   IMPLICIT NONE
   INTEGER  :: ji, jj, jk                         ! dummpy loop variables
-  REAL(8) :: zsdmp, zbdmp                     ! Surface and bottom damping coeff
+  REAL(wp) :: zsdmp, zbdmp                     ! Surface and bottom damping coeff
   CHARACTER(LEN=200) :: meshfile = 'mesh_mask.nc'   ! mesh file
-  CHARACTER(LEN=200) :: outfile = 'dmp_mask.nc'     ! output file
-  REAL(8) :: zlat, zlat2, zlat0
+  CHARACTER(LEN=200) :: outfile = 'resto.nc'     ! output file
+  REAL(wp) :: zlat, zlat2, zlat0
 
   ! Read namelist
@@ -55,11 +56,11 @@
 
   !Calculate surface and bottom damping coefficients
-  zsdmp = 1. / ( pn_surf * rday )
-  zbdmp = 1. / ( pn_bot  * rday )
+  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.
+     resto(:,:) = 0._wp
      
      IF (.NOT. (jk == 1 .AND. ln_zero_top_layer) ) THEN 
@@ -83,7 +84,7 @@
                     zlat = ABS(gphit(ji,jj))
                     IF ( nn_hdmp <= zlat .AND. zlat <= zlat2 ) THEN
-                       resto(ji,jj) = resto(ji,jj) * 0.5 * (  1. - COS( rpi*(zlat-nn_hdmp)/zlat0 ) )
+                       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.
+                       resto(ji,jj) = 0._wp
                     ENDIF
                  END DO
@@ -98,6 +99,6 @@
 
         ! Damping in Med/Red Seas (or local modifications if full field is set)
-        IF (ln_med_red_seas .AND. (cp_cfg == 'orca')) THEN
-           CALL med_red_dmp(resto)
+        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
 
@@ -105,5 +106,9 @@
               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
 
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90	(revision 4739)
@@ -8,5 +8,5 @@
    CONTAINS  
 
-   SUBROUTINE med_red_dmp(presto)
+   SUBROUTINE med_red_dmp(presto, jk, ln_31_lev)
       !!------------------------------------
       !!    **ROUTINE: med_red_dmp
@@ -17,13 +17,22 @@
       !!-----------------------------------
       INTEGER :: ij0,ij1,ii0,ii1,ji,jj     
-      REAL(8), DIMENSION(:,:), ALLOCATABLE :: zmrs
-      REAL(8) :: zhfac, zsdmp, zbdmp
-      REAL(8), DIMENSION(jpi,jpj), INTENT(inout) :: presto
+      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.
+         zmrs(:,:) = 0._wp
          !
          SELECT CASE ( jp_cfg )
@@ -43,27 +52,27 @@
             !                                        ! =======================
             ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea
-            ii0 = 157   ;   ii1 = 181   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1.
+            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.
+            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.
+            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.
-            ii0 = 142   ;   ii1 = 142   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. / 90.
-            ii0 = 143   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40
-            ii0 = 144   ;   ii1 = 144   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75
+            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.
+            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
+            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
+            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. / 90.
+            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.
+            ii0 = 160   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0._wp
             !
             !                                        ! =======================
@@ -71,23 +80,23 @@
             !                                        ! =======================
             ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea
-            ij0 = 324   ;   ij1 = 333   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1.
+            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.
+            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.
+            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. / 90.
-            ii0 = 566   ;   ii1 = 566   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40
-            ii0 = 567   ;   ii1 = 567   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75
+            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.
+            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 * ABS( FLOAT(ji - ii1) )
+               zmrs( ji , ij0:ij1 ) = 0.1_wp * ABS( FLOAT(ji - ii1) )
             END DO 
             !                                       ! ========================
@@ -100,13 +109,22 @@
          END SELECT
 
-         ! Note that the original "online" code had a dependency on model levels
-         ! here (as opposed to depth)
-         ! This has been removed but can be reproduced using the "custom" module
-         ! if required
-         zsdmp = 1. / ( pn_surf * rday )
-         zbdmp = 1. / ( pn_bot  * rday )
-         zhfac = (  zbdmp + (zsdmp-zbdmp) * EXP( -gdept(1,1)/pn_dep )  )
+         zsdmp = 1._wp / ( pn_surf * rday )
+         zbdmp = 1._wp / ( pn_bot  * rday )
 
-         presto(:,:) = zmrs(:,:) * zhfac + ( 1. - zmrs(:,:) ) * presto(:,:)
+         ! 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 )         
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90	(revision 4739)
@@ -6,4 +6,5 @@
   PUBLIC
   
+  INTEGER, PUBLIC, PARAMETER   :: dp=8 , sp=4, wp=dp
   INTEGER :: tmask_id, umask_id, vmask_id, fmask_id
   INTEGER :: gdept_id
@@ -13,19 +14,19 @@
   INTEGER :: jpi, jpj, jpk                      ! Size of domain
   INTEGER  :: ncin, ncout                              ! File handles for netCDF files
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphit, glamt
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphiu, glamu
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphiv, glamv
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphif, glamf
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: tmask, umask, vmask, fmask
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gdept
-  REAL(8), DIMENSION(:,:),   ALLOCATABLE :: resto
+  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(8),PARAMETER :: rday = 86400           ! seconds in a day
-  REAL(8),PARAMETER ::   rpi = 3.141592653589793
-  REAL(8),PARAMETER ::   rad = 3.141592653589793/180.
-  REAL(8),PARAMETER ::   ra =  6371229.
+  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
@@ -42,9 +43,11 @@
   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_coast, ln_zero_top_layer, ln_custom, &
+                          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
 
@@ -106,5 +109,5 @@
      CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) )
 
-     CALL check_nf90( nf90_def_var(ncout, 'resto', nf90_float, (/id_x,id_y,id_z/), resto_id ) )
+     CALL check_nf90( nf90_def_var(ncout, 'resto', nf90_double, (/id_x,id_y,id_z/), resto_id ) )
      CALL check_nf90( nf90_enddef(ncout) )
 
Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90	(revision 4738)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90	(revision 4739)
@@ -18,10 +18,14 @@
       !! ** Action  : - resto, the damping coeff. for T and S
       !!----------------------------------------------------------------------
-      REAL(8), DIMENSION(jpi,jpj), INTENT(inout)  ::   presto   ! restoring coeff. (s-1)
-      REAL(8), DIMENSION(jpi,jpj), INTENT(in)  ::   mask   ! restoring coeff. (s-1)
+      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(8) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar
-      REAL(8), DIMENSION(6)  ::   zfact               ! 1Dworkspace
+      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
       !!----------------------------------------------------------------------
       !
@@ -29,14 +33,18 @@
       !
       
+      ! Read namelist
+      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
+      READ( numnam, nam_dmp_create )
+      CLOSE( numnam )
 
-      zfact(1) =  1.
-      zfact(2) =  1.
-      zfact(3) = 11. / 12.
-      zfact(4) =  8. / 12.
-      zfact(5) =  4. / 12.
-      zfact(6) =  1. / 12.
-      zfact(:) = zfact(:) / ( 5. * rday )    ! 5 days max restoring time scale
+      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.
+      presto(:,:) = 0._wp
 
       ! damping along the forced closed boundary over 6 grid-points
@@ -51,15 +59,15 @@
       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,*)
+         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.
-         zlat0 = 10.                     ! zlat0 : latitude strip where resto decreases
-         zlat1 = 30.                     ! zlat1 : resto = 1 before zlat1
+         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. / ( 5. * rday )   ! z1_5d : 1 / 5days
+         z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days
 
          DO jj = 1, jpj
@@ -67,5 +75,5 @@
                zlat = ABS( gphit(ji,jj) )
                IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
-                  presto(ji,jj) = 0.5 * z1_5d * (  1. - COS( rpi*(zlat2-zlat)/zlat0 )  ) 
+                  presto(ji,jj) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  ) 
                ELSEIF( zlat < zlat1 ) THEN
                   presto(ji,jj) = z1_5d
