Changeset 5269


Ignore:
Timestamp:
2015-05-13T15:28:30+02:00 (5 years ago)
Author:
timgraham
Message:

Removed use of 3D variables.
Compiles and runs with no errors but still need to check output

Location:
branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90

    r5257 r5269  
    286286         WRITE(*,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  & 
    287287            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 
    288 !! helsinki 
     288       
     289      ! Create the output file 
     290      CALL make_coord_file() 
    289291 
    290292      !                                            ! ======================= 
     
    293295      ! 
    294296      ! non-dimensional "sigma" for model level depth at w- and t-levels 
     297 
    295298 
    296299!======================================================================== 
     
    298301! Siddorn and Furner 2012 (ln_sf12=T) 
    299302! or  tanh function       (both false)                     
     303! To reduce memory loop over jpk and write each level to file 
    300304!======================================================================== 
    301305      IF      ( ln_s_sh94 ) THEN  
     
    305309      ELSE                  
    306310                           CALL s_tanh() 
    307       ENDIF  
    308  
    309       ! 
    310       where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    311       where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    312       where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    313       where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    314       where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    315       where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    316       where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
    317  
    318 #if defined key_agrif 
    319       ! Ensure meaningful vertical scale factors in ghost lines/columns 
    320       ! TODO: How can we deal with this in offline tool?  
    321       IF( .NOT. Agrif_Root() ) THEN 
    322          !   
    323          IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    324             e3u_0(1,:,:) = e3u_0(2,:,:) 
    325          ENDIF 
    326          ! 
    327          IF((nbondi ==  1).OR.(nbondi == 2)) THEN 
    328             e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 
    329          ENDIF 
    330          ! 
    331          IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    332             e3v_0(:,1,:) = e3v_0(:,2,:) 
    333          ENDIF 
    334          ! 
    335          IF((nbondj ==  1).OR.(nbondj == 2)) THEN 
    336             e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 
    337          ENDIF 
    338          ! 
    339311      ENDIF 
    340 #endif 
    341 !! 
    342       ! HYBRID :  
    343       DO jj = 1, jpj 
    344          DO ji = 1, jpi 
    345             DO jk = 1, jpk-1 
    346                IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    347             END DO 
    348             IF( scobot(ji,jj) == 0.               )   mbathy(ji,jj) = 0 
    349          END DO 
    350       END DO 
    351       WRITE(*,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    352  
    353       ! min max values over the domain 
    354       WRITE(*,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    355       WRITE(*,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
    356          &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) ) 
    357       WRITE(*,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   & 
    358          &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   & 
    359          &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   & 
    360          &                          ' w ', MINVAL( e3w_0  (:,:,:) ) 
    361  
    362       WRITE(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   & 
    363          &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) ) 
    364       WRITE(*,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   & 
    365          &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   & 
    366          &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   & 
    367          &                          ' w ', MAXVAL( e3w_0  (:,:,:) ) 
    368        
    369          ! selected vertical profiles 
    370       WRITE(*,*) 
    371       WRITE(*,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 
    372       WRITE(*,*) ' ~~~~~~  --------------------' 
    373       WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    374       WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     & 
    375          &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 
    376       DO jj = 20, 20 
    377          DO ji = 20, 20 
    378             WRITE(*,*) 
    379             WRITE(*,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    380             WRITE(*,*) ' ~~~~~~  --------------------' 
    381             WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    382             WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
    383                &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
    384          END DO 
    385       END DO 
    386       DO jj = 74,74 
    387          DO ji = 100, 100 
    388             WRITE(*,*) 
    389             WRITE(*,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    390             WRITE(*,*) ' ~~~~~~  --------------------' 
    391             WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    392             WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
    393                &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
    394          END DO 
    395       END DO 
    396  
    397 !================================================================================ 
    398 ! check the coordinate makes sense 
    399 !================================================================================ 
    400       DO ji = 1, jpi 
    401          DO jj = 1, jpj 
    402  
    403             IF( hbatt(ji,jj) > 0.) THEN 
    404                DO jk = 1, mbathy(ji,jj) 
    405                  ! check coordinate is monotonically increasing 
    406                  IF (e3w_0(ji,jj,jk) <= 0. .OR. e3t_0(ji,jj,jk) <= 0. ) THEN 
    407                     WRITE(*,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    408                     WRITE(*,*) 'e3w',e3w_0(ji,jj,:) 
    409                     WRITE(*,*) 'e3t',e3t_0(ji,jj,:) 
    410                     STOP 1 
    411                  ENDIF 
    412                  ! and check it has never gone negative 
    413                  IF( gdepw_0(ji,jj,jk) < 0. .OR. gdept_0(ji,jj,jk) < 0. ) THEN 
    414                     WRITE(*,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
    415                     WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) 
    416                     WRITE(*,*) 'gdept',gdept_0(ji,jj,:) 
    417                     STOP 1 
    418                  ENDIF 
    419                  ! and check it never exceeds the total depth 
    420                  IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    421                     WRITE(*,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    422                     WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) 
    423                     STOP 1 
    424                  ENDIF 
    425                END DO 
    426  
    427                DO jk = 1, mbathy(ji,jj)-1 
    428                  ! and check it never exceeds the total depth 
    429                 IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    430                     WRITE(*,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    431                     WRITE(*,*) 'gdept',gdept_0(ji,jj,:) 
    432                     STOP 1 
    433                  ENDIF 
    434                END DO 
    435  
    436             ENDIF 
    437  
    438          END DO 
    439       END DO 
    440       ! 
    441       ! Write output file 
    442       CALL write_coord_file() 
    443       ! 
     312 
     313      CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) 
     314      CALL check_nf90( nf90_close(ncout) ) 
     315 
     316  
    444317      ! 
    445318END PROGRAM SCOORD_GEN 
     
    462335      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
    463336      ! 
    464       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
    465       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3  
    466       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
    467  
    468       ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) ) 
    469       ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk))  
    470       ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk))  
    471       ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 
     337      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     338      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 
     339      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3  
     340      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
     341 
     342      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) 
     343      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj))  
     344      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj))  
     345      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) 
     346      ALLOCATE( z_gsi3w3m1(jpi,jpj) ) 
     347  
     348      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0. 
     349      z_esigt3  = 0.   ;   z_esigw3  = 0.  
     350      z_esigt3p1  = 0. ;   z_esigw3p1  = 0.  
     351      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0. 
     352      z_esigwu3 = 0.   ;   z_esigwv3 = 0. 
     353 
     354      DO jk = 1,jpk 
     355         DO ji = 1, jpi 
     356            DO jj = 1, jpj 
     357 
     358               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
     359                     z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) 
     360                     z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb ) 
     361                     z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp)    , rn_bb ) 
     362               ELSE ! shallow water, uniform sigma 
     363                     z_gsigw3(ji,jj) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
     364                     z_gsigw3p1(ji,jj) =   REAL(jk,wp)            / REAL(jpk-1,wp) 
     365                     z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) 
     366               ENDIF 
     367               ! 
     368             !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk 
     369               IF( jk .EQ. 1) THEN 
     370                  z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) ) 
     371                  z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) 
     372               ELSE 
     373                  z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) 
     374                  z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) 
     375               ENDIF 
     376               IF( jk .EQ. jpk) THEN 
     377                  z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) ) 
     378               ELSE 
     379                  z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) 
     380               ENDIF 
     381               ! 
     382 
     383               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 
     384               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 
     385               gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft ) 
     386               gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw ) 
     387               gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft ) 
     388              ! 
     389            END DO   ! for all jj's 
     390         END DO    ! for all ji's 
     391 
     392         DO ji = 1, jpim1 
     393            DO jj = 1, jpjm1 
     394               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) )   & 
     395                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     396               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) )   & 
     397                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     398               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj)     & 
     399                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) )   & 
     400                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     401               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) )   & 
     402                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     403               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) )   & 
     404                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     405               ! 
     406               e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     407               e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     408               e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     409               e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     410               ! 
     411               e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     412               e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     413               e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 
     414            END DO 
     415         END DO 
     416 
     417         z_gsigt3m1 = z_gsigt3 
     418         z_gsiw3m1 = z_gsiw3 
     419  
     420         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0 
     421         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0 
     422         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0 
     423         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0 
     424         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0 
     425         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0 
     426         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0 
     427         
     428         CALL write_netcdf_vars(jk) 
     429         DO jj = 1, jpj 
     430            DO ji = 1, jpi 
     431                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk ) 
     432                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0 
     433            END DO 
     434         END DO 
     435      END DO !End of loop over jk 
     436 
     437      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) 
     438      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
     439 
     440   END SUBROUTINE s_sh94 
     441 
     442   SUBROUTINE s_sf12 
     443 
     444      !!---------------------------------------------------------------------- 
     445      !!                  ***  ROUTINE s_sf12 ***  
     446      !!                      
     447      !! ** Purpose :   stretch the s-coordinate system 
     448      !! 
     449      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
     450      !!                mixed S/sigma/Z coordinate 
     451      !! 
     452      !!                This method allows the maintenance of fixed surface and or 
     453      !!                bottom cell resolutions (cf. geopotential coordinates)  
     454      !!                within an analytically derived stretched S-coordinate framework. 
     455      !! 
     456      !! 
     457      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 
     458      !!---------------------------------------------------------------------- 
     459      ! 
     460      USE utils 
     461      REAL(wp) ::   zsmth               ! smoothing around critical depth 
     462      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     463      ! 
     464      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     465      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 
     466      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3  
     467      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
     468 
     469      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) 
     470      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj))  
     471      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj))  
     472      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) 
     473      ALLOCATE( z_gsi3w3m1(jpi,jpj) ) 
    472474 
    473475      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0. 
     
    475477      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0. 
    476478      z_esigwu3 = 0.   ;   z_esigwv3 = 0. 
    477  
    478       DO ji = 1, jpi 
    479          DO jj = 1, jpj 
    480  
    481             IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
    482                DO jk = 1, jpk 
    483                   z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) 
    484                   z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
    485                END DO 
    486             ELSE ! shallow water, uniform sigma 
    487                DO jk = 1, jpk 
    488                   z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
    489                   z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) 
    490                   END DO 
    491             ENDIF 
    492             ! 
    493             DO jk = 1, jpk-1 
    494                z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 
    495                z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 
    496             END DO 
    497             z_esigw3(ji,jj,1  ) = 2. * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) ) 
    498             z_esigt3(ji,jj,jpk) = 2. * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 
    499             ! 
    500             ! Coefficients for vertical depth as the sum of e3w scale factors 
    501             z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 
    502             DO jk = 2, jpk 
    503                z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
    504             END DO 
    505             ! 
    506             DO jk = 1, jpk 
    507                zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 
    508                zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 
    509                gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    510                gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    511                gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    512             END DO 
    513            ! 
    514          END DO   ! for all jj's 
    515       END DO    ! for all ji's 
    516  
    517       DO ji = 1, jpim1 
    518          DO jj = 1, jpjm1 
    519             DO jk = 1, jpk 
    520                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
    521                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    522                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    523                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    524                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
    525                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    526                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    527                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    528                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    529                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    530                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     479       
     480 
     481      DO jk = 1,jpk 
     482         DO ji = 1, jpi 
     483            DO jj = 1, jpj 
     484             IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 
     485               
     486                zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
     487                                                     ! could be changed by users but care must be taken to do so carefully 
     488                zzb = 1.0-(zzb/hbatt(ji,jj)) 
     489             
     490                zzs = rn_zs / hbatt(ji,jj)  
     491               
     492                IF (rn_efold /= 0.0) THEN 
     493                   zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
     494                ELSE 
     495                   zsmth = 1.0  
     496                ENDIF 
     497                   
     498                z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp) 
     499                z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp) 
     500                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
     501                z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth  ) 
     502                z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth  ) 
     503                z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth  ) 
     504  
     505             ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma 
     506 
     507                z_gsigw3(ji,jj) =  REAL(jk-1,wp)     /REAL(jpk-1,wp) 
     508                z_gsigw3p1(ji,jj) =  REAL(jk,wp)     /REAL(jpk-1,wp) 
     509                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
     510 
     511             ELSE  ! shallow water, z coordinates 
     512 
     513               z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     514               z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     515               z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     516 
     517             ENDIF 
     518 
     519             !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk 
     520             IF( jk .EQ. 1) THEN 
     521                z_esigw3(ji,jj  ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj )) 
     522                z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) 
     523             ELSE 
     524                z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) 
     525                z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) 
     526             ENDIF 
     527             IF( jk .EQ. jpk) THEN 
     528                z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj)) 
     529             ELSE 
     530                z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) 
     531             ENDIF 
     532 
     533             gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj) 
     534             gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj) 
     535             gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj) 
     536 
     537            ENDDO   ! for all jj's 
     538         ENDDO    ! for all ji's 
     539 
     540         DO ji=1,jpi-1 
     541            DO jj=1,jpj-1 
     542 
     543               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / & 
     544                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     545               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / & 
     546                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     547               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) +  & 
     548                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / & 
     549                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     550               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / & 
     551                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     552               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / & 
     553                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     554 
     555               e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj) 
     556               e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj) 
     557               e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj) 
     558               e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj) 
    531559               ! 
    532                e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    533                e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    534                e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    535                e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    536                ! 
    537                e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    538                e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    539                e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 
    540             END DO 
    541         END DO 
    542       END DO 
    543  
    544       DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3) 
    545       DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    546  
    547    END SUBROUTINE s_sh94 
    548  
    549    SUBROUTINE s_sf12 
    550  
    551       !!---------------------------------------------------------------------- 
    552       !!                  ***  ROUTINE s_sf12 ***  
    553       !!                      
    554       !! ** Purpose :   stretch the s-coordinate system 
    555       !! 
    556       !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
    557       !!                mixed S/sigma/Z coordinate 
    558       !! 
    559       !!                This method allows the maintenance of fixed surface and or 
    560       !!                bottom cell resolutions (cf. geopotential coordinates)  
    561       !!                within an analytically derived stretched S-coordinate framework. 
    562       !! 
    563       !! 
    564       !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 
    565       !!---------------------------------------------------------------------- 
    566       ! 
    567       USE utils 
    568       REAL(wp) ::   zsmth               ! smoothing around critical depth 
    569       REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
    570       ! 
    571       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
    572       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3  
    573       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
    574       ! 
    575       ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) ) 
    576       ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk))  
    577       ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk))  
    578       ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 
    579  
    580       z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0. 
    581       z_esigt3  = 0.   ;   z_esigw3  = 0.  
    582       z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0. 
    583       z_esigwu3 = 0.   ;   z_esigwv3 = 0. 
    584  
    585       DO ji = 1, jpi 
    586          DO jj = 1, jpj 
    587  
    588           IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 
    589                
    590               zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
    591                                                      ! could be changed by users but care must be taken to do so carefully 
    592               zzb = 1.0-(zzb/hbatt(ji,jj)) 
    593              
    594               zzs = rn_zs / hbatt(ji,jj)  
    595                
    596               IF (rn_efold /= 0.0) THEN 
    597                 zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
    598               ELSE 
    599                 zsmth = 1.0  
    600               ENDIF 
    601                 
    602               DO jk = 1, jpk 
    603                 z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp) 
    604                 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
    605               ENDDO 
    606               z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  ) 
    607               z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  ) 
     560               e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj) 
     561               e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj) 
     562               e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj) 
     563            ENDDO 
     564         ENDDO 
     565         ! Keep some arrays for next level 
     566         z_gsigt3m1 = z_gsigt3 
     567         z_gsiw3m1 = z_gsiw3 
    608568  
    609           ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma 
    610  
    611             DO jk = 1, jpk 
    612               z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp) 
    613               z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
    614             END DO 
    615  
    616           ELSE  ! shallow water, z coordinates 
    617  
    618             DO jk = 1, jpk 
    619               z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
    620               z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
    621             END DO 
    622  
    623           ENDIF 
    624  
    625           DO jk = 1, jpk-1 
    626              z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 
    627              z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 
    628           END DO 
    629           z_esigw3(ji,jj,1  ) = 2.0 * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  )) 
    630           z_esigt3(ji,jj,jpk) = 2.0 * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 
    631  
    632           ! Coefficients for vertical depth as the sum of e3w scale factors 
    633           z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 
    634           DO jk = 2, jpk 
    635              z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
    636           END DO 
    637  
    638           DO jk = 1, jpk 
    639              gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 
    640              gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 
    641              gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 
    642           END DO 
    643  
    644         ENDDO   ! for all jj's 
    645       ENDDO    ! for all ji's 
    646  
    647       DO ji=1,jpi-1 
    648         DO jj=1,jpj-1 
    649  
    650           DO jk = 1, jpk 
    651                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
    652                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    653                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    654                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    655                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
    656                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    657                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    658                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    659                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    660                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    661                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    662  
    663              e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
    664              e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 
    665              e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 
    666              e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 
    667              ! 
    668              e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 
    669              e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 
    670              e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 
    671           END DO 
    672  
    673         ENDDO 
    674       ENDDO 
    675       ! 
    676       !                                               ! ============= 
    677  
    678       DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3) 
     569         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0 
     570         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0 
     571         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0 
     572         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0 
     573         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0 
     574         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0 
     575         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0 
     576          
     577         WRITE(*,*) 'Writing level ',jk,' to file' 
     578         CALL write_netcdf_vars(jk) 
     579         WRITE(*,*) 'Written level ',jk,' to file' 
     580      ENDDO ! End of loop over jk 
     581 
     582      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) 
    679583      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    680584 
     
    728632!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    729633      DO jk = 1, jpk 
    730          zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 
    731          zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 
    732          gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 
    733          gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 
    734          gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 
    735634      END DO 
    736635!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
    737       DO jj = 1, jpj 
    738          DO ji = 1, jpi 
    739             DO jk = 1, jpk 
    740               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 
    741               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 
    742               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 
    743               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) 
     636      DO jk = 1, jpk 
     637         DO jj = 1, jpj 
     638            DO ji = 1, jpi 
     639              zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 
     640              zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 
     641              gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft ) 
     642              gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw ) 
     643              gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft ) 
     644              e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 
     645              e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 
     646              e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 
     647              e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) 
    744648              ! 
    745               e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 
    746               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 
    747               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 
     649              e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 
     650              e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 
     651              e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 
     652              IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk ) 
     653              IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0 
    748654            END DO 
    749655         END DO 
    750       END DO 
     656         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0 
     657         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0 
     658         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0 
     659         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0 
     660         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0 
     661         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0 
     662         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0 
     663  
     664         CALL write_netcdf_vars(jk) 
     665      ENDDO ! End of loop over jk 
     666 
    751667 
    752668      DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w ) 
     
    826742      !!---------------------------------------------------------------------- 
    827743      USE utils, ONLY : jpk,jk,wp 
    828       REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate 
    829       REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate 
     744      REAL(wp), INTENT(in   ) ::   pk1           ! continuous "k" coordinate 
     745      REAL(wp)                ::   p_gamma       ! stretched coordinate 
    830746      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth 
    831747      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth 
    832       REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter 
    833       REAL(wp)                ::   za1,za2,za3    ! local variables 
    834       REAL(wp)                ::   zn1,zn2        ! local variables 
    835       REAL(wp)                ::   za,zb,zx       ! local variables 
     748      REAL(wp), INTENT(in   ) ::   psmth         ! Smoothing parameter 
     749      REAL(wp)                ::   za1,za2,za3   ! local variables 
     750      REAL(wp)                ::   zn1,zn2       ! local variables 
     751      REAL(wp)                ::   za,zb,zx      ! local variables 
    836752      !!---------------------------------------------------------------------- 
    837753      ! 
     
    849765      zx = 1.0-za/2.0-zb 
    850766  
    851       DO jk = 1, jpk 
    852         p_gamma(jk) = za*(pk1(jk)*(1.0-pk1(jk)/2.0))+zb*pk1(jk)**3.0 +  & 
    853                     & zx*( (rn_alpha+2.0)*pk1(jk)**(rn_alpha+1.0)- & 
    854                     &      (rn_alpha+1.0)*pk1(jk)**(rn_alpha+2.0) ) 
    855         p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0-psmth) 
    856       ENDDO  
     767      p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 +  & 
     768                  & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- & 
     769                  &      (rn_alpha+1.0)*pk1**(rn_alpha+2.0) ) 
     770      p_gamma = p_gamma*psmth+pk1*(1.0-psmth) 
    857771 
    858772      ! 
  • branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90

    r5257 r5269  
    1313   !! All coordinates 
    1414   !! --------------- 
    15    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m) 
    16    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m) 
    17    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f 
    18    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m) 
    19    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw 
    20    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m) 
     15   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m) 
     16   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m) 
     17   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f 
     18   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m) 
     19   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw 
     20   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m) 
    2121 
    2222   !! s-coordinate and hybrid z-s-coordinate 
     
    4545   INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    4646   INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation 
    47    INTEGER  ::   numnam                   ! File handle for namelist  
     47   INTEGER,PARAMETER  ::   numnam=8       ! File handle for namelist  
    4848   REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    4949   REAL(wp) ::   zrfact 
     
    6060                        rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 
    6161 
     62  ! IDs for output netcdf file 
     63  INTEGER :: id_x, id_y, id_z 
     64  INTEGER :: ncout 
     65  INTEGER, DIMENSION(11) :: var_ids  !Array to contain all variable IDs 
     66 
    6267   CONTAINS 
    6368 
     
    7176         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) ) 
    7277         ! 
    73       ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) ,                         & 
    74          &      gdept_0(jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         & 
    75          &      gdepw_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(2) ) 
     78      ALLOCATE( gdep3w_0(jpi,jpj) , e3v_0(jpi,jpj) , e3f_0(jpi,jpj) ,                         & 
     79         &      gdept_0(jpi,jpj) , e3t_0(jpi,jpj) , e3u_0 (jpi,jpj) ,                         & 
     80         &      gdepw_0(jpi,jpj) , e3w_0(jpi,jpj) , e3vw_0(jpi,jpj) , e3uw_0(jpi,jpj) , STAT=ierr(2) ) 
    7681         ! 
    7782         ! 
     
    97102 
    98103     ! Find the size of the input bathymetry 
    99      CALL dimlen(ncin, 'x', jpi)     
    100      CALL dimlen(ncin, 'y', jpj)     
     104     CALL dimlen(ncin, 'lon', jpi)     
     105     CALL dimlen(ncin, 'lat', jpj)     
    101106      
    102107     ALLOCATE( bathy(jpi, jpj) ) 
    103108      
    104109     ! Read the bathymetry variable from file 
    105      CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', var_id ), 'Cannot get variable ID for bathymetry') 
     110     CALL check_nf90( nf90_inq_varid( ncin, 'Bathymetry', var_id ), 'Cannot get variable ID for bathymetry') 
    106111     CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) ) 
    107112 
     
    125130   
    126131  
    127    SUBROUTINE write_coord_file() 
    128      ! Write out variables to the a netcdf coordinates file 
     132   SUBROUTINE make_coord_file() 
     133     ! Create new coordinates file and define dimensions and variables ready for 
     134     ! writing 
    129135      
    130      INTEGER :: id_x, id_y, id_z 
    131      INTEGER :: ncout 
    132      INTEGER, DIMENSION(11) :: var_ids  !Array to contain all variable IDs 
    133136 
    134137     !Create the file 
     
    141144     ! 
    142145     !Define variables 
    143      CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_x/), var_ids(1)) ) 
    144      CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(2)) ) 
    145      CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(3)) ) 
    146      CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_x/), var_ids(4)) ) 
    147      CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_x/), var_ids(5)) ) 
    148      CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_x/), var_ids(6)) ) 
    149      CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_x/), var_ids(7)) ) 
    150      CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(8)) ) 
    151      CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(9)) ) 
    152      CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(10)) ) 
    153      CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y,id_x/), var_ids(11)) ) 
     146     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_z/), var_ids(1)) ) 
     147     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(2)) ) 
     148     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(3)) ) 
     149     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_z/), var_ids(4)) ) 
     150     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_z/), var_ids(5)) ) 
     151     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_z/), var_ids(6)) ) 
     152     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_z/), var_ids(7)) ) 
     153     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(8)) ) 
     154     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(9)) ) 
     155     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(10)) ) 
     156     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y/), var_ids(11)) ) 
    154157      
    155158     ! End define mode 
    156159     CALL check_nf90( nf90_enddef(ncout) ) 
     160      
     161     WRITE(*,*) 'Opened coord_zgr.nc file and defined variables' 
    157162 
    158      !Write variables to file 
    159      CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0) ) 
    160      CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0) ) 
    161      CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0) ) 
    162      CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0) ) 
    163      CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0) ) 
    164      CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0) ) 
    165      CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0) ) 
    166      CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0) ) 
    167      CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0) ) 
    168      CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0) ) 
    169      CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) 
    170       
    171      CALL check_nf90( nf90_close(ncout) ) 
     163   END SUBROUTINE make_coord_file 
    172164 
    173    END SUBROUTINE write_coord_file 
     165   SUBROUTINE write_netcdf_vars(kk) 
     166   ! Write  variables to the netcdf file at level kk 
     167     INTEGER, INTENT(in) :: kk 
     168 
     169     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     170     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     171     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     172     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     173     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     174     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     175     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     176     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     177     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     178     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 
     179 
     180   END SUBROUTINE write_netcdf_vars 
    174181 
    175182   SUBROUTINE check_nf90( istat, message ) 
Note: See TracChangeset for help on using the changeset viewer.