Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 2415)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 2416)
@@ -17,24 +17,24 @@
 #if   defined key_tradmp   ||   defined key_esopa
    !!----------------------------------------------------------------------
-   !!   key_tradmp                                         internal damping
+   !!   'key_tradmp'                                       internal damping
    !!----------------------------------------------------------------------
-   !!   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
+   !!   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 dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables
-   USE trdmod_oce         ! ocean space and time domain variables
-   USE trdtra         ! ocean space and time domain variables
-   USE zdf_oce         ! ocean vertical physics
-   USE phycst          ! Define parameters for the routines
-   USE dtatem          ! temperature data
-   USE dtasal          ! salinity data
-   USE zdfmxl          ! mixed layer depth
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing
-   USE prtctl          ! Print control
+   USE oce            ! ocean: variables
+   USE dom_oce        ! ocean: domain variables
+   USE trdmod_oce     ! ocean: trend variables
+   USE trdtra         ! active tracers: trends
+   USE zdf_oce        ! ocean: vertical physics
+   USE phycst         ! physical constants
+   USE dtatem         ! data: temperature
+   USE dtasal         ! data: salinity
+   USE zdfmxl         ! vertical physics: mixed layer depth
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -43,6 +43,6 @@
    PUBLIC   tra_dmp      ! routine called by step.F90
    PUBLIC   tra_dmp_init ! routine called by opa.F90
-   PUBLIC   dtacof       ! routine called by tradmp.F90 and trcdmp.F90
-   PUBLIC   dtacof_zoom  ! routine called by tradmp.F90 and trcdmp.F90
+   PUBLIC   dtacof       ! routine called by in both tradmp.F90 and trcdmp.F90
+   PUBLIC   dtacof_zoom  ! routine called by in both tradmp.F90 and trcdmp.F90
 
 #if ! defined key_agrif
@@ -52,14 +52,14 @@
 #endif
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   strdmp   !: damping salinity trend (psu/s)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ttrdmp   !: damping temperature trend (Centigrade/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ttrdmp   !: damping temperature trend (Celcius/s)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   resto    !: restoring coeff. on T and S (s-1)
    
-   !                             !!* Namelist namtra_dmp : T & S newtonian damping *
-   INTEGER  ::   nn_hdmp =   -1   ! = 0/-1/'latitude' for damping over T and S
-   INTEGER  ::   nn_zdmp =    0   ! = 0/1/2 flag for damping in the mixed layer
-   REAL(wp) ::   rn_surf =   50.  ! surface time scale for internal damping        [days]
-   REAL(wp) ::   rn_bot  =  360.  ! bottom time scale for internal damping         [days]
-   REAL(wp) ::   rn_dep  =  800.  ! depth of transition between rn_surf and rn_bot [meters]
-   INTEGER  ::   nn_file =    2   ! = 1 create a damping.coeff NetCDF file 
+   !                                !!* Namelist namtra_dmp : T & S newtonian damping *
+   INTEGER  ::   nn_hdmp =   -1      ! = 0/-1/'latitude' for damping over T and S
+   INTEGER  ::   nn_zdmp =    0      ! = 0/1/2 flag for damping in the mixed layer
+   REAL(wp) ::   rn_surf =   50._wp  ! surface time scale for internal damping        [days]
+   REAL(wp) ::   rn_bot  =  360._wp  ! bottom time scale for internal damping         [days]
+   REAL(wp) ::   rn_dep  =  800._wp  ! depth of transition between rn_surf and rn_bot [meters]
+   INTEGER  ::   nn_file =    2      ! = 1 create a damping.coeff NetCDF file 
 
    !! * Substitutions
@@ -69,7 +69,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -94,15 +93,9 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      REAL(wp) ::   zta, zsa       ! temporary scalars
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
-      !!----------------------------------------------------------------------
-
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
-      ENDIF
-
-      SELECT CASE ( nn_zdmp ) 
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zta, zsa     ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( nn_zdmp )     !==    type of damping   ==!
       !
       CASE( 0 )                   !==  newtonian damping throughout the water column  ==!
@@ -114,6 +107,5 @@
                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
-                  ! save the salinity trend (used in asmtrj)
-                  strdmp(ji,jj,jk) = zsa
+                  strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj)
                   ttrdmp(ji,jj,jk) = zta
                END DO
