Changeset 1099


Ignore:
Timestamp:
2008-06-09T10:58:04+02:00 (12 years ago)
Author:
ctlod
Message:

trunk: code review stuff, see tickets: #205, #191, #130

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/DOM/domzgr.F90

    r1083 r1099  
    44   !! Ocean initialization : domain initialization 
    55   !!============================================================================== 
    6  
     6   !! History :  OPA  !  1995-12  (G. Madec)  Original code : s vertical coordinate 
     7   !!                 !  1997-07  (G. Madec)  lbc_lnk call 
     8   !!                 !  1997-04  (J.-O. Beismann)  
     9   !!            8.5  !  2002-09 (A. Bozec, G. Madec)  F90: Free form and module 
     10   !!             -   !  2002-09 (A. de Miranda)  rigid-lid + islands 
     11   !!  NEMO      1.0  !  2003-08  (G. Madec)  F90: Free form and module 
     12   !!             -   !  2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function 
     13   !!            2.0  !  2006-04  (R. Benshila, G. Madec)  add zgr_zco 
     14   !!            3.0  !  2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style 
    715   !!---------------------------------------------------------------------- 
    8    !!   dom_zgr     : defined the ocean vertical coordinate system 
     16 
     17   !!---------------------------------------------------------------------- 
     18   !!   dom_zgr          : defined the ocean vertical coordinate system 
    919   !!       zgr_bat      : bathymetry fields (levels and meters) 
    1020   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain 
     
    1525   !!       zgr_sco      : s-coordinate 
    1626   !!       fssig        : sigma coordinate non-dimensional function 
     27   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!) 
    1728   !!--------------------------------------------------------------------- 
    18    !! * Modules used 
    1929   USE oce             ! ocean dynamics and tracers 
    2030   USE dom_oce         ! ocean space and time domain 
     
    2232   USE lib_mpp         ! distributed memory computing library 
    2333   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    24    USE closea 
    25    USE solisl 
     34   USE closea          ! closed seas 
     35   USE solisl          ! solver - island in rigid-lid 
    2636   USE c1d             ! 1D configuration 
    2737 
     
    2939   PRIVATE 
    3040 
    31    !! * Routine accessibility 
    32    PUBLIC dom_zgr        ! called by dom_init.F90 
    33  
    34    !! * Module variables 
    35       REAL(wp) ::          & !!: Namelist nam_zgr_sco 
    36       sbot_min =  300.  ,  &  !: minimum depth of s-bottom surface (>0) (m) 
    37       sbot_max = 5250.  ,  &  !: maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    38       theta    =    6.0 ,  &  !: surface control parameter (0<=theta<=20) 
    39       thetb    =    0.75,  &  !: bottom control parameter  (0<=thetb<= 1) 
    40       r_max    =    0.15      !: maximum cut-off r-value allowed (0<r_max<1) 
    41  
    42  
     41   PUBLIC   dom_zgr    ! called by dom_init.F90 
     42 
     43!!gm   DOCTOR name for the namelist parameter... 
     44   !                                 !!! ** Namelist nam_zgr_sco ** 
     45   REAL(wp) ::   sbot_min =  300.     ! minimum depth of s-bottom surface (>0) (m) 
     46   REAL(wp) ::   sbot_max = 5250.     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     47   REAL(wp) ::   theta    =    6.0    ! surface control parameter (0<=theta<=20) 
     48   REAL(wp) ::   thetb    =    0.75   ! bottom control parameter  (0<=thetb<= 1) 
     49   REAL(wp) ::   r_max    =    0.15   ! maximum cut-off r-value allowed (0<r_max<1) 
    4350  
    4451   !! * Substitutions 
     
    4653#  include "vectopt_loop_substitute.h90" 
    4754   !!---------------------------------------------------------------------- 
    48    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     55   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     56   !! $Id:$ 
     57   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4958   !!---------------------------------------------------------------------- 
    5059 
     
    5867      !!      vertical scale factors. 
    5968      !! 
    60       !! ** Method  reference vertical coordinate 
    61       !! 
    62       !! ** Action :  
    63       !! 
    64       !! History : 
    65       !!   9.0  !  03-08  (G. Madec)  original code 
    66       !!   9.0  !  05-10  (A. Beckmann)  modifications for hybrid s-ccordinates 
    67       !!---------------------------------------------------------------------- 
    68       INTEGER ::   ioptio = 0      ! temporary integer 
    69  
     69      !! ** Method  : - reference 1D vertical coordinate (gdep._0, e3._0) 
     70      !!              - read/set ocean depth and ocean levels (bathy, mbathy) 
     71      !!              - vertical coordinate (gdep., e3.) depending on the  
     72      !!                coordinate chosen : 
     73      !!                   ln_zco=T   z-coordinate   (forced if lk_zco) 
     74      !!                   ln_zps=T   z-coordinate with partial steps 
     75      !!                   ln_zco=T   s-coordinate  
     76      !! 
     77      !! ** Action  :   define gdep., e3., mbathy and bathy 
     78      !!---------------------------------------------------------------------- 
     79      INTEGER ::   ioptio = 0   ! temporary integer 
     80      !! 
    7081      NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco 
    7182      !!---------------------------------------------------------------------- 
    7283 
    73       ! Read Namelist nam_zgr : vertical coordinate' 
    74       ! --------------------- 
    75       REWIND ( numnam ) 
     84      REWIND ( numnam )                ! Read Namelist nam_zgr : vertical coordinate' 
    7685      READ   ( numnam, nam_zgr ) 
    7786 
    78       ! Parameter control and print 
    79       ! --------------------------- 
    80       ! Control print 
    81       IF(lwp) THEN 
     87      IF(lwp) THEN                     ! Control print 
    8288         WRITE(numout,*) 
    8389         WRITE(numout,*) 'dom_zgr : vertical coordinate' 
     
    8995      ENDIF 
    9096 
    91       ! Check Vertical coordinate options 
    92       ioptio = 0 
     97      ioptio = 0                       ! Check Vertical coordinate options 
    9398      IF( ln_zco ) ioptio = ioptio + 1 
    9499      IF( ln_zps ) ioptio = ioptio + 1 
    95100      IF( ln_sco ) ioptio = ioptio + 1 
    96       IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 
    97  
     101      IF ( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' ) 
    98102      IF( lk_zco ) THEN 
    99103          IF(lwp) WRITE(numout,*) '          z-coordinate with reduced incore memory requirement' 
    100           IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' ) 
     104          IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' reduced memory with zps or sco option is impossible' ) 
    101105      ENDIF 
    102106 
    103107      ! Build the vertical coordinate system 
    104108      ! ------------------------------------ 
    105  
    106109                     CALL zgr_z        ! Reference z-coordinate system (always called) 
    107  
    108110                     CALL zgr_bat      ! Bathymetry fields (levels and meters) 
    109  
    110111      IF( ln_zco )   CALL zgr_zco      ! z-coordinate 
    111  
    112112      IF( ln_zps )   CALL zgr_zps      ! Partial step z-coordinate 
    113  
    114113      IF( ln_sco )   CALL zgr_sco      ! s-coordinate or hybrid z-s coordinate 
    115  
    116 !!bug gm control print: 
    117       IF( nprint == 1 .AND. lwp )   THEN 
    118          WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    119          WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   & 
    120             &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) ) 
    121          WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  & 
    122             &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  & 
    123             &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   & 
    124             &                   ' w ',   MINVAL( fse3w(:,:,:) ) 
    125  
    126          WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   & 
    127             &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) ) 
    128          WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  & 
    129             &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  & 
    130             &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   & 
    131             &                   ' w ',   MAXVAL( fse3w(:,:,:) ) 
    132       ENDIF 
    133 !!!bug gm 
    134  
     114      ! 
    135115   END SUBROUTINE dom_zgr 
    136116 
     
    153133      !! 
    154134      !! ** Action  : - gdept_0, gdepw_0 : depth of T- and W-point (m) 
    155       !!              -  e3t_0, e3w_0    : scale factors at T- and W-levels (m) 
    156       !! 
    157       !! Reference : 
    158       !!      Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 
    159       !! 
    160       !! History : 
    161       !!   9.0  !  03-08  (G. Madec)  F90: Free form and module 
    162       !!---------------------------------------------------------------------- 
    163       !! * Local declarations 
     135      !!              - e3t_0  , e3w_0   : scale factors at T- and W-levels (m) 
     136      !! 
     137      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 
     138      !!---------------------------------------------------------------------- 
    164139      INTEGER  ::   jk                     ! dummy loop indices 
    165140      REAL(wp) ::   zt, zw                 ! temporary scalars 
    166       REAL(wp) ::   & 
    167       zsur , za0, za1, zkth, zacr,      &  ! Values set from parameters in 
    168       zdzmin, zhmax                        ! par_CONFIG_Rxx.h90 
     141      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in 
     142      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90 
    169143      !!---------------------------------------------------------------------- 
    170144 
     
    177151      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 
    178152      ! 
    179        IF(  ppa1  == pp_to_be_computed  .AND.  & 
     153      IF(   ppa1  == pp_to_be_computed  .AND.  & 
    180154         &  ppa0  == pp_to_be_computed  .AND.  & 
    181155         &  ppsur == pp_to_be_computed           ) THEN 
    182          za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) )          & 
    183              / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) & 
    184              &                         *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      & 
    185              &                             - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  ) 
    186  
    187          za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 
    188          zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  ) 
    189  
    190        ELSE 
     156         ! 
     157         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      & 
     158            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      & 
     159            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  ) 
     160         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr ) 
     161         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  ) 
     162      ELSE 
    191163         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur 
    192        ENDIF 
    193  
    194  
    195       ! Parameter print 
    196       ! --------------- 
    197       IF(lwp) THEN 
     164      ENDIF 
     165 
     166      IF(lwp) THEN                         ! Parameter print 
    198167         WRITE(numout,*) 
    199168         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates' 
     
    222191      ! ====================== 
    223192      IF( ppkth == 0.e0 ) THEN            !  uniform vertical grid        
    224  
    225193         za1 = zhmax / FLOAT(jpk-1)  
    226194         DO jk = 1, jpk 
     
    232200            e3t_0  (jk) =  za1 
    233201         END DO 
    234  
    235       ELSE 
    236  
     202      ELSE                                ! Madec & Imbard 1996 function 
    237203         DO jk = 1, jpk 
    238204            zw = FLOAT( jk ) 
    239205            zt = FLOAT( jk ) + 0.5 
    240             gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) )  ) 
    241             gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) )  ) 
    242             e3w_0  (jk) =          za0      + za1        * TANH(      (zw-zkth)/zacr   ) 
    243             e3t_0  (jk) =          za0      + za1        * TANH(      (zt-zkth)/zacr   ) 
    244          END DO 
    245          gdepw_0(1) = 0.e0      ! force first w-level to be exactly at zero 
    246  
    247       ENDIF 
    248  
    249       ! Control and  print 
    250       ! ================== 
    251  
    252       IF(lwp) THEN 
     206            gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  ) 
     207            gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  ) 
     208            e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   ) 
     209            e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
     210         END DO 
     211         gdepw_0(1) = 0.e0                     ! force first w-level to be exactly at zero 
     212      ENDIF 
     213 
     214      IF(lwp) THEN                        ! control print 
    253215         WRITE(numout,*) 
    254216         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:' 
     
    256218         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk ) 
    257219      ENDIF 
    258  
    259       DO jk = 1, jpk 
    260          IF( e3w_0(jk)  <= 0. .OR. e3t_0(jk)  <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' ) 
    261          IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' ) 
    262       END DO 
    263  
     220      DO jk = 1, jpk                      ! control positivity 
     221         IF( e3w_0  (jk) <= 0.e0 .OR. e3t_0  (jk) <= 0.e0 )   CALL ctl_stop( ' e3w or e3t =< 0 '    ) 
     222         IF( gdepw_0(jk) <  0.e0 .OR. gdept_0(jk) <  0.e0 )   CALL ctl_stop( ' gdepw or gdept < 0 ' ) 
     223      END DO 
     224      ! 
    264225   END SUBROUTINE zgr_z 
    265226 
     
    298259      !! ** Action  : - mbathy: level bathymetry (in level index) 
    299260      !!              - bathy : meter bathymetry (in meters) 
    300       !! 
    301       !! History : 
    302       !!   9.0  !  03-08  (G. Madec)  Original code 
    303       !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates 
    304       !!---------------------------------------------------------------------- 
    305       !! * Modules used 
     261      !!---------------------------------------------------------------------- 
    306262      USE iom 
    307  
    308       !! * Local declarations 
    309       INTEGER ::   ji, jj, jl, jk       ! dummy loop indices 
    310       INTEGER ::   inum                 ! temporary logical unit 
    311       INTEGER  ::   & 
    312          ii_bump, ij_bump, ih           ! bump center position 
    313       INTEGER , DIMENSION(jpidta,jpjdta) ::   & 
    314          idta                           ! global domain integer data 
    315       REAL(wp) ::   & 
    316          r_bump, h_bump, h_oce,      &  ! bump characteristics  
    317          zi, zj, zh                     ! temporary scalars 
    318       REAL(wp), DIMENSION(jpidta,jpjdta) ::   & 
    319          zdta                           ! global domain scalar data 
     263      !! 
     264      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
     265      INTEGER  ::   inum                      ! temporary logical unit 
     266      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position 
     267      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics  
     268      REAL(wp) ::   zi     , zj     , zh      ! temporary scalars 
     269      INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data 
     270      REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data 
    320271      !!---------------------------------------------------------------------- 
    321272 
     
    324275      IF(lwp) WRITE(numout,*) '    ~~~~~~~' 
    325276 
    326       ! ======================================== 
    327       ! global domain level and meter bathymetry (idta,zdta) 
    328       ! ======================================== 
    329       !                                               ! =============== !  
    330       IF( ntopo == 0 .OR. ntopo == -1 ) THEN          ! defined by hand ! 
    331          !                                            ! =============== ! 
    332  
     277      !                                               ! ================== !  
     278      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  ! 
     279         !                                            ! ================== ! 
     280         !                                            ! global domain level and meter bathymetry (idta,zdta) 
     281         ! 
    333282         IF( ntopo == 0 ) THEN                        ! flat basin 
    334  
    335283            IF(lwp) WRITE(numout,*) 
    336284            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin' 
    337  
    338             idta(:,:) = jpkm1                            ! flat basin  
    339             zdta(:,:) = gdepw_0(jpk) 
     285            idta(:,:) = jpkm1                            ! before last level 
     286            zdta(:,:) = gdepw_0(jpk)                     ! last w-point depth 
    340287            h_oce     = gdepw_0(jpk) 
    341  
    342          ELSE                                         ! bump 
     288         ELSE                                         ! bump centered in the basin 
    343289            IF(lwp) WRITE(numout,*) 
    344290            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump' 
    345  
    346             ii_bump = jpidta / 2           ! i-index of the bump center 
    347             ij_bump = jpjdta / 2           ! j-index of the bump center 
    348             r_bump  = 50000.e0             ! bump radius (meters)        
    349             h_bump  = 2700.e0              ! bump height (meters) 
    350             h_oce   = gdepw_0(jpk)         ! background ocean depth (meters) 
     291            ii_bump = jpidta / 2                           ! i-index of the bump center 
     292            ij_bump = jpjdta / 2                           ! j-index of the bump center 
     293            r_bump  = 50000.e0                             ! bump radius (meters)        
     294            h_bump  = 2700.e0                              ! bump height (meters) 
     295            h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters) 
    351296            IF(lwp) WRITE(numout,*) '            bump characteristics: ' 
    352297            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump 
     
    354299            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index' 
    355300            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters' 
    356             ! zdta : 
    357             DO jj = 1, jpjdta 
     301            !                                         
     302            DO jj = 1, jpjdta                              ! zdta : 
    358303               DO ji = 1, jpidta 
    359304                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 
     
    362307               END DO 
    363308            END DO 
    364  
    365             ! idta : 
    366             IF( ln_sco ) THEN      ! s-coordinate (zsc       ): idta()=jpk 
     309            !                                              ! idta : 
     310            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk 
    367311               idta(:,:) = jpkm1 
    368             ELSE                   ! z-coordinate (zco or zps): step-like topography 
     312            ELSE                                                ! z-coordinate (zco or zps): step-like topography 
    369313               idta(:,:) = jpkm1 
    370314               DO jk = 1, jpkm1 
     
    377321            ENDIF 
    378322         ENDIF 
    379  
    380          ! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio) 
     323         !                                            ! set GLOBAL boundary conditions  
     324         !                                            ! Caution : idta on the global domain: use of jperio, not nperio 
    381325         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 
    382326            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1.e0 
     
    396340         ENDIF 
    397341 
    398          ! ======================================= 
    399          ! local domain level and meter bathymetry (mbathy,bathy) 
    400          ! ======================================= 
    401           
    402          mbathy(:,:) = 0                                 ! set to zero extra halo points 
    403          bathy (:,:) = 0.e0                              ! (require for mpp case) 
    404           
    405          DO jj = 1, nlcj                                 ! interior values 
     342         !                                            ! local domain level and meter bathymetries (mbathy,bathy) 
     343         mbathy(:,:) = 0                                   ! set to zero extra halo points 
     344         bathy (:,:) = 0.e0                                ! (require for mpp case) 
     345         DO jj = 1, nlcj                                   ! interior values 
    406346            DO ji = 1, nlci 
    407347               mbathy(ji,jj) = idta( mig(ji), mjg(jj) ) 
     
    409349            END DO 
    410350         END DO 
    411  
    412          !                                            ! =============== ! 
    413       ELSEIF( ntopo == 1 ) THEN                       !   read in file  ! 
    414          !                                            ! =============== ! 
    415  
    416          IF( ln_zco )   THEN 
    417             CALL iom_open ( 'bathy_level.nc', inum )   ! Level bathymetry 
     351         ! 
     352         !                                            ! ================ ! 
     353      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain) 
     354         !                                            ! ================ ! 
     355         ! 
     356         IF( ln_zco )   THEN                          ! zco : read level bathymetry  
     357            CALL iom_open( 'bathy_level.nc', inum )   
    418358            CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 
    419359            CALL iom_close (inum) 
    420360            mbathy(:,:) = INT( bathy(:,:) ) 
    421361         ENDIF 
    422  
    423          IF( ln_zps .OR. ln_sco )   THEN 
    424             CALL iom_open ( 'bathy_meter.nc', inum )   ! meter bathymetry 
     362         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry 
     363            CALL iom_open( 'bathy_meter.nc', inum )  
    425364            CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 
    426365            CALL iom_close (inum) 
     
    429368      ELSE                                            !      error      ! 
    430369         !                                            ! =============== ! 
    431          WRITE(ctmp1,*) '          parameter , ntopo = ', ntopo 
     370         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo 
    432371         CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) ) 
    433372      ENDIF 
    434  
    435       ! ======================= 
    436       ! NO closed seas or lakes 
    437       ! ======================= 
    438  
    439       IF( nclosea == 0 ) THEN 
    440          DO jl = 1, jpncs 
     373      ! 
     374      !                                               ! =========================== ! 
     375      IF( nclosea == 0 ) THEN                         !   NO closed seas or lakes   ! 
     376         DO jl = 1, jpncs                             ! =========================== ! 
    441377            DO jj = ncsj1(jl), ncsj2(jl) 
    442378               DO ji = ncsi1(jl), ncsi2(jl) 
    443                   mbathy(ji,jj) = 0                   ! suppress closed seas 
    444                   bathy (ji,jj) = 0.e0                ! and lakes 
    445 !!bug gm          bathy (ji,jj) = 300.e0                ! and lakes 
     379                  mbathy(ji,jj) = 0                   ! suppress closed seas and lakes from bathymetry 
     380                  bathy (ji,jj) = 0.e0                 
    446381               END DO 
    447382            END DO 
    448383         END DO 
    449384      ENDIF 
    450  
    451385#if defined key_orca_lev10 
    452       ! 10 time the vertical resolution 
    453       mbathy(:,:) = 10 * mbathy(:,:) 
    454       IF(lwp) WRITE(numout,*) ' ATTENTION: 300 niveaux avec bathy levels "vraie?"' 
     386      !                                               ! ================= ! 
     387      mbathy(:,:) = 10 * mbathy(:,:)                  !   key_orca_lev10  ! 
     388      !                                               ! ================= ! 
     389      IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' ) 
    455390#endif 
    456       ! =========== 
    457       ! Zoom domain  
    458       ! =========== 
    459  
    460       IF( lzoom )   CALL zgr_bat_zoom 
    461  
    462       ! ================ 
    463       ! Bathymetry check 
    464       ! ================ 
    465  
    466       CALL zgr_bat_ctl 
    467  
     391      !                                               ! =============== ! 
     392      IF( lzoom        )   CALL zgr_bat_zoom          !   Zoom domain   ! 
     393      !                                               ! =============== ! 
     394 
     395      !                                               ! =================== ! 
     396      IF( .NOT. lk_c1d )   CALL zgr_bat_ctl           !   Bathymetry check  ! 
     397      !                                               ! =================== ! 
    468398   END SUBROUTINE zgr_bat 
    469399 
     
    479409      !! 
    480410      !! ** Action  : - update mbathy: level bathymetry (in level index) 
    481       !! 
    482       !! History : 
    483       !!   9.0  !  03-08  (G. Madec)  Original code 
    484       !!---------------------------------------------------------------------- 
    485       !! * Local variables 
     411      !!---------------------------------------------------------------------- 
    486412      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers 
    487413      !!---------------------------------------------------------------------- 
    488  
     414      ! 
    489415      IF(lwp) WRITE(numout,*) 
    490416      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain' 
    491417      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~' 
    492  
     418      ! 
    493419      ! Zoom domain 
    494420      ! =========== 
    495  
     421      ! 
    496422      ! Forced closed boundary if required 
    497       IF( lzoom_w )   mbathy( mi0(jpizoom):mi1(jpizoom) , :  ) = 0 
    498       IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) ) = 0 
    499       IF( lzoom_e )   mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0 
    500       IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0 
    501  
     423      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0 
     424      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , : ) = 0 
     425      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0 
     426      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0 
     427      ! 
    502428      ! Configuration specific domain modifications 
    503429      ! (here, ORCA arctic configuration: suppress Med Sea) 
     
    508434            !                                     ! ======================= 
    509435            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea' 
    510  
    511436            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices 
    512437            ij0 =  98   ;   ij1 = 110 
     
    515440            !                                     ! ======================= 
    516441            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea' 
    517    
    518442            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe 
    519443            ij0 = 314   ;   ij1 = 370  
    520  
    521444         END SELECT 
    522445         ! 
     
    524447         ! 
    525448      ENDIF 
    526  
     449      ! 
    527450   END SUBROUTINE zgr_bat_zoom 
    528451 
     
    549472      !! ** Action  : - update mbathy: level bathymetry (in level index) 
    550473      !!              - update bathy : meter bathymetry (in meters) 
    551       !! 
    552       !! History : 
    553       !!   9.0  !  03-08  (G. Madec)  Original code 
    554       !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates 
    555       !!---------------------------------------------------------------------- 
    556       !! * Local declarations 
    557       INTEGER ::   ji, jj, jl           ! dummy loop indices 
    558       INTEGER ::   & 
    559          icompt, ibtest, ikmax          ! temporary integers 
    560       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    561          zbathy                         ! temporary workspace 
     474      !!---------------------------------------------------------------------- 
     475      INTEGER ::   ji, jj, jl                    ! dummy loop indices 
     476      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers 
     477      REAL(wp), DIMENSION(jpi,jpj) ::   zbathy   ! temporary workspace 
    562478      !!---------------------------------------------------------------------- 
    563479 
     
    566482      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~' 
    567483 
    568       ! ================ 
    569       ! Bathymetry check 
    570       ! ================ 
    571  
    572       IF( .NOT. lk_c1d )   THEN 
    573  
    574          ! Suppress isolated ocean grid points 
    575  
     484      !                                          ! Suppress isolated ocean grid points 
     485      IF(lwp) WRITE(numout,*) 
     486      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points' 
     487      IF(lwp) WRITE(numout,*)'                   -----------------------------------' 
     488      icompt = 0 
     489      DO jl = 1, 2 
     490         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
     491            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west 
     492            mbathy(jpi,:) = mbathy(  2  ,:) 
     493         ENDIF 
     494         DO jj = 2, jpjm1 
     495            DO ji = 2, jpim1 
     496               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   & 
     497                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  ) 
     498               IF( ibtest < mbathy(ji,jj) ) THEN 
     499                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   & 
     500                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 
     501                  mbathy(ji,jj) = ibtest 
     502                  icompt = icompt + 1 
     503               ENDIF 
     504            END DO 
     505         END DO 
     506      END DO 
     507      IF( icompt == 0 ) THEN 
     508         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points' 
     509      ELSE 
     510         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed' 
     511      ENDIF 
     512      IF( lk_mpp ) THEN 
     513         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     514         CALL lbc_lnk( zbathy, 'T', 1. ) 
     515         mbathy(:,:) = INT( zbathy(:,:) ) 
     516      ENDIF 
     517 
     518      !                                          ! East-west cyclic boundary conditions 
     519      IF( nperio == 0 ) THEN 
     520         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 
     521         IF( lk_mpp ) THEN 
     522            IF( nbondi == -1 .OR. nbondi == 2 ) THEN 
     523               IF( jperio /= 1 )   mbathy(1,:) = 0 
     524            ENDIF 
     525            IF( nbondi == 1 .OR. nbondi == 2 ) THEN 
     526               IF( jperio /= 1 )   mbathy(nlci,:) = 0 
     527            ENDIF 
     528         ELSE 
     529            IF( ln_zco .OR. ln_zps ) THEN 
     530               mbathy( 1 ,:) = 0 
     531               mbathy(jpi,:) = 0 
     532            ELSE 
     533               mbathy( 1 ,:) = jpkm1 
     534               mbathy(jpi,:) = jpkm1 
     535            ENDIF 
     536         ENDIF 
     537      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN 
     538         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio 
     539         mbathy( 1 ,:) = mbathy(jpim1,:) 
     540         mbathy(jpi,:) = mbathy(  2  ,:) 
     541      ELSEIF( nperio == 2 ) THEN 
     542         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio 
     543      ELSE 
     544         IF(lwp) WRITE(numout,*) '    e r r o r' 
     545         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio 
     546         !         STOP 'dom_mba' 
     547      ENDIF 
     548 
     549      ! Set to zero mbathy over islands if necessary  (lk_isl=F) 
     550      IF( .NOT. lk_isl ) THEN    ! No island 
    576551         IF(lwp) WRITE(numout,*) 
    577          IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points' 
    578          IF(lwp) WRITE(numout,*)'                   -----------------------------------' 
    579  
    580          icompt = 0 
    581  
    582          DO jl = 1, 2 
    583  
    584             IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    585                mbathy( 1 ,:) = mbathy(jpim1,:) 
    586                mbathy(jpi,:) = mbathy(  2  ,:) 
    587             ENDIF 
    588             DO jj = 2, jpjm1 
    589                DO ji = 2, jpim1 
    590                   ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),   & 
    591                      &          mbathy(ji,jj-1),mbathy(ji,jj+1) ) 
    592                   IF( ibtest < mbathy(ji,jj) ) THEN 
    593                      IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   & 
    594                         &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 
    595                      mbathy(ji,jj) = ibtest 
    596                      icompt = icompt + 1 
    597                   ENDIF 
    598                END DO 
    599             END DO 
    600  
    601          END DO 
    602          IF( icompt == 0 ) THEN 
    603             IF(lwp) WRITE(numout,*)'     no isolated ocean grid points' 
    604          ELSE 
    605             IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed' 
    606          ENDIF 
    607          IF( lk_mpp ) THEN 
     552         IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands' 
     553         IF(lwp) WRITE(numout,*) '         ----------------------------' 
     554         ! 
     555         mbathy(:,:) = MAX( 0, mbathy(:,:) ) 
     556         ! 
     557         !  Boundary condition on mbathy 
     558         IF( .NOT.lk_mpp ) THEN  
     559!!gm        !!bug ???  think about it ! 
     560            !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab 
    608561            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    609562            CALL lbc_lnk( zbathy, 'T', 1. ) 
    610563            mbathy(:,:) = INT( zbathy(:,:) ) 
    611564         ENDIF 
    612  
    613          ! 3.2 East-west cyclic boundary conditions 
    614  
    615          IF( nperio == 0 ) THEN 
    616             IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west',   & 
    617                ' boundary: nperio = ', nperio 
    618             IF( lk_mpp ) THEN 
    619                IF( nbondi == -1 .OR. nbondi == 2 ) THEN 
    620                   IF( jperio /= 1 )   mbathy(1,:) = 0 
    621                ENDIF 
    622                IF( nbondi == 1 .OR. nbondi == 2 ) THEN 
    623                   IF( jperio /= 1 )   mbathy(nlci,:) = 0 
    624                ENDIF 
    625             ELSE 
    626                IF( ln_zco .OR. ln_zps ) THEN 
    627                   mbathy( 1 ,:) = 0 
    628                   mbathy(jpi,:) = 0 
    629                ELSE 
    630                   mbathy( 1 ,:) = jpkm1 
    631                   mbathy(jpi,:) = jpkm1 
    632                ENDIF 
    633             ENDIF 
    634          ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN 
    635             IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions',   & 
    636                ' on mbathy: nperio = ', nperio 
    637             mbathy( 1 ,:) = mbathy(jpim1,:) 
    638             mbathy(jpi,:) = mbathy(  2  ,:) 
    639          ELSEIF( nperio == 2 ) THEN 
    640             IF(lwp) WRITE(numout,*) '   equatorial boundary conditions',   & 
    641                ' on mbathy: nperio = ', nperio 
    642          ELSE 
    643             IF(lwp) WRITE(numout,*) '    e r r o r' 
    644             IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio 
    645             !         STOP 'dom_mba' 
    646          ENDIF 
    647  
    648          ! Set to zero mbathy over islands if necessary  (lk_isl=F) 
    649          IF( .NOT. lk_isl ) THEN    ! No island 
    650             IF(lwp) WRITE(numout,*) 
    651             IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands' 
    652             IF(lwp) WRITE(numout,*) '         ----------------------------' 
    653  
    654             mbathy(:,:) = MAX( 0, mbathy(:,:) ) 
    655  
    656             !  Boundary condition on mbathy 
    657             IF( .NOT.lk_mpp ) THEN  
    658                !!bug ???  y reflechir! 
    659                !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab 
    660                zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    661                CALL lbc_lnk( zbathy, 'T', 1. ) 
    662                mbathy(:,:) = INT( zbathy(:,:) ) 
    663             ENDIF 
    664  
    665          ENDIF 
    666  
    667565      ENDIF 
    668566 
    669567      ! Number of ocean level inferior or equal to jpkm1 
    670  
    671568      ikmax = 0 
    672569      DO jj = 1, jpj 
     
    675572         END DO 
    676573      END DO 
    677       !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ??? 
    678  
     574!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ??? 
    679575      IF( ikmax > jpkm1 ) THEN 
    680576         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1' 
     
    685581      ENDIF 
    686582 
     583      ! control print 
    687584      IF( lwp .AND. nprint == 1 ) THEN 
    688585         WRITE(numout,*) 
    689          WRITE(numout,*) ' bathymetric field ' 
     586         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels ' 
    690587         WRITE(numout,*) ' ------------------' 
    691          WRITE(numout,*) ' number of non-zero T-levels ' 
    692          CALL prihin( mbathy, jpi, jpj, 1, jpi,   & 
    693                       1     , 1  , jpj, 1, 3  ,   & 
    694                       numout ) 
    695          WRITE(numout,*) 
    696       ENDIF 
    697  
     588         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout ) 
     589         WRITE(numout,*) 
     590      ENDIF 
     591      ! 
    698592   END SUBROUTINE zgr_bat_ctl 
    699593 
     
    706600      !! 
    707601      !! ** Method  :   set 3D coord. arrays to reference 1D array if lk_zco=T 
    708       !! 
    709       !! History : 
    710       !!        !  06-04  (R. Benshila, G. Madec)  Original code  
    711602      !!---------------------------------------------------------------------- 
    712603      INTEGER  ::   jk 
    713604      !!---------------------------------------------------------------------- 
    714  
     605      ! 
    715606      IF( .NOT.lk_zco ) THEN 
    716607         DO jk = 1, jpk 
     
    727618         END DO 
    728619      ENDIF 
    729  
     620      ! 
    730621   END SUBROUTINE zgr_zco 
    731622 
     
    787678      !!         - - - - - - -   gdept, gdepw and e3. are positives 
    788679      !!       
    789       !!  Reference : 
    790       !!     Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 
    791       !! 
    792       !!  History : 
    793       !!    8.5  !  02-09 (A. Bozec, G. Madec)  F90: Free form and module 
    794       !!    9.0  !  02-09 (A. de Miranda)  rigid-lid + islands 
    795       !!---------------------------------------------------------------------- 
    796       !! * Local declarations 
    797       INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    798       INTEGER  ::   & 
    799          ik, it            ! temporary integers 
    800        
    801       REAL(wp) ::   &   
    802          ze3tp, ze3wp,    &  ! Last ocean level thickness at T- and W-points 
    803          zdepwp,          &  ! Ajusted ocean depth to avoid too small e3t 
    804          zdepth,          &  !    "         " 
    805          zmax, zmin,      &  ! Maximum and minimum depth 
    806          zdiff               ! temporary scalar 
    807  
    808       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    809          zprt  !    "           " 
    810  
    811       LOGICAL ::  ll_print                          ! Allow  control print for debugging 
    812  
     680      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 
     681      !!---------------------------------------------------------------------- 
     682      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     683      INTEGER  ::   ik, it           ! temporary integers 
     684      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     685      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     686      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     687      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
     688      REAL(wp) ::   zdiff            ! temporary scalar 
     689      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace 
    813690      !!--------------------------------------------------------------------- 
    814       !! OPA8.5, LODYC-IPSL (2002) 
    815       !!--------------------------------------------------------------------- 
    816        
    817       ! Local variable for debugging 
    818       ll_print=.FALSE. 
    819 !!!   ll_print=.TRUE. 
    820        
    821       ! Initialization of constant 
    822       zmax = gdepw_0(jpk) + e3t_0(jpk) 
    823       zmin = gdepw_0(4) 
    824        
    825       ! Ocean depth 
    826       IF(lwp .AND. ll_print) THEN  
    827          WRITE(numout,*) 
    828          WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
    829          CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
    830       ENDIF 
    831691 
    832692      IF(lwp) WRITE(numout,*) 
     
    835695      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    836696 
     697      ll_print=.FALSE.                     ! Local variable for debugging 
     698!!    ll_print=.TRUE. 
     699       
     700      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
     701         WRITE(numout,*) 
     702         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     703         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     704      ENDIF 
     705 
    837706 
    838707      ! bathymetry in level (from bathy_meter) 
    839708      ! =================== 
    840  
    841       ! initialize mbathy to the maximum ocean level available 
    842       mbathy(:,:) = jpkm1 
    843  
    844       ! storage of land and island's number (zera and negative values) in mbathy 
     709      zmax = gdepw_0(jpk) + e3t_0(jpk)     ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 
     710      zmin = gdepw_0(4)                    ! minimum depth = 3 levels 
     711 
     712      mbathy(:,:) = jpkm1                  ! initialize mbathy to the maximum ocean level available 
     713 
     714      !                                    ! storage of land and island's number (zera and negative values) in mbathy 
     715      WHERE( bathy(:,:) <= 0. )   mbathy(:,:) = INT( bathy(:,:) ) 
     716 
     717      ! bounded value of bathy 
     718!!gm  bathy(:,:) = MIN(  zmax,  MAX( bathy(:,:), zmin )  )     !!gm : simpler   as zmin is > 0 
    845719      DO jj = 1, jpj 
    846720         DO ji= 1, jpi 
    847             IF( bathy(ji,jj) <= 0. )   mbathy(ji,jj) = INT( bathy(ji,jj) ) 
    848          END DO 
    849       END DO 
    850  
    851       ! bounded value of bathy 
    852       ! minimum depth == 3 levels 
    853       ! maximum depth == gdepw_0(jpk)+e3t_0(jpk)  
    854       ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk) 
    855       DO jj = 1, jpj 
    856          DO ji= 1, jpi 
    857             IF( bathy(ji,jj) <= 0. ) THEN 
    858                bathy(ji,jj) = 0.e0 
    859             ELSE 
    860                bathy(ji,jj) = MAX( bathy(ji,jj), zmin ) 
    861                bathy(ji,jj) = MIN( bathy(ji,jj), zmax ) 
     721            IF( bathy(ji,jj) <= 0. ) THEN   ;   bathy(ji,jj) = 0.e0 
     722            ELSE                            ;   bathy(ji,jj) = MIN(  zmax,  MAX( bathy(ji,jj), zmin )  ) 
    862723            ENDIF 
    863724         END DO 
     
    877738      END DO 
    878739 
    879      ! Scale factors and depth at T- and W-points 
    880       
    881      ! intitialization to the reference z-coordinate 
    882      DO jk = 1, jpk 
    883         gdept(:,:,jk) = gdept_0(jk) 
    884         gdepw(:,:,jk) = gdepw_0(jk) 
    885         e3t  (:,:,jk) = e3t_0  (jk) 
    886         e3w  (:,:,jk) = e3w_0  (jk) 
    887      END DO 
    888      hdept(:,:) = gdept(:,:,2 ) 
    889      hdepw(:,:) = gdepw(:,:,3 )      
    890       
    891      !  
    892      DO jj = 1, jpj 
    893         DO ji = 1, jpi 
    894            ik = mbathy(ji,jj) 
    895            ! ocean point only 
    896            IF( ik > 0 ) THEN 
    897               ! max ocean level case 
    898               IF( ik == jpkm1 ) THEN 
    899                  zdepwp = bathy(ji,jj) 
    900                  ze3tp  = bathy(ji,jj) - gdepw_0(ik) 
    901                  ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 
    902                  e3t(ji,jj,ik  ) = ze3tp 
    903                  e3t(ji,jj,ik+1) = ze3tp 
    904                  e3w(ji,jj,ik  ) = ze3wp 
    905                  e3w(ji,jj,ik+1) = ze3tp 
    906                  gdepw(ji,jj,ik+1) = zdepwp 
    907                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp 
    908                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 
    909                  ! standard case 
    910               ELSE 
    911 !!alex 
    912                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN 
    913                     gdepw(ji,jj,ik+1) = bathy(ji,jj) 
    914                  ELSE 
    915 !!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1) 
    916                     gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 
    917                  ENDIF 
    918 !!Alex 
    919 !!Alex           gdepw(ji,jj,ik+1) = bathy(ji,jj) 
    920 !!bug gm  verifier les gdepw_0 
    921                  gdept(ji,jj,ik  ) =  gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) )   & 
    922                     &                             * ((gdept_0(      ik  ) - gdepw_0(ik) )   & 
    923                     &                             / ( gdepw_0(      ik+1) - gdepw_0(ik) )) 
    924                  e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &  
    925                     &                      / ( gdepw_0(      ik+1) - gdepw_0(ik) )  
    926                  e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   & 
    927                     &                * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 
    928                  !       ... on ik+1 
    929                  e3w(ji,jj,ik+1) = e3t(ji,jj,ik) 
    930                  e3t(ji,jj,ik+1) = e3t(ji,jj,ik) 
    931                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 
    932               ENDIF 
    933            ENDIF 
    934         END DO 
    935      END DO 
    936  
    937      it = 0 
    938      DO jj = 1, jpj 
    939         DO ji = 1, jpi 
    940            ik = mbathy(ji,jj) 
    941            ! ocean point only 
    942            IF( ik > 0 ) THEN 
    943               ! bathymetry output 
    944               hdept(ji,jj) = gdept(ji,jj,ik  ) 
    945               hdepw(ji,jj) = gdepw(ji,jj,ik+1) 
    946               e3tp (ji,jj) = e3t(ji,jj,ik  ) 
    947               e3wp (ji,jj) = e3w(ji,jj,ik  ) 
    948               ! test 
    949               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  ) 
    950               IF( zdiff <= 0. .AND. lwp ) THEN  
    951                  it=it+1 
    952                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    953                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    954                  WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1),   & 
    955                                  ' zdiff   = ', zdiff 
    956                  WRITE(numout,*) ' e3tp    = ', e3t(ji,jj,ik  ), ' e3wp    = ', e3w(ji,jj,ik  ) 
    957               ENDIF 
    958            ENDIF 
    959         END DO 
     740      ! Scale factors and depth at T- and W-points 
     741      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     742         gdept(:,:,jk) = gdept_0(jk) 
     743         gdepw(:,:,jk) = gdepw_0(jk) 
     744         e3t  (:,:,jk) = e3t_0  (jk) 
     745         e3w  (:,:,jk) = e3w_0  (jk) 
     746      END DO 
     747      hdept(:,:) = gdept(:,:,2 ) 
     748      hdepw(:,:) = gdepw(:,:,3 )      
     749      !  
     750      DO jj = 1, jpj 
     751         DO ji = 1, jpi 
     752            ik = mbathy(ji,jj) 
     753            IF( ik > 0 ) THEN               ! ocean point only 
     754               ! max ocean level case 
     755               IF( ik == jpkm1 ) THEN 
     756                  zdepwp = bathy(ji,jj) 
     757                  ze3tp  = bathy(ji,jj) - gdepw_0(ik) 
     758                  ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 
     759                  e3t(ji,jj,ik  ) = ze3tp 
     760                  e3t(ji,jj,ik+1) = ze3tp 
     761                  e3w(ji,jj,ik  ) = ze3wp 
     762                  e3w(ji,jj,ik+1) = ze3tp 
     763                  gdepw(ji,jj,ik+1) = zdepwp 
     764                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp 
     765                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 
     766                  ! 
     767               ELSE                         ! standard case 
     768                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj) 
     769                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 
     770                  ENDIF 
     771!gm Bug?  check the gdepw_0 
     772                  !       ... on ik 
     773                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
     774                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   & 
     775                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) )) 
     776                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &  
     777                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) )  
     778                  e3w  (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   & 
     779                     &                  * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 
     780                  !       ... on ik+1 
     781                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik) 
     782                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik) 
     783                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 
     784               ENDIF 
     785            ENDIF 
     786         END DO 
     787      END DO 
     788      ! 
     789      it = 0 
     790      DO jj = 1, jpj 
     791         DO ji = 1, jpi 
     792            ik = mbathy(ji,jj) 
     793            IF( ik > 0 ) THEN               ! ocean point only 
     794               hdept(ji,jj) = gdept(ji,jj,ik  ) 
     795               hdepw(ji,jj) = gdepw(ji,jj,ik+1) 
     796               e3tp (ji,jj) = e3t(ji,jj,ik  ) 
     797               e3wp (ji,jj) = e3w(ji,jj,ik  ) 
     798               ! test 
     799               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  ) 
     800               IF( zdiff <= 0. .AND. lwp ) THEN  
     801                  it = it + 1 
     802                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     803                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     804                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff 
     805                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  ) 
     806               ENDIF 
     807            ENDIF 
     808         END DO 
    960809      END DO 
    961810 
    962811      ! Scale factors and depth at U-, V-, UW and VW-points 
    963  
    964       ! initialisation to z-scale factors 
    965       DO jk = 1, jpk 
     812      DO jk = 1, jpk                        ! initialisation to z-scale factors 
    966813         e3u (:,:,jk)  = e3t_0(jk) 
    967814         e3v (:,:,jk)  = e3t_0(jk) 
     
    969816         e3vw(:,:,jk)  = e3w_0(jk) 
    970817      END DO 
    971  
    972      ! Computed as the minimum of neighbooring scale factors 
    973      DO jk = 1,jpk 
    974         DO jj = 1, jpjm1 
    975            DO ji = 1, fs_jpim1   ! vector opt. 
    976               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 
    977               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 
    978               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 
    979               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 
    980            END DO 
    981         END DO 
    982      END DO 
     818      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     819         DO jj = 1, jpjm1 
     820            DO ji = 1, fs_jpim1   ! vector opt. 
     821               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 
     822               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 
     823               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 
     824               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 
     825            END DO 
     826         END DO 
     827      END DO 
     828      !                                     ! Boundary conditions 
     829      CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. ) 
     830      CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. ) 
     831      ! 
     832      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     833         WHERE( e3u (:,:,jk) == 0.e0 )   e3u (:,:,jk) = e3t_0(jk) 
     834         WHERE( e3v (:,:,jk) == 0.e0 )   e3v (:,:,jk) = e3t_0(jk) 
     835         WHERE( e3uw(:,:,jk) == 0.e0 )   e3uw(:,:,jk) = e3w_0(jk) 
     836         WHERE( e3vw(:,:,jk) == 0.e0 )   e3vw(:,:,jk) = e3w_0(jk) 
     837      END DO 
     838       
     839      ! Scale factor at F-point 
     840      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     841         e3f (:,:,jk) = e3t_0(jk) 
     842      END DO 
     843      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     844         DO jj = 1, jpjm1 
     845            DO ji = 1, fs_jpim1   ! vector opt. 
     846               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 
     847            END DO 
     848         END DO 
     849      END DO 
     850      CALL lbc_lnk( e3f, 'F', 1. )          ! Boundary conditions 
     851      ! 
     852      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     853         WHERE( e3f(:,:,jk) == 0.e0 )   e3f(:,:,jk) = e3t_0(jk) 
     854      END DO 
     855!!gm  bug ? :  must be a do loop with mj0,mj1 
     856      !  
     857      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     858      e3w(:,mj0(1),:) = e3w(:,mj0(2),:)  
     859      e3u(:,mj0(1),:) = e3u(:,mj0(2),:)  
     860      e3v(:,mj0(1),:) = e3v(:,mj0(2),:)  
     861      e3f(:,mj0(1),:) = e3f(:,mj0(2),:)  
     862 
     863      ! Control of the sign 
     864      IF( MINVAL( e3t  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' ) 
     865      IF( MINVAL( e3w  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' ) 
     866      IF( MINVAL( gdept(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' ) 
     867      IF( MINVAL( gdepw(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' ) 
    983868      
    984      ! Boundary conditions 
    985      CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. ) 
    986      CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. ) 
    987       
    988      ! set to z-scale factor if zero (i.e. along closed boundaries) 
    989      DO jk = 1, jpk 
    990         DO jj = 1, jpj 
    991            DO ji = 1, jpi 
    992               IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk) 
    993               IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk) 
    994               IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk) 
    995               IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk) 
    996            END DO 
    997         END DO 
    998      END DO 
    999       
    1000      ! Scale factor at F-point 
    1001       
    1002      ! initialisation to z-scale factors 
    1003      DO jk = 1, jpk 
    1004         e3f (:,:,jk) = e3t_0(jk) 
    1005      END DO 
    1006       
    1007      ! Computed as the minimum of neighbooring V-scale factors 
    1008      DO jk = 1, jpk 
    1009         DO jj = 1, jpjm1 
    1010            DO ji = 1, fs_jpim1   ! vector opt. 
    1011               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 
    1012            END DO 
    1013         END DO 
    1014      END DO 
    1015      ! Boundary conditions 
    1016      CALL lbc_lnk( e3f, 'F', 1. ) 
    1017       
    1018      ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1019      DO jk = 1, jpk 
    1020         DO jj = 1, jpj 
    1021            DO ji = 1, jpi 
    1022               IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk) 
    1023            END DO 
    1024         END DO 
    1025      END DO 
    1026 !!bug ? gm:  must be a do loop with mj0,mj1 
    1027  
    1028      ! we duplicate factor scales for jj = 1 and jj = 2 
    1029      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)  
    1030      e3w(:,mj0(1),:) = e3w(:,mj0(2),:)  
    1031      e3u(:,mj0(1),:) = e3u(:,mj0(2),:)  
    1032      e3v(:,mj0(1),:) = e3v(:,mj0(2),:)  
    1033      e3f(:,mj0(1),:) = e3f(:,mj0(2),:)  
    1034  
    1035  
    1036       
    1037      ! Compute gdep3w (vertical sum of e3w) 
    1038       
    1039      gdep3w   (:,:,1) = 0.5 * e3w (:,:,1) 
    1040      DO jk = 2, jpk 
    1041         gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w (:,:,jk)  
    1042      END DO 
     869      ! Compute gdep3w (vertical sum of e3w) 
     870      gdep3w(:,:,1) = 0.5 * e3w(:,:,1) 
     871      DO jk = 2, jpk 
     872         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)  
     873      END DO 
    1043874           
    1044      ! Control print 
    1045  9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ') 
    1046  9610 FORMAT(10x,i4,4f9.2) 
    1047       IF(lwp .AND. ll_print) THEN 
     875      !                                               ! ================= ! 
     876      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     877         !                                            ! ================= ! 
    1048878         DO jj = 1,jpj 
    1049879            DO ji = 1, jpi 
    1050                ik = MAX(mbathy(ji,jj),1) 
    1051                zprt(ji,jj) = e3t(ji,jj,ik) 
    1052             END DO 
    1053          END DO 
    1054          WRITE(numout,*) 
    1055          WRITE(numout,*) 'domzgr e3t(mbathy)' 
    1056          CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1057          DO jj = 1,jpj 
    1058             DO ji = 1, jpi 
    1059                ik = MAX(mbathy(ji,jj),1) 
    1060                zprt(ji,jj) = e3w(ji,jj,ik) 
    1061             END DO 
    1062          END DO 
    1063          WRITE(numout,*) 
    1064          WRITE(numout,*) 'domzgr e3w(mbathy)' 
    1065          CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1066          DO jj = 1,jpj 
    1067             DO ji = 1, jpi 
    1068                ik = MAX(mbathy(ji,jj),1) 
    1069                zprt(ji,jj) = e3u(ji,jj,ik) 
    1070             END DO 
    1071          END DO 
    1072  
    1073          WRITE(numout,*) 
    1074          WRITE(numout,*) 'domzgr e3u(mbathy)' 
    1075          CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1076          DO jj = 1,jpj 
    1077             DO ji = 1, jpi 
    1078                ik = MAX(mbathy(ji,jj),1) 
    1079                zprt(ji,jj) = e3v(ji,jj,ik) 
    1080             END DO 
    1081          END DO 
    1082          WRITE(numout,*) 
    1083          WRITE(numout,*) 'domzgr e3v(mbathy)' 
    1084          CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1085          DO jj = 1,jpj 
    1086             DO ji = 1, jpi 
    1087                ik = MAX(mbathy(ji,jj),1)  
    1088                zprt(ji,jj) = e3f(ji,jj,ik) 
    1089             END DO 
    1090          END DO 
    1091  
    1092          WRITE(numout,*) 
    1093          WRITE(numout,*) 'domzgr e3f(mbathy)' 
    1094          CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1095          DO jj = 1,jpj 
    1096             DO ji = 1, jpi 
    1097                ik =  MAX(mbathy(ji,jj),1) 
    1098                zprt(ji,jj) = gdep3w(ji,jj,ik) 
    1099             END DO 
    1100          END DO 
    1101          WRITE(numout,*) 
    1102          WRITE(numout,*) 'domzgr gdep3w(mbathy)' 
    1103          CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1104  
     880               ik = MAX( mbathy(ji,jj), 1 ) 
     881               zprt(ji,jj,1) = e3t   (ji,jj,ik) 
     882               zprt(ji,jj,2) = e3w   (ji,jj,ik) 
     883               zprt(ji,jj,3) = e3u   (ji,jj,ik) 
     884               zprt(ji,jj,4) = e3v   (ji,jj,ik) 
     885               zprt(ji,jj,5) = e3f   (ji,jj,ik) 
     886               zprt(ji,jj,6) = gdep3w(ji,jj,ik) 
     887            END DO 
     888         END DO 
     889         WRITE(numout,*) 
     890         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     891         WRITE(numout,*) 
     892         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     893         WRITE(numout,*) 
     894         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     895         WRITE(numout,*) 
     896         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     897         WRITE(numout,*) 
     898         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     899         WRITE(numout,*) 
     900         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1105901      ENDIF   
    1106  
    1107  
    1108       DO jk = 1,jpk 
    1109          DO jj = 1, jpj 
    1110             DO ji = 1, jpi  
    1111                IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN 
    1112                   IF(lwp) THEN 
    1113                      WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 ' 
    1114                      WRITE(numout,*) ' =========           --------------- ' 
    1115                      WRITE(numout,*) 
    1116                   ENDIF 
    1117                   STOP 'domzgr.psteps' 
    1118                ENDIF 
    1119                IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN 
    1120                   IF(lwp) THEN 
    1121                      WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 ' 
    1122                      WRITE(numout,*) ' =========        ------------------ ' 
    1123                      WRITE(numout,*) 
    1124                   ENDIF 
    1125                   STOP 'domzgr.psteps' 
    1126                ENDIF 
    1127             END DO 
    1128          END DO 
    1129       END DO   
    1130  
    1131    IF(lwp) THEN 
    1132       WRITE(numout,*) ' e3t lev 21 ' 
    1133       CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1134       WRITE(numout,*) ' e3w lev 21  ' 
    1135       CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1136       WRITE(numout,*) ' e3u lev 21  ' 
    1137       CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1138       WRITE(numout,*) ' e3v lev 21  ' 
    1139       CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1140       WRITE(numout,*) ' e3f lev 21  ' 
    1141       CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1142       WRITE(numout,*) ' e3t lev 22 ' 
    1143       CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1144       WRITE(numout,*) ' e3w lev 22  ' 
    1145       CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1146       WRITE(numout,*) ' e3u lev 22  ' 
    1147       CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1148       WRITE(numout,*) ' e3v lev 22  ' 
    1149       CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1150       WRITE(numout,*) ' e3f lev 22  ' 
    1151       CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
    1152    ENDIF 
    1153  
    1154       ! =========== 
    1155       ! Zoom domain  
    1156       ! =========== 
    1157  
    1158       IF( lzoom )   CALL zgr_bat_zoom 
    1159  
    1160       ! ================ 
    1161       ! Bathymetry check 
    1162       ! ================ 
    1163  
    1164       CALL zgr_bat_ctl 
    1165  
     902        
     903      !                                               ! =============== ! 
     904      IF( lzoom )   CALL zgr_bat_zoom                 !   Zoom domain   ! 
     905      !                                               ! =============== ! 
     906 
     907      !                                               ! =================== ! 
     908      IF( .NOT. lk_c1d )   CALL zgr_bat_ctl           !   Bathymetry check  ! 
     909      !                                               ! =================== ! 
    1166910   END SUBROUTINE zgr_zps 
    1167911 
     
    1223967      !!      schemes. 
    1224968      !! 
    1225       !! Reference : 
    1226       !!      Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
    1227       !! 
    1228       !! History : 
    1229       !!        !  95-12  (G. Madec)  Original code : s vertical coordinate 
    1230       !!        !  97-07  (G. Madec)  lbc_lnk call 
    1231       !!        !  97-04  (J.-O. Beismann)  
    1232       !!   8.5  !  02-06  (G. Madec)  F90: Free form and module 
    1233       !!   9.0  !  05-10  (A. Beckmann)  new stretching function 
    1234       !!---------------------------------------------------------------------- 
    1235       !! * Local declarations 
    1236       INTEGER  ::   ji, jj, jk, jl 
    1237       INTEGER  ::   iip1, ijp1, iim1, ijm1 
    1238       REAL(wp) ::   zrmax, ztaper 
    1239  
    1240       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    1241          zenv, ztmp, zmsk, zri, zrj, zhbat 
    1242  
     969      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     970      !!---------------------------------------------------------------------- 
     971      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument 
     972      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
     973      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars 
     974      REAL(wp), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk    ! 2D workspace 
     975      REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     - 
     976      !! 
    1243977      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max 
    1244978      !!---------------------------------------------------------------------- 
    1245979 
    1246       ! Read Namelist nam_zgr_sco : sigma-stretching parameters 
    1247       ! ------------------------- 
    1248       REWIND( numnam ) 
     980      REWIND( numnam )                        ! Read Namelist nam_zgr_sco : sigma-stretching parameters 
    1249981      READ  ( numnam, nam_zgr_sco ) 
    1250982 
    1251       IF(lwp) THEN 
     983      IF(lwp) THEN                            ! control print 
    1252984         WRITE(numout,*) 
    1253985         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate' 
     
    1262994      ENDIF 
    1263995 
    1264       ! ??? 
    1265       ! ------------------ 
    1266       hift(:,:) = sbot_min 
     996      hift(:,:) = sbot_min                     ! set the minimum depth for the s-coordinate 
    1267997      hifu(:,:) = sbot_min 
    1268998      hifv(:,:) = sbot_min 
    1269999      hiff(:,:) = sbot_min 
    1270  
    1271  
    1272       ! set maximum ocean depth  
    1273       ! ----------------------- 
    1274        DO jj = 1, jpj 
    1275           DO ji = 1, jpi 
    1276              bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) ) 
    1277           END DO 
    1278        END DO 
    1279  
    1280       ! Define the envelop bathymetry 
    1281       ! ============================= 
     1000      !                                        ! set maximum ocean depth  
     1001      bathy(:,:) = MIN( sbot_max, bathy(:,:) ) 
     1002 
     1003      !                                        ! ============================= 
     1004      !                                        ! Define the envelop bathymetry   (hbatt) 
     1005      !                                        ! ============================= 
    12821006      ! Smooth the bathymetry (if required)  
    1283  
    12841007      scosrf(:,:) = 0.e0              ! ocean surface depth (here zero: no under ice-shelf sea) 
    12851008      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth  
    1286  
    1287   
     1009      ! 
    12881010      ! use r-value to create hybrid coordinates 
    1289   
    12901011      DO jj = 1, jpj 
    12911012         DO ji = 1, jpi 
     
    12931014         END DO 
    12941015      END DO 
    1295  
    12961016      jl = 0 
    12971017      zrmax = 1.e0 
     
    13001020         !                                                  ! ================ ! 
    13011021         jl = jl + 1 
    1302  
    13031022         zrmax = 0.e0 
    13041023         zmsk(:,:) = 0.e0 
    1305   
    13061024         DO jj = 1, nlcj 
    13071025            DO ji = 1, nlci 
     
    13191037         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain 
    13201038         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 
    1321          ztmp(:,:) = zmsk(:,:) 
    1322          CALL lbc_lnk( zmsk, 'T', 1. ) 
     1039         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1. ) 
    13231040         DO jj = 1, nlcj 
    13241041            DO ji = 1, nlci 
     
    13261043            END DO 
    13271044         END DO 
    1328   
    1329          WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 
    1330  
     1045         !  
     1046         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 
     1047         ! 
    13311048         DO jj = 1, nlcj 
    13321049            DO ji = 1, nlci 
     
    13481065            END DO 
    13491066         END DO 
    1350   
     1067         !  
    13511068         DO jj = 1, nlcj 
    13521069            DO ji = 1, nlci 
     
    13541071            END DO 
    13551072         END DO 
    1356  
     1073         ! 
    13571074         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero 
    1358          ztmp(:,:) = zenv(:,:) 
    1359          CALL lbc_lnk( zenv, 'T', 1. ) 
     1075         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1. ) 
    13601076         DO jj = 1, nlcj 
    13611077            DO ji = 1, nlci 
     
    13661082      END DO                                                !     End loop     ! 
    13671083      !                                                     ! ================ ! 
    1368  
    1369       ! save envelop bathymetry in hbatt 
     1084      ! 
     1085      !                                        ! envelop bathymetry saved in hbatt 
    13701086      hbatt(:,:) = zenv(:,:)  
    1371  
    1372 !!bug gm  :   CAUTION:  tapering hard coded !!!!  orca2 only 
    1373 !!bug gm                NOT valid in mpp ===> taper have been changed 
    1374  
    1375        DO jj = 1, jpj 
    1376           DO ji = 1, jpi 
    1377 !!bug mpp    ztaper = EXP( -(FLOAT(jj-74)/10.)**2 ) 
    1378              ztaper = EXP( -(gphit(ji,jj)/8)**2 ) 
    1379              hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper) 
    1380           END DO 
    1381        END DO 
    1382  
    1383       DO jj = 1, jpj 
    1384          zrmax  = EXP( -(gphit(10,jj)/8)**2 ) 
    1385          ztaper = EXP( -(FLOAT(jj-74)/10.)**2 ) 
    1386          IF( nprint == 1 .AND. lwp ) WRITE(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax 
    1387       END DO 
    1388  
    1389       IF( nprint == 1 .AND. lwp )   THEN 
    1390          WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 
    1391          WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 
    1392       ENDIF 
    1393  
    1394       ! Control print 
     1087      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN 
     1088         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
     1089         DO jj = 1, jpj 
     1090            DO ji = 1, jpi 
     1091               ztaper = EXP( -(gphit(ji,jj)/8)**2 ) 
     1092               hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper) 
     1093            END DO 
     1094         END DO 
     1095      ENDIF 
     1096      ! 
     1097      IF(lwp) THEN                             ! Control print 
     1098         WRITE(numout,*) 
     1099         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 
     1100         WRITE(numout,*) 
     1101         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout ) 
     1102         IF( nprint == 1 )   THEN         
     1103            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 
     1104            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 
     1105         ENDIF 
     1106      ENDIF 
     1107 
     1108      !                                        ! ============================== 
     1109      !                                        !   hbatu, hbatv, hbatf fields 
     1110      !                                        ! ============================== 
    13951111      IF(lwp) THEN 
    13961112         WRITE(numout,*) 
    1397          WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 
    1398          WRITE(numout,*) 
    1399          CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout ) 
    1400       ENDIF 
    1401  
    1402       ! 1. computation of hbatu, hbatv, hbatf fields 
    1403       ! -------------------------------------------- 
    1404  
     1113         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min 
     1114      ENDIF 
    14051115      hbatu(:,:) = sbot_min 
    14061116      hbatv(:,:) = sbot_min 
    14071117      hbatf(:,:) = sbot_min 
    1408  
    1409       IF(lwp) THEN 
    1410          WRITE(numout,*) 
    1411          WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min 
    1412          WRITE(numout,*) 
    1413       ENDIF 
    1414  
    14151118      DO jj = 1, jpjm1 
    14161119        DO ji = 1, jpim1 
    1417            hbatu(ji,jj)= 0.5 *( hbatt(ji  ,jj)+hbatt(ji+1,jj  ) ) 
    1418            hbatv(ji,jj)= 0.5 *( hbatt(ji  ,jj)+hbatt(ji  ,jj+1) ) 
    1419            hbatf(ji,jj)= 0.25*( hbatt(ji  ,jj)+hbatt(ji  ,jj+1)   & 
    1420                                +hbatt(ji+1,jj)+hbatt(ji+1,jj+1) ) 
     1120           hbatu(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) ) 
     1121           hbatv(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) ) 
     1122           hbatf(ji,jj) = 0.25* ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   & 
     1123              &                 + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 
    14211124        END DO 
    14221125      END DO 
    1423   
     1126      !  
    14241127      ! Apply lateral boundary condition 
    1425       ! CAUTION: retain non zero value in the initial file 
    1426       !          this should be OK for orca cfg, not for EEL 
    1427       zhbat(:,:) = hbatu(:,:) 
    1428       CALL lbc_lnk( hbatu, 'U', 1. ) 
     1128!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 
     1129      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1. ) 
    14291130      DO jj = 1, jpj 
    14301131         DO ji = 1, jpi 
     
    14351136         END DO 
    14361137      END DO 
    1437  
    1438       zhbat(:,:) = hbatv(:,:) 
    1439       CALL lbc_lnk( hbatv, 'V', 1. ) 
     1138      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1. ) 
    14401139      DO jj = 1, jpj 
    14411140         DO ji = 1, jpi 
     
    14461145         END DO 
    14471146      END DO 
    1448  
    1449       zhbat(:,:) = hbatf(:,:) 
    1450       CALL lbc_lnk( hbatf, 'F', 1. ) 
     1147      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1. ) 
    14511148      DO jj = 1, jpj 
    14521149         DO ji = 1, jpi 
     
    14581155      END DO 
    14591156 
    1460  
    14611157!!bug:  key_helsinki a verifer 
    14621158      hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
     
    14661162 
    14671163      IF( nprint == 1 .AND. lwp )   THEN 
    1468          WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift(:,:) ), ' f ', MAXVAL( hiff(:,:) ),  & 
    1469             &                        ' u ',   MAXVAL( hifu(:,:) ), ' v ', MAXVAL( hifv(:,:) ) 
    1470          WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift(:,:) ), ' f ', MINVAL( hiff(:,:) ),  & 
    1471             &                        ' u ',   MINVAL( hifu(:,:) ), ' v ', MINVAL( hifv(:,:) ) 
     1164         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  & 
     1165            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 
     1166         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  & 
     1167            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 
    14721168         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  & 
    14731169            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 
     
    14771173!! helsinki 
    14781174 
    1479       ! 2. Computation of gsig and esig fields 
    1480       ! -------------------------------------- 
    1481  
    1482       ! 2.1 Coefficients for model level depth at w- and t-levels 
    1483  
     1175      !                                            ! ======================= 
     1176      !                                            !   s-ccordinate fields     (gdep., e3.) 
     1177      !                                            ! ======================= 
     1178      ! 
     1179      ! non-dimensional "sigma" for model level depth at w- and t-levels 
    14841180      DO jk = 1, jpk 
    14851181        gsigw(jk) = -fssig( FLOAT(jk)-0.5 ) 
    14861182        gsigt(jk) = -fssig( FLOAT(jk)     ) 
    14871183      END DO 
    1488  
    14891184      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
    1490  
    1491       ! 2.2 Coefficients for vertical scale factors at w-, t- levels 
    1492  
    1493 !!bug gm :  define it from analytical function, not like juste bellow.... 
    1494 !!          or betteroffer the 2 possibilities.... 
    1495       DO jk=2,jpk 
    1496         esigw(jk)=gsigt(jk)-gsigt(jk-1) 
    1497       END DO 
    1498       DO jk=1,jpkm1 
    1499         esigt(jk)=gsigw(jk+1)-gsigw(jk) 
    1500       END DO 
    1501         esigw(1)=esigw(2) 
    1502         esigt(jpk)=esigt(jpkm1) 
    1503  
    1504       ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors 
    1505  
     1185      ! 
     1186      ! Coefficients for vertical scale factors at w-, t- levels 
     1187!!gm bug :  define it from analytical function, not like juste bellow.... 
     1188!!gm        or betteroffer the 2 possibilities.... 
     1189      DO jk = 1, jpkm1 
     1190         esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
     1191         esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
     1192      END DO 
     1193      esigw( 1 ) = esigw(  2  ) 
     1194      esigt(jpk) = esigt(jpkm1) 
     1195 
     1196!!gm  original form 
     1197!!org DO jk = 1, jpk 
     1198!!org    esigt(jk)=fsdsig( FLOAT(jk)     ) 
     1199!!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 
     1200!!org END DO 
     1201!!gm 
     1202      ! 
     1203      ! Coefficients for vertical depth as the sum of e3w scale factors 
    15061204      gsi3w(1) = 0.5 * esigw(1) 
    15071205      DO jk = 2, jpk 
    1508         gsi3w(jk) = gsi3w(jk-1)+ esigw(jk) 
    1509       END DO 
    1510  
    1511 !!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
     1206         gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
     1207      END DO 
     1208!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    15121209      DO jk = 1, jpk 
    1513 !! change the solver stat.... 
    1514 !!       zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 
    1515 !!       gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 
    1516          gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*(FLOAT(jk)-0.5)/FLOAT(jpkm1)) 
    1517 !!       zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)  
    1518 !!       gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 
    1519 !!       gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw) 
    1520          gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1)) 
    1521          gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1)) 
    1522       END DO 
    1523  
    1524 !!bug gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
     1210         zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 
     1211         zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)  
     1212         gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 
     1213         gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 
     1214         gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw) 
     1215      END DO 
     1216!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
    15251217      DO jj = 1, jpj 
    15261218         DO ji = 1, jpi 
     
    15301222              e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 
    15311223              e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1)) 
    1532  
     1224              ! 
    15331225              e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1)) 
    15341226              e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 
     
    15371229         END DO 
    15381230      END DO 
    1539  
    1540 ! HYBRID 
     1231      ! 
     1232      ! HYBRID :  
    15411233      DO jj = 1, jpj 
    15421234         DO ji = 1, jpi 
     
    15471239         END DO 
    15481240      END DO 
    1549  
    1550       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    1551  
    1552  
    1553       ! =========== 
    1554       ! Zoom domain  
    1555       ! =========== 
    1556  
    1557       IF( lzoom )   CALL zgr_bat_zoom 
    1558  
    1559       ! 2.4 Control print 
    1560  
    1561       IF(lwp) THEN 
     1241      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
     1242         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
     1243 
     1244 
     1245      !                                            ! =========== 
     1246      IF( lzoom )   CALL zgr_bat_zoom              ! Zoom domain  
     1247      !                                            ! =========== 
     1248 
     1249      !                                            ! ============= 
     1250      IF(lwp) THEN                                 ! Control print 
     1251         !                                         ! ============= 
    15621252         WRITE(numout,*)  
    15631253         WRITE(numout,*) ' domzgr: vertical coefficients for model level' 
    1564          WRITE(numout,9400) 
    1565          WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk) 
    1566       ENDIF 
    1567  9400 FORMAT(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w') 
    1568  9410 FORMAT(10x,i4,5f11.4) 
    1569  
    1570       IF(lwp) THEN 
     1254         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" ) 
     1255         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk ) 
     1256      ENDIF 
     1257      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain 
     1258         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     1259         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   & 
     1260            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) ) 
     1261         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   & 
     1262            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   & 
     1263            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   & 
     1264            &                          ' w ', MINVAL( fse3w (:,:,:) ) 
     1265 
     1266         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   & 
     1267            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) ) 
     1268         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   & 
     1269            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   & 
     1270            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   & 
     1271            &                          ' w ', MAXVAL( fse3w (:,:,:) ) 
     1272      ENDIF 
     1273      ! 
     1274      IF(lwp) THEN                                  ! selected vertical profiles 
    15711275         WRITE(numout,*) 
    15721276         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 
    15731277         WRITE(numout,*) ' ~~~~~~  --------------------' 
    1574          WRITE(numout,9420) 
    1575          WRITE(numout,9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk),     & 
    1576                              fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk) 
     1278         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')") 
     1279         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     & 
     1280            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk ) 
    15771281         DO jj = mj0(20), mj1(20) 
    15781282            DO ji = mi0(20), mi1(20) 
     
    15801284               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    15811285               WRITE(numout,*) ' ~~~~~~  --------------------' 
    1582                WRITE(numout,9420) 
    1583                WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),     & 
    1584                     &                 fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk) 
     1286               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')") 
     1287               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     & 
     1288                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 
    15851289            END DO 
    15861290         END DO 
     
    15901294               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    15911295               WRITE(numout,*) ' ~~~~~~  --------------------' 
    1592                WRITE(numout,9420) 
    1593                WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),     & 
    1594                     &                 fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk) 
    1595             END DO 
    1596          END DO 
    1597       ENDIF 
    1598  
    1599  9420 FORMAT(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ') 
    1600  9430 FORMAT(10x,i4,4f9.2) 
    1601  
    1602 !!bug gm   no more necessary?  if ! defined key_helsinki 
     1296               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')") 
     1297               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     & 
     1298                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 
     1299            END DO 
     1300         END DO 
     1301      ENDIF 
     1302 
     1303!!gm bug?  no more necessary?  if ! defined key_helsinki 
    16031304      DO jk = 1, jpk 
    16041305         DO jj = 1, jpj 
    16051306            DO ji = 1, jpi 
    16061307               IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN 
    1607                   IF(lwp) THEN 
    1608                      WRITE(numout,*) 
    1609                      WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 ' 
    1610                      WRITE(numout,*) ' =========           --------------- ' 
    1611                      WRITE(numout,*) 
    1612                      WRITE(numout,*) '             point (i,j,k)= ',ji,jj,jk 
    1613                      WRITE(numout,*) 
    1614                   ENDIF 
    1615                   STOP 'domzgr' 
     1308                  WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1309                  CALL ctl_stop( ctmp1 ) 
    16161310               ENDIF 
    16171311               IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN 
    1618                   IF(lwp) THEN 
    1619                      WRITE(numout,*) 
    1620                      WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 ' 
    1621                      WRITE(numout,*) ' =========        ------------------ ' 
    1622                      WRITE(numout,*) 
    1623                      WRITE(numout,*) '             point (i,j,k)= ', ji, jj, jk 
    1624                      WRITE(numout,*) 
    1625                   ENDIF 
    1626                   STOP 'domzgr' 
     1312                  WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     1313                  CALL ctl_stop( ctmp1 ) 
    16271314               ENDIF 
    16281315            END DO 
    16291316         END DO 
    16301317      END DO 
    1631 !!bug gm    #endif 
    1632  
     1318!!gm bug    #endif 
     1319      ! 
    16331320   END SUBROUTINE zgr_sco 
    16341321 
    16351322#endif 
    1636  
    16371323 
    16381324   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.