New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15278 – NEMO

Changeset 15278


Ignore:
Timestamp:
2021-09-22T18:44:26+02:00 (15 months ago)
Author:
dbruciaferri
Message:

reorganising code in more general structure and adding partial-steps for MEs

Location:
NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_cfg

    r15131 r15278  
    1515!----------------------------------------------------------------------- 
    1616   ! 
    17    ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens.  
     17   ln_e3_dep   = .false.   ! =T : e3=dk[depth] in discret sens.  
    1818   !                       !      ===>>> will become the only possibility in v4.0 
    1919   !                       ! =F : e3 analytical derivative of depth function 
     
    3434&namzgr        !   vertical coordinate 
    3535!----------------------------------------------------------------------- 
    36    ln_mes      = .true.    !  Multi-Envelope s-coordinate 
     36   ln_zps      = .true.    !  z-coordinate with partial steps 
    3737   ln_linssh   = .true.    !  linear free surface 
    38 / 
    39 !----------------------------------------------------------------------- 
    40 &namzgr_mes    !   MEs-coordinate                      
    41 !----------------------------------------------------------------------- 
    42    ln_envl     =   .TRUE. , .TRUE. , .TRUE. , .TRUE., .FALSE.  ! (T/F) If the envelope is used 
    43    nn_strt     =     1    ,   1    ,    1   ,   1    ,   1     ! Stretch. funct.: Madec 1996 (0) or 
    44                                                                ! Song & Haidvogel 1994 (1) or                     
    45                                                                ! Siddorn & Furner 2012 (2) 
    46    nn_slev     =     3    ,   37   ,   17   ,   18   ,    0    ! number of s-lev between env(n-1) 
    47                                                                ! and env(n) 
    48    rn_e_hc     =     0.0  ,    0.0 ,   0.0  ,   0.0  ,   0.0   ! critical depth for transition to         
    49                                                                ! stretch. coord. 
    50    rn_e_th     =     1.0  ,    1.2 ,   2.0  ,   0.0  ,   0.0   ! surf. control param.:                   
    51                                                                ! SH94 or MD96: 0<=th<=20 
    52                                                                ! SF12: thickness surf. cell 
    53    rn_e_bb     =     0.9  ,    0.9 ,   0.85 ,   0.0  ,   0.0   ! bot. control param.: 
    54                                                                ! SH94 or MD96: 0<=bb<=1 
    55                                                                ! SF12: offset for calculating Zb 
    56    rn_e_al     =     0.0  ,    0.0 ,   0.0  ,   0.0  ,   0.0   ! alpha stretching param with SF12 
    57    rn_e_ba     =     0.0  ,    0.0 ,   0.0  ,   0.0  ,   0.0   ! SF12 bathymetry scaling factor for  
    58                                                                ! calculating Zb 
    59    rn_bot_min  = 10.0       ! minimum depth of the ocean bottom (>0) (m) 
    60    rn_bot_max  = 5900.0     ! maximum depth of the ocean bottom (= ocean depth) (>0) (m) 
    61  
    62    ln_loc_mes  = .TRUE. 
    6338/ 
    6439!----------------------------------------------------------------------- 
     
    6742   nn_bathy    =  1 
    6843   nn_msh      =  1 
    69    jphgr_msh   =  0 
    70    ldbletanh   = .TRUE. 
    71    ppglam0     =  999999.d0 
    72    ppgphi0     =  999999.d0 
    73    ppe1_deg    =  999999.d0 
    74    ppe2_deg    =  999999.d0 
    75    ppe1_m      =  999999.d0 
    76    ppe2_m      =  999999.d0                  
    77    ppa0        =  103.9530096000000 
    78    ppa1        =  2.415951269000000 
    79    ppa2        =  100.7609285000000 
    80    ppacr       =  7.0 
    81    ppacr2      =  13.0 
    82    ppdzmin     =  999999.0 
    83    pphmax      =  999999.0 
    84    ppkth       =  15.35101370000000 
    85    ppkth2      =  48.02989372000000 
    86    ppsur       = -3958.951371276829 
    8744   rn_atfp     =  0.1 
    8845   rn_e3zps_min=  25.0 
     
    9047   rn_hmin     = -8.0 
    9148   rn_rdt      =  1350.0 
     49   jphgr_msh   =       0               !  type of horizontal mesh 
     50   ppglam0     =  999999.0             !  longitude of first raw and column T-point (jphgr_msh = 1) 
     51   ppgphi0     =  999999.0             ! latitude  of first raw and column T-point (jphgr_msh = 1) 
     52   ppe1_deg    =  999999.0             !  zonal      grid-spacing (degrees) 
     53   ppe2_deg    =  999999.0             !  meridional grid-spacing (degrees) 
     54   ppe1_m      =  999999.0             !  zonal      grid-spacing (degrees) 
     55   ppe2_m      =  999999.0             !  meridional grid-spacing (degrees) 
     56   ppsur       =   -3958.951371276829  !  ORCA r4, r2 and r05 coefficients 
     57   ppa0        =    103.9530096000000  ! (default coefficients) 
     58   ppa1        =      2.415951269000000  ! 
     59   ppkth       =     15.35101370000000 ! 
     60   ppacr       =       7.0             ! 
     61   ppdzmin     =  999999.0             !  Minimum vertical spacing 
     62   pphmax      =  999999.0             !  Maximum depth 
     63   ldbletanh   =   .TRUE.              !  Use/do not use double tanf function for vertical coordinates 
     64   ppa2        =     100.7609285000000 !  Double tanh function parameters 
     65   ppkth2      =      48.02989372000000  ! 
     66   ppacr2      =      13.              ! 
    9267/ 
    9368!----------------------------------------------------------------------- 
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_ref

    r15129 r15278  
    107107   ln_isfcav   = .false.   !  ice shelf cavity 
    108108   ln_linssh   = .false.   !  linear free surface 
     109   ln_loc_zgr  = .false.   !  to localise the chosen vert. coord. system 
     110                           !  If TRUE, the user must provide: 
     111                           !     1) a bathy_meter.nc file including 
     112                           !        *) s2z_msk: 2D mask for s(=2), s2z(=1) and z(=0) areas 
     113                           !        *) s2z_wgt: 2D field with 1 in s-areas, distance-based 
     114                           !                    weights (0<wgt<1) in s2z-areas and 0 elsewhere. 
     115                           !     2) a domain_cfg_global.nc file including all the geometrical 
     116                           !        information describing the vertical grid used globally. 
    109117/ 
    110118!----------------------------------------------------------------------- 
     
    160168   rn_bot_max  = 2201.5    ! maximum depth of the ocean bottom (= ocean depth) (>0) (m) 
    161169 
    162    ln_loc_mes  = .FALSE.   ! To use MEs-coordinates only in localised regions.  
    163                            ! If TRUE, the user must provide: 
    164                            !     1) a bathy_meter.nc file including 
    165                            !        *) s2z_msk: 2D mask for s(=2), s2z(=1) and z(=0) areas 
    166                            !        *) s2z_wgt: 2D field with 1 in s-areas, distance-based  
    167                            !                    weights (0<wgt<1) in s2z-areas and 0 elsewhere. 
    168                            !     2) a domain_cfg_global.nc file including all the geometrical 
    169                            !        information describing the vertical grid used globally.   
     170   ln_pst_mes   = .FALSE.  ! To use partial-steps algorithm where MEs-coordinates are used. 
     171   ln_pst_l2g   = .FALSE.  ! To use partial-steps algorithm in the transition area  
     172                           ! (requires ln_loc_zgr = .TRUE.) 
     173   rn_e3pst_min = 20.      ! partial step thickness is set larger than the minimum of 
     174   rn_e3pst_rat = 0.1      ! rn_e3pst_min and rn_e3pst_rat*e3t, with 0<rn_e3pst_rat<1 
    170175/ 
    171176!----------------------------------------------------------------------- 
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/dom_oce.f90

    r15271 r15278  
    178178   LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF  
    179179   LOGICAL, PUBLIC ::   ln_linssh    !: variable grid flag 
     180   LOGICAL, PUBLIC ::   ln_loc_zgr   !: To localise (.TRUE.) or not (.FALSE.) the  
     181                                     !  chosen vertical coordinate system 
    180182 
    181183   !                                                        !  ref.   ! before  !   now   ! after  ! 
     
    245247   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    246248 
    247    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mes_msk      ! mask for local-ME coordinates (ln_loc_mes=TRUE)   
    248                                                                        !    MEs-area (=1), 
    249                                                                        !    z-area (=0) 
    250  
     249   ! LOCAL VERTICAL COORDINATE SYSTEM  
     250   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: l2g_msk ! mask for identifying areas using 
     251                                                                  ! -) the local coord. system (=2), 
     252                                                                  ! -) transitioning from local to global (=1), 
     253                                                                  ! -) using global coord. system (=0) 
     254                                                                  ! If ln_loc_zgr=.FALSE. then l2g_msk(:,:)=2 
    251255 
    252256   !!---------------------------------------------------------------------- 
     
    370374 
    371375      ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       & 
    372          &     tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                       &  
    373          &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 
    374          &     mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , mes_msk(jpi,jpj) , STAT=ierr(9) ) 
     376         &      tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                      & 
     377         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj), & 
     378         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) ,                   & 
     379         &      l2g_msk(jpi,jpj) , STAT=ierr(9) ) 
    375380 
    376381! (ISF) Allocation of basic array    
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domwri.f90

    r15271 r15278  
    1616   !!---------------------------------------------------------------------- 
    1717   USE dom_oce         ! ocean space and time domain 
    18    USE mes, ONLY : ln_loc_mes 
    1918   USE in_out_manager  ! I/O manager 
    2019   USE iom             ! I/O library 
     
    596595         END DO 
    597596      END DO 
    598       IF( ln_loc_mes ) zx1 = zx1 * mes_msk  
     597      IF( ln_loc_zgr .AND. ln_mes ) WHERE (l2g_msk(:,:) == 0.0)  zx1(:,:) = 0.0  
    599598      CALL lbc_lnk( zx1, 'T', 1. ) 
    600599      ! 
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domzgr.f90

    r15152 r15278  
    113113      INTEGER ::   ios 
    114114      ! 
    115       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_mes, ln_isfcav, ln_linssh 
     115      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_mes, ln_isfcav, ln_linssh, ln_loc_zgr 
    116116      !!---------------------------------------------------------------------- 
    117117      ! 
     
    141141 
    142142      IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time' 
     143      IF( ln_loc_zgr.AND. lwp) WRITE(numout,*) '   the chosen vert. coord. system is implemented LOCALLY' 
    143144 
    144145      ioptio = 0                       ! Check Vertical coordinate options 
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90

    r15271 r15278  
    4444   REAL(wp)           :: rn_e_al(max_nn_env)     ! Array specifing stretching parameter rn_alpha of SF12 for 
    4545                                                 ! each vertical sub-zone 
    46    LOGICAL, PUBLIC    :: ln_loc_mes              ! To use localised MEs (.TRUE.) or not (.FALSE. 
     46   LOGICAL, PUBLIC    :: ln_pst_mes              ! To apply partial steps to MEs (.TRUE.) or not (.FALSE.). 
     47   LOGICAL, PUBLIC    :: ln_pst_l2g              ! To apply partial steps to the transition zone (.TRUE.)  
     48                                                 ! or not (.FALSE.) when ln_loc_zgr = .TRUE.  
     49   REAL, PUBLIC       :: rn_e3pst_min            ! partial steps parameters (require ln_pst_mes = .TRUE.) 
     50   REAL, PUBLIC       :: rn_e3pst_rat            ! partial steps parameters (require ln_pst_mes = .TRUE.) 
    4751   ! 
    4852   REAL(wp), POINTER, DIMENSION(:,:,:) :: envlt  ! array for the envelopes 
     
    669673      REAL(wp)                     :: rn_bot_max ! max depth of ocean bottom (= ocean depth) (>0) (m) 
    670674 
    671       NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl, & 
    672                           nn_strt    , nn_slev    , rn_e_hc, & 
    673                           rn_e_th    , rn_e_bb    , rn_e_ba, & 
    674                           rn_e_al    , ln_loc_mes 
     675      NAMELIST/namzgr_mes/rn_bot_min  , rn_bot_max  , ln_envl   , & 
     676          &               nn_strt     , nn_slev     , rn_e_hc   , & 
     677          &               rn_e_th     , rn_e_bb     , rn_e_ba   , & 
     678          &               rn_e_al     , ln_pst_mes  , ln_pst_l2g, & 
     679          &               rn_e3pst_min, rn_e3pst_rat 
    675680      !!---------------------------------------------------------------------- 
    676681      ! 
     
    720725      WRITE(ctlmes,*) 'TOT number of levels (jpk) IS DIFFERENT from sum over nn_slev(:)' 
    721726      IF ( SUM(nn_slev(:)) /= jpk ) CALL ctl_stop( ctlmes ) 
     727      ! 
     728      IF( .NOT.(ln_loc_zgr) .AND. ln_pst_l2g) ln_pst_l2g = .FALSE. 
    722729      ! 
    723730      ! 3) Reading Bathymetry and envelopes 
     
    810817            WRITE(numout,*) '          ------------------------------------------------------------------' 
    811818         ENDDO 
     819         IF( ln_loc_zgr) THEN 
     820           WRITE(numout,*) '' 
     821           WRITE(numout,*) '          LOCALISED MEs ' 
     822           WRITE(numout,*) '          ============= ' 
     823         ENDIF 
     824         IF( ln_pst_mes .OR. ln_pst_l2g) THEN 
     825           WRITE(numout,*) '' 
     826           WRITE(numout,*) '          APPLYING PARTIAL-STEPS to ' 
     827           WRITE(numout,*) '          ========================= ' 
     828           WRITE(numout,*) '' 
     829           WRITE(numout,*) '            a) MEs area       : ', ln_pst_mes 
     830           IF( ln_loc_zgr) WRITE(numout,*) '            b) Transition area: ', ln_pst_l2g 
     831           WRITE(numout,*) '' 
     832           WRITE(numout,*) '            Partial step thickness is set larger than the minimum of' 
     833           WRITE(numout,*) '' 
     834           WRITE(numout,*) '                  rn_e3pst_min = ', rn_e3pst_min 
     835           WRITE(numout,*) '            and' 
     836           WRITE(numout,*) '                  rn_e3pst_rat = ', rn_e3pst_rat 
     837           WRITE(numout,*) '' 
     838           WRITE(numout,*) '            with 0 < rn_e3pst_rat < 1' 
     839         ENDIF 
    812840      ENDIF 
    813841 
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90

    r15271 r15278  
    1212   USE closea            ! closed seas 
    1313   USE mes               ! MEs-coordinates 
     14   USE zgrloc            ! localised vert. coor. system  
    1415   ! 
    1516   USE in_out_manager    ! I/O manager 
     
    2425 
    2526   PUBLIC   zgr_mes        ! called by domzgr.F90 
    26    ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.) 
    27    ! 
    28    REAL(wp), POINTER, DIMENSION(:)     :: e3t_1d_in  , e3w_1d_in   !: ref. scale factors (m) 
    29    REAL(wp), POINTER, DIMENSION(:)     :: gdept_1d_in, gdepw_1d_in !: ref. depth (m) 
    30    ! 
    31    REAL(wp), POINTER, DIMENSION(:,:)   :: mbathy_in  
    32    ! 
    33    REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m] 
    34    REAL(wp), POINTER, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in         !:   -      - 
    35    !  
    36    REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_in, gdepw_in               !: depths [m] 
    3727 
    3828   !! * Substitutions 
     
    5040      !!              MEs-coordinates systems 
    5141      !!---------------------------------------------------------------------- 
    52       ! 
    53       ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.) 
    54       ! 
     42      INTEGER  ::   ji, jj, jk     ! dummy loop index 
    5543      !----------------------------------------------------------------------- 
    5644      ! 
     45      ! Initialise the mask to use MEs-coord. globally 
     46      l2g_msk(:,:) = 2.0 
     47      
    5748      ! Generating a global MEs vertical grid 
     49      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    5850      CALL mes_build 
    59  
     51       
    6052      ! Local MEs 
    61       IF( ln_loc_mes ) THEN 
    62          ! 
    63          CALL wrk_alloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in) 
    64          CALL wrk_alloc( jpi, jpj, mbathy_in) 
    65          CALL wrk_alloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in) 
    66          CALL wrk_alloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in) 
    67          ! 
    68          ! Reading the input vertical grid that will be used globally 
    69          CALL zgr_dom_read 
    70          ! Creating the local MEs within the global input vertical grid 
    71          CALL zgr_mes_local 
    72          ! 
    73          CALL wrk_dealloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in) 
    74          CALL wrk_dealloc( jpi, jpj, mbathy_in) 
    75          CALL wrk_dealloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in) 
    76          CALL wrk_dealloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in) 
    77          ! 
     53      ! ~~~~~~~~~ 
     54      IF( ln_loc_zgr )  CALL zgr_loc 
     55       
     56      ! MEs with partial steps 
     57      ! ~~~~~~~~~~~~~~~~~~~~~~ 
     58      IF(ln_pst_mes .OR. ln_pst_l2g)  CALL zgr_pst_mes 
     59      ! 
     60      ! Computing gdep3w 
     61      gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 
     62      DO jk = 2, jpk 
     63         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     64      END DO 
     65      ! 
     66      ! From here equal to sco code - domzgr.F90 line 2183 
     67      ! 
     68      CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 
     69      CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 
     70      CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 
     71      CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 
     72      CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 
     73      CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 
     74      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     75      ! 
     76      WHERE (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
     77      WHERE (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
     78      WHERE (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
     79      WHERE (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
     80      WHERE (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
     81      WHERE (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
     82      WHERE (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     83      ! 
     84      IF ( ln_loc_zgr ) THEN 
     85         IF ( .NOT. ln_pst_l2g ) THEN 
     86            ! Only in the transition zone since MEs 
     87            ! zone has been already taken care of 
     88            IF ( lwp ) WRITE(numout,*) 'Refine mbathy' 
     89            DO jj = 1, jpj 
     90               DO ji = 1, jpi 
     91                  IF ( l2g_msk(ji,jj) == 1._wp ) THEN 
     92                     DO jk = 1, jpkm1 
     93                        IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
     94                     END DO 
     95                     IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
     96                  END IF 
     97               END DO 
     98            END DO 
     99         END IF 
    78100      END IF 
    79101 
     
    82104! ===================================================================================================== 
    83105 
    84    SUBROUTINE zgr_dom_read 
     106   SUBROUTINE zgr_pst_mes 
    85107     !!--------------------------------------------------------------------- 
    86      !!              ***  ROUTINE zgr_dom_read  *** 
     108     !!              ***  ROUTINE zgr_pst_mes  *** 
    87109     !! 
    88      !! ** Purpose :   Read the vertical information in the domain configuration file 
     110     !! ** Purpose : Apply partial steps to MEs coordinates  
     111     !!              (and transition zone if using local MEs) 
    89112     !! 
    90113     !!---------------------------------------------------------------------- 
    91      INTEGER  ::   jk     ! dummy loop index 
    92      INTEGER  ::   inum   ! local logical unit 
    93      REAL(WP) ::   z_cav 
    94      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace 
     114     INTEGER                             ::   ji, jj, jk, jp, jl   ! dummy loop index 
     115     INTEGER                             ::   jbku, jbkv, jbkf     ! dummy loop index 
     116     REAL(wp)                            ::   zdepth, e3psm, e3psr ! local variables 
     117     REAL(wp)                            ::   zgdept, zgdepw, zgdepwp1 
     118     REAL(wp)                            ::   ze3t, ze3w, e3u_j, e3u_jp1 
     119     REAL(wp)                            ::   e3t_i, e3t_ip1, e3w_i, e3w_ip1 
     120     REAL(wp)                            ::   e3t_j, e3t_jp1, e3w_j, e3w_jp1 
     121     REAL(wp)                            ::   ps_i, ps_ip1, ps_j, ps_jp1 
     122     REAL(wp)                            ::   ps_uj, ps_ujp1, ps_p, ps_s 
     123     ! 
     124     REAL(wp), POINTER, DIMENSION(:,:)   :: pst_msk              !: mask for areas where to 
     125                                                                 !: apply partial steps 
     126     REAL(wp), POINTER, DIMENSION(:,:)   :: apst                 !: tracking points where 
     127                                                                 !: we apply partial steps 
     128     ! ARRAYS TO BACKUP ORIGINAL MEs COORDINATE SYSTEM  
     129     REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_mes, e3w_mes     !: scale factors [m] 
     130     REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_mes, gdepw_mes !: depths [m] 
    95131     !!---------------------------------------------------------------------- 
    96      ! 
    97      IF(lwp) THEN 
    98        WRITE(numout,*) 
    99        WRITE(numout,*) '   zgr_dom_read : read the vertical coordinates in domain_cfg_global.nc file' 
    100        WRITE(numout,*) '   ~~~~~~~~' 
    101      ENDIF 
    102      ! 
    103      CALL iom_open( 'domain_cfg_global.nc', inum ) 
    104      ! 
    105      !                          !* ocean cavities under iceshelves 
    106      CALL iom_get( inum, 'ln_isfcav', z_cav ) 
    107      IF( z_cav == 0._wp ) THEN   ;   ln_isfcav = .false.   ;   ELSE   ;   ln_isfcav = .true.   ;   ENDIF 
    108      ! 
    109      !                          !* vertical scale factors 
    110      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d_in ) ! 1D reference coordinate 
    111      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d_in ) 
    112      ! 
    113      CALL iom_get( inum, jpdom_global, 'bottom_level' , mbathy_in )     ! 2D mbathy 
    114      ! 
    115      CALL iom_get( inum, jpdom_global, 'e3t_0'  , e3t_in  )     ! 3D coordinate 
    116      CALL iom_get( inum, jpdom_global, 'e3u_0'  , e3u_in  ) 
    117      CALL iom_get( inum, jpdom_global, 'e3v_0'  , e3v_in  ) 
    118      CALL iom_get( inum, jpdom_global, 'e3f_0'  , e3f_in  ) 
    119      CALL iom_get( inum, jpdom_global, 'e3w_0'  , e3w_in  ) 
    120      CALL iom_get( inum, jpdom_global, 'e3uw_0' , e3uw_in ) 
    121      CALL iom_get( inum, jpdom_global, 'e3vw_0' , e3vw_in ) 
    122      ! 
    123      !                          !* depths 
    124      !                                   !- old depth definition (obsolescent feature) 
    125      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  & 
    126         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  & 
    127         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  & 
    128         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN 
    129         CALL ctl_warn( 'zgr_dom_read : old definition of depths and scale factors used ', &  
    130            &           '           depths at t- and w-points read in the domain configuration file') 
    131         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d_in )    
    132         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d_in ) 
    133         CALL iom_get( inum, jpdom_global , 'gdept_0' , gdept_in ) 
    134         CALL iom_get( inum, jpdom_global , 'gdepw_0' , gdepw_in ) 
    135         ! 
    136      ELSE                                !- depths computed from e3. scale factors 
    137         CALL e3_to_depth( e3t_1d_in, e3w_1d_in, gdept_1d_in, gdepw_1d_in )    ! 1D reference depth 
    138         CALL e3_to_depth( e3t_in   , e3w_in   , gdept_in   , gdepw_in    )    ! 3D depths 
    139         IF(lwp) THEN 
    140           WRITE(numout,*) 
    141           WRITE(numout,*) '              GLOBAL reference 1D z-coordinate depth and scale factors:' 
    142           WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" ) 
    143           WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d_in(jk), gdepw_1d_in(jk), e3t_1d_in(jk), e3w_1d_in(jk), jk = 1, jpk ) 
    144         ENDIF 
    145      ENDIF 
    146      ! 
    147      CALL iom_close( inum ) 
    148      ! 
    149    END SUBROUTINE zgr_dom_read 
    150  
    151 ! ===================================================================================================== 
    152  
    153    SUBROUTINE zgr_mes_local 
    154      !!--------------------------------------------------------------------- 
    155      !!              ***  ROUTINE zgr_dom_merge  *** 
    156      !! 
    157      !! ** Purpose :   Create a local MEs grid within a gloabal grid  
    158      !!                using different vertical coordinates. 
    159      !! 
    160      !!---------------------------------------------------------------------- 
    161      INTEGER                      ::   ji, jj, jk ! dummy loop index 
    162      INTEGER                      ::   inum       ! local logical unit 
    163      REAL(wp), DIMENSION(jpi,jpj) ::   s2z_msk    ! mask for MEs-area (=2), 
    164                                                   !          transition area (=1), 
    165                                                   !          z-area (=0) 
    166      REAL(wp), DIMENSION(jpi,jpj) ::   s2z_wgt    ! weigths for computing model 
    167                                                   ! levels in transition area 
    168      !!---------------------------------------------------------------------- 
    169  
    170      IF(lwp) THEN 
    171        WRITE(numout,*) 
    172        WRITE(numout,*) '   zgr_mes_merge : read the vertical coordinates in domain_cfg_global.nc file' 
    173        WRITE(numout,*) '   ~~~~~~~~' 
    174      ENDIF 
    175      ! 
    176      CALL iom_open( 'bathy_meter.nc', inum ) 
    177      CALL iom_get( inum, jpdom_data, 's2z_msk', s2z_msk) 
    178      CALL iom_get( inum, jpdom_data, 's2z_wgt', s2z_wgt) 
    179  
    180      mes_msk(:,:) = 1. 
    181      WHERE (s2z_msk(:,:).eq.0.0)  mes_msk(:,:) = 0.0 
    182  
    183      DO jj = 1,jpj 
    184         DO ji = 1,jpi 
    185            SELECT CASE (INT(s2z_msk(ji,jj))) 
    186              CASE (0) ! global zps area 
    187                   mbathy (ji,jj  ) = mbathy_in(ji,jj) 
    188                   gdept_0(ji,jj,:) = gdept_in(ji,jj,:) 
    189                   gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:) 
    190                   e3t_0  (ji,jj,:) = e3t_in (ji,jj,:) 
    191                   e3w_0  (ji,jj,:) = e3w_in (ji,jj,:) 
    192                   e3u_0  (ji,jj,:) = e3u_in (ji,jj,:) 
    193                   e3v_0  (ji,jj,:) = e3v_in (ji,jj,:) 
    194                   e3f_0  (ji,jj,:) = e3f_in (ji,jj,:) 
    195                   e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:) 
    196                   e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:) 
    197              CASE (1) ! MEs to zps transition area  
    198                   gdept_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdept_0(ji,jj,:) + & 
    199                     &                ( 1._wp - s2z_wgt(ji,jj) ) * gdept_in(ji,jj,:) 
    200                   gdepw_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdepw_0(ji,jj,:) + & 
    201                     &                ( 1._wp - s2z_wgt(ji,jj) ) * gdepw_in(ji,jj,:) 
    202              CASE (2) ! MEs area 
    203                   CYCLE      
    204            END SELECT 
    205         END DO 
    206      END DO 
    207      ! 
    208      ! e3t, e3w for transition zone 
    209      ! as finite differences 
     132 
     133     e3psm = rn_e3pst_min 
     134     e3psr = rn_e3pst_rat 
     135 
     136     CALL wrk_alloc( jpi, jpj, pst_msk, apst) 
     137     CALL wrk_alloc( jpi, jpj, jpk, gdept_mes, gdepw_mes, e3t_mes, e3w_mes) 
     138 
     139     ! Initialise the mask to identify where to apply partial steps 
     140     pst_msk(:,:) = 0.0 
     141     IF(ln_pst_mes)  WHERE (l2g_msk(:,:) == 2.0)  pst_msk(:,:) = 1.0 
     142     IF(ln_pst_l2g)  WHERE (l2g_msk(:,:) == 1.0)  pst_msk(:,:) = 1.0 
     143 
     144 
     145     ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     146     ! find the number of ocean levels such that the last level thickness 
     147     ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 
    210148     DO jj = 1, jpj 
    211149        DO ji = 1, jpi 
    212            IF ( s2z_msk(ji,jj) == 1._wp ) THEN 
    213               DO jk = 1,jpkm1 
    214                  e3t_0(ji,jj,jk)   = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk) 
    215                  e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk) 
    216               ENDDO 
    217               ! Surface 
    218               jk = 1 
    219               e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1)) 
     150           IF ( pst_msk(ji,jj) == 1.0 .AND. bathy(ji,jj) > 0._wp) THEN 
     151              DO jk = jpkm1, 1, -1 
     152                 zdepth = gdepw_0(ji,jj,jk) + MIN( e3psm, e3t_0(ji,jj,jk)*e3psr ) 
     153                 IF ( bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1 
     154              END DO 
     155            END IF 
     156        END DO 
     157     END DO 
     158 
     159     ! =====================================================  
     160     ! Scale factors and depth at T- and W- points 
     161     ! ===================================================== 
     162     ! 
     163     ! Backup to the reference MEs-coordinate 
     164     e3t_mes  (:,:,:) = e3t_0  (:,:,:) 
     165     e3w_mes  (:,:,:) = e3w_0  (:,:,:) 
     166     gdept_mes(:,:,:) = gdept_0(:,:,:) 
     167     gdepw_mes(:,:,:) = gdepw_0(:,:,:) 
     168      
     169     ! Mask identifying points where we apply partial steps: 
     170     !  bottom cell is     a partial cell = 0  
     171     !  bottom cell is not a partial cell = 1 
     172     apst(:,:) = 1._wp 
     173      
     174     DO jj = 1, jpj 
     175        DO ji = 1, jpi 
     176           jk = mbathy(ji,jj) 
     177           ! 
     178           ! We use partial steps only at ocean points where  
     179           ! envelopes differ from the bottom topography  
     180           IF ( pst_msk(ji,jj) == 1.0 .AND. jk > 0 ) THEN     
    220181              ! 
    221               ! Bottom 
    222               jk = jpk 
    223               e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk)) 
     182              IF ( gdepw_0(ji,jj,jk+1) /= bathy(ji,jj) ) apst(ji,jj) = 0      
     183              ! 
     184              IF ( jk == jpkm1 ) THEN          ! max ocean level case 
     185                 zgdepw = bathy(ji,jj) 
     186                 ze3t   = bathy(ji,jj) - gdepw_mes(ji,jj,jk) 
     187                 ze3w   = 0.5_wp * e3w_mes(ji,jj,jk) * ( 1._wp + ( ze3t/e3t_mes(ji,jj,jk) ) ) 
     188                 e3t_0(ji,jj,jk  ) = ze3t 
     189                 e3w_0(ji,jj,jk  ) = ze3w 
     190                 gdepw_0(ji,jj,jk+1) = zgdepw 
     191                 gdept_0(ji,jj,jk  ) = gdept_mes(ji,jj,jk-1) + ze3w 
     192                 gdept_0(ji,jj,jk+1) = gdept_0(ji,jj,jk) + ze3t 
     193                 ! 
     194              ELSE                             ! standard case 
     195                 zgdepw   = gdepw_mes(ji,jj,jk  ) 
     196                 zgdepwp1 = gdepw_mes(ji,jj,jk+1) 
     197                 zgdept   = gdept_mes(ji,jj,jk  ) 
     198                 ze3t     = e3t_mes  (ji,jj,jk  ) 
     199                 ze3w     = e3w_mes  (ji,jj,jk  ) 
     200                 IF ( bathy(ji,jj) <= zgdepwp1 ) gdepw_0(ji,jj,jk+1) = bathy(ji,jj) 
     201                 !gm Bug?  check the gdepw_1d 
     202                 !       ... on jk 
     203                 gdept_0(ji,jj,jk) = zgdepw + ( gdepw_0(ji,jj,jk+1) - zgdepw      ) & 
     204                    &                       * ((zgdept-zgdepw) / (zgdepwp1-zgdepw)) 
     205                 e3t_0  (ji,jj,jk) = ze3t   * (gdepw_0(ji,jj,jk+1) - zgdepw) / (zgdepwp1 - zgdepw) 
     206                 e3w_0  (ji,jj,jk) = 0.5_wp * (gdepw_0(ji,jj,jk+1) + zgdepwp1 - 2._wp * zgdepw) & 
     207                    &                       * (ze3w / (zgdepwp1 - zgdepw)) 
     208                 !       ... on jk+1 
     209                 gdept_0(ji,jj,jk+1) = gdept_0(ji,jj,jk) + e3t_0(ji,jj,jk) 
     210              ENDIF 
     211           ENDIF 
     212        ENDDO 
     213     ENDDO      
     214     ! 
     215     ! set value at mbathy+1 
     216     DO jj = 1, jpj 
     217        DO ji = 1, jpi 
     218           jk = mbathy(ji,jj) 
     219           IF ( pst_msk(ji,jj) == 1.0 .AND. jk > 0) THEN 
     220               e3t_0  (ji,jj,jk+1) = e3t_0  (ji,jj,jk) 
     221               e3w_0  (ji,jj,jk+1) = e3t_0  (ji,jj,jk) 
    224222           END IF 
    225223        END DO 
    226224     END DO 
    227225     ! 
    228      ! MBATHY transition zone 
    229      DO jj = 1, jpj 
    230         DO ji = 1, jpi 
    231            IF ( s2z_msk(ji,jj) == 1._wp ) THEN 
    232               DO jk = 1, jpkm1 
    233                  IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 
    234                  IF( scobot(ji,jj) == 0.e0              ) mbathy(ji,jj) = 0 
    235                  IF( scobot(ji,jj) < 0.e0               ) mbathy(ji,jj) = INT( scobot(ji,jj)) ! do we need it? 
    236               END DO 
    237            END IF 
    238         END DO 
    239      END DO 
    240      ! 
    241      ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0  
    242      ! for transition zone 
    243      ! 
     226     ! ===========================================  
     227     ! Scale factors at U-, V-, UW- and VW- points 
     228     ! =========================================== 
    244229     DO jj = 1, jpjm1 
    245230        DO ji = 1, jpim1 
    246            IF ( s2z_msk(ji,jj) == 1._wp ) THEN 
    247               DO jk = 1, jpk 
    248                  e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         & 
    249                                   MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) /         & 
    250                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  ) 
    251  
    252                  e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         & 
    253                                   MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) /         & 
    254                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  ) 
    255  
    256                  e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         & 
    257                                    MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) /         & 
    258                                    MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  ) 
    259  
    260                  e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         & 
    261                                    MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) /         & 
    262                                    MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  ) 
    263  
    264                  e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)         +       & 
    265                                   MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk)      +       & 
    266                                   MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) +       & 
    267                                   MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) )    /       & 
    268                                   MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  & 
    269                                 +   MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))  ) 
     231           ! 
     232           DO jp = 0,1 
     233              !   
     234              ! e3u at the bottom 
     235              ! ~~~~~~~~~~~~~~~~~ 
     236              jbku = mbathy(ji+jp,jj) 
     237              ! Only ocean points   
     238              IF (pst_msk(ji,jj) == 1.0 .AND. jbku > 0) THEN 
     239                 ! work at bottom topography and level below 
     240                 DO jk = jbku, jbku+1 
     241                    ps_i   = apst(ji  ,jj  )   ! is bottom cell at T(ji  ,jj) point a partial cell or not 
     242                    ps_ip1 = apst(ji+1,jj  )   ! is bottom cell at T(ji+1,jj) point a partial cell or not 
     243                    ! correction for level jk 
     244                    IF ( (jk /= mbathy(ji  ,jj)) .AND. (jk /= mbathy(ji  ,jj)+1) ) ps_i   = 1 
     245                    IF ( (jk /= mbathy(ji+1,jj)) .AND. (jk /= mbathy(ji+1,jj)+1) ) ps_ip1 = 1 
     246                    ! 
     247                    ps_p = ps_i * ps_ip1 
     248                    ps_s = ps_i + ps_ip1 
     249                    ! 
     250                    e3t_i   = e3t_0(ji  ,jj,jk) 
     251                    e3t_ip1 = e3t_0(ji+1,jj,jk) 
     252                    e3w_i   = e3w_0(ji  ,jj,jk) 
     253                    e3w_ip1 = e3w_0(ji+1,jj,jk) 
     254                    ! 
     255                    e3u_0(ji,jj,jk)  = MIN( e3u_0(ji,jj,jk) ,                              & 
     256                     &                        ps_p    * e3u_0(ji,jj,jk)                    & ! no zps cell 
     257                     &                     + (1-ps_p) * ( ps_ip1 * e3t_i + ps_i * e3t_ip1  & ! one of neighb is zps 
     258                     &                                   + (1-ps_s) * MIN(e3t_i,e3t_ip1) ) & ! both neighb are zps 
     259                     &                    ) 
     260                      
     261                    e3uw_0(ji,jj,jk) = MIN(e3uw_0(ji,jj,jk),                               & 
     262                     &                        ps_p    * e3uw_0(ji,jj,jk)                   & ! no zps  
     263                     &                     + (1-ps_p) * ( ps_ip1 * e3w_i + ps_i * e3w_ip1  & ! one of neighb is zps 
     264                     &                                   + (1-ps_s) * MIN(e3w_i,e3w_ip1) ) & ! both side are zps 
     265                     &                    ) 
     266                 END DO 
     267              END IF 
     268              !  
     269              ! e3v at the bottom 
     270              ! ~~~~~~~~~~~~~~~~~ 
     271              jbkv = mbathy(ji,jj+jp) 
     272              ! Only ocean points   
     273              IF (pst_msk(ji,jj) == 1.0 .AND. jbkv > 0) THEN 
     274                 ! work at bottom topography and level below 
     275                 DO jk = jbkv, jbkv+1 
     276                    ps_j   = apst(ji, jj  )   ! is bottom cell at T(ji,jj  ) point a partial cell or not 
     277                    ps_jp1 = apst(ji, jj+1)   ! is bottom cell at T(ji,jj+1) point a partial cell or not 
     278                    ! correction for level jk 
     279                    IF ( (jk /= mbathy(ji  ,jj)) .AND. (jk /= mbathy(ji  ,jj)+1) ) ps_j   = 1 
     280                    IF ( (jk /= mbathy(ji+1,jj)) .AND. (jk /= mbathy(ji+1,jj)+1) ) ps_jp1 = 1 
     281                    ! 
     282                    ps_p = ps_j * ps_jp1 
     283                    ps_s = ps_j + ps_jp1 
     284                    ! 
     285                    e3t_j   = e3t_0(ji,jj  ,jk) 
     286                    e3t_jp1 = e3t_0(ji,jj+1,jk) 
     287                    e3w_j   = e3w_0(ji,jj  ,jk) 
     288                    e3w_jp1 = e3w_0(ji,jj+1,jk) 
     289                    ! 
     290                    e3v_0(ji,jj,jk)  = MIN( e3v_0(ji,jj,jk) ,                              & 
     291                     &                        ps_p    * e3v_0(ji,jj,jk)                    & ! no zps cell 
     292                     &                     + (1-ps_p) * ( ps_jp1 * e3t_j + ps_j * e3t_jp1  & ! one of neighb is zps 
     293                     &                                   + (1-ps_s) * MIN(e3t_j,e3t_jp1) ) & ! both neighb are zps 
     294                     &                    ) 
     295 
     296                    e3vw_0(ji,jj,jk) = MIN(e3vw_0(ji,jj,jk),                               & 
     297                     &                        ps_p    * e3vw_0(ji,jj,jk)                   & ! no zps  
     298                     &                     + (1-ps_p) * ( ps_jp1 * e3w_j + ps_j * e3w_jp1  & ! one of neighb is zps 
     299                     &                                   + (1-ps_s) * MIN(e3w_j,e3w_jp1) ) & ! both side are zps 
     300                     &                    ) 
     301                 END DO 
     302              END IF 
     303           END DO 
     304        END DO 
     305     END DO 
     306 
     307     ! =========================  
     308     ! Scale factors at F-points 
     309     ! =========================    
     310     DO jj = 1, jpjm1 
     311        DO ji = 1, jpim1 
     312           DO jp = 0,1 
     313              DO jl = 0,1 
     314                 jbkf = mbathy(ji+jl,jj+jp) 
     315                 ! Only ocean points   
     316                 IF (pst_msk(ji,jj) == 1.0 .AND. jbkf > 0) THEN 
     317                    ! work at bottom topography and level below 
     318                    DO jk = jbkf, jbkf+1 
     319                       ps_uj   = MIN(apst(ji,jj  ),apst(ji+1,jj  )) ! is bot cell at U(ji,jj  ) pnt a part. cell or not 
     320                       ps_ujp1 = MIN(apst(ji,jj+1),apst(ji+1,jj+1)) ! is bot cell at U(ji,jj+1) pnt a part. cell or not 
     321                       ! correction for level jk 
     322                       IF (      (jk /= mbathy(ji  ,jj  )) .AND. (jk /= mbathy(ji  ,jj  )+1) & 
     323                         & .AND. (jk /= mbathy(ji+1,jj  )) .AND. (jk /= mbathy(ji+1,jj  )+1) ) ps_uj   = 1 
     324                       IF (      (jk /= mbathy(ji  ,jj+1)) .AND. (jk /= mbathy(ji  ,jj+1)+1) & 
     325                         & .AND. (jk /= mbathy(ji+1,jj+1)) .AND. (jk /= mbathy(ji+1,jj+1)+1) ) ps_ujp1 = 1 
     326                       ! 
     327                       ps_p = ps_uj * ps_ujp1 
     328                       ps_s = ps_uj + ps_ujp1 
     329                       ! 
     330                       e3u_j   = e3u_0(ji,jj  ,jk) 
     331                       e3u_jp1 = e3u_0(ji,jj+1,jk) 
     332                       ! 
     333                       e3f_0(ji,jj,jk) = MIN( e3f_0(ji,jj,jk) ,                                & 
     334                        &                       ps_p    * e3f_0(ji,jj,jk)                      & ! no zps cell 
     335                        &                    + (1-ps_p) * ( ps_ujp1 * e3u_j + ps_uj * e3u_jp1  & ! one of neighb is zps 
     336                        &                                  + (1-ps_s) * MIN(e3u_j,e3u_jp1)   ) & ! both neighb are zps 
     337                        &                   ) 
     338                    END DO 
     339                 END IF 
    270340              END DO 
    271            END IF 
     341           END DO 
    272342        END DO 
    273343     END DO 
    274      ! Computing gdep3w 
    275      gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 
    276      DO jk = 2, jpk 
    277         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
    278      END DO 
    279      ! 
    280      ! From here equal to sco code - domzgr.F90 line 2079 
    281      ! Lateral B.C. we apply only for transition zone 
    282      ! since for s- and zps- areas has already been applied 
    283  
    284      CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 
    285      CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 
    286      CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 
    287      CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 
    288      CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 
    289      CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 
    290      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    291  
    292      WHERE (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    293      WHERE (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    294      WHERE (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    295      WHERE (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    296      WHERE (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    297      WHERE (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    298      WHERE (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
    299  
    300      IF( lwp ) WRITE(numout,*) 'Refine mbathy' 
    301      DO jj = 1, jpj 
    302         DO ji = 1, jpi 
    303            IF ( s2z_msk(ji,jj) == 1._wp ) THEN 
    304               DO jk = 1, jpkm1 
    305                  IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    306               END DO 
    307               IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
    308            END IF 
    309         END DO 
    310      END DO 
    311  
    312    END SUBROUTINE zgr_mes_local 
     344     ! 
     345     ! set to z-scale factor if zero (i.e. along closed boundaries) 
     346     DO jk = 1, jpk 
     347        WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_0(ji,jj,jk) 
     348     END DO 
     349     ! 
     350     e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     351     e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
     352     e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
     353     e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
     354     e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
     355 
     356     ! Control of the sign 
     357     IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   e3t_0 <= 0' ) 
     358     IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   e3w_0 <= 0' ) 
     359     IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   gdept_0 < 0' ) 
     360     IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   gdepw_0 < 0' ) 
     361 
     362     CALL wrk_dealloc( jpi, jpj, pst_msk, apst) 
     363     CALL wrk_dealloc( jpi, jpj, jpk, gdept_mes, gdepw_mes, e3t_mes, e3w_mes) 
     364 
     365   END SUBROUTINE zgr_pst_mes 
    313366 
    314367END MODULE zgrmes 
Note: See TracChangeset for help on using the changeset viewer.