Index: branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90
===================================================================
--- branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 5255)
+++ branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 5257)
@@ -44,7 +44,7 @@
       !!----------------------------------------------------------------------
       !
+      USE utils
       IMPLICIT NONE
 
-      USE utils
 
      !!----------------------------------------------------------------------
@@ -54,30 +54,30 @@
       READ( numnam, namzgr_sco )
       CLOSE( numnam )
-      jpk = rn_jpk
-
-      WRITE(numout,*)
-      WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
-      WRITE(numout,*) '~~~~~~~~~~~'
-      WRITE(numout,*) '   Namelist namzgr_sco'
-      WRITE(numout,*) '     stretching coeffs '
-      WRITE(numout,*) '        Number of levels                             rn_jpk        = ',rn_jpk
-      WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
-      WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
-      WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
-      WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
-      WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
-      WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
-      WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
-      WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
-      WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
-      WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
-      WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
-      WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
-      WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
-      WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
-      WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
-      WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
-      WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
-      WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
+      jpk = INT(rn_jpk)
+
+      WRITE(*,*)
+      WRITE(*,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
+      WRITE(*,*) '~~~~~~~~~~~'
+      WRITE(*,*) '   Namelist namzgr_sco'
+      WRITE(*,*) '     stretching coeffs '
+      WRITE(*,*) '        Number of levels                             rn_jpk        = ',rn_jpk
+      WRITE(*,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
+      WRITE(*,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
+      WRITE(*,*) '        Critical depth                               rn_hc         = ',rn_hc
+      WRITE(*,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
+      WRITE(*,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
+      WRITE(*,*) '        Song and Haidvogel 1994 stretching coefficients'
+      WRITE(*,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
+      WRITE(*,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
+      WRITE(*,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
+      WRITE(*,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
+      WRITE(*,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
+      WRITE(*,*) '        Siddorn and Furner 2012 stretching coefficients'
+      WRITE(*,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
+      WRITE(*,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
+      WRITE(*,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
+      WRITE(*,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
+      WRITE(*,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
+      WRITE(*,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
 
       
@@ -102,5 +102,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
+           IF( bathy(ji,jj) > 0. )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
          END DO
       END DO
@@ -114,5 +114,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-           IF( bathy(ji,jj) == 0._wp ) THEN
+           IF( bathy(ji,jj) == 0. ) THEN
              iip1 = MIN( ji+1, jpi )
              ijp1 = MIN( jj+1, jpj )
@@ -120,5 +120,5 @@
              ijm1 = MAX( jj-1, 1 )
              IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
-        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
+        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. ) THEN
                zenv(ji,jj) = rn_sbot_min
              ENDIF
@@ -128,13 +128,13 @@
       ! 
       ! smooth the bathymetry (if required)
-      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
+      scosrf(:,:) = 0.             ! ocean surface depth (here zero: no under ice-shelf sea)
       scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
       !
       jl = 0
-      zrmax = 1._wp
+      zrmax = 1.
       !   
       !     
       ! set scaling factor used in reducing vertical gradients
-      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
+      zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax )
       !
       ! initialise temporary evelope depth arrays
@@ -145,30 +145,48 @@
       !
       ! initialise temporary r-value arrays
-      zri(:,:) = 1._wp
-      zrj(:,:) = 1._wp
+      zri(:,:) = 1.
+      zrj(:,:) = 1.
       !                                                            ! ================ !
-      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
+      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) !  Iterative loop  !
          !                                                         ! ================ !
          jl = jl + 1
-         zrmax = 0._wp
+         zrmax = 0.
          ! we set zrmax from previous r-values (zri and zrj) first
          ! if set after current r-value calculation (as previously)
          ! we could exit DO WHILE prematurely before checking r-value
          ! of current zenv
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
+!         DO jj = 1, nlcj
+!            DO ji = 1, nlci
+         DO jj = 1, jpi !jpi or jpim1?
+            DO ji = 1, jpj
                zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
             END DO
          END DO
-         zri(:,:) = 0._wp
-         zrj(:,:) = 0._wp
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
-               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
-               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
+         zri(:,:) = 0.
+         zrj(:,:) = 0.
+!         DO jj = 1, nlci
+ !           DO ji = 1, nlcj
+ !              iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
+ !              ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
+ !              IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN
+ !                 zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+ !              END IF
+ !              IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN
+ !                 zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+ !              END IF
+ !              IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
+ !              IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
+ !              IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
+ !              IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
+ !           END DO
+ !        END DO
+         DO jj = 1, jpi-1
+            DO ji = 1, jpj-1
+               iip1 = ji+1      
+               ijp1 = jj+1      
+               IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN
                   zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
                END IF
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
+               IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN
                   zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
                END IF
@@ -182,10 +200,9 @@
          WRITE(*,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
          !
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
+         DO jj = 1, jpi
+            DO ji = 1, jpj
                zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
             END DO
          END DO
-         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
          !                                                  ! ================ !
       END DO                                                !     End loop     !
@@ -198,28 +215,29 @@
       !
       ! Envelope bathymetry saved in hbatt
+      ! TODO - get this section to work
       hbatt(:,:) = zenv(:,:) 
-      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
-         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
-               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
-            END DO
-         END DO
-      ENDIF
+!      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) THEN
+ !        CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
+ !        DO jj = 1, jpj
+ !           DO ji = 1, jpi
+ !              ztaper = EXP( -(gphit(ji,jj)/8.)**2. )
+ !              hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper )
+ !           END DO
+ !        END DO
+ !     ENDIF
       !
       !                                        ! ==============================
       !                                        !   hbatu, hbatv, hbatf fields
       !                                        ! ==============================
-      WRITE(numout,*)
-      WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+      WRITE(*,*)
+      WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
       hbatu(:,:) = rn_sbot_min
       hbatv(:,:) = rn_sbot_min
       hbatf(:,:) = rn_sbot_min
-      DO jj = 1, jpjm1
-        DO ji = 1, jpim1   ! NO vector opt.
-           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
-           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
-           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
+      DO jj = 1, jpj-1
+        DO ji = 1, jpi-1   ! NO vector opt.
+           hbatu(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
+           hbatv(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
+           hbatf(ji,jj) = 0.25 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
               &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
         END DO
@@ -229,7 +247,7 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( hbatu(ji,jj) == 0._wp ) THEN
-               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
-               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
+            IF( hbatu(ji,jj) == 0. ) THEN
+               IF( zhbat(ji,jj) == 0. )   hbatu(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0. )   hbatu(ji,jj) = zhbat(ji,jj)
             ENDIF
          END DO
@@ -238,7 +256,7 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( hbatv(ji,jj) == 0._wp ) THEN
-               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
-               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
+            IF( hbatv(ji,jj) == 0. ) THEN
+               IF( zhbat(ji,jj) == 0. )   hbatv(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0. )   hbatv(ji,jj) = zhbat(ji,jj)
             ENDIF
          END DO
@@ -247,7 +265,7 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( hbatf(ji,jj) == 0._wp ) THEN
-               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
-               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
+            IF( hbatf(ji,jj) == 0. ) THEN
+               IF( zhbat(ji,jj) == 0. )   hbatf(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0. )   hbatf(ji,jj) = zhbat(ji,jj)
             ENDIF
          END DO
@@ -260,11 +278,11 @@
       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
 
-         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
+         WRITE(*,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
             &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
-         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
+         WRITE(*,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
             &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
-         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
+         WRITE(*,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
             &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
-         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
+         WRITE(*,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
             &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
 !! helsinki
@@ -326,23 +344,23 @@
          DO ji = 1, jpi
             DO jk = 1, jpk-1
-               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
-            END DO
-            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
-         END DO
-      END DO
-      WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
+            END DO
+            IF( scobot(ji,jj) == 0.               )   mbathy(ji,jj) = 0
+         END DO
+      END DO
+      WRITE(*,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
 
       ! min max values over the domain
-      WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
-      WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
+      WRITE(*,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
+      WRITE(*,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
          &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
-      WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
+      WRITE(*,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
          &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
          &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
          &                          ' w ', MINVAL( e3w_0  (:,:,:) )
 
-      WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
+      WRITE(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
          &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
-      WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
+      WRITE(*,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
          &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
          &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
@@ -350,27 +368,27 @@
       
          ! selected vertical profiles
-      WRITE(numout,*)
-      WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
-      WRITE(numout,*) ' ~~~~~~  --------------------'
-      WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
-      WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
+      WRITE(*,*)
+      WRITE(*,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
+      WRITE(*,*) ' ~~~~~~  --------------------'
+      WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+      WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
          &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
-      DO jj = mj0(20), mj1(20)
-         DO ji = mi0(20), mi1(20)
-            WRITE(numout,*)
-            WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
-            WRITE(numout,*) ' ~~~~~~  --------------------'
-            WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
-            WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
+      DO jj = 20, 20
+         DO ji = 20, 20
+            WRITE(*,*)
+            WRITE(*,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
+            WRITE(*,*) ' ~~~~~~  --------------------'
+            WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+            WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
                &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
          END DO
       END DO
-      DO jj = mj0(74), mj1(74)
-         DO ji = mi0(100), mi1(100)
-            WRITE(numout,*)
-            WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
-            WRITE(numout,*) ' ~~~~~~  --------------------'
-            WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
-            WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
+      DO jj = 74,74
+         DO ji = 100, 100
+            WRITE(*,*)
+            WRITE(*,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
+            WRITE(*,*) ' ~~~~~~  --------------------'
+            WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+            WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
                &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
          END DO
@@ -383,25 +401,25 @@
          DO jj = 1, jpj
 
-            IF( hbatt(ji,jj) > 0._wp) THEN
+            IF( hbatt(ji,jj) > 0.) THEN
                DO jk = 1, mbathy(ji,jj)
                  ! check coordinate is monotonically increasing
-                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
-                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
-                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
-                    CALL STOP 1
+                 IF (e3w_0(ji,jj,jk) <= 0. .OR. e3t_0(ji,jj,jk) <= 0. ) THEN
+                    WRITE(*,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(*,*) 'e3w',e3w_0(ji,jj,:)
+                    WRITE(*,*) 'e3t',e3t_0(ji,jj,:)
+                    STOP 1
                  ENDIF
                  ! and check it has never gone negative
-                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
-                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
-                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
-                    CALL STOP 1
+                 IF( gdepw_0(ji,jj,jk) < 0. .OR. gdept_0(ji,jj,jk) < 0. ) THEN
+                    WRITE(*,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:)
+                    WRITE(*,*) 'gdept',gdept_0(ji,jj,:)
+                    STOP 1
                  ENDIF
                  ! and check it never exceeds the total depth
-                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
-                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
-                    CALL STOP 1
+                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                    WRITE(*,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:)
+                    STOP 1
                  ENDIF
                END DO
@@ -409,8 +427,8 @@
                DO jk = 1, mbathy(ji,jj)-1
                  ! and check it never exceeds the total depth
-                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
-                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
-                    CALL STOP 1
+                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                    WRITE(*,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(*,*) 'gdept',gdept_0(ji,jj,:)
+                    STOP 1
                  ENDIF
                END DO
@@ -421,5 +439,9 @@
       END DO
       !
-      !
+      ! Write output file
+      CALL write_coord_file()
+      !
+      !
+END PROGRAM SCOORD_GEN
 
 !!======================================================================
@@ -437,5 +459,5 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      USE utils
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
       !
@@ -449,8 +471,8 @@
       ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
 
-      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
-      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
-      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
-      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
+      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
+      z_esigt3  = 0.   ;   z_esigw3  = 0. 
+      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
+      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
 
       DO ji = 1, jpi
@@ -459,5 +481,5 @@
             IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
                DO jk = 1, jpk
-                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
+                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5, rn_bb )
                   z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
                END DO
@@ -465,5 +487,5 @@
                DO jk = 1, jpk
                   z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
-                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
+                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp)
                   END DO
             ENDIF
@@ -473,9 +495,9 @@
                z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
             END DO
-            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
-            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
+            z_esigw3(ji,jj,1  ) = 2. * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
+            z_esigt3(ji,jj,jpk) = 2. * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
             !
             ! Coefficients for vertical depth as the sum of e3w scale factors
-            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
+            z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
             DO jk = 2, jpk
                z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
@@ -483,8 +505,8 @@
             !
             DO jk = 1, jpk
-               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp)
-               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp)
-               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
-               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
+               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
+               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
                gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
             END DO
@@ -543,5 +565,5 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      USE utils
       REAL(wp) ::   zsmth               ! smoothing around critical depth
       REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
@@ -556,8 +578,8 @@
       ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
 
-      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
-      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
-      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
-      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
+      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
+      z_esigt3  = 0.   ;   z_esigw3  = 0. 
+      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
+      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
 
       DO ji = 1, jpi
@@ -568,22 +590,22 @@
               zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
                                                      ! could be changed by users but care must be taken to do so carefully
-              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
+              zzb = 1.0-(zzb/hbatt(ji,jj))
             
               zzs = rn_zs / hbatt(ji,jj) 
               
-              IF (rn_efold /= 0.0_wp) THEN
+              IF (rn_efold /= 0.0) THEN
                 zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
               ELSE
-                zsmth = 1.0_wp 
+                zsmth = 1.0 
               ENDIF
                
               DO jk = 1, jpk
                 z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
-                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
+                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
               ENDDO
               z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
               z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
  
-          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
+          ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma
 
             DO jk = 1, jpk
@@ -596,5 +618,5 @@
             DO jk = 1, jpk
               z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
-              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
             END DO
 
@@ -605,6 +627,6 @@
              z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
           END DO
-          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
-          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
+          z_esigw3(ji,jj,1  ) = 2.0 * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
+          z_esigt3(ji,jj,jpk) = 2.0 * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
 
           ! Coefficients for vertical depth as the sum of e3w scale factors
@@ -615,6 +637,6 @@
 
           DO jk = 1, jpk
-             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
-             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
+             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
+             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
              gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
           END DO
@@ -671,5 +693,5 @@
       !!----------------------------------------------------------------------
 
-      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      USE utils
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
 
@@ -680,9 +702,9 @@
       ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
 
-      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
-      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
+      z_gsigw  = 0.   ;   z_gsigt  = 0.   ;   z_gsi3w  = 0.
+      z_esigt  = 0.   ;   z_esigw  = 0. 
 
       DO jk = 1, jpk
-        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
+        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 )
         z_gsigt(jk) = -fssig( REAL(jk,wp)        )
       END DO
@@ -696,9 +718,9 @@
          z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
       END DO
-      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
-      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
+      z_esigw( 1 ) = 2. * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
+      z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) )
       !
       ! Coefficients for vertical depth as the sum of e3w scale factors
-      z_gsi3w(1) = 0.5_wp * z_esigw(1)
+      z_gsi3w(1) = 0.5 * z_esigw(1)
       DO jk = 2, jpk
          z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
@@ -706,8 +728,8 @@
 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
       DO jk = 1, jpk
-         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp)
-         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp)
-         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
-         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
+         zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
+         zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
+         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
          gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
       END DO
@@ -721,5 +743,5 @@
               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
               !
-              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
@@ -744,13 +766,14 @@
       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
       !!----------------------------------------------------------------------
+      USE utils, ONLY : wp
       REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
       REAL(wp)             ::   pf   ! sigma value
       !!----------------------------------------------------------------------
       !
-      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpk-1) + rn_thetb )  )   &
+      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb )  )   &
          &     - TANH( rn_thetb * rn_theta                                )  )   &
          & * (   COSH( rn_theta                           )                      &
-         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
-         & / ( 2._wp * SINH( rn_theta ) )
+         &     + COSH( rn_theta * ( 2. * rn_thetb - 1. ) )  )              &
+         & / ( 2. * SINH( rn_theta ) )
       !
    END FUNCTION fssig
@@ -768,4 +791,5 @@
       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
       !!----------------------------------------------------------------------
+      USE utils, ONLY : wp
       REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
       REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
@@ -774,9 +798,9 @@
       !
       IF ( rn_theta == 0 ) then      ! uniform sigma
-         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpk-1 )
+         pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 )
       ELSE                        ! stretched sigma
-         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
-            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpk-1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
-            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
+         pf1 =   ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
+            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta )  )  &
+            &        / ( 2. * TANH( 0.5 * rn_theta ) )  )
       ENDIF
       !
@@ -801,4 +825,5 @@
       !! Reference  :   Siddorn and Furner, in prep
       !!----------------------------------------------------------------------
+      USE utils, ONLY : jpk,jk,wp
       REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
       REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
@@ -809,5 +834,4 @@
       REAL(wp)                ::   zn1,zn2        ! local variables
       REAL(wp)                ::   za,zb,zx       ! local variables
-      integer                 ::   jk
       !!----------------------------------------------------------------------
       !
@@ -816,18 +840,18 @@
       zn2  =  1. -  zn1
 
-      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
-      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
-      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
+      za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) 
+      za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0)
+      za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1)
      
       za = pzb - za3*(pzs-za1)-za2
-      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
-      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
-      zx = 1.0_wp-za/2.0_wp-zb
+      za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) )
+      zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1)
+      zx = 1.0-za/2.0-zb
  
       DO jk = 1, jpk
-        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
-                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
-                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
-        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
+        p_gamma(jk) = za*(pk1(jk)*(1.0-pk1(jk)/2.0))+zb*pk1(jk)**3.0 +  &
+                    & zx*( (rn_alpha+2.0)*pk1(jk)**(rn_alpha+1.0)- &
+                    &      (rn_alpha+1.0)*pk1(jk)**(rn_alpha+2.0) )
+        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0-psmth)
       ENDDO 
 
@@ -835,3 +859,2 @@
    END FUNCTION fgamma
 
-END PROGRAM SCOORD_GEN
Index: branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90
===================================================================
--- branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 5255)
+++ branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 5257)
@@ -1,13 +1,13 @@
 MODULE utils
 
