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 5268)
+++ branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 5269)
@@ -286,5 +286,7 @@
          WRITE(*,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
             &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
-!! helsinki
+      
+      ! Create the output file
+      CALL make_coord_file()
 
       !                                            ! =======================
@@ -293,4 +295,5 @@
       !
       ! non-dimensional "sigma" for model level depth at w- and t-levels
+
 
 !========================================================================
@@ -298,4 +301,5 @@
 ! Siddorn and Furner 2012 (ln_sf12=T)
 ! or  tanh function       (both false)                    
+! To reduce memory loop over jpk and write each level to file
 !========================================================================
       IF      ( ln_s_sh94 ) THEN 
@@ -305,141 +309,10 @@
       ELSE                 
                            CALL s_tanh()
-      ENDIF 
-
-      !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
-
-#if defined key_agrif
-      ! Ensure meaningful vertical scale factors in ghost lines/columns
-      ! TODO: How can we deal with this in offline tool? 
-      IF( .NOT. Agrif_Root() ) THEN
-         !  
-         IF((nbondi == -1).OR.(nbondi == 2)) THEN
-            e3u_0(1,:,:) = e3u_0(2,:,:)
-         ENDIF
-         !
-         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
-            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
-         ENDIF
-         !
-         IF((nbondj == -1).OR.(nbondj == 2)) THEN
-            e3v_0(:,1,:) = e3v_0(:,2,:)
-         ENDIF
-         !
-         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
-            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
-         ENDIF
-         !
       ENDIF
-#endif
-!!
-      ! HYBRID : 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            DO jk = 1, jpk-1
-               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(*,*) ' 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(*,*) ' 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(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
-         &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_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  (:,:,:) ),   &
-         &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
-      
-         ! selected vertical profiles
-      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 = 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 = 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
-      END DO
-
-!================================================================================
-! check the coordinate makes sense
-!================================================================================
-      DO ji = 1, jpi
-         DO jj = 1, jpj
-
-            IF( hbatt(ji,jj) > 0.) THEN
-               DO jk = 1, mbathy(ji,jj)
-                 ! check coordinate is monotonically increasing
-                 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( 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( 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
-
-               DO jk = 1, mbathy(ji,jj)-1
-                 ! and check it never exceeds the total depth
-                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
-
-            ENDIF
-
-         END DO
-      END DO
-      !
-      ! Write output file
-      CALL write_coord_file()
-      !
+
+      CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) )
+      CALL check_nf90( nf90_close(ncout) )
+
+ 
       !
 END PROGRAM SCOORD_GEN
@@ -462,12 +335,141 @@
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
       !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
-
-      ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) )
-      ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 
-      ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 
-      ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
+      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
+      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
+      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
+      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
+ 
+      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
+      z_esigt3  = 0.   ;   z_esigw3  = 0. 
+      z_esigt3p1  = 0. ;   z_esigw3p1  = 0. 
+      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
+      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
+
+      DO jk = 1,jpk
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+
+               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
+                     z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb )
+                     z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb )
+                     z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp)    , rn_bb )
+               ELSE ! shallow water, uniform sigma
+                     z_gsigw3(ji,jj) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
+                     z_gsigw3p1(ji,jj) =   REAL(jk,wp)            / REAL(jpk-1,wp)
+                     z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp)
+               ENDIF
+               !
+             !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
+               IF( jk .EQ. 1) THEN
+                  z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) )
+                  z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
+               ELSE
+                  z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
+                  z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
+               ENDIF
+               IF( jk .EQ. jpk) THEN
+                  z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) )
+               ELSE
+                  z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
+               ENDIF
+               !
+
+               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
+               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
+               gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft )
+               gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw )
+               gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft )
+              !
+            END DO   ! for all jj's
+         END DO    ! for all ji's
+
+         DO ji = 1, jpim1
+            DO jj = 1, jpjm1
+               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj)     &
+                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               !
+               e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               !
+               e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+            END DO
+         END DO
+
+         z_gsigt3m1 = z_gsigt3
+         z_gsiw3m1 = z_gsiw3
+ 
+         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
+         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
+         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
+         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
+         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
+         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
+         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
+        
+         CALL write_netcdf_vars(jk)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
+                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
+            END DO
+         END DO
+      END DO !End of loop over jk
+
+      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
+      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+
+   END SUBROUTINE s_sh94
+
+   SUBROUTINE s_sf12
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_sf12 *** 
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
+      !!                mixed S/sigma/Z coordinate
+      !!
+      !!                This method allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !!
+      !!
+      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
+      !!----------------------------------------------------------------------
+      !
+      USE utils
+      REAL(wp) ::   zsmth               ! smoothing around critical depth
+      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
+      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
+      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
+      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
+      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
 
       z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
@@ -475,206 +477,108 @@
       z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
       z_esigwu3 = 0.   ;   z_esigwv3 = 0.