@@ -125,15 +117,14 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( avt(ji,jj,jk) <= 5.e-4 ) THEN
+                  IF( avt(ji,jj,jk) <= 5.e-4_wp ) THEN
                      zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
                      zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
                   ELSE
-                     zta = 0.e0
-                     zsa = 0.e0   
+                     zta = 0._wp
+                     zsa = 0._wp  
                   ENDIF
                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
-                  ! save the salinity trend (used in asmtrj)
-                  strdmp(ji,jj,jk) = zsa
+                  strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj)
                   ttrdmp(ji,jj,jk) = zta
                END DO
@@ -149,12 +140,10 @@
                      zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
                   ELSE
-                     zta = 0.e0
-                     zsa = 0.e0   
+                     zta = 0._wp
+                     zsa = 0._wp  
                   ENDIF
-                  ! add the trends to the general tracer trends
                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
-                  ! save the salinity trend (used in asmtrj)
-                  strdmp(ji,jj,jk) = zsa
+                  strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj)
                   ttrdmp(ji,jj,jk) = zta
                END DO
@@ -163,11 +152,8 @@
          !
       END SELECT
-
+      !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_dmp, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_dmp, ztrds )
-         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_dmp, ttrdmp )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_dmp, strdmp )
       ENDIF
       !                           ! Control print
@@ -191,5 +177,6 @@
       REWIND ( numnam )                  ! Read Namelist namtra_dmp : temperature and salinity damping term
       READ   ( numnam, namtra_dmp )
-      IF( lzoom )   nn_zdmp = 0           ! restoring to climatology at closed north or south boundaries
+      
+      IF( lzoom )   nn_zdmp = 0          ! restoring to climatology at closed north or south boundaries
 
       IF(lwp) THEN                       ! Namelist print
@@ -226,6 +213,6 @@
          &   CALL ctl_stop( 'no temperature and/or salinity data define key_dtatem and key_dtasal' )
 
-      strdmp(:,:,:) = 0.e0       ! internal damping salinity trend (used in asmtrj)
-      ttrdmp(:,:,:) = 0.e0
+      strdmp(:,:,:) = 0._wp       ! internal damping salinity trend (used in asmtrj)
+      ttrdmp(:,:,:) = 0._wp
       !                          ! Damping coefficients initialization
       IF( lzoom ) THEN   ;   CALL dtacof_zoom( resto )
@@ -250,21 +237,20 @@
       !! ** Action  : - resto, the damping coeff. for T and S
       !!----------------------------------------------------------------------
-      !! * Arguments
-      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   ! temporary scalar
-      REAL(wp), DIMENSION(6)  ::   zfact        ! temporary workspace
-      !!----------------------------------------------------------------------
-
-      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
-
-      presto(:,:,:) = 0.e0
+      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
+      !!----------------------------------------------------------------------
+
+      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
@@ -285,8 +271,9 @@
          !
          !                          ! Initialization : 
-         presto(:,:,:) = 0.e0
-         zlat0 = 10.                     ! zlat0 : latitude strip where resto decreases
-         zlat1 = 30.                     ! zlat1 : resto = 1 before zlat1
-         zlat2 = zlat1 + zlat0           ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
+         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
@@ -295,7 +282,7 @@
                   zlat = ABS( gphit(ji,jj) )
                   IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
-                     presto(ji,jj,jk) = 0.5 * ( 1./(5.*rday) ) * ( 1. - cos(rpi*(zlat2-zlat)/zlat0) ) 
+                     presto(ji,jj,jk) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  ) 
                   ELSEIF( zlat < zlat1 ) THEN
-                     presto(ji,jj,jk) = 1./(5.*rday)
+                     presto(ji,jj,jk) = z1_5d
                   ENDIF
                END DO
@@ -325,20 +312,19 @@
       USE iom
       USE ioipsl