-   IMPLICIT NONE
    USE netcdf
 
+   IMPLICIT NONE
    PUBLIC             ! allows the acces to par_oce when dom_oce is used
    !                  ! exception to coding rules... to be suppressed ???
 
-   PUBLIC dom_oce_alloc
-   PUBLIC read_bathy
+!   PUBLIC dom_oce_alloc
 
+   INTEGER, PARAMETER   :: dp=8 , sp=4, wp=dp
 
    !! All coordinates
@@ -45,26 +45,33 @@
    INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
    INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation
+   INTEGER  ::   numnam                   ! File handle for namelist 
    REAL(wp) ::   zrmax, ztaper   ! temporary scalars
    REAL(wp) ::   zrfact
    !
-   REAL(wp), DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
-   REAL(wp), DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zri, zrj, zhbat
+
+   !Namelist variables
+   REAL(wp) :: rn_jpk, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax, rn_theta
+   REAL(wp) :: rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
+   LOGICAL :: ln_s_sh94, ln_s_sf12, ln_sigcrit
 
    NAMELIST/namzgr_sco/rn_jpk, ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
                         rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
 
+   CONTAINS
 
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(12) :: ierr
+      INTEGER, DIMENSION(4) :: ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
       !
       ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), &
-         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )
+         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) )
          !
