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

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

File:
1 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      ! 
Note: See TracChangeset for help on using the changeset viewer.