-      !! * Arguments
-      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
+      !!
+      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 or TRC (tracer 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           !    -          -
-      INTEGER ::   inum0                        ! logical unit for file restoring damping term
-      INTEGER ::   icot                         ! logical unit for file distance to the coast
-      REAL(wp) ::   zinfl, zlon                 ! temporary scalars
-      REAL(wp) ::   zlat, zlat0, zlat1, zlat2   !    -         -
-      REAL(wp) ::   zsdmp, zbdmp                !    -         -
+      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                !   -      -
       REAL(wp), DIMENSION(jpk)         ::   zhfac
       REAL(wp), DIMENSION(jpi,jpj)     ::   zmrs
@@ -347,19 +333,18 @@
       !!----------------------------------------------------------------------
 
-      ! ====================================
-      !  ORCA configuration : global domain
-      ! ====================================
-
+      !                                   ! ====================
+      !                                   !  ORCA configuration : global domain
+      !                                   ! ====================
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) '              dtacof : Global domain of ORCA'
       IF(lwp) WRITE(numout,*) '              ------------------------------'
-
-      ! ... Initialization : 
-      presto(:,:,:) = 0.e0
+      !
+      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.'
+         IF(lwp) WRITE(numout,*) '              Damping poleward of ', kn_hdmp, ' deg.'
          !
          CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )
@@ -373,6 +358,6 @@
 
          !                            ! Compute arrays resto 
-         zinfl = 1000.e3                   ! distance of influence for damping term
-         zlat0 = 10.                       ! latitude strip where resto decreases
+         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|
@@ -382,7 +367,7 @@
                zlat = ABS( gphit(ji,jj) )
                IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
-                  presto(ji,jj,1) = 0.5 * ( 1. - cos(rpi*(zlat-zlat1)/zlat0 ) )
+                  presto(ji,jj,1) = 0.5_wp * (  1._wp - COS( rpi*(zlat-zlat1)/zlat0 )  )
                ELSEIF ( zlat > zlat2 ) THEN
-                  presto(ji,jj,1) = 1.
+                  presto(ji,jj,1) = 1._wp
                ENDIF
             END DO
@@ -393,7 +378,7 @@
                DO ji = 1, jpi
                   zlat = gphit(ji,jj)
-                  zlon = MOD( glamt(ji,jj), 360. )
-                  IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45. < zlon .AND. zlon < 100. ) THEN
-                     presto(ji,jj,1) = 0.e0
+                  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
@@ -401,6 +386,6 @@
          ENDIF
 
-         zsdmp = 1./(pn_surf * rday)
-         zbdmp = 1./(pn_bot  * rday)
+         zsdmp = 1._wp / ( pn_surf * rday )
+         zbdmp = 1._wp / ( pn_bot  * rday )
          DO jk = 2, jpkm1
             DO jj = 1, jpj
@@ -408,7 +393,7 @@
                   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 * ( 1. - COS( rpi*zdct(ji,jj,jk)/zinfl) )
+                  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) )
+                  presto(ji,jj,jk) = presto(ji,jj,jk) * (  zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)  )
                END DO
             END DO
@@ -417,98 +402,93 @@
       ENDIF
 
-
+      !                                  ! =========================
+      !                                  !  Med and Red Sea damping    (ORCA configuration only)
+      !                                  ! =========================
       IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN
-
-         !                                         ! =========================
-         !                                         !  Med and Red Sea damping
-         !                                         ! =========================
          IF(lwp)WRITE(numout,*)
          IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas'
-
-
-         zmrs(:,:) = 0.e0                             ! damping term on the Med or Red Sea
-
+         !
+         zmrs(:,:) = 0._wp
+         !
          SELECT CASE ( jp_cfg )
          !                                           ! =======================
          CASE ( 4 )                                  !  ORCA_R4 configuration 
             !                                        ! =======================
-            ! Mediterranean Sea
-            ij0 =  50   ;   ij1 =  56
-            ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
+            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.e0
+            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.e0
+            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*( 1.- COS( rpi*(jk-1)/16. ) ) / rday
+               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./rday
+               zhfac (jk) = 1._wp / rday
             END DO
             !                                        ! =======================
          CASE ( 2 )                                  !  ORCA_R2 configuration 
             !                                        ! =======================