-      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         &
-         &      gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         &
-         &      gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) )
+      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) ,                         &
+         &      gdept_0(jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         &
+         &      gdepw_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(2) )
          !
          !
@@ -74,18 +81,18 @@
          &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
          &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
-         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) )
+         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(3) )
 
-      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(9) )
-      !
+      ALLOCATE( mbathy(jpi,jpj) , STAT=ierr(4) )
+     !
       dom_oce_alloc = MAXVAL(ierr)
       !
    END FUNCTION dom_oce_alloc
-
+   
 
    SUBROUTINE read_bathy()
      !! Read bathymetry from input netcdf file
-     INTEGER :: var_id
+     INTEGER :: var_id, ncin
 
-     CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening mesh_mask file' )
+     CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening bathy.nc file' )
 
      ! Find the size of the input bathymetry
@@ -96,6 +103,8 @@
      
      ! Read the bathymetry variable from file
-     CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', tmask_id ), 'Cannot get variable ID for bathymetry')
+     CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', var_id ), 'Cannot get variable ID for bathymetry')
      CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+     CALL check_nf90( nf90_close(ncin), 'Error closing bathy.nc file' )
 
    END SUBROUTINE read_bathy
@@ -113,9 +122,52 @@
      CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len))
 
-  END SUBROUTINE dimlen
+   END SUBROUTINE dimlen
+  
+ 
+   SUBROUTINE write_coord_file()
+     ! Write out variables to the a netcdf coordinates file
+     
+     INTEGER :: id_x, id_y, id_z
+     INTEGER :: ncout
+     INTEGER, DIMENSION(11) :: var_ids  !Array to contain all variable IDs
 
+     !Create the file
+     CALL check_nf90( nf90_create('coord_zgr.nc', NF90_CLOBBER, ncout), 'Could not create output file')
+     !
+     !Define dimensions
+     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) )
+     !
+     !Define variables
+     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_x/), var_ids(1)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(2)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(3)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_x/), var_ids(4)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_x/), var_ids(5)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_x/), var_ids(6)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_x/), var_ids(7)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(8)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(9)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(10)) )
+     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y,id_x/), var_ids(11)) )
+     
+     ! End define mode
+     CALL check_nf90( nf90_enddef(ncout) )
 
-   
-   SUBROUTINE write_coord_file()
+     !Write variables to file
+     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) )
+     
+     CALL check_nf90( nf90_close(ncout) )
 
    END SUBROUTINE write_coord_file
@@ -127,6 +179,6 @@
 
       IF (istat /= nf90_noerr) THEN
-         WRITE(numerr,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
-         IF ( PRESENT(message) ) THEN ; WRITE(numerr,*) message ; ENDIF
+         WRITE(*,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
+         IF ( PRESENT(message) ) THEN ; WRITE(*,*) message ; ENDIF
          STOP
       ENDIF