-
-      DO ji = 1, jpi
-         DO jj = 1, jpj
-
-            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, rn_bb )
-                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
-               END DO
-            ELSE ! shallow water, uniform sigma
-               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 ) / REAL(jpk-1,wp)
-                  END DO
-            ENDIF
-            !
-            DO jk = 1, jpk-1
-               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
-               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. * ( 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 * 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)
-            END DO
-            !
-            DO jk = 1, jpk
-               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
-           !
-         END DO   ! for all jj's
-      END DO    ! for all ji's
-
-      DO ji = 1, jpim1
-         DO jj = 1, jpjm1
-            DO jk = 1, jpk
-               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
-                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+      
+
+      DO jk = 1,jpk
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+             IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
+              
+                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-(zzb/hbatt(ji,jj))
+            
+                zzs = rn_zs / hbatt(ji,jj) 
+              
+                IF (rn_efold /= 0.0) THEN
+                   zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
+                ELSE
+                   zsmth = 1.0 
+                ENDIF
+                  
+                z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
+                z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)
+                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+                z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth  )
+                z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(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
+
+                z_gsigw3(ji,jj) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
+                z_gsigw3p1(ji,jj) =  REAL(jk,wp)     /REAL(jpk-1,wp)
+                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+
+             ELSE  ! shallow water, z coordinates
+
+               z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+               z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+               z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+
+             ENDIF
+
+             !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
+             IF( jk .EQ. 1) THEN
+                z_esigw3(ji,jj  ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ))
+                z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
+             ELSE
+                z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
+                z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
+             ENDIF
+             IF( jk .EQ. jpk) THEN
+                z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj))
+             ELSE
+                z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
+             ENDIF
+
+             gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj)
+             gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj)
+             gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj)
+
+            ENDDO   ! for all jj's
+         ENDDO    ! for all ji's
+
+         DO ji=1,jpi-1
+            DO jj=1,jpj-1
+
+               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) +  &
+                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+
+               e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj)
+               e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj)
+               e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj)
+               e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj)
                !
-               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-               !
-               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
-            END DO
-        END DO
-      END DO
-
-      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3)
-      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
-
-   END SUBROUTINE s_sh94
-
-   SUBROUTINE s_sf12
-
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE s_sf12 *** 
-      !!                     
-      !! ** Purpose :   stretch the s-coordinate system
-      !!
-      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
-      !!                mixed S/sigma/Z coordinate
-      !!
-      !!                This method allows the maintenance of fixed surface and or
-      !!                bottom cell resolutions (cf. geopotential coordinates) 
-      !!                within an analytically derived stretched S-coordinate framework.
-      !!
-      !!
-      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
-      !!----------------------------------------------------------------------
-      !
-      USE utils
-      REAL(wp) ::   zsmth               ! smoothing around critical depth
-      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
-      !
-      ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) )
-      ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 
-      ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 
-      ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
-
-      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
-         DO jj = 1, jpj
-
-          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
-              
-              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-(zzb/hbatt(ji,jj))
-            
-              zzs = rn_zs / hbatt(ji,jj) 
-              
-              IF (rn_efold /= 0.0) THEN
-                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
-              ELSE
-                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)/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  )
+               e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj)
+               e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj)
+               e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj)
+            ENDDO
+         ENDDO
+         ! Keep some arrays for next level
+         z_gsigt3m1 = z_gsigt3
+         z_gsiw3m1 = z_gsiw3
  
-          ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma
-
-            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)/REAL(jpk-1,wp)
-            END DO
-
-          ELSE  ! shallow water, z coordinates
-
-            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)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
-            END DO
-
-          ENDIF
-
-          DO jk = 1, jpk-1
-             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
-             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 * (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
-          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)
-          END DO
-
-          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)
-             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
-          END DO
-
-        ENDDO   ! for all jj's
-      ENDDO    ! for all ji's
-
-      DO ji=1,jpi-1
-        DO jj=1,jpj-1
-
-          DO jk = 1, jpk
-                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
-                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-
-             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
-             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
-             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
-             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
-             !
-             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
-             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
-             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
-          END DO
-
-        ENDDO
-      ENDDO
-      !
-      !                                               ! =============
-
-      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3)
+         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
+         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
+         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
+         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
+         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
+         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
+         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
+         
+         WRITE(*,*) 'Writing level ',jk,' to file'
+         CALL write_netcdf_vars(jk)
+         WRITE(*,*) 'Written level ',jk,' to file'
+      ENDDO ! End of loop over jk
+
+      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
       DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
 
@@ -728,25 +632,37 @@
 !!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 ) / 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
 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            DO jk = 1, jpk
-              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
-              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
-              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
-              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+              zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
+              zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
+              gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft )
+              gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw )
+              gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft )
+              e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
+              e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
+              e3f_0(ji,jj) = ( (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) )
-              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) )
+              e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
+              e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
+              IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
+              IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
             END DO
          END DO