-            ! Mediterranean Sea
-            ij0 =  96   ;   ij1 = 110
-            ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
+            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.e0
+            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.e0
-            ! Decrease before Gibraltar Strait
-            ij0 = 101   ;   ij1 = 102
-            ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
-            ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
-            ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
-            ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75e0
-            ! Red Sea
-            ij0 =  87   ;   ij1 =  96
-            ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
-            ! Decrease before Bab el Mandeb Strait
-            ij0 =  91   ;   ij1 =  91
-            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80e0
+            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.40e0
+            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.e0 / 90.e0
+            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.e0
+            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*( 1.- COS( rpi*(jk-1)/16. ) ) / rday
+               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./rday
+               zhfac (jk) = 1._wp / rday
             END DO
             !                                        ! =======================
          CASE ( 05 )                                 !  ORCA_R05 configuration
             !                                        ! =======================
-            ! Mediterranean Sea
-            ii0 = 568   ;   ii1 = 574 
-            ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
+            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.e0
-            ! Black Sea (remaining part
-            ii0 = 641   ;   ii1 = 651
-            ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
-            ! Decrease before Gibraltar Strait
-            ij0 = 324   ;   ij1 = 333
-            ii0 = 565   ;   ii1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
-            ii0 = 566   ;   ii1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
-            ii0 = 567   ;   ii1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75
-            ! Red Sea
-            ii0 = 641   ;   ii1 = 665
-            ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
-            ! Decrease before Bab el Mandeb Strait
-            ii0 = 666   ;   ii1 = 675
+            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 * ABS( FLOAT(ji - mi1(ii1)) )
+               zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) )
             END DO 
-            zsdmp = 1./(pn_surf * rday)
-            zbdmp = 1./(pn_bot  * rday)
+            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) )
+               zhfac(jk) = (  zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep )  )
             END DO
             !                                       ! ========================
@@ -520,11 +500,11 @@
 
          DO jk = 1, jpkm1
-            presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1. - zmrs(:,:) ) * presto(:,:,jk)
+            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.e0
-         presto(:,:,jpk) = 0.e0
+         presto(:,:, 1 ) = 0._wp
+         presto(:,:,jpk) = 0._wp
          !                         !--------------------!
       ELSE                         !     No damping     !
@@ -594,8 +574,8 @@
            &               '         create the "dist.coast.nc" file using IDL' )
 
-      pdct(:,:,:) = 0.e0
-      zxt(:,:) = cos( rad * gphit(:,:) ) * cos( rad * glamt(:,:) )
-      zyt(:,:) = cos( rad * gphit(:,:) ) * sin( rad * glamt(:,:) )
-      zzt(:,:) = sin( rad * gphit(:,:) )
+      pdct(:,:,:) = 0._wp
+      zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )
+      zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )
+      zzt(:,:) = SIN( rad * gphit(:,:) )
 
 
@@ -610,7 +590,7 @@
                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. ) 
-               llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1. ) 
-               llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )
+               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
@@ -678,5 +658,4 @@
          ! Compute cartesian coordinates of coastline points
          ! and the number of coastline points
-
          icoast = 0
          DO jj = 1, jpj
@@ -704,9 +683,8 @@
 
          ! Distance for the T-points
-
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( tmask(ji,jj,jk) == 0. ) THEN
-                  pdct(ji,jj,jk) = 0.
+               IF( tmask(ji,jj,jk) == 0._wp ) THEN
+                  pdct(ji,jj,jk) = 0._wp
                ELSE
                   DO jl = 1, icoast
@@ -727,6 +705,6 @@
       ! ----------------------------------------------------------    
       clname = 'dist.coast'
-      itime = 0
-      CALL ymds2ju( 0     , 1      , 1     , 0.e0 , zdate0 )
+      itime  = 0
+      CALL ymds2ju( 0     , 1      , 1     , 0._wp , zdate0 )
       CALL restini( 'NONE', jpi    , jpj   , glamt, gphit ,   &
          &          jpk   , gdept_0, clname, itime, zdate0,   &
@@ -734,5 +712,5 @@
       CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
       CALL restclo( icot )
-
+      !
    END SUBROUTINE cofdis
 