-      END DO
+         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
+         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
+         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
+         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
+         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
+         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
+         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
+ 
+         CALL write_netcdf_vars(jk)
+      ENDDO ! End of loop over jk
+
 
       DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
@@ -826,12 +742,12 @@
       !!----------------------------------------------------------------------
       USE utils, ONLY : jpk,jk,wp
-      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
-      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
+      REAL(wp), INTENT(in   ) ::   pk1           ! continuous "k" coordinate
+      REAL(wp)                ::   p_gamma       ! stretched coordinate
       REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
       REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
-      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
-      REAL(wp)                ::   za1,za2,za3    ! local variables
-      REAL(wp)                ::   zn1,zn2        ! local variables
-      REAL(wp)                ::   za,zb,zx       ! local variables
+      REAL(wp), INTENT(in   ) ::   psmth         ! Smoothing parameter
+      REAL(wp)                ::   za1,za2,za3   ! local variables
+      REAL(wp)                ::   zn1,zn2       ! local variables
+      REAL(wp)                ::   za,zb,zx      ! local variables
       !!----------------------------------------------------------------------
       !
@@ -849,10 +765,8 @@
       zx = 1.0-za/2.0-zb
  
-      DO jk = 1, jpk
-        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 
+      p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 +  &
+                  & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- &
+                  &      (rn_alpha+1.0)*pk1**(rn_alpha+2.0) )
+      p_gamma = p_gamma*psmth+pk1*(1.0-psmth)
 
       !
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 5268)
+++ branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 5269)
@@ -13,10 +13,10 @@
    !! All coordinates
    !! ---------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
 
    !! s-coordinate and hybrid z-s-coordinate
@@ -45,5 +45,5 @@
    INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
    INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation
-   INTEGER  ::   numnam                   ! File handle for namelist 
+   INTEGER,PARAMETER  ::   numnam=8       ! File handle for namelist 
    REAL(wp) ::   zrmax, ztaper   ! temporary scalars
    REAL(wp) ::   zrfact
@@ -60,4 +60,9 @@
                         rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
 
+  ! IDs for output netcdf file
+  INTEGER :: id_x, id_y, id_z
+  INTEGER :: ncout
+  INTEGER, DIMENSION(11) :: var_ids  !Array to contain all variable IDs
+
    CONTAINS
 
@@ -71,7 +76,7 @@
          &      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(2) )
+      ALLOCATE( gdep3w_0(jpi,jpj) , e3v_0(jpi,jpj) , e3f_0(jpi,jpj) ,                         &
+         &      gdept_0(jpi,jpj) , e3t_0(jpi,jpj) , e3u_0 (jpi,jpj) ,                         &
+         &      gdepw_0(jpi,jpj) , e3w_0(jpi,jpj) , e3vw_0(jpi,jpj) , e3uw_0(jpi,jpj) , STAT=ierr(2) )
          !
          !
@@ -97,11 +102,11 @@
 
      ! Find the size of the input bathymetry
-     CALL dimlen(ncin, 'x', jpi)    
-     CALL dimlen(ncin, 'y', jpj)    
+     CALL dimlen(ncin, 'lon', jpi)    
+     CALL dimlen(ncin, 'lat', jpj)    
      
      ALLOCATE( bathy(jpi, jpj) )
      
      ! Read the bathymetry variable from file
-     CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', var_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 /) ) )
 
@@ -125,10 +130,8 @@
   
  
-   SUBROUTINE write_coord_file()
-     ! Write out variables to the a netcdf coordinates file
+   SUBROUTINE make_coord_file()
+     ! Create new coordinates file and define dimensions and variables ready for
+     ! writing
      
-     INTEGER :: id_x, id_y, id_z
-     INTEGER :: ncout
-     INTEGER, DIMENSION(11) :: var_ids  !Array to contain all variable IDs
 
      !Create the file
@@ -141,35 +144,39 @@
      !
      !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)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_z/), var_ids(1)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(2)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(3)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_z/), var_ids(4)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_z/), var_ids(5)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_z/), var_ids(6)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_z/), var_ids(7)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(8)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(9)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(10)) )
+     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y/), var_ids(11)) )
      
      ! End define mode
      CALL check_nf90( nf90_enddef(ncout) )
+     
+     WRITE(*,*) 'Opened coord_zgr.nc file and defined variables'
 
-     !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 make_coord_file
 
-   END SUBROUTINE write_coord_file
+   SUBROUTINE write_netcdf_vars(kk)
+   ! Write  variables to the netcdf file at level kk
+     INTEGER, INTENT(in) :: kk
+
+     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+
+   END SUBROUTINE write_netcdf_vars
 
    SUBROUTINE check_nf90( istat, message )